U–Th and ¹⁰Be constraints on sediment recycling in proglacial settings, Lago Buenos Aires, Patagonia

The estimation of sediment transfer times remains a challenge to our understanding of sediment budgets and the relationships between erosion and climate. Uranium (U) and thorium (Th) isotope disequilibria offer a means of more robustly constraining sediment transfer times. Here, we present new uranium and thorium disequilibrium data for a series of nested moraines around Lago Buenos Aires in Argentine Patagonia. The glacial chronology for the area is constrained using in situ cosmogenic ¹⁰Be analysis of glacial outwash. Sediment transfer times within the periglacial domain were estimated by comparing the deposition ages of moraines to the theoretical age of sediment production, i.e., the comminution age inferred from U disequilibrium data and recoil loss factor estimates. Our data show first that the classical comminution age approach must include weathering processes accounted for by measuring Th disequilibrium. Second, our combined data suggest that the pre-deposition history of the moraine sediments is not negligible, as evidenced by the large disequilibrium of the youngest moraines despite the equilibrium of the corresponding glacial flour. Monte Carlo simulations suggest that weathering was more intense before the deposition of the moraines and that the transfer time of the fine sediments to the moraines was on the order of 100–200 kyr. Long transfer times could result from a combination of long sediment residence times in the proglacial lake (recurrence time of a glacial cycle) and the remobilization of sediments from moraines deposited during previous glacial cycles. ¹⁰Be data suggest that some glacial cycles are absent from the preserved moraine record (seemingly every second cycle), supporting a model of reworking moraines and/or fluctuations in the extent of glacial advances. The chronological pattern is consistent with the U–Th disequilibrium data and the 100–200 kyr transfer time. This long transfer time raises the question of the proportion of freshly eroded sediments that escape (or not) the proglacial environments during glacial periods.

1 Introduction

The sedimentary cycle incorporates the erosion of rocks followed by the transport and deposition of sediments. While rates of erosion and deposition can be accurately documented, tracing the history of sediments between production and deposition remains challenging. Importantly, mechanisms of transfer and the alteration of sediments during transport play a key role in the evolution of basins and feedbacks between erosion and climate, especially because the age of sediment strongly controls its susceptibility to weathering (Vance et al., 2009; White and Brantley, 2003). This is particularly the case in glacial settings because glaciers are highly efficient at eroding landscapes (Hallet et al., 1996; Koppes and Montgomery, 2009) and produce highly reactive and easily weathered sediments (Anderson, 2005; Anderson et al., 1997; White and Brantley, 2003). Moreover, glaciers can create large overdeepenings that are subsequently filled with sediments isolated from interactions with surface processes. The efficiency of evacuating those sediments post-deposition is poorly known.

Silicate weathering is an important surface parameter for cooling climate over geologic and glacial–interglacial timescales because it consumes CO₂ (Ebelmen, 1845; Walker et al., 1981). But its role in controlling CO₂ concentrations and climate variations over glacial–interglacial cycles is contentious (Cogez et al., 2015; Foster and Vance, 2006; Lupker et al., 2013; Vance et al., 2009; von Blanckenburg et al., 2015). Understanding whether weathering varied over these cycles requires robust determinations of sediment transfer times. For example, long transport times can cause a lag and/or a damping of the response of weathering to climate or erosion forcing and bias reconstructions of erosion and weathering intensity and variations through time.

Geochemical tools offer a means of measuring sediment transport times. Uranium series isotopes are particularly useful because the diversity of chemical elements in the radioactive decay chain (U, Th and Ra) leads to disequilibrium during surface process fractionation. Moreover, the half-lives of these isotopes range from 1500 to 250 000 years, corresponding to sediment transport process times. The timescales of weathering processes have been particularly well documented using U–Th–Ra disequilibria in soils, sediments and river waters (Ackerer et al., 2016; Chabaux et al., 2008, 2012, 2013; Dosseto and Schaller, 2016; Dosseto et al., 2010, 2012, 2014; Granet et al., 2010; Keech et al., 2013; Ma et al., 2013). An alternative approach uses the fine fraction of silicates to date the time since the physical erosion of sediments (DePaolo et al., 2006, 2012). The method is based on the α recoil of uranium during radioactive decay, triggering the loss of a fraction of the daughter isotopes compared to the parent in small (≤50 µm) grains. The time since comminution can theoretically be estimated, despite difficulties and limitations discussed in previous studies (Dosseto et al., 2010; Handley et al., 2013a, b; Lee et al., 2010; Maher et al., 2006). Here we use this method and also take into account weathering, as this was shown to be an important driver of U-series disequilibria in a boreal environment (Andersen et al., 2013).

In this study, we sampled sequences of nested moraines around Lago Buenos Aires in Argentine Patagonia, which range in age from 0 to 1 Ma (Kaplan et al., 2005; Singer et al., 2004). The aim was to constrain the pre-deposition history of sediments within the moraines. First we refined the deposition age of five of the moraines using in situ cosmogenic ¹⁰Be exposure and depth profile dating. To constrain the sediment transport times, we used the comminution age approach combined with ²³⁰Th disequilibrium measurements. We show that taking weathering into account can help to resolve previous issues with the comminution age methodology. Finally, we are able to estimate that the pre-deposition history of the fine sediment (4–50 µm) is likely on the order of 100–200 kyr, including sediment recycling and chemical weathering in the proglacial system. Our findings have implications for the understanding of the couplings between erosion and climate.

2 Settings and methods

2.1 The Lago Buenos Aires moraines

Lago Buenos Aires in Argentina (Lago General Carreras in Chile) is a large proglacial lake, around 100 km long and 5–20 km wide; it is oriented east–west at 46^∘ S (Fig. 1). The present climate of the area is temperate, with an average annual precipitation of 100 mm yr⁻¹ and temperatures between 4 and 14 ^∘C, producing open steppe vegetation. On the eastern edge of the lake, a series of frontal moraines are nested from the youngest close to the lake to the oldest 50 km further east. The chronology of these moraines has been studied: the five innermost moraines (Fenix 1–5) are Last Glacial Maximum (LGM) in age (Douglass et al., 2006; Kaplan et al., 2004). Singer et al. (2004) dated lava flows interbedded with six of the intermediary moraines (Moreno 1–3 and Deseado 1–3) and showed that they range in age between 109 and 760 ka. Finally, six moraines (Telken 1–6) in the outermost part of the system were likely deposited between 760 and 1016 ka. Mercer and Sutter (1982) showed that the oldest glaciogenic sediments we find in the area are 6 Ma and consist of tills interbedded with lava flows that preserved them from erosion. Cosmogenic nuclide exposure dating of boulders supports the assertion that erosion and degradation of older moraines in this region yields erroneously young ages (Hein et al., 2009; Kaplan et al., 2005). Recently, Hein et al. (2017) used exposure dating of cobbles on outwash related to the Moreno moraines to show that they were constructed at ca. 260–270 ka during marine isotope stage (MIS) 8. For this study, we further refined the chronology by providing direct age constraints for the Deseado 1, Deseado 2 and Telken 5 moraines (see Sect. 2.2).

Figure 1 Map of the study area adapted from Bendle et al. (2017). (a) Southern South America showing the location of Lago Buenos Aires. The three main contemporary ice masses in Patagonia are also shown: the Northern (NPI) and Southern Patagonian Ice Field (SPI) and the Cordillera Darwin Ice Field (CDI). (b) Inset shows a west–east transect through the study area from the Patagonian Andes to the frontal moraine systems of the proglacial lake (adapted from Kaplan et al., 2009); sample locations as per (a). Panels (c) and (d) show the series of frontal moraines sampled in this study. Silt samples are shown by red circles and cobbles and profiles for cosmogenic analysis are shown in blue; the glacial flour sample at the mouth of Los Exploradores glacier is shown, and the extent of the present day ice field is shown in black. Successive glaciations have formed the large overdeepening in which the proglacial lake has formed.

