Lateglacial Shifts in Seasonality Reconcile Conflicting North Atlantic Temperature Signals

The accelerating flux of glacial meltwater to the oceans due to global warming is a potential trigger for future climate disturbance. Past disruption of Atlantic Ocean circulation, driven by melting of land‐based ice, is linked in models to reduced ocean‐atmosphere heat transfer and abrupt cooling during stadial events. The most recent stadial, the Younger Dryas (YD), is traditionally viewed as a severe cooling centered on the North Atlantic but with hemispheric influence. However, indications of summer warmth question whether YD cooling was truly year‐round or restricted to winter. Here, we present a beryllium‐10‐dated glacier record from the north‐east North Atlantic, coupled with 2‐D glacier‐climate modeling, to reconstruct Lateglacial summer air temperature patterns. Our record reveals that, contrary to the prevailing model, the last glacial advance in Scotland did not occur during the YD but predated the stadial, while the YD itself was characterized by warming‐driven deglaciation. We argue that these apparently paradoxical findings can be reconciled with regional and global climate events by invoking enhanced North Atlantic seasonality—with anomalously cold winters but warming summers—as an intrinsic response to globally increased poleward heat fluxes.

. We also note that several modeling studies have proposed that stadial cooling was dominated by wintertime (Renssen & Isarin, 1997;Schenk et al., 2018) and that YD summers were comparatively warm (Flückinger et al., 2008;Schenk et al., 2018), consistent with the seasonality hypothesis (Denton et al., 2005).While such reports offer tantalizing insight on the dynamics of abrupt climate change, the full expression of YD summer climate in the North Atlantic has not yet been tested from a terrestrial glacial vantage.Thus, the dominant view of the YD remains one of year-round Northern Hemisphere cooling accompanied, via the bipolar seesaw concept, by Southern Hemisphere warming (Cheng et al., 2020;Mangerud, 2021;Rea et al., 2020).
Capitalizing on refined cosmogenic 10 Be production rates, we employed the geologic record of mountain glaciation in Scotland to evaluate Lateglacial ablation-season climate in the mid-latitude North Atlantic basin, with emphasis on the YD and preceding millennia.Recognizing that temperature and precipitation play a role in mass balance, we utilized 2-D glacier-climate modeling to test glacier sensitivity to both climatic variables and identify the likeliest thermal scenarios represented by our 10 Be-dated geologic record.By focusing on terrestrial glaciers, we avoid the potential marine influences of calving and sub-surface melting, thereby delivering a conceptually straightforward record of the cryospheric response to abrupt climate change in the pivotal North Atlantic region.

Geographical Context and Site Descriptions
Following the last glacial maximum (LGM: ∼27-19°ka; Clark et al., 2009), when Scotland was inundated by the British ice sheet (BIS), deglaciation was interrupted by a temporary regrowth of ice referred to locally as the Loch Lomond Readvance (LLR), during which myriad mountain glaciers and the larger West Highland ice field (WHIF) were nourished in upland regions.Representing a major shift in mean North Atlantic climate, the LLR has a long association with the YD.However, this relationship has yet to be confirmed via reproducible, direct dating of glacial deposits.We sampled moraines deposited by three small cirque glaciers in NW Scotland during the LLR, all three of which were first mapped by Sissons (1977) (Cadh' a' Mhoraire, Glas Tholl, Toll an Lochan; Figure 1).Small glaciers, such as those restricted to cirques and mountain valleys, are ideal for this study due to their short response time-years rather than decades-to climate forcing relative to larger glaciers and ice sheets (Stone et al., 1998).To resolve the timing of post-LLR deglaciation, we sampled erratics on the cirque floors up-valley of the terminal moraines (Figures 2-4 and Figure S1 in Supporting Information S1).At Cadh' a' Mhoraire and Glas Tholl we also measured 10 Be in perched erratics outside the LLR moraines to determine when the BIS ice margin retreated to its LLR position.Analytical procedures are detailed in the Methods section.
Our northernmost site, Cadh' a' Mhoraire (57.99°, −5.20°), is a north-facing cirque incised into the NW flank of Ben Mor Coigach and exhibiting a headwall of ∼350 m relief (Figures 1 and 2).Local bedrock is uniformly coarse-grained feldspathic sandstone (arkose) of the Precambrian Torridonian Group, containing rounded quartz pebbles of up to ∼5 cm diameter.Twenty-five km south of Cadh' a' Mhoraire, Toll an Lochan (57.77°, −5.32°) is a NE-facing cirque situated between the peaks of Beinn Dearg Mor and Beinn Dearg Bheag and draining into Strath na Sealga (Figures 1 and 3).The lower (100-200 m elevation) and upper (300-400 m elevation) portions of the cirque are separated by a higher-angle slope of ice-molded bedrock, while the cirque headwall, ∼170-300 m in relief, is characterized by colluviated till, talus, and cliffs.Local lithology is dominated both by sandstones and conglomerates of the Torridonian Group and the underlying Lewisian gneiss, which outcrops on the lower slopes of Beinn Dearg Bheag.Accordingly, glacial erratics analyzed in this study include both lithologies.Our third site, Glas Tholl (57.81°, 5.23°), is situated on the east side of the An Teallach massif, immediately below the 1,062 m-high summit of Bidien a'Ghlas Thuill (Figures 1 and 4).On its southern and western margins, the cirque is bounded by precipitous ∼500 m-high cliffs and talus; the northern cirque margin comprises steep colluviated till and soliflucted weathering debris.As at Cadh' a' Mhoraire, the underlying bedrock-and all boulders sampled at Glas Tholl-is coarse-grained sandstone/conglomerate of the Torridon Group.this study; red: recalculated sites of Ballantyne and Stone (2012)).Also shown is the approximate extent of the WHIF and larger satellite ice masses.

