A cosmogenic nuclide-derived chronology of pre-Last Glacial Cycle glaciations during MIS 8 and MIS 6 in northern Patagonia

. The precise environmental mechanisms controlling Quaternary glacial cycles remain ambiguous. To address this problem, it is critical to better comprehend the drivers of spatio-temporal variability in ice-sheet evolution by establishing reliable chronologies of former outlet-glacier advances. When spanning multiple glacial cycles, such chronologies have the capacity to contribute to knowledge on the topic of interhemispheric phasing of glaciations and climate events. In southern 20 Argentina, reconstructions of this kind are achievable, as Quaternary expansions of the Patagonian Ice Sheet have emplaced a well-preserved geomorphological record covering several glacial cycles. Moreover, robust ice-sheet reconstructions from Patagonia are powerful barometers of former climate change, as Patagonian glaciers are influenced by the Southern Westerly Winds and its coupled Antarctic Circumpolar Current. Former shifts in these circulation mechanisms are essential to better constrain, as they may have played a critical role in pacing regional and possibly global Quaternary climate change. Here, we 25 present a new set of cosmogenic 10 Be and 26 Al exposure ages from pre-Last Glacial Cycle moraine boulder, glaciofluvial outwash cobble and bedrock samples. This dataset constitutes the first direct chronology dating pre-LGM glacier advances in northern Patagonia, and completes our effort to date the entire preserved moraine record of the Río Corcovado valley system (43°S, 71°W). We find the outermost margins of the study site depict at least three distinct pre-Last Glacial Cycle stadials occurring around 290-270 ka, 270-245 ka, and 130-150 ka. Combined with the local LGM chronology, we discover that a 30 minimum of four distinct Pleistocene stadials occurred during Marine Isotope Stages eight, six, and two in northern Patagonia. Evidence for Marine Isotope Stage four and three deposits were not found at the study site. This may illustrate former longitudinal and latitudinal asynchronies in Patagonian Ice Sheet mass balance during these Marine Isotope Stages.


Introduction
Understanding the environmental mechanisms driving the Quaternary build-up and collapse cycles of Earth's major ice sheets remains a focus of ongoing research.Key to this research is the establishment of detailed chronologies offering direct dating of major ice-sheet expansion and recession events (Darvill et al., 2016).Numerical glacier chronologies are most informative when presenting unambiguous dating of ice-sheet margins spanning the entirety of the last glacial cycle (LGC: ~125-15 ka) and ideally earlier Quaternary glacial cycles (Kaplan et al., 2010;Schaefer et al., 2015).To the east of the Patagonian Andes, the Argentinian foreland contains one of the most complete and well-preserved sequences of Quaternary glacial deposits in the world (Clapperton, 1993).This unique geomorphological record provides the opportunity to reconstruct and date glacier fluctuations of the formerly 2500 km-long (from north to south) Patagonian Ice Sheet (PIS) over several glacial cycles (Mercer, 1976, Hein et al., 2011).Furthermore, the PIS occupied a key mid-latitude position in the ocean-dominated Southern Hemisphere (SH) and was the only ice mass to fully intersect the precipitation-bearing southern westerly wind belt, a fundamental feature of the southern mid-latitude climate system.Hence, reconstructing the timing of PIS expansion and recession events during several glacial cycles can provide rare insight into southern mid-latitude Quaternary climate evolution (Davies et al., 2020).
Palaeo-climate proxy records such as those retrieved from ice-and ocean-sediment cores indicate that some of the most potent Quaternary cooling events occurred during the Middle Pleistocene (Parrenin et al., 2013;Shakun et al., 2015), a period (~800 -130 ka; Hughes et al., 2020) characterized by 100-ka pacing of glacial-interglacial cycles.Paradoxically, there is a dearth of directly-dated terrestrial glacial records for these glacial events (Hughes et al., 2020).In Patagonia, bracketing ages for pre-LGC glaciations have been produced in some places (e.g.near Lago Argentino; 50.5°S) with 40 Ar/ 39 Ar dating, K-Ar dating, and palaeo-magnetic measurements on lava flows interbedded between Pliocene/Pleistocene glacial-till units (Mercer, 1976;Sylwan et al., 1991;Meglioli, 1992;Singer et al., 2004).For instance, Mercer (1976) used such methods to constrain the likely timing of the Greatest Patagonian Glaciation to between 1.2 Ma and 1.0 Ma in the Río Gallegos (52°S) and Lago Buenos Aires (46.5°S) regions.While such studies have been instrumental in establishing the broad chronological framework for Patagonian glaciations (Rabassa and Coronato, 2009), direct numerical dating of individual outlet-glacier advances is now required to develop detailed knowledge of pre-LGM and pre-LGC ice-sheet activity in Patagonia.In Patagonia, direct numerical dating of pre-LGC glacial advances, achievable using terrestrial cosmogenic nuclide (TCN) surface exposure dating, has been limited to two valleys of eastern central Patagonia (Kaplan et al., 2005;Hein et al., 2009;2011;2017;Cogez et al., 2018).Throughout the rest of the southern mid-latitudes (23 -66 °S), terrestrial glacial records yielding TCN exposure ages that pre-date the LGC have been reported using moraine boulders in New Zealand (Putnam et al., 2013), the southern central Andes (Terrizzano et al., 2017) and Tasmania (Barrows et al., 2002;Kiernan et al. 2004;2010;2014;2017;Colhoun andBarrows, 2011, Augustinus et al., 2017).Because pre-LGC moraines and boulders are often poorly preserved, the majority of these investigations have identified the approximate timing of individual pre-LGC glacial advances with small (< n = 5 per moraine) and scattered sets of ages.Consequently, robust and direct numerical chronologies of pre-LGC glacial expansion events remain particularly scarce in the southern mid-latitudes.Patagonian investigations have shown that TCN exposure dating can yield direct ages for the deposition of pre-LGC moraine-outwash complexes preserved in the Argentinian steppe (e.g.Hein et al., 2017).While moraine boulders are well-suited to dating glacier advances of the LGC, Hein et al. (2009) revealed that targeting outwash terraces in eastern Patagonia can be more appropriate for dating pre-LGC glacial margins, because their low-gradient surfaces are less prone to degradation via gravity-driven diffusion when compared to steep-sided and unconsolidated moraines.On pre-LGC moraines, such degradation more frequently causes moraine surface lowering and boulder exhumation post-deposition, leading to significant exposure-age scatter and underestimations of moraine formation age.In contrast, surface cobbles from well-preserved outwash plains tend to display comparatively less post-depositional scatter and are thought to date the glacial event more closely (Hein et al., 2009).This approach has been used successfully to constrain the timing of glacial advances in several locations of central and southern Patagonia (Darvill et al., 2015;Hein et al., 2009;2011;2017;Mendelová et al., 2020).However, to date, no pre-LGC glacial events have been dated using this technique in northern Patagonia.Consequently, here we adopt this approach and target both moraine boulders and outwash surface cobbles to constrain the timing of pre-LGC glacial advances in northeastern Patagonia.
We focus our investigation on the Río Corcovado (RC) valley system (43°S; 71°W), formerly host to a major outlet glacier of the PIS in northeastern Patagonia (Caldenius, 1932).We present a multi-glacial-cycle 10 Be and 26 Al TCN exposure-age chronology from the three stratigraphically oldest glacial moraine/outwash complexes deposited by the RC outlet glacier.
This study builds on our earlier work (Leger et al., 2020(Leger et al., , 2021a,b) ,b) to compile a geochronological reconstruction that spans the entire preserved moraine-outwash record of the RC valley system.We use our cosmogenic nuclide-derived chronology to test whether major middle Pleistocene advances previously dated in central Patagonia (Hein et al., 2009;2017) coincide with moraines found in its northeastern sector (this study), and thus whether such advances depict strong regional climate signals expressed across much of the palaeo-ice front.We then compare the different magnitudes of middle Pleistocene Patagonian glaciations relative to other formerly glaciated regions of the world, and thereby assess the presence of former interhemispheric and inter-regional a/synchronies in ice-sheet behaviour.Finally, we evaluate the synchronicity of PIS expansion events with other global and regional palaeo-climate proxy records, and explore the potential role of various insolation parameters on controlling middle-to-late Pleistocene climate variations in Patagonia and the southern midlatitudes.