U and Th samples were taken from silty beds inside the moraines several meters below the surface in order to avoid potential post-depositional modifications, such as weathering or dust inputs (see picture in the Supplement). In addition to moraine samples, a glacial flour sample was taken at the front of the Los Exploradores glacier in Laguna San Rafael National Park to estimate the initial U–Th composition of silts found in the moraines (Fig. 1). This glacier covers part of the total catchment draining into Lago Buenos Aires. Consequently, it is unlikely to be representative of the entire basin but provides an estimate of the initial composition of silts after comminution.

2.2 ¹⁰Be methodology

2.2.1 Sampling

The cosmogenic nuclide exposure dating methodology followed the sampling recommendations of Hein et al. (2009) and Darvill et al. (2015). We sampled quartzite cobbles for ¹⁰Be analysis from the surfaces of outwash grading to the Deseado 1, Deseado 2 and Telken 5 moraines. Similar samples were taken for Moreno 1 and 3 moraines and are presented in Hein et al. (2009). Using outwash cobbles has been shown to overcome issues of post-depositional erosion that affect moraine boulders in this area (Hein et al., 2009), but relies on finding locations where there is a clear geomorphic relationship between the glacial outwash and moraines. At suitable locations, cobbles were targeted based on quartz composition, size (6–15 cm and 300–1 000 g to contain sufficient ¹⁰Be) and preservation (sub-rounded shape, little aeolian erosion and partially buried, suggesting they were in situ). Surface lowering or up-freezing of cobbles can result in vertical movement through the outwash unit. In this region, it is more likely that cobbles will underestimate the true age of the unit (Hein et al., 2009, 2017). To circumvent this issue three or four cobbles were analyzed at each sampling location, and the oldest was taken as the best estimate of moraine deposition age.

In addition to cobbles from the outwash surfaces, ¹⁰Be depth profiles through exposures within the outwash units were used to estimate surface erosion rates and nuclide inheritance. We sampled through two outwash profiles relating to the Deseado 2 and Telken 5 moraines. Each of the samples making up the depth profiles consisted of 50 to 100 quartz-rich gravel clasts taken from horizontal lines 2–4 cm wide at various depths below the surface. Five depths were sampled per profile between 50 and 240 cm, which is sufficient depth to account for nuclide attenuation from the surface and reliably estimate any cosmogenic ¹⁰Be inheritance. The top 50 cm of the profiles were carbonated due to post-depositional processes but the rest of the profiles were of relatively consistent density. See the Supplement for pictures.

2.2.2 Analytical methods

Surface cobbles were analyzed individually as independent estimates of exposure age. The gravels at each depth in the profiles were amalgamated to produce an average nuclide concentration for that depth. Samples were crushed and sieved to obtain the 125–250 µm and the 250–500 µm fractions, which were then purified using a Frantz magnetic separator to isolate non-magnetic fractions. Following dissolution in concentrated HF, the residual sample material was leached with 10 mL H₂O to extract Be(OH)₂ (Stone, 1998). A combination of anion and cation resins and precipitation was used to purify Be (Norton et al., 2008; von Blanckenburg et al., 2004). The oxidized BeO was mixed with Nb powder, pressed into cathodes and measured on the 0.5 MeV Tandy accelerator at ETH Zürich (Müller et al., 2010). The resulting ¹⁰Be∕⁹Be ratios were normalized to S2007N (Christl et al., 2013; Kubik and Christl, 2010), with a ¹⁰Be half-life of 1.387±0.012 Ma (Chmeleff et al., 2010; Korschinek et al., 2010). Measured ¹⁰Be∕⁹Be ratios were 30–500 times higher than the procedure blank ($\mathrm{9.60}\times {\mathrm{10}}^{-\mathrm{15}}\pm \mathrm{4.77}\times {\mathrm{10}}^{-\mathrm{15}}$). Blank-corrected ¹⁰Be concentrations decrease exponentially with depth from greater than 36×10⁵ at g⁻¹ for near-surface samples to less than 2×10⁵ at g⁻¹ at 2 m of depth.

2.2.3 Deposition age calculations

