the Creative Commons Attribution 3.0 License.

the Creative Commons Attribution 3.0 License.

# Numerical modelling of landscape and sediment flux response to precipitation rate change

### Alexander C. Whittaker

### Mustapha Zakari

### Benjamin Campforts

Laboratory-scale experiments of erosion have
demonstrated that landscapes have a natural (or intrinsic) response time to a
change in precipitation rate. In the last few decades there has been growth
in the development of numerical models that attempt to capture landscape
evolution over long timescales. However, there is still an uncertainty
regarding the
validity of the basic assumptions of mass transport that are made in deriving
these models. In this contribution we therefore return to a principal
assumption of sediment transport within the mass balance for surface
processes; we explore the sensitivity of the classic end-member landscape
evolution models and the sediment fluxes they produce to a change in
precipitation rates. One end-member model takes the mathematical form of a
kinetic wave equation and is known as the stream power model, in which sediment
is assumed to be transported immediately out of the model domain. The second
end-member model is the transport model and it takes the form of a diffusion
equation, assuming that the sediment flux is a function of the water flux and
slope. We find that both of these end-member models have a response time that
has a proportionality to the precipitation rate that follows a negative power
law. However, for the stream power model the exponent on the water flux term
must be less than one, and for the transport model the exponent must be
greater than one, in order to match the observed concavity of natural
systems. This difference in exponent means that the transport model generally
responds more rapidly to an increase in precipitation rates, on the order of
10^{5} years for post-perturbation sediment fluxes to return to within
50 % of their initial values, for theoretical landscapes with a scale of
100×100 km. Additionally from the same starting conditions, the
amplitude of the sediment flux perturbation in the transport model is
greater, with much larger sensitivity to catchment size. An important finding
is that both models respond more quickly to a wetting event than a drying
event, and we argue that this asymmetry in response time has significant
implications for depositional stratigraphies. Finally, we evaluate the extent
to which these constraints on response times and sediment fluxes from simple
models help us understand the geological record of landscape response to
rapid environmental changes in the past, such as the Paleocene–Eocene thermal
maximum (PETM). In the Spanish Pyrenees, for instance, a relatively rapid (10
to 50 kyr) duration of the deposition of gravel is observed for a climatic shift
that is thought to be towards increased precipitation rates. We suggest that the
rapid response observed is more easily explained through a diffusive
transport model because (1) the model has a faster response time, which is consistent
with the documented stratigraphic data, (2) there is a high-amplitude spike
in sediment flux, and (3) the assumption of instantaneous transport is
difficult to justify for the transport of large grain sizes as an alluvial
bedload. Consequently, while these end-member models do not reproduce all
the complexity of processes seen in real landscapes, we argue that variations
in long-term erosional dynamics within source catchments can fundamentally
control when, how, and where sedimentary archives can record past
environmental change.

How river networks form and how landscapes erode remains a basic research
question despite more than a century of experimentation and study. At a
fundamental level, the root of the problem is a lack of an equation of motion
for erosion derived from first principles (Dodds and Rothman, 2000). A range
of heuristic erosion equations have, however, been proposed from stochastic
models (Banavar et al., 1997; Pastor-Satorras and Rothman, 1998) to
deterministic models based on the St. Venant shallow water equations
(Izumi and Parker, 1995; Smith, 2010; Smith and Bretherton, 1972), diffusive
*transport-limited* conditions (Whipple and Tucker, 2002), or the
stream power law (e.g. Howard and Kerby, 1983; Whipple and Tucker, 2002; Willett et al., 2014,
among many others). These models, in various forms, have been explored to try
to understand how landscape evolves and responds to tectono-environmental
change. In general terms, numerical studies have found that landscapes
typically recover from a shift in tectonic uplift after 10^{5} to 10^{6}
years (reviewed in Romans et al., 2016). These apparently long
response timescales to tectonic perturbations have been supported by field
observations of landscapes upstream of active faults
(Cowie et al., 2008; Whittaker and Boulton, 2012; Whittaker et al., 2007), although
the precise appropriateness of any time-integrated erosion law to specific
field sites is not always easy to establish. Sediment flux response times for
the advective stream power law have been previously characterized by
Whipple (2001) and Baldwin et al. (2003), and for the transport model
they have been studied by Armitage et al. (2011) and
Armitage et al. (2013), but not systematically or using 2-D models.
Furthermore, to our knowledge no comparison between the transport models has
been previously made.

The response of landscapes and sediment routing systems to a change in the
magnitude or timescale of precipitation rates is expected to depend on the
long-term erosion law implemented
(Armitage et al., 2011, 2013; Castelltort and Van Den Dreissche, 2003).
Some numerical modelling studies based on treating erosion as a length-dependent diffusive problem suggest that landscape responses to a change in
rainfall are also on the order of 10^{5} to 10^{6} years, similar to
tectonic perturbations, although they produce diagnostically different
stratigraphic signatures from the latter
(Armitage et al., 2011). However, other modelling
contributions with different assumptions suggest that response times to a
precipitation change may be more rapid (Simpson and Castelltort, 2012), although field
data sets remain equivocal (see Demoulin et al., 2017, for a recent
review). In laboratory studies, a series of experiments in which granular piles
of a length scale of the order of centimetres are eroded due to surface water
have demonstrated that a change in precipitation rate leads a period of
adjustment of the landscape topography until a new steady state is achieved
(Bonnet and Crave, 2003; Rohais et al., 2011). These experiments use a mixture
of granular silica of a mean diameter between 10 and 20 µm that
is eroded by water released from a fine sprinkler system above. Given the
complexity of these experiments, unfortunately there have been insufficient
different precipitation rates studied to fully understand how the recovery
timescale varies as a function of precipitation or other parameters.

It has been increasingly recognized over the last 2 decades that many basic geomorphic measures of catchments, such as the scaling between channel slopes and catchment drainage areas, are typically unable to distinguish the erosional processes behind their formation (Dodds and Rothman, 2000; Tinkler and Whol, 1998; Tucker and Whipple, 2002; Whipple, 2004). Erosion and transport can be described by equations that encompass both advective and diffusive processes (Smith and Bretherton, 1972) and at topographic steady state, it is very well established that fluvial erosion models based on either of these two end-members can produce very similar river longitudinal profiles (Tucker and Whipple, 2002; van der Beek and Bishop, 2003).

