Research ArticleGEOLOGY

Transient glacial incision in the Patagonian Andes from ~6 Ma to present

See allHide authors and affiliations

Science Advances  12 Feb 2020:
Vol. 6, no. 7, eaay1641
DOI: 10.1126/sciadv.aay1641


We report a mountain-scale record of erosion rates in the central Patagonian Andes from >10 million years (Ma) ago to present, which covers the transition from a fluvial to alpine glaciated landscape. Apatite (U-Th)/He ages of 72 granitic cobbles from alpine glacial deposits show slow erosion before ~6 Ma ago, followed by a two- to threefold increase in the spatially averaged erosion rate of the source region after the onset of alpine glaciations and a 15-fold increase in the top 25% of the distribution. This transition is followed by a pronounced decrease in erosion rates over the past ~3 Ma. We ascribe the pulse of fast erosion to local deepening and widening of valleys, which are characteristic features of alpine glaciated landscapes. The subsequent decline in local erosion rates may represent a return toward a balance between rock uplift and erosion.


Late Cenozoic cooling is marked by the onset of widespread periodic glaciations, starting in polar regions at 25 to 10 million years (Ma) ago, and expanding to midlatitude mountains and continental regions at ~6 to 2.5 Ma ago. There has been a long debate focused on the onset of the late Cenozoic “icehouse” climate and how it may have changed the size and shape of mountain ranges around the world and affected the delivery of continental-derived sediment to the oceans [e.g., (14)]. This debate has been strongly influenced by studies of modern erosion, which show that, at local scales and short time intervals, temperate glaciers are usually more erosive than rivers (1). However, erosion faster than rock uplift would decrease the height of the glacial catchment relative to the equilibrium line altitude (ELA) and subsequently reduce the amount of ice discharge (5, 6). This negative feedback operates at the scale of the glacial network and will therefore have a response time that is longer than that for the erosion processes operating along the bed of the glacier. This relationship suggests that the onset of mountain glaciations might include an initial phase of fast erosion, followed by a return to rates that are balanced with the local rates of rock uplift. The history of alpine glaciation in Fiordland, New Zealand over the past 2 Ma provides a possible example of this transient response, where the headward carving of deep valley troughs is accompanied by a flattening of the upland landscape (7), which would have reduced the amount of ice discharge.

Low-temperature thermochronology provides a way to study erosion at a regional scale and on a time scale of millions of years [e.g., (4, 711)]. The “cooling age” measures the transit time of a sample from a closure depth to the Earth’s surface. Low-temperature thermochronometers are well suited for this kind of study because the closure depth is fairly close to the Earth’s surface. For example, the (U-Th)/He apatite (AHe) system has an approximate closure temperature of 65°C (10, 12) and an approximate closure depth of 2.5 km, assuming a surface temperature of ~0°C (as expected in a mountain setting) and an average continental thermal gradient of ~25°C/km.

Past studies of this kind have focused on bedrock cooling ages, which use bedrock samples collected from the modern landscape surface. Here, we report detrital cooling ages, which supply samples from older bedrock surfaces, and thus provide a longer record of transit times. The transit time is estimated by the lag time (13), which is the difference between the cooling age of the detrital sample and the age of the deposit containing the sample. This approach carries the assumption that transport from bedrock source to the site of deposition is short relative to the transit time from the closure depth to the Earth’s surface. Mountain landscapes generally lack settings with substantial residence time, as judged by the age of upland deposits [see discussion in (14)]. Exceptions to this, such as the intermontane basins of the Central Andes, are easy to recognize and avoid. The lag time concept is best applied in settings adjacent to the piedmont transition, where the erosive flux of the bedrock landscape is captured by aggradation in the piedmont region, a conclusion that holds for both fluvial and glacial transport processes.

Our study focuses on distinctive granitic cobbles, which are a widespread component of glacial deposits along the eastern flank of the Patagonian Andes and are known to have been sourced from the Patagonian batholith, exposed in the core of the range (Fig. 1). These cobbles are used to estimate lag times, as defined by their AHe ages and the ages of the glacial deposits in which they are found. We report lag times from 72 cobbles, collected from four glacial deposits with ages from ~6 Ma to 18.5 thousand years (ka), and compile 51 published AHe ages from the modern bedrock landscape. These data are used to construct a record of regional-scale erosion rates, extending back to 10 Ma ago and earlier. This record provides direct evidence about the transition of this midlatitude mountain range from a fluvial landscape to a largely glacial one.