Setting
Our study site focuses on the ice-marginal environment of the former PIS during Quaternary ice-age maxima, at the latitude of its major eastward-flowing Palena outlet glacier (43°S; Leger et al., 2020;2021a).This up to 100 km-long outlet glacier diverged towards the eastern mountain front into three glaciers occupying the Río Frío valley to the northeast, the Río Huemul (RH) valley to the east, and the RC valley to the south (Fig. 1; Leger et al., 2021a).The local valleys are characterised by a westward-sloping bed topography such that the former outlet glaciers flowed up along reverse-sloping beds.The surface geology of the study area is characterised by well-preserved and distinct sequences of ice-contact glaciogenic deposits interspersed among their associated glaciolacustrine and glaciofluvial sediment-landform assemblages (Leger et al., 2020).Terminal moraine complexes indicate that the RC and RH outlet glaciers re-converged to deposit moraine sequences located approximately 60 km to the east of the Argentinian town of Corcovado (43°54′S; 71°46′W), into an environment presently climatically defined as the semi-arid Patagonian steppe (Fig. 1; Garreaud et al., 2013).Local aridity (~570 mm a -1 ; Fick and Hijmans, 2017; Fig. 1a) results from the Patagonian Andes acting as a powerful orographic barrier to the dominant Southern Westerly Winds (Fig. 1a), causing a potent west-east rain-shadow effect (Garreaud et al., 2013;Fig. 1a).Along the flowlines of the RC outlet glacier, ice flow eroded bedrock mainly sourced from the Mesozoic north Patagonian batholith towards the west (Hervé et al., 2017) along with formations of Jurassic/Cretaceous volcanic and sedimentary units located farther east (Haller et al., 2003).Moraine-outwash complexes studied here thus contain a mixture of rock types including quartz-bearing leuco-monzonite, granite and white granodiorite that are well-suited to TCN exposure dating using in situ 10 Be and 26 Al radionuclides.
Ice sheet-wide geomorphological mapping of glacial sediment-landform assemblages has been conducted by Caldenius (1932), Glasser and Jansson (2008), and Davies et al. (2020) while the first detailed glacial geomorphological map specific to the RC and RH valleys was generated by Leger et al. (2020).These maps enabled the identification of at least eight stratigraphically-distinct moraine-outwash complexes deposited by the RC outlet glacier (Fig. 1c).This study focuses on dating results from the three stratigraphically oldest glacial margins: i.e. the RC 0 (outermost), RC I and RC II (Fig. 1c) moraine-outwash complexes.Directly inboard of these margins, the RC outlet glacier deposited another series of at least five younger moraine complexes (RC III-VII moraines), previously dated to between 26.4 ± 1.4 ka and 19.9 ± 1.1 ka using TCN exposure dating and Bayesian age modelling (Leger et al., 2021a).The latter moraine belts reflect the local Last Glacial Maximum (LGM) expansions of the PIS.