Exposure ages for the outwash surface cobbles were calculated from ¹⁰Be concentrations using version 2.3 of the CRONUS online calculator (Balco et al., 2008). The shielding factor was measured in the field (>0.999999 for every sample), a density of 2.65 g cm⁻³ was used for every sample and no erosion rate was applied since there were no visible signs of surface erosion. We used the global ¹⁰Be production rate of Borchers et al. (2016) (4 $\mathrm{at}{\mathrm{g}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$) and the “Lm” scaling scheme (as described in Balco et al., 2008, Lal, 1991, and Stone, 2000, with paleomagnetic corrections from Nishiizumi et al., 1989). However, the scatter between the different scaling schemes is comparable to the external uncertainties. The production rate from Kaplan et al. (2011), calculated from a site very close to LBA, gives 3.6 $\mathrm{at}{\mathrm{g}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ (Borchers et al., 2016). The difference in calculated exposure ages between these two production rates remains within uncertainties.

We modeled the most probable deposition age, surface erosion rate and nuclide inheritance for each depth profile using the Monte Carlo simulation approach of Hidy et al. (2010). A site production rate of 6.2 $\mathrm{at}{\mathrm{g}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ was used, taken from an average of the scaling schemes in the CRONUS calculator. Changing this production rate by ±0.3 $\mathrm{at}{\mathrm{g}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ does not significantly alter the results. This range of site production rates (5.9–6.5 $\mathrm{at}{\mathrm{g}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$) corresponds to SLHL production rates of 3.6–4 $\mathrm{at}{\mathrm{g}}^{-\mathrm{1}}{\mathrm{yr}}^{-\mathrm{1}}$ using the Lm scaling scheme. We used an average profile density of 2.5 g cm⁻³ (based on field observations) and a shielding factor of 1. The a priori values assigned to the Bayesian Monte Carlo simulation were the following: 0 to 1 cm ka⁻¹ erosion rate with a maximum erosion threshold of 1000 cm and uniform distribution (based on erosion rates in Hein et al., 2009, for the Hatcher profile in nearby Lago Pueyrredón and Darvill et al., 2015); 0 to 10⁵ at g⁻¹ inheritance with uniform distribution (based on concentrations in the lower levels of the profiles); and 200 to 800 ka age range for Deseado 2 and 700 to 1200 ka age range for Telken 5 with uniform distribution. These ages were chosen based on the argon ages of Singer et al. (2004) and the Moreno age refinement of Hein et al. (2017).

2.3 U–Th methodology

2.3.1 Principles and theory

U–Th-series disequilibria have been used for a few decades in order to determine the time constant of weathering. Refer to the Supplement for more references. The comminution age theory introduced by DePaolo et al. (2006) is a promising approach to estimating the age of sediment production. ²³⁸U decays into ²³⁴Th with a half-life of 4.5 Ga. During this α emission, the daughter isotope recoils over a distance of several tens of nanometers depending on the type of mineral (e.g., 30 nm for a feldspar). ²³⁴Th quasi-instantaneously produces ²³⁴U, with a longer half-life of 245 kyr, which is long enough to be measurable and usable on geologic timescales. If the ²³⁸U decay occurs on the edge of a sediment grain, the daughter isotope is ejected from the grain and a disequilibrium in the decay chain starts as the ²³⁴U∕²³⁸U ratio of the grain starts to decrease. If the grain is small enough with a large surface area compared to its volume, the disequilibrium can be measured. Steady state is reached after a time characterized by the half-life of ²³⁴U. During this transient state, the magnitude of the disequilibrium in a grain equates to the time since it was produced. The steady-state disequilibrium value depends on the fraction of ²³⁴U ejected out of the grain due to α recoil. The process can be modeled using the following two equations describing the evolution of ²³⁸U and ²³⁴U:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle \frac{\mathrm{d}{N}_{\mathrm{238}}}{\mathrm{d}t}}=-{\mathit{\lambda }}_{\mathrm{238}}{N}_{\mathrm{238}}\\ {\displaystyle }& {\displaystyle \frac{\mathrm{d}{N}_{\mathrm{234}}}{\mathrm{d}t}}=-{\mathit{\lambda }}_{\mathrm{234}}{N}_{\mathrm{234}}+(\mathrm{1}-f_{\mathit{\alpha }}^{\mathrm{234}}){\mathit{\lambda }}_{\mathrm{238}}{N}_{\mathrm{238}}\end{array}$$

λ₂₃₄ and λ₂₃₈ are the decay constants of ²³⁴U and ²³⁸U, respectively. N₂₃₄ and N₂₃₈ are the number of nuclides for each isotope, and $f_{\mathit{\alpha }}^{\mathrm{234}}$ is the recoil loss fraction of ²³⁴U (the proportion of ²³⁴U in the grain ejected during α decays). This parameter is critical, since it represents the process that determines the possibility of tracing the comminution age. When the effect of this process starts to be measurable, the clock starts. We discuss in more detail how this parameter can be evaluated in Sect. 2.3.4. The evolution of the U∕Th ratio over time is shown in Fig. 3 (blue curves). The steady-state disequilibrium model is based on the assumption that the sediment grain has not been weathered over time, or at least that the loss of nuclide due to weathering was negligible compared to loss by α recoil. The assumption may be valid for very arid areas, but is otherwise unlikely to hold true. However, the influence of weathering on the comminution age model has not been investigated.

In order to decouple the respective influence of α recoil and weathering, we measured ²³⁰Th concentration. The solubility of thorium is much lower than uranium. Therefore the ratio of ²³⁰Th∕²³⁴U can inform weathering rates. ²³⁰Th is produced by ²³⁴U decay with a half-life of 75 kyr. It can also be ejected out of the grain due to α recoil. The recoil loss fraction of ²³⁰Th, $f_{\mathit{\alpha }}^{\mathrm{230}}$, is related to $f_{\mathit{\alpha }}^{\mathrm{234}}$ (Hashimoto et al., 1985; Neymark, 2011).

$$\begin{array}{}\text{(1)}& {\displaystyle }{\displaystyle }{f}_{\mathit{\alpha }}^{\mathrm{230}}=f_{\mathit{\alpha }}^{\mathrm{234}}\times {\displaystyle \frac{{\mathit{\alpha }}_{\mathrm{230}}{M}_{\mathrm{234}}}{{\mathit{\alpha }}_{\mathrm{234}}{M}_{\mathrm{230}}}}=\mathrm{1.176}\times f_{\mathit{\alpha }}^{\mathrm{234}}\end{array}$$

α₂₃₀ and α₂₃₄ are the decay energy of ²³⁸U to ²³⁴U and of ²³⁴U to ²³⁰Th, respectively (71.79 and 82.96 keV, respectively). M₂₃₄ and M₂₃₀ are the atomic mass of ²³⁴U and ²³⁰Th. This gives $f_{\mathit{\alpha }}^{\mathrm{230}}=\mathrm{1.176}\times f_{\mathit{\alpha }}^{\mathrm{234}}$ (in practice $f_{\mathit{\alpha }}^{\mathrm{230}}$ is probably larger since ²³⁴U is on damaged sites, so it increases the probability for ²³⁰Th to be ejected from the grain during a decay). The evolution of the number of nuclides ²³⁰Th is then given by

$$\begin{array}{}\text{(2)}& {\displaystyle }{\displaystyle \frac{\mathrm{d}{N}_{\mathrm{230}}}{\mathrm{d}t}}=-{\mathit{\lambda }}_{\mathrm{230}}{N}_{\mathrm{230}}+(\mathrm{1}-f_{\mathit{\alpha }}^{\mathrm{230}}){\mathit{\lambda }}_{\mathrm{234}}{N}_{\mathrm{234}}.\end{array}$$

Since $f_{\mathit{\alpha }}^{\mathrm{230}}>f_{\mathit{\alpha }}^{\mathrm{234}}$ and λ₂₃₀>λ₂₃₄, $({}^{\mathrm{230}}\mathrm{Th}/{}^{\mathrm{238}}\mathrm{U})<({}^{\mathrm{230}}\mathrm{Th}/{}^{\mathrm{234}}\mathrm{U})<({}^{\mathrm{234}}\mathrm{U}/{}^{\mathrm{238}}\mathrm{U})$ is always verified. However, weathering would increase the ratio of ²³⁰Th∕²³⁴U and ²³⁰Th∕²³⁸U because Th is less soluble than U, and likely decrease the ratio of ²³⁴U∕²³⁸U because radioactive decay places ²³⁴U on a more fragile site compared to ²³⁸U so that it could be more easily weathered (DePaolo et al., 2012; Handley et al., 2013a). Consequently, Th isotopes can help to identify and quantify the effects of α recoil and weathering, respectively. New equations can be derived for the evolution of the different nuclide contents following Chabaux et al. (2003, 2008), Dosseto et al. (2008) and Dosseto and Schaller (2016):

$$\begin{array}{ll}{\displaystyle }& {\displaystyle \frac{\mathrm{d}{N}_{\mathrm{238}}}{\mathrm{d}t}}=-k_{\mathrm{238}}{N}_{\mathrm{238}}-{\mathit{\lambda }}_{\mathrm{238}}{N}_{\mathrm{238}}\\ {\displaystyle }& {\displaystyle \frac{\mathrm{d}{N}_{\mathrm{234}}}{\mathrm{d}t}}=-k_{\mathrm{234}}{N}_{\mathrm{234}}-{\mathit{\lambda }}_{\mathrm{234}}{N}_{\mathrm{234}}+(\mathrm{1}-f_{\mathit{\alpha }}^{\mathrm{234}}){\mathit{\lambda }}_{\mathrm{238}}{N}_{\mathrm{238}}\\ \text{(3)}& {\displaystyle }& {\displaystyle \frac{\mathrm{d}{N}_{\mathrm{230}}}{\mathrm{d}t}}=-k_{\mathrm{230}}{N}_{\mathrm{230}}-{\mathit{\lambda }}_{\mathrm{230}}{N}_{\mathrm{230}}+(\mathrm{1}-f_{\mathit{\alpha }}^{\mathrm{230}}){\mathit{\lambda }}_{\mathrm{234}}{N}_{\mathrm{234}},\end{array}$$

where k₂₃₈, k₂₃₄ and k₂₃₀ are the leaching coefficients of ²³⁸U, ²³⁴U and ²³⁰Th, respectively. The solution is shown in Fig. 3.

The steady-state activity ratios can also be derived.

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\left({}^{\mathrm{234}}\mathrm{U}/{}^{\mathrm{238}}\mathrm{U}\right)}_{\text{steady}}={\displaystyle \frac{{\mathit{\lambda }}_{\mathrm{234}}\left(\mathrm{1}-f_{\mathit{\alpha }}^{\mathrm{234}}\right)}{{\mathit{\lambda }}_{\mathrm{234}}+k_{\mathrm{234}}-k_{\mathrm{238}}}}\\ {\displaystyle }& {\displaystyle }{\left({}^{\mathrm{230}}\mathrm{Th}/{}^{\mathrm{234}}\mathrm{U}\right)}_{\text{steady}}={\displaystyle \frac{{\mathit{\lambda }}_{\mathrm{230}}\left(\mathrm{1}-f_{\mathit{\alpha }}^{\mathrm{230}}\right)}{{\mathit{\lambda }}_{\mathrm{230}}+k_{\mathrm{230}}-k_{\mathrm{238}}}}\\ {\displaystyle }& {\displaystyle }{\left({}^{\mathrm{230}}\mathrm{Th}/{}^{\mathrm{238}}\mathrm{U}\right)}_{\text{steady}}\\ {\displaystyle }& {\displaystyle }={\displaystyle \frac{{\mathit{\lambda }}_{\mathrm{230}}{\mathit{\lambda }}_{\mathrm{234}}\left(\mathrm{1}-f_{\mathit{\alpha }}^{\mathrm{230}}\right)\left(\mathrm{1}-f_{\mathit{\alpha }}^{\mathrm{234}}\right)}{\left({\mathit{\lambda }}_{\mathrm{230}}+k_{\mathrm{230}}-k_{\mathrm{238}}\right)\left({\mathit{\lambda }}_{\mathrm{234}}+k_{\mathrm{234}}-k_{\mathrm{238}}\right)}}\end{array}$$

Here, parentheses are used to represent activity ratios corresponding to the isotopic ratio multiplied by the ratio of the respective decay constants of each isotope.

2.3.2 Sample preparation

Silt samples were collected from moraines. Water has likely percolated through these silts since moraine deposition, triggering the precipitation of secondary carbonates, oxides, organic matter and clays. These secondary fractions have different U and Th disequilibria than the primary silicate and it is therefore important to remove them before analysis without altering the primary silicate fraction. We used the protocol established by Lee (2009), with slight modifications to help overcome the low solubility of Th that can make it difficult to remove during successive leaching, and allow Th to adsorb onto mineral surfaces after ejection from the grain. Our protocol is as follows.

1. The samples were sieved at 50 µm in order to collect the ≤50 µm fraction.

2. Approximately 10 g of this fraction was then weighed in a quartz crucible and heated for 4 h at 550 ^∘C to burn organic matter. The Lee (2009) optimal protocol uses H₂O₂ to dissolve organic matter, but other studies have shown that this reagent is not always fully selective (Tessier et al., 1979). The ashing of organic matter yields similar results, but is normally conducted at the end of the procedure (Lee, 2009). By ashing first, we were able to process larger amounts of silt, limiting furnace contamination and aiding the process of dissolving precipitates (mainly oxides formed during heating) in the following steps. Oxides turn the samples a brown to reddish color.

3. Sample grains were dispersed in 0.1 N sodium oxalate and the ≤4 µm fraction was removed using Stokes settling. The ≤4 µm fraction consists of clays and also primary silicates with larger U–Th disequilibria.

4. For each sample, 4 g of sediment was weighed in a 50 mL centrifuge tube and leached twice with 32 mL 1 M magnesium nitrate to remove the exchangeable fraction and any residues of the ashing process.

5. 32 mL of acetic acid buffered to pH 5 with 1 M sodium acetate was then added and agitated at room temperature during 5 h to dissolve the carbonates.

6. 40 mL of 0.04 M hydroxylamine hydrochloride in 25 % (v∕v) acetic acid was added and heated at 95 ^∘C for 6 h twice to dissolve oxides. After this procedure, the samples still had a brown to reddish color from oxides formed in the furnace. These amorphous oxides can be resistant to the hydroxylamine hydrochloride treatment, as shown by Gontier (2014), who recommended using oxalic acid in ammonium oxalate (agitated in the dark at room temperature for 4 h) (Leleyter and Probst, 1999). This procedure preserves silicate from dissolution (Gontier, 2014). Note that for the last three steps, the leaching residues were centrifuged and rinsed twice with MQ water.

At this stage the samples should contain only primary silicates between 4 and 50 µm. To evaluate whether sample preparation fully removed secondary phases (e.g., carbonates, oxides, clays) without affecting primary silicates, we made scanning electron microscopy (SEM) images at the University of Strasbourg, as shown in Sect. 3.2. Our observations showed the method to be effective (see the Supplement for a further discussion on the efficacy of Th removal in this protocol) and the samples were deemed ready for disequilibrium analysis.

Two duplicates of this pretreatment protocol were measured (see Table 3). They are remarkably consistent within errors for (²³⁰Th∕²³⁸U) and (²³⁰Th∕²³⁴U) and only slightly different for (²³⁴U∕²³⁸U) (only twice the external uncertainty). Since the surface area data are also slightly different, this small difference could be associated with heterogeneities in the powder, fractionation of the grain size during the sampling of the powder or biases in the leaching process.

2.3.3 U–Th disequilibrium analysis

A first aliquot of the samples was dedicated to U–Th analysis. This aliquot was powdered in an agate ball mill to optimize sample homogeneity and acid dissolution. Approximately 100 mg of this powder was spiked with ²³³U–²²⁹Th tracer and digested with concentrated HF–HNO₃–HClO₄ acids following Pelt et al. (2013). Separation and purification of U and Th followed standard anionic resin chromatography (Dequincey et al., 2002; Granet et al., 2007; Pelt et al., 2008). U and Th isotopic ratios and concentrations were determined by standard sample bracketing (SSB) using an MC-ICPMS Neptune and following an optimized procedure upgraded from previously published protocols (Bosia et al., 2016; Ma et al., 2012; Pelt et al., 2008). The ²³³U–²²⁹Th spike was calibrated against the gravimetric NIST SRM 3164 and 3159 U and Th pure solutions. IRMM-184 and IRMM-035 isotopic standard solutions were spiked and used as bracketing solutions for U and Th SSB analysis, with the ²³³U∕²³⁴U and ²²⁹Th∕²³⁰Th ratios of these mixtures calibrated by TIMS Triton. We used the consensus value from Sims et al. (2008) for the ²³²Th∕²³⁰Th ratio of the IRMM-035 standard and our own absolute TIMS Triton measurements for the ²³⁴U∕²³⁸U ratio of the IRMM-184 standard instead of the certified values. The peak tailing of the ²³²Th over the ²³⁰Th was corrected using an exponential law fitted by signals measured in 229.6 and 230.6 masses. The precision and accuracy of the data were estimated from regular analysis of international pure solutions: HU1 for (²³⁴U∕²³⁸U) at 1.0008±0.0005 (2σ, N=5), IRMM-036 for ²³²Th∕²³⁰Th at 329 064±1313 (2σ, N=5) and rock basalts (BCR-2, AThO and BE-N). Results are consistent with published data and the ²³⁸U–²³⁴U–²³⁰Th secular equilibrium assumed for old (several million years) BCR-2 and BE-N “unweathered” basalts within errors. Based on these pure solutions and homogeneous basalt powders we estimate an uncertainty (2 SD) of 0.2 % for (²³⁴U∕²³⁸U), 0.5–1 % for (²³⁰Th∕²³²Th) and 1–1.5 % for (²³⁸U∕²³²Th), (²³⁰Th∕²³⁸U) and (²³⁰Th∕²³⁴U). A duplicate of the protocol was measured and agrees within the internal measurement uncertainties for (²³⁴U∕²³⁸U), (²³⁰Th∕²³⁸U), (²³⁰Th∕²³⁴U) and a rock standard BCR-2 (see Table 3). Usual blanks in the lab are in the range of 20–50 pg for U and 100–200 pg for Th and are therefore negligible here.

2.3.4 Estimation of the recoil loss factor

The recoil loss parameter (f_α) is a critical parameter in comminution age theory as presented in Sect. 2.3.1. Different methods have been proposed to estimate this parameter, as summarized in Maher et al. (2006) and Lee et al. (2010), but discrepancies up to an order of magnitude can be found between the different methods (Handley et al., 2013a, b). Here we used the N₂ gas absorption technique, which consists of surface area and fractal dimension measurements. Other methods are more subjective in that they involve visual estimation and assumptions for the aspect ratio and surface roughness coefficient. We measured the specific surface area for all samples and assessed the efficiency of this method to determine f_α and the relation between f_α and U–Th disequilibria.

A second aliquot of the samples (at least 1 g sample⁻¹) was prepared for these analyses, which were conducted at the Institut de Physique du Globe de Strasbourg on a Carlo Erba Sorptomatic 1990 machine. Details of the procedure are given in the Supplement. Fractal dimensions were also measured for each sample and the recoil loss factor f_α was calculated following Bourdon et al. (2009) after Semkow (1991):

$$\begin{array}{}\text{(4)}& {\displaystyle }{\displaystyle }{f}_{\mathit{\alpha }}={\displaystyle \frac{\mathrm{1}}{\mathrm{4}}}\left[{\displaystyle \frac{{\mathrm{2}}^{D-\mathrm{1}}}{\mathrm{4}-D}}{\left[{\displaystyle \frac{a}{R}}\right]}^{D-\mathrm{2}}\right]R\cdot S\cdot {\mathit{\rho }}_s,\end{array}$$

where S is the specific surface area, D is the fractal dimension, a is the diameter of the adsorbate molecule (0.35 nm), R is the recoil length (30 nm) and ρ_s is the density of the solid (2650 kg m⁻³).

2.3.5 Determination of the recycling time and weathering intensities: Monte Carlo analysis

A Monte Carlo approach was used to determine weathering intensities during both the pre- and post-deposition histories of the different moraine silts, and the duration of this pre-deposition history (called “recycling time” from this point forward).

¹⁰Be exposure ages provide information on the deposition age of the moraines. The age extracted from U–Th data provide information on the time since the sediment was produced, so it is equal to the sum of sediment recycling time (time elapsed between comminution and deposition of the moraine) and the deposition age. The problem as formulated in Sect. 2.3.1 has a large number of unknowns compared to the amount of data in each equation, so the problem is underdetermined. However, if geomorphology and climate did not vary too much between each glacial period, we can assume that recycling times and weathering intensities are comparable between the different samples. In fact, the glacial maxima observed since the mid-Pleistocene transition are quite comparable. We can then assume that the processes (especially weathering) occurring during a glacial cycle are comparable in intensity from one cycle to another. Regarding the recycling time, the other main parameter (apart from climate) influencing it is the morphology (mainly valleys and overdeepening). The Great Patagonian Glaciation occurred before 1 Ma (Mercer and Sutter, 1982); these glaciations are most probably responsible for carving the U-shaped valleys and the overdeepening in which the lake is currently. The results of the numerical experiments of Kaplan et al. (2009), regarding the extent of ice related to the elevation of the accumulation area and the depth of the overdeepening, suggest that this overdeepening was most probably already carved before 1 Ma. Recently, Christeleit et al. (2017) argued using thermochronometric data that the relief of Patagonia at the latitudes of our study was probably achieved between 10 and 5 Ma. On that basis, one can assume that the recycling time probably remained constant over the period covered by the samples. Hence all samples are assumed to share the same weathering intensity and the same pre-depositional history duration.

We used the theoretical scheme described in Sect. 2.3.1 to estimate these weathering intensities and recycling times. Known or informed parameters in the Monte Carlo simulations were the activity ratios (²³⁴U∕²³⁸U) and (²³⁰Th∕²³⁸U), the initial activity ratio (from the glacial flour sample), the recoil loss factor f_α and the deposition ages of the moraines (the post-depositional duration). Unknown parameters were the weathering intensities of ²³⁸U, ²³⁴U and ²³⁰Th for the pre- and post-deposition history. The relative weathering intensity of ²³⁴U compared to ²³⁸U, and ²³⁰Th compared to ²³⁸U, were assumed to be the same before and after moraine deposition. Consequently, we needed to determine five unknowns: the respective weathering coefficients $k_{\mathrm{238}}^{\text{pre}}$, $k_{\mathrm{238}}^{\text{post}}$, k₂₃₄∕k₂₃₈, k₂₃₀∕k₂₃₈ and the recycling time t_recycl using Monte Carlo simulations. We chose a large, arbitrary number of values for the unknowns and calculated the activity ratios associated with these sets of values. Then, we calculated a misfit M for the data as the difference between the measured and modeled values weighted by the uncertainty. The lower the misfit, the larger the probability of this set of values (we show the likelihood L for better visibility):

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }M=(m-o{)}^T{\mathit{\sigma }}^{-\mathrm{1}}(m-o)\\ \text{(5)}& {\displaystyle }& {\displaystyle }L=\exp (-M/\mathrm{2}),\end{array}$$

where m is the vector of modeled parameters (the five unknowns above), o is the observed values (the activity ratios (²³⁴U∕²³⁸U) and (²³⁰Th∕²³⁸U)) and σ is the diagonal matrix of uncertainties on the data. In theory, an optimal solution can be determined corresponding to the lower misfit (larger likelihood). However, no solutions were found because the number of unknowns was too great. We also know the relative weathering intensities of ²³⁸U, ²³⁴U and ²³⁰Th since Th is less soluble than U and ²³⁴U is on more fragile mineralogic sites. Hence, k₂₃₀ is lower than k₂₃₈, and k₂₃₄ is slightly larger than k₂₃₈. This information was used to further constrain the Monte Carlo simulations. For different couples of k₂₃₄∕k₂₃₈ and k₂₃₀∕k₂₃₈, we ran a Monte Carlo optimization to find an optimal solution (lower misfit). All optimal solutions for $k_{\mathrm{238}}^{\text{pre}}$, $k_{\mathrm{238}}^{\text{post}}$ and the recycling time are shown in Fig. 6, and we use the optimal solutions to characterize the pre-depositional duration of these sediments.

To test whether exchangeable Th might not be efficiently leached by our chemical protocol, we also made a simulation with $f_{\mathit{\alpha }}^{\mathrm{230}}=\mathrm{0}$ and $f_{\mathit{\alpha }}^{\mathrm{230}}={\displaystyle \frac{\mathrm{1.176}\times f_{\mathit{\alpha }}^{\mathrm{234}}}{\mathrm{2}}}$ (half of the theoretical value described in Sect. 2.3.1). If Th was totally readsorbed onto the mineral surfaces after ejection from the grain (and not subsequently removed during leaching), then it is equivalent to having no ejection ($f_{\mathit{\alpha }}^{\mathrm{230}}=\mathrm{0}$). The case with $f_{\mathit{\alpha }}^{\mathrm{230}}={\displaystyle \frac{\mathrm{1.176}\times f_{\mathit{\alpha }}^{\mathrm{234}}}{\mathrm{2}}}$ is an intermediary between total removal of the adsorbed Th and no removal at all, as discussed further in the Supplement.

3 Results

3.1 Cosmogenic ¹⁰Be and deposition ages

The results of ¹⁰Be analyses are shown in Fig. 2 and Tables 2 and 1. Refer to the Supplement for more details on the profile age results.

Figure 2 ¹⁰Be exposure ages from outwash surface cobbles (circles) and modeled depth profiles (squares) for five of the Lago Buenos Aires moraines. Depth profile uncertainties are 1σ. Symbols and colors are the same as in Fig. 1. The δ¹⁸O benthic stack record from Lisiecki and Raymo (2005) is also plotted, and gray shading shows glacial stages.

Figure 3 Solutions to Eq. 3 that we use as the theoretical framework to analyze our data for different values of f_α and k₂₃₈ (we consider $k_{\mathrm{234}}/k_{\mathrm{238}}=\mathrm{1.4}$ and $k_{\mathrm{230}}/k_{\mathrm{238}}=\mathrm{0.5}$). The time evolution of ²³⁴U∕²³⁸U (a), ²³⁰Th∕²³⁸U (b) and ²³⁰Th∕²³⁴U (c) activity ratios. The vertical dashed line represents an event through which the weathering intensity varies: this corresponds in our case to the time of sediment remobilization. Our ¹⁰Be exposure ages give the time elapsed since this event (shown by the dashed vertical line), whereas the U–Th data indicate the duration of the first phase since comminution and remobilization (from the x axis origin to the dashed vertical line). (d, e and f) The co-evolution of the ratios in (a), (b) and (c), respectively. Different scenarios are plotted corresponding to different f_α values and weathering intensities to illustrate potential patterns in the comminution weathering model (see legend and Sect. 2.3.1). Observing these figures helps in understanding the range of variations in the different activity ratios and the effect of each parameter (α recoil and weathering) on these ratios before and after the remobilization of sediments.

Table 1 Cosmogenic ¹⁰Be nuclide data for outwash surface cobbles and depth profile samples. Moreno 1 and 3 data are already published in Hein et al. (2017).

Table 2 Modeled exposure ages for outwash surface cobbles and depth profile samples (see main text and the Supplement for details on profile modeling). Moreno 1 and 3 ages have been already published in Hein et al. (2017). Bold values are the oldest samples for each outwash that we consider to be the most representative deposition age of the moraine (see Sect. 3.1 for details).

Eight surface cobbles from outwash relating to the Moreno 1 and 3 moraines yielded ¹⁰Be exposure ages ranging from 168–269 ka (Hein et al., 2017). To these we add six more surface exposure ages and two depth profile ages from outwash related to the Deseado and Telken moraines. Three cobbles from Deseado 1 produced relatively tightly clustered exposure ages of 430–468 ka. A cobble from the ice-distal Deseado 2 outwash yielded a much younger, stratigraphically inconsistent age of 293 ka (sample D2-T12), similar to the Moreno ages. Two further cobbles from the Deseado 2 outwash yielded ages of 520 and 618 ka, the oldest of which agrees within errors with a modeled depth profile age of 600 ka ($+\mathrm{70}/-\mathrm{35}$ kyr) from the same outwash unit. The agreement between older ages from surface cobbles and the Deseado 2 depth profile is unsurprising given that modeling suggests that surface erosion (<0.2 cm ka⁻¹) and inheritance within the outwash sediments (<2.10⁴ atom g⁻¹) was relatively low (see the Supplement). There are a number of factors that can result in some age scatter within outwash surface cobbles, including surface deflation and up-freezing of clasts (Darvill et al., 2015; Hein et al., 2017). Moreover, Hein et al. (2009, 2017) suggested that the oldest surface exposure ages from outwash in this region will most closely relate to deposition age. The D2-T12 sample is >200 kyr younger than the other Deseado 2 cobbles, and it is likely that this cobble was either deposited at a much later time (perhaps during Moreno-stage glaciation, although the difference in altitude makes this unlikely) or was not in situ. The close agreement between other Deseado 2 ages supports the removal of D2-T12 as an outlier. Finally, a modeled depth profile through outwash relating to the Telken 5 moraine produced an age of 780 ka ($+\mathrm{170}/-\mathrm{80}$ kyr, 1σ confidence interval), with relatively low surface erosion (<0.1 cm ka⁻¹) and inheritance (<6.10⁴ atomg g⁻¹). This age is within the lower range of 760–1016 ka known from radiometric ages (Singer et al., 2004). Sensitivity tests showed that modifications of the a priori values of erosion rates and inheritance lead to insignificant changes in the resulting age.

3.2 Evaluation of the U–Th chemical procedure with SEM images

For three samples taken randomly, we looked at SEM images after 50 µm sieving and preparation (Fig. 4) to assess the effects of our chemical and mechanical treatment. We were particularly interested in whether secondary phases could remain after processing and if the primary silicates could have been altered by this treatment. The minerals observed are mainly quartz, feldspars and micas corresponding to the lithology of the Andes in these latitudes. We list below the main observations.

• No trace of carbonate, organic matter or oxides were observed in samples that underwent the entire protocol. Prior to treatment, we only note the presence of oxides (Fig. 4b). This does not, however, preclude the presence of carbonate or organic matter since a difference in mass before and after steps 2 and 5 of the protocol was observed. Therefore, the protocol seems efficient at eliminating carbonates, oxides and organic matter.

• Samples were generally cleaner after preparation. This is because ≤4 µm fractions containing mainly clays were removed by Stokes settling. However, clays are still observed in the samples after preparation (see Fig. 4b and c), mainly agglomerated around primary grains as clay pellets. In one sample, up to one-third of the grains have clay pellets. These clays could have been precipitated directly at the surface of the grains during weathering or agglomerated subsequently. It is difficult to fully eliminate the clays and a potential bias could be introduced if too many clays remain in the sample. This potentially limits the application of comminution ages to samples with a very low amount of clays, since their complete removal is currently not possible.

• Most micas show surfaces having experienced weathering (Fig. 4d). The shapes observed are typical of weathering both before and after the protocol. This shows that our silt samples have experienced weathering and that it could potentially have affected U–Th disequilibrium.

• Following the preparation protocol, two feldspar grains in one sample displayed evidence of corrosion (Fig. 4e and f). While this damage could be associated with alteration of the primary silicates during leaching, these are the only two corroded grains visible in the samples we tested. Since the corrosion pits are smaller than 4 µm and all the grains are larger than 4 µm in our samples, ubiquitous damage would be observable on other grains. As such, we assume that primary mineral corrosion during the protocol is minimal.

Figure 4 Scanning electron microscopy (SEM) images of selected samples. (a) Overview of a sample after the preparation process described in Sect. 2.3.2. Grains vary in size and shape between 4 and 50 µm and are angular to subangular. White grains are heavy minerals such as zircon or oxides. (b) The same sample before 50 µm sieving. The only notable difference is the presence of clays. (c) Zoomed image of a mineral grain covered in clay pellets that were not removed by the preparation process. (d) Mica grains with cauliflower facies, which is characteristic of weathering. (e, f) Feldspar grains with corroded surfaces. It is unclear whether these corrosion features result from geological processes or the sampling procedure, but they were only observed twice.

3.3 U–Th – f_α – age

U–Th results are shown in Fig. 5 and Table 3. (²³⁴U∕²³⁸U) values from moraines display little variation, ranging from 0.95 to 0.97, while glacial flour has a composition of 1.007 and is close to equilibrium. For the moraine samples, a general decrease in (²³⁴U∕²³⁸U) is observed with time, though with significant noise. (²³⁰Th∕²³⁸U) values mostly vary between 0.94 and 0.97 with the youngest moraine sample (F1-U-M2, Fenix moraine) having a composition around 1. Similarly, (²³⁰Th∕²³⁴U) values mostly range from 0.99 and 1.01, with the youngest moraine around 1.03. Again, an overall decrease in (²³⁰Th∕²³⁴U) is observed with time. The glacial flour is close to equilibrium compositions with $({}^{\mathrm{230}}\mathrm{Th}/{}^{\mathrm{238}}\mathrm{U})=\mathrm{1.02}$ and $({}^{\mathrm{230}}\mathrm{Th}/{}^{\mathrm{234}}\mathrm{U})=\mathrm{1.01}$. Specific surface areas vary from 1.5 to 7 m² g⁻¹, which is typical for this kind of sample, and the fractal dimensions range from 2.45 to 2.6. The recoil loss factors f_α calculated using Eq. 4 are between 0.005 and 0.025, with no observable trend with time.

Figure 5 2-D cross plots of ²³⁴U∕²³⁸U and ²³⁰Th∕²³⁸U activity ratios, associated recoil loss factors ($\mathrm{1}-f_{\mathit{\alpha }}^{\mathrm{234}}$) and $(\mathrm{1}-f_{\mathit{\alpha }}^{\mathrm{234}})\times (\mathrm{1}-f_{\mathit{\alpha }}^{\mathrm{230}}$), and ¹⁰Be exposure ages. Symbol colors represent the third parameter in each plot. The sorting observed in the symbol colors helps to explain the noise associated with the other two parameters. In particular, the scatter in the decrease of the activity ratios over time can be largely explained by the variability in the recoil loss factor (top left plot). The plots also let us assume that the variability in the weathering coefficients could be low; otherwise a residual dispersion would be observed. In all cases, error bars are smaller than the size of the dots.

Table 3 U–Th concentrations, activity ratios and BET data (specific surface areas, fractal dimensions and calculated recoil loss factors). Duplicates are also shown, two for the whole process (including sieving, leachings, HF digestion and chromatography) and one for HF digestion and chromatography. The samples are sorted by deposition age (F1: Fenix, M1: Moreno 1, M3: Moreno3, D1: Deseado 1, D2: Deseado 2, TI: Telken 1, TK: Telken 5).

Handley et al. (2013a) analyzed only a few samples with the gas adsorption technique to characterize the specific surface area and f_α. We analyzed every sample with this technique, highlighting the consistency of the noise in (²³⁴U∕²³⁸U) and the dispersion in f_α derived from specific surface area (Fig. 5). The relationship between (²³⁴U∕²³⁸U), f_α and age is described by a surface in three dimensions. In other words, the noise observed in the (²³⁴U∕²³⁸U) age diagram around the general decreasing pattern is in large part associated with the dispersion of f_α. This suggests that specific surface area data are consistent and reliable or that, if a bias exists in these data, it is systematic. Such a systematic bias would arise from a rather random event due to sample processing, so this latter explanation seems unreasonable. We conclude that the specific surface area determination based on the gas adsorption technique is valid to evaluate the recoil loss factor in the U–Th comminution age theory.

We observe that (²³⁰Th∕²³⁴U) is larger than (²³⁴U∕²³⁸U) and (²³⁰Th∕²³⁸U) has comparable compositions to (²³⁴U∕²³⁸U). As mentioned in Sect. 2.3.1, this observation is incompatible with a simple comminution age model (only α recoil). An enrichment of ²³⁰Th compared to ²³⁴U could be the imprint of weathering on silts found in the moraines. Similarly, the incompatibilities of (²³⁴U∕²³⁸U) and f_α in the framework of the simple comminution age model suggest that weathering must be considered. In the comminution age theory as described by DePaolo et al. (2006), 1−f_α must always be lower than (²³⁴U∕²³⁸U), and for a sample old enough to have reached steady state, $({}^{\mathrm{234}}\mathrm{U}/{}^{\mathrm{238}}\mathrm{U})=\mathrm{1}-f_{\mathit{\alpha }}$. This is not what we observe, similarly to Handley et al. (2013a). However, as described in Sect. 2.3.1 such an observation can be explained if weathering is considered and k₂₃₄≥k₂₃₈.

Given the large disequilibrium of the youngest samples in (²³⁴U∕²³⁸U) and in (²³⁰Th∕²³⁴U) (Table 3 and Fig. 5) and the fact that the glacial flour (the most probable initial composition for the moraine silts) is close to equilibrium, we cannot ignore the pre-depositional history (i.e., before deposition in the moraines). Based on a starting value close to equilibrium, the disequilibrium values measured here cannot be reached in only 20 or even 200 kyr. Since our estimations of the f_α seem consistent and reliable, this requires consideration of a pre-depositional history involving weathering with a different intensity than during the post-deposition history.

3.4 Pre-deposition history and Monte Carlo analysis

We attempt to constrain the pre-depositional history using a Monte Carlo analysis. The results show an inverse relationship, with the greatest probability for $k_{\mathrm{238}}^{\text{pre}}/k_{\mathrm{238}}^{\text{post}}$ between 2 and 3 and a recycling time between 100 and 200 kyr. The optimal values of k₂₃₄∕k₂₃₈ and k₂₃₀∕k₂₃₈ associated with this solution describe an anticorrelation with $\mathrm{1.2}\le k_{\mathrm{234}}/k_{\mathrm{238}}\le \mathrm{2}$ and $\mathrm{0.01}\le k_{\mathrm{230}}/k_{\mathrm{238}}\le \mathrm{0.8}$ (Fig. 6). These parameters are, however, poorly constrained by the inversion process.

Figure 6 Results of Monte Carlo simulations for the U–Th data. The same Monte Carlo simulations were conducted for different fixed values of k₂₃₄∕k₂₃₈ and k₂₃₀∕k₂₃₈ (see Sect. 2.3.5). For each simulation we estimated the optimal solution for $k_{\mathrm{238}}^{\text{pre}}$, $k_{\mathrm{238}}^{\text{post}}$ and recycling time (pre-deposition duration, after comminution), i.e., the solution with the lowest misfit between measured and calculated activity ratio data. We used this misfit to calculate a likelihood between 0 and 1 (see Eq. 5), as shown by the color map (red is lower misfit, blue is higher misfit; stretched to fit our results). The larger the likelihood, the closer to the data the solution is, so the more likely the set of parameter values. (a) The ratio $k_{\mathrm{238}}^{\text{pre}}/k_{\mathrm{238}}^{\text{post}}$ as a function of the recycling time. Each dot is one Monte Carlo solution. The optimal solution (red dots) is for 3 times stronger weathering before deposition and a recycling time of 100 kyr. (b) Relative weathering intensities of ²³⁴U, ²³⁸U and ²³⁰Th. The relative weathering intensities are poorly constrained, but a relationship can be extracted between their ratios: the larger k₂₃₄∕k₂₃₈, the smaller k₂₃₀∕k₂₃₈ must be. (c) A relation between $k_{\mathrm{238}}^{\text{post}}$ and k₂₃₄∕k₂₃₈ can be extracted, even if both parameters are poorly constrained. The estimation of sediment transfer times using U–Th disequilibria series could be improved with better constraints on the relative mobility of ²³⁴U and ²³⁸U.

Despite the numerous parameters to be constrained, the model converges for key results as shown above. We show in particular that, in the framework of the modified comminution age theory incorporating weathering, as described in Sect. 2.3.1, the silts from the LBA moraines were likely eroded 100 to 200 kyr before being deposited in the moraines and that on average they experienced around 2–3 times more intense weathering during this interval than after deposition (Fig. 6). The best fit corresponding to a recycling time of 180 kyr, a ratio $k_{\mathrm{238}}^{\text{pre}}/k_{\mathrm{238}}^{\text{post}}=\mathrm{2.35}$, $k_{\mathrm{234}}/k_{\mathrm{238}}=\mathrm{1.4}$ and $k_{\mathrm{230}}/k_{\mathrm{238}}=\mathrm{0.6}$, is shown in the Supplement Fig. S5.

4 Discussion

4.1 Glacial chronology and perspectives on the control of ice extent in the Patagonian Andes

The new ¹⁰Be exposure ages in this study help clarify the timing of the deposition of several LBA moraines. The oldest outwash cobbles are taken as closest to the age of deposition because depth profiles show that average nuclide inheritance is relatively low and outwash surfaces show evidence of deflation (Hein et al., 2011, 2017). It is possible that even the oldest cobbles underestimate the age of deposition if they have also been exhumed, and we apply no erosion correction to our exposure ages. As discussed in Hein et al. (2017), outwash surface cobbles imply that the Moreno moraines were deposited at ca. 260–270 ka during marine isotope stage 8.

Unlike the Moreno system, surface cobbles from the Deseado moraines appear to date from two different glacial cycles. Exposure ages from Deseado 1 outwash yield relatively tightly clustered ages of 430–470 ka, suggesting that the limit was deposited during MIS 12. The older of two published moraine boulder exposure ages from the Deseado 1 moraine (Kaplan et al., 2005) is consistent with the cobble ages in this study. In contrast, surface cobbles from Deseado 2 outwash yield ages of 520 and 618 ka (excluding the anomalously young D2-T12 age of 293 ka). The dating would be less conclusive without the accompanying depth profile age of 600 ka $+\mathrm{70}/-\mathrm{35}$ kyr. Slight changes in scaling factors (see Sect. 2.2.3) do not affect this finding. Taken together, the oldest surface cobble and depth profile suggest that the Deseado 2 limit was deposited at ca. 600–620 ka. Within errors, the limit may relate to MIS 16. The scatter in ages is interesting given that luminescence dating yielded an even younger age of 123±18 ka for the Deseado 2 limit (Smedley et al., 2016). This scatter may imply greater sediment reworking or a more complex relationship between moraines and outwash in this sequence. The Telken 5 depth profile gave an age of 780 ka $+\mathrm{170}/-\mathrm{80}$ kyr that is stratigraphically consistent within our chronological dataset, but less helpful in determining when the limit was deposited. The error range spans MIS 18–24 (and radiometric datings gave a lower age of 760 ka for the Telken series), with larger probability around MIS 20, so it is possible that the moraine was deposited during MIS 20. Since the entire Telken moraine system must be older than 760 ka (Singer et al., 2004), this implies that Telken 1 to 4 are also that age and that the Telken moraines represent different glacial advances during the same glacial cycle.

In summary, there is a clear pattern in the timing of moraine deposition in which alternate glacial cycles are represented in our chronology: Telken 5 during MIS 20; Deseado 2 during MIS 16; Deseado 1 during MIS 12; and the Moreno moraines during MIS 8. Previous work has shown that the innermost Fenix moraines relate to the Last Glacial Maximum during MIS 2 (Douglass et al., 2006; Kaplan et al., 2004; Smedley et al., 2016). We strongly caution that scatter in surface cobble ages and error ranges in depth profiles may complicate this pattern, particularly for the older limits. Moreover, we did not analyze the enigmatic Deseado 3 moraine and so cannot say whether this relates to the counterparts dated here or intervening glacial stages. However, the dating of the Moreno and Deseado 1 and 2 moraines does imply that the intervening glacial cycle is absent from the record.

The pattern in the timing of moraine deposition around Lago Buenos Aires implies that there are a few glacial cycles that are either not recorded in this area or have been erased or removed. Either the ice lobe did not advance during alternate glacial cycles, or it advanced to similar or less extensive positions so that moraines and outwash were then removed by the following advance. A similar pattern was observed by Hein et al. (2009, 2011) for the Lago Pueyrredón ice lobe 100 km to the south, so there could be a regional driver of alternate glacial advances. One possibility is that erosion over successive glacial cycles caused entrenchment of ice lobes within large basins on the eastern side of the Patagonian Ice Sheet (Kaplan et al., 2009). This erosion model has been linked to the pattern of nested limits seen across the former ice sheet (Anderson et al., 2012; Kaplan et al., 2009) but may be more complex if only alternate cycles are represented in the moraine record. Alternatively, a purely climatic forcing could have caused alternate strong advances, although there is no simple relationship between available climate records and alternate glacial advances at Lago Buenos Aires. A complex erosion–climate feedback mechanism may determine when or how far glacial advances occurred in this region, but more detailed glacial models are required to further test such a model.

The chronology presented here suggests that moraine reworking has probably occurred in the area and/or that glacial erosion in the Andes may feed back into the glacial lobe advance to produce this observed pattern of absent intervening glacial cycles. Regardless, this chronology allows us to broadly constrain deposition ages of moraines targeted for U–Th analysis, which may in turn inform sediment recycling times between preserved moraines.

4.2 Implications for sediment recycling, chemical weathering, climate and their interactions in proglacial systems

Following the conclusions of Handley et al. (2013a), the residence time of 100–200 kyr could be an artifact due to an addition of old dust. In our area, when a new moraine is being formed, the older sediments are on the east, whereas the dominant winds come from the west. It is therefore unlikely that these winds could add older material to that being deposited.

The long time estimated using the whole set of data and the Monte Carlo simulations, along with the absence of inheritance in the ¹⁰Be data, suggest that the sediment was not exposed at the surface in the proglacial system before being deposited in its moraine. This implies that the sediment must have been buried, either in the proglacial lake or within a sedimentary pile (moraine, till, channel, etc., at least a few meters below the surface) before final deposition.

The dated sediment was most probably eroded during glacial periods, as erosion rates are likely much larger than during an interglacial. Likewise, deposition in an end moraine occurs during a glacial period. This means that the sediment deposited in the moraine may spend on average as long as 100–200 kyr in the proglacial system (till, lake sediments, etc.) before being deposited in the moraines. Hence the sediment eroded during a glacial period is on average deposited in the frontal moraine during the next glacial cycle. Using postglacial sediment budget and reservoir theory, Hoffmann and Hillebrand (2016) modeled the residence time of sediments in a periglacial system in the Canadian Rocky Mountains and also found a time of 100 kyr. This average residence time implies that a fraction of the sediment deposited in a moraine may have been produced during the same glacial cycle with another fraction being much older.

This time appears to be long. However, our measurements of cosmogenic exposure ages using ¹⁰Be nuclide concentrations in outwash cobbles and profiles suggest that the preserved moraines represent every other glacial cycle, indicating an advance over the proglacial sediment of one or two previous advances. Moreover, it has been shown that proglacial lakes occupying overdeepenings are filled following deglaciation (Eyles et al., 1991; Houbolt and Jonker, 1968). The lakes are filled with sediments during interglacial periods and emptied during the subsequent glacial periods. These processes imply sediment erosion and deposition over a full glacial–interglacial cycle, which is consistent with both ¹⁰Be and U–Th disequilibrium data.

We also obtain a weathering rate that is 2–3 times higher during the 100–200 kyr of the pre-depositional phase than after deposition. Based on the argument above, these sediments are exposed to water in the proglacial system before deposition in the moraine. Authors such as Anderson (2005) have shown that proglacial environments favor weathering. Since weathering rate is a time-dependent parameter (White and Brantley, 2003), the calculated recycling time of 100–200 kyr likely has a non-negligible impact on weathering fluxes.

Understanding how sediments are evacuated and transported to the oceans and when they experience weathering would help to quantify the relationships between erosion and climate. In particular, if 100 to 200 kyr are needed to escape the periglacial area, important lags could be observed between an erosional forcing and the weathering and climate response. This suggests that measured weathering variations could have occurred over the last glacial–interglacial cycles or could have been smoothed and/or damped because of a lag in the sediment transport. Estimating these lag times would help in understanding the relationships between erosion, climate and marine biogeochemical cycles. Erosion rates can be determined using thermochronology or cosmogenic isotopes. Sedimentation rates in the ocean can also be estimated using marine core dating methods. The link between erosion and sedimentation is poorly known (Sadler, 1981), especially because the transport times are poorly constrained. The transport time of 100 to 200 kyr presented here suggests that the pathway between initial erosion and deposition is potentially complex. An effort to constrain these transport times appears to be potentially fruitful to reveal the actual link between erosion and sedimentation rates and climate.

4.3 Perspectives on comminution age method and its applications to characterize sediment transfer

Our data let us estimate only a mean recycling time for moraine sediment (the fine silty fraction). Here we have discussed it in terms of the recycling of the previous glacial cycle moraine, but it could also be interpreted as a mixing of much older, deeper, reworked sediment with new freshly eroded sediment. In other words, we do not quantify the amount of sediment escaping the proglacial system or the amount of sediment rapidly deposited in the moraine after erosion. It would be beneficial to think about age distributions and not only mean ages. However, this necessitates being able to measure U–Th disequilibrium on single grains, which remains an analytical challenge (Bosia et al., 2018).

Quantifying or reducing the effects of weathering remains a major challenge. Pure primary minerals with no clays would minimize these secondary processes. To this end, working on pure zircon grains could be an option. This would require improved mineral separation methods for fractions smaller than 50 µm. One would also have to assume that comminution effectively occurred on such a heavy mineral and that there is no initial disequilibrium (a problem being that zircons are much more enriched in U than the surrounding minerals).

5 Conclusions

Our study is a first attempt to quantify long-term sediment transfer times in a proglacial area. We take advantage of the particularly well-preserved series of nested moraines of the Lago Buenos Aires in Patagonia. Our approach involves ¹⁰Be exposure dating of these moraines combined with U–Th disequilibrium measurement on the fine fraction of moraine sediment, within a “modified” comminution theoretical framework, to characterize the pre-deposition and post-deposition histories of the sediments.

We show that weathering cannot be neglected when determining the comminution age. Measuring Th isotopes helps constrain the weathering process. Even with this complication, it is possible to calculate comminution ages if the studied samples have well-constrained deposition ages and have experienced the same pre-depositional history. In this case, both the weathering intensities and the pre-depositional history duration are the same for all the samples. We also show that specific surface area measurements based on the gas adsorption technique for estimation of the recoil loss factor is a reliable method and should be applied to all samples. One caveat is that clay minerals may not be removed from the samples, especially clay pellets agglomerated at some mineral surfaces, which could be a strong limitation of the use of this method to calculate transport times.

Our new chronology shows that some glacial cycles are either not preserved or not recorded in the LBA area; i.e., the associated moraines may have been removed by following glacial advances or may have been deposited farther upstream. Hein et al. (2011) found a similar pattern in the Lago Pueyrredón area, suggesting that the mechanisms driving glacial advance may have a regional extent. Future work may reveal whether erosion feedbacks in the Andes could be responsible for this pattern.

The absence of moraine records for a few glacial cycles in the LBA area can be explained by our U–Th data showing that there has been reworking and/or recycling of the sediments. Using a Monte Carlo approach we estimate that the silts from the Lago Buenos Aires frontal moraine system have a 100 kyr residence time in the proglacial system (transported in the lake sediments, remobilized from the previous moraines or stored in intermediate reservoirs) and experienced 3 times more intense weathering before deposition than in the moraine. Our data represent a step forward in the effort to constrain the pathways and timescales over which sediments are transported from source (erosion in mountain belts) to sink (ocean sediments). Although considered a daunting challenge by Sadler and Jerolmack (2015) and requiring additional development, this approach can contribute to our understanding of basin-scale sediment budgets and how erosion, weathering and sedimentation evolve through time.