Fig. 1 Map of study area.

Study area with the modern exposure of the Patagonian batholith, faults, and sample locations. LBA, Lago Buenos Aires; MLBA, Meseta del LBA; NPI, North Patagonian Icefield; CTJ, Chile Triple Junction; LOFZ, Liquiñe-Ofqui fault zone. Geologic units other than the batholith and the glacial deposits are not shown. The extent of the Guivel and Mercer till exposures is not visible at this scale, as the meters-thick units are intercalated with basalt flows and outcrop in narrow bands that follow the contour of the meseta edge. The four sampling locations are marked by stars. Published bedrock AHe data (single samples or sample suites) are shown by squares (15), circles (20), diamonds (8), and pentagons (16). Large black symbols indicate data that appear in the top of Fig. 3 and in table S1. Small white symbols indicate other published bedrock ages, which are shown for reference but not used for our analysis. Faults in the region are shown in orange (15).

Conventional studies of bedrock cooling ages are limited to samples currently exposed at the bedrock surface. The vertical sampling method can be used to estimate past erosion rates [e.g., (15, 16)], but it generally provides a short record of erosion rates and it is known to produce biased results, as surface topography–induced vertical deflections of the closure temperature isotherm will depend on the rate at which surface relief is increasing or decreasing (17). In contrast, detrital ages provide a more comprehensive sampling of AHe cooling ages at the regional scale, and the interpretation is focused on the transit time from the mean closure depth to the mean surface elevation. As a result, this approach is not influenced by the change-of-topography biases associated with age-elevation transects.


Glacial deposits are located along the entire eastern margin of the Patagonian Andes and generally coincide with the termini of glacial advances (Fig. 1). Their exceptional preservation is due to the aridity and low topographic gradient in the eastern Patagonian foreland and, in some cases, due to burial by basalt flows, which are widespread south of 46°S. These deposits are typically conglomeratic, with rounded clasts set in fine-grained matrix. The clasts include volcanic, metamorphic, and granitic material, all of which can be traced to sources in the west within the higher part of the range. The proportion of granitic clasts is 5 to 20% in these deposits. The present surface exposures of the batholith have igneous ages of 155 to 115 Ma and igneous pressures indicating initial emplacement at depths of ~10 km (18). There are minor intrusions with ages as young as 10 Ma (18), although these were emplaced well below the AHe closure depth. As a result, the AHe ages from these granites are related to cooling during erosion and are not influenced by postmagmatic cooling. Figure 1 shows ice flow lines, which were determined from a reconstruction of glacial topography during the Last Glacial Maximum (LGM) (Fig. 1). At its previous LGM size and position, the Patagonian ice sheet transported material from the exposed batholith to the locations of the sampled deposits.

Our study site is located to the east of the Chile Triple Junction (CTJ) (Fig. 1), where the Nazca ridge, which separates the Nazca and Antarctic plates, subducts below the South American plate. Hypotheses regarding the tectonic and thermal effects of the CTJ include “collision deformation” by the Nazca ridge as it subducts (19), heating and/or dynamic uplift caused by an “asthenospheric window” (20), and strike-slip slivering of the forearc north of the triple junction (15). As of yet, there is little direct evidence to support these ideas. Recent studies (15, 21) infer that there are many local strike-slip faults east of the CTJ, but there is no direct exposure of these structures. One of these faults, located at the southern end of the Liquiñe-Ofqui fault zone (LOFZ), has direct evidence of offset (22). Most other features ascribed to the LOFZ are based primarily on the linear appearance of fjords.

In contrast to proposed CTJ-associated tectonic and thermal effects, a recent paleotopography study (18) shows that the topography around 46°S, the latitude of the CTJ, has been steady over the past 60 Ma. In another nearby study, the authors demonstrate that the region east of CTJ has been a site of relatively slow erosion and uplift over the past 6 Ma, except for a period of fast valley incision associated with the onset of glaciation (16). The evidence for glacial incision is widespread and most markedly demonstrated by the deep fjords and overdeepened lakes that characterize this region, with maximum depths extending to 1468 m below sea level (16). Current models for glacial erosion show that fjords and overdeepening can form only where the background rate of rock uplift is low (23). The Patagonian Andes are crisscrossed by a dense network of fjords, and our study area is flanked to the east by Lago Buenos Aires, an overdeepened glacial lake with a maximum depth of 586 m. Last, regional-scale thermochronologic data show that samples south of the CTJ yield older cooling ages relative to those to the north [e.g., (8)]. Proposals for ridge collision and dynamic topography predict a northward propagating region of young uplift and erosion, which would yield the opposite of what has been observed. This summary is meant to highlight an ongoing debate in this area about the role of the CTJ.