Methods
Preliminary identification and mapping of major glacial sediment-landform assemblages studied here was carried out using remotely-sensed imagery and topographic data.These geomorphological interpretations were then ground-checked during a 5-week field season during the 2020 austral summer.Details of the geomorphological mapping methodology and criteria, and the complete map of the study area, are given by Leger et al. (2020).To date the RC 0 -II glacier advances, we sampled both outwash surface cobbles and moraine boulders.Specifically, we measured 10 Be in five cobbles from the outermost proglacial outwash surface sampled, here termed the "Tecka outwash", which we correlate stratigraphically with the RC 0 advance (Fig. 1c).We sampled five and six cobbles respectively from the RC I and RC II outwash terraces, six moraine boulders located along a single broad ridge of the RC II moraine complex, and one ice-moulded bedrock surface located between the RC I and II margins .Additionally, we measured 26 Al in five of these 23 samples to determine whether they present simple or complex exposure/burial histories.The chosen moraine boulders are large (90-200 cm high) and rounded to sub-rounded granodiorite and granite erratics embedded in the RC II moraine crest (Fig. 3).Where found, the top 2-5 cm of boulder surfaces exhibiting smooth rock fragments protruding from more eroded surrounding surfaces were sampled using hammer, chisel and angle grinder.
To date the formation of proglacial outwash plains, quartz-bearing and fluvially-rounded cobbles (4-10 cm diameter; crushed whole) were sampled from well-preserved terraces characterised by low surface gradients (~0.5°-1°), following the procedure established by Hein et al. (2009;2011) (Fig. 2a-d).Surface ventifacts and desert varnish, indicative of prolonged cobble exposure, were considered a required sampling criterion.Sample preparation and cosmogenic isotope ratio measurements were conducted at three different laboratories: the cosmogenic isotope analysis and Accelerator Mass Spectrometer (AMS) facilities of the Scottish Universities Environmental Research Centre (SUERC) (East Kilbride, UK), the French National Cosmogenic Nuclides Laboratory (LN2C) and the ASTER AMS facility of the European Centre for Research and Teaching in Environmental Geosciences (CEREGE, Aix-en-Provence, France), and the University of Edinburgh's Cosmogenic Nuclide Laboratory (Edinburgh, UK).Sample details and nuclide concentrations are displayed in Table 1.Details relevant to sample location, photographs, preparation, and wet chemistry are provided in the supplementary materials.
10 Be and 26 Al exposure ages were calculated using the online calculator formerly known as the CRONUS-Earth online calculator version 3 (Balco et al., 2008).Rock density is assumed to be 2.65 g cm -3 and elevation flag is STD.For 10 Be exposure-age calculations, we used the NIST_27900 standardisation (equivalent to 07KNSTD within rounding error) and the central Patagonia production rate (Kaplan et al., 2011) re-calculated from the ICE-D online database (http://calibration.iced.org/).This local production rate yields a current sea-level high-latitude production rate of 3.96 ± 0.24 10 Be atoms g -1 yr -1 using Lm scaling, versus 4.09 ± 0.19 10 Be atoms g -1 yr -1 for the global mean of Borchers et al. (2016).Details regarding specific AMS standards, nominal 10 Be/ 9 Be ratios and 10 Be half-life employed are included in the supplementary materials.
26 Al exposure ages were calculated using the Purdue Z92-0222 AMS standardisation (equivalent to KNSTD within rounding error) and the default CRONUS-Earth online calculator's global production rate of Borchers et al. (2016), as the central Patagonia production rate does not feature 26 Al/ 27 Al measurements.We however consider it reasonable to use both the central Patagonia and global production rates when comparing 10 Be and 26 Al ages, as their respective 10 Be calibrations agree within error, with reported LSDn scaling-fitting parameter values of 0.823 ± 0.067 and 0.846 ± 0.040, respectively (calibration: ICE-D online database).All exposure ages are displayed in Table 2. Below, we discuss exposure ages calculated using the time-dependent LSDn scaling scheme of Lifton et al. (2014) with 1σ analytical uncertainties.Exposure ages are between 6-9% and 0.6-1% older when using the St and Lm scaling schemes of Lal (1991) and Stone (2000), respectively.The 10 Be New Zealand production rate (Putnam et al., 2010), often used for Patagonian TCN exposure-age datasets due to its southern midlatitude location (44°S), decreases our exposure ages by between 1.7% and 1.5% using LSDn scaling, which is less than 1σ analytical uncertainties (3.5-1.6%).

Geomorphology
In the RC valley, the three most distal preserved moraine complexes (RC 0-II) are morphologically distinct from the younger LGM moraines (RC III-VII).They display greater relief (70-90 m) than the RC III -VII moraines (10-40 m; Leger et al., 2021a), are more subdued, and have broader-crested ridges.Along our sampling transect, the three moraine complexes are 14 km (RC 0), 12 km (RC I) and 5 km (RC II) more distal than the outermost-dated LGM moraine (RC III).The sampled RC II moraine is a ~1 km-wide multi-ridge complex presenting a flat-topped, subdued crest surface elevated up to 45 m above its ice-distal outwash plain.Ice-distal slope gradient is variable but does not exceed 5°, a relatively low grade compared to the slopes of the younger, better preserved LGM moraines which feature slope gradients of 9-19° (Leger et al., 2021a).The RC I moraine complex features slope gradients (>5°) and surface geomorphologies for individual moraine ridges comparable to the RC II moraine belt, but is much wider in places (up to 4.6 km) and exhibits a greater concentration of distinct hummocky ridges and moraines.The subdued moraine geomorphologies characteristic of the RC 0-II margins are suggestive of progressive lateral slope downwasting and post-depositional moraine surface erosion of fine material (Clapperton, 1993;Glasser and Jansson, 2008;Hein et al., 2017).Geomorphological observations thus suggest the RC II and older moraine complexes were deposited significantly prior to the local LGM.In the adjacent RH valley, the outermost glacial margin consists of a ~2 km-wide belt of kettle-kame topography displaying a chaotic hummocks-and-hollows configuration, and was formerly termed the "Tecka drift" (Haller et al., 2003;Leger et al., 2020).The sampled Tecka outwash surface marks the distal limit of this kettle-kame deposit (Fig. 1c).Based on geomorphological mapping and moraine-complex traceability across the RC and RH valleys, the Tecka limit (RH valley) is assumed here to correspond stratigraphically to the RC 0 advance, for which no well-preserved outwash plain surface could be found in the RC valley (Fig. 1c).
The outwash surfaces sampled in this study all present long-term preservation properties (Fig. 2).They display low-gradient surfaces (0.5-1°), desert pavements, sparse vegetation cover, and preserved braided palaeo-meltwater channels suggestive of limited surface erosion/deflation (Fig. 2a,d; Leger et al., 2020).These glaciofluvial terraces are composed of sorted and fluvially-rounded sands and gravel deposits.Near the locations of outwash cobble sampling, the outwash surfaces display no detectable signs of post-deposition activation and/or incision by post-glacial hydrology.Where checked, modern soil thickness varied between 0 and 20 cm.This is comparable to observations made further south in the Lago Pueyrredón (47.3°S) and Lago Belgrano (47.5°S) valleys (Hein et al., 2009;Mendelová et al., 2020).At the locations of sampling, these terraces occur adjacent to moraine complexes and mark their distal limits (Fig. 1c), thus enabling the unambiguous stratigraphic association of an outwash deposit with a specific glacial advance.

TCN exposure-age chronology
The exposure-age results for the 23 samples are presented in Table 2. Four of the five cobbles sampled on the outermost Tecka/RC 0 outwash surface (Fig. 1c) range from 252.0 ± 6.3 ka to 284.1 ± 7.0 ka and afford a mean 10 Be exposure age of 268.7 ± 14.4 ka (arithmetic mean ± 1σ standard deviation).One older age (19RHS12: 311.5 ± 8.1 ka) lies outside the remaining population's 95% (2σ) confidence level and is rejected as a statistical outlier.The remaining exposure-age population is relatively well-clustered given the landform age, but still displays a >1 mean square weighted deviation (MSWD) value of 4.38 (>k = 2.63).Such a MSWD value is indicative of greater exposure-age scatter than can be predicted by 1σ analytical uncertainties alone (Jones et al., 2019).
The five 10 Be ages from the RC I outwash surface cobbles range from 185.0 ± 3.0 ka to 257.4 ± 6.6 ka.Although the age population presents no obvious statistical or stratigraphical outliers, the exposure ages are poorly clustered, and present a high MSWD value of 33.6 (> k = 2.41), indicating significantly more exposure-age scatter than solely predicted by analytical uncertainties.From the 26 Al measured in two samples (RC20-13 and RC20-15) from this population, the resulting 26 Al/ 10 Be concentration ratios are 6.4 ± 0.3 and 6.0 ± 0.6, respectively.These ratios are consistent with the canonical 26 Al/ 10 Be surface spallation production ratio of ~6.75, the value currently used in the CRONUS Earth online calculator (Balco et al., 2008), and are here interpreted as indicating a single, continuous exposure history (within uncertainty) post erosion (Granger and Muzikar, 2001;Balco and Rovey, 2008; supplementary materials figure 1).
The ice-moulded bedrock surface (RC20-01), located inboard of the RC I moraine and outboard of the RC II moraine, yields a 10 Be exposure age of 185.7 ± 5.0 ka, an 26 Al exposure age of 181.1 ± 8.5 ka, and an 26 Al/ 10 Be concentration ratio of 6.8 ± 0.3, which is consistent with continuous exposure post erosion of the surface.
The six 10 Be ages from the RC II outwash surface cobbles range from 120.6 ± 2.0 ka to 146.6 ± 4.3 ka with a mean exposure age of 131.3 ± 11.1 ka.The age population is tightly clustered considering analytical uncertainties associated with TCN exposure dating of pre-LGC landforms, and thus features no obvious statistical or stratigraphical outliers.However, the population still yields a MSWD value of 11.4 (> k = 2.26), diagnostic of some non-analytical exposure-age scatter in the 10 Be age distribution.Three of the six 10 Be ages from the RC II moraine boulders range from 125.5 ± 3.4 ka to 161.9 ± 3.6 ka and afford a mean exposure age of 145.3 ± 18.4 ka.Younger ages from the remaining three boulders (RC20-03, 10, 11) plot outside the remaining boulder population's 95% confidence level and the 2σ envelope associated with the youngest RC II outwash cobble (RC20-07).We consider 10 Be inheritance an unlikely source of exposure-age scatter compared to the high potential for boulder exhumation causing young outliers.Moreover, the RC II outwash and RC II moraine belt are geomorphologically likely to represent the same glacier expansion event.These three younger boulder ages, which furthermore would indicate a glaciation occuring during the warmer MIS 5 interglacial, a less probable alternative, were thus rejected as stratigraphical and statistical outliers (Table 2).The three remaining moraine boulder ages still exhibit a high MSWD value of 28.7 (>k = 3.0), indicating that the remaining dataset displays substantially more exposure-age scatter than can be predicted solely by analytical uncertainties.The two oldest sampled boulders, RC20-04 (161.9 ± 3.6 ka) and RC20-05 (148.6 ± 3.0 ka), yield 26 Al exposure ages of 143.6 ± 5.4 ka and 137.5 ± 5.4 ka, corresponding to 26 Al/ 10 Be concentration ratios of 6.3 ± 0.3 and 6.5 ± 0.3, respectively.Such ratios are here interpreted as indicating a continuous exposure post erosion of the sampled boulders.

TCN exposure-age scatter sources
We consider post-depositional shielding of TCN exposure dating samples by vegetation and/or snow to be negligible given the semi-arid and windy conditions (annual 1970-2000 mean wind speed of ~5.3 m s -1 ; WorldClim 2.1 data; Fick and Hijmans., 2017) that characterise the study site and more generally the Argentinian steppe foreland (Garreaud et al., 2013;Hein et al., 2009;2017).Modern daily winter (June-August) precipitation is estimated at 2.2 mm near the sampling site (Fick and Hijmans, 2017), which would represent between 20 and 50 cm thick snowfalls given various snow densities.Snow covers of between 75-150 cm thick persisting for four months of the year are typically required to reduce nuclide production by ~5% (snow density range: 0.16-0.33g cm -3 ; Dunai, 2010).Snow accumulation at the study site is thus considered too low to impact exposure ages beyond analytical uncertainties on 10 5 yr timescales.Cosmogenic nuclide inheritance from previous exposure has been shown to be negligible in other eastern Patagonian valleys (e.g.Douglass et al., 2007;Hein et al., 2009) and is considered unlikely for our samples based on the long (>80 km) distance separating glacial deposits and source regions.Such transport distance should have enabled efficient glacial erosion of transported clasts.For moraine boulder samples, we expect the sources of exposure-age scatter to originate mainly from rock surface erosion and boulder exhumation through moraine surface deflation and lateral creep (Putkonen and Swanson, 2003;Briner et al., 2005;Hein et al., 2017).
For all outwash surface cobbles sampled, total rock-surface erosion is considered negligible due to same reasons as described for moraine boulder samples, but also due to the fluvially-rounded and polished nature of target cobbles.Such interpretation is further supported by the analysis of 10 Be concentrations in a proglacial outwash depth-profile of MIS 8 -old sediments deposited more than 65 km east of the closest bedrock source region, in an eastern Patagonian setting similar to our study site (Hein et al., 2009).Results from this analysis indicate that nuclide inheritance is negligible in outwash deposits of the Río Blanco and Hatcher units, in the Lago Pueyrredón valley (47.5°S).We do, however, expect exposure-age scatter to reflect cobble exhumation via a combination of outwash surface deflation (Hein et al., 2009;Darvill et al., 2015) and near surface turbation (e.g.cryoturbation) caused by possible local development of soil permafrost during cold intervals (Trombotto, 2008).We thus assume post-depositional disturbance to predominantly cause cobble exhumation and young apparent exposure ages (Phillips et al., 1990;Hein et al., 2011), and thus consider the oldest cobble exposure age a better minimum-age estimate of outwash-plain stabilisation following glacier-front retreat, with the exception of obvious statistical outliers (Table 2, Fig. 4).Surface cobbles from the Tecka/RC 0 outwash (Fig. 1c) yield a mean exposure age of 268.7 ± 14.4 ka (n=4), excluding one statistical outlier with an age of 311.5 ± 8.1 ka, and display a discernible amount of exposure-age scatter (MSWD: 4.38 >k).
At the sampling sites, the outwash terrace surface displays preserved braided meltwater channels that suggest minimal outwash surface deflation post deposition.We therefore assess the possible impact of cryoturbation alone on exposure-age underestimation, and estimate broadly the potential magnitude of former cobble exhumation through soil.To do so, we model the constant exhumation of all sampled cobbles (excluding outliers) through a given soil horizon (density: 1.3 g cm -3 ) and calculate the resulting exposure-age bias caused by cosmic-ray attenuation with depth (Gosse and Phillips, 2001; see supplementary materials).We acknowledge that cobble upfreezing is unlikely to have been continuous and of similar magnitude for all cobbles.We assume that some of this variability is considered by reporting a population mean exposure age post-simulation that matches the original age of the oldest cobble, although this assumption undoubtedly yields uncertainties.Using such simulations, we calculated that a soil thickness of ~12 cm would have been required for the Tecka/ RC 0 cobbles to yield a mean 10 Be exposure age (283.1 ± 15.0 ka: +5.7%) that approximately matches that of the oldest original surface cobble.This estimated magnitude of cobble upfreezing is plausible given our observations of modern soil thicknesses (0-20 cm) near to the sampling locations.Therefore, despite local semi-arid conditions, cobble exhumation due to cryoturbation through a realistic soil thickness can potentially explain the observed exposure-age scatter (Fig. 5).Such results may also support the hypothesis that outwash surface erosion and deflation was minimal at this location.
In summary, our interpretations suggest that the exposure ages from the Tecka outwash place the timing of the RC 0 advance

The RC I outwash cobbles and the RC I outermost advance
Surface cobbles from the RC I outwash are poorly clustered (MSWD: 33.6 > k) and afford a mean exposure age of 214.0 ± 29.5 ka and an oldest cobble age of 257.4 ± 6.6 ka (Table 2, Fig. 4).As previously justified (Sect.4.3.),we consider the population's oldest cobble (257.4 ± 6.6 ka) as the closest minimum-age estimate for the RC I outwash deposition.
Simulations of constant cobble upfreezing through soil (Sect.5.1.1.)demonstrate that a mean exposure age (257.2 ± 35.3 ka: + 20%) similar to that of the oldest original cobble is obtained when simulating constant exhumation through a ~39 cm-thick soil horizon.Although reasonable, this average soil thickness is greater than modern observations (0-20 cm).For these cobbles, three distinct outwash surface locations were sampled.They display near identical surface morphologies, vegetation cover and apparent modern soil thicknesses (Fig. 2d).It is thus challenging to assess whether spatially variable cryoturbation could have contributed to the observed scatter.This may indicate that, additionally to cobble upfreezing, outwash surface deflation played, in this instance, an important role in causing significant exhumation of certain RC I cobbles post deposition.This hypothesis is supported by the observation that RC I cobble exposure ages display three distinct age groups that correspond with the three different sampling locations (Table 1, SM figure 3).The first group (sample RC20-16) yields an exposure age of 257.4 ± 6.6 ka, while the second (RC20-12 and 13) and third (RC20-14 and 15) groups display mean exposure ages of 222.4 ± 20.2 ka and 191.5 ± 9.3 ka, respectively.We note that the 1σ standard deviations associated with these age groups do not overlap.Thus, it seems plausible that outwash surface deflation was greater at the locations of the second and third sample groups and, combined with clast upfreezing, could therefore explain the substantial exposure-age scatter observed in the 10 Be age distribution of the RC I outwash cobbles.
An alternative hypothesis to explain the large RC I exposure-age scatter would be to argue that the RC I outwash could be a composite glaciofluvial deposit formed by meltwater runoff associated with numerous distinct RC outlet-glacier advances.
While we commonly expect younger glaciofluvial sedimentation to bury previous deposits, certain cobbles associated with the oldest advances could have been remobilised by either meltwater or ice during younger advances and re-deposited towards the outwash surface.The overlap between the youngest RC I outwash cobble age (185.0 ± 3.0 ka) and the deglacial 10 Be age (185.7 ± 5.0 ka) obtained from the ice-moulded bedrock surface sampled directly inboard of the RC I limit could support this hypothesis.Moreover, the RC I moraine belt is a wider (~4 km) complex constituted of more ice-contact hummocks and hummocky ridges suggestive of a greater glaciogenic sediment volume than the RC 0 and RC II moraine complexes (Fig. 6a).This geomorphological distinction could signal that a higher number of individual outlet-glacier advances have reached the RC I margin.
Establishing which of the previous hypotheses are more likely to be valid would require a more detailed chronological and sedimentological analysis through a depth profile of the RC I outwash deposit.However, regardless of which mechanism is responsible for the observed exposure-age scatter, we argue the oldest exposure age from the RC I outwash suggests a separate extensive advance of the RC outlet glacier occurred at around 270-245 ka, thus during the latest stage of the MIS 8 interval (Figs. 4,6).

The sampled bedrock surface and the innermost RC I advances
The sampled ice-moulded bedrock's position and elevation (1026 m a.s.l.) suggests it was exposed after deglaciation from the RC I moraine-outwash complex, and prior to the RC II advance (Fig. 1c, 6a).Its 10 Be age of 185.7 ± 5.0 ka is thus interpreted as a minimum-age estimate of ice-front retreat following the innermost RC outlet-glacier advance reaching the RC I moraine complex.The older RC I outwash-cobble mean exposure age (214.0 ± 29.5 ka) brings further evidence for this stratigraphic relationship.The 26 Al/ 10 Be ratio (6.8 ± 0.3) from this surface bedrock sample does not indicate a complex exposure-burial history (within uncertainty), thus suggesting the RC glacier did not advance beyond and bury the bedrock for a prolonged and >100 ka period (SM Fig. 1).This ratio also demonstrates efficient subglacial erosion of the bedrock surface, which created a fresh surface free of inherited nuclides when last ice-covered.
Unlike other samples, we noticed signs of surface erosion on the bedrock outcrop (Fig. 2e).The surface demonstrated a lack of fluvio-glacial polish preservation and showed signs of homogeneous granular disintegration of a depth that proved challenging to quantify in the field.The coarser grained and flat nature of this rock surface may have caused it to be more subject to frost wedging than rounded moraine boulders.We thus suspect a certain degree of exposure-age underestimation from this bedrock sample (RC20-01).Applying the rock-surface erosion rate of 0.2 mm ka -1 estimated for semi-arid central Patagonia (Douglass et al., 2007;Hein et al., 2017) increases its age by 3% to 191.3 ± 5.3 ka, which is within the 1σ analytical uncertainty of the exposure age.Subsequently, this minimum deglacial exposure age suggests the youngest outletglacier advance to have reached the RC I moraine complex had to occur prior to ~190 ka.This interpretation is coherent with RC I outwash cobble exposure ages, and together indicate the RC I complex was most likely deposited prior to the MIS 6 (191-130 ka;Lisiecki and Raymo, 2005) cold interval.Because palaeo-climate proxy records such as Antarctic atmospheric temperature and global SST suggest a return to warmer, near-interglacial conditions between the MIS 7d interstadial (~220-230 ka) and ~190 ka (Fig. 8; Parrenin et al., 2013;Shakun et al., 2015), it is possible that the innermost RC I advance occurred during or prior to the MIS 7d interstadial.
To summarise, outwash cobble and bedrock ages suggest that outermost RC I glaciogenic deposits were likely deposited by a MIS 8 glacier advance (~270-245 ka), as portrayed by the oldest outwash cobble age (Figs.6b,7).On the other hand, the innermost and youngest outlet-glacier advance to have reached the multi-ridge, 4 km-wide RC I moraine complex (Fig. 1c, 6a) may have been significantly younger, and perhaps occurred during the MIS 7d interstadial.However, the latter part of this interpretation remains unclear, due to the lack of chronological data for the innermost RC I deposits combined with the analytical uncertainties and the large scatter displayed by outwash cobble exposure ages from this margin.More chronological constraints are thus required to test this hypothesis.

The RC II exposure ages
As for the older RC I and Tecka outwash cobble samples, we argue the observed scatter (MSWD: 11.4 >k) in the RC II outwash cobbles exposure ages (mean: 131.3 ± 11.1 ka) likely originates from cobble exhumation and minor outwashsurface deflation (Fig. 5).This interpretation is reinforced by cobble ages displaying a younger mean age and a tighter cluster than moraine boulders from the stratigraphically-equivalent glacial margin (RC II; Table 2, Fig. 4).Consequently, the oldest surface cobble (146.7 ± 4.3 ka) is deemed a better minimum-age estimate for the RC II glacial advance.We find that a mean exposure age (145.5 ± 12.6 ka: + 11%) similar to that of the oldest cobble can be obtained when simulating constant exhumation of all cobbles (see Sect. 5.1.1.)through a 24 cm-thick soil horizon.This estimate, although highly uncertain, is consistent with modern soil thickness observations near sampling locations (0-20 cm).This simulation suggests that, despite local semi-arid conditions, cobble upfreezing through a realistic soil thickness can cause the observed exposure-age scatter leading to the age underestimation of outwash plain formation.In summary, the RC II outwash cobble ages seem to suggest a minimum-age estimate of ~145 ka.Therefore, we argue the RC II limit reflects an extensive expansion of the local PIS outlet glacier occurring towards the latest stage of MIS 6 (Figs.4,6,7).
Following outlier removal (n = 3), the RC II moraine boulder mean exposure age (n = 3; 145.3 ± 18.4 ka) overlaps the RC II mean outwash cobble age within 1σ analytical uncertainties.The boulder population also displays two ages that are consistent with the set of outwash cobbles (125.5 ± 3.4 ka and 148.6 ± 3.0 ka).These results support our stratigraphic interpretation that these two landforms were deposited by the same outlet-glacier advance.Because we sampled rounded boulders and specifically targeted raised rock fragments presenting polished surfaces, we here consider boulder surface erosion to be negligible, in contrast to the sampled ice-moulded bedrock surface, where granular disintegration was ubiquitous.The population's oldest boulder (RC20-04: 161.9 ± 3.6 ka) displays a ~15 ka older exposure age than the oldest RC II outwash surface cobble.Cosmogenic nuclide inheritance from pre-glacial-transport exposure is considered unlikely given the >80 km distance to source region and the sub-rounded morphology of the boulders sampled indicative of subglacial clast erosion (Fig. 3b,c).We instead propose this older boulder may reflect at least one earlier glacial advance to the same terminal position, or may have been re-deposited during the RC II advance following exposure during previous glacier advances reaching proximal inboard ice-front positions.The 26 Al/ 10 Be concentration ratio from this sample (6.3 ± 0.3) does not suggest a prolonged period of boulder burial (SM Fig. 1).We therefore hypothesise previous PIS advances might have locally occurred during earlier MIS 6 cold intervals (e.g.MIS 6c,e), while the youngest penultimate MIS 6 glacial maximum is more accurately dated by the RC II surface outwash cobbles, as the glaciofluvial terraces likely remained active until the end of the glaciation.
Overall, our results from the RC II limit (n = 12 ages) indicate that outwash surface cobble ages display less postdepositional scatter than moraine boulder ages, with MSWD values of ~11 vs ~30, respectively.This comparison thus agrees with results of Hein et al. (2009Hein et al. ( , 2017) ) suggesting that well-rounded cobbles sampled from the surface of former proglacial outwash plains are better-suited to numerically date pre-LGC glaciations in semi-arid, undisturbed environments such as the Argentinian steppe of Patagonia.However, even if semi-arid conditions usually limit soil thickness to <20 cm on outwash plains of the Argentinian steppe, our simulations show that potential cobble exhumation due to cryoturbation alone can cause meaningful landform-age underestimations for pre-LGC glacial deposits.This source of post-depositional exposure-age scatter is thus important to consider when interpreting outwash surface cobble exposure ages.Another advantage of targeting outwash cobbles is their abundance in formerly active glaciofluvial environments, while moraine boulders on pre-LGC margins are often sparse and where found, are highly eroded.This enables rigorous selection of cobbles that are most suitable to TCN exposure dating.On the other hand, obtaining exposure ages from both moraine boulders and outwash cobbles can provide additional insight on former glacier activity that may not be evident without the combined approach.

The MIS 8 glaciations
Our chronology suggests that two of the most extensive PIS advances preserved at the study site (RC 0 and RC I outermost) occurred during MIS 8, at around 290-270 ka, and 270-245 ka, respectively (Figs. 6,7).These expansion events are coeval, within dating uncertainties, with glacial advances recorded by the two largest outlet glaciers of the former PIS in central Patagonia.There, the penultimate glacial advance in both the Lago Buenos Aires valley (46.5°S;Hein et al., 2017;Fig. 1) and the Lago Pueyrredón valley (47.5°S;Hein et al., 2009), are thought to have occurred during MIS 8 (Fig. 8i,h).
Altogether, these three pre-LGC Patagonian TCN chronologies indicate that the MIS 8 glacial events were some of the PIS's most extensive Pleistocene glaciations.These advances likely affected the former PIS across a large latitudinal range, as they can now be traced across multiple valleys from central to northern Patagonia.The geographical ubiquity of these events is diagnostic of a strong climatic forcing, while their timing coincides with some of the lowest middle Pleistocene δDeuterium values recorded over Antarctica at around 280-270 ka, associated with maximum MIS 8 Antarctic cooling of -8.5°C relative to pre-industrial modern temperatures (Parrenin et al., 2013;Fig. 8).The RC 0 advance may also be coeval with a severe peak in dust flux recorded in the EPICA Dome C Antarctic ice core at ~270 ka (Lambert et al., 2008;Fig. 8).Isotopic composition analyses have demonstrated that Antarctic ice-core dust concentrations are positively correlated to outwashplain activity in southern Patagonia, and thus may directly reflect PIS outlet-glacier expansion events (Sugden et al., 2009; During MIS 8, Northern Hemisphere (NH) mid-latitude summer insolation intensity reached two distinct minima at 280 ka and 255 ka (Berger and Loutre, 1991;Fig. 8).These troughs are synchronous with peaks in the SH mid-latitude seasonality curve indicative of colder winters and warmer summers in Patagonia at those times (Fig. 8d; Darvill et al., 2016).
Importantly, the signal of NH mid-latitude summer insolation intensity also mirrors orbitally-controlled SH seasonal duration (Fig. 8c; Denton et al., 2021).Hence, Patagonia was experiencing longer, colder winters and shorter, warmer summers around 280 ka and 255 ka.In constrast, SH mid-latitude summer insolation intensity reached maximum values at 280 ka and 255 ka, but minimum values at 292 ka and 268 ka, which are within exposure-age scatter and analytical uncertainties associated with dating of the RC 0 and outermost RC I advances (this study), and of the MIS 8 advances of the Lago Buenos Aires and Lago Pueyrredón outlet glaciers also (Hein et al., 2009;2017).Therefore, due to current TCN exposure dating uncertainties, it remains a challenge to determine whether the MIS 8 Patagonian glaciation was coeval with a SH, or a NH mid-latitude summer insolation intensity minimum.However, more accurate dating of younger PIS outletglacier advances occurring during the global LGM (MIS 2; e.g.Peltier et al., 2021), including the RC outlet glacier (RC III-VII; Leger et al., 2021a), has shown that these later SH glacier advances, and also SSTs and Antarctic atmospheric temperature depressions, were anti-phased with SH summer insolation intensity, and were instead coeval with high SH seasonality and winter duration.Consequently, it is plausible that the PIS advanced at ~280 ka (RC 0) and ~255 ka (RC I), when SH summer insolation intensity was high, but when SH winters were cold and long.The proposed timing for the RC 0 advance is consistent with other SH palaeo-climate proxy records which suggest that peak MIS 8 cooling occurred relatively early during the glacial cycle (i.e.280-270 ka; Fig. 8), in contrast with the later MIS 6 and MIS 5d to MIS 2 glaciations, which are characterised by maximum global cooling occurring just prior to sudden glacial terminations (Hughes et al., 2020).
Paradoxically, MIS 8 displays a relatively weak global cooling signal in several records such as benthic δ 18 O (Hughes et al., 2020), and is thought to have been characterised by the lowest sea-level depression relative to other middle Pleistocene 100 ka glacial intervals (-93 m; Spratt and Lisiecki, 2016).Moreover, glacial terrestrial records from different regions demonstrate contrasting ice-sheet and glacier behaviours during MIS 8. Indeed, although stratigraphic evidence for MIS 8 glaciations is often elusive, some has been reported for the Fennoscandian and British-Irish ice sheets (Beets et al., 2005;Davies et al., 2012;White et al., 2010;2017), former ice sheets and glaciers overprinting Poland (Krznanian glaciations; Lindner and Marks;1999), Siberia and Russia (Astakhov et al., 2011;2016), Central and southern Europe (e.g.Preusser et al., 2011), and the West Coast Range of Tasmania (Colhoun and Barrows, 2011;Kiernan et al., 2010).In North America, chronological constraints for stratigraphic evidence of pre-MIS 6 glacial deposits are, in most cases, too imprecise to determine the Marine Isotope Stage to which they relate (Hughes et al., 2020).The review by Hughes et al. (2020) concludes that robust chronological evidence of MIS 8 glaciations is rare, and, where it is found, is often characterised by margins located inboard of MIS 6 limits.In general, direct chronological evidence of MIS 8 glaciations more extensive than MIS 6-2 advances has only been found in Russia (east of the Urals; Astakhov et al., 2016) and Patagonia.Therefore, the extensive MIS 8 glaciation of Patagonia implies a strong regional cooling signal (Hein et al., 2017).
It has been noted in previous work that the major eastern Patagonian valleys, including the RC valley (Figs.6,7), feature records of progressively decreasing glacier extents as the glacier advances get younger (Singer et al., 2004;Kaplan et al., 2009;Cogez et al., 2018).This pattern may not be entirely climate-related.Instead, it might be caused by a combination of long-term erosion of the accumulation area and focused subglacial erosion along former glacial valleys that resulted in glacial overdeepening (Kaplan et al., 2009).During Quaternary glacial cycles, selective subglacial erosion, shown to be most effective along the flowlines of warm-based PIS outlet glaciers and near the centre of the ice sheet where ice was thickest (Clapperton, 1983), has caused rapid and zonally asymmetric overdeepening of main Patagonian valleys (Rabassa and Clapperton, 1990).In the past ~1.2Ma, this erosive process has been estimated to cause >1000 m of bedrock elevation lowering in certain Patagonian valleys (Singer et al., 2004).Kaplan et al. (2009) argued that such overdeepening would likely have overcome the effects of tectonic-and glacial isostatic adjustment-related uplift.For a given location and identical climate forcing, erosion-driven bed lowering may have led to progressively more negative glacier mass balances during each major PIS expansion event, causing less extensive glacier advances with time.Moreover, if more effective towards the former ice divide, subglacial erosion would also have increased retrograde basal slope gradients and thus augmented basal shear stress during subsequent glacier advances.These processes could therefore partly explain the observed pattern of decreasing outlet-glacier extent through time in eastern Patagonia (Fig. 6).Consequently, although the geographically widespread MIS 8 advances of the PIS were most likely the result of strong climate forcing, their outboard positions relative to MIS 6-2 deposits do not necessarily indicate colder regional conditions during MIS 8 than during the younger advances.

The MIS 6 glaciation
TCN exposure ages from the RC II deposit suggest that a major expansion of the PIS occurred during MIS 6, and towards the latter part of the glacial interval, at around 130-150 ka (Figs. 4,6,7).Our chronology suggests that maximum PIS expansion occurred just prior to the rapid global warming leading to MIS 5e, one of the warmer middle-to-late Pleistocene interglacial periods (Fig. 8).This is consistent with global sea-level reconstructions suggesting that the second greatest global ice volume of the past 800 ka, associated with a sea-level minimum of -125 m, was reached towards the end of MIS 6, at ~140 -135 ka (Spratt and Lisiecki, 2016).At that time, mean atmospheric temperatures over Antarctica are estimated to have reached their second lowest values of the past 800 ka (-9.3 °C relative to pre-industrial period; Fig. 8).Global cooling during MIS 6 was more persistent than during other middle Pleistocene glacial cycles.Indeed, this cycle is associated with four major minima in Antarctic atmospheric temperatures spread over ~50 ka, and occurring at around ~185, ~165, ~155 and ~140 ka (Parrenin et al., 2013;Fig. 8).As observed for the MIS 8 interval, these peaks in atmospheric cooling occurred a few ka after NH summer insolation intensity minima, associated with SH seasonality and winter duration maxima.Global SSTs and atmospheric CO 2 concentrations remained consistently low for the same 50 ka time-window, between ~190 and ~140 ka (Bereiter et al., 2015;Shakun et al., 2015).The persistence of the cold phase throughout MIS 6 may have facilitated multiple expansions of the PIS at this time.In our 10 Be chronology, this may be reflected by the potential recycling and redeposition of one RC II moraine boulder sample, yielding a 10 Be exposure ages of 162 ± 4 ka, during the younger penultimate advance dated by the RC II outwash cobbles (~130-150 ka).
Evidence for several distinct MIS 6 glaciations has been reported from other regions of the world.For the European Ice Sheet, for instance, stratigraphic investigations have led to the dating, and naming, of two glaciations within MIS 6, respectively the Late Saalian Drenthe glaciation at 175-155 ka, and the Warthe Stadial between 150 ka and 135 ka (Toucanne et al., 2009;Margari et al., 2014).More generally, MIS 6 glacial deposits are ubiquitous near former glacial margins and are documented by a wider body of evidence than other pre-LGC glaciations (e.g.Rinterknecht et al. 2012;Putnam et al., 2013;Kiernan et al., 2014;Evans et al., 2019;Fernandes et al., 2021), suggesting that the greatest global ice volume of the middle Pleistocene was likely reached during that interval (Hughes et al., 2020).In Patagonia, however, our chronology is amongst the first published datasets offering direct and robust dating of a MIS 6 glacial expansion event.
Indeed, the well-studied Lago Pueyrredón and Lago Buenos Aires pre-LGC moraine-outwash records do not provide clear chronological evidence for MIS 6 advances (Hein et al., 2009;2017).However, this may just imply that the MIS 2 advances of the Lago Pueyrredón and Lago Buenos Aires outlet glaciers were more extensive than any MIS 6 glacial events.Therefore, as for the MIS 8 glaciation, MIS 6 expansions of outlet glaciers may have been ubiquitous through the former PIS.However, the current lack of robust pre-LGC glacier chronologies in Patagonia prevents this hypothesis from being tested adequately.

Synthesis: The timing of Patagonian glaciations
Along with companion publications focusing on the LGC (Leger et al., 2020;2021a), our findings illustrate that the preserved sequence of moraine-outwash complexes in the RC valley system records a minimum of four major PIS expansion events occurring during MIS 8, 6, and 2. In Antarctic atmospheric paleo-temperature proxy records, global SST records, and global atmospheric CO 2 records, the strongest minima between today and 330 ka are all recorded by the RC moraineoutwash complexes, with the exception of the 70-60 ka (MIS 4) cold interval, for which no preserved moraine-outwash complex could be found at the study site.However, we note that a major advance of the PIS near this time did occur on the western side of the Andes at 42°S (García et al., 2021;Gómez et al., 2022).There, Garcia et al. (2021) dated glacial sediments on the Isle of Chiloe to 57.8 ± 4.7 ka, or early MIS 3, with a somewhat smaller advance occurring during MIS 2 at 26.0 ± 2.9 ka.The lack of a corresponding MIS 3/4 moraine east of the mountains suggests that the PIS behaved asynchronously between its western and eastern margins at this time.In southeastern Patagonia, more-extensive-than MIS 2 advances were also dated to MIS 4 and 3 using TCN exposure dating (e.g.Darvill et al., 2016;García et al., 2018;Peltier et al., 2021).We argue that while the RC outlet glacier most likely advanced during the MIS 4 and 3 cooling events, these intervals were relatively short-lived, impeding the PIS from becoming locally thick enough to allow the RC outlet glacier to advance along its retrograde bed slope and generate an extensive advance.Very likely, in addition, all glaciogenic deposits associated with MIS 3/4 advances may have been eroded during the younger MIS 2 advances.This may suggest that steady cooling over periods >15-20 ka is an important pre-requisite for northeastern PIS outlet glaciers to reach a mass balance positive enough to enable highly extensive advances.

Climatic and orbital drivers of SH glaciations -hypotheses
Our chronology implies that the most extensive middle Pleistocene PIS expansion events (MIS 8, 6 and 2) preserved match the strongest minima in Antarctic atmospheric temperatures (Parrenin et al., 2013), global oceanic SSTs (Shakun et al., 2015), and global atmospheric CO 2 concentrations (Bereiter et al., 2015) (Fig. 8).They also appear to occur around the timing of minima in NH summer insolation intensity (60°N) and maxima in SH seasonality, while being out-of-phase with mid-latitude SH summer insolation intensity (Fig. 8c).However, one must note that this statement can only be advanced with confidence for the local MIS 2 expansions of the PIS.For the local MIS 8 and MIS 6 glaciations, this observation is based on current knowledge of 10 Be production rates and the assumptions made in this paper, and does not take into account the full exposure-age range covered by dating analytical uncertainties.The first observations of the asynchrony between SH mid-latitude glacier advances/recessions and local summer insolation intensity (e.g.Mercer, 1976;Denton et al., 1999) encouraged the hypothesis that major glacial-termination-inducing warming signals were primarily controlled by NH summer insolation intensity driving NH glacial mass balance, and were propagated from the northern to the southern hemisphere through oceanic and atmospheric circulation shifts (e.g.Broecker, 1998;Denton et al., 2010).Some of these interhemispheric climate-transmission mechanisms, such as the impact of the abrupt weakening of the Atlantic meridional overturning circulation in the north Atlantic following the last deglaciation, have been well documented (Barker et al., 2009;Denton et al., 2010).However, there is still a great deal of debate regarding the dominant trigger mechanisms responsible for global ice-sheet expansion/retreat cycles (Putnam et al., 2013).Indeed, the above hypothesis entails the dilemma that while NH summer insolation intensity would be a key driver of NH ice-sheet mass balance, it would instead have little impact on mass fluctuations of SH mid-latitude ice sheets and glaciers (Denton et al., 2021).This problem invites the consideration of alternative hypotheses, such as the potential roles of SH seasonality and seasonal duration on controlling global ice-sheet mass fluctuations.García et al. (2018) noted that the extensive advances of the Torres del Paine (51°S) and the Última Esperanza (51.5°S) outlet glaciers at 48 ka (MIS 3) coincide with a minimum in obliquity-modulated high-latitude winter insolation intensity.
The authors therefore hypothesise that colder SH winters during periods of high obliquity may favour SH mid-latitude glaciation by promoting local atmospheric cooling, Southern Ocean stratification and expanded Antarctic sea ice.Since we observe that middle-to-late Pleistocene PIS expansion events appear to occur during periods of longer, colder winters and shorter, warmer summers, we suggest that other PIS chronologies (including this study) may support García et al. (2018)'s hypothesis.
Alternatively, the "Zealandia Switch" hypothesis (Denton et al., 2021)  Ocean (Lorrey and Bostock, 2017).A reduced Aghulhas Leakage has implications for the Atlantic salt input resulting in a potential weakening of Atlantic meridional overturning circulation (Caley et al., 2012;De Boer et al., 2013).The authors thus argue that such a southern oceanic and atmospheric circulation configuration could promote a global glacial mode (Denton et al., 2021) that would greatly affect SH and thus Patagonian glacier mass balance.Moreover, variations in the intensity and latitudinal position of the Southern Westerly Winds influence Southern Ocean upwelling of deep-stored CO 2 and nutrients (Anderson et al., 2009;Ai et al., 2020).This mechanism is thought to represent a major contribution to variations in global atmospheric CO 2 concentrations over millennial timescales (Sarmiento and Toggweiler, 1984).Former equatorward shifts and weakening of the Southern Westerly Winds could have decreased Southern Ocean upwelling and enhanced CO 2 storage in the deep ocean, thus causing delayed global cooling amplification (Sime et al., 2013).Conversely, a poleward migration and strengthening of the Southern Westerly Winds is thought to amplify global warming via increasing Southern Ocean upwelling, CO 2 outgassing and increasing CO 2 atmospheric concentrations.This hypothetical CO 2 capacitor storage-and-release mechanism may explain the highly amplified warming observed during major glacial terminations (Anderson et al., 2009;Denton et al., 2021).This positive feedback mechanism may partly justify why middle-to-late Pleistocene global oceanic and atmospheric temperature variations often reached their most extreme values with some notable delay (a few ka) relative to NH mid-latitude summer insolation intensity.
To summarise, empirical data including this study suggest that major middle-to-late Pleistocene cycles of SH midlatitude glaciers and ice sheets appear to be coeval with NH glaciers and ice sheets, but out-of-phase with local mid-latitude summer insolation intensity (Putnam et al., 2013;Peltier et al., 2021), traditionally thought to be the primary driver of middle-to-late Pleistocene glacial expansion/demise cycles (Mercer, 1984).However, recent hypotheses suggest that this paradox might be explained by the potentially more dominant climate-forcing effect of the seasonality and seasonal-duration expressions of SH insolation, through their impact on atmospheric temperature, the position of the Southern Westerly Winds belt and Subtropical Front, Antarctic sea ice extent, and Southern Ocean stratification and CO 2 storage-and-release feedback mechanisms.Testing the above hypotheses, and determining which of seasonality versus seasonal duration played a primary role in driving SH climate and glacial variations during the middle-to-late Pleistocene, remains a major challenge and represents a key avenue for future research.

Conclusions
• We provide geomorphological and direct chronological evidence for at least three extensive, stratigraphically distinct pre-Last Glacial Cycle advances of the Río Corcovado glacier, a major outlet of the former northern Patagonian Ice Sheet.TCN surface exposure dating of outwash surface cobbles yields best age range estimates for the RC 0, RC I and RC II advances of 290-270 ka, 270-245 ka, and 130-150 ka, respectively.
• Together with evidence we have reported in Leger et al. (2021a, b), our complete chronology reveals a minimum of four distinct expansion events of northeastern Patagonian Ice Sheet outlet glaciers that occurred during the MIS 8, MIS 6, and MIS 2 intervals.
• Our dataset presents the first robust direct dating evidence of MIS 6 glaciations in Patagonia.
• We did not find evidence for MIS 4 and 3 deposits in the RC valley geomorphological record.However, extensive glaciations did occur at these times on the former western margin of the ice sheet, and in other, more southern Patagonian regions.This highlights longitudinal and latitudinal asynchronies in former Patagonian Ice Sheet mass balance during these stages.We suggest that MIS 4 and 3 outlet-glacier advances likely occurred in northeastern Patagonia, but were of relatively shorter extent, and their associated deposits were likely eroded by the younger, more extensive MIS 2 glaciation.
• We interpret that major middle-to-late Pleistocene Patagonian Ice Sheet expansions may have occurred in synchrony with the greatest cooling events recorded in global SST and Antarctic atmospheric temperatures, with the most prominent peaks in Antarctic ice-core dust concentrations, and with major minima in global CO 2 atmospheric concentrations.We observe a particularly good correlation between Antarctic atmospheric temperature signals and southern mid-latitude icesheet volume fluctuations, over several glacial cycles.
• Our results suggest that Patagonian glaciations were more extensive during MIS 8 than during MIS 6-2, while global ice volumes and sea level depression were not particularly high during MIS 8.This implies an interhemispheric asynchrony in the relative volume and extent of Earth's major ice sheets during MIS 8.However, we acknowledge that a larger MIS 8 eastern Patagonian Ice Sheet is not necessarily indicative of a more potent local climate-forcing event than during younger glaciations (MIS 6-2) as glacial erosion, which greatly influences the local topography, may play a role in controlling ice-sheet extent.
• Major middle-to-late Pleistocene advances of the Patagonian Ice Sheet appear to have occurred out-of-phase with local summer insolation intensity but in synchrony with orbitally-controlled periods of longer and colder winters.We hypothesise that seasonality and seasonal duration may exert more control over southern mid-latitude ice-sheet mass balance than mid-latitude summer insolation intensity over millennial timescales.This may be a consequence of the possible effect of seasonal duration variability on the strength and migration of the Southern Westerly Winds, and the impact of these migrations on global climate and Southern Ocean circulation and CO 2 upwelling, as proposed by the "Zealandia Switch" hypothesis.
• Our findings corroborate work by Hein et al. (2009;2017) and indicate that, in eastern Patagonia, pre-LGC outwash surface cobbles yield exposure ages displaying less post-depositional scatter than moraine boulders from the same margins, and hence facilitate more accurate estimations of the approximate timing of pre-LGC glacial advances.

Figure 2 .
Figure 2. Field photographs of: (A-D) proglacial outwash plains and outwash surface cobbles sampled for TCN exposure dating and; (E, F) the ice-moulded bedrock surface sampled (RC20-01).Ages displayed are 10 Be exposure ages ± 1σ analytical uncertainties.Further sample coordinates and characteristics are presented inTable 1.

Figure 3 .
Figure 3. Field photographs of moraine granite and granodiorite boulders sampled for TCN exposure dating.Ages displayed are 10 Be exposure ages ± 1σ analytical uncertainties.
Table footnotes.All cobblesamples were crushed whole, without prior cutting.Sample preparation and wet chemistry was conducted at three different laboratories: ϒ samples: The University of Edinburgh's Cosmogenic Nuclide Laboratory (Edinburgh, UK), * samples: The Scottish Universities Environmental Research Centre (SUERC) (East Kilbride, UK), Δ samples: The French National Cosmogenic Nuclides Laboratory (LN2C) of the European Centre for Research and Teaching in Environmental Geosciences (CEREGE, Aix-en-Provence, France).Nuclide ratios in ϒ and * samples were measured at the SUERC AMS facility (East Kilbride, UK), while nuclide ratios in Δ samples were measured at the ASTER AMS facility (CEREGE, Aix-en-Provence, France).SiO2 is the chemical formula of silica (quartz).# Total 27 Al mass from both carrier ( 27 Al concentration: 982 ± 0.1 µg g -1 ) and sample, determined by ICP-OES.Be Carrier solution used at LN2C has a 9 Be concentration of 3025 ± 9 µg .g - while carrier solutions used at SUERC and the University of Edinburgh's Cosmogenic Nuclide Laboratory have 9 Be concentrations of 849 ± 12 µg .g - and 1000 µg .g - , respectively.Conversions of isotope ratios to 10 Be and 26 Al concentrations were conducted following standard equations, as described by Balco (2006).

Table footnotes. 10
Be and 26 Al exposure ages were calculated using "the online calculator formerly known as the CRONUS-Earth online calculator version 3"(Balco et al., 2008).Rock density is assumed to be 2.65 g .cm - and elevation flag is STD.All samples were collected in 2020. 10Be ages reported here are calculated using the central Patagonia production rate(Kaplan et al., 2011: recalculated from the ICE-D online database: http://calibration.ice-d.org/)while 26 Al ages are here calculated using the global 26 Al production rate of Borchers et al. (2016).AMS standardizations employed for calculator input data are NIST_27900 ( 10 Be) and Purdue Z92-0222 ( 26 Al).Scaling schemes: St is the time-independent version of Lal (1991) and Stone (2000), Lm is the timedependent version of Lal (1991) and Stone (2000), and LSDn is the time-dependent scheme of Lifton et al. (2014).Ages are reported with 1σ analytical and external uncertainties, the latter including production rate and scaling uncertainties in addition to analytical ones.Summary statistics were calculated for each dated landform using only 10 Be exposure ages and LSDn scaling.This includes arithmetic means with 1σ standard deviations (std), 1σ propagated (from individual ages) analytical uncertainties, and propagated 1σ analytical plus production rate uncertainty (PR%).Summary statistics for sets of outwash surface cobble ages also display the oldest cobble exposure age; here considered a better minimum-age estimate of outwash deposit formation.Summary statistics also include uncertainty-weighted means and 1σ standard deviations, MSWDs and Peak ages; calculated using standard equations as described by Jones et al. (2019).Ages in bold are here interpreted as the most appropriate summary ages per landform and are the ones used throughout the paper for discussion (see justification in section 4.3).

Figure 4 .
Figure 4. (A): Kernel density plots, adapted from IceTEA tools for exposure dating outputs (Jones et al., 2019) and summary statistics for dated landforms.Thick red curve represent the summed probability distribution for the exposure-age population, excluding outliers, while thin red curves depict Gaussian curves for individual samples.Outliers are denoted by grey dashed curves.Vertical black and purple dashed lines denote the population's arithmetic mean and peak ages, respectively.(B): Same kernel density plots as (A) but visualised on single time scale to enable better comparison of different landform ages and relative exposure-age scatter between surface outwash cobble and moraine boulder samples.Panel B also features the global LR04 stack of benthic δ 18 O records of Lisiecki and Raymo (2005); commonly used to define the timing of Marine Isotope Stages, here also denoted by alternating white and blue vertical bands.

Figure 5 .
Figure 5. Illustration of the impact of combined outwash-plain surface deflation and near-surface turbation on exposure-age scatter for outwash surface cobbles sampled.The relationship between depth and nuclide concentration is illustrated by the decay curve (purely schematic) on the bottom left-hand corner of the diagram, which demonstrates that radionuclide production in rock minerals is reduced with depth below surface due to cosmic-ray attenuation.The attenuation effect causes cobbles that have experienced some exhumation post-deposition, due to the effect of upfreezing (dashed vertical arrows; scenario 1), to display exposure-age scatter causing landform deposition age underestimation.Cobbles that remained near the outwash surface since deposition (scenario 2) will display less/no scatter and better estimate the true minimum outwash deposition age.Cobble upfreezing is considered to mainly occur within the soil horizon (brown-coloured) due to the impact of potential permafrost formation.Soil thickness is locally estimated to reach a maximum of ~20 cm.The outwash surface photograph focuses on one of

Figure 6 .
Figure 6.(A): Glacial geomorphological map focused on the RC 0 -RC II moraine-outwash sequence around the location of TCN exposure dating samples, highlighted by white dots (with the exception of Tecka/RC 0 outwash samples (Fig. 1c)).Exposure ages displayed are 10 Be ages ± 1σ analytical uncertainties, calibrated with LSDn scaling and central Patagonia production rate (Kaplan et al., 2011).Ages in red were interpreted as outliers."Bed" stands for "bedrock" and depicts the ice moulded bedrock surface sampled.Elevation data is provided by the 30-m resolution DEM from the ALOS World 3D missions (AW3D30, version 2.2; JAXA; https://www.Eorc.jaxa.jp/ALOS/en/aw3d30/)with a shaded relief background and elevation contour lines at 15-m intervals.(B): Simplified and smoothed elevation profile graph (DEM data: AW3D30) drawn along the black line denoted in panel (A), with TCN exposure dating results ( 10 Be) from the RC 0 -III moraine-outwash dataset after interpretation of the best minimum-age estimate for each dated landform.Palaeo-ice surface elevation and slope gradients are hypothetical and purely schematic.*: Stratigraphic order of events (1: oldest, 5: youngest).

Figure 7 .
Figure 7. Palaeo-glacial reconstructions for the RC, RH, RF and Lago Palena/General Vintter (LP/LGV) valleys, for each advances/stillstands and deglacial events interpreted in this investigation (A-C).Panel D displays the ice extent associated with the RC III moraine complex, described in detail in companion publication (Leger et al., 2021a).The latter reconstruction includes the formation of glaciolacustrine phase one (here termed PL1; Leger et al., 2021a).Ice-sheet and mountain glacier models were digitised manually in ArcGIS.Apart from when delineated by confidently-mapped moraine limits (e.g. the RC, RH and LP/LGV moraine sequences), the geographical location of ice margins are inferred from topography.Topographic data is from the AW3D30 DEM."OC" and "MB" denote the type of rock samples used for TCN exposure dating of the represented event, and stand for "outwash cobbles", and "moraine boulders", respectively.Yellow stars indicate the approximate sample locations for TCN exposure dating of each represented event, while each panel features the number of 10 Be exposure ages produced (n = x) on the bottom left-hand corner (this excludes 26 Al ages).There are no previously published glacier chronologies in the RH, RF and LP/LGV valleys.Relative ice extents in those neighbouring valleys are thus inferred from our RC chronology and cross-valley comparisons of moraine set numbers, preservation and morphostratigraphy.Hence, these inferences yield some uncertainties.Proglacial lake volume estimates were computed from DEM data (AW3D30).

Figure 8 .
Figure 8. Vertical plot comparing this investigation's TCN 10 Be chronology to other Patagonian outlet-glacier TCN chronologies, to insolation intensity variations, and to palaeo-climate proxy records over the past 350 ka.(A) Stable isotope δ Deuterium record (orange curve) from the EPICA Dome C (EDC) ice core, 75°S, east-Antarctica and Antarctic atmospheric temperature (overlying dashed black curve) from a stack of isotopic temperature reconstructions from five different ice cores (EDC, Vostok, Dome Fuji, TALDICE, and EDML) (Parrenin et al., 2013), with most potent cooling events of the past 350 ka highlighted by black circles.TCN exposure ages from the RC valley interpreted as best minimum-age estimate for each reconstructed glacial event are displayed (red dots) with analytical (black horizontal bars) and external uncertainties (red horizontal bars).(B): Global stack (n = 49 cores) of paired planktonic δ 18 O-derived Sea Surface Temperature (SST) data with Marine Isotope Stage numbers indicated (Shakun et al., 2015).(C): Dark green line: summer (mid-June) insolation intensity for 60°N (Berger and Loutre, 1991).Dashed black curve: seasonal duration at 60°S (insolation threshold = 300 W m -2 ) after Huybers and Denton (2008) and Denton et al. (2021).(D): Southern Hemisphere seasonality at 30°C, 40°S, and 60°S calculated by subtracting June from December insolation intensity values (Berger and Loutre, 1991), such that increasing seasonality indicates colder winters and warmer summers (after Darvill et al., 2016).(E): Atmospheric CO2 concentrations composite record from several Antarctic ice cores including EDC, WAIS, Vostock, Siple Dome, TALDICE, EDML, and LAW Dome (Bereiter et al., 2015).(F): Dust flux record from the EDC ice core (Lambert et al., 2008).(G-I): Kernel density distribution curves for TCN exposure ages from the RC O -RC II moraineoutwash records (this study: G), the Hatcher moraine-outwash complex, in the Lago Cochrane/Pueyrredón valley (Hein et al., 2009; H), and the Moreno 2/3 moraine-outwash limits, in the Lago General Carerra / Buenos Aires valley (Hein et al., 2017; I).
focuses on the SH seasonal-duration aspect of orbitally-controlled insolation as a potential driver of late-Quaternary glaciation.By pacing the seasonal cycle of weakened and northward-migrated Southern Westerly Winds during winter versus strengthened and southward-migrated Southern Westerly Winds during summer, this hypothesis implies that SH seasonal duration might be a key driver of global climate change during middle-to-late Pleistocene glacial/interglacial cycles.Using global climate model simulations, the authors show that a reduced SH summer duration and increased SH winter duration could cause millennial-scale equatorward migrations of the Southern Westerly Winds and its coupled Antarctic Circumpolar Current.Such zonal wind configuration would induce a northward migration of the subtropical front and a latitudinally-restricted subtropical gyre system in the Indo-Pacific Ocean basins (ibid).The bathymetry of the Australia/Zealandia continent is such that a restricted Indo-Pacific subtropical gyre could decrease Indonesian Throughflow and restrict westward Agulhas Leakage into the South Atlantic