# Milankovitch-paced erosion in the southern Central Andes

G. Burch Fisher, Lisa V. Luna, William H. Amidon, Douglas W. Burbank, Bas de Boer, Lennert B. Stap, Bodo Bookhagen, Vincent Godard, Michael E. Oskin, Ricardo N. Alonso, et al.

► To cite this version:

G. Burch Fisher, Lisa V. Luna, William H. Amidon, Douglas W. Burbank, Bas de Boer, et al. Milankovitch-paced erosion in the southern Central Andes. Nature Communications, 2023, 14 (1), pp.424. 10.1038/s41467-023-36022-0. insu-04131497

HAL Id: insu-04131497

https://insu.hal.science/insu-04131497

Submitted on 16 Jun 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
<figure><img src="image_2.png" /></figure>

# Milankovitch-paced erosion in the southern Central Andes

Received: 16 November 2020

Accepted: 11 January 2023

Published online: 26 January 2023

<figure><img src="image_3.png" /></figure>

Check for updates

G. Burch Fisher Ⓓ1,2 Ⓒ, Lisa V. Luna Ⓓ3, William H. Amidon Ⓓ4, Douglas W. Burbank Ⓓ2, Bas de Boer Ⓓ5, Lennert B. Stap Ⓓ6, Bodo Bookhagen Ⓓ7, Vincent Godard8,9, Michael E. Oskin Ⓓ10, Ricardo N. Alonso Ⓓ11, Erik Tuenter12 & Lucas J. Lourens Ⓓ13

It has long been hypothesized that climate can modify both the pattern and magnitude of erosion in mountainous landscapes, thereby controlling morphology, rates of deformation, and potentially modulating global carbon and nutrient cycles through weathering feedbacks. Although conceptually appealing, geologic evidence for a direct climatic control on erosion has remained ambiguous owing to a lack of high-resolution, long-term terrestrial records and suitable field sites. Here we provide direct terrestrial field evidence for long-term synchrony between erosion rates and Milankovitch-driven, 400-kyr eccentricity cycles using a Plio-Pleistocene cosmogenic radionuclide paleo-erosion rate record from the southern Central Andes. The observed climate-erosion coupling across multiple orbital cycles, when combined with results from the intermediate complexity climate model CLIMBER-2, are consistent with the hypothesis that relatively modest fluctuations in precipitation can cause synchronous and nonlinear responses in erosion rates as landscapes adjust to ever-evolving hydrologic boundary conditions imposed by oscillating climate regimes.

The hypothesis that climate exerts a first-order control on landscape form and function through physical and chemical erosion dates back more than a century<fnref n="1" />. Nonetheless, quantifying the influence of climate on erosion rate and morphology in any given landscape has proven difficult. Whereas numerical models provide clear rationales for how precipitation foci and oscillations may impact deformation and sediment-flux regimes in orogenic systems²⁻⁵, empirical evidence has remained equivocal with convincing evidence both for and against strong coupling between climate, erosion, and tectonics⁶. This ambiguity partially derives from temporal and spatial mismatches in both

the datasets and proxies used to infer such linkages. For example, many interpretations rely on comparing present-day precipitation patterns with erosion and/or tectonic rates measured over 10²⁻ to 10⁶⁻ year timescales. Further complications arise in deciphering the primary driver of erosion in orogenic systems due to the common spatial coincidence between precipitation, topographic relief maxima, and deformational foci as a consequence of orographic effects⁷.

Acknowledging these complications, researchers have sought to exploit sedimentary archives in order to extract temporal records of erosion rates for direct comparison to known oscillations and

<figure><img src="image_4.png" /></figure>

e-mail: gbf@ucsb.edu

<footnote n="1">¹Jackson School of Geosciences, University of Texas at Austin, Austin, TX 78712, USA. ²Earth Research Institute, University of California, Santa Barbara, CA 93106, USA. ³Institute of Environmental Science and Geography, University of Potsdam, Potsdam, Germany. ⁴Department of Earth and Climate Sciences, Middlebury College, Middlebury, VT 05753, USA. ⁵Earth and Climate Cluster, Faculty of Science, Vrije Universiteit Amsterdam, Amsterdam, the Netherlands. ⁶Institute for Marine and Atmospheric Research, Utrecht University, Utrecht, the Netherlands. ⁷Institute of Geoscience, University of Potsdam, Potsdam, Germany. ⁸Aix-Marseille Univ., CNRS, IRD, INRAE, CEREGE, Aix-en-Provence, France. ⁹Institut Universitaire de France, Paris, France. ¹⁰Department of Earth and Planetary Sciences, University of California, Davis, CA, USA. ¹¹Universidad Nacional de Salta, Salta, Argentina. ¹²Royal Netherlands Meteorological Institute (KNMI), De Bilt, the Netherlands. ¹³Department of Earth Sciences, Faculty of Geosciences, Utrecht University, Utrecht, the Netherlands.</footnote>
excursions in climate⁸⁻²⁵. Unfortunately, these studies are rarely straightforward, and commonly rely on less-than-ideal field exposures, temporal controls, and source-area and cosmogenic radionuclide systematic constraints. In addition, many studies lack the sampling frequencies necessary to unambiguously link paleo-erosion-rate trends with climatic changes. For example, in the Tien Shan region, the addition of new sediment archives resulted in the revised conclusion that tectonics was predominantly responsible for controlling erosion in the region during the Plio-Pleistocene¹⁶ rather than the onset of Quaternary glaciation as previously proposed¹¹. Likewise, in the south-central Andes of Argentina, two studies developed paleo-erosion rate records along the same section¹⁴,¹⁸. Denser sampling and improved temporal controls used by the later study¹⁸ allowed identification of two pronounced reductions in paleo-erosion rates not identified by the earlier study and consistent with regional and marine climate proxies.

Despite the complexities involved with extracting high-quality erosion-rate data from terrestrial sediment records, such data have the potential to provide unique insights into how landscapes respond to diverse climatic and tectonic perturbations across a range of forcing frequencies. To date, considerable research effort has focused on identifying Milankovitch cycles, commonly referred to as precession (-23 kyr cyclicity), obliquity (-41 kyr), and eccentricity (-100 kyr and -400 kyr), within terrestrial sedimentary records. These orbital periodicities, which change the latitudinal distribution of solar insolation by season, have been widely acknowledged across marine and atmospheric records as key global climate drivers. Yet, linkage of sediment

fluxes or erosion rates to these climate drivers over multiple cycles has remained largely theoretical, especially in active tectonic settings²,⁴,⁵. Thus, any direct linkage between climate oscillations and sediment flux from the continents could have far-reaching implications for landscape evolution theory, chemical weathering regimes, and ecosystem and evolutionary trajectories²⁶,²⁷.

To address the hypothesis that long-term variations in climate can drive consequent responses in erosion rate, we present 49 ¹⁰Be and 23 ²⁶Al cosmogenic radionuclide-derived paleo-erosion rates sampled along the Río Iruya canyon in northwest Argentina spanning the period from 6-1.8 Ma (Fig. 1 and Supplementary Figs. 1 to 3). This approach to calculating paleo-erosion rates is similar to the well-established practice of calculating erosion rates (10²–10⁴-year timescales) from present-day cosmogenic radionuclide concentrations in river sediments, which quantifies the combined signal (also referred to as denudation) from both physical erosion and chemical weathering within the contributing catchment²⁸. The technique relies on the fact that cosmogenically produced nucleons attenuate exponentially with depth through the upper few meters of rock and regolith, where they spall O and Si atoms in quartz into ¹⁰Be (half-life = 1.387 Myr)²⁹ and ²⁶Al (half-life = 0.705 Myr)³⁰, respectively, which accumulate over time. This behavior permits radionuclide concentrations in surface samples to be converted into erosion rates, whereby the concentration is inversely proportional to the rate³¹. Assuming that the fluvial system reliably integrates sediment from throughout the catchment, a sample of river sediment provides an accurate estimate of the mean catchment erosion rate²⁸. If these sediments are then rapidly buried and progressively

<figure><img src="image_5.png" /><figcaption>**Fig. 1 | Geographic and climatic setting of the Río Iruya study area.** a Major modern South American climate and topographic components including modern mean annual rainfall derived from the Tropical Rainfall Monitoring Mission (TRMM) 2B31 dataset⁷, the South American low-level jet (SALL), and the internally-drained, high-elevation Altiplano and Puna plateaus (see Supplementary Fig. 1 for detailed climatic and topographic characteristics of the study area). The white star locates the study site in a and b. b Plot of modern peak mean annual rainfall (±2σ envelope) (blue) extracted from 115, 50-km-wide by 1000-km-long, orogen-perpendicular swath profiles along the Andes (see⁷) indicating a transitional zone (gray) between -18–25°S latitude. Linearly-interpolated mean annual rainfall derived from CLIMBER-2 model results (magenta), extracted from a longitudinal swath spanning 95–115°W at 10° latitude intervals from 10°N to 35°S in the model (see Methods for discussion</figcaption></figure>

of the geographic offset of the model), shows maximum precipitation values from a precession maximum at 2.670 Ma and minimum values from a precession minima at 2.487 Ma. Magenta percentages indicate the maximum percent increase in mean annual precipitation between the plotted precession minimum and maximum values, with the mean differences in parentheses. c Overview of the Río Iruya watershed geography and Eastern Cordillera topography with predominant thrust faults labeled. Locations of modern cosmogenic radionuclide (CRN) erosion-rate samples (IR-3M and IR-6M) are shown as white stars. d Close-up of the Río Iruya canyon section showing the pre- and post-diversion Río Iruya flow directions (hashed and solid white arrows, respectively), CRN paleo-erosion rate sample locations (white circles) with select stratigraphic ages (white numbers), and total estimated incision in the canyon over the -150 years post diversion (blue-red color gradient).
<figure><img src="image_6.png" /></figure>