On the basis of an assessment of the available information, we conclude that the design of our experiment and the analysis and interpretation of the results are not influenced by the issues outlined above. Our study area lies ~200 km east of the CTJ and is east of the easternmost proposed strike-slip fault (Cachet fault; Fig. 1) (15). Thermal and thermochronologic modeling (16, 21, 24) indicates that the hot conditions associated with subduction of young lithosphere in the vicinity of the CTJ is too deep and too recent to have affected the shallow thermal field (<65°C) associated with the AHe cooling ages used in this study. Some may assume that thermochronologic ages in the Patagonian Andes can be used as a direct record of tectonic uplift (15, 20, 21), yet it is known that cooling ages record exhumation, not uplift (25). Surface processes, such as fluvial incision, landsliding, and glacial abrasion, control erosional exhumation. Tectonic processes can play a role in maintaining topography but are generally not important over the short time scale and modest amounts of erosion documented in our study.

We sampled four glacial deposits in the vicinity of Lago Buenos Aires and compiled published AHe bedrock ages (Fig. 1). The two oldest deposits are interbedded within a sequence of tabular basaltic flows that underlie the Meseta del Lago Buenos Aires. The oldest, herein referred to as the “Mercer” till (informal name), is exposed on the northwest margin of the Meseta del Lago Buenos Aires and is bracketed by basalt flows with ages of 4.7 and 7.3 Ma (26, 27). The “Guivel” till [informal name; 3.3 ± 0.1 Ma (28)] crops out on the southern margin of the meseta. Cobbles from these deposits were collected well below (>10 m) the overlying basalt flow and with no evidence of nearby feeder dikes or sills, which ensures that the AHe ages from these samples were not thermally reset.

The next youngest deposit is Telken VII (1.016 ± 0.01 Ma), which is exposed as a till-covered hill and marks the largest glacial advance in the area, known as the Great Patagonian Glaciation (29). The youngest glacial deposit is Fenix I (18.5 ± 0.4 ka), exposed as a sharp-crested moraine with striated and glacially polished cobbles (30, 31). Note that the four deposits span approximately 75 km north to south and were therefore likely derived from similar source areas within the core of the range (see Fig. 1). Examples of the sampled till and moraine morphology are shown in fig. S1. A collection of 51 published bedrock AHe ages (table S1; locations shown in Fig. 1) was used here to estimate the distribution of AHe lag times for the modern landscape and extend the lag time record to the present day.

To better resolve the depositional age of the Mercer till, we use a deep-ocean temperature record (Fig. 2) (32), which is based on the global benthic foraminifera δ18O record and corrected for the associated evolution of polar ice volumes. Deep-ocean temperature is a record of time-averaged, high-latitude, sea-surface temperature (33) and provides a useful proxy for cool and warm events in our study area. The three well-dated glacial deposits coincide with extreme cold events in this record. We therefore infer that the Mercer deposit was also associated with an extreme cold event. Given the bracketing ages provided by the basalt flows, the likely event would have been at ~5.7 Ma ago (Fig. 2).

Fig. 2 Regional climate record and deposit ages.

The black curve is the deep-ocean water temperature from an ice volume–corrected model based on the global benthic foraminifera δ18O record (32). The colored rectangles indicate the depositional age constraints on the sampled glacial deposits. The circles indicate intersection of the lowest temperature associated with the depositional age range for each deposit.

Figure 1 shows the modern water divide for this portion of the Andes as well as the ice divide and ice flow lines associated with the LGM ice sheet, as determined from (34). With the modern topography, deeply incised glacial valleys allow westward drainage across the full width of the range, far to the east of the highest peaks. The older sampled glacial deposits may have been deposited in association with a more western water divide. However, the water divide likely reached its present position by the time that the younger two sample locations were deposited. Therefore, cobbles in the younger deposits could only have been transported from their western sources during glaciations. This conclusion likely holds for all four sampling locations, as each is a glacial deposit. The ice divide depends on the glacier shape and size, but we use the LGM ice divide as an approximation for the location of older ice divides. On this basis, we estimate that the glacial catchments for our sample locations were ~10,000 km2 for the two older deposits and 5000 km2 for the younger two. The bedrock sample locations straddle the ice divide and cover an area of 6500 km2.