Methods
For mid-latitude glaciers with melt zones, mass balance is dominated by ablation season (summer) air temperature (Anderson & Mackintosh, 2012;Larsen et al., 2015;Oerlemans, 2001Oerlemans, , 2005;;Rupper & Roe, 2008), as verified by the widespread retreat of glaciers today due to global warming (Roe et al., 2020;Zemp et al., 2019).Coupled with accurate mapping and cosmogenic-nuclide surface-exposure dating, relict moraines provide a robust terrestrial record of past changes primarily in summer air temperature.Specifically, surface-exposure ages of erratics on moraine crests indicate the culmination of glacier stability and onset of retreat due to warming; transects of perched erratics reflect periods of net retreat.

Glacial-Geomorphic Mapping and Cosmogenic 10 Be Geochronology
The glacial geomorphology of our three sites was mapped over the course of three field surveys (2015)(2016)(2017), during which we also identified and sampled suitable erratic boulders for cosmogenic 10 Be surface-exposure dating.Specifically, we selected glacially molded boulders of Torridonian sandstone and Lewisian gneiss located in stable positions, with no obvious indication of post-depositional movement.Samples were collected from the top surfaces of boulders using a hammer and chisel; sample locations and elevations were measured with Trimble   then rinsed in water and boiled in 6M HCl to remove metal oxides, and quartz isolated by froth-flotation mineral separation.All samples were treated to successive etches in 5% HF/5% HNO 3 and 2% HF/2% HNO 3 until they comprised pure quartz as established via Inductively Coupled Plasma Optical Emission Spectrometry (ICP-OES).
Beryllium extractions were conducted at the University of Maine cosmogenic nuclide laboratory, where pure quartz was first weighed, spiked with a low-10 Be background beryllium carrier, and dissolved in hot concentrated HF.Beryllium was isolated using ion-chromatography following the procedures developed by the University of Washington laboratory (http://depts.washington.edu/cosmolab/chem.shtml).BeO samples were then mixed with niobium powder and packed into stainless steel targets, after which 10 Be/ 9 Be ratios were measured at the Center for Accelerator Mass Spectrometry (CAMS), Lawrence Livermore National Laboratory, relative to the 07KNSTD3110 standard ( 10 Be/ 9 Be = 2.85 × 10 −12 ).Fifty-two samples were measured in 10 batches between 2015 and 2020.Ten procedural blanks provide 10 Be values ranging from 9,420 to 38,168 10 Be atoms per blank (arithmetic mean = 23,080 ± 9,297 10 Be atoms), corresponding to ∼1.5% of the 10 Be concentration of a typical sample measurement (∼1,000,000 10 Be atoms).Reported concentration uncertainties include the reported analytical uncertainty (1σ) propagated with machine background, procedural blank, and boron correction uncertainties (Table S1 in Supporting Information S1).Population ages are reported with the standard error of the mean propagated with the 2.7% production rate uncertainty (Putnam et al., 2019).
We calculated surface-exposure ages using version 3 of the University of Washington online calculator (https:// hess.ess.washington.edu), in conjunction with the recent local 10 Be production rate and time-independent "St" scaling (Balco et al., 2008;Putnam et al., 2019).Our choice of production rate is based on both the proximity of the Rannoch Moor calibration site to our field area (<160 km) and the straightforward geomorphic and stratigraphic principles on which that calibration was made (Putnam et al., 2019).Being located at the same geomagnetic latitude and longitude as Rannoch Moor, and corresponding to the same Lateglacial period, our three target cirques have experienced a similar intensity of cosmic radiation as the calibration site.Second, the Rannoch Moor rate was calibrated against landforms for which robust independent age control is provided by replicable bracketing 14 C ages (Bromley et al., 2014(Bromley et al., , 2018;;Putnam et al., 2019), resulting in a high degree of similarity with other rates globally that employed robustly dated calibration surfaces (e.g., Kaplan et al., 2011;Kelly et al., 2015;Putnam et al., 2010).In contrast, three earlier published estimates of the Scottish 10 Be production rate are based not on independent dating but on assumed landform ages inferred from traditional interpretations of ice core data (Ballantyne & Stone, 2012;Borchers et al., 2016;Small & Fabel, 2015).Not only do those rates stand out as unusually high relative to the global data set (Phillips et al., 2016) (see Section 5.2), the age assumptions involved negate our ability to test the seasonality hypothesis via 10 Be-dated glacial records.Accordingly, the Rannoch Moor production rate is the most appropriate for this study.
Our choice of scaling reflects the fact that St provides good agreement among cosmogenic-derived ages and independent chronological constraints, both at Rannoch Moor and elsewhere (see Putnam et al. (2019) and references therein).Using any one of the three most commonly applied schemes (St, Lm, LSDn) causes our age calculations to vary by less than 0.5%, reflecting the geographic and temporal proximity of our study area to the calibration site.Thus, our glacial-climatic interpretations are independent of scaling choice.
All calculated ages have been corrected for horizon shielding (http://stoneage.ice-d.org/math/skyline/skyline_in.html).We have not corrected for erosion, however, despite evidence for some boulder surfaces having lost as much as 2 cm to post-depositional weathering (e.g., resistant quartz pebbles protruding from less resistant rock matrix: Figure S2 in Supporting Information S1), because we cannot confirm that erosion rates are uniform among our three sites and multiple samples.Nonetheless, we note that accounting for ∼2 cm of surface erosion makes our calculated exposure ages ∼2% older than reported here.All sample and nuclide measurement data are given in Table 1 and calculated ages in Table 2; statistics for dated landforms are given in Table S1 in Supporting Information S1.

Glacier-Climate Modeling
We employed the two-dimensional, finite-difference numerical glacier model of Kessler et al. (2006) to evaluate plausible climatic scenarios consistent with glacier expansion to mapped and dated landforms.The model considers ice flow by ice deformation and sliding, as well as an optional component that accounts for avalanching from headwalls.Ice thickness for each grid cell is determined via the continuity equation together with the underlying  topographic grid of the model domain.The model defines the climate for each scenario by determining the net mass balance (b z ), which is calculated using an assigned mass-balance gradient (∇b z ), maximum accumulation rate (B z max ), and equilibrium-line altitude (ELA).We specified a ∇b z of 7 m w.e.km −1 based on the median value determined from a set of eight modern glaciers in south-western Norway (∼60° latitude) (Rasmussen & Andreassen, 2005), where the temperate maritime climate is similar to that of the north-western Scottish Highlands.We assigned a modern B z max value of 2 m w.e.yr −1 based on average precipitation rates for NW Scotland determined from the ERA5 Reanalysis (Hersbach et al., 2020), which are also similar to accumulation rates measured from the Norwegian glaciers (Rasmussen & Andreassen, 2005).Recognizing that reconstructed palaeo-precipitation rates for the British Isles are highly variable, the sensitivity tests described below encompass a range of lower B Z max values.Here, we focused our modeling efforts on the Cadh' a' Mhoraire, Toll an Lochan, and Glas Tholl field areas.For each glacier, we performed a set of simulations in which the ELA was lowered in 10-m increments from the top of the cirque headwall until the modeled glacier margin gave a best fit to the targeted landform.Each simulation was run for 1,000-model years, allowing the modeled glacier to achieve equilibrium (i.e., constant volume) several hundred years before the simulation was complete.Most simulated glacier volumes achieved equilibrium within 300 model years.
Following the approach of Kessler et al. (2006), we also conducted a series of sensitivity tests to explore relationships among the three Scottish glacier systems and the past climatological parameters: ∇b z , B z max (i.e., precipitation change; we also refer to this as ∆P) and temperature (∆T).First, we determined the amount of ELA change required to compensate for a change in ∇b z .Results indicated that a minor ELA adjustment of ±10 m is required to compensate for the effects of changing the ∇b z by ±1 m w.e.km −1 (i.e., consistent with the range represented by glaciers of SW Norway) and maintain a good fit between the modeled ice margin and the targeted lateral and terminal moraines.We note that this range of ∇b z produced ice-surface profiles that gave a good match to lateral moraine profiles.
To evaluate sensitivity to ΔP, we repeated the sets of simulations described above but with B z max set to 75%, 50%, 25%, and 10% of the modern value to derive two specific outcomes.First, for each simulation set, we determined the ∆ELA from the B z max = 100% (i.e., modern precipitation) best-fit result that would be required to counteract the reduction in precipitation and produce a good fit with the targeted glacial landform.For each of these model results, any change in ELA must be due to the effect of temperature change (∆T) on glacier mass balance, because accumulation is held constant.These simulations produced a subset of best-fit results (i.e., in which the glacier extends to a targeted mapped and dated landform) from which we produced a set of curves (Figures S3-S5 in Supporting Information S1) exhibiting the non-linear sensitivity of the selected glacier systems to changes in temperature and precipitation (discussed below; Table S2 in Supporting Information S1).Second, we determined how much glacier recession would result from a given decrease in precipitation with the ELA held constant (Figures S3-S5 in Supporting Information S1).The range of prescribed precipitation values used here is admittedly extreme relative to Lateglacial reconstructions for the British Isles; existing estimates for the YD, for example, are equivocal, ranging from marginally drier to slightly wetter than today (Section 5.1).Nonetheless, as detailed in Section 4.1, this conservative approach confirms that our conclusions are not impacted by changes in precipitation.Finally, we also tested whether the effects of avalanching could have a significant impact on the results.The avalanching function included in the Kessler et al. (2006) model transports any accumulation on slopes >30° down to less steep terrain, resulting in a redistribution of mass on the upper part of the glacier.We found that the effects of avalanching have negligible impacts on modeled glacier geometry for the targeted paleo-glacier systems.
Resulting ELA reconstructions have been converted to temperature (∆T) using an adiabatic lapse rate of 6.5°C km −1 , which we feel is the most suitable approach given the absence of paleo-lapse rate estimates for the British Isles.We incorporated the effects of changes in ∇b z and a 50% reduction in ∆P into attendant ∆T uncertainties (Figure 5).

Model Sensitivity Tests
In line with modern glaciological observations, our model sensitivity tests indicate that mass balance of the three reconstructed LLR glaciers was dominated by melting, and thus air temperature.As shown in Figures S3-S5 in Supporting Information S1, the negative impact on mass balance of a 50% reduction in precipitation is balanced by a cooling of approximately 0.1°C, while even extremely dry conditions-10% of current values, equivalent to the modern Sahel-are offset by a ∼0.5°C cooling (Table S2 in Supporting Information S1).We also note that, for all three cirques, sizable glaciers remain under all five precipitation scenarios when temperatures are held stable (Figures S3-S5 in Supporting Information S1).Together, these results confirm that ΔELA reconstructed from the moraine record is a reliable reflection of atmospheric ΔT; to reflect the full range of measured precipitation impacts, we nevertheless apply a conservative ±0.2°C uncertainty to all temperature estimates.

Cadh' a' Mhoraire
Moraines corresponding to the LLR form three principal arcuate ridges deposited by north-flowing ice originating in Cadh' a' Mhoraire and have been described previously by Chandler and Lukas (2017).Outside the LLR moraines, the landscape is dominated by glacially molded bedrock with a thin cover of peat and low vegetation, and perched erratic boulders (e.g., 57,58); moraine forms corresponding to the decaying ice sheet are absent.The outer LLR moraine is the most prominent and complete and forms the southern shore of Lochan Tuath.Average moraine relief is 2-4 m and moraine crests are mantled with abundant embedded boulders derived from the cirque up-valley.Inboard of the inner terminal moraine, the sparsely vegetated cirque floor is dominated by ice-molded bedrock mantled with erratic boulders (e.g., ANT-17-53, 55); moraine ridges and mounds are absent.Physical weathering of the rock surfaces (e.g., erratic boulders, glacially molded bedrock) since the end of the LLR is indicated by the abundance of protruding quartz pebbles (Figure S2 in Supporting Information S1), which have a maximum observed relief of ∼2 cm.
We measured 10 Be in 17 samples from Cadh' a' Mhoraire and the immediate vicinity.Ages of three boulders perched on bedrock outside the Lateglacial moraines record the disappearance of the BIS by 14.0 ± 0.3 ka (Figure 2).The subsequent LLR is represented by three arcuate moraines describing a glacier of 0.6 km 2 at its maximum (Figure 2).The outer moraine gives a mean age of 13.5 ± 0.3 ka (n = 4) and the inner moraine a mean age of 12.7 ± 0.3 (n = 5) (Figure 2).Three samples from the middle moraine include one that returned an age of 13.3 ± 0.4 ka (ANT-17-45), which is stratigraphically concordant with the outer and inner moraines, and two (ANT-17-44, 46) that are younger than all ages on the stratigraphically younger inner moraine and thus we reject as outliers (Figure 2).Accordingly, we take 13.3 ka is the approximate age of the middle moraine.Erratics on the cirque floor 400 m inside and 70 m above the inner moraine give a mean of 12.5 ± 0.3 ka (n = 2) for complete cirque deglaciation (Figure 2).Based on our model runs, the gradual, active retreat of the Cadh' a' Mhoraire terminus from the outer to inner LLR moraines represents an ΔELA of ∼+30 m, driven by a 0.2°C warming, after which abandonment of the inner moraine and terminus retreat into the inner cirque reflects an ΔELA (ΔT) of +50 m (+0.3°C).In all, the complete melting of the Cadh' a' Mhoraire glacier (ΔELA = +100 m) corresponds to a warming of 0.65°C (Figure 5).

Toll an Lochan
The glacial geology of Toll an Lochan is dominated by a suite of three arcuate terminal moraines deposited by a 1.3 km 2 glacier during the LLR (Figure 3).These landforms overlie thick, fluvially incised units of ground moraine presumably deposited by the last ice sheet and lacking any clear ice-marginal limits.Following retreat of the Toll an Lochan LLR glacier, the terminal moraines have been incised by river erosion, resulting in extensive exposures of the ridges' internal structure.The largest of these sedimentologic sections, ∼2-3 m high, is located immediately upstream of the middle terminal moraine suite (which comprises two closely nested moraine ridges) at 57.784°, −5.305° and reveals a sequence of lacustrine sands and clays deposited in a shallow pond/lake impounded by the moraine (Figure 3 and Figure S6 in Supporting Information S1).The back-tilted and, in places, faulted nature of these lacustrine sediments indicates post-depositional disturbance by overriding ice, which also deposited the overlying basal till during a minor readvance of the cirque glacier terminus.Together, the presence of glacio-lacustrine sediments and interbedding of these units with glacial tills, along with the widespread glacial molding of bedrock and erratic boulders at our three sites, describe a dynamic, warm-based glacier regime during the LLR.A low-relief (<1.5 m), bouldery lateral moraine ridge bounding the south side of the high tarn represents the final stage of the much-diminished Toll an Lochan glacier (Figure 3).This landform was not dated and thus is not included in the current chronology.
We measured beryllium in 15 samples collected from erratics embedded in moraine crests and perched on bedrock in Toll an Lochan.The outer and inner moraines give mean ages of 13.3 ± 0.3 ka (n = 5; Figure 3) and 13.0 ± 0.3 ka (n = 4), respectively.One boulder on the middle moraine gives an age of 13.0 ± 0.3 ka (ANT-17-09), which is stratigraphically concordant with the inner and outer moraines; two additional boulders (ANT-17-10, 12) give relatively young ages that are inconsistent with the local chrono-stratigraphy and thus are considered outliers and removed from consideration.The gradual contraction of the Toll an Lochan terminus from the outer to inner LLR moraines corresponds to an ΔELA of ∼+50, driven by a 0.2°C warming (Figure 5).Two erratics on bedrock 200 m higher and 0.75 km up-valley from the inner moraine indicate extensive deglaciation by 12.8 ± 0.4 ka (Figure 3), reflecting a 60 m rise in ELA and warming of 0.4°C.A third deglacial age (13.7 ± 0.3 ka; ANT-17-23) predates its two neighbors by almost a millennium and also predates the three LLR moraines located in stratigraphically older positions at the mouth of Toll an Lochan, and thus is pruned as an outlier.From the outer LLR moraine, complete deglaciation of Toll an Lochan (ΔELA = +220 m) corresponds to a total warming of 1.4°C (Figure 5).

Glas Tholl
The LLR is defined at Glas Tholl by a single composite moraine ridge deposited by a 1 km 2 glacier (Figure 4).The moraine, which in its southern half comprises a high-relief composite ridge, visibly overlies an extensive pavement of glacially molded bedrock riven with relict meltwater channels and mantled with perched boulders deposited by the decaying ice sheet prior to the LLR (Figure S1 in Supporting Information S1).As at the other two cirques, there are no conspicuous moraine landforms preserved outside the Glas Tholl LLR moraine.Inboard of the moraine, the cirque floor is characterized by a thin till sheet, which is conspicuously fluted in its southern portion, and abundant perched erratics (including samples 10,13,and 14).That cirque glaciation in Glas Tholl was wet-based is reflected by the highly molded nature of the bedrock and erratics, and reinforced by the relatively low concentrations of cosmogenic beryllium-10 in bedrock samples measured by Stone et al. (1998) (Figure 4), which indicate effective removal of any nuclide inheritance by basal erosion (7.84-8.45× 10 4 atoms/g).
We measured beryllium in 20 samples from Glas Tholl, including five from erratics on bedrock outside the LLR moraine that confirm retreat of the BIS from this area by 15.0 ± 0.4 ka (n = 5) (Figure 4).Eleven ages from the LLR moraine itself range from 13.6 ± 0.3 ka to 11.6 ± 0.5 ka (n = 11), giving a non-normal bimodal distribution (Figure 4) that is indicative of likely geologic outliers.Finally, four ages from erratic boulders on the cirque floor up-valley of the LLR moraine provide an internally consistent grouping of 12.8 ± 0.3 ka (Figure 4), indicating significant post-LLR deglaciation by that time.We note this age is consistent with three lower-resolution 10 Be dates for the deglaciation of Glas Tholl reported by Stone et al. (1998) (recalculated mean 12.6 ± 0.5 ka; Figure 4).
The range exhibited by the LLR moraine ages presents a challenge to establishing the landform's true age.While the most conservative approach is simply to conclude that moraine deposition occurred between 13.6 and 11.6 ka, the tight grouping afforded by the four up-valley ages (12.8 ± 0.3 ka) provides a valuable measure of stratigraphic control with which we can reduce this uncertainty and arrive at a more realistic age estimate.Assuming that the moraine's true depositional age lies within the range of the 11 measured samples, any exposure age on this moraine that is younger than all ages on a stratigraphically younger surface must be erroneous (see Balter-Kennedy et al., 2020).Following this treatment, we identify six exposure ages on the LLR moraine that are younger than all four ages of recessional erratics up-valley and are thus chrono-stratigraphically discordant.Pruned of these outliers, the population gives a stratigraphically concordant age range of 12.7-13.6ka.Given the voluminous nature of the ridge, we posit that the broad age range might reflect the influence of post-depositional slumping and/or melt-out of internal ice.
Upon the culmination of the LLR, retreat of the Glas Tholl terminus from the moraine to the inner cirque represents a ∼+50 m rise in ELA, driven by a 0.3°C warming; complete deglaciation of Glas Tholl (ΔELA = +155 m) corresponds to a total warming of 1.0°C (Figure 5).

Nature and Timing of Glacier Fluctuations in Scotland
The glacial geomorphology of our study area supports the longstanding model of the LLR as a conspicuous regrowth of ice, rather than a pause in the final stages of ice sheet deglaciation (Ballantyne, 2012).At each of the three cirques, the moraine stratigraphy describes separate, locally nourished glaciers that were unrelated to the distal, pre-LLR geomorphology of earlier ice sheet decay.The distribution of recessional ice sheet moraines throughout our study area does not support a Lateglacial retreat of ice back to the cirques (Chandler & Lukas, 2017;Sissons, 1977).Further support for the LLR as a discrete event is provided by 10 Be ages of pre-LLR erratics at Glas Tholl (Figure 4), which indicate a period up to 1,500 years between ice sheet deglaciation and the subsequent LLR, and by the wealth of sedimentologic data documenting the Lateglacial reoccupation of western fjords by WHIF outlet glaciers following a period of open marine conditions after the LGM (Bromley et al., 2018;Gray & Brooks, 1972;Peacock et al., 1989;Sissons, 1967).
Geomorphic indicators of wet-based glaciation and meltwater confirm that, like the concomitant WHIF, the cirque glaciers investigated here were temperate in nature and thus sensitive primarily to melt season (summer) temperature.Based on our numerical modeling results, as well as on the close coupling between modern temperate glaciers and summer temperature (Zemp et al., 2019), we postulate the LLR occurred in response to relatively cool summers, and hence lower melt rates, conducive to survival of winter snowpack and ice accumulation.
Our 10 Be-dated moraine record suggests the LLR occurred as early as 13.4 ka, which is the average of outer moraines at Cadh' a' Mhoraire and Toll an Lochan, and probably no later than 13.1 ka, which represents the average of all LLR moraines (minus outliers).Replicated at Cadh' a' Mhoraire and Toll an Lochan, and the most stratigraphically plausible scenario at Glas Tholl, we propose that this geologic signal shows that Scotland's last glaciers culminated during the Bølling-Allerød (B-A) period and not the YD (Figure 6).This pattern aligns with 10 Be ages from boulders located on drift deposited by an LLR cirque glacier in Coire an Lochain when recalculated in the same manner (Kirkbride et al., 2014) (n = 9; peak age = 13.8 ka; mean age = 14.3 ± 1.0 ka; Figure S7 in Supporting Information S1).
According to our data, the LLR predates previous estimates by as much as 1.5 kyr (Lowe et al., 2019;MacLeod et al., 2011;Palmer et al., 2020;Small & Fabel, 2016).Nonetheless, we note that the last expansion of Scottish glaciers aligns with events reported from Arctic (Wittmeier et al., 2020), tropical (Bromley et al., 2011;Jomelli et al., 2014;Martin et al., 2020), and southern mid latitudes (Putnam et al., 2010;Sagredo et al., 2018), where major Lateglacial advances predated the YD.A plausible mechanism for this widespread return to positive mass balance is explored in Section 5.3.
'Steady-state' ELAs for the three cirques ranged from 385 to 695 m, in broad agreement with earlier reconstructions (Sissons, 1977).Regional variability in absolute ELAs potentially reflect added influences of aspect, topography, avalanching, and windblown snow on the glacier mass balance and sensitivity (Barr & Lovell, 2014;Nesje & Dahl, 2000).For instance, Chandler and Lukas (2017) reported that the Cadh' a' Mhoraire glacier experienced topographically enhanced snow accumulation, thereby depressing the ELA below the "true" regional snowline.Such variability in palaeo-ELA values parallel the glaciologic variability observed today among neighboring glaciers in, for example, the Southern Alps, European Alps, and Norway (Chinn et al., 2005;Rabatel et al., 2013;Torsnes et al., 1993).Following a period during which glacier termini fluctuated within 200-300 m of their maximum LLR extents, the three cirque glaciers underwent marked retreat during the early YD (Figures 5 and 6) and, at Cadh' a' Mhoraire, had largely melted away by 12.6 ka.This pattern is mirrored by recalculated 10 Be ages from three nearby cirques-Maol Chean Dearg, Coire Fearchair, and Coire nan Arr (Figure 1)-where erratics up-valley of the LLR moraines represent final cirque deglaciation (Figures S8-S10 in Supporting Information S1; Ballantyne & Stone, 2012).We include these three sites in our assessment for two specific reasons: first, each cirque affords an internally consistent data set with minimal scatter, making them comparable to our own ages; second, the original data were generated via a methodology-sample preparation, AMS measurement and standardization-comparable to our own, thus enabling a true apples-to-apples age comparison.Recalculated in the same manner as our data, we obtain mean ages for deposition of 12.5 ± 0.3 ka (Maol Chean Dearg), 12.5 ± 0.3 ka (Coire Fearchair) and 12.7 ± 0.3 ka (Coire nan Arr) (Figures S8-10 in Supporting Information S1), revealing that deglaciation occurred ∼0.3 kyr after the YD onset and ∼1 kyr earlier than the end of the YD.Reflecting the variable depression of LLR ELAs described above, the magnitude of ΔELA (and thus ΔT) required for full deglaciation also varies among cirques (Figure 5).
Altogether, our 10 Be data indicate coeval moraine construction and retreat among the multiple cirques described here, leading us to conclude that our composite chronology describes a regional pattern of pre-YD moraine construction and deglaciation throughout the stadial itself.This behavior mirrors that of the larger, concurrent WHIF reconstructed from bracketing 14 C ages (Bromley et al., 2014(Bromley et al., , 2018)), suggesting that Scottish ice masses large and small fluctuated broadly in unison during the Lateglacial.As supported by our model experiments, we implicate rising summer temperature as the cause of ice retreat, plausibly linked to the resumption of global warming at the YD onset (Barker et al., 2009;Bereiter et al., 2018;Denton et al., 2021;Marcott et al., 2014).Recognizing that precipitation also plays a role in mass balance, our model sensitivity tests indicate that the magnitude of desiccation required to drive deglaciation to completion (more than 90% below modern precipitation values) is wholly inconsistent with published YD estimates for Scotland (Ballantyne, 2002;Benn & Ballantyne, 2005;Benn & Lukas, 2006;Chandler & Lukas, 2017;Karger et al., 2021).Thus, we consider that the Scottish mountain glacier record primarily reflects changes in ablation season temperature.
Glaciologically, our Scottish cirque record aligns with behavior of the adjacent WHIF (Bromley et al., 2014(Bromley et al., , 2018;;Putnam et al., 2019) and with mountain glaciers elsewhere (Jomelli et al., 2014;Martin et al., 2020;Putnam  Summed probabilities and mean ages for all moraine samples (blue curve; minus outliers) from our three cirques and all deglacial erratics (red curve) from our sites plus the three recalculated CRONUS sites.Curves normalized to Younger Dryas onset (12,850 years: Rasmussen et al., 2006).Blue and white/black vertical lines are mean and standard error of each population, respectively.et al., 2010;Sagredo et al., 2018).In the North Atlantic region, our record also mirrors the YD retreat of Arctic glaciers (Carlson et al., 2021;Funder et al., 2021;Lakeman et al., 2018;Levy et al., 2016;Wittmeier et al., 2020) and coincided with the resulting influx of meltwater to the ocean (Bard et al., 2000;McManus et al., 2004;Murton et al., 2010).However, while prior studies evoked marine mechanisms to explain deglaciation under apparently cold conditions-including bathymetry and channel morphology and sub-aqueous melting (Lakeman et al., 2018;Funder et al., 2021)-the cirque glaciers in our study were entirely terrestrial and thus immune to hypothesized marine effects.We reiterate that the likeliest mechanism forcing glacier retreat was atmospheric warming during the melt season.From a glacial vantage, events in Scotland bear a striking resemblance to the more global pattern of Lateglacial advance during, but not necessarily linked to, the Antarctic Cold Reversal followed by YD retreat (Figure 6).

Comparison With Existing Models of the Scottish Lateglacial
Acknowledging that our interpretation of the LLR contrasts with the traditional view, in which glaciers reformed in response to depressed YD temperatures and reached their maxima either early or late in the stadial (Ballantyne et al., 2007(Ballantyne et al., , 2012;;Lowe et al., 2019;Macleod et al., 2011;Palmer et al., 2020;Small & Fabel, 2016), we now appraise the empirical basis for that long-standing model in the context of our data.Amassed over more than half a century, previous attempts to pin the LLR in time can be sorted into two broad groups: those employing direct dating methods, such as cosmogenic nuclide geochronology, and those based on indirect, and in some cases assumed, stratigraphic relationships.The former primarily comprises 10 Be surface-exposure ages of erratic boulders and scoured bedrock produced by the action of LLR glaciers themselves, and is potentially the most valuable in establishing the timing of the glacial event directly.Since the application of 10 Be dating in Scotland spans more than two decades, however, the cumulative published data set exhibits considerable variability arising from methodological advances and our evolving understanding of nuclide production rates.For example, earlier studies applying 10 Be production rates ∼20% higher than current assessments reported exposure ages that are now known to be significantly younger than the corresponding landform's true age (Table S3 in Supporting Information S1).Accordingly, upon recalculation with a more accurate rate, many glacial landforms originally assigned to the YD now predate the stadial; notable examples include LLR chronologies from Orkney (Ballantyne et al., 2007), Beinn Inverveigh (Golledge et al., 2007;see Putnam et al., 2019), and the Cairngorms (Kirkbride et al., 2014) (Figure S7 and Table S3 in Supporting Information S1).Similarly, recalculated 10 Be ages of perched erratics and scoured bedrock on cirque floors signify that deglaciation occurred not at the close of the YD, as originally assumed (Ballantyne & Stone, 2012;Stone et al., 1998), but early in the stadial (Putnam et al., 2019;Figures S8-10 in Supporting Information S1).
The pronounced shift in apparent landform age has key ramifications for our understanding of the LLR and its climatic drivers, underscoring the importance of robust computation of production rates.While several attempts have been made to constrain production rates in Scotland (Ballantyne & Stone, 2012;Small & Fabel, 2015;Stone et al., 1998), the Rannoch Moor calibration used in this study remains the only published work to have calculated the in-situ production of cosmogenic 10 Be from an independently dated landform and without underlying assumptions of how glacier activity is tied to the Greenland ice core signal (Putnam et al., 2019).
Where dating control has been sought, indirect approaches to correlating the LLR with the YD have focused largely on radiocarbon ages of organic macrofossils in glacial contexts (Lowe et al., 2019;Macleod et al., 2011;Rose et al., 1988), tephrochronology (Macleod et al., 2015;Palmer et al., 2020), and varve chronologies from ice-dammed lakes (Palmer et al., 2010(Palmer et al., , 2012(Palmer et al., , 2020)).We note, however, that one such recent study essentially replicated the minimum-limiting 14 C ages reported by Bromley et al. (2014) for the K2 core site on Rannoch Moor (Lowe et al., 2019).Those ages comprise part of a larger data set that not only indicates early YD deglaciation of Rannoch Moor, but which also forms the basis for the 10 Be production rate employed in this study (Putnam et al., 2019).
Elsewhere, a growing body of sedimentary data documents the former existence of lakes impounded by LLR glaciers, the best known of which-the Glen Roy-Spean-Laggan lake system-was dammed by the WHIF (Sissons, 1978).According to the Lochaber Master Varve Chronology (Macleod et al., 2015;Palmer et al., 2010Palmer et al., , 2012Palmer et al., , 2020)), this lake existed for approximately 500 years, during which Icelandic tephras became incorporated into the varved lake sediments.Pertinent to our study is the assertion by Palmer et al. (2020) that cryptotephra preserved in varves at the Loch Laggan West (LLW) site corresponds to the Vedde Ash, which erupted from the Katla volcano approximately 12.1 ka (Lane et al., 2012), in which case both the lake and damming ice field existed during the YD (Palmer et al., 2020).We note, however, that the geochemistry of the LLW tephra also matches that of the R1 tephra, which erupted from Katla approximately 13.6 ka and is also evident in the North Atlantic region (Thornalley et al., 2011).In the absence of conclusive dating control (see Lane et al., 2012), therefore, the LLW tephra might just as readily be reinterpreted to suggest an Allerød age for the Glen Roy-Spean-Laggan lake system, consistent both with the Rannoch Moor deglacial 14 C chronology (Bromley et al., 2014;Putnam et al., 2019) and with our new 10 Be-dated moraine record.
Finally, a 14 C-dated stratigraphic section at Croftamie, near the terminus of the former Loch Lomond glacier, is widely cited as robust evidence for the LLR culminating late in the YD (Benn, 2021;Lowe et al., 2019;Macleod et al., 2011).This interpretation assumes that the overlying "Gartocharn till" was deposited as lodgement till by the advancing glacier (Coope & Rose, 2008;Rose, 1981;Rose et al., 1988), in which case the ten 14 C dates on underlying plant remains afford a maximum-limiting age for the advance (Macleod et al., 2011;Rose et al., 1988).Yet a demonstrably (sub)glacial origin for the Gartocharn unit remains elusive and, as elaborated by Bromley et al. (2018), this hypothesis is undermined by key glacial-geomorphic inconsistencies that preclude the Croftamie stratigraphy from defining the timing of the LLR.
Beyond Scotland, the pattern of Allerød advance and stadial retreat described by our record is also evident in neighboring Scandinavia, where Lateglacial deposits have been mapped extensively (see Mangerud et al., 2022).Moraine chronologies from both southern and Arctic Norway, constructed using comparable methods to those employed in this study, indicate culmination of Lateglacial advance prior to (Wittmeier et al., 2020) and at the YD transition (Andersen et al., 1995), with net recession during the stadial.Farther east, bracketing 10 Be ages place construction of the Middle Swedish End Moraine not within the YD, as previously assumed, but during the Bølling-Allerød (Anjar et al., 2014;Cuzzone et al., 2016), while 10 Be ages from the Salpausselkä I moraines in southern Finland constrain abandonment to the late Allerød (see Cuzzone et al., 2016, and references therein).That the ice sheet margin retreated during the YD is confirmed by 10 Be ages from the younger Salpausselkä II moraine (Cuzzone et al., 2016).Collectively, these data agree with the Scottish chronologies presented here and are further supported by relative sea-level evidence for early Allerød expansion of the Scandinavian ice sheet (Lohn et al., 2007).

Implications for Regional, Hemispheric, and Global Stadial Climate
Mechanisms proposed to explain the apparent paradox of North Atlantic deglaciation during stadial conditions include marine instability (Callard et al., 2020;Small et al., 2017), anomalous warm sub-surface currents (Jennings et al., 2017;Knutz et al., 2011), and enhanced aridity (Allard et al., 2021).Each of these mechanisms invokes specific dispensations that, potentially plausible on a local basis, become irreconcilable with the basinwide picture of stadial deglaciation.An alternative view hypothesizes that North Atlantic stadials, as a predictable by-product of climate warming, can only be sustained by elevated meltwater from fringing ice masses, and thus are a symptom of strong summertime warming (Denton et al., 2021(Denton et al., , 2022)).By this model, the intensified discharge of meltwater into the North Atlantic documented during stadials generally (e.g., Boswell et al., 2019;Martinez-Lamas et al., 2020;Pellegrini et al., 2018;Soulet et al., 2013;Toucanne et al., 2009Toucanne et al., , 2015) ) serves to weaken AMOC, foster expansion of seasonal sea ice, and halt the ocean-atmosphere heat exchange that ordinarily supports mild winter temperatures downwind.The local result is anomalous wintertime cooling.Moreover, because the cold anomaly must be sustained by elevated summertime melting, a strong seasonality signal emerges in typically low-seasonality maritime regions.If true for the Lateglacial, the profound retreat of Scottish glaciers during the YD stadial was a product of global warming that, regionally, also produced concomitant winter cooling.
New research suggests that the energy driving stadial deglaciation may have originated in the tropics, where invigorated transfer of water vapor and latent heat to the troposphere delivered a rapid and global rise in radiative forcing (Denton et al., 2022).Building upon that model, we propose that the LLR, as part of a global pattern of Lateglacial resurgence, reflects the temporary cessation of the widespread stadial warming trend that also coincided with flatlining greenhouse gas forcing (Marcott et al., 2014), a pause in CO 2 outgassing (Anderson et al., 2009), declining mean ocean temperature (Bereiter et al., 2018), and the Antarctic Cold Reversal.Within the North Atlantic basin, the reduction in meltwater production during the LLR permitted reinvigoration of AMOC and reinstated a maritime climate of low thermal seasonality in the British Isles.
From this global glacial vantage, the resumption of stadial conditions after 12.9 ka was by no means unique.The documented rise in radiative forcing and poleward energy fluxes (Denton et al., 2022) renewed the influx of meltwater to the North Atlantic (Hillaire-Marcel et al., 2013), impeding North Atlantic deep water and weakening AMOC for the duration of the YD (Brauer et al., 2008;Ezat et al., 2014).During winter, anomalous expansion of sea ice imposed a quasi-continental winter climate on the North Atlantic margin (Buizert et al., 2014); the YD expansion of permafrost and collapse of floral and faunal communities reflects this winter climate anomaly on land (Isarin, 1997).This view is supported by recent TraCE-21ka modeling, which, despite depicting summertime YD cooling, suggests AMOC primarily affects winter climate via sea ice anomalies (Buizert et al., 2018;Denton et al., 2021).Nonetheless, the overall warming trend continued to dominate North Atlantic summers, maintaining the meltwater flux, and pushing deglaciation of Scotland to completion in the early YD.
The apparent discordance between our chronology and existing proxy signatures might be reconciled by considering the importance of seasonality in abrupt climate change (Björck et al., 2002;Denton et al., 2005).By the model described here, the strong cooling traditionally associated with the YD was a winter anomaly sustained by the summer meltdown of Northern Hemisphere terrestrial ice and focused on the North Atlantic, where ocean circulation is sensitive to salinity-buoyancy disturbance.If true, the proxies used regularly to inform the prevailing view of YD summer cooling are potentially skewed toward winter temperature (Brauer et al., 2008;Denton et al., 2005Denton et al., , 2021Denton et al., , 2022)); it is also conceivable that some climate indicators present a more complex response to shifting seasonality than is currently understood (Luoto & Nevalainen, 2013;Schenk et al., 2020;Tóth et al., 2012).

Conclusions
We measured cosmogenic 10 Be in glacially transported boulders in three Scottish cirques to reconstruct the timing of the LLR and subsequent deglaciation.Together, our glacial chronologies indicate that the LLR culminated late in the Allerød period, and not during the YD stadial as previously inferred, and that the YD itself was characterized by moraine abandonment and deglaciation.Glacier-climate modeling indicates that the mass balance of these temperate mountain glaciers was most sensitive to ablation-season air temperature, leading us to conclude that YD recession was driven by atmospheric warming during the summer despite anomalous wintertime cooling that defined the North Atlantic region.These results are consistent with the emerging picture of stadials as periods of enhanced thermal seasonality resulting from the response of North Atlantic ocean-atmosphere heat transfer to global tropospheric warming.

Figure 1 .
Figure 1.Cirque locations (a) NW Scotland study area; yellow rectangle is the area in Panel (b) White shading depicts the approximate extent of the northernmost West Highland ice field (WHIF).(b) Cirque locations (yellow: this study; red: recalculated sites of Ballantyne and Stone (2012)).Also shown is the approximate extent of the WHIF and larger satellite ice masses.

Figure 2 .
Figure 2. Glacial geomorphology and 10 Be chronology for Cadh' a' Mhoraire.Geomorphology and 10 Be ages for Cadh' a' Mhoraire.All statistics in Table S1 in Supporting Information S1.Red curves are summed probabilities for each group.Solid black lines are individual samples; dashed black lines are stratigraphic outliers (italics on map.) Yellow columns represent 1σ standard error of the mean.Imagery from Landmap (2014).

Figure 3 .
Figure 3. Glacial geomorphology and 10 Be chronology for Toll an Lochan.Red star denotes location of the sediment exposure detailed in Figure S3 in Supporting Information S1.Imagery from Landmap (2014).
GNNS.Full sample details, including metadata, are stored online in the NOAA NCEI Paleoclimate data repository viaBromley et al. (2022).Samples were prepared to the clean quartz stage at the University of Maine and University of Galway.First, mass-weighted sample thicknesses were measured using calipers, after which samples were crushed, pulverized, and sieved to retrieve the 250-710 μm size fraction.Samples were

Figure 4 .
Figure 4. Glacial geomorphology and 10 Be chronology for Glas Tholl.Red circles are sample sites of Stone et al. (1998).The likeliest age range for the LLR moraine is denoted by the horizontal black bar and vertical gray shading, while black and blue lines represent retained and pruned ages, respectively.Imagery from Landmap (2014).

Figure 5 .
Figure 5.Time versus relative temperature.Comparison of glacier-inferred temperature change for the NW Scotland cirques to the NGRIP δ 18 O record (NGRIP North GRIP Members, 2004).Temperature reconstructions are based on the moist adiabatic lapse rate (6.5°C/km), while the thickness of the temperature curves (purple) incorporates the ±0.2°C uncertainty applied to our reconstructions based upon the sensitivity experiments.Green triangles (and 1σ error bars) represent 10 Be-dated moraines, the ages of which mark the culmination of glacial equilibrium and the onset of moraine abandonment.Red circles are perched erratics on bedrock and represent uninterrupted retreat of a glacier in disequilibrium with climate.Yellow shading denotes the period encompassing the Younger Dryas stadial.

Figure 6 .
Figure 6.Timing of Lateglacial maxima and recession in NW Scotland.Summed probabilities and mean ages for all moraine samples (blue curve; minus outliers) from our three cirques and all deglacial erratics (red curve) from our sites plus the three recalculated CRONUS sites.Curves normalized to Younger Dryas onset (12,850 years:Rasmussen et al., 2006).Blue and white/black vertical lines are mean and standard error of each population, respectively.

Table 2
(Putnam et al., 2019)ace-Exposure Ages Calculated Using the Rannoch Moor Production Rate(Putnam et al., 2019)and Three Scaling Models

Table 2 Continued
Note.Values in bold denote ages calculated with our prefered scaling scheme.