# Thickness of Fluvial Deposits Records Climate Oscillations

X. Yuan, L. Guerit, J. Braun, D. Rouby, C. Shobe

► To cite this version:

X. Yuan, L. Guerit, J. Braun, D. Rouby, C. Shobe. Thickness of Fluvial Deposits Records Climate Oscillations. Journal of Geophysical Research: Solid Earth, 2022, 127 (4), pp.e2021JB023510. 10.1029/2021JB023510. insu-03747061

HAL Id: insu-03747061

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

Submitted on 8 Aug 2022

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>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

JGR Solid Earth

RESEARCH ARTICLE

10.1029/2021JB023510

**Key Points:**

* We use a numerical model and a new analytical solution to quantify a physical link between fluvial deposits and climate oscillations

* Our method provides a theoretical framework for extracting information on past climate variations from fluvial terrace deposits

* Our results explain time lag of 20%–25% of forcing period commonly observed between the timing of maximum rainfall and erosion

**Supporting Information:**

Supporting Information may be found in the online version of this article.

**Correspondence to:**

X. P. Yuan,
yuanxiaoping@cug.edu.cn

**Citation:**

Yuan, X. P., Guerit, L., Braun, J., Rouby, D., & Shobe, C. M. (2022). Thickness of fluvial deposits records climate oscillations. *Journal of Geophysical Research: Solid Earth*, 127, e2021JB023510. https://doi.org/10.1029/2021JB023510

Received 27 OCT 2021
Accepted 19 MAR 2022

**Author Contributions:**

**Conceptualization:** X. P. Yuan, L. Guerit, J. Braun

**Funding acquisition:** J. Braun, D. Rouby

**Methodology:** X. P. Yuan, J. Braun

**Validation:** X. P. Yuan, L. Guerit, J. Braun, C. M. Shobe

**Writing – original draft:** X. P. Yuan

**Writing – review & editing:** L. Guerit, J. Braun, D. Rouby, C. M. Shobe

© 2022 The Authors.
This is an open access article under
the terms of the Creative Commons
Attribution-NonCommercial License,
which permits use, distribution and
reproduction in any medium, provided the
original work is properly cited and is not
used for commercial purposes.

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

Thickness of Fluvial Deposits Records Climate Oscillations

X. P. Yuan<sup>1,2</sup> <sup>●</sup>, L. Guerit<sup>3,4</sup> <sup>●</sup>, J. Braun<sup>2,5</sup>, D. Rouby<sup>3</sup> <sup>●</sup>, and C. M. Shobe<sup>2,6</sup>

<sup>1</sup>Hubei Key Laboratory of Critical Zone Evolution, School of Earth Sciences, China University of Geosciences, Wuhan, China, <sup>2</sup>Helmholtz Centre Potsdam, German Research Centre for Geosciences (GFZ), Potsdam, Germany, <sup>3</sup>Géosciences Environnement Toulouse, CNRS-IRD-Université de Toulouse, Toulouse, France, <sup>4</sup>Now at University of Rennes, CNRS, Géosciences Rennes, Rennes, France, <sup>5</sup>Institute of Geosciences, University of Potsdam, Potsdam, Germany, <sup>6</sup>Department of Geology and Geography, West Virginia University, Morgantown, WV, USA

**Abstract** Fluvial deposits offer Earth's best-preserved geomorphic record of past climate change over geological timescales. However, quantitatively extracting this information remains challenging in part due to the complexity of erosion, sediment transport and deposition processes and how each of them responds to climate. Furthermore, sedimentary basins have the potential to temporarily store sediments, and rivers subsequently rework those sediments. This may introduce time lags into sedimentary signals and obscure any direct correlation with climate forcing. Here, using a numerical model that combines all three processes—and a new analytical solution—we show that the thickness of fluvial deposits at the outlet of a mountain river can be linked to the amplitude and period of rainfall oscillations but is modulated by the mountain uplift rate. For typical uplift rates of a few mm/yr, climate oscillations at Milankovitch periods lead to alluvial sediment thickness of tens of meters as observed in nature. We also explain the time lag of the order of 20%–25% of the forcing period that is commonly observed between the timing of maximum rainfall and erosion. By comparing to field datasets, our predictions for the thickness and time lag of fluvial deposits are broadly consistent with observations despite the simplicity of our modeling approach. These findings provide a new theoretical framework for quantitatively extracting information on past rainfall variations from fluvial deposits.

**Plain Language Summary** Climate influences the evolution of terrestrial landscapes through the amount of precipitation, which provides water to erode rocks and transport sediment in rivers. At the outlets of mountain ranges, rivers can deposit part of their sediment load; the shape of the deposits is influenced by the amount of flow in the rivers. If the climate changes such that the precipitation rate increases, rivers can cut into their own previous deposits. The remaining deposits are then abandoned above the riverbed. On the contrary, if precipitation decreases, rivers tend to deposit more sediment, leading to increases in the thickness of sediments at the outlets of mountain rivers. Thus, there is a relationship between the amount of precipitations and the thickness of sediments deposited at river outlets. We study this with a computer model that allows us to relate the variations in precipitation rates to variations in thickness of fluvial terrace deposits. This work can be used to better understand how rivers respond to climatic changes, and also to reconstruct climatic variations of the past from observed river deposits.

1. Introduction

Landscapes constantly evolve in response to tectonic and climatic perturbations via fluvial erosion and sediment transport that remove mass from eroding mountain ranges and deposit it in sedimentary basins. Under the influence of changing external controls, rivers can also incise into previously deposited sediments (e.g., Hancock & Anderson, 2002; Humphrey & Heller, 1995; Tofelde et al., 2017). Cycles of sediment aggradation and incision, which can lead for example to the formation of fluvial fill terraces in proximal basins, have been interpreted to reflect variations in climate (Dey et al., 2016; Tofelde et al., 2017), tectonics (Yanites et al., 2010), or base level (Blum & Törnqvist, 2000). Geomorphic cycles as represented by fill terraces are thus often used to constrain the variations of external controls of the past (e.g., Bookhagen et al., 2006; Dey et al., 2016; Tofelde et al., 2017).

However, such aggradation and incision cycles can also arise from autogenic processes (e.g., Malatesta et al., 2017; Yuan, Braun, Guerit, Rouby, & Cordonnier, 2019). Furthermore, drainage systems have the potential to store sediments for tens to hundreds of thousands of years (Castellort & Van Den Driessche, 2003), and rivers subsequently rework those sediments. This may introduce time lags and signal modulation into sedimentary signals and obscure any direct correlation with climate forcing (Duller et al., 2019; Romans et al., 2016;
Tofelde et al., 2021). Consequently, the utility of geomorphic records for quantitatively inferring the amplitude and frequency of climate oscillations, as well as the extent to which climatic signals are filtered by sedimentary systems, remain highly debated (e.g., Armitage et al., 2013; Braun et al., 2015; Castelltort & Van Den Driessche, 2003; Romans et al., 2016; Straub et al., 2020). Furthermore, it is still unclear how climate oscillations are quantitatively preserved in terrestrial landscapes over a range of geological timescales (e.g., thousands of years to tens of millions of years; Braun et al., 2015; Castelltort & Van Den Driessche, 2003; Godard et al., 2013), and whether or not these oscillations always lead to cycles of fluvial aggradation and incision, because of their competing and interacting effects on river discharge, sediment supply, and transport capacity (Dey et al., 2016; Hancock & Anderson, 2002; Tofelde et al., 2017), and the influence of grain size, rock erodibility, and other variables (Armitage et al., 2013).