**Fig. 2 |** *Advantageous characteristics of the Río Iruya canyon section for cosmogenic radionuclide paleo-erosion rate analysis.* a Stratigraphy, volcanic ashes with high-resolution ages (red stars), magnetostratigraphic (PMAG) sample sites, magnetic polarity (black = normal polarity; white = reversed polarity), and overall chronologic constraints for the Río Iruya section³²,³⁴ (see Supplementary Fig. 2 and Supplementary Dataset 1). b Calculated sediment-accumulation rates for the Río Iruya section¹² (mean in black and maximum and minimum envelope in gray). Dashed red line is the overall mean accumulation rate for the section. Rate uncertainties are greatest for brief magnetochrons where the precise positions of polarity boundaries are not well resolved relative to the length of the magnetochron (see Supplementary Dataset 1). Mean exhumation rates for cosmogenic radionuclide samples (blue dots and trend line) are based on the sample depth below the pre-diversion river level and assuming 150 years since the diversion. The

<figure><img src="image_7.png" /></figure>

associated error envelope assumes a+/- 15 m error on the reconstructed pre-diversion surface elevation. c An error-weighted two-nuclide test for complex transport and burial history showing minimal diversion from the expected decay-corrected ²⁶Al:¹⁰Be ratio (7:1 in dashed blue line) for all ²⁶Al samples with detectable concentrations and AMS errors less than 35%, minus the outlier at 2.49 Ma (95% confidence intervals are shown by the gray envelope; n = 8; mean ratio = 7.4:1). Large ²⁶Al errors are related to dilution of the ²⁶Al/²⁷Al ratio by considerable inherent Al in the sample quartz. Note ²⁶Al and ¹⁰Be errors are based on the paleo-erosion concentration component, N<sub>E</sub> (i.e., corrected for exhumation, decay, and burial components), with 2σ errors reported (¹⁰Be errors fall within the markers). d Plot of ²⁶Al errors with sediment age and concentration showing the non-systematic distribution of high AMS error samples (>75%).

preserved over millions of years, a time-resolved record of paleo-erosion rates can be reconstructed by measuring ¹⁰Be or ²⁶Al in quartz from sediment in a well-dated depositional sequence¹⁰,¹¹,¹⁴,¹⁷,¹⁸.

The Río Iruya canyon, located on the eastern flank of the northern Argentinian Andes and with a watershed bordering the internally-drained central Andean Plateau (Puna Plateau), provides an ideal location to apply a cosmogenic radionuclide paleo-erosion rate technique for three reasons (Fig. 1). First, the close proximity of the depositional foreland basin to a steep source area (<100 km) minimizes the likelihood of prolonged sediment storage, which may alter ¹⁰Be concentrations acquired in the target catchment during source-to-sink transport. Second, rapid foreland-basin deposition rates (~1 m/kyr)³² suggest that ¹⁰Be concentrations at deposition have been minimally affected by post-depositional nuclide accumulation during burial (Fig. 2). Lastly, the Río Iruya canyon is the result of a large-scale flood-diversion project completed in 1865 that led to the rerouting and shortening of the mainstem Río Iruya (Fig. 1c). Subsequent vertical incision of ~100 m has recently exposed >8 km of continuous Miocene-to-Pleistocene-aged foreland sediments known as the Oran Group³²⁻³⁴. The rapid exhumation of the canyon (0.5–1 m/yr) minimizes ¹⁰Be production due to recent exposure and, when combined with a continuous, high-resolution magnetochronology³² and a simple pre-burial transport history (Fig. 2 and Supplementary Fig. 2), permits extraction of a robust record of erosion-rate variability that minimizes confounding and commonly unconstrained sources of post-depositional ¹⁰Be and ²⁶Al.

The Río Iruya watershed drains the eastern edge of the Puna Plateau, which is defined by a series of mountain ranges known as the Eastern Cordillera. In the Iruya watershed, the Eastern Cordillera exhibits nearly 4000 m of relief, plunging steeply downward into the adjacent foreland basin, which is defined by a N-S striking fold and thrust belt (Fig. 1). The Eastern Cordillera is likely growing over a footwall ramp, which extends eastward as a detachment surface with splay faults connecting upwards to several foreland faults and folds³³,³⁵. In the study region, uplift of the eastern Puna Plateau was likely underway by the Oligocene³⁶,³⁷, creating significant relief, deep basins and orographic barriers by ca. 10–8 Ma³⁵,³⁸,³⁹. Facies migration within regional foreland basin sequences suggests the onset of thrusting in the Eastern Cordillera at roughly 15 Ma¹⁰ consistent with the onset of deposition of the Oran Group between ca. 13.5–8 Ma³³,³⁴,⁴¹,⁴². Uplift and erosion of the Eastern Cordillera has continued until today, resulting in progressive incision of the Río Iruya into the metamorphic core of the Eastern Cordillera and with its headwaters tapping into basins of the eastern Puna Plateau (Fig. 3)⁴³. A detrital fission-track study along the Río Iruya suggests that areas of peak incision (e.g., the Río Astillero) have experienced roughly 6 km of exhumation over the last 10 Myr, yielding a long-term average rate of roughly 0.6 mm/yr⁴³.

Amidon et al.³² have argued that the Río Iruya watershed experienced major provenance shifts at roughly 4.0 Ma and 2.4 Ma, each of which corresponds temporally with changes in the spatial pattern of thrust faulting and surface uplift³². The first provenance shift is recorded by the addition of Neogene zircons as well as a very restricted range of
<figure><img src="image_8.png" /></figure>

**Fig. 3** | Schematic illustration of late Miocene-to-Pleistocene tectonics, climate, and deposition in the Río Iruya catchment. Pink rectangle indicates the approximate location of the Río Iruya section of this study. a Latest Miocene to mid-Pliocene interval featuring ongoing uplift of the Eastern Cordillera and early stages of foreland development. Largely transverse drainages dominate the depositional system and Río Iruya watershed may not have been directly connected to the eastern Puna Plateau. b Mid-Pliocene-to-Early Pleistocene period characterized by ongoing tectonic uplift in the eastern Cordillera (Sierra Santa Victoria and Sierra

Tilcarca). The Río Iruya watershed may have been connected to the eastern Puna Plateau based on the strong proportion of Neogene-aged zircons in the sediment. c Hinterland Puna Plateau drainages such as the Quebrada de Humahuaca are fully diverted to the south by sustained range growth and increasing aridification. The Río Iruya may have become disconnected from the Puna Plateau and begun to establish its modern extent based on the near disappearance of Neogene zircons from the sediment.

single-grain quartz chemistries occurring between 4.3 and 3.5 Ma³². The observed increase of Neogene-aged zircons is interpreted to result from the watershed gaining connection to a source of sediments from the volcanic provinces of the eastern Puna Plateau. The restricted quartz chemistry is interpreted as a more localized sediment source. Together, these observations suggest a change in the extent of the Río Iruya

watershed during this period. This shift is roughly contemporaneous with out-of-sequence deformation in the foreland between 4.5 and 4.0 Ma as thrusting stepped backwards onto the Pescado Thrust (Fig. 1)³³,³⁴. Meanwhile, in the Humahuaca region just to the southwest, active thrusting along the eastern edge of the Puna Plateau disconnected the basin from the foreland, rerouting fluvial networks into a more N-S
axial geometry and triggering deposition of the fluvial Tilcara Formation by roughly 4.2 Ma (Fig. 3)⁴⁴. Likewise, the nearby Sierra Alta experienced a pulse of uplift at around 4.3 Ma (Fig. 3)⁴⁵.

The second provenance shift occurred between 2.4 and 2.1 Ma as the proportion of Neogene zircons drops by ~90% and the distribution of single-grain quartz chemistry regains its full complexity³². This provenance change corresponds to an important stratigraphic boundary as the Terciario Subandino sequence is unconformably overlain by the coarse El Simbolar conglomerates³³,³⁴. This provenance shift is likely related to ongoing tectonic deformation in the region, which eventually caused another change in the watershed geometry. For example, foreland deformation stepped out to the easternmost Aguaragüe fault at roughly 2.7 Ma³³. In the Humahuaca region (Fig. 3), increasing aridification caused denudation rates to drop by an order of magnitude from ~4 Ma to 2.7 Ma²¹, and deposition of the fluvial Tilcara formation ceased, giving way to coarse alluvial conglomerates by 2.4 Ma⁴⁴.

Climatically, the modern Río Iruya region has a semi-arid sub-tropical climate in which precipitation is heavily controlled by topography, with mean annual precipitation varying from <0.5 m/yr on the arid Puna Plateau to ~2.5 m/yr on high-relief, east-facing mountain slopes⁷ (Fig. 1, Supplementary Fig. 1). Nearly all precipitation is derived from easterly or northeasterly sources due to the topographic barrier of the Puna Plateau and hyper-arid Atacama Desert immediately to the west. Roughly 80% of precipitation in the study region falls in the Austral summer months between November and February⁴⁶ when easterly winds bring moisture off the Atlantic Ocean, fueling the South American Monsoon System (SAMS) over Brazil⁴⁷. Some of this moisture is channeled down the Andean Front to the Río Iruya watershed via the South American Low Level jet (SALL): a low-level atmospheric flow drawn toward the Chaco high pressure system near 25° S (Fig. 1a). Given its location in the continental interior and the highly seasonal precipitation, precipitation in the Río Iruya region is sensitive to the strength and position of the South American monsoon and the SALL, which are in turn governed by patterns of global oceanic and atmospheric circulation⁴⁸.

Although several studies suggest the SALL system was in place by the mid-Miocene³⁹,⁴⁹, very little is known about the terrestrial climatic conditions in the Andes during the Pliocene-Pleistocene time given the lack of high-resolution lacustrine or glacial records during this period. We thus assess the role of climatic forcing using the low-resolution, intermediate complexity CLIMBER-2 model⁵⁰. The CLIMBER-2 model, although coarse in spatial resolution (latitudinal resolution of 10° and longitudinal resolution of ~51°) and insensitive to regional orographic forcings and climate systems driven by the topography of the Andes (e.g., the SALL), provides one of the only robust modeling tools with which to perform long-term climate system simulations over millions of years. To this end, the model provides a broad understanding of the magnitude and frequency of hemispheric moisture transport to the study region during the Plio-Pleistocene.

In this work, we document the direct erosional response of a tectonically-active southern Central Andes watershed (Río Iruya) to imposed climate oscillations during the Pliocene, as well as to a major climate transition at the Plio-Pleistocene boundary. We propose that the 400-kyr eccentricity-paced erosion-rate signal observed during the Pliocene arises from nonlinear fluvial incision and hillslope erosion thresholds, namely those associated with landsliding in the catchment: a signal resolvable in the sedimentary record because of the erosional sensitivity of the study region to relatively small perturbations in precipitation. We posit that the observed synchrony between erosion rates and climate oscillations in this unique study area is likely restricted to those regions dominated by transitional climate regimes (i.e., insolation-driven monsoons) and where geomorphic response times and processes are sensitive enough to record the imposed climatic oscillations (e.g.,

arid to semi-arid climates). These findings should newly inform ongoing debates concerning the role of climate cyclicity in driving landscape evolution, as well as the extent to which nonlinearities in the geomorphic response of a landscape may destroy, dampen, or enhance climate signals across diverse climate-erosion-tectonic couplings and those recorded in sedimentary archives.

## Results

Paleo-erosion rates spanning approximately 6–1.8 Ma can be broadly grouped into three periods of behavior spanning 6–4 Ma, 4–2.4 Ma, and 2.4–1.8 Ma (Fig. 4 and Fig. 5). Although these intervals are defined by the nature of the erosion-rate record, transitions between them are roughly contemporaneous with previously documented tectonic and sediment provenance shifts at ~4.0 and ~2.4 Ma described above. All paleo-erosion concentrations ($N_E$) and paleo-erosion rates ($R_E$) presented have been corrected for exhumation ($N_X$), decay ($N_D$), and burial components ($N_B$) (see Methods for details). Outside of the decay component, the burial component provided the largest modification of the paleo-erosion concentrations used to infer erosion rates ($N_E$) with values ranging from 1–36% (mean = 10%). Concentration modification associated with exhumation was negligible (<0.05%) as expected given the recent, rapid incision of the Río Iruya canyon.

### Paleo-erosion Rates from 6 Ma to 4 Ma

Prior to ~4 Ma, erosion rates vary widely from ~0.05 to 0.42 mm/yr based on the assumptions of the preferred catchment-uplift and production-rate evolution scenario (Fig. 5, Supplementary Fig. 4). Despite the large number of paleo-erosion rate estimates ($n = 21$), it is difficult to attribute a consistent pattern or frequency to the erosion rates during this time interval, especially given that many of the erosion-rate peaks are defined by single data points. It seems unlikely that this ambiguous erosion-rate pattern is due to uncertainty in stratigraphic age given that the assigned stratigraphic ages are tightly constrained by 8 magnetic reversals spread throughout the period from 4.2 to 5.3 Ma and anchored by four volcanic ash beds (Fig. 2, Supplementary Fig. 2). Although the geometry of the watershed prior to 4 Ma is uncertain, field observations and provenance studies indicate deposition during this period occurred in a larger and more distal depositional system characterized by finer-grained sediments and a poorly-constrained source area (Fig. 3a)³².

### Paleo-erosion Rates from 4 Ma to 2.4 Ma

At ~4 Ma,¹⁰Be concentrations record the onset of a pronounced periodicity in erosion rates from ~4 to 2.4 Ma (Fig. 4 and Fig. 5). This period is characterized by up to three complete paleo-erosion rate cycles, displaying >10-fold rate changes (~0.02 to 0.23 mm/yr) that are largely in-phase with predicted Milankovitch long eccentricity cycles at the 400-kyr frequency (Fig. 4 to Fig. 6), a correlation that is discussed in detail below.

### Paleo-erosion Rates from 2.4 Ma to 1.8 Ma

The final period of the erosion-rate record from ~2.4 to 1.8 Ma departs from the 400-kyr cyclicity and is characterized instead by a large increase in erosion rates, which fluctuate about a mean consistent with the modern rate of ~0.34 mm/yr (Fig. 5).

### CLIMBER-2 Climate Model

Results from the CLIMBER-2 model suggest that, although precipitation oscillations in the region are dominated by the 23-kyr precession and 41-kyr obliquity-frequency bands during this period (Fig. 7), 400-kyr eccentricity cycles strongly modulate changes in amplitude of the 23-kyr precipitation cycles (Fig. 4e, f). When eccentricity is high, for example, precession cycles cause ~20% swings in summer insolation at 23°S, cycles that are damped to ~2% during periods of low eccentricity
<figure><img src="image_9.png" /><figcaption>**Fig. 4** | Cosmogenic radionuclide concentrations and climate datasets. a Global compilation of benthic δ¹⁸O values⁵¹ (VDPB, the Vienna Pee-Dee belemnite standard) and b mean austral summer insolation (Dec.-Mar.) calculated at the 23°S latitude of the Río Iruya study area⁸³. c ¹⁰Be and d ²⁶Al paleo-erosion concentrations (Nₑ - i.e., corrected for exhumation, decay, and burial components) (2σ errors are contained within the black circles for ¹⁰Be) on reversed axes where decreasing concentration is indicative of increasing erosion rate, assuming a constant mean catchment-production rate. e Modeled mean annual precipitation rate for the Río Iruya region (taking the grid box that covers the region of interest, i.e., from 116.4-65°W and 30-20°S) derived from the CLIMBER-2 model with f the Gaussian-filtered 23-kyr (magenta), 41-kyr (green), and 400-kyr (black) CLIMBER-2 precipitation components. Note: the 400-kyr eccentricity curve is shown in light gray and repeated across all panels for reference. N.H. = Northern Hemisphere; DZs = detrital zircons; ENSO = El Niño Southern Oscillation; Obl. = obliquity; Prec. = precession.</figcaption></figure>

(Fig. 4b). These model results are notable because paleo-erosion rates correlate with the amplitude of CLIMBER-2 precipitation oscillations from ~4.0 to 2.4 Ma (Fig. 6), which themselves are largely driven by insolation.

Although spatially coarse, several lines of evidence suggest that the CLIMBER-2 data are, broadly, a reasonable predictor of Plio-Pleistocene climatic variability. For example, strong 400-kyr Milan-kovitch pacing is observed in Antarctic ice volume, sea surface
<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /></figure>

temperature, and global carbon-cycle dynamics during the Pliocene⁵¹,⁵². Likewise, the dominance of the 23-kyr frequency oscillation in the CLIMBER-2 precipitation results is in agreement with late Quaternary hydrologic records across South America linking SAMS precipitation intensity to precessionally-driven summer-insolation maxima⁵³,⁵⁴. Finally, the Plio-Pleistocene CLIMBER-2 model results mimic modern north-south hemispheric precipitation gradients,

indicating adequate sensitivity to latitudinal moisture transport in the model (Fig. 1b and Supplementary Fig. 5).

## Discussion

The most striking result of this study is the correlation between erosion rate and the 400-kyr amplitude modulation of 23-kyr summer insolation oscillations from ~4.0 to 2.4 Ma (Fig. 5a, Fig. 6b). Given the
**Fig. 5** | ¹⁰Be-derived paleo-erosion rates from the Rio Iruya and nearby Quebrada de Humahuaca from 6 Ma to 1.8 Ma. **a** Paleo-erosion rates calculated for the Rio Iruya catchment based on the preferred uplift regime and mean catchment production rate evolution scenario (see Supplementary Fig. 4b and Supplementary Dataset 2). Black circles indicate samples, and the blue piecewise cubic spline is interpolated between sample errors (±2σ envelope). Long-eccentricity cycles (400 kyr) are shown in the background for reference (light gray) with the modern Rio Iruya catchment erosion rate (0.34 mm/yr) shown in dark gray. **b** Paleo-erosion rates (±1σ bars)(red) and c δD values of volcanic glass shards sampled throughout

the adjacent intermontane Quebrada de Humahuaca (purple) (±1σ on age and δD values) (Fig. 3, Supplementary Fig. 3)²¹,⁸⁶. Data from the Humahuaca Basin document sustained surface uplift of the Humahuaca Basin and surrounding ranges in the Eastern Cordillera from -6 Ma to 3.5 Ma resulting in the development of a significant orographic barrier along the Eastern Cordillera by 3.5 Ma. Continued orographic isolation and aridification of the basin followed, driving an order of magnitude reduction in basin denudation rates despite the presence of considerable topography. V-SMOW - Vienna standard mean ocean water.

<figure><img src="image_13.png" /><figcaption>**Fig. 6** | Pearson's correlation between erosion rates and climatic parameters.</figcaption></figure>

Erosion rates compared to **a**, **b** 400-kyr eccentricity, **c**, **d** insolation at 23° S latitude, and **e**, **f** the CLIMBER-2 precipitation dataset for three different time periods: 1.8–2.48 Ma (red), 2.49–3.8 Ma (orange), and 3.8–6 Ma (or 5 Ma for the CLIMBER-2 dataset) (blue). Pearson's correlation coefficients and *p*-values are presented for

three periods of time below each plot with 95% confidence interval envelopes shown for each correlation. Note the strong correlation and low *p*-values for the relationships between erosion rate and 400-kyr eccentricity and CLIMBER-2 precipitation for the period from 2.49–3.8 Ma (orange) in **a** and **e**. Also note the different x-axis scalings between the upper and lower panels.

synchrony of observed erosion-rate cycles with the amplitude of precipitation changes predicted by the CLIMBER-2 model (Fig. 4b–d and Fig. 8a), we hypothesize that changes in erosion rates during this period are driven by changes in summer monsoon precipitation in the study area. Although we argue for precipitation changes as the primary driver of erosion rates, the preservation of orbitally-paced erosion-rate cyclicity during this limited time window (~4-2.4 Ma) was made possible by uniquely favorable tectonic and geomorphic conditions within the watershed. For example, the onset of the clear cyclicity at ~4.0 Ma corresponds to a major shift in sediment provenance described above, which is interpreted to record a more restricted sediment source, and the possible addition of sediments from the interior margin of the Puna Plateau³²(Supplementary Fig. 3).

Although the observed orbital pacing of erosion rates suggests a direct climatic control on long-term erosion rates, any model to explain this observed correlation must answer several key questions.

First, given that erosion rates are generally rapid enough to record changes at 23-kyr precessional timescales, why is this signal, or noise associated with undersampling it, not observed? Relatedly, if such a high-frequency signal is simply being 'smoothed' by the geomorphic system, why then should erosion rates rise and fall in 400-kyr cycles when modeled precipitation is essentially constant when averaged over multiple 23-kyr precessional cycles? Finally, how do we reconcile order-of-magnitude changes in observed erosion rates with apparently modest changes in modeled mean annual precipitation (~20%) during this period (Fig. 1b)?

Recognizing these questions, our preferred hypothesis to explain the observed 400-kyr erosion-rate pattern calls upon nonlinear threshold hillslope behavior driven by landsliding coupled to channel incision (e.g., lateral or vertical)⁵⁵,⁵⁶, causing erosion rates to nonlinearly track precipitation amplitudes above some threshold (Fig. 8)³⁷. During periods of increased aridity, sediment production and
<figure><img src="image_14.png" /><figcaption>**Fig. 7** | CLIMBER-2 wavelet and power spectra results. **a** CLIMBER-2 precipitation and **b** temperature wavelet analyses and power spectra. The grid cell of CLIMBER-2 covers the area from 30–20°S and 65–116.4°W. Precipitation shows strong 23-kyr periodicities from 4.5 to 3 Ma, with strong 41-kyr periodicities from 3 to 2 Ma. In contrast, the temperature spectra are largely dominated by 400-kyr periodicities throughout both periods with 100-kyr and 41-kyr periodicities increasing in power</figcaption></figure>

at <3 Ma. The continuous black line in the wavelet analysis indicates areas of the spectrum susceptible to distortion due to edge effects of the time series. The dashed red line in the power spectrum and the black outlined areas in the wavelet plots indicate the 95%-confidence interval. Horizontal dashed lines in the power spectrum indicate the dominant orbital frequencies at 23, 41, 100, and 400-kyr periodicities.

transport would diminish significantly, thereby allowing long-time-
scale, slower hillslope diffusion processes to dominate the erosion-rate
signal. During wetter periods, erosion rates would be dominated by
both an increased frequency of landsliding and greater exceedance of
transport and incision thresholds, volumetrically overwhelming and
diluting cosmogenic radionuclide signals generated during preceding
transport-limited periods: similar to model predictions of punctuated
denudation driven by increasing runoff intensity². Although steady-
state exploration of stochastic stream power models in the study area
during the Pliocene is beyond the scope of this study, previous work
indicates that climatic perturbations are most likely to influence ero-
sion rates in semiarid landscapes characterized by steep channels,
highly variable flood discharges, and sustained flood variability with
increasing mean annual runoff⁸. These characteristics as well as the
importance of landsliding frequency and landscape adjustment in
modulating sediment flux in the study region is well documented and
has been previously implicated in explaining order-of-magnitude,
erosion-rate variability observed for millennial-scale wet-arid climate
transitions during the late Quaternary¹³,⁵⁷,⁵⁹. In addition, the magnitude
of paleo-erosion rates calculated in this study agrees with empirical
landscape thresholds documented in the study region, whereby mean
catchment specific stream power values become asymptotic with
increasing catchment erosion rates between 0.1–0.5 mm/yr⁷, strongly
implicating landslide dynamics as a primary driver of erosion rates in
this landscape. The preferred hypothesis is also indirectly supported
by results from higher-resolution, regional paleoclimate modeling in
the study area that predicts greater interannual climatic variability
along the eastern Andes during the Pliocene compared to more recent
time periods (-21 ka, 6 ka, pre-industrial)⁶⁰, providing a potential cli-
matic mechanism for punctuated exceedance of erosional thresholds
and sustained runoff variability in the Río Iruya region.

Whereas our hypothesized mechanism provides a compelling
explanation that satisfies the observed results and is well-documented
in the modern study region, other potential explanations are certainly
possible. For example, the simplest explanation for the observed 400-
kyr pacing of erosion rates is that Pliocene moisture transport into the
study area was driven directly by 400-kyr eccentricity cycles instead of
modulating precession amplitudes. Although evidence exists for
strong pacing by long-eccentricity cycles in select marine sea-surface
temperature and carbon-isotope records during the Pliocene², such a
scenario contradicts the precession-dominated CLIMBER-2 precipita-
tion results during this time (Fig. 4 and Fig. 7). This explanation also
contradicts empirical evidence from late-Quaternary continental
records in the Andes documenting the dominance of higher-frequency

precipitation oscillations largely driven by precession and/or 100-kyr
eccentricity cycles<sup>53,54</sup>, though we acknowledge the potential dis-
connect between these geologically recent records and the climate
state of the Pliocene.

Alternatively, the 400-kyr cycles driving erosion rates may derive
from temperature oscillations instead of precipitation. The 400-kyr
frequency dominates the CLIMBER-2 temperature signal during the
Pliocene (Fig. 7b) and the effect of temperature oscillations would have
been especially pronounced if large glaciers and/or periglacial pro-
cesses were at work in the catchment during this period<sup>23,61</sup>. However,
there is sparse evidence to suggest that the Río Iruya catchment was
glaciated even during the Last Glacial Maximum (LGM) (-20 ka), with
most regional glaciers during the LGM small and localized to only the
highest elevations (> 4000 m)<sup>62</sup>. In addition, modern mean Austral
winter temperatures at ~2800 m elevation in the headwaters of the
study area currently average ~7–8 °C with average low temperatures of
~0–1 °C (town of Iruya). The Pliocene is known to have had a mean
global temperature similar, if not slightly warmer, to present day<sup>63</sup>,
making it unlikely that periglacial/glacial processes were responsible
for the observed erosion rate variability during the period of study.

Another possible mechanism for the observed erosion-rate
cyclicity is that catchment erosion rates stayed relatively constant
throughout the period from 4–2.4 Ma and that the observed signal
arises from variations in sediment-transport dynamics that impacted
¹⁰Be concentrations, namely sediment provenance concentrations
and/or stochastic burial during transport. One such scenario is that the
onset of cyclicity at 4 Ma was triggered by the Iruya headwaters tap-
ping into a previously restricted basin(s) on the eastern Puna Plateau
and rapidly evacuating its sedimentary fill into the foreland. Cyclic
changes in monsoonal strength could then extract pulses of low-
concentration sediments into the foreland (e.g., the Oran group)
during high-amplitude cycles, giving the impression of faster erosion.
Alternatively, cycles of enhanced fluvial storage at shallow burial
depths would allow additional ¹⁰Be to build up during the transport
process. Cyclic reworking and final burial of these sediments could
create periods of apparently slower erosion. These explanations
appear unlikely based on results from end-member modeling for both
scenarios (Supplementary Fig. 6). Periods of apparently rapid erosion
would require roughly 50–90% of sediment be derived from a source
with zero ¹⁰Be concentration to produce the observed cyclicity in an
otherwise steadily eroding landscape. Likewise, periods of apparently
slow erosion rates would require shallow sediment storage (<1 m
depth) over timescales of 10–30 kyr, which seems unlikely given the
proximity of the relatively steep Río Iruya watershed to the
<figure><img src="image_15.png" /><figcaption>**Fig. 8 |** **Paleo-erosion rates and conceptual mechanism for observed 400-kyr cycles.** a Calculated $^{10}$Be-derived paleo-erosion rates (±2σ envelope) from 1.8-3.8 Ma compared to modern $^{10}$Be-derived erosion rates (0.34 ± 0.01 mm/yr) and 400-kyr eccentricity pacing. Erosion-rate values shown are based on the preferred uplift scenario (Supplementary Fig. 4b). b Conceptual model illustrating how a 400-kyr erosion-rate signal with an order-of-magnitude variation in erosion rates might be produced from a nonlinear erosional response to precession-forced precipitation variability (in purple; Fig. 4f) above a hydrologic threshold for transport/landsliding. Black dots represent erosion rates that track the precipitation amplitude above the threshold (dashed black line) for each 23-kyr precession cycle. Note the nonlinear erosion-rate axis. SAMS = South American Monsoon System.</figcaption></figure>

depositional loci and the general lack of accommodation space in the modern watershed (e.g., narrow canyons and lack of paleo-terraces in all but the highest reaches) (Supplementary Fig. 7). Finally, cyclic tapping of long-buried sediments seems unlikely given that $^{26}$Al/$^{10}$Be ratios are mostly within uncertainty of their expected production rate ratios (Fig. 2C).

A fourth potential explanation is that the 400-kyr erosion-rate signal is the product of a climatically driven 23-kyr erosion-rate cycle that has been filtered by geomorphic processes within the watershed. For example, varying timescales of sediment transport and storage could act to overprint the underlying hillslope erosion rate. Much of our current understanding of how source watersheds respond to and convey climate signals to depositional areas relies on numerical models with conflicting arguments that include shredding, dampening, and even amplification of sediment-flux signals by autogenic geomorphic processes$^{64,65}$. Empirical results from nearby paleo-erosion rate studies$^{13,19}$ have proposed that catchment architecture and transport distance may dictate which periodicities are preserved in the sediment record, whereby preserved periodicities are proportional to the square of the river length above the site of deposition$^{19}$. Whereas this regional observation is compelling, the Río Iruya dataset

represents an integrated geochemical signature of catchment erosion that is largely independent from the autogenic geomorphic processes that drive fluvial terrace or alluvial fan preservation. Evidence from experimental studies, however, indicates that catchments undergoing pronounced aggradation/incision cycles, and when coupled to significant long-term sediment storage, may result in mixing temporally disparate erosion-rate signals during periods of incision as modern $^{10}$Be signals are progressively mixed with stratigraphically older, cosmogenically-decayed signals$^{66}$. Although some temporal mixing of cosmogenic radionuclide signals by geomorphic processes is necessary to allow precession amplitudes to closely track 400-kyr cycles based on our hypothesized mechanism, the storage times necessary to decay a shielded $^{10}$Be signal to produce the observed 400-kyr cyclicity in our record would be prohibitive (»10$^6$ years based on $^{10}$Be half-lives).

In summary, we posit that, during the Pliocene, precipitation was controlled by the strength of the SAMS, which oscillated on a 23-kyr timescale at an amplitude controlled by longer 400-kyr eccentricity cycles. The Río Iruya watershed, favorably conditioned by tectonically steepened and likely poorly vegetated threshold hillslopes (similar to the modern catchment), responded nonlinearly via landslide dynamics to this precipitation regime, allowing relatively small perturbations in precipitation (±20%) to produce order-of-magnitude changes in erosion rates in concert with 400-kyr pacing (Fig. 8).

Another interesting feature of the paleo-erosion rate dataset is the overall increase in erosion rate and loss of the 400-kyr cyclicity starting at ~2.4 Ma (Fig. 5a). One possible explanation is that the loss of 400 kyr erosion-rate cyclicity may have been driven by a watershed reorganization in the headwaters of the Río Iruya associated with the ~2.4 Ma provenance shift. Alternatively, the loss of erosion-rate cyclicity may have resulted from a change in precipitation conditions associated with the continued growth of the Eastern Cordillera, eventually concentrating precipitation along the eastern range front as the orographic barrier became more pronounced. This explanation is supported during this time by an order-of-magnitude decrease in erosion rates on the leeward side of the growing Eastern Cordillera in the Quebrada de Humahuaca$^{21}$ (Fig. 5 and Supplementary Fig. 3). The ongoing tectonism and increased precipitation on the windward flank of the Eastern Cordillera, when combined with a potentially shorter, steeper catchment would have conspired to drive the more rapid erosion rates observed following the Plio-Pleistocene transition. However, it remains unclear why the newly reorganized catchment would not continue to record precipitation cyclicity during this period, given the relatively faithful tracking of 400-kyr cyclicity over the preceding ~1.5 Myr period (Fig. 5a).

A third possible explanation for the loss of erosion-rate cyclicity at ~2.4 Ma is a change in the global climate state driven by the onset of major Northern Hemisphere glaciation at the Plio-Pleistocene boundary around 2.6 Ma. Growth of the Northern Hemisphere ice sheets ushered in a climate regime largely paced by obliquity (~41 kyr) as the earth progressively cooled from the Pliocene into the Pleistocene:$^{51,67}$ a finding supported by the CLIMBER-2 precipitation results as well (Fig. 4f and Fig. 7). The CLIMBER-2 data also suggest that, as obliquity begins to dominate the precipitation signal between 3–2.5 Ma, the amplitude becomes modulated by eccentricity mega-cycles (~1.2 Myr), producing both higher amplitude oscillations than the precession cycles and decreased variability between individual oscillations due to the protracted length of the modulating cycle (Fig. 4f). If erosion rates continued to track amplitude modulations of the dominant obliquity signal instead of precession during this time, it is plausible that the more gradual changes in amplitude associated with the eccentricity mega-cycles failed, or were too subtle, to be recorded by erosion rates within the catchment. Alternatively, the strengthening of the obliquity signal in combination with a diminishing precession signal during the early Pleistocene may have produced a more chaotic temporal precipitation regime that, when coupled to a newly evolving basin
architecture, prevented any dominant climate forcing frequency to emerge in the erosion-rate data. More broadly, the SAMS may also have undergone significant changes during this period as El Niño Southern Oscillation conditions supplanted those of the permanent El Niño state⁶⁸,⁶⁹ and the Intertropical Convergence Zone was pushed southward due to increased interhemispheric temperature gradients driven by Northern Hemisphere ice-sheet growth⁷⁰.

## Methods

### Sample age model and uncertainty

Detailed timing in the Río Iruya section was established using a densely sampled magnetostratigraphic transect combined with several precisely dated, intercalated volcanic ash beds yielding robust age control in the section that spans 1.94 to 6.49 Ma³² (Supplementary Fig. 2 and Supplementary Dataset 1). Calculation of age uncertainties on paleo-erosion rate samples was performed using published magnetostratigraphic data and allowed the polarity reversal boundaries to be flexible with respect to their confining samples or radiometric ash ages. Assuming a linear accumulation rate between polarity reversals, the maximum and minimum age of each cosmogenic radionuclide sample was determined by moving the reversal boundaries as far forward (making the sample older) and as far back (making the sample younger) as possible in the stratigraphy (Supplementary Fig. 8). This method allows for proper quantification of sample-age uncertainty and highlights the asymmetric nature of age uncertainty associated with unevenly spaced sampling and polarity reversal boundaries. Midpoint, maximum, and minimum age uncertainties for each sample are reported in Supplementary Dataset 2, with a median maximum and minimum divergence from the midpoint age of 0.033 Myr (n = 41) and 0.034 Myr (n = 39), respectively. For samples in the lower part of the section defined only by ashes (5.5–6.49 Ma), the midpoint age was calculated using a linear interpolation based on stratigraphic location between the ashes. Uncertainties for these samples were cautiously set as the 75th percentile of maximum and minimum uncertainties from the magnetostratigraphy component of the dataset (0.048 and 0.049 Ma, respectively). In the case of the youngest sample (IR22c), which lacked lower bound age control, the midpoint age was calculated using the average sedimentation rate for the entire section (930 m/Ma), yielding an age of 1.83 Ma with uncertainties again assumed to be within the overall 75th percentile. In cases where a maximum age uncertainty could be defined but not a minimum one (IR20.1c and IR21c), due to the topmost reversal boundary lying down section from the cosmogenic radionuclide sample, we assumed uncertainty symmetry about the midpoint and used the maximum uncertainty value for both uncertainties.

### Cosmogenic radionuclide sampling strategy and analytical procedures

Quartz-rich, cobble conglomerate beds -1-m thick were excavated in the Oran Group section along the Río Iruya and sieved in the field to obtain the <1 mm interstitial size fraction. Samples were then further sieved at Middlebury College to obtain the 0.25- to 1-mm fraction. This fraction then followed a standard quartz cleaning protocol⁷¹ whereby samples were subjected to carbonate dissolution using a heated 1:1 hydrochloric acid mixture and repeated 2% hydrofluoric acid leachings performed using commercial hot dog rollers and ultrasonic baths to provide heat and mechanical abrasion. Sample purities were checked following three rounds of 2% hydrofluoric acid leaching and continued until quartz purity targets were met. Quartz dissolution (ranging from 95 to 621 g per sample) using hydrofluoric acid, ¹⁰Be and ²⁶Al extraction by ion-exchange chromatography, and sample-target preparation were all performed at the University of California at Santa Barbara using established and publicly available protocols⁵⁷,⁷¹. During digestion of the quartz samples, a low-ratio ⁹Be spike with a ¹⁰Be/⁹Be ratio value of 3 × 10⁻¹⁵ was added. In addition, during each digestion batch an

empty teflon beaker was treated with approximately the same maximum volume of hydrofluoric acid as used in that given batch of samples and was spiked with the same low-ratio ⁹Be spike. These blanks were then processed through the same protocols as samples, yielding ¹⁰Be/⁹Be values for the 5 blanks ranging from 3.5 to 13 × 10⁻¹⁵ (Supplementary Dataset 2 and 3). No ²⁷Al spike was added, given the high concentrations of Al in the quartz mineral grain lattice, with no ²⁶Al contamination detected in the processed blanks during the chemical processing steps. ¹⁰Be and ²⁶Al targets were then measured via accelerator mass spectrometry at the Purdue Rare Isotope Measurement Laboratory (PRIME) using the 07KNSTD standardization and assuming a ¹⁰Be half-life of 1.36 ± 0.07 Myr⁷². Reference standard KNSTD was used for ²⁶Al assuming a half-life of 0.705 Myr³⁰. In total, 49 ¹⁰Be and 23 ²⁶Al samples were processed and analyzed with ²⁶Al measurements restricted to samples with ages <3.5 Ma (-5 ²⁶Al half-lives).

### Modern erosion rates from ¹⁰Be

Modern erosion rates were calculated using ¹⁰Be for catchment samples taken -3 km upstream of the sampled section of the Río Iruya canyon (IR-3M) as well as further upstream in the headwaters of the Río Iruya (IR-6M) (Supplementary Fig. 3). The production of ¹⁰Be in quartz in the upper few meters of the Earth's surface arises predominantly from high-energy neutron bombardment of O and Si nuclei, with negative muon capture and fast muon reactions dominant at depths greater than -3 meters³¹,⁷³,⁷⁴. Production P of ¹⁰Be can be modeled as an exponential decay function with depth from the surface for each production pathway where

$$P = P_0 e^{-zP/\Lambda} \quad (1)$$

and P₀ (at g⁻¹ yr⁻¹) is the surface production rate, z (cm) is the depth below the surface, ρ the density of the regolith (g cm⁻³), and Λ the attenuation mean free path (g cm⁻²)³¹. Because P₀ varies with latitude, altitude, and topographic shielding, modern mean effective production rates (topographic shielding factor x production rate) across the Río Iruya watershed were estimated using the 90-meter Shuttle Radar Topography Mission digital elevation model and the CRONUScalc Matlab codes⁷⁵ assuming a sea-level high-latitude (SLHL) production rate of 4.01 atoms g⁻¹ yr⁻¹ and based on the time-independent Lal-Stone scaling (St)³¹,⁷⁶,⁷⁷. Assuming rapid transport of sediments from the hillslopes into the fluvial network and proper amalgamation in the fluvial system of an integrated catchment signal, the concentration of ¹⁰Be in a river sand sample can be used to calculate the erosion rate, R_E, of the upstream catchment²⁴ where

$$R_E = \frac{\Lambda}{\rho} \left( \frac{P_0}{N_E} - \lambda \right) \quad (2)$$

and Λ is the radioactive decay constant for ¹⁰Be (yr⁻¹) and N_E is the ¹⁰Be concentration of a modern river sample (atoms g⁻¹).

### Paleo-erosion rate calculations and uncertainty propagation

Paleo-erosion rate calculations were performed using a simplified version of the detailed methodology presented in¹⁷, which tracks and compounds uncertainties across the different cosmogenic radionuclide ingrowth components. Here we omit any treatment of the excavated depth of sampling or paleochannel thickness, due to the rapid exhumation rates (-1 m/yr for 150 years) in the Río Iruya and the lack of paleochannel-derived samples that might warrant such nuanced treatment. The major ingrowth and decay components of the measured cosmogenic radionuclide signal (N_A) are modeled by the linear combination:

$$N_A = N_E + N_T + N_B - N_D + N_X \quad (3)$$
where N<sub>A</sub> is the total measured concentration of each sample, which encompasses the contribution from erosion in the watershed (N<sub>E</sub> – the value of interest), transport to the site of deposition (N<sub>T</sub>), spallation-induced and muogenic production during burial (N<sub>B</sub>), radioactive decay during burial (N<sub>D</sub>), and lastly, exhumation to the surface where the sediment is re-exposed to cosmogenic bombardment (N<sub>X</sub>) and then sampled<sup>17</sup>. Production from spallation and fast and slow muons is integrated throughout both the burial and exhumation terms (cf. Supplementary Table 1), though production during exhumation in this study was negligible due to the rapid incision of the canyon. Each component and associated uncertainty were constrained and inverted progressively back through time (N<sub>X</sub> to N<sub>E</sub>) utilizing a Monte Carlo approach (100,000 iterations) to compound uncertainty in the erosion-rate component (N<sub>E</sub>) and calculated erosion rate (R<sub>E</sub>) for each sample based on Eq. 2 (Supplementary Dataset 2 and 3 and Supplementary Table 1). Resulting asymmetric uncertainty distributions derive from the exponential depth dependence of cosmogenic radionuclide production, the time dependence of decay, and asymmetry in age uncertainties (Supplementary Fig. 8). Whereas uncertainties were not assigned to production rate values in the source area (P<sub>E</sub>) and foreland depocenter (P<sub>B</sub>), several scenarios of production-rate evolution through time were modeled to account for a wide range of possible catchment uplift scenarios (see Supplementary Fig. 4 and Supplementary Dataset 2). These differing scenarios, it should be noted, only result in minor changes to the absolute magnitude of calculated paleo-erosion rates from 1.8 Ma to 4 Ma and not the observed patterns or cyclicity. The transport contribution (N<sub>T</sub>) to measured cosmogenic radionuclide concentrations in the analysis is assumed to be negligible based on: (1) the close proximity of the depositional basin to a relatively steep source area starting at ~4 Ma;<sup>32,44</sup> and (2) good agreement between <sup>26</sup>Al:<sup>10</sup>Be sample ratios (for <sup>26</sup>Al AMS errors <35%) with expected <sup>26</sup>Al:<sup>10</sup>Be ratios assuming no significant burial during transport from ~2–3 Ma (Fig. 2c). Our calculations assume a relative contribution by nucleon spallation, negative muon capture, and fast muon reactions at the surface equal to 98.15%, 1.2%, and 0.65%, respectively<sup>73,74,78</sup>. To calculate the decay of <sup>10</sup>Be during burial, a half-life of 1.387 Myr<sup>29</sup> was used as compared to the half-life of 1.36 Myr<sup>72</sup> assumed during analytical measurement. All additional calculated, assumed, and measured parameters are detailed in Supplementary Table 1. Note that no attempt was made to derive paleo-erosion rates from <sup>26</sup>Al due to the larger analytical uncertainties resulting from the dilution of the <sup>26</sup>Al/<sup>27</sup>Al ratio by abundant <sup>27</sup>Al quartz impurities in the sample material. A detailed presentation and discussion of the uncertainty workflow and equations used in the paleo-erosion rate analysis, as well as open-source R code, is available in<sup>17</sup>.

CLIMBER-2 climate model simulations and data

The Earth system model of intermediate complexity, (EMIC) CLIMBER-2, was chosen to test the hypothesis that climate perturbations in the Rio Iruya region largely result from global-scale forcings during the Plio-Pleistocene. The low-resolution, intermediate complexity model seeks to incorporate a large set of processes and feedbacks normally reserved for comprehensive climate models, but at a low enough spatial resolution to avoid their computational costs, thereby yielding fast turnaround times for long time-scale paleoclimate reconstructions. Whereas the low-spatial resolution of the model is incapable of resolving orographic and vegetative forcings associated with the Andes (and thereby associated regional systems like the South American low level jet), CLIMBER-2 has been shown to reliably replicate both large-scale climate reconstructions from the Last Glacial Maximum, as well as more complex model simulations for the Pliocene<sup>79,80</sup>.

Here we implemented version 2.3 of CLIMBER-2 that includes atmospheric, ocean, and terrestrial components, including a dynamic vegetation model and thermodynamic sea-ice model<sup>50,79</sup>. The CLIMBER-2 atmospheric component is represented by a 2.5-dimensional statistical-dynamical model with latitudinal resolution of 10° and

longitudinal resolution of ~51° (7 zones) using a daily time step. Circulation, temperature, and humidity are calculated across 10 vertical levels in the model and 16 levels for long-wave radiation. Vegetation potential in each grid cell is recorded at a yearly time-step based on temperature and precipitation using the terrestrial model VECODE<sup>81</sup>. The ocean model integrates the Pacific, Atlantic, and Indian oceans via connection through the Southern Ocean and calculates zonally averaged flow between grid cells, as well as temperature and salinity changes, at a time step of 5 days, a latitudinal resolution of 2.5°, and across 20 unevenly distributed vertical levels<sup>82</sup>.

Model simulations were performed using a transient 5-million-year simulation which utilizes an established solution for astronomical parameters<sup>83</sup> to calculate insolation, along with prescribed changes of ice-sheet topography and atmospheric CO<sub>2</sub> concentrations<sup>79</sup>. Ice-sheet topography change is adopted from previous simulations of 3-D ice-sheet models<sup>52</sup> driven by the LR04 benthic δ18O stack record<sup>67</sup>. Temperature derived from the simulation with the 3-D ice-sheet model is used to reconstruct CO<sub>2</sub> concentration using previously described methods<sup>84</sup>. Whereas precipitation and temperature are recorded daily for each cell in the model, outputs in this paper were derived from 1000-year time-step averages to filter out sub-millennial model variability not pertinent to orbital time scales. Frequency analysis was performed using AnalySeries software version 2.0.8. Both precession and long-eccentricity frequency filtering components were acquired using Gaussian shaped filters at 0.044 and 0.0025 with a bandwidth of 0.008 and 0.001, respectively. Wavelet analyses were performed using previously described methods<sup>85</sup>. Due to the simple spatial representation of landmass in the CLIMBER-2 model, a geographic offset was accounted for in assessing data for South America (e.g., Fig. 1b) as the continent is shifted westward by ~35–45 degrees of longitude in the model (see Supplementary Fig. 5). CLIMBER-2 outputs used in this paper are provided in Supplementary Dataset 4.

Data availability

All cosmogenic radionuclide and CLIMBER-2 data generated in this study are provided in the Supplementary Information.

Code availability

All computational analyses rely on previously published models and code directly cited within the methods of the manuscript.

References

1. Gilbert, G. K. Rhythms and geologic time. *Science* **11**, 1001–1012 (1900).

2. Tucker, G. E. & Slingerland, R. Drainage basin responses to climate change. *Water Resour. Res.* **33**, 2031–2047 (1997).

3. Willett, S. D. Orogeny and orography: The effects of erosion on the structure of mountain belts. *J. Geophys. Res.-Solid Earth* **104**, 28957–28981 (1999).

4. Godard, V., Tucker, G. E., Burch Fisher, G., Burbank, D. W. & Bookhagen, B. Frequency-dependent landscape response to climatic forcing. *Geophys. Res. Lett.* **40**, 859–863 (2013).

5. Braun, J., Voisin, C., Gourlan, A. T. & Chauvel, C. Erosional response of an actively uplifting mountain belt to cyclic rainfall variations. *Earth Surf. Dyn.* **3**, 1–14 (2015).

6. Whipple, K. X. The influence of climate on the tectonic evolution of mountain belts. *Nat. Geosci.* **2**, 97–104 (2009).

7. Bookhagen, B. & Strecker, M. R. Orographic barriers, high-resolution TRMM rainfall, and relief variations along the eastern Andes. *Geophys. Res. Lett.* **35**, L06403, https://doi.org/10.1029/2007GL032011 (2008).

8. Willenbring, J. K. & von Blanckenburg, F. Long-term stability of global erosion rates and weathering during late-Cenozoic cooling. *Nature* **465**, 211–214 (2010).
9. Molnar, P. Late Cenozoic increase in accumulation rates of terrestrial sediment: how might climate change have affected erosion rates? *Annu. Rev. Earth Planet. Sci.* **32**, 67–89 (2004).

10. Schaller, M. et al. Paleoerosion rates from cosmogenic 10Be in a 1.3 Ma terrace sequence: Response of the River Meuse to changes in climate and rock uplift. *J. Geol.* **112**, 127–144 (2004).

11. Charreau, J. et al. Paleo-erosion rates in Central Asia since 9Ma: A transient increase at the onset of Quaternary glaciations? *Earth Planet. Sci. Lett.* **304**, 85–92 (2011).

12. Hidy, A. J., Gosse, J. C., Blum, M. D. & Gibling, M. R. Glacial-interglacial variation in denudation rates from interior Texas, USA, established with cosmogenic nuclides. *Earth Planet. Sci. Lett.* **390**, 209–221 (2014).

13. Schildgen, T. F. et al. Landscape response to late Pleistocene climate change in NW Argentina: Sediment flux modulated by basin geometry and connectivity. *J. Geophys. Res.: Earth Surf.* **121**, 392–414 (2016).

14. Val, P., Hoke, G. D., Fosdick, J. C. & Wittmann, H. Reconciling tectonic shortening, sedimentation and spatial patterns of erosion from 10Be paleo-erosion rates in the Argentine Precordillera. *Earth Planet. Sci. Lett.* **450**, 173–185 (2016).

15. Marshall, J. A., Roering, J. J., Gavin, D. G. & Granger, D. E. Late Quaternary climatic controls on erosion rates and geomorphic processes in western Oregon, USA. *Geol. Soc. Am. Bull.* **129**, 715–731 (2017).

16. Puchol, N. et al. Limited impact of Quaternary glaciations on denudation rates in Central Asia. *Geol. Soc. Am. Bull.* **129**, 479–499 (2017).

17. Oskin, M. E. et al. Steady 10 Be-derived paleoerosion rates across the Plio-Pleistocene climate transition, Fish Creek-Vallecito basin, California. *J. Geophys. Res.: Earth Surf.* **122**, 1653–1677 (2017).

18. Amidon, W. H. et al. Mio-Pliocene aridity in the south-central Andes associated with Southern Hemisphere cold periods. *Proc. Natl Acad. Sci.* **114**, 6474–6479 (2017).

19. Tofelde, S. et al. 100 kyr fluvial cut-and-fill terrace cycles since the Middle Pleistocene in the southern Central Andes, NW Argentina. *Earth Planet. Sci. Lett.* **473**, 141–153 (2017).

20. Madella, A., Delunel, R., Akçar, N., Schlunegger, F. & Christl, M. ¹⁰Be- inferred paleo-denudation rates imply that the mid-Miocene western central Andes eroded as slowly as today. *Sci. Rep.* **8**, 2299 (2018).

21. Pingel, H., Schildgen, T., Strecker, M. R. & Wittmann, H. Pliocene-Pleistocene orographic control on denudation in north-west Argentina. *Geology* **47**, 359–362 (2019).

22. Lenard, S. J. P. et al. Steady erosion rates in the Himalayas through late Cenozoic climatic changes. *Nat. Geosci.* **13**, 448–452 (2020).

23. Mariotti, A. et al. Nonlinear forcing of climate on mountain denudation during glaciations. *Nat. Geosci.* **14**, 16–22 (2021).

24. Charreau, J. et al. A 6 Ma record of palaeodenudation in the central Himalayas from in situ cosmogenic 10 Be in the Surai section. *Basin Res.* **33**, 1218–1239 (2021).

25. Mandal, S. K., Scherler, D. & Wittmann, H. Tectonic accretion controls erosional cyclicity in the Himalaya. *AGU Adv.* **2**, e2021AV000487 (2021).

26. Rohrmann, A. et al. Miocene orographic uplift forces rapid hydrological change in the southern central Andes. *Sci. Rep.* **6**, 35678, https://doi.org/10.1038/srep35678 (2016).

27. Hoorn, C. et al. Amazonia through time: andean uplift, climate change, landscape evolution, and biodiversity. *Science* **330**, 927–931 (2010).

28. von Blanckenburg, F. The control mechanisms of erosion and weathering at basin scale from cosmogenic nuclides in river sediment. *Earth Planet. Sci. Lett.* **242**, 224–239 (2006).

29. Chmeleff, J., von Blanckenburg, F., Kossert, K. & Jakob, D. Determination of the 10Be half-life by multicollector ICP-MS and liquid scintillation counting. *Nucl. Inst. Methods Phys. Res., B* **268**, 192–199 (2010).

30. Nishiizumi, K. Preparation of 26Al AMS standards. *Nucl. Instrum. Methods Phys. Res. Sect. B: Beam Interact. Mater. At.* **223-224**, 388–392 (2004).

31. Lal, D. & Labeling, Cosmic-Ray of Erosion Surfaces—Insitu nuclide production-rates and erosion models. *Earth Planet. Sci. Lett.* **104**, 424–439 (1991).

32. Amidon, W. H. et al. Provenance and tectonic implications of Orán Group foreland basin sediments, Rio Iruya canyon, NW Argentina (23° S). *Basin Res.* **29**, 96–112 (2015).

33. Echavarria, L., Hernández, R., Allmendinger, R. & Reynolds, J. Subandean thrust and fold belt of northwestern Argentina: Geometry and timing of the Andean evolution. *AAPG Bull.* **87**, 965–985 (2003).

34. Hernández, R. M., Reynolds, J. & Di, A. Salvo, Análisis tectosedimentario y ubicación geocronológica del Grupo Orán en el río Iruya. *Boletín de. Informaciones Petroleras* **45**, 80–93 (1996).

35. Starck, D. & Schulz, A. La configuración estructural del límite Cordillera Oriental – Sierras Subandinas en el extremo norte de la República Argentina. *XII Congr. Geol.ógico Boliviano, Mem.* **32**, 221–229 (1996).

36. Hongn, F. et al. Middle Eocene deformation and sedimentation in the Puna-Eastern Cordillera transition (23°–26°S): Control by pre-existing heterogeneities on the pattern of initial Andean shortening. *Geology* **35**, 271 (2007).

37. Insel, N. et al. Paleozoic to early Cenozoic cooling and exhumation of the basement underlying the eastern Puna plateau margin prior to plateau growth. *Tectonics* **31**, 1–17 (2012).

38. Kleinert, K. & Strecker, M. R. Climate change in response to orographic barrier uplift: Paleosol and stable isotope evidence from the late Neogene Santa Maria basin, northwestern Argentina. *Geol. Soc. Am. Bull.* **113**, 728–742 (2001).

39. Mulch, A., Uba, C. E., Strecker, M. R., Schoenberg, R. & Chamberlain, C. P. Late Miocene climate variability and surface elevation in the central Andes. *Earth Planet. Sci. Lett.* **290**, 173–182 (2010).

40. DeCelles, P. G., Carrapa, B., Horton, B. K. & Gehrels, G. E. Cenozoic foreland basin system in the central Andes of northwestern Argentina: Implications for Andean geodynamics and modes of deformation. *Tectonics* **30**, 1–30 (2011).

41. Reynolds, J. H. et al. Middle Miocene tectonic development of the transition zone, Salta Province, northwest Argentina: Magnetic stratigraphy from the Metan Subgroup, Sierra de Gonzalez. *Geol. Soc. Am. Bull.* **112**, 1736–1751 (2000).

42. Reynolds, J. H., Hernández, R. M., Galli, C. I. & Idleman, B. D. Magnetostatigraphy of the Quebrada La Porcelana section, Sierra de Ramos, Salta Province, Argentina: age limits for the neogene Orán Group and uplift of the southern Sierras Subandinas. *J. South Am. Earth Planet.* **14**, 681–692 (2001).

43. Rahl, J. M., Harbor, D. J., Galli, C. I. & O'Sullivan, P. Foreland basin record of uplift and exhumation of the Eastern Cordillera, Northwest Argentina. *Tectonics* **37**, 4173–4193 (2018).

44. Pingel, H., Strecker, M. R., Alonso, R. N. & Schmitt, A. K. Neotectonic basin and landscape evolution in the Eastern Cordillera of NW Argentina, Humahuaca Basin (-24°S). *Basin Res.* **25**, 554–573 (2013).

45. Streit, R. L. et al. Controls on intermontane basin filling, isolation and incision on the margin of the Puna Plateau, NW Argentina (-23°S). *Basin Res.* **29**, 131–155 (2015).

46. Strecker, M. R. et al. Tectonics and climate of the Southern Central Andes. *Annu. Rev. Earth Planet. Sci.* **35**, 747–787 (2007).

47. Vuille, M. et al. A review of the South American monsoon history as recorded in stable isotopic proxies over the past two millennia. *Climate* **8**, 1309–1321 (2012).

48. Garreaud, R. D., Molina, A. & Farias, M. Andean uplift, ocean cooling and Atacama hyperaridity: A climate modeling perspective. *Earth Planet. Sci. Lett.* **292**, 39–50 (2010).
49. Uba, C. E., Strecker, M. R. & Schmitt, A. K. Increased sediment accumulation rates and climatic forcing in the central Andes during the late Miocene. *Geology* **35**, 979 (2007).

50. Petoukhov, V. et al. CLIMBER-2: a climate system model of intermediate complexity. Part I: model description and performance for present climate. *Clim. Dyn.* **16**, 1–17 (2000).

51. Zachos, J., Pagani, M., Sloan, L., Thomas, E. & Billups, K. Trends, rhythms, and aberrations in global climate 65 Ma to present. *Science* **292**, 686 (2001).

52. de Boer, B., Lourens, L. J. & van de Wal, R. S. W. Persistent 400,000-year variability of Antarctic ice volume and the carbon cycle is revealed throughout the Plio-Pleistocene. *Nat. Commun.* **5**, 2999 (2014).

53. Cheng, H. et al. Climate change patterns in Amazonia and biodiversity. *Nat. Commun.* **4**, 1411 (2013).

54. Baker, P. A. & Fritz, S. C. Nature and causes of Quaternary climate variation of tropical South America. *Quat. Sci. Rev.* **124**, 31–47 (2015).

55. Burbank, D. W. et al. Bedrock incision, rock uplift and threshold hillslopes in the northwestern Himalayas. *Nature* **379**, 505–510 (1996).

56. Larsen, I. J. & Montgomery, D. R. Landslide erosion coupled to tectonics and river incision. *Nat. Geosci.* **5**, 468–473 (2012).

57. Bookhagen, B. & Strecker, M. R. Spatiotemporal trends in erosion rates across a pronounced rainfall gradient: Examples from the southern Central Andes. *Earth Planet. Sci. Lett.* **327–328**, 97–110 (2012).

58. DiBiase, R. A. & Whipple, K. X. The influence of erosion thresholds and runoff variability on the relationships among topography, climate, and erosion rate. *J. Geophys. Res.* **116**, 1–17 (2011).

59. Trauth, M. H., Bookhagen, B., MARWAN, N. & Strecker, M. R. Multiple landslide clusters record Quaternary climate changes in the northwestern Argentine Andes. *Palaeogeogr., Palaeoclimatol., Palaeoecol.* **194**, 109–121 (2003).

60. Mutz, S. G. et al. Estimates of late Cenozoic climate change relevant to Earth surface processes in tectonically active orogens. *Earth Surf. Dyn.* **6**, 271–301 (2018).

61. Marshall, J. A. et al. Frost for the trees: Did climate increase erosion in unglaciated landscapes during the late Pleistocene? *Sci. Adv.* **1**, e1500715–e1500715 (2015).

62. D'Arcy, M. et al. Timing of past glaciation at the Sierra de Aconquija, northwestern Argentina, and throughout the Central Andes. *Quat. Sci. Rev.* **204**, 37–57 (2019).

63. Fedorov, A. V. et al. Patterns and mechanisms of early Pliocene warmth. *Nature* **496**, 43–49 (2013).

64. Straub, K. M., Duller, R. A., Foreman, B. Z. & Hajek, E. A. Buffered, incomplete, and shredded: the challenges of reading an imperfect stratigraphic record. *J. Geophys. Res.: Earth Surf.* **125**, 1–44 (2020).

65. Romans, B. W., Castellort, S., Covault, J. A., Fildani, A. & Walsh, J. P. Environmental signal propagation in sedimentary systems across timescales. *Earth Sci. Rev.* **153**, 7–29 (2016).

66. Tofelde, S., Savi, S., Wickert, A. D., Bufe, A. & Schildgen, T. F. Alluvial channel response to environmental perturbations: fill-terrace formation and sediment-signal disruption. *Earth Surf. Dyn.* **7**, 609–631 (2019).

67. Lisiecki, L. E., Raymo, M. E. & Pliocene-Pleistocene, A. stack of 57 globally distributed benthic δ18O records. *Paleoceanography* **20**, 1–17 (2005).

68. Lease, R. O. & Ehlers, T. A. Incision into the Eastern Andean Plateau During Pliocene Cooling. *Science* **341**, 774–776 (2013).

69. Fedorov, A. V. et al. The Pliocene paradox (mechanisms for a permanent El Nino). *Science* **312**, 1485–1489 (2006).

70. Schneider, T., Bischoff, T. & Haug, G. H. Migrations and dynamics of the intertropical convergence zone. *Nature* **513**, 45–53 (2014).

71. Kohl, C. P. & Nishiizumi, K. Chemical isolation of quartz for measurement of insitu-produced cosmogenic nuclides. *Geochimica et Cosmochimica Acta* **56**, 3583–3587 (1992).

72. Nishiizumi, K. et al. Absolute calibration of 10Be AMS standards. *Nucl. Instrum. Methods Phys. Res. Sect. B: Beam Interact. Mater.* **At. 258**, 403–413 (2007).

73. Heisinger, B. et al. Production of selected cosmogenic radionuclides by muons 1. Fast muons. *Earth Planet. Sci. Lett.* **200**, 345–355 (2002).

74. Heisinger, B. et al. Production of selected cosmogenic radionuclides by muons: 2. Capture of negative muons. *Earth Planet. Sci. Lett.* **200**, 357–369 (2002).

75. Marrero, S. M. et al. Cosmogenic nuclide systematics and the CRONUScalc program. *Quat. Geochronol.* **31**, 160–187 (2016).

76. Stone, J. O. Air pressure and cosmogenic isotope production. *J. Geophys. Res.-Solid Earth* **105**, 23753–23759 (2000).

77. Borchers, B. et al. Geological calibration of spallation production rates in the CRONUS-Earth project. *Quat. Geochronol.* **31**, 188–198 (2016).

78. Braucher, R., Brown, E. T., Bourlès, D. L. & Colin, F. In situ produced 10Be measurements at great depths: implications for production rates by fast muons. *Earth Planet. Sci. Lett.* **211**, 251–258 (2003).

79. Stap, L. B. et al. Modeled influence of land ice and CO₂ on polar amplification and paleoclimate sensitivity during the past 5 million years. *Paleoceanogr. Paleoclimatology* **33**, 381–394 (2018).

80. Willeit, M., Ganopolski, A., Calov, R. & Brovkin, V. Mid-Pleistocene transition in glacial cycles explained by declining CO₂ and regolith removal. *Sci. Adv.* **5**, eaav7337 (2019).

81. Brovkin, V., Ganopolski, A. & Svirezhev, Y. A continuous climate-vegetation classification for use in climate-biosphere studies. *Ecol. Model.* **101**, 251–261 (1997).

82. Stocker, T. F., Mysak, L. A. & Wright, D. G. A zonally averaged, coupled ocean-atmosphere model for paleoclimate studies. *J. Clim.* **5**, 773–797 (1992).

83. Laskar, J. et al. A long-term numerical solution for the insolation quantities of the Earth. *Astron. Astrophys.* **428**, 261–285 (2004).

84. Van de Wal, R. S. W., de Boer, B., Lourens, L. J., Köhler, P. & Bintanja, R. Reconstruction of a continuous high-resolution CO₂ record over the past 20 million years. *Climate* **7**, 1459–1469 (2011).

85. Grinsted, A., Moore, J. C. & Jevrejeva, S. Application of the cross wavelet transform and wavelet coherence to geophysical time series. *Nonlinear Process. Geophysics* **11**, 561–566 (2004).

86. Pingel, H. et al. Pliocene orographic barrier uplift in the southern Central Andes. *Geology* **42**, 691–694 (2014).

## Acknowledgements

This work was supported by NSF grant EAR-1148233 to W.H.A and D.W.B and a Jackson Distinguished Postdoctoral Fellowship to G.B.F at the University of Texas at Austin. B.d.B was funded by NWO Earth and Life Sciences (ALW)(#863.15.019) while the CLIMBER-2 dataset has been generated by E.T. within the framework of the NWO-ALW funded Vici project (#865.10.001) of L.J.L. We thank Louise McCarren for assistance in the field and Taylor Schildgen and Heiko Pingel for invaluable discussions regarding the manuscript.

## Author contributions

D.W.B and W.H.A conceived the project. G.B.F, W.H.A, D.W.B, L.V.L., and R.N.A. conducted all fieldwork and topographic data analyses. G.B.F, W.H.A, L.V.L., and B.B. performed all chemical separation and target preparation of cosmogenic radionuclide samples. G.B.F, W.H.A, L.V.L., B.B., and M.E.O analyzed the cosmogenic radionuclide data and performed error analysis. G.B.F, W.H.A, L.V.L., D.W.B., B.B., and V.G. developed conceptual model for observed erosion-rate cyclicity. B.d.B, L.B.S., E.T., and L.J.L. performed, processed, and analyzed the CLIMBER-2 simulations. G.B.F. wrote the manuscript with insights and contributions from all authors.

## Competing interests

The authors declare no competing interests.
**Additional information**

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s41467-023-36022-0.

**Correspondence** and requests for materials should be addressed to G. Burch Fisher.

**Peer review information** *Nature Communications* thanks Martin Clausen, Isaac Larsen, Nadine McQuarrie and the other, anonymous, reviewer(s) for their contribution to the peer review of this work.

**Reprints and permissions information** is available at http://www.nature.com/reprints

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2023