Lag time results for all samples are shown in Fig. 3 (see Materials and Methods for details). Vertical red ticks show AHe ages for each cobble and modern bedrock sample, and the blue curves show probability density distributions for each glacial deposit, as estimated by the kernel density method (35). Figure 3 and Table 1 include estimates of the mean and first quartile lag times for each of the distributions, along with their 2 SE uncertainties. We use the harmonic mean for the detrital samples, given that they are sampled by yield and are therefore affected by local variations in erosion rates. The arithmetic mean is used for the bedrock samples given that they are sampled in space (see Materials and Methods for details). In addition, note that the means include all lag times in the distribution, including some lag times that are greater than the range shown in the plots. Table S2 and fig. S2 provide a full report of all AHe grain ages.

Fig. 3 AHe lag times for samples in this study.

Vertical red ticks indicate the lag times of individual cobbles reported here and published modern bedrock samples (8, 15, 20). Blue curves in the bottom four panels are lag time probability density plots (35), as estimated by AHe minimum ages for each cobble. (Top) Solid blue lines represent predicted probability density curves for modern bedrock using Leones and Nef age-elevation relationships (AERs) in (15); dashed blue lines represent predicted probability density curves for modern bedrock using DES and LL AERs in (20) (see Results and Materials and Methods for details; AER locations are shown in Fig. 1). Vertical dashed lines indicate relevant mean values, and black arrows indicate first quartile values for the lag time distributions. The uncertainties are ±2 SE, and the uncertainty for first quartile values is determined numerically using a bootstrap method.

Table 1 Mean and first quartile lag time and erosion rates.

View this table:

The lag time distributions show a simple evolution with decreasing age. This pattern suggests that, before the onset of alpine glaciation (~6 Ma ago), rates of erosion were slow. Erosion rates increased by 3.3 Ma ago and then returned toward slow at 1.02 Ma ago, 18.5 ka ago, and present day. The oldest deposit (Mercer) has a mean lag time of 24 Ma and a first quartile lag time of 12 Ma. In contrast, the distribution for the next deposit, Guivel, has mean and first quartile lag times of 5 and 0.5 Ma, respectively. The Telken VII distribution has mean and first quartile lag times of 10 and 3 Ma, and the distribution from Fenix I has mean and first quartile lag times of 8 and 4.5 Ma. The distribution of the bedrock samples has mean and first quartile lag times of 6 and 3.5 Ma. The bedrock sample mean is significantly lower than that for the Fenix I deposit, likely because bedrock samples tend to be collected at lower elevations where access is easiest but where AHe cooling ages are youngest.

We use four published age-elevation transects of modern bedrock samples (15, 20) to provide a more direct comparison with the detrital samples. These four transects provide local estimates of the dependence of AHe cooling ages on modern elevation. We used a digital elevation dataset paired with sub-icefield bedrock elevation information (36) to extract elevation data for the portion of each ice catchment underlain by batholith (see Fig. 1). Similar to (37), we multiplied the resulting hypsometry by each age-elevation relationship and then smoothed and normalized it to generate probability density curves of AHe lag times for these regions, as shown in the top of Fig. 3. In comparison, the individual bedrock samples and the cobble samples from the Fenix I deposit show a larger range, both toward smaller and larger lag times, than indicated by these estimated modern bedrock probability density curves. This result suggests that there is a significant spatial variation in erosion rates across the study region.

The black line in Fig. 4 (and in fig. S3 at a more optimal scale) shows the evolution of the spatially averaged erosion rates through time. This curve was estimated by first converting the lag time estimates into an average erosion rate for the age interval covered by that lag time. This conversion was done using a modified version of the age2edot program (38), which finds a simultaneous solution for the thermal structure of the upper crust and the erosion rate at the surface as a function of the observed lag time and the thermally sensitive diffusion properties of the AHe thermochronometer (see Materials and Methods for details). Sensitivity testing indicates that the erosion rates vary by about ±15% over the plausible range of the thermal and diffusive parameters used in the age2edot calculation. Note, however, that this source of uncertainty would shift the entire curve. That is, the relative variations in the erosion rates shown in Fig. 4 come from the lag time data alone and not from the age2edot calculation. Figure S4 shows the erosion rates estimated for all of the cobble lag time data as a function of geologic age. The spatially averaged erosion rate curve (Fig. 4 and fig. S3) was determined by averaging the erosion rate estimates from fig. S4 at 0.5-Ma steps along the age axis. The spatially averaged erosion rate curve includes a fair amount of smoothing due to the fact that the individual erosion rate estimates are averages over the duration of the lag time interval. The smoothing decreases with decreasing lag times and increasing erosion rate. As a result, the temporal resolution of this erosion rate curve tends to increase with decreasing age.

Fig. 4 Summary of the evolution of erosion rates in the source region sampled by the AHe cobble ages.

The black curve shows the evolution of the averaged erosion rate within the source region (fig. S3 shows a larger plot of this curve). The color bars show the fastest rates, as represented by the first quartile value of the erosion rate distribution for each deposit. The horizontal extent of each color bar shows the lag interval (AHe closure to deposition) for each of these “fastest erosion” estimates. The gray curve shows, in schematic fashion, the evolution of the fastest eroding part of the source region. The right-hand axis shows the correspondence between lag time and erosion rate, as determined from the age2edot program (see fig. S6A).

The spatially averaged erosion rate curve (Fig. 4 and fig. S3) shows that, before the onset of alpine glaciation at ~6 Ma ago, the spatially averaged erosion rate for the source region was ~0.06 to 0.1 km/Ma. After the start of glaciations, this erosion rate increased to a steady value of ~0.23 km/Ma. This estimate for the mean erosion rate over the past 6 Ma matches well with previous estimates (0.2 to 0.3 km/Ma) for this time interval, as determined from regional studies of bedrock cooling ages (8, 11).

To gain further insight into the variable rates of erosion across the landscape, we plot the evolution of the fast-eroded cobbles from the source region. For this purpose, we use the value at the first quartile for each lag time distribution (arrows in Fig. 3). Figure 4 shows colored bars that indicate the age interval for the first quartile lag time for each deposit. The gray curve shows an interpretation of the evolution of this “fastest-eroding” component: These estimates begin with a ~0.2 km/Ma erosion rate before the onset of alpine glaciation, increase to ~3 km/Ma at 3.3 Ma ago, and then decline to <0.8 km/Ma between 3.3 and 1.02 Ma ago and to 0.5 km/Ma between 1.02 Ma ago and 18.5 ka ago. The modern bedrock samples indicate a similar spatially averaged erosion rate of 0.5 km/Ma for the first quartile lag time (Table 1). These data indicate a nearly 15-fold increase in the fastest-eroding component with the onset of glaciation, followed by a return to rates similar to those before the onset of glaciation.


This study provides a record of the evolution of erosion rates during the transition from fluvial to glacial conditions in the Patagonian Andes. From these observations, we infer that glaciation did not result in uniformly fast erosion but rather varied spatially across the landscape. We infer that early glaciations would have deepened the preexisting, fluvial, low-concavity valley profiles. Headward propagation of valley incision along longitudinal profiles is known to cause transient increases in local elevation gradients and relief and will lead to ice sliding velocities driving higher incision rates (7). However, as erosion propagates headward, the area of the ice catchment above the ELA is reduced, resulting in slower subsequent glacial erosion (57). Our conclusions are in agreement with (16), which used 4He/3He apatite data and thermokinetic modeling to estimate the formation timing of the exceptionally deep valleys in the Central Patagonian Andes and concludes that valley incision probably occurred shortly after the onset of glaciations in the Patagonian Andes. Their study (16) indicates incision sometime between 10 and 6 Ma ago, while the results of this study indicate sometime between 6 and 3 Ma ago.

If our interpretation is correct, then one might expect the initial pulse of fast incision to be followed by decay to rates that balance with rates of rock uplift. In geomorphic systems, a return to steady state commonly occurs in an exponential fashion. Fitting an exponential to the first quartile erosion rates over the past ~4 Ma (Fig. 4, solid gray curve) indicates an exponential time constant of ~2 Ma. If correct, this relationship would predict that 95% progress to steady state would take ~6 Ma (three times the exponential time constant). As a result, we might expect that the Patagonian Andes are still in transition to a steady state. In comparison, other midlatitude mountain ranges with more recent glacial onset [e.g., the Alps at 2.6 Ma (39)] may be a few millions of years from balance between erosion and rock uplift. Our finding is consistent with recent modeling results that place landscape response to glaciation between a few tens of thousands and a few millions of years (40).

The fastest rates of erosion in our study, ~3 km/Ma, are among the highest rates of glacial erosion observed in the geologic record using AHe thermochronometry (4, 41). Comparably high rates and magnitudes of glacier valley incision occurred in both the Northern and Southern Hemispheres but during the Pleistocene: British Columbia, Canada (51°N), >5 km/Ma at ~1.8 Ma ago (42); Rhône Valley, Switzerland (46°N), 1 to 1.5 km/Ma at ~1 Ma ago (43); and Fiordland, New Zealand (45°S), ~5 km/Ma between ~2 and 1 Ma ago (7). Rapid incision occurred earlier at higher latitudes in both hemispheres: An increase in glacial erosion rates ~200 km south of this study area occurred between 10 and 5 Ma ago (16). The direct relationship between latitude and glacial onset appears to broadly hold in the Southern Hemisphere, as shown in the Andes and continuing onto the Antarctic Peninsula, where glacial onset was recorded at 33.5 Ma ago (44). Study of this long-lived, northward transition to icy conditions over the late Cenozoic helps develop an understanding of the complex interactions between climate, erosion, and tectonics both in the Southern Hemisphere and globally.


The dataset presented in this study resolves systematic changes in mountain-scale erosion rates over a ~6-Ma response to glacial conditions in a midlatitude mountain range. We demonstrate that, in one location, glacial erosion takes on a range of rates over space and time and this may not typically be captured by bedrock studies due to sparse sampling. In this study location, a 15-fold increase in the highest (i.e., first quartile) erosion rates during major topographic adjustment, followed by a decrease in these erosion rates, indicates that there is a measurable transient landscape response to the onset of glaciation. The measured erosion rate is not simply due to the presence or absence of actively sliding alpine glaciers but is a function of the relief and shape of the valley profiles over time and the magnitude of the ice flux. The fastest erosion occurs when flowing ice initially appears on a landscape and the transition time to an equilibrated glacial landscape is on the order of 4 to 6 Ma. As a result, assuming comparable conditions, mountainous landscapes that have been more recently glaciated may still be in a phase of incision and topographic adjustment and may require several million more years of periodic glacial conditions to reach balance between erosion and rock uplift.


AHe thermochronometry

From four deposits, we measured AHe ages for 6 to 29 cobbles per deposit. We collected cobble-sized samples (6 to 10 kg each) of granite or granodiorite from the deposits to ensure (i) provenance from the Patagonian batholith and (ii) the occurrence of apatite for AHe dating. The Mercer till yielded only six dated cobbles (because of poor apatite quality), while the other deposits each contained >15 dated cobbles. The uneven sample count might influence our ability to detect the onset of glacial erosion, but we inferred that the lack of young lag times (<~7.5 Ma) in the oldest deposit (Mercer) provides a reasonable sampling of slow erosion rates before the onset of glaciation. That is, the faster erosion associated with glaciation ensures that the onset of glacial erosion should be well represented by the sampled cobbles.

We isolated apatite crystals using conventional mineral separation methods, and individual crystals were selected for suitability for AHe analysis (euhedral crystals, free of visible inclusions). Molar abundances of U, Th, Sm, and 4He were measured using isotope dilution. A total of 206 crystals from 72 cobbles were analyzed at the University of California Santa Cruz and the Berkeley Geochronology Center (table S2). Measured AHe ages can yield variance between aliquots that exceeds analytical uncertainty, primarily because of undetected inclusions of U- and Th-bearing minerals, variations in crystal size, and zonation of the parent nuclides. These discordance problems generally introduce an old-side bias to the measured AHe age. When possible, we measured replicate AHe ages to screen out older discordant ages. On average, we measured three replicate ages per cobble: Some cobbles have up to six replicate ages, and 22% of the cobbles have only one AHe age. The cobbles with replicated ages indicate that discordance was present in only 1 of 5 cobbles; we therefore believe that issues related to discordance are limited to only about 4% of our cobble ages. Furthermore, AHe data from previous studies in the same region [e.g., (15, 16)] show little discordance in grain ages. For those cobbles with multiple replicate ages, we use a maximum likelihood method to estimate a minimum age (45), defined as the age of the youngest fraction of grain ages that are statistically concordant as defined by analytical errors. This calculation assumes that the replicate AHe ages are a two-component mixture, consisting of young grain ages free from biasing effects and older grain ages that are randomly affected by biasing effects.

Helium diffusion in apatite is estimated using numerical integration with time-invariant diffusion parameters (46). We have not accounted for the potential influence of radiation damage on He diffusivity (more detail on this point below).

Probability density estimates

The probability density curves in Fig. 3 and fig. S2 were estimated using the kernel density method (35) with the kernel set to 2.4*SE, where SE is the standard error for the represented discrete observations (e.g., grain ages and minimum ages). The recommended kernel size is 0.6*SE value (35), but we selected a larger size for our study to emphasize the general features of the probability density distributions. A probability density curve should integrate to unit probability mass. In keeping with this constraint, each probability curve was normalized so that it has the same integrated area as the others.

Figure S5 compares two different estimates of the probability density curves. The blue curve, which is the one shown in Fig. 3, was based on the minimum ages from the replicate ages for each of the cobble samples. The gray curve was determined using all of the AHe replicate ages with a weighting applied to ensure unit weight for each cobble sample. The comparison shows that these two methods give similar results.

Converting lag time to erosion rate

We used the age2edot program to convert lag times (i.e., the AHe age minus the deposition age) into erosion rates (38). The age2edot model does not use a prescribed age-elevation relationship but instead determines an average erosion rate as a function of the cooling age, the closure properties of the AHe thermochronometer, and the one-dimensional thermal structure of the upper 30 km of Earth’s crust, all of which are sensitive to the erosion rate. For this calculation, we assumed that the depth to the closure isotherm can be treated in a quasi-steady fashion, which is consistent with the fast response time (<~1 Ma) of the AHe system to changes in erosion rate (13). This calculation also accounts for the thermal-magmatic structure of the region, as guided by a recent magmatic arc model (47), which postulates that the thermal structure of the arc is strongly controlled by the emplacement of mantle-derived basaltic melt into the lower crust of the arc. Arc volcanoes occur on a ~100-km spacing across the region, but models of subduction magmatism indicate that there is likely much more widespread melt within the lower crust at depths greater than 20 to 30 km (13, 47). This effect was accounted for in the age2edot calculation. One might anticipate that the shallow crust might be strongly influenced by feeder dikes associated with volcanoes. Thermal analysis, however, shows that, in the shallow crust, potential resetting around a feeder dike would be limited to a region on the order of the thickness of the dike (48). We therefore conclude that thermal resetting is relatively rare in the shallow crust and our AHe lag time data are primarily a result of cooling during erosion and can be used to estimate erosion rates.

The age2edot program represents the thermal structure of the upper crust using an infinite layer with fixed boundary temperatures. The thermal diffusivity and thermal conductivity of crust are set using a new relationship (49) that accounts for the temperature sensitivity of these properties in typical crustal rocks. The volumetric heat production was set to a uniform value of 2 μW/m3, which is the average for the Sierra Nevada batholith (50). The upper boundary corresponds to the local mean elevation of the land surface (~1000 m) and was set to a temperature equal to the long-term atmospheric temperature at that elevation (~0°C). The lower boundary was set at 20-km depth and ~800°C, the approximate solidus temperature for a granodioritic crust at 20 km (48). The crust beneath the source region is likely mainly composed of Patagonian batholith, hence the choice of the granodioritic solidus temperature. For comparison, the more pelitic composition of the schist belts exposed on the flanks of the range would have a solidus temperature of ~700°C at 20-km depth (48). Erosion is represented by a vertical velocity through the layer. The presence of melt below the basal boundary ensures that the material crossing that boundary is always at the solidus temperature. We solved for the lag time of the AHe closure surface as a function of the thermal properties of the crust, the diffusion properties of the AHe system, and the erosion rate.

The depth to the top of the lower crustal melt zone is controlled by the flux rate of the mantle-derived melt, which, in turn, is controlled by the subduction velocity and the corner flow velocity within the supra-slab mantle (51). The top of the melt zone, which marks the shallowest region in the crust with coexisting melt and crust, should remain at a fixed temperature, as defined by the selected solidus curve. The depth of this boundary is mainly controlled by conductive heat transport through the upper crust. Thus, surface erosion will cause the crust to move toward the surface, but the top of the melt zone should remain at a steady depth. This situation was correctly represented in the age2edot model by a fixed basal boundary condition, which ensures that as the material rises through that boundary, it enters into the model domain at the temperature set for the boundary.

Figure S6A shows that the temperature and depth of the basal boundary condition have little influence on the predicted relationship between erosion rate and AHe cooling age. Figure S6B shows the dependence of closure depth on erosion rate. The estimated closure depth for our cobble samples is between 2.4 km for long lag times and slow erosion and 1.1 km for short lag times and fast erosion. To verify the validity of a quasi-steady state solution for this setting, we run a full transient calculation of the temperature history of a sample and the evolution of the sample AHe age (fig. S7). The steadiness of the closure depth was measured by the velocity ratio (vertical axis of fig. S7), defined as the ratio of the velocity of the closure surface relative to the vertical material velocity (equal to the erosion rate). The velocity ratio shows high values following the onset of erosion, which indicates unsteady migration of the closure surface, but the ratio drops back down to nearly zero within 2 Ma.

We estimate helium diffusion in apatite using time-invariant diffusion parameters (46). Radiation damage is unlikely to be a significant source of variance in diffusion parameters given that most AHe ages are relatively young, and all samples are likely to have only experienced cooling (i.e., no reheating) through geologic time. To test this, we used the cooling paths estimated from the first quartile lag times, the alpha damage annealing model (ADAM) (52), and minimum and maximum measured U and Th concentrations (table S2). In the case producing the largest difference (i.e., the Mercer deposit with the slowest apparent erosion rates and the lowest U and Th concentrations; table S2), the corrected AHe age would be no more than 30% less relative to the nominal age. This bias would preferentially influence our estimates of the lowest erosion rates, for example, shifting them upward from 0.2 to 0.3 km/Ma. In contrast, our estimated cooling paths for the young Guivel cobble ages would yield a revised erosion rate of 3.2 km/Ma instead of 3 km/Ma. Because these differences are relatively small, we concluded that our use of time-invariant diffusion parameters from (46) is a sufficient approximation in this setting and eliminates additional computational expense of time-variant diffusivity for each of the 206 crystals. Comparable assumptions were made for bedrock ages measured in the same region (15).

Estimating spatially averaged erosion rates

Each of the cobbles in this study represents an individual bedrock sample from an upland granitic source region “collected” by a glacier in the geologic past. As a result, areas with faster erosion will yield more cobble samples per unit area of the source region than those areas with slower erosion. We used a method that accounts for this bias and provides spatially averaged erosion rates for the source region. Consider a randomly sampled distribution of erosion rates, ei, where i = 1 to n, that are determined, in some fashion, from sediment materials eroded from the source region. The sample mean is defined bye¯=wieiwi(1)where wi are weights for each sample. The conventional mean uses uniform weights for the samples (wi = 1). For our problem, this estimate would give mean erosion rate as sampled by the sediment yield. Our objective is to estimate the mean erosion rate as sampled by area. To do so, we set the weights as wi = 1/ei, which removes the bias due to spatially varying erosion rates. Substitution into Eq. 1 givese¯=((1/ei)n)1(2)

This derivation shows that the mean erosion rate by area is simply the reciprocal of the mean of the reciprocal erosion rate. This estimator is called the harmonic mean and is known to be useful when averaging certain kinds of rate measurements, including the use of detrital quartz 10Be measurements to estimate the spatially averaged erosion rate of a river catchment. In the same way, we used the harmonic mean to estimate spatially averaged erosion rates from our cobble lag times.


Supplementary material for this article is available at

Fig. S1. Example photos of till and moraine deposit morphology.

Fig. S2. Lag time plot showing all AHe ages (including all replicate ages).

Fig. S3. Average erosion rate curve from Fig. 4.

Fig. S4. Color bars showing lag intervals for all AHe cobble ages in our study, plotted separately for each deposit.

Fig. S5. Simplified version of Fig. 3.

Fig. S6. Age2edot output.

Fig. S7. Evolution of the closure depth in response to an instantaneous increase in erosion rate.

Table S1. Published bedrock ages shown in Fig. 3.

Table S2. Crystal AHe measurements and ages.

Data S1. Google Earth file (.kml) of sample names and locations.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We thank A. Bryk for assistance with the hypsometry datasets. We also greatly appreciate the constructive feedback of three anonymous reviewers and B. Schoene’s useful feedback and thoughtful handling of the manuscript. Funding: This work was supported by the Yale College Dean’s, Von Damm, and Stiles Mellon Forum Research grants, and NSF Graduate Research Fellowship DGE 1752814 (to C.D.W.), NSF EAR-1650313 (to M.T.B.), and the Ann and Gordon Getty Foundation and the Esper Larsen Research Grant (to D.L.S.). Author contributions: All authors conceptualized and developed the project and collected, processed, and analyzed samples. C.D.W., M.T.B., and D.L.S. interpreted data, developed the thermal model, and wrote the manuscript with coauthor support. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. The age2edot code is available by request from M.T.B. (mark.brandon{at} Additional data related to this paper may be requested from the authors.
View Abstract

Stay Connected to Science Advances

Navigate This Article