Non-uniqueness or equifinality is a common problem when comparing the morphology generated from landscape evolution models (Hancock et al., 2016). Consequently, we aim to explore whether the sediment flux responses of fluvial systems to a precipitation perturbation may be diagnostically different for the two end-member deterministic models across a range of parameter space. This issue is pertinent because within sedimentary basins, a change in the erosional dynamics upstream could be recorded by changes in the total sediment volumes stored in sedimentary basins (Allen et al., 2013; Michael et al., 2014), in sediment delivery or sediment accumulation rates linked to landscape response times (Armitage et al., 2015; Foreman et al., 2012), and/or in the grain sizes deposited as a function of sediment flux output (Armitage et al., 2011; D'Arcy et al., 2016; Paola et al., 1992; Whittaker et al., 2011).

In this article we make a comparative study between the transport and stream power model to further explore the potential differences between these two end-member hypothetical landscape evolution models. We will focus on the transient period of adjustment to a perturbation in precipitation rates, and using end-member numerical models we attempt to evaluate how the response time varies as a function of the model forcing. To this end we aim to find the model parameters that generate similar landscape morphologies such that we can subsequently explore how the same end-member models respond to a change in precipitation rate. We believe that the results of this study have implications for understanding the responses of landscapes to past changes in climate and could potentially be compared with and tested against further laboratory experiments.

## 2.1 Erosion within a single dimension system

We aim to understand the effects of the most basic assumptions of mass
transport in landscape evolution on the sediment flux record. In other words,
how do the response times vary for the advective stream power law and the
diffusive transport model? To this end we derive the two models from first
principles to demonstrate clearly how, from the same starting point, the
fundamental assumptions made about mass transport initially give rise to very
different model equations. We use this framework as a context for our
investigation of an eroding system responding to precipitation change. We
first define a one-dimensional system from which the basic equations can be
assembled. Following Dietrich et al. (2003) we define a landscape of
elevation *z* composed of bedrock, thickness *η* (units of metres), and a
surface layer of sediment with thickness *h* (units of metres; see
Fig. 1). This landscape is forced externally through uplift
rate *U* (units of m yr^{−1}). The bedrock is transferred into sediment
through erosion at a rate *E* (units of m yr^{−1}) and the sediment is
transported across the system with a flux *q*_{s} (units of
m^{2} yr^{−1}). Assuming that the densities of sediment and bedrock are
equal, then the change in bedrock thickness is

and the rate of change in sediment thickness is

It then follows that the rate of change in landscape elevation is

It is important to realize that to solve Eq. (3), we are required to make some assumptions that fundamentally affect the erosional dynamics of the modelled system, and we illustrate this below.

One basic assumption to make is that there is always a supply of transportable sediment; we can then follow through with the summation in Eq. (3), giving

This may be appropriate when modelling the transport of sediment along the
riverbed and when considering the formation of alluvial fans
(Guerit et al., 2014; Paola et al., 1992; Whipple and Tucker, 2002). In the absence
of surface water we can assume that sediment flux is simply a function of
local slope ${q}_{\mathrm{s}}=-\mathit{\kappa}{\partial}_{x}z$, where *κ* is the hill
slope diffusion coefficient. In the presence of flowing water, the
sediment flux is a function of the flowing water and local slope
${q}_{\mathrm{s}}=-c{q}_{\mathrm{w}}^{\mathit{\delta}}{\left({\partial}_{x}z\right)}^{\mathit{\gamma}}$,
where *c* is the transport coefficient (units (m^{2} yr^{−1})^{1−n}),
*q*_{w} is the water flux per unit width (units m^{2} yr^{−1}),
and the exponents *δ*>1 and *γ*≥1 are dependent on how
sediment grains are transported along the bed
(Paola et al., 1992; Smith and Bretherton, 1972). Furthermore, *δ*>1 is required to
create concentrated flow (Smith and Bretherton, 1972). The change in landscape
elevation is then given by

which can be written as

Equation (6) is non-linear in the case that *γ*≠1.
In deriving this equation of elevation change due to sediment transport we
have simply summed the two terms for sediment flux, the linear and
potentially non-linear slope-dependent terms. This summation has been done as
it is the simplest way to generate landscape profiles that have the desired
convex and concave elements observed in natural landscapes
(Smith and Bretherton, 1972).

To solve this equation in one dimension we assume that the water flux is a
function of the precipitation transported down the river network. The water
collected is taken from the upstream drainage area, *a*, which is related to
the main stream length, *l*, by *l*∝*a*^{h}, where *h* is the exponent
taken from the empirical Hack's law (Hack, 1957). The main stream length
is related to the longitudinal length of the catchment by *l*∝*x*^{d},
where $\mathrm{1}\le d\le \mathrm{1.1}$ (Maritan et al., 1996; Tarboton et al., 1990).
Therefore, we can write $x\propto {a}^{h/d}$, and the water flux is the
precipitation rate, *α* (m yr^{−1}), multiplied by the length
of the drainage system,

where *k*_{w} is a constant of proportionality with units m^{1−p}
and $p=d/h$. Furthermore it is observed that river catchments are typically
longer than they are wide, and so *p*<2 (Dodds and Rothman, 2000). Therefore given
that $\mathrm{0.5}<h<\mathrm{0.7}$ (Rigon et al., 1996), then $\mathrm{1.4}<p<\mathrm{2}$, and
the transport model (Eq. 6) becomes

For simplicity, we will also assume *k*_{w}=1 and vary *c* when
exploring the model behaviour.

However, returning to Eq. (3), it is clear that the transport model is not the only solution. If we assume that the rate of change in sediment thickness is zero over geological timescales, which is to say all sediment created is transported out of the model domain, then Eq. (3) becomes

This assumption has been made previously when studying small mountain catchments where the river may erode directly into the bedrock. However, recent numerical studies, such as Rudge et al. (2015), have expanded this model to cover continent-scale landscapes.

It is clearly plausible to suppose that erosion is primarily due to flowing
water, so the assumption of geologically instantaneous transport may well be
valid for mass that is transported as suspended load within the water column.
Such an assumption is less clear for bedload transport. We can assume that
the speed at which suspended loads travel down-system is a function of the
height achieved within each hop, which is a function of the water depth,
settling velocity, and flow velocity. For small grains <1 mm, the settling
velocity is given by the force balance between the weight of the grain and
the viscous drag given by the Stokes law (Dietrich, 1982). For a particle of
diameter $\mathrm{1}\times {\mathrm{10}}^{-\mathrm{4}}$ m and density 2800 kg m^{−3} the settling
velocity is ∼0.01 m s^{−1}. Therefore the distance travelled
assuming a flow velocity of 1 km h^{−1} and an elevation of suspension of
1 m is roughly 3 km. Using a similar argument the travel distance of a
sediment grain typical of the Bengal Fan is estimated to be ∼10^{4} m
(Ganti et al., 2014). This suggests that rapid transport of sediment
across a continent is possible.

The percentage of mass transported in suspension may also be quite significant. For a small Alpine braided river it was found that the majority of mass was transported as suspended load (Meunier et al., 2006), and for the river systems draining the Tian Shan, China, 70 % of mass is transported as suspended and dissolved load (Liu et al., 2011). Therefore significant mass may be transported rapidly, geologically instantaneously, down-system, suggesting that the assumption that ${\partial}_{\mathrm{t}}h\sim \mathrm{0}$ may be valid in some circumstances.

Assuming surface flow is the primary driver of landscape erosion and that
positive *x* is in the downstream direction, then erosion, *E*, as a function
of the power of the flow to detach particles of rock per unit width can be
written as

where *k*_{b} is a dimensional constant that parameterizes bedrock
erodibility (Howard and Kerby, 1983; units
(m^{2} yr^{−1})^{1−m} yr kg^{−1}), *ρ*_{w} is water
density, *g* is gravity, and *m* and *n* are constants. The exponent *m*∼0.5, as it is a function of how the stream flow width is proportional to the
water flux (Lacey, 1930; Leopold and Maddock, 1953; Whittaker et al., 2007). The
exponent *n*>0 acts upon the slope.

In two dimensions the change in elevation is then given by

where the constant *k* lumps together the other constants (units
m^{−1} (m^{2} yr^{−1})^{1−m}), and if *n*≠1
Eq. (11) becomes non-linear. Using a version of
Eq. (11) to invert river profiles for uplift histories, it
is argued by some authors that *n* is close to unity (Rudge et al., 2015).
However, certain river profiles may arguably be indicative of *n*>1
(Lague, 2014). Furthermore if *n*>1, Eq. (11)
becomes non-linear and the model response to precipitation rate change will
become a function of both uplift and precipitation rates
(Whipple, 2001).

To solve Eq. (11) in 1-D, as before we will assume that
*q*_{w}=*k*_{w}*α**x*^{p} where $\mathrm{1.4}<p<\mathrm{2}$. The stream
power law for landscape erosion in 1-D is then

where *k*_{p}=*k**k*_{w}*ρ*_{w}*g* (units
m^{−p} (m^{2} yr^{−1})^{1−m}).

We have demonstrated two different fundamental equations for change in elevation in 2-D (Eqs. 6 and 11) and the equivalent 1-D forms (Eqs. 8 and 12). These two models of elevation change differ in that Eq. (11) is an advection equation and Eq. (6) is a diffusion equation. This means that the time evolution of Eq. (11) would be a migrating wave of erosion travelling either up or down the catchment (Braun et al., 2015). This wave could also potentially take the form of a shock wave, in which due to the change in gradient the lower reaches of the migrating wave could travel faster than the upper reaches, creating a breaking wave (Pritchard et al., 2009; Smith et al., 2000). The time evolution of Eq. (6) is very different because here the evolution is dominated by diffusive processes. The diffusion coefficient is a function of the down-system collection of water, which can lead to the concentration of flow and the creation of realistic morphologies (Smith and Bretherton, 1972). It is not, however, completely established how the transport model responds differently to changes in precipitation forcing in comparison to the stream power model.

## 2.2 Linear and non-linear solutions

If *n*=1 (Eqs. 11 and 12) and *γ*=1 (Eqs. 6 and 8) then the models are
linear, and we can solve the equations both analytically and in 1-D and 2-D
numerical schemes. For the stream power model we use an implicit finite-difference scheme (Braun and Willett, 2013) and for the transport model we use an
explicit finite-element scheme with linear elements (Simpson and Schlunegger, 2003). If
*n*≠1 and if *γ*≠1 the equations become non-linear. In this case
the numerical solutions can become unstable for simple explicit schemes and
may suffer from too much numerical diffusion for implicit schemes, unless the
size of the time step is limited by the appropriate Courant–Friedrichs–Lewy
(CFL) condition (Campforts and Covers, 2015). Given the short time steps required
to obtain an accurate solution, we explore the non-linear solutions for
erosion down a river-long profile in 1-D. We solve for the stream power model
(Eq. 12) using an explicit total variation diminishing
scheme with the appropriate CFL condition (Campforts and Covers, 2015). For the
transport model (Eq. 8) we use an explicit finite-element
model with quadratic elements and the appropriate CFL condition to find a
stable solution.

## 2.3 Generalizing to a two-dimensional system

To solve Eqs. (6) and (11) over a 2-D domain requires an algorithm to route surface flow down the landscape. In our case, to explore how a model landscape responds to a change in precipitation rate we will make the simplest assumption available: that water flows down the steepest slope. We then solve for Eq. (11) using the numerical model Fastscape (Braun and Willett, 2013), with a resolution of 1000 by 1000 nodes for a 100 by 100 km domain, giving a spatial resolution of 100 m. Erosion by sediment transport in 2-D is solved following the MATLAB model of Simpson and Schlunegger (2003), which is available from Simpson (2017). We solve Eq. (6) on a triangular grid with a resolution of 316 by 316 nodes for a 100 by 100 km domain, giving a spatial resolution of the order of 300 m. We also explored how the models evolve for a domain that is 500 by 500 km in size. The time step used for both models is 10 kyr.

We will explore how an idealized landscape evolves under uniform uplift at a
rate of 0.1 mm yr^{−1}. The initial condition is of a flat surface with a
small amount of noise added to create a roughness. The boundary conditions
are fixed elevation at the left and right sides and no flow at the
sides. To explore the response of the two models to a change in precipitation
rate we start the model with an initial precipitation rate of *α*_{0}=1 m yr^{−1}. For the linear models we then increase or decrease the
precipitation rate to a new value, *α*_{1}, after 10 Myr of model run
time. This is to be sure that the steady state has been reached before
applying the perturbation. For the non-linear models (Sects. 3.3 and 3.4),
the precipitation rate is changed after 5 Myr as in this case steady state
was reached earlier. As the coefficients *c* and *k* have units that are
related to the exponents *δ* and *m* in Eqs. (6) and
(11), respectively
(Armitage et al., 2013; Whipple and Meade, 2006), when modelling
increasing values of *δ* and *m* the coefficients are likewise
increased.

The response time for the transport model scales by the effective diffusivity and can be given by

where *L* is the model length scale (in this case the length of the domain).
For the stream power model the response time is a function of the velocity at
which the kinematic wave travels up the catchment
(Whipple, 2001; Whipple and Tucker, 1999). The response time is therefore
given by the time it takes for this wave of incision to travel up the
catchment length, *l*_{c},

Therefore we expect the response time to be a function of the choice of both
the constants *c* and *k* and the exponents *δ* and *m* within both
models. The effect of varying the coefficients *m* and *δ* independently
has been previously explored
(Armitage et al., 2013; Whipple and Meade, 2006), and we therefore will
not do so in detail again here. Instead we aim to compare the two models and
therefore search for the values of *c*, *k*, *m*, and *δ* that generate
similar topography at steady state. This steady state is then perturbed by a
change in precipitation rate.

## 2.4 Generating similar landscapes

It has been previously demonstrated that both end-member models can generate
convex-up long profiles
(Crosby et al., 2007; Kirkby, 1971; Smith and Bretherton, 1972; Smith et al., 2000; Whipple and Tucker, 2002).
From solving both Eqs. (8) and (12),
where *γ*=1 and *n*=1, we find that in the range $\mathrm{1}<\mathit{\delta}\le \mathrm{1.5}$ and
$\mathrm{0.3}\le m\le \mathrm{0.7}$ the two end-member models are comparable (see Appendix
A). Given the possible additional degree of freedom introduced if we also
vary *γ* and *n*, it is clear that river-long profiles are not a unique
identifier of erosional processes. However, in order to compare how the
end-member models respond to a change in precipitation rate, it is preferable
to perturb catchments of a similar morphology. We will subsequently explore
how the models in their linear and non-linear forms respond to a change in
precipitation rates within the Results section.

### 2.4.1 Erosion by sediment transport

Six models have been run without a change in precipitation to find the steady-state topography. The models explored are first a set of three with varying
*δ* and constant *c*, i.e. *δ*=1.1, 1.3, and 1.5 with
$c={\mathrm{10}}^{-\mathrm{4}}$ (m^{2} yr${}^{-\mathrm{1}}{)}^{\mathrm{1}-\mathit{\delta}}$ (Fig. 2a),
and a set of three in which *δ* and *c* co-vary, i.e. *δ*=1.1 with
$c={\mathrm{10}}^{-\mathrm{2}}$ (m^{2} yr${}^{-\mathrm{1}}{)}^{\mathrm{1}-\mathit{\delta}}$, *δ*=1.3 with
$c={\mathrm{10}}^{-\mathrm{3}}$ (m^{2} yr${}^{-\mathrm{1}}{)}^{\mathrm{1}-\mathit{\delta}}$, and *δ*=1.5 with
$c={\mathrm{10}}^{-\mathrm{4}}$ (m^{2} yr${}^{-\mathrm{1}}{)}^{\mathrm{1}-\mathit{\delta}}$ (Fig. 2b).
These values are chosen because they generate response times within the range
of observations from normal fault-bounded sedimentary systems that have
responded to changes in slip rate
(Armitage et al., 2011; Densmore et al., 2007).

When the transport coefficient *c* is the same for the three values of the
exponent *δ*, the model wind-up time increases with decreasing *δ*
and takes several million years when *δ*<1.5
(Fig. 2a). Steady-state sediment flux is greater for
increasing *δ* when *c* is kept constant. The dimensions (units) of *c*
depend on *δ*, which means that the value of the coefficient *c* must be
adjusted when *δ* is changed to yield the same unit erosion rate per
water flux regardless of *δ* (Armitage et al., 2013).
Consequently, when *c* is suitably adjusted the model can reach a steady
state in a similar time for all three values of *δ*
(Fig. 2b).

We subsequently analyse the topography for the relationship between trunk
river slope and drainage area, as shown in Fig. 3, using TopoToolbox2
(Schwanghart and Scherler, 2014). For the case in which *δ*=1.5 the scaling between
channel slopes and catchment drainage areas, the slope–area exponent
*θ*, is equal to −0.42, and for *δ*=1.3, *θ* is equal to
−0.23 (Fig. 3b). The same value is calculated using the
spatial transformation described in Perron and Royden (2012), commonly referred
to as *χ* analysis (Table 1). Given the reduction in
*θ* from *δ*=1.5 to 1.3, we did not analyse the case for
*δ*=1.1 as the slope–area relationship will clearly lie below the
observed range ($-\mathrm{0.7}<\mathit{\theta}<-\mathrm{0.35}$;
e.g. Snyder et al., 2000; Wobus et al., 2006). Therefore, for river
networks defined by routing water down the steepest slope of descent, the
transport model can create catchment morphologies that have a concavity
similar to that observed in nature if *δ*∼1.5.

### 2.4.2 Comparison to erosion by stream power

In order to provide a comparison for the morphology of the transport model we
explore how the stream power model evolves to a steady state. The landscape
derived from the stream power model, as shown in Eq. (11), evolves
towards a steady state with a slightly different behaviour in comparison to
the transport model (Fig. 4). As before we ran six models
for which in this case the first set of three are *m*=0.3, 0.5, and 0.7 with
$k={\mathrm{10}}^{-\mathrm{5}}$ m^{−1} (m^{2} yr${}^{-\mathrm{1}}{)}^{\mathrm{1}-m}$
(Fig. 4a). The second set of three are of *m*=0.3 with
$k={\mathrm{10}}^{-\mathrm{4}}$ m^{−1} (m^{2} yr${}^{-\mathrm{1}}{)}^{\mathrm{1}-m}$, *m*=0.5 with
$k={\mathrm{10}}^{-\mathrm{5}}$ m^{−1} (m^{2} yr${}^{-\mathrm{1}}{)}^{\mathrm{1}-m}$, and *m*=0.7 with
$k={\mathrm{10}}^{-\mathrm{6}}$ m^{−1} (m^{2} yr${}^{-\mathrm{1}}{)}^{\mathrm{1}-m}$
(Fig. 4b). This range of *m* is chosen as it spans the
range of observed concavities within catchments. As with the transport model
the coefficient *k* can be adjusted along with *m* as they are related, and
increasing *k* reduces the model wind-up time (Fig. 4b).
Decreasing the exponent *m* increases the timescale taken to reach a steady
state (Fig. 4a); however, by varying *k* by a factor of
100, the steady-state sediment flux is reached within 3 Myr for the three
values of *m* (Fig. 4b).

Following the previous examples, we analyse the topography for the
relationship between trunk river slope and drainage area
(Fig. 5). Both the transport model and the stream power
model can create landscapes with similar slope–area relationships calculated
using the *χ*-analysis approach (Table 1). For both
models, the values of the intercept, *k*_{s}, and the gradient,
*θ*, are of similar magnitudes for *δ*=1.5 and *m*=0.5.
Absolute elevation for the model shown in Fig. 5a is
higher than the transport model example due to the larger value of *k*
relative to *c*. However, importantly, both models can create similar
landscape morphologies at steady state.

The stream power and transport models can both fit the observed slope–area
relationships of the present day landscape morphology, for example *θ* ranging
from −0.35 to −0.7 (Snyder et al., 2000; Wobus et al., 2006), when the
water flux exponent is *m*∼0.5 or *δ*∼1.5 for the stream power
and transport model, respectively. Therefore, both models may be a reasonable
representation of how, on a gross scale, a landscape erodes. We therefore
keep the exponents in the range $\mathrm{0.3}\le m\le \mathrm{0.7}$ and $\mathrm{1.3}\le \mathit{\delta}\le \mathrm{1.5}$ and explore how the models in their linear and non-linear forms
respond to a change in precipitation rates.

## 3.1 Response to precipitation rate reduction

When the model is perturbed by a change in precipitation rate the sediment
flux output will first change as the erosive power changes (e.g.
Fig. 6). The model will subsequently return to the steady-state
output, as the slope of the fluvial system will adjust to the new
precipitation rate, and the landscape will re-achieve the same steady state.
In Fig. 6a we display the response of erosion for the transport
model in terms of sediment flux out of the model domain for a reduction in
precipitation rate from 1 to 0.5 m yr^{−1} at 10 Myr of model evolution.
We explore how the transport model responds for *δ*=1.5, $c={\mathrm{10}}^{-\mathrm{4}}$ and
*δ*=1.3, $c={\mathrm{10}}^{-\mathrm{3}}$ as these two values of *δ* generate reasonable
slope–area relationships (Fig. 3b,
Table 1). The response to a reduction in precipitation is
similar for the two model parameter sets, with the flux initially reducing by
half and then recovering to within 10 % of steady-state values within
∼ 2 Myr (Fig. 6a; see Table 2). Changing
the transport coefficient, *c*, does not affect the predicted gradient of
catchment slope versus catchment area (see Appendix A,
Fig. A2). However, changing *c* changes the model
elevation (Fig. A2). Furthermore, the larger the value
of *c* the faster the response (Eq. 13; see
Armitage et al., 2013). A small increase in the exponent *δ*
will strongly reduce response times, as it will increase the water flux term
(Eq. 13). Therefore an order-of-magnitude decrease in *c* counters
the change in *δ* for the two model sets (Fig. 6a). For the
values chosen both models respond at a similar rate to the change in
precipitation (Fig. 6a; see Table 2).

The response of the stream power model to an identical reduction in
precipitation at a model time of 10 Myr takes a similar form, with an
initial decrease in sediment flux out followed by a gradual recovery
(Fig. 6b). In a similar manner as the transport model, response
is a function of the exponent *m* and the coefficient *k*
(Eq. 14). We have modelled three parameter sets: *m*=0.3 and
$k={\mathrm{10}}^{-\mathrm{4}}$, *m*=0.5 and $k={\mathrm{10}}^{-\mathrm{5}}$, and *m*=0.7 and $k={\mathrm{10}}^{-\mathrm{6}}$
(Fig. 6b). The response time to achieve a return to 10 % of the
steady-state sediment flux varies from 3 Myr in the case of *m*=0.3 to less
than 1 Myr when *m*=0.7. In addition to response time being longer for smaller
values of *m*, the peak magnitude of the flux response is reduced for smaller
values of *m* (Fig. 6b).

The magnitude of the response for all the runs is greater for the transport
model when compared to the stream power model (Fig. 6).
Consequently, response time, while being a function of the transport
coefficients *c* and *k*, may still systematically differ
between the two models: the transport model with *δ*=1.5 and $c={\mathrm{10}}^{-\mathrm{4}}$
generates a maximum model elevation of ∼240 m, and the stream power
model with *m*=0.5 and $k={\mathrm{10}}^{-\mathrm{5}}$ generates a maximum elevation of ∼180 m. These two models have a similar slope–area relationship at steady
state (Table 1) and are therefore comparable, suggesting a
faster response to a reduction in precipitation rates for the stream power
model (Fig. 6).

To explore how the difference in response time and magnitude is expressed in
the landscape, we extract the river profiles of the main trunk systems for
models in which *δ*=1.5 and *m*=0.5 during the response to the reduction in
precipitation rate, while the uplift rate is constant
(Figs. 7 and 8). For the
transport model in which *δ*=1.5 and $c={\mathrm{10}}^{-\mathrm{4}}$, the catchment
elevation increases to a new steady state that has an elevation that is
roughly 2.6 times higher than the steady-state elevation after 10 Myr
(Fig. 7). Just under half of this new topographic
elevation is achieved within the first 500 kyr. In contrast, for the stream
power model in which *m*=0.5 and $k={\mathrm{10}}^{-\mathrm{5}}$, the steady-state topography is
achieved within a fraction of the time when compared to the transport model.
This is in line with the more rapid response of this model to a relative
drying of the climate using these parameters (compare Fig. 6a
and b). Furthermore the increase in elevation due to the reduced surface
water flux is only a factor of ∼1.2, which is less than half of the
increase for the transport model. Our results confirm that two different
end-member erosion models encompassing advective and diffusive phenomena
can produce landscapes with similar morphologies if particular parameter
sets are selected accordingly.

## 3.2 Response to different magnitudes of precipitation rate change

The response time of the transport model is known to be a function of the
transport coefficient and the magnitude of the precipitation rate
(Armitage et al., 2013). This behaviour is displayed in
Fig. 9a, in which the response of the transport model with
*δ*=1.5 and $c={\mathrm{10}}^{-\mathrm{4}}$ for a change in precipitation from 1 to 0.25,
0.5, 0.75, and 2 m yr^{−1} is plotted. The response time, measured as the
time for the sediment flux to recover by half and by 90 % to the steady-state value, is shown additionally in Fig. 10 as black solid
and dashed lines, respectively, and in Table 2. For a
reduction to 0.25 m yr^{−1} the prediction is for a long response time of
6.07 Myr, while for an increase to 2 m yr^{−1} the prediction is a for
rapid response time of 310 kyr for 90 % recovery towards previous sediment
flux values. The equivalent half-life, recovery by 50 % towards previous
sediment flux values, is 1.42 Myr and 90 kyr.

The stream power model likewise has a response time that is a function of
precipitation rate (Fig. 9b). For a reduction to
0.25 m yr^{−1} the prediction is for a response time of 1.66 Myr, while
for an increase to 2 m yr^{−1} the prediction is for a recovery time of
600 kyr for 90 % recovery (Table 2). The equivalent half-life is 0.98 Myr and 340 kyr (Table 2). The stream power
model is therefore faster to recover for a reduction in precipitation rate
yet slower to respond to an increase in precipitation rate. This is because
the response time of the stream power model is more weakly a function of
precipitation rate. Importantly, these results therefore suggest that there is a
fundamental asymmetry in the response timescale to a climate perturbation.
The models suggest that it takes longer for surface processes to recover from
a drying event compared to a wetting event.

Both models display a response time that is a function of the precipitation rate (Figs. 9 and 10). The relationship between precipitation rate and the transport model response can be expressed as

where in this case *δ*=1.5. This proportionality is in agreement with
our numerical model results, in which the slope of trend for the transport model
in the log–log plot is −1.5 (Fig. 10).

In contrast, the response time of the stream power model is not as strongly
inversely dependent on the precipitation rate (Fig. 10). For
this model, the response time is a function of the velocity at which the wave
of incision travels upstream. This velocity is directly related to the
inverse of the water flux, ${q}_{\mathrm{w}}^{m}$, which is in turn a
function of the drainage length and precipitation rate, *α*. Therefore
for the stream power model we can write that response time is

This proportionality, which is in agreement with the approximate analytical solutions of Whipple (2001), is likewise in agreement with our numerical model results, in which the slope of trend for the stream power model in the log–log plot is −0.5 (Fig. 10). Consequently, for these two models, which were derived from the same starting point (Fig. 1) and applied to catchments of similar topography and morphology, we find that above a certain magnitude of precipitation rate change, the transport model responds more rapidly than the stream power model and vice versa.

The position of the critical point at which the stream power model responds more
rapidly than the transport model is a function of the water flux and the
collection of coefficients. In the model comparison developed here, we have
compared two model catchments that have a similar slope–area exponent, *θ*
between −0.4 and −0.5 (*δ*=1.5 and *m*=0.5) and model domain length
of *L*=100 km, giving catchments of roughly 50 km length. In this case the
90 % recovery of the sediment flux signal is predicted to be more rapid for
the transport model when compared to the stream power model for an increase
in precipitation rate (Fig. 10). If, however, the model domain
is increased to *L*=500 km then it takes twice as long for the transport
model to recover from an increase in precipitation rate from 1 to
2 m yr^{−1}: 0.63 Myr compared to 0.31 Myr for *L*=100 km
(Fig. 11a and Table 2).

The stream power model is insensitive to the size of the model domain because
of the particular choice of *m*=0.5 and the shape of drainage network that
forms under the assumptions of routing water down the steepest slope of
descent (Fig. 11b). Taking the drainage length to be directly
proportional to the catchment area, *l*_{d}∝*a*, and given that
catchment length is proportional to drainage area raised to the Hack
exponent, *h*, we can re-write Eq. (14) as

Therefore, in the case that *h*=0.5 and *m*=0.5, as in the numerical
model here, the response time becomes independent of system length
(Whittaker and Boulton, 2012). If *h*<*m* then response times would decrease with
increasing drainage basin size, and if *h*>*m* then response times would
increase with drainage basin size. There is good empirical evidence for
$\mathrm{0.5}<h<\mathrm{0.7}$ (Rigon et al., 1996), which fundamentally controls the
plan view shape of catchments, yet there is not a complete consensus on the
value of *m* (Lague, 2014; Temme et al., 2017).

A final key difference between the transient sediment flux responses of the
two models is that the peak magnitude of system response to a change in
precipitation rate is systematically larger for the transport model
(Fig. 9). For an increase in precipitation rates from 1 to
2 m yr^{−1}, the sediment flux increases from 1×10^{6} m^{3} to
2.5×10^{6} m^{3} for erosion by sediment transport. This is
3 times greater than the equivalent increase for the stream power model. The
reduction in sediment flux is likewise larger for the transport model
(Fig. 9). Therefore, although response time is a function of
precipitation rate, the magnitude of change is consistently larger for the
transport model.

## 3.3 Non-linear response timescales

Up to this point we have compared how the models respond to a precipitation
rate change when the solutions are linear. However, there is reasonable
debate as to the value of the slope exponent *n* in the stream power model
(Croissant and Braun, 2014; Lague, 2014; Rudge et al., 2015) and likewise within
the transport model it is plausible that the slope exponent *γ*>1. The
response time for the stream power model for various values of *n* has been
explored within Baldwin et al. (2003). Here we expand on this by exploring
the equivalent response times for the transport model. To explore the
implications of the non-linearity introduced by relaxing the constraint that
*n*=1 and *γ*=1 for both models, we solve Eqs. (8)
and (12) for *p*=1.1, *δ*=1.5, and $c=\mathrm{5}\times {\mathrm{10}}^{-\mathrm{5}}$
and *m*=0.5 and $k={\mathrm{10}}^{-\mathrm{4}}$, respectively, with different uplift rates. We have
modelled the response due to an uplift rate between 0.1 and
1.0 mm yr^{−1} for the case in which *γ*=1.2 and *n*=1.2 in
Eqs. (8) and (12)
(Fig. 12).

We find that for both the transport and stream power model, when the slope
exponent is greater than one, the model response time is a function of uplift
rate. The faster the rate of uplift, the faster the system responds to a
change in precipitation rate. If the response time for a system recovery to
steady state by 50 or 10 % is plotted on a log–log plot against uplift rate
we find that the response time is proportional to the uplift rate raised to a
negative power (Fig. 13). In the case of *n*=1.2 or
*γ*=1.2 the slope of trend is −0.1667, and for *n*=2 or *γ*=2 the
slope of trend is −0.5 (Fig. 13). These slopes are in
agreement with the approximate analytical solutions of Whipple (2001)
and numerical models of Baldwin et al. (2003); i.e. the stream power
response time *τ*_{sp} has a proportionality,

and equivalently we infer from our numerical model (Fig. 13) the transport model response time as

This implies that both models have the same form of response dependency on uplift rates. Therefore, regardless of the rate of uplift we should expect the transport model to respond more rapidly to a large increase in precipitation rate and the stream power model to respond more rapidly to a reduction in precipitation rate (Fig. 10). Our results are also consistent with the field-based findings of Whittaker and Boulton (2012), who showed that landscape response times for rivers close to the detachment-limited end-member were shorter for terrain uplifted by faster-slipping active faults.

## 3.4 Response time as a function of the initial precipitation rate

Up until this point we have only explored how the numerical models respond to
an increase or decrease in precipitation rate by keeping the initial
precipitation rate fixed at *α*_{0}=1 m yr^{−1} and varying the
final precipitation rate *α*_{1}. In this final section we will
instead keep the final precipitation rate fixed at *α*_{1}=1 m yr^{−1} and vary the initial precipitation rate *α*_{0} from
values of 0.5 to 1.5 m yr^{−1}. We will focus again on the 1-D models and
look at the linear and non-linear cases with *n*=1.2 and *γ*=1.2.

For the linear and non-linear transport model we find that if the initial
precipitation is less than the final precipitation (*α*_{0}<*α*_{1}) then the response time is not very sensitive to the initial
precipitation rate (Fig. 14a). If *α*_{0}>*α*_{1} then
the response time is a function of the initial precipitation rate, but the
relationship cannot be explained by a simple power law (Fig. 14a).
The change in response time as a function of the initial precipitation rate
is, however, small compared to the change in response time as a function of the
final precipitation rate.

In the case of the linear and non-linear stream power model, the response time has a no dependence on the initial precipitation rate and is only a function of the final precipitation rate (Fig. 14b). With all other parameters being held constant, the initial precipitation rate will set up the topography and hence the slope of the pre-perturbation landscape. Elevations will be lower for higher precipitation rates, and the topographic gradient will be smaller. For the case of the stream power model, the change in erosion rates migrates up the catchment and so the old topography does not impact the response time. For the transport model, however, the remnant topography does have a small effect on the response time, but only if the previous precipitation rate was higher than the new post-perturbation precipitation rate.

In deriving the two end-member models to describe landscape evolution, we showed that if the rate of transport of sediment were assumed to be instantaneous (i.e. all sediment is transported out of the model domain) then the stream power model would be appropriate to describe catchment erosion. However, if it is instead assumed that the rate of sediment transport is not instantaneous, then we arrive at a model in which erosion scales with the rate of change of sediment flux, which itself is dependent on both linear and potentially non-linear slope-dependent terms. These two end-members can produce similar steady-state landscapes, as noted by a number of previous studies (Tucker and Whipple, 2002; Whipple and Tucker, 2002). However, as we demonstrate above, when perturbed by a change in conditions such as rainfall rate, they can produce very different landscape responses, which vary in terms of their style, magnitude, and tempo. We explore the nature and implications of these responses below.

It is also important to stress that the catchment responses and the predicted sediment fluxes out of these two model domains might be variously relevant to different erosional and depositional domains (Lague, 2014; Temme et al., 2017). A model of instantaneous sediment transport might be more relevant for suspended sedimentary loads, for which transport times can be very small, while the transport model might be more appropriate for bedload-dominated systems, even in cases in which bedrock is clearly incised (Paola et al., 1992; Valla et al., 2010). Furthermore, given that these two models have different response times, it is possible that fine-grained deposits might record signals of climate change differently from, for example, the gravel deposits within alluvial valleys. Below, we will therefore first discuss how the two model responses compare in terms of their response time and place our results in the wider context of sediment routing system response to environmental change. Second, we will compare the model results with three records of change in sediment deposition during the Paleocene–Eocene thermal maximum (PETM), a known and well-constrained period of rapid climate change. Finally, we will summarize the key implications from our results.

## 4.1 Response times as a function of model choice

Under certain parameter sets it is relatively straightforward to generate two
landscapes eroded by the transport or stream power model that have similar
elevation, slope, and area metrics (Figs. 3 and
5). To find a path to break the apparent non-uniqueness
of these solutions we have explored the transient sediment flux response out
of the model domain for two end-member solutions to erosion. The first
observation is that both models respond at a first order in a broadly
similar way to a precipitation rate (climate) driver (Figs. 9
and 10). Both models have a response that is an inverse
function of the magnitude of precipitation rate change. Both models have a
response that is related to uplift in an identical manner
(Fig. 13). However, the responses for catchments that are
comparable in slope–area relationship and maximum elevation, but which are
governed by different erosional dynamics defined by *c*, *k*, *m*, and
*δ*, actually display different response times by almost 1 order of
magnitude (Figs. 2 and 4).

We have demonstrated that models limited by their ability to transport
sediment tend to have shorter response times to an increase in rainfall rate
and thus re-achieve pre-perturbation sediment flux values more rapidly
compared to stream-power-dominated systems, particularly when catchment
length scales are small (e.g. <100 km, Fig. 10). These
model observations suggest that the sediment fluxes from small alluvial
catchments, even when captured in downstream depocentres, may be difficult to
tie to changing climate parameters unless depositional chronologies are
exceptionally well constrained (D'Arcy et al., 2017). Conversely,
catchments whose erosional dynamics lie close to the stream power end-member
model may be well placed to record longer-term climate shifts, but may be
buffered to very high-frequency variations in the climate driver
(Armitage et al., 2013; Simpson and Castelltort, 2012). It is important to
stress that the trend in response is asymmetric, by which we mean that both
models show a faster response for a precipitation increase relative to a
precipitation decrease (Fig. 10). This is an important
outcome, which has to date not been widely recognized or investigate in field
scenarios. In particular, it raises the prospect that for
glacial–interglacial cycles characterized by wetter, cooler stadial
periods and dryer, warmer interstadials, the rapid climate recovery from peak
glacial conditions typically seen in *δ*^{18}O records might be mediated
by a longer landscape response time to this change. Conversely, physically
slower boundary condition changes towards wetter conditions may give rise to
faster landscape response times. We suggest that an exploration of these
differences may be a promising avenue of future research.

Given that the response time is a function of the water flux exponent (*m* or
*δ*) and that the water flux exponent for the transport model is
greater than that for the stream power model, there will be a cross-over
point at which the stream power model responds faster than the transport model.
This cross-over point is a function of the erodibility coefficient *k* and
the transport coefficient *c*. In the scenario in which we have tried to
initiate the perturbation in precipitation rates from similar catchments, we
find that this cross-over point is towards large reductions in precipitation
rates (Fig. 10). This implies that the transport model
generally responds faster than the stream power model (10^{5} to
10^{6} yr) for examples in which the parameter combinations used here
produce similar steady-state landscapes.

For such conditions, the stream power model predicts a landscape response
time to a change in precipitation of the order of 10^{6} yr, and this time
is related to the precipitation rate to the inverse power of *m*
(Fig. 10). The transport model predicts a wider range of
response times of the order of 10^{6} to 10^{5} yr that is related to the
precipitation rate to the inverse power of *δ*; in this case the
response time is also length dependent (Fig. 10 and
Table 2). It has been suggested that a transition from a
landscape controlled by detachment-limited erosion (stream power model) to
sediment transport at longer system lengths may explain the longevity of
mountain ranges (Baldwin et al., 2003). This hypothesis is somewhat backed
up by the analysis of response times for the transport model, as the response
time increases with system length (Table 2) unlike the stream
power model, which has a response that is only slightly modified by system
length (Baldwin et al., 2003; Whipple, 2001). To date, physical constraints
on landscape and sediment flux response times to climate changes in the
geologic past are relatively scarce
(Ganti et al., 2014; Romans et al., 2016; Temme et al., 2017) because real
systems are complex. They include internal dynamics, such as vegetation and
autogenic behaviour, which are often omitted from model studies, and because
of the need for stratigraphic archives to be complete with well-established
chronologies (Allen et al., 2013; Forman and Straub, 2017). In principle, however, the
dominant long-term incision process governing catchment behaviour
fundamentally determines the sediment flux response and may itself help
identify catchment erosional dynamics; we explore this question in Sect. 4.2.

Finally, it is worth noting that the model response time has implications for
the inverse modelling of river profiles. When river-long profiles are inverted
for uplift, erosion is typically assumed to be captured by the stream power
model (Pritchard et al., 2009). Studies of continent-scale
inversion have found that the best fit value of *k* for the stream power
model increases by 2 orders of magnitude to fit river profiles in Africa
relative to Australia (Rudge et al., 2015). Such a large change in *k*
would result in a highly significant difference in response time from
continent to continent, which in itself would imply that tectonic and
climatic signals are preserved in landscapes and sediment archives over
vastly different time periods (Demoulin et al., 2017). Such an
outcome may reflect fundamental differences in bedrock erodiblity
(Roy et al., 2015), but alternatively could be satisfactorily explained by
differing long-term erosional dynamics and sediment transport. These
differences are enhanced in the case in which *n*>1 in the stream power
erosion model.

## 4.2 Relevance of model responses to sediment records of climate change

To what extent do these model results, which start from similar steady-state topographies, help us to understand whether stratigraphic records of sediment accumulation through time do or do not reflect the effects of climatic change on sediment routing systems governed by differing long-term erosional dynamics? One motivation for this study has come from the increasing number of field and stratigraphic investigations of terrestrial sedimentary deposits, apparently contemporaneous with (and taken to record) known past climate perturbations, such as the Palaeocene–Eocene thermal maximum (PETM), a hyperthermal event that occurred around 56 Ma. Stratigraphers often correlate changing stratigraphic characteristics with changing environmental boundary conditions in a qualitative way (Allen, 2017; Romans et al., 2016). However, to evaluate quantitatively how sediment routing systems respond to climate with reference to real examples, it is imperative to consider systems in which the timescales of erosion (or as a proxy, deposition) are known, stratigraphic sections are complete, and the driving mechanisms well documented (Allen et al., 2013; D'Arcy et al., 2017).

To compare our model predictions with observations it is clear that we have to use the depositional record. Therefore, there is an implicit assumption that stratigraphy is a faithful recorder of erosion. It is, however, possible that climatic change will also alter processes that control sediment deposition, for example by altering how sediment partitions from transport into stratigraphy. By using estimates of the total volume of sediment deposited within the Escanillia Eocene sedimentary system in the Spanish Pyrenees, it has been demonstrated that climatic change can recreate observed changes in grain size deposition (Armitage et al., 2015). This example of a close model-to-stratigraphic-observation prediction might be evidence that the stratigraphic record is a faithful record of a change in sediment flux delivery to the depositional environment.

The PETM is arguably the most rapid global warming event of the Cenozoic,
with a rise in global surface temperatures by 5 to 9 ^{∘}C
(Dunkley Jones et al., 2010), forming a clear step change in climate for
which depositional records are well constrained in a number of basins
worldwide (Foreman et al., 2012). It is therefore a good example for
high-level comparison with our model outputs. While the large-magnitude
glacial–interglacial cycles of the past 1 Myr are also plausible candidates
to investigate these links in principle
(Armitage et al., 2013), we note that many terrestrial records
of sedimentation over ca. 100 kyr, such as fluvial terraces and alluvial
fans, have depositional chronologies that are often incomplete or reworked
(D'Arcy et al., 2017; Demoulin et al., 2017).

The initial warming associated with the PETM event occurred at ca. 55.5 Ma
and may have been as abrupt as 20 kyr, with a duration of 100 to 200 kyr
based on the synthesis of *δ*^{13}C and *δ*^{18}O records
(Foreman et al., 2012; Schmitz and Pujalte, 2007). The event has been associated
with clear changes to global weather patterns; for instance, hydrogen isotope
records suggest increased moisture delivery towards the poles at the onset of
the PETM, consistent with predictions of storm track migrations during global
warming (Sluijs et al., 2006). This event has also been argued by an
increasing number of authors to have produced a significant geomorphic and
erosional impact based on sedimentary evidence and its apparent effect on
the global hydrological cycle and catchment run-off
(Foreman, 2014; Foreman et al., 2012).

A clear response to the PETM is recorded within both the Spanish Pyrenees and the western US; however, the responses are arguably not the same. At the onset of the PETM there is strong evidence for the contemporaneous increase in precipitation rates and the deposition of coarse gravels known as the Claret Conglomerate (Schmitz and Pujalte, 2007) in the Tremp Basin of the Spanish Pyrenees. In the western US the PETM is marked by the deposition of well-documented channel sandstone bodies in the Piceance Creek and Bighorn basins (Foreman, 2014; Foreman et al., 2012). In the US cases, the deposits include coarse channelized sands, marked by upper flow regime bed forms, some of which are consistent with a synchronous increase in both water and sediment discharge. At Claret, where the style of sedimentation abruptly changes from a semi-arid alluvial plain to an extensive braid plain or megafan deposit, the conglomerate has a thickness of ∼10 m, while the total carbon isotope excursion (CIE) in the same section measures ∼35 m (Manners et al., 2013).

### 4.2.1 Claret Conglomerate, Spanish Pyrenees

The Claret Conglomerate was likely deposited rapidly, representing a fast
response to climate change. If we assume a constant rate of deposition, then
the Claret Conglomerate accounts for roughly 30 % of the total duration of
deposition for the CIE (170 kyr; Röhl et al., 2007), suggesting
that deposition occurred over a duration of up to 50 kyr. Indeed,
Schmitz and Pujalte (2007) argue that the deposition of this unit may have been
markedly quicker than the conservative estimate above, perhaps taking less
than 10 kyr, based on their detailed comparison of *δ*^{13}C and
*δ*^{18}O records at the field site compared to marine records of the
excursion. Therefore unless there is a major unconformity within the CIE, the
implication is that the erosional system responded rapidly at this particular
field site, in 10 to 50 kyr, to a significant shift in climatic conditions.
These values suggest sedimentation rates of up to 1 mm yr^{−1}. If such a
sedimentation rate had been sustained for the duration of the deposition of
the Tremp Group (Maastrichtian to the end of the Palaeocene), the sediment thickness
would be >15 km. This is an order of magnitude more than actually observed
(Cuevas, 1992) and would therefore suggest that sediment fluxes
increased dramatically at the PETM.

Erosional source catchment areas were likely <100 km in length given the
palaeo-geography of the Pyrenees at the time (Manners et al., 2013). The
very short duration of the erosional response, which is required for the
sediments to be transported and deposited on a timescale of ca. 10^{4}
years, is therefore difficult to model within a stream power (advective)
end-member model for catchments of this scale (Table 2),
although a version of such a model has been recently used to explore the
controls on the evolution of later Miocene megafans in the northern Pyrenees
(Mouchené et al., 2017). To use the stream power model would
require a significant increase in the bedrock erodibility parameter, *k*, within
the model (by greater than 1 order of magnitude), implying slopes and
topography in the palaeo-Pyrenees that were highly subdued. In contrast, the
sediment transport model more easily reproduces the documented response
timescales given an increase in precipitation; it is also consistent with the
volumetrically significant export of bedload-transported gravel clasts and
therefore honours the independent field data more effectively. We also note
that the transport model displays a response time that has a stronger
dependence on precipitation rate change and has a greater amplitude of
perturbation (e.g. Fig. 9). We therefore suggest that the erosional
pulse that led to the deposition of the Claret Conglomerate is most
appropriately modelled as a diffusive system response to a sharp increase in
precipitation over the source catchments of the developing Pyrenean mountain
chain at that time.

### 4.2.2 Sandstone bodies in the Piceance Creek and Bighorn basins, western US

The time-equivalent sections in the Bighorn and Piceance Creek basins of the western US also provide clear evidence of anomalous sedimentation at the PETM; however, in this case the duration of deposition is somewhat longer, >100 kyr (Foreman, 2014; Foreman et al., 2012). Here the deposits are of smaller grain sizes, with the boundary sandstone sequence in the Bighorn Basin being made up of fine to coarse sand with little gravel (Foreman, 2014). In the Piceance Creek basin, the PETM section documents the rapid progradation of coarse-grained sands, which is consistent with greater discharges, over silty underlying strata and in that sense these observations also match sediment transport model predictions for rapid increases in sediment flux driven by enhanced precipitation (Foreman et al., 2012). However, it is notable that the documented changes in fluvial style persisted beyond the PETM and we therefore suggest that the fast response of the system to the increase in precipitation, but the persistence of coarser-grain sedimentation as the climate presumably dried and cooled, may indeed reflect the marked asymmetry in sediment flux responses to wetting and drying noted in Fig. 9.

In contrast, the Bighorn Basin boundary sandstone sediments are contained within the PETM time period and indicate uniform flow depths and widths during this time, while also being coarser than the underlying horizons. Moreover, proxy data suggest a net decrease rather than an increase in precipitation (Foreman, 2014; Foreman et al., 2012). While the progradation of such coarse-grained facies could be represented as a diffusive process driven by increasing rainfall-driven discharge (Armitage et al., 2011; Paola et al., 1992), this is apparently inconsistent with the sedimentological characteristics of the deposit. Although sediment fluxes are not explicitly reconstructed in this work, this response apparently requires greater volumes and grain sizes of sediment delivered despite lowered rainfall conditions and is thus difficult to capture in either of the end-member models used here. Foreman et al. (2012) and Foreman (2014) argue for the preferential removal of fine-grained floodplain deposits speculatively linked to changing vegetation and the reduced cohesion of overbank sediments. Consequently, while two of the PETM sections considered here are broadly consistent with landscape responses governed by a sediment transport model, some depositional stratigraphies are not immediately consistent with either end-member model and may reflect important complexity, such as the effects of vegetation, lacking from simple model solutions.

## 4.3 Summary and model limitations

In this section we consider the implications of our model
outputs, both generally for interpreting sediment routing system response to
boundary condition change and specifically in the context of the well-studied
PETM event. While the sediment flux response of the models to a change in
precipitation are at a first-order level broadly similar, there are four key
differences to highlight. First, starting from the same initial conditions,
the sediment transport model appears to be more sensitive to precipitation
change than the equivalent stream power model. It is therefore a good
candidate for which rapid catchment-wide responses are recorded to, for example, a climate
change event, as we argued for the PETM Claret Conglomerate in the Spanish
Pyrenees. Second, we note that in both model cases there is a quicker response
to a wetting than a drying event, something which has not been well established
or demonstrated from field observations. Nevertheless we argue that field
data sets, including PETM studies, may already have recorded this asymmetry,
although it may not have been recognized as such. Third, the sediment
transport model has a greater magnitude of peak sediment flux and is
particularly sensitive to catchment size. Finally, we note that response time
in both models is a function of uplift rate for *n*>1, which means that in such
cases, perhaps counter-intuitively, the more perturbed the system the
faster it responds (Whittaker and Boulton, 2012).

However, it is important to recognize that in deriving these two classic
end-member models we have simplified landscape evolution considerably. We
acknowledge that change in the model parameters, *c*, *k*, *m*, and *δ*,
will alter the response times depicted here
(Armitage et al., 2013). However, in order to compare the two
models we have specifically used values of *c*, *k*, *m*, and *δ* that
generate comparable model landscapes, and we then changed the precipitation
rate to understand the form of the model response. No model incorporates all
the complexities that characterize sediment routing systems from source to
sink (Allen, 2017) and the act of simplification inherent in
considering erosional end-member models necessitates that in arguing for the
applicability of one over the other, we simply consider the broad styles of
behaviour suggested by model outputs. We do not, for example, consider
autogenic behaviours (Forman and Straub, 2017), nor do we consider coupled
issues of vegetation turnover in response to climate change, which may a
play an important role in examples such as the Bighorn Basin considered here
(Foreman, 2014). Nonetheless, a significant finding of this work
has been the clear asymmetry in response time of these end-member models in
terms of a wetting event (faster) compared to a drying event (slower). This
implies that aridification events are harder to preserve in the sedimentary
record, not only because they are typically associated with reduced sediment
fluxes, but also because the timescale of landscape response may be >10^{6} years.

Deterministic numerical models of landscape evolution rest on fundamental assumptions on how sediment is transported down-system. The stream power law is based on the assumption that all sediment generated is transported instantaneously out of the landscape. Transport models assume that there is an endless supply of sediment to be transported. The existence of knickpoints within river-long profiles, assumed to be produced by a system perturbation such as a base level, has been used to provide evidence in support of the stream power law in upland areas (Snyder et al., 2000; Whipple and Tucker, 1999; Whittaker et al., 2008). Knickpoints, however, can likewise be a result of changes in lithology (Grimaud et al., 2014; Roy et al., 2015) and are certainly not a unique indicator of erosion dynamics (Grimaud et al., 2016; Tucker and Whipple, 2002; Valla et al., 2010). In this contribution we therefore attempted to understand how the sediment flux signal out of the eroding catchment may generate a distinguishable difference between the end-member models in terms of a response to a change in run-off. This idea is motivated from field observations of past landscape responses to climate excursions, such as the PETM, which are manifested in the rapid deposition of coarse sedimentary packages in terrestrial depocentres (Armitage et al., 2011; Foreman et al., 2012).

Both models suggest that the response time of landscape to a change in precipitation rate has a proportionality of the form of a negative power law (Eqs. 15 and 16). The key difference is in the value of the exponent. For the stream power model, the exponent must be less than one in order to match the observed concavity of river profiles. In contrast, for the transport model the exponent on the precipitation rate must be greater than one in order to generate a river network (Smith and Bretherton, 1972) and to generate the observed concavity of river profiles. This results in the transport model responding more rapidly to an increase in precipitation rate in comparison to the stream power law model (Fig. 10). In contrast, the stream power model is faster to respond to a reduction in rainfall rate. This is fundamentally because the response time of this model is more weakly a function of precipitation than the sediment transport model. Significantly, therefore, our results show that there is a fundamental asymmetry in the response of both models to a climatic perturbation, with the response time to a drying event longer than that to an increase in rainfall. In general terms, the magnitude of the response to a change in precipitation rate appears greater across the range of model space investigated here for the sediment transport (diffusive) model solutions, while for the stream power (advective) model, the magnitude of the sediment flux perturbation is smaller, but is more localized within the catchment with respect to knickpoint retreat.

While this study does not address whether or not these sediment flux signals
will be preserved in the stratigraphic record, a problem that fundamentally
rests on the availability of accommodation to capture the eroded sediment
(Allen et al., 2013; Whittaker et al., 2011), it does suggest that
landscapes governed by these simple erosional end-members should be sensitive
to climate change. Moreover, there are some important diagnostic
differences between their sediment flux responses to an identical
perturbation, including the amplitude, timescale, and locus of the erosional
response. Using published stratigraphic examples, we suggest that the
timescales and magnitude of coarse sediment deposition in the Spanish
Pyrenees at the time of the PETM are best described using the diffusive
transport model end-member. Moreover, we argue that these model end-members
allow us to constrain the range of likely sediment flux scenarios that
precipitation changes may generate and that numerical models, in conjunction
with a range of field and independently constrained proxy data sets, are best
placed to tease apart when and in what circumstances climate signals are
likely to have been generated in erosional catchment systems, which
fundamentally determines whether they can be subsequently *captured*
in sedimentary depocentres downstream.

The 1-D solution to the transport model is available from John Armitage (armitage@ipgp.fr). The 1-D solution to the stream power model is available from Benjamin Campforts (benjamin.compforts@kuleuven.be). Fastscape is available from Jean Braun (GFZ Potsdam) by request. The 2-D solution to the transport model was developed by Guy Simpson (University of Geneva) and is available as part of Simpson (2017).

The solution to the one-dimensional stream power law (Eq. 12) assuming that at the end of the catchment at *x*=*L*
elevation *z*=0 and *m**p*≠1 is

and for the case in which *m**p*=1 this simplifies to

For the transport model (Eq. 8) there is an exact solution
for the case that *δ**p*=2, which assumes that at *x*=0, ${\partial}_{x}z=\mathrm{0}$ and at *x*=*L*, *z*=0:

where

For other values of *δ**p* the steady-state solution is solved
numerically; Eq. (8) is solved using the finite-element method with linear weight functions. We use a non-uniform 1-D nodal
spacing, for which the spatial resolution is increased with increasing gradient.
The numerical model is bench-marked against the analytical solution for the
case in which *n**p*=2.

The steady-state solutions are plotted in the case that $\mathit{\delta}=p=\sqrt{\mathrm{2}}$ and for reference the stream power model solution with *m*=0.5
and $p=\sqrt{\mathrm{2}}$ (Fig. A1). Such a value of *p* assumes
that *h*∼0.7, which is towards the higher end for observed Hack
exponents, and that the river catchment is very elongated. When plotting the
logarithm of the model slope against drainage area
(Fig. A1b), for which area is given by $a={x}^{p}/{k}_{\mathrm{w}}$ and assuming *k*_{w}=1, for the simple stream
power law derived here the slope–area exponent $\mathit{\theta}=-m$. The value of
the dimensional constant *k* has no impact on the slope–area exponent as
expected. The transport model likewise creates river-long profiles that have
on average a negative curvature. For small values of *x* there is, however, a
region of positive curvature in which *κ*>*c**k*_{w}*α*^{δ}*L*^{δp}. For the slope–area analysis this
leads to a positive gradient in the trend for small catchment
areas. This relationship subsequently has a negative slope for larger
catchments. The point of inflection is dependent on the value of
*D*_{e}; for smaller values of *κ* the region of positive
gradient is reduced. There is therefore a critical catchment area that is
dependent on the diffusive term *κ*. After this critical point the slope–area relationship becomes negative. At distances down-system, where the
upstream area is greater than this critical area, the gradient $\mathit{\theta}=-\mathrm{0.88}$; *θ* is insensitive to the coefficient *c* as would be expected.

The range of gradients found for river catchments for this type of slope–area
analysis, usually referred to as concavity, generally lies within the range
$\mathit{\theta}=-\mathrm{0.35}$ to −0.70 (Snyder et al., 2000; Wobus et al., 2006). It is
trivial to find the values of *m* for the steady-state solution to the stream
power law that fit such values of *θ*. To further explore how
*θ* depends on *δ* and *p* within the transport model we solve
Eq. (8) numerically for *δ*=1, 1.5, and 2 while
keeping *h*=0.7 or 0.6 (Fig. A2). The result is that
*θ* varies from −0.3 for the case of *δ*=1 to −1.31 for
*δ*=2. The values of the gradient for the slope–area analysis for $\mathrm{1.4}<p<\mathrm{2}$, in which we assume *d*=1 and hence $p=\mathrm{1}/h$, are displayed in Table
A1. For the transport model the slope is dependent on
both *δ* and *p*.

Clearly there is a combination of *δ**p* values that is equally capable of
fitting the observed river-long profile. Furthermore, for the transport model
the slope is a function of the Hack exponent *h* (and therefore *p*) and the
choice of *δ*. This because of the diffusivity term that leads to
positive curvature and rounded 1-D profiles (Figs. A1b and
A2). The magnitude of the water flux term within the
transport equation (Eq. 8) is dependent on how much water
the river network captures, which is in turn a function of how elongated the
catchment is.

The positive slope–area relationship for the transport model for small
catchment areas (see Figs. A1b and
A2b) has been previously explored in
Willgoose et al. (1991). The gradient of the relationship between slope
and catchment area is dominantly a function of the exponent *δ* within
Eq. (6). The value of this exponent is likely within the
range of $\mathrm{1}<\mathit{\delta}<\mathrm{2}$ depending on the bedload transport law assumed
(Armitage et al., 2013). If the observations of trunk river slope
against catchment area are representative of a landscape at steady state,
then for the smaller range of $\mathrm{1}<\mathit{\delta}\le \mathrm{1.5}$, a realistic catchment
topography can be generated.

The authors declare that they have no conflict of interest.

We would like to thank Guy Simpson (University of Geneva) for sharing his
numerical model that solves the Smith and Bretherton (1972) equations in two
dimensions. We thank Tom Dunkley Jones (University of Birmingham) for
discussions on the duration of conglomeratic deposition during the PETM,
Gareth Roberts (Imperial College London) for discussions on slope exponents,
and Jean Braun (GFZ Potsdam) for sharing his numerical model Fastscape. This
work was initiated under funding by the Royal Astronomical Society through a
research fellowship awarded to John Armitage while he was at Royal Holloway
University of London.

Edited by: Jean Braun

Reviewed by: Laure Guerit and one anonymous referee

Allen, P. A.: Sediment Routing Systems: The Fate of Sediment from Source to Sink, Cambridge University Press, 2017. a, b

Allen, P. A., Armitage, J. J., Carter, A., Duller, R. A., Michael, N. A.,
Sinclair, H. D., Whitchurch, A. L., and Whittaker, A. C.: The *Q*_{s}
problem: Sediment volumetric balance of proximal foreland basin systems,
Sedimentology, 60, 102–130, https://doi.org/10.1111/sed.12015, 2013. a, b, c, d

Armitage, J. J., Duller, R. A., Whittaker, A. C., and Allen, P. A.: Transformation of tectonic and climatic signals from source to sedimentary archive, Nat. Geosci., 4, 231–235, https://doi.org/10.1038/ngeo1087, 2011. a, b, c, d, e, f, g

Armitage, J. J., Dunkley Jones, T., Duller, R. A., Whittaker, A. C., and Allen, P. A.: Temporal buffering of climate-driven sediment flux cycles by transient catchment response, Earth Planet. Sci. Lett., 369–370, 200–210, https://doi.org/10.1016/j.epsl.2013.03.020, 2013. a, b, c, d, e, f, g, h, i, j, k

Armitage, J. J., Allen, P. A., Burgess, P. M., Hampson, G. J., Whittaker, A. C., Duller, R. A., and Michael, N. A.: Physical stratigraphic model for the Eocene Escanilla sediment routing system: Implications for the uniqueness of sequence stratigraphic architectures, J. Sediment. Res., 85, 1510–1524, https://doi.org/10.2110/jsr.2015.97, 2015. a, b

Baldwin, J. A., Whipple, K. X., and Tucker, G. E.: Implications of the shear stress river incision model for the timescale of postorogenic decay of topography, J. Geophys. Res., 108, 2158, https://doi.org/10.1029/2001JB000550, 2003. a, b, c, d, e

Banavar, J. R., Colaiori, F., Flammini, A., Giacometi, A., MAritan, A., and Rinaldo, A.: Sculpting of a fractal river basin, Phys. Rev. Lett., 23, 4522–4525, 1997. a

Bonnet, S. and Crave, A.: Landscape response to climate change: Insights from experimental modeling and implications for tectonic versus climatic uplift of topograph, Geology, 31, 123–126, 2003. a

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

Braun, J., Voisin, C., Gourlan, A. T., and Chauvel, C.: Erosional response of an actively uplifting mountain belt to cyclic rainfall variations, Earth Surf. Dynam., 3, 1–14, https://doi.org/10.5194/esurf-3-1-2015, 2015. a

Campforts, B. and Covers, G.: Keeping the edge: A numerical method that avoids knickpoint smearing when solving the stream power law, J. Geophys. Res., 120, 1189–1205, https://doi.org/10.1002/2014JF003376, 2015. a, b

Castelltort, S. and Van Den Dreissche, J.: How plausible are high-frequency sediment supply-driven cycles in the stratigraphic record?, Sediment. Geol., 157, 3–13, https://doi.org/10.1016/S0037-0738(03)00066-6, 2003. a

Cowie, P. A., Whittaker, A. C., Attal, M., Roberts, G., Tucker, G. E., and Ganas, A.: New constraints on sediment-flux–dependent river incision: Implications for extracting tectonic signals from river profiles, Geology, 36, 535–538, https://doi.org/10.1130/G24681A.1, 2008. a

Croissant, T. and Braun, J.: Constraining the stream power law: a novel approach combining a landscape evolution model and an inversion method, Earth Surf. Dynam., 2, 155–166, https://doi.org/10.5194/esurf-2-155-2014, 2014. a

Crosby, B. T., Whipple, K. X., Gasparini, N. M., and Wobus, C. W.: Formation of fluvial hanging valleys: Theory and simulation, J. Geophys. Res., 112, F03510, https://doi.org/10.1029/2006JF000566, 2007. a

Cuevas, J. L.: Estratigrafia del “Garumniense” de la Conca de Tremp, Prepirineo de Lerida, Acata Geological Hispanica, 27, 95–108, 1992. a

D'Arcy, M., Whittaker, A. C., and Roda-Boluda, D. C.: Measuring alluvial fan sensitivity to past climate changes using a self-similarity approach to grain-size fining, Death Valley, California, Sedimentology, 64, 388–424, https://doi.org/10.1111/sed.12308, 2016. a

D'Arcy, M., Whittaker, A. C., and Roda-Boluda, D. C.: Measuring alluvial fan sensitivity to past climate changes using a self-similarity approach to grain-size fining, Death Valley, California, Sedimentology, 64, 388–424, https://doi.org/10.1111/sed.12308, 2017. a, b, c

Demoulin, A., Mather, A., and Whittaker, A.: Fluvial archives, a valuable record of vertical crustal deformation, Quaternary Sci. Rev., 166, 10–37, https://doi.org/10.1016/j.quascirev.2016.11.011, 2017. a, b, c

Densmore, A. L., Allen, P. A., and Simpson, G.: Development and response of a coupled catchment fan system under changing tectonic and climatic forcing, J. Geophys. Res., 112, F01002, https://doi.org/10.1029/2006JF000474, 2007. a

Dietrich, W. E., Bellugi, D. G., Sklar, L. S., Srock, J. D., Heimsath, A. M., and Roering, J. J.: Geomorphic transport laws for predicting landscape form and dynamics, in: Prediction in Geomorphology, edited by: Wilcock, P. R. and Iverson, R. M., Geoph. Monog. Series, 135, 1–30, https://doi.org/10.1029/135GM09, 2003. a

Dietrich, W. M.: Settling velocity of natural particles, Water Resour. Res., 18, 1615–1626, 1982. a

Dodds, P. S. and Rothman, D. H.: Geometry of river networks. I. Scaling, fluctuations and deviations, Phys. Rev. E, 63, 016115, https://doi.org/10.1103/PhysRevE.63.016115, 2000. a, b, c

Dunkley Jones, T., Ridgwell, A., Lunt, D. J., Maslin, M. A., Schmidt, D. N., and Valdes, P. J.: A Palaeogene perspective on climate sensitivity and methane hydrate instability, P. T. Roy. Soc. Pt. A, 368, 2395–2415, https://doi.org/10.1098/rsta.2010.0053, 2010. a

Foreman, B. Z.: Climate-driven generation of a fluvial sheet sand body at the Paleocene-Eocene boundary in north-west Wyoming (USA), Basin Res., 26, 225–241, https://doi.org/10.1111/bre.12027, 2014. a, b, c, d, e, f

Foreman, B. Z., Heller, P. L., and Clementz, M. T.: Fluvial response to abrupt global warming at the Palaeocene/Eocene boundary, Nature, 491, 92–95, https://doi.org/10.1038/nature11513, 2012. a, b, c, d, e, f, g, h, i

Forman, B. Z. and Straub, K. M.: Autogenic geomorphic processes determine the resolution and fidelity of terrestrial paleoclimate records, Sci. Adv., 3, e1700683, https://doi.org/10.1126/sciadv.1700683, 2017. a, b

Ganti, V., Lamb, M. P., and McElroy, B.: Quantitative bounds on morphodynamics and implications for reading the sedimentary record, Nat. Comun., 5, 3298, https://doi.org/10.1038/ncomms4298, 2014. a, b

Grimaud, J. L., Chardon, D., and Beauvais, A.: Very long-term incision of big rivers, Earth Planet. Sci. Lett., 405, 74–84, https://doi.org/10.1016/j.epsl.2014.08.021, 2014. a

Grimaud, J. L., Paola, C., and Voller, V.: Experimental migration of knickpoints: influence of style of base-level fall and bed lithology, Earth Surf. Dynam., 4, 11–23, https://doi.org/10.5194/esurf-4-11-2016, 2016. a

Guerit, L., Métivier, F., Devauchelle, O., Lajeunesse, E., and Barrier, L.: Laboratory alluvial fans in one dimension, Phys. Rev. E, 90, 022203, https://doi.org/10.1103/PhysRevE.90.022203, 2014. a

Hack, J. T.: Studies of longitudingal profiles in Virginia and Maryland, US Geological Survey Proffesional Papaer, 294-B, 1957. a

Hancock, G. R., Coulthard, T. J., and Lowry, J. B. C.: Predicting uncertainty in sediment transport and landscape evolution - the influence of initial surface conditions, Comput. Geosci., 90, 117–130, https://doi.org/10.1016/j.cageo.2015.08.014, 2016. a

Howard, A. D. and Kerby, G.: Channel changes in badlands, Geol. Soc. Am. Bull., 94, 739–752, 1983. a, b

Izumi, N. and Parker, G.: Inception of channelization and drainage basin formation: upstream-driven theory, J. Fluid Mech., 283, 341–363, 1995. a

Kirkby, M. J.: Hillslope process-response models based on the continuity equation, Special Publication Institute of British Geographers, 3, 15–30, 1971. a

Lacey, G.: Stable channels in alluvium, in: Minutes of the Proceedings, edited by: Grierson, W. W., Institution of Civil Engineers Publishing, vol. 229, 259–292, 1930. a

Lague, D.: The stream power river incision model: evidence, theory and beyond, Earth Surf. Proc. Land., 39, 38–61, https://doi.org/10.1002/esp.3462, 2014. a, b, c, d

Leopold, L. B. and Maddock, T.: The hydraulic geometry of stream channels and some physiographic implications, Tech. Rep. 252, US Geological Survey Proffesional Papers, 57 pp., 1953. a

Liu, Y., Métivier, F., Gaillerdet, J., Ye, B., Meunier, P., Narteau, C., Lajeunesse, E., Han, T., and Malverti, L.: Erosion rates deduced from seasonal mass balance along the upper Urumqi River in Tianshan, Solid Earth, 2, 283–301, https://doi.org/10.5194/se-2-283-2011, 2011. a

Manners, H. R., Grimes, S. T., Sutton, P. A., Domingo, L., Leng, M. J., Twitchett, R. J., Hart, M. B., Dunkley Jones, T., Pancost, R. D., Duller, R., and Lopez-Martinez, N.: Magnitude and profile of organic carbon isotope records from the Paleocene – Eocene Thermal Maximum: Evidence from northern Spain, Earth Planet. Sci. Lett., 376, 220–230, https://doi.org/10.1016/j.epsl.2013.06.016, 2013. a, b

Maritan, A., Rinaldo, A., Rigon, R., Giacometti, A., and Rordriguez-Iturbe, I.: Scaling laws for river networks, Phys. Rev. E, 53, 1510–1515, 1996. a

Meunier, P., Métivier, F., Lajeunesse, E., Mériaux, A. S., and Faure, J.: Flow pattern and sediment transport in a braided river: The ”torrent de St Pierre” (French Alps), J. Hydrol., 330, 496–505, https://doi.org/10.1016/j.jhydrol.2006.04.009, 2006. a

Michael, N. A., Whittaker, A. C., Carter, A., and Allen, P. A.: Volumetric budget and grain-size fractionation of a geological sediment routing system: Eocene Escanilla Formation, south-central Pyrenees, Geol. Soc. Am. Bull., 126, 585–599, https://doi.org/10.1130/B30954.1, 2014. a

Mouchené, M., van der Beek, P., Carretier, S., and Mouthereau, F.: Autogenic versus allogenic controls on the evolution of a coupled fluvial megafan-mountainous catchment system: numerical modelling and comparison with the Lannemezan megafan system (northern Pyrenees, France), Earth Surf. Dynam., 5, 125–143, https://doi.org/10.5194/esurf-5-125-2017, 2017. a

Paola, C., Heller, P. L., and Angevine, C. L.: The large-scale dynamics of grain-size variation in alluvial basins, 1: Theory, Basin Res., 4, 73–90, https://doi.org/10.1111/j.1365-2117.1992.tb00145.x, 1992. a, b, c, d, e

Pastor-Satorras, R. and Rothman, D. H.: Stochastic Equation for the Erosion of Inclined Topography, Phys. Rev. Lett., 80, 4349–4352, https://doi.org/10.1103/PhysRevLett.80.4349, 1998. a

Perron, J. T. and Royden, L.: An integral approach to bedrock river profile analysis, Earth Surf. Proc. Land., 38, 570–576, https://doi.org/10.1002/esp.3302, 2012. a, b

Pritchard, D., Roberts, G. G., White, N. J., and Richardson, C. N.: Uplift histories from river profiles, Geophys. Res. Lett., 36, L24301, https://doi.org/10.1029/2009GL040928, 2009. a, b

Rigon, R., Rodriguez-Iturbe, I., Maritan, A., Giacometti, A., Tarboton, D. G., and Rinaldo, A.: On Hack's law, Water Resour. Res., 32, 3367–3374, 1996. a, b

Rohais, S., Bonnet, S., and Eschard, R.: Sedimentary record of tectonic and climatic erosional perturbations in an experimental coupled catchment-fan system, Basin Res., 23, 1–15, https://doi.org/10.1111/j.1365-2117.2011.00520.x, 2011. a

Röhl, U., Westerhold, T., Bralower, T. J., and Zachos, J. C.: On the duration of the Paleocene-Eocene thermal maximum (PETM), Geochem. Geophy. Geosy., 8, Q12002, https://doi.org/10.1029/2007GC001784, 2007. a

Romans, B. W., Castelltort, S., Covault, J. A., Fildani, A., and Walsh, J. P.: Environmental signal propagation in sedimentary systems across timescales, Earth-Sci. Rev., 153, 7–29, https://doi.org/10.1016/j.earscirev.2015.07.012, 2016. a, b, c

Roy, S. G., Koons, P. O., Upton, P., and Tucker, G. E.: The influence of crustal strength fields on the patterns and rates of fluvial incision, J. Geophys. Res., 120, 275–299, https://doi.org/10.1002/2014JF003281, 2015. a, b

Rudge, J. F., Roberts, G. G., White, N. J., and Richardson, C. N.: Uplift histories of Africa and Australia from linear inverse modeling of drainage inventories, J. Geophys. Res., 120, 894–914, https://doi.org/10.1002/2014JF003297, 2015. a, b, c, d

Schmitz, B. and Pujalte, V.: Abrupt increase in seasonal extreme precipitation at the Paleocene-Eocene boundary, Geology, 35, 215–218, https://doi.org/10.1130/G23261A.1, 2007. a, b, c

Schwanghart, W. and Scherler, D.: Short Communication: TopoToolbox 2 – MATLAB-based software for topographic analysis and modeling in Earth surface sciences, Earth Surf. Dynam., 2, 1–7, https://doi.org/10.5194/esurf-2-1-2014, 2014. a

Simpson, G.: Practical Finite Element Modeling in Earth Science using Matlab, Wiley-Blackwell, 2017. a, b

Simpson, G. and Castelltort, S.: Model shows that rivers transmit high-frequency climate cycles to the sedimentary record, Geology, 40, 1131–1134, https://doi.org/10.1130/G33451.1, 2012. a, b

Simpson, G. and Schlunegger, F.: Topographic evolution and morphology of surfaces evolving in response to coupled fluvial and hillslope sediment transport, J. Geophys. Res., 108, 2300, https://doi.org/10.1029/2002JB002162, 2003. a, b

Sluijs, A., Schouten, S., Pagani, M., Woltering, M., Brinkhuis, H., Sinninghe Damsté, J. S., Dickens, G. R., Huber, M., Reichart, G. J., Stein, R., Matthiessen, J., Lourens, L. J., Pedentchouk, N., Backman, J., Moran, K., and the Expedition 302 scientists: Subtropical Arctic Ocean temperatures during the Palaeocene/Eocene thermal maximum, Nature, 411, 610–613, https://doi.org/10.1038/nature04668, 2006. a

Smith, T. R.: A theory for the emergence of channelized drainage, J. Geophys. Res., 115, F02023, https://doi.org/10.1029/2008JF001114, 2010. a

Smith, T. R. and Bretherton, F. P.: Stability and conservation of mass in drainage basin evolution, Water Resour. Res., 8, 1506–1529, https://doi.org/10.1029/WR008i006p01506, 1972. a, b, c, d, e, f, g, h, i

Smith, T. R., Merchant, G. E., and Birnir, B.: Transient attractors: towards a theory of the graded stream for alluvial and bedrock channels, Comput. Geosci., 26, 541–580, https://doi.org/10.1016/S0098-3004(99)00128-4, 2000. a, b

Snyder, N. P., Whipple, K. X., Tucker, E., and Merrits, D. J.: Landscape response to tectonic forcing: Digital elevation model analysis of stream profiles in the Mendocino triple junction region, northern California, Geol. Soc. Am. Bull., 112, 1250–1263, https://doi.org/10.1130/0016-7606(2000)112<1250:LRTTFD>2.0.CO;2, 2000. a, b, c, d

Tarboton, D. G., Bras, R. L., and Rodriguez-Iturbe, I.: Comment on “On the fractal dimension of stream network” by Paolo La Barbera and Renzo Rosso, Water Resour. Res., 26, 2243–2244, https://doi.org/10.1029/WR026i009p02243, 1990. a

Temme, A. J. A. M., Armitage, J. J., Attal, M., van Gorp, W., Coulthard, T. J., and Schoorl, J. M.: Choosing and using landscape evolution models to inform field stratigraphy and landscape reconstruction studies, Earth Surf. Proc. Land., 42, 2167–2183, https://doi.org/10.1002/esp.4162, 2017. a, b, c

Tinkler, K. J. and Whol, E. E.: Rivers over Rock: fluvial processes in bedrock channels, American Geophysical Union, Washington, USA, Geophys. Monog. Series, 107, https://doi.org/10.1029/GM107, 1998. a

Tucker, G. E. and Whipple, K. X.: Topographic outcomes predicted by stream erosion models: Sensitivity analysis and intermodel comparison, J. Geophys. Res., 107, 2179, https://doi.org/10.1029/2001JB000162, 2002. a, b, c, d

Valla, P. G., van der Beek, P. A., and Lague, D.: Fluvial incision into bedrock: Insights from morphometric analysis and numerical modeling of gorges incising glacial hanging valleys (Western Alps, France), J. Geophys. Res., 116, F02010, https://doi.org/10.1029/2008JF001079, 2010. a, b

van der Beek, P. and Bishop, P.: Ceonzoic river profile development in the Upper Lachlan catchment (SE Australia) as a test of quantative fluvial incision models, J. Geophys. Res., 108, 2309, https://doi.org/10.1029/2002JB002125, 2003. a

Whipple, K. X.: Fluvial landscape response time: how plausible is steady-state denudation, Am. J. Sci., 301, 313–325, 2001. a, b, c, d, e, f

Whipple, K. X.: Bedrock rivers and the geomorphology of active orogens, Ann. Rev. Earth Planet. Sci., 32, 151–185, https://doi.org/10.1146/annurev.earth.32.101802.120356, 2004. a

Whipple, K. X. and Meade, B. J.: Orogen response to changes in climatic and tectonic forcing, Earth Planet. Sci. Lett., 243, 218–228, https://doi.org/10.1016/j.epsl.2005.12.022, 2006. a, b

Whipple, K. X. and Tucker, G. E.: Dynamics of the stream-power river incision model: Implications for height limits of mountain ranges, landscape response times, and research needs, J. Geophys. Res., 104, 17661–17674, 1999. a, b

Whipple, K. X. and Tucker, G. E.: Implications of sediment-flux-dependent river incision models for landscape evolution, J. Geophys. Res., 107, https://doi.org/10.1029/2000JB000044, 2002. a, b, c, d, e

Whittaker, A. C. and Boulton, S. J.: Tectonic and climatic controls on knickpoint retreat rates and landscape response times, J. Geophys. Res., 117, F02024, https://doi.org/10.1029/2011JF002157, 2012. a, b, c

Whittaker, A. C., Cowie, P. A., Attal, M., Tucker, G. E., and Roberts, G. P.: Bedrock channel adjustment to tectonic forcing: Implications for predicting river incision rates, Geology, 35, 103–106, https://doi.org/10.1130/G23106A.1, 2007. a, b

Whittaker, A. C., Attal, M., Cowie, P. A., Tucker, G. E., and Roberts, G.: Decoding temporal and spatial patterns of fault uplift using transient river long-profiles, Geomorphology, 100, 506–526, https://doi.org/10.1016/j.geomorph.2008.01.018, 2008. a

Whittaker, A. C., Duller, R. A., Springett, J., Smithells, R. A., Whitchurch, A. L., and Allen, P. A.: Decoding downstream trends in stratigraphic grain size as a function of tectonic subsidence and sediment supply, Geol. Soc. Am. Bull., 123, 1363–1382, 2011. a, b

Willett, S. D., McCoy, S. W., Perron, J. T., Goren, L., and Chen, C.: Dynamic Reorganization of River Basins, Science, 343, 6175, https://doi.org/10.1126/science.1248765, 2014. a

Willgoose, G., Bras, R. L., and Rodriguez-Iturbe, I.: A physical explanation of an observed link area-slope relationship, Water Resour. Res., 27, 1697–1702, 1991. a

Wobus, C., Whipple, K. X., Kirby, E., Snyder, N., Johnson, J., Spyropolou, K., Crosby, B., and Sheehan, D.: Tectonics from topography: Procedures, promise, and pitfalls, Geol. Soc. Am. Sp., 398, 55–74, 2006. a, b, c