In recent years, a parameterization of long-term fluvial erosion and sediment transport by rivers (Davy & Lague, 2009) has gained wide acceptance among geomorphologists as it allows the description of fluvial systems as a continuum from detachment-limited behavior (erosion is limited by the capacity of rivers to incise their bed) to transport-limited behavior (erosion is limited by the capacity of rivers to transport their sediment; Braun, 2021; Shobe et al., 2017; Yuan, Braun, Guerit, Rouby, & Cordonnier, 2019; Yuan et al., 2022). In this work, we use this parameterization to study the first-order behaviors of sediment aggradation and incision cycles: their thickness and timing, and how they relate to the amplitude and timing of precipitation variations. In Section 2, we show the parameterization of landscape evolution model, a simple model setup of mountain and basin, and the imposed cyclic precipitation variations. In Section 3, we analyze sediment transport behaviors in basins under cyclic precipitation variations. In Section 4, we use the modeling to relate quantitatively the thickness of fluvial aggradation and incision cycles in basins to the observed amplitude and frequency of precipitation variations in nature. We further suggest a new theoretical framework based on numerical and analytical solutions to this model that may be used to quantitatively extract the amplitude and timing of precipitation variations from the thickness of fluvial fill terrace deposits.

## 2. Method

### 2.1. Fluvial Erosion-Deposition Landscape Evolution Model

The numerical model parameterization rests on the assumption that fluvial erosion is in proportion to the stream power, and that sediment deposition is in proportion to the concentration of sediment in the water column (i.e., the ratio of sediment flux $Q_s$ to water discharge $Q_w$; Davy & Lague, 2009). We use a landscape evolution model which solves the parameterization with an efficient method (Yuan, Braun, Guerit, Rouby, & Cordonnier, 2019). The rate of topographic change $\partial h/\partial t$, in response to tectonic uplift, fluvial erosion and sediment deposition, is given by Yuan, Braun, Guerit, Rouby, and Cordonnier (2019):

$$ \frac{\partial h}{\partial t} = U - K_f \tilde{p}^m A^m S^n + \frac{GQ_s}{\tilde{p}A}, \quad \text{with} \quad Q_s = \int_A \left( U - \frac{\partial h}{\partial t} \right) dA, \qquad (1) $$

where $h$ is elevation, $t$ is time, $U$ is tectonic uplift rate, $K_f$ is erodibility, $A$ is drainage area, $S$ is local slope in the steepest-descent direction of water flow (D8 rule; Braun & Willett, 2013), and $m$ and $n$ are the area and the slope exponents, respectively. Dimensionless $\tilde{p}$ represents any spatial or temporal variation in precipitation $p$ relative to the mean precipitation rate $p_0$. The water discharge is $Q_w = pA$. $G = d^*v_s/p_0$ is a dimensionless deposition coefficient controlling the relative efficiency of deposition versus erosion (Davy & Lague, 2009; Yuan, Braun, Guerit, Rouby, & Cordonnier, 2019), with $d^*$ the dimensionless sediment concentration ratio in the water column and $v_s$ the net settling velocity of sediment. However, such parameters are not explicit in our model, and therefore, $G$ values are derived from morphological analysis at the scale of a landscape (Guerit et al., 2019).

In this work, we use $m = 0.4$ and $n = 1$ to keep the system linear (Hilley et al., 2019; Stock & Montgomery, 1999), and impose $K_f = 2 \times 10^{-5} m^{0.2}/yr$, and $G = 1$ (Yuan, Braun, Guerit, Rouby, & Cordonnier, 2019), close to the estimation of ~1.6 with $n = 1$ based on various natural and experimental examples ($G = 1.6n^{-1.1}$; Guerit et al., 2019). We also explore later the sensitivity of the model to the value of $K_f$, and to varying $G$ from 1 to 2. Our solution scheme for Equation 1 (Yuan, Braun, Guerit, Rouby, & Cordonnier, 2019) is implicit in time and gives first-order accuracy. The potential influence of using higher-accuracy solution schemes is discussed in Section 4.4.
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Simulation of a sediment-routing system and associated riverbed elevation during climate oscillations. Numerical setup with an uplifted domain of 100 × 50 km coupled to a basin of 100 × 50 km. Model resolution: Δx = Δy = 1 km. The landscape evolves by fluvial erosion and deposition, sediment can be stored in the basin and build alluvial deposits (Movie S1). The precipitation rate is constant during 15 Myr and the cyclic variations are then imposed (see text for details).</figcaption></figure>

The landscape evolution model also incorporates a linear diffusion term to approximate hillslope processes (Culling, 1960), that is,

$$ \frac{\partial h}{\partial t} = K_d \nabla^2 h, \qquad (2) $$

where $K_d$ is the hillslope diffusion coefficient. This work assumes a constant diffusion coefficient of $K_d = 0.01 \text{ m}^2/\text{yr}$ (Armitage et al., 2013; Densmore et al., 2007).

## 2.2. Model Setup

We use the model to study the first-order behaviors of sediment aggradation and incision cycles, and the geomorphic records of fluvial deposits using a typical continental sediment-routing system (Densmore et al., 2007; Yuan, Braun, Guerit, Rouby, & Cordonnier, 2019; Carretier et al., 2020), under cyclic precipitation variations. This system comprises a source area where uplift is assumed steady and uniform as well as a proximal basin where uplift is zero and sediments are either deposited or transported out of the model system (Figure 1). The sediment flux out of the system represents the flux that would be transported to a more distal sediment sink like an ocean basin. In the modeled basin, sediment can be stored both within and outside of channels (Carretier et al., 2020). No flexural subsidence is involved in this work, but we test the influence of basin subsidence on modeling results in the Supplementary Information. We also assume for simplicity that the contact between the uplifted domain and the proximal basin edge is a vertical fault (e.g., Armitage et al., 2013; Carretier et al., 2020; Densmore et al., 2007).

The numerical simulations are composed of a region of dimensions 100 × 50 km, uplifting at a constant rate $U$ (mm/yr) and coupled to a basin ($U = 0$) of the same dimensions (Figure 1). The model system is discretized into 10,000 cells with cell size of $\Delta x = \Delta y = 1 \text{ km}$. We also test the influence of model resolutions ($\Delta x = \Delta y = 1 \text{ km}$, 500 m, and 100 m; Figures S2 and S3 in Supporting Information S1) on modeling results. The initial topography has random noise elevation of up to 1 m. The base level, defined by the basin edge, is fixed at an elevation of $h = 0 \text{ m}$ during the simulation. It is an open boundary that enables the export of sediment beyond the model system. The other three boundaries are characterized by a no-flux leaving boundary condition. The model is run for a total time of 20 Myr (20,000 time steps of 1,000 years).

We chose a timestep of 1,000 years and a spatial resolution of 1,000 m to balance model accuracy and efficiency. It would have been possible to use significantly larger timesteps given the stability of the implicit model (Braun &
Willett, 2013; Yuan, Braun, Guerit, Rouby, & Cordonnier, 2019), but model local truncation error increases with timestep such that we choose not to use a larger timestep. Likewise, accuracy could be increased by using smaller timesteps and finer spatial discretizations at the cost of model efficiency. We show in Supporting Information S1 that our results persist for finer model grid resolutions.

## 2.3. Cyclic Variations in Precipitation Rates

The piedmont is the most proximal deposition area outside of the mountain range and it is thus likely to record sedimentary signals without major transformation due to transport and reworking in the floodplain (Romans et al., 2016). The response of the system to cyclic variations in precipitation rates is thus studied at the transition between the mountain range and the basin (i.e., the outlets of mountain rivers; Figure 1).

During a simulation, the initial flat landscape first evolves until it reaches topographic and flux steady states (i.e., the average topography of the system and the average sediment flux out of the domain no longer evolve) under a constant uplift rate ($U = 1$ mm/yr) and constant mean precipitation rate (dimensionless $\tilde{p} = 1$) for $t \le 15$ Myr (Figure 2a). The climatic variations are introduced after 15 Myr by varying the precipitation rate with a sinusoidal perturbation of period $P$ and amplitude $\delta p$. The variation in precipitation relative to the mean precipitation is thus:

$$ \tilde{p} = 1, \quad t \le 15 \text{ Myr}, \\ \tilde{p} = 1 + \delta p \sin(2\pi t/P), \quad t > 15 \text{ Myr}. \qquad (3) $$

The response of topography and sediment fluxes to periodic oscillations in precipitation rate of amplitude $\delta p = 0.5$ is analyzed for different periods $P = 1$ Myr, 400 kyr, and 100 kyr (e.g., Armitage et al., 2013) (upper panels in Figures 2b–2d).

# 3. Results

During the first phase of the simulation ($t \le 15$ Myr), the sediment fluxes exported outside of the mountain range ($Q_{\text{mountain}}$) and eroded in the basin (absolute value of $Q_{\text{basin}}$) increase and reach steady state (Figure 2a). As a consequence, the total sediment flux exported out of the numerical domain ($Q_{\text{system}}$) increases and reaches steady state during this phase (Figure 2a; Movie S1). With the parameterization and for a constant mean precipitation rate (i.e., $\tilde{p} = 1$) and a constant uplift rate of $U = 1$ mm/yr, the steady state is reached within 15 Myr (Figure 2a) with an e-folding response time of $\tau \approx 6$ Myr, consistent with the theoretical response time (Yuan, Braun, Guerit, Rouby, & Cordonnier, 2019) of

$$ \tau = (1 + G/\tilde{p}) U^{1/n-1} K_f^{-1/n} \tilde{p}^{-m/n} k^{-m/n} (1 - bm/n)^{-1} L^{1-bm/n}, \qquad (4) $$

where $L = 50$ km is the width of mountain range, $k \approx 0.5$ and $b \approx 2$ are Hack's law coefficients (Hack, 1957), using $K_f = 2 \times 10^{-5}$, $m^{0.2}/\text{yr}$, $m = 0.4$, $n = 1$, and $G = 1$ (see their definitions in Section 2.1).

## 3.1. Sediment Flux Adjusts Immediately to Climatic Oscillations but Topography Does Not

After reaching steady state, the system responds to the imposed cyclic perturbations in precipitation rate (Figures 2b–2d; Movie S1). Cyclic changes in precipitation rate induce oscillations in the sediment flux out of the mountain belt and out of the numerical domain (black and red curves, respectively, in the middle panels in Figures 2b–2d), and the local sediment flux out of the outlets of mountain rivers ($Q_{\text{outlet}}$; Figure S1 in Supporting Information S1). The model shows that sediment fluxes out of the mountain and out of the numerical domain increase during wet periods (with respect to the reference period) and decrease during dry periods, in general agreement with field observations (e.g., Bookhagen et al., 2005; Fuller et al., 2009). This results from enhanced erosion in the mountains with higher precipitation rates (as the erosion rate is proportional to $\tilde{p}^m$ in Equation 1) together with a decrease in sediment deposition in the basin due to higher water discharge (gray curves in Figures 2b–2d, as the deposition rate is inversely proportional to $\tilde{p}$ in Equation 1).

From the three scenarios explored here, the amplitude of the sediment flux response appears to be independent of the frequency of rainfall oscillations (Figures 2b–2d). The three model-wide sediment fluxes adjust immediately
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Response of the sediment flux and of the average riverbed elevation (see Movie S2) to cyclic precipitation variations, at the outlets of mountain rivers. (a) First phase of the numerical simulation with a constant precipitation rate (i.e., $\delta p = 0$). As the mountain grows and sediment accumulate, the total erosion flux of the numerical domain ($Q_{\text{system}}$, red) increases and reaches steady. The total flux can be divided into the erosion flux in the uplifting region ($Q_{\text{mountain}} = \int_{\text{mountain}} (U - \partial h / \partial t) dA$, black) plus the erosion flux in the basin ($Q_{\text{basin}} = \int_{\text{basin}} (U - \partial h / \partial t) dA$, gray). Similarly, the average elevation along section S-S' (Figure 1) increases and reaches steady. After this first phase, a varying precipitation rate (variation amplitude $\delta p = 0.5$) is imposed with (b) a period of 1 Myr, (c) a period of 400 kyr, and (d) a period of 100 kyr. In each scenario, the sediment fluxes vary in phase with the precipitation rate while the riverbed elevation increases and decreases with a temporal lag. The amplitude of the riverbed elevation variation $\Delta h$ is related to the period of the precipitation variation (see text for details).</figcaption></figure>

to the changes in precipitation rate (Figures 2b–2d), that is, there is no time lag between the forcing (precipitation change) and the response in fluxes (out of the mountain, the basin, and the system). This is because the timescale of the climatic perturbations (i.e., 100 kyr, 400 kyr, and 1 Myr) is small with respect to the response time of the mountain range ($\tau \sim 6$ Myr for these simulations, see Equation 4; Braun et al., 2015; Godard et al., 2013).
<figure><img src="image_6.png" /><figcaption>**Figure 3.** Evolution of the riverbed elevation at the outlet of the mountain range (section S-S' in Figure 1) in response to a precipitation cycle with a period of 1 Myr (Figure 2b). (a) Section at 15.9 Myr showing the aggradation (dry) period between 15.4 and 15.9 Myr. The time interval between light gray lines is 0.1 Myr. (b and c) Sections at 16.1 and 16.2 Myr, respectively, showing the incision (wet) period. Sedimentary deposits of significant thickness are stored above the riverbed as a result of strong contrasts between high channel mobility during aggradation period and low channel mobility during incision period.</figcaption></figure>

However, the modeling results reveal the existence of a time lag between the maximum change in precipitation rate and the maximum incision or aggradation at the outlets of mountain rivers (Figures 2b–2d), corresponding to the local sediment flux out of the outlets of mountain rivers $Q_{\text{outlet}} = \int_{\text{outlet}} (-\partial h/\partial t) dA = 0$ (Figure S1 in Supporting Information S1). This lag occurs because the change in topography (either incision or aggradation) is driven by differences between the perturbed precipitation rate and the steady-state (mean) precipitation rate. In consequence, incision (or aggradation) persists until the precipitation rate is back to its initial value: maximum depth of incision (or thickness of aggradation) records the end of a climatic perturbation, not its climax. This implies that short period perturbations will lead to low total incision/deposition and therefore to limited changes in riverbed elevation. Longer-period climate oscillations are therefore more likely to be recorded in the geomorphic record than shorter ones. We demonstrate that the modeling results are independent of model grid resolutions by testing $\Delta x = \Delta y = 1$ km, 500 m, and 100 m (Figure S2 in Supporting Information S1).

## 3.2. Transient Storage in the Basin Amplifies the Erosion Signal From the Mountain Range

The aggradation and incision cycles are shown by plotting the average riverbed elevation at the transition points through time (the third panels in Figures 2b–2d; Movie S2) and along a cross-section (Figure 3). During wet
<figure><img src="image_7.png" /><figcaption>**Figure 4.** Strong variations in channel mobility between the aggradation and incision phases. Model resolution: $\Delta x = \Delta y = 100$ m. (a) Zoom along one channel in the basin during an aggradation phase (15.4–15.9 Myr), see Movie S3. During this phase, the channel is quite mobile as it moves within the white area. (b) Zoom along the same channel during an incision phase (15.9–16.4 Myr), see Movie S3. The channel is less mobile and almost locked during this phase. This is also illustrated by the relative lateral position of the channel through time at (c) 3 km and (d) 6 km from the mountain front. Blue and yellow areas correspond to wetter and drier periods, respectively, 0 km on the y-axis corresponds to the position of the channel at 15 Myr. Evolution of the riverbed elevation at the outlet of the mountain range is shown in Figure S3 in Supporting Information S1.</figcaption></figure>

periods, the increased water discharge $Q_w$ causes increased erosion efficiency and reduced deposition, both in the mountains and in the basins (black and gray curves in the middle panels in Figures 2b–2d). As a consequence, incision occurs in the basin together with an increase of the sediment flux exported from the numerical domain. On the contrary, during dry periods, erosion rate is reduced in the mountains and sediment is deposited in the basin. In other words, changes in transport capacity (increasing with water discharge $Q_w$) within the basin dominate over changes in sediment supply $Q_s$ from the mountains such that net deposition takes place during drier periods and net incision takes place during wetter periods. This is in good agreement with field observations that episodes of incision, for example within alluvial fans, are associated with wetter periods whereas aggradation phases are associated with drier periods (Bookhagen et al., 2006; Malatesta et al., 2017; Poisson & Avouac, 2004).

We propose that these cycles of aggradation and incision could be recorded in the form of fill terraces in natural settings. As the model does not incorporate channel width dynamics and lateral river erosion processes, it cannot reproduce the details of fill terrace formation *per se*. However, as shown in Figures 2b and 3, and Figure S3 in Supporting Information S1 (higher model resolution), sedimentary deposits of significant thickness (i.e., up to several tens of meters) are stored above the riverbed. This results from the strong contrast in mobility of the riverbed during aggradation or incision phases (Figure 4; Movie S3). During aggradation phases (Figures 4a and 4c), the channel is laterally mobile and the model replicates the behavior of aggrading rivers (Tofelde et al., 2019)—here, lateral channel mobility arises from rapid sediment deposition that changes the flow routing direction. During incision phases (Figures 4b and 4d), the channel is less laterally mobile as it is locked within its banks.
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Modeled thickness of aggradation/incision cycles at the transition points between the uplifted domain and the basin. (a) Sensitivity analyses of oscillation period *P*, amplitude *δp*, uplift rate *U*, and deposition parameter *G* on the thickness *Δh* based on Equation 5. (b) Modeled thicknesses (symbols) compared to the analytical results (curves) for various deposition and precipitation parameters using the uplift rate *U* = 1 mm/yr. The black and red curves represent the estimated thickness from Equation 5 for *G* = 1, *δp* = 0.5 and *G* = 1, *δp* = 0.8, respectively. The difference between the modeled and analytical results for large oscillation periods is partially because (a) we use, in the modeling, the average thickness at the transition points between the uplifted domain and the basin; and (b) several assumptions in the derivation of the analytical result.</figcaption></figure>

Several studies have argued that climate-driven sedimentary signals are likely to be damped when they leave the mountain front (Armitage et al., 2013; Braun et al., 2015) or the sedimentary basin (Métivier & Gaude-mer, 1999; Castelltort & Van Den Driessche, 2003), while others have suggested that climatic-driven fluctuations in sediment flux may be fully transferred to the marine domain (Simpson & Castelltort, 2012). Here, the relative amplitude (i.e., signal to mean) of the climate forcing is not simply transmitted to the basin and out of the model domain but is, in fact, slightly amplified in the sediment flux record (Figures 2b–2d), as observed in previous studies (Godard et al., 2013; Romans et al., 2016; Simpson & Castelltort, 2012). This signal amplification results from the dynamics of the basin, which temporarily stores sediment during drier periods and releases it during wetter periods, leading to phases of aggradation (dry periods) and incision (wet periods) in the basin. The amplitude of these aggradation and incision cycles is greatest at the transition points between the basin and the uplifted domain, in agreement with previous studies (Humphrey & Heller, 1995).

# 4. Discussion

## 4.1. Thickness of Fluvial Deposits as a Quantitative Marker of Climatic Cycles

To quantify the relationship between the aggradation and incision thickness and the forcing amplitude and period, we have derived an approximate analytical solution (Text S1 in Supporting Information S1) for the thickness of aggradation and incision cycles at the transition points between the mountain and the basin (i.e., the outlet of mountain river), in response to a cyclic variation of the precipitation rate:

$$ \Delta h = \frac{PUG(1+m)}{2} - \frac{PUG(1+m)}{\pi \sqrt{1 - \delta p^2/(G+1)^2}} \left[ \frac{\pi}{2} - \tan^{-1} \left( \frac{\delta p/(G+1)}{\sqrt{1 - \delta p^2/(G+1)^2}} \right) \right]. \quad (5) $$
<figure><img src="image_9.png" /><figcaption>**Figure 6.** Parametric studies on sediment aggradation and incision cycles. (a and b) Riverbed elevation as a function of time for various uplift rates *U* (using a constant *K*<sub>f</sub> = 2 × 10<sup>-5</sup> m<sup>0.2</sup>/yr), and for various erodibilities *K*<sub>f</sub> (using a constant uplift rate *U* = 1 mm/yr), respectively. (c and d) Riverbed elevation as a function of time for various mountain lengths (using a constant basin length of 50 km), and for various basin lengths (using a constant mountain length of 50 km), respectively. In (a), the thickness of aggradation and incision cycles is proportional to the uplift rate. In (b and d), the thicknesses of aggradation and incision cycles do not depend on the erodibility *K*<sub>f</sub>, nor on the lengths of the mountain and the basin.</figcaption></figure>

Equation 5 implies that the aggradation and incision thicknesses are proportional to the uplift rate *U* (Figures 5 and 6a), to the amplitude *δp* and the period *P* of the climate oscillations (Figure 5a), but only weakly dependent on the deposition coefficient *G* in the range of 1–2 (Figure 5a), around the estimation of *G* ~1.6 with *n* = 1 (Guerit et al., 2019). The analytical results are consistent with modeled thickness of aggradation and incision cycles based on numerical simulations (Figure 5b). Interestingly, the amplitude of aggradation and incision thickness at the transition points does not depend on the efficiency of fluvial erosion as encompassed in the erodibility parameter *K*<sub>f</sub> (Figure 6b), nor on the lengths of the mountain and the basin (Figures 6c and 6d). The basin subsidence rate affects the average riverbed elevation at the outlets of mountain rivers, but only weakly influences the amplitude of aggradation and incision thickness (Figure S4 in Supporting Information S1). Based on this analytical solution and on the simulation results, the thickness of aggradation and incision preserved in basins may be exploited as a record of past climate oscillations because the thickness (a) is directly proportional to the amplitude and period of climate oscillations, and (b) does not depend on weakly constrained erosional and transport parameters (i.e., *K*<sub>f</sub> and *G*). However, inverting the amplitude of climate oscillations from the preserved sediment thickness requires knowledge of the uplift rate in the adjacent mountain range.

For uplift rates characteristic of tectonically active areas (i.e., a few mm/yr) and a rainfall oscillation of 0.2–0.8, within the reported ranges of variations (~0.1–1) in the Quaternary (Maher & Thompson, 1995), the model predicts that the amplitude of aggradation and incision cycles at Milankovitch periods (20–100 kyr) is of the order of tens of meters. Such a thickness is of the same order of magnitude as typical fill terraces preserved in the field (Dey et al., 2016; Hancock & Anderson, 2002; Tofelde et al., 2017).

## 4.2. Comparison to Observations

To further explore the possibility that the thickness of aggradation and incision cycles records climatic signals, we obtain the sediment aggradation thicknesses and the aggradation durations of fluvial fill terraces along
the Musone River in the Apennines, Italy (Wegmann & Pazzaglia, 2009; Figure 7, Figure S5 and Table S1 in Supporting Information S1). The example contains four paired middle Pleistocene-to-Holocene fill terraces (Qt2, Qt3, Qt4–5, and Qt7) in the frontal part of the mountains and faults (i.e., the Marche and Cingoli Ridges and Thrusts, Figure S5a in Supporting Information S1). The fill terraces of Qt3 and Qt4-5 are approximately at the Musone River outlet in front of the Cingoli Ridge, and the terraces of Qt2 and Qt7 are ~5–10 km downstream from this outlet (Wegmann & Pazzaglia, 2009). The formation of these fill terraces is interpreted to result from climatic cycles since the late Quaternary (Wegmann & Pazzaglia, 2009). The dry and wet periods correspond to the phases of fluvial aggradation and incision, respectively, in agreement with the modeling results. The sediment aggradation durations are relatively well constrained, so we replace the oscillation period P in Equation 5 with the aggradation duration $D_a = P/2$, as observed in the model. The aggradation thickness ($\Delta h$) and the uncertainty of the predicted thickness ($\Delta(\Delta h)$) are

$$ \Delta h = D_a U G(1 + m) - \frac{D_a U G(1 + m)}{\sqrt{1 - \delta p^2/(G + 1)^2}} \left[ 1 - \frac{2}{\pi} \tan^{-1} \left( \frac{\delta p/(G + 1)}{\sqrt{1 - \delta p^2/(G + 1)^2}} \right) \right], \quad (6) $$
and $\Delta(\Delta h) = \sqrt{(\partial \Delta h / \partial D_a \Delta D_a)^2 + (\partial \Delta h / \partial U \Delta U)^2}$,

where $\Delta D_a$ and $\Delta U$ are the uncertainties of aggradation duration and uplift rate, respectively.

The climatic oscillation amplitude $\delta p$ can be extracted from field observation of preserved sediment thickness by comparing to the predicted aggradation thickness using Equation 6. We use $G = 1.5$ to simulate a transport-limited landscape (Guerit et al., 2019) and the constrained uplift rates and their uncertainties ($U$ and $\Delta U$), the aggradation durations and their uncertainties ($D_a$ and $\Delta D_a$) for various oscillation amplitudes $\delta p$ (Figure 7). The climatic aggradation duration ($D_a$) and the uncertainty of aggradation duration ($\Delta D_a$) are suggested in the literature (Table S1 in Supporting Information S1). The uplift rate ($U$) and the uncertainty of uplift rate ($\Delta U$) can be obtained as follows. Low-temperature thermochronometry data show long-term exhumation rates of 1–2 mm/yr in the Apennines (Balestrieri et al., 2003; Simoni et al., 2003). These exhumation rates can be assumed as a proxy for rock uplift rates, within the observed rock uplift rates of 1–3 mm/yr along the Apennines (Devoti et al., 2011; Serpelloni et al., 2013). Thus, we obtain the uplift rate of $U = 2 \pm 1$ mm/yr (the average uplift rate of $U = 2$ mm/yr and the uncertainty of $\Delta U = 1$ mm/yr) for this area. Despite the relatively large uncertainties on thickness predictions arising from uncertainties $\Delta U$, the model and observed values are broadly consistent with each other for a proper rainfall oscillation amplitude $\delta p$. Fitting the observed thickness in the Apennines requires an amplitude of rainfall variations of $\delta p \sim 0.6–0.8$ (Figure 7).

We present the above example of field data that might show deposits consistent with our model, but our comparison is tentative and preliminary because there exist many complexities in the field examples, and our modeling setup requires comparison against observations of fill terraces at the outlets of mountain rivers. These datasets are rarely found in the literature. For further progress in model-data comparison, it would be valuable to have more field measurements of sediment aggradation thickness and duration (i.e., the initial and final timings of deposits) of fluvial fill terraces at the outlets of mountain rivers, and constrained uplift rates in the source areas.

### 4.3. Time Lag Between Maximum Rainfall and Erosion

The model predicts a time lag between maximum rainfall at $P/4$ and erosion at the outlets of mountain rivers (Figures 2b and 8a) due to the fact that changes in topography are driven by the change in precipitation rate with respect to the initial steady-state value. Consequently, incision persists until the precipitation rate returns to its initial value (i.e., $\tilde{p} = 1$) at $\sim P/2$ (Figure 2b). Thus, the time lag between maximum rainfall and erosion is empirically close to 25% of the oscillation period $P$. We further perform a number of simulations with various oscillation periods after steady-state conditions, and monitor the associated timings of maximum erosion at the outlets of mountain rivers. We find that, when normalized by the oscillation period, this time lag decreases from 25% to a few percent in a predictable manner as a function of the oscillation period normalized by the response time of the mountain range ($P/\tau$ increases from 0.001 to 10; Figure 9). This time lag is of the order of 20%–25% of the forcing period at typical Milankovitch periods (Figure 9) assuming an uplift rate of the order of 1 mm/yr, typical of active orogenic settings. Such a time lag has been documented in the NW Himalayas (Bookhagen
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Predicted over observed aggradation thicknesses of the fill terraces along the Musone River in the Apennines. (a–d) Comparisons with $U = 2 \pm 1$ mm/yr and with $\delta p = 0.2, 0.4, 0.6$, and $0.8$, respectively. Best fit is obtained for $\delta p \sim 0.6–0.8$ (see text for details). Observed datasets are from Figure S5 and Table S1 in Supporting Information S1.</figcaption></figure>

et al., 2006; Dey et al., 2016), in Western Colorado (Duller et al., 2019; Foreman et al., 2012), and in Northern Spain (Duller et al., 2019), with a typical offset of several thousand to tens of thousands of years (Figures 8 and 9; Table S2 in Supporting Information S1). For the examples in the NW Himalayas, the time lag of 4–5 kyr is obtained between the maximum rainfall and the minimum riverbed elevation (Figures 8b and 8c; Table S2 in Supporting Information S1).

When the period of rainfall oscillation is large relative to the orogenic response timescale (e.g., $P/\tau > 0.1$), the time lag decreases but remains on the order of several to 20 percent of the forcing period. This implies that the timing of aggradation cessation deduced from cosmogenic isotope exposure data is likely to lag behind by a measurable, and potentially predictable, amount of time from the maximum precipitation rate. The results help explain, for example, why the minimum riverbed elevation at ~3 ka in the NW Himalayas (Bookhagen et al., 2006; Dey et al., 2016; Figures 8b and 8c) lags behind the mid-Holocene (i.e., 7–8 ka) climatic optimum. The observed time lag of 4–5 kyr (22%–28% of the forcing period, Table S2 in Supporting Information S1) is consistent with the predicted time lag (~25%) between maximum rainfall and maximum topographic response.

### 4.4. Model Limitations

Several limitations preclude precise comparisons between model predictions and field observations. First, we use a simple sediment-routing system evolving under constant and uniform controlling factors (i.e., uplift rate, amplitude of climate oscillations). An increase of one of these controlling factors (e.g., $U$ and $\delta p$) during model simulation will increase the amplitude of incision cycles, and may lead to the destruction of former geomorphic records.
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Modeled time lag and inferred time lags for three different natural systems. (a) The modeled time lag from Figure 2b. (b) Along the Sutlej valley (NW Himalayas), we propose that the minimum riverbed elevation (observed ~3 ka) lags behind the mid-Holocene climatic optimum (7–8 ka) by 4–5 kyr (after Bookhagen et al., 2006). (c) The time lag of ~4 ka can be proposed in the Kangra valleys (NW Himalayas) (after Dey et al., 2016).</figcaption></figure>

Thus, field records are likely recording stronger-than-average changes of the amplitude of the controlling factors (Dott, 1983). Second, the model captures the changing fluvial erosion dynamics at the transition between aggrada- tion and incision cycles, thus the preserved sediment thickness may be comparable to nature. However, the averaged, long-term model cannot explicitly simulate water flow, channel geometry (width), or lateral erosion/ channel migration, so it does not simulate fill terrace formation *per se*. Models of shallow-water flow dynamics (e.g., Coulthard et al., 2013; Davy et al., 2017) and accompanying erosion laws may overcome these issues. Third, for the sake of simplicity, we did not take the threshold of erosion (or critical shear stress) required to entrain sediment into account, because it depends on various parameters, such as grain size, relative roughness, and discharge variability (e.g., Buffington & Montgomery, 1997; DiBiase & Whipple, 2011). The threshold of erosion, which spans a relatively large range and is difficult to constrain over landscape evolution timescales (Theodoratos & Kirchner, 2020), likely influences model behavior. Fourth, our modeling uses the D8 rule: flow

is routed into the steepest-descent cell of the eight surrounding cells (Braun & Willett, 2013). Future studies can consider how channel dynamics might change if more continuous (e.g., route-to-many) flow routing algorithms are used (e.g., Shelef & Hilley, 2013). This could be important when studying the detailed landscape dynamics at the mountain front, including the modeled lateral channel mobility (Yuan, Braun, Guerit, Rouby, & Cordonnier, 2019). Fifth, modeled changes in channel path (avulsions) might also depend on the accuracy of the model solution scheme. Avulsions might arise from numerical “over-deposition” when the model deposits more sediment than it would with a more accurate solution method. Using higher-accuracy explicit approaches or adaptive time-stepping methods might therefore change the frequency of modeled avulsions and therefore affect lateral channel mobility. Last but not least, our simple approach may miss some complicated changes in the hydrologic and geomorphic processes that accompany climate oscil- lations. For example, outwash terraces produced by glacial recession will involve increased discharge due to ice melting, as well as the liberation of hillslope material due to accelerated mass wasting in the wake of glacial recession (e.g., Moon et al., 2012). Also, the deposition coefficient *G* in our modeling is proportional to sediment settling velocity, related to parti- cle size in transport. We use a constant *G* value in each simulation, which implicitly assumes a constant grain-size distribution of sediment delivered

<figure><img src="image_12.png" /><figcaption>**Figure 9.** Time lag $\theta$ between the maximum precipitation rate and the maximum riverbed erosion with respect to the period of climate oscillations *P*. The time lag $\theta$ in the y-axis is normalized by the period of climatic oscillations *P*, and the period *P* in the x-axis is normalized by the response time of the system ($\tau$; Equation 4). Black and red curves represent the time lag predicted by the modeling for the erodibilities of $K_f = 2 \times 10^{-5} m^{0.2}/yr$ and $1 \times 10^{-4} m^{0.2}/yr$, respectively. The gray box indicates a system with a response time of 1–10 Myr submitted to typical Milankovitch climatic variations of 20–100 kyrs.</figcaption></figure>
to channels. The distribution of upland processes might cause *G* to change with precipitation as the grain-size distribution delivered to channels may change during glacial-interglacial fluctuations. Further work can investigate the complicated changes by introducing variable parameters on *K*<sub>d</sub> (which controls on the amount of hillslope material delivered to channels) and *G* (which sets the sediment transport and deposition behaviors) due to climate oscillations.

## 5. Conclusions

This work explores and quantifies potential physical links between climate oscillations and fluvial sediment aggradation and incision cycles at the outlets of mountain rivers. In contrast to autogenic sediment signals observed far from the source area with the same model (Yuan, Braun, Guerit, Rouby, & Cordonnier, 2019), we document a clear link between modeled climate oscillations and sediment signals in proximal basins. In particular, the thickness of fluvial aggradation is proportional to the amplitude and period of variations in precipitation rate but is modulated by the source area uplift rate, which partly controls the topography of the source and its rate of erosion. The model predictions for sediment thickness are broadly consistent with field datasets; the model may be useful for estimating the amplitude of climatic forcing from observed aggradation phases in basins. This is particularly relevant for fill terraces located at the transition between a mountain range and its adjacent basin. Modeling results also help explain the time lag between the maximum change in rainfall and the maximum change in topography as observed in nature.

The modeling also shows that, unlike diffusion models (Métivier & Gaudemer, 1999; Castelltort & Van Den Driessche, 2003), variations in sediment flux resulting from short-period (tens to hundreds of thousands of years) climate-driven variations in source area erosion are not buffered (Braun, 2021) but on the contrary slightly amplified as they cross basins. Marine sedimentary basins may therefore record sediment flux cycles resulting from climate oscillations (Simpson & Castelltort, 2012; Yuan, Braun, Guerit, Simon, et al., 2019). Our results provide a framework to help untangle the complexity of erosion, sediment transport and deposition processes, how each of them is affected by climate, and how past climate variations are recorded in the geomorphic record.

## Data Availability Statement

The codes used for the simulations are available in Yuan, Braun, Guerit, Rouby, and Cordonnier (2019) and https://doi.org/10.5281/zenodo.3833983 website (Bovy & Braun, 2020). Figures were made using ParaView, InkScape, and Matplotlib.

## Acknowledgments

The work is part of the COLORS project funded by TOTAL. C.M.S. received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie Grant agreement number 833132. The authors thank John Armitage and Rob Duller for their feedbacks on an early draft. The authors also thank François Guillocheau, Joel Scheingross, Kimberly Huppert, Stefanie Tofelde and Taylor Schildgen for their helpful discussions. The authors are grateful to the AE Peter G. DeCelles, and two reviewers George Hilley and Chris Paola for their very constructive comments. Open access funding enabled and organized by Projekt DEAL.

## References

Armitage, J. J., Jones, T. D., Duller, R. A., Whittaker, A. C., & Allen, P. A. (2013). Temporal buffering of climate-driven sediment flux cycles by transient catchment response. *Earth and Planetary Science Letters*, **369**, 200–210. https://doi.org/10.1016/j.epsl.2013.03.020

Balestrieri, M., Bernet, M., Brandon, M. T., Picotti, V., Reiners, P., & Zattin, M. (2003). Pliocene and Pleistocene exhumation and uplift of two key areas of the Northern Apennines. *Quaternary International*, **101**, 67–73. https://doi.org/10.1016/s1040-6182(02)00089-7

Blum, M. D., & Törnqvist, T. E. (2000). Fluvial responses to climate and sea-level change: A review and look forward. *Sedimentology*, **47**, 2–48. https://doi.org/10.1046/j.1365-3091.2000.00008.x

Bookhagen, B., Fleitmann, D., Nishiizumi, K., Strecker, M. R., & Thiede, R. C. (2006). Holocene monsoonal dynamics and fluvial terrace formation in the northwest Himalaya, India. *Geology*, **34**(7), 601–604. https://doi.org/10.1130/g22698.1

Bookhagen, B., Thiede, R. C., & Strecker, M. R. (2005). Late Quaternary intensified monsoon phases control landscape evolution in the north-west Himalaya. *Geology*, **33**(2), 149–152. https://doi.org/10.1130/g20982.1

Bovy, B., & Braun, J. (2020). *Fastscape-lem/fastscapelib-fortran*. Zenodo. https://doi.org/10.5281/zenodo.3833983

Braun, J. (2021). Comparing the transport-limited and ξ-q models for sediment transport. *Earth Surface Dynamics Discussions*, **1–39**.

Braun, J., Voisin, C., Gourlan, A., & Chauvel, C. (2015). Erosional response of an actively uplifting mountain belt to cyclic rainfall variations. *Earth Surface Dynamics*, **3**(1), 1–14. https://doi.org/10.5194/esurf-3-1-2015

Braun, J., & Willett, S. D. (2013). A very efficient O(n), implicit and parallel method to solve the stream power equation governing fluvial incision and landscape evolution. *Geomorphology*, **180**, 170–179. https://doi.org/10.1016/j.geomorph.2012.10.008

Buffington, J. M., & Montgomery, D. R. (1997). A systematic analysis of eight decades of incipient motion studies, with special reference to gravel-bedded rivers. *Water Resources Research*, **33**(8), 1993–2029. https://doi.org/10.1029/96wr03190

Carretier, S., Guerit, L., Harries, R., Regard, V., Maffre, P., & Bonnet, S. (2020). The distribution of sediment residence times at the foot of mountains and its implications for proxies recorded in sedimentary basins. *Earth and Planetary Science Letters*, **546**, 116448. https://doi.org/10.1016/j.epsl.2020.116448

Castelltort, S., & Van Den Driessche, J. (2003). How plausible are high-frequency sediment supply-driven cycles in the stratigraphic record? *Sedimentary Geology*, **157**(1–2), 3–13. https://doi.org/10.1016/s0037-0738(03)00066-6
Coulthard, T. J., Neal, J. C., Bates, P. D., Ramirez, J., Almeida, G. A., & Hancock, G. R. (2013). Integrating the LISFLOOD-FP 2D hydrodynamic model with the CAESAR model: Implications for modelling landscape evolution. *Earth Surface Processes and Landforms, 38*(15), 1897–1906. https://doi.org/10.1002/esp.3478

Culling, W. (1960). Analytical theory of erosion. *The Journal of Geology*, **68**(3), 336–344. https://doi.org/10.1086/626663

Davy, P., Croissant, T., & Lague, D. (2017). A precipitation method to calculate river hydrodynamics, with applications to flood prediction, landscape evolution models, and braiding instabilities. *Journal of Geophysical Research: Earth Surface*, **122**(8), 1491–1512. https://doi.org/10.1002/2016jf004156

Davy, P., & Lague, D. (2009). Fluvial erosion/transport equation of landscape evolution models revisited. *Journal of Geophysical Research*, **114**, F03007. https://doi.org/10.1029/2008JF001146

Densmore, A. L., Allen, P. A., & Simpson, G. (2007). Development and response of a coupled catchment fan system under changing tectonic and climatic forcing. *Journal of Geophysical Research*, **112**(F1). https://doi.org/10.1029/2006JF000474

Devoti, R., Esposito, A., Pietrantonio, G., Pisani, A. R., & Riguzzi, F. (2011). Evidence of large scale deformation patterns from GPS data in the Italian subduction boundary. *Earth and Planetary Science Letters*, **311**(3–4), 230–241. https://doi.org/10.1016/j.epsl.2011.09.034

Dey, S., Thiede, R. C., Schildgen, T. F., Wittmann, H., Bookhagen, B., Scherler, D., et al. (2016). Climate-driven sediment aggradation and incision since the late Pleistocene in the NW Himalaya, India. *Earth and Planetary Science Letters*, **449**, 321–331. https://doi.org/10.1016/j.epsl.2016.05.050

DiBiase, R. A., & Whipple, K. X. (2011). The influence of erosion thresholds and runoff variability on the relationships among topography, climate, and erosion rate. *Journal of Geophysical Research*, **116**(F4). https://doi.org/10.1029/2011JF002095

Dott, R. H. (1983). Episodic sedimentation; how normal is average? How rare is rare? Does it matter? *Journal of Sedimentary Research*, **53**(1), 5–23. https://doi.org/10.1306/212f8148-2b24-11d7-8648000102c1865d

Duller, R. A., Armitage, J. J., Manners, H. R., Grimes, S., & Jones, T. D. (2019). Delayed sedimentary response to abrupt climate change at the Paleocene–Eocene boundary, northern Spain. *Geology*, **47**(2), 159–162. https://doi.org/10.1130/g45631.1

Foreman, B. Z., Heller, P. L., & Clementz, M. T. (2012). Fluvial response to abrupt global warming at the Palaeocene/Eocene boundary. *Nature*, **491**(7422), 92–95. https://doi.org/10.1038/nature11513

Fuller, T. K., Perg, L. A., Willenbring, J. K., & Lepper, K. (2009). Field evidence for climate-driven changes in sediment supply leading to strath terrace formation. *Geology*, **37**(5), 467–470. https://doi.org/10.1130/g25487a.1

Godard, V., Tucker, G. E., Burch Fisher, G., Burbank, D. W., & Bookhagen, B. (2013). Frequency-dependent landscape response to climatic forcing. *Geophysical Research Letters*, **40**(5), 859–863. https://doi.org/10.1002/grl.50253

Guerit, L., Yuan, X.-P., Carretier, S., Bonnet, S., Rohais, S., Braun, J., & Rouby, D. (2019). Fluvial landscape evolution controlled by the sediment deposition coefficient: Estimation from experimental and natural landscapes. *Geology*, **47**(9), 853–856. https://doi.org/10.1130/g46356.1

Hack, J. T. (1957). *Studies of Longitudinal Stream Profiles in Virginia and Maryland* (Vol. 294). US Government Printing Office.

Hancock, G. S., & Anderson, R. S. (2002). Numerical modeling of fluvial strath-terrace formation in response to oscillating climate. *Geological Society of America Bulletin*, **114**(9), 1131–1142. https://doi.org/10.1130/0016-7606(2002)114<1131:nmofst>2.0.co;2

Hilley, G. E., Porder, S., Aron, F., Baden, C. W., Johnstone, S. A., Liu, F., et al. (2019). Earth's topographic relief potentially limited by an upper bound on channel steepness. *Nature Geoscience*, **12**, 828–832. https://doi.org/10.1038/s41561-019-0442-3

Humphrey, N. F., & Heller, P. L. (1995). Natural oscillations in coupled geomorphic systems: An alternative origin for cyclic sedimentation. *Geology*, **23**(6), 499–502. https://doi.org/10.1130/0091-7613(1995)023<0499:noicgs>2.3.co;2

Maher, B. A., & Thompson, R. (1995). Paleorainfall reconstructions from pedogenic magnetic susceptibility variations in the Chinese loess and paleosols. *Quaternary Research*, **44**(3), 383–391. https://doi.org/10.1006/qres.1995.1083

Malatesta, L. C., Prancevic, J. P., & Avouac, J.-P. (2017). Autogenic entrenchment patterns and terraces due to coupling with lateral erosion in incising alluvial channels. *Journal of Geophysical Research: Earth Surface*, **122**(1), 335–355. https://doi.org/10.1002/2015jf003797

Métivier, F., & Gaudemer, Y. (1999). Stability of output fluxes of large rivers in South and East Asia during the last 2 million years: Implications on floodplain processes. *Basin Research*, **11**(4), 293–303. https://doi.org/10.1046/j.1365-2117.1999.00101.x

Moon, T., Joughin, I., Smith, B., & Howat, I. (2012). 21st-century evolution of Greenland outlet glacier velocities. *Science*, **336**(6081), 576–578. https://doi.org/10.1126/science.1219985

Poisson, B., & Avouac, J.-P. (2004). Holocene hydrological changes inferred from alluvial stream entrenchment in north Tian Shan (northwestern China). *The Journal of Geology*, **112**(2), 231–249. https://doi.org/10.1086/381659

Romans, B. W., Castellort, S., Covault, J. A., Fildani, A., & Walsh, J. (2016). Environmental signal propagation in sedimentary systems across timescales. *Earth-Science Reviews*, **153**, 7–29. https://doi.org/10.1016/j.earscirev.2015.07.012

Serpelloni, E., Faccenna, C., Spada, G., Dong, D., & Williams, S. D. (2013). Vertical GPS ground motion rates in the Euro-Mediterranean region: New evidence of velocity gradients at different spatial scales along the Nubia-Eurasia plate boundary. *Journal of Geophysical Research: Solid Earth*, **118**(11), 6003–6024. https://doi.org/10.1002/2013jb010102

Shelef, E., & Hilley, G. E. (2013). Impact of flow routing on catchment area calculations, slope estimates, and numerical simulations of landscape development. *Journal of Geophysical Research: Earth Surface*, **118**(4), 2105–2123. https://doi.org/10.1002/jgrf.20127

Shobe, C. M., Tucker, G. E., & Barnhart, K. R. (2017). The SPACE 1.0 model: A Landlab component for 2-D calculation of sediment transport, bedrock erosion, and landscape evolution. *Geoscientific Model Development*, **10**(12), 4577–4604. https://doi.org/10.5194/gmd-10-4577-2017

Simoni, A., Elmi, C., & Picotti, V. (2003). Late Quaternary uplift and valley evolution in the Northern Apennines: Lamone catchment. *Quaternary International*, **101**, 253–267. https://doi.org/10.1016/s1040-6182(02)00106-4

Simpson, G., & Castellort, S. (2012). Model shows that rivers transmit high-frequency climate cycles to the sedimentary record. *Geology*, **40**(12), 1131–1134. https://doi.org/10.1130/g33451.1

Stock, J. D., & Montgomery, D. R. (1999). Geologic constraints on bedrock river incision using the stream power law. *Journal of Geophysical Research*, **104**(B3), 4983–4993. https://doi.org/10.1029/98jb02139

Straub, K. M., Duller, R. A., Foreman, B. Z., & Hajek, E. A. (2020). Buffered, incomplete, and shredded: The challenges of reading an imperfect stratigraphic record. *Journal of Geophysical Research: Earth Surface*, **125**(3), e2019JF005079. https://doi.org/10.1029/2019JF005079

Theodoratos, N., & Kirchner, J. W. (2020). Dimensional analysis of a landscape evolution model with incision threshold. *Earth Surface Dynamics*, **8**(2), 505–526. https://doi.org/10.5194/esurf-8-505-2020

Tofelde, S., Bernhardt, A., Guerit, L., & Romans, B. W. (2021). Times associated with source-to-sink propagation of environmental signals during landscape transience. *Frontiers in Earth Science*, **9**, 227. https://doi.org/10.3389/feart.2021.628315

Tofelde, S., Savi, S., Wickert, A. D., Bufe, A., & Schildgen, T. (2019). Alluvial channel response to environmental perturbations: Fill-terrace formation and sediment-signal disruption. *Earth Surface Dynamics*, **7**, 609–631. https://doi.org/10.5194/esurf-7-609-2019
Tofelde, S., Schildgen, T. F., Savi, S., Pingel, H., Wickert, A. D., Bookhagen, B., et al. (2017). 100 kyr fluvial cut-and-fill terrace cycles since the Middle Pleistocene in the southern Central Andes, NW Argentina. *Earth and Planetary Science Letters, 473*, 141–153. https://doi.org/10.1016/j.epsl.2017.06.001

Wegmann, K. W., & Pazzaglia, F. J. (2009). Late Quaternary fluvial terraces of the Romagna and Marche Apennines, Italy: Climatic, lithologic, and tectonic controls on terrace genesis in an active orogen. *Quaternary Science Reviews*, 28(1–2), 137–165. https://doi.org/10.1016/j.quascirev.2008.10.006

Yanites, B. J., Tucker, G. E., Mueller, K. J., & Chen, Y.-G. (2010). How rivers react to large earthquakes: Evidence from central Taiwan. *Geology, 38*(7), 639–642. https://doi.org/10.1130/g30883.1

Yuan, X. P., Braun, J., Guerit, L., Rouby, D., & Cordonnier, G. (2019). A new efficient method to solve the stream power law model taking into account sediment deposition. *Journal of Geophysical Research: Earth Surface, 124*(6), 1346–1365. https://doi.org/10.1029/2018jf004867

Yuan, X. P., Braun, J., Guerit, L., Simon, B., Bovy, B., Rouby, D., et al. (2019). Linking continental erosion to marine sediment transport and deposition: A new implicit and O(N) method for inverse analysis. *Earth and Planetary Science Letters*, 524, 115728. https://doi.org/10.1016/j.epsl.2019.115728

Yuan, X. P., Huppert, K., Braun, J., Shen, X., Liu-Zeng, J., Guerit, L., et al. (2022). Propagating uplift controls on high-elevation, low-relief landscape formation in the southeast Tibetan Plateau. *Geology*, 50(1), 60–65. https://doi.org/10.1130/g49022.1