the Creative Commons Attribution 4.0 License.

the Creative Commons Attribution 4.0 License.

# Statistical modelling of co-seismic knickpoint formation and river response to fault slip

### Thomas Croissant

### Edwin Baynes

### Dimitri Lague

Most landscape evolution models adopt the paradigm of constant and uniform uplift. It results that the role of fault activity and earthquakes on landscape building is understood under simplistic boundary conditions. Here, we develop a numerical model to investigate river profile development subjected to fault displacement by earthquakes and erosion. The model generates earthquakes, including mainshocks and aftershocks, that respect the classical scaling laws observed for earthquakes. The distribution of seismic and aseismic slip can be partitioned following a spatial distribution of mainshocks along the fault plane. Slope patches, such as knickpoints, induced by fault slip are then migrated at a constant rate upstream a river crossing the fault. A major result is that this new model predicts a uniform distribution of earthquake magnitude rupturing a river that crosses a fault trace and in turn a negative exponential distribution of knickpoint height for a fully coupled fault, i.e. with only co-seismic slip. Increasing aseismic slip at shallow depths, and decreasing shallow seismicity, censors the magnitude range of earthquakes cutting the river towards large magnitudes and leads to less frequent but higher-amplitude knickpoints, on average. Inter-knickpoint distance or time between successive knickpoints follows an exponential decay law.

Using classical rates for fault slip (15 mm year^{−1}) and knickpoint
retreat (0.1 m year^{−1}) leads to high spatial densities of knickpoints. We
find that knickpoint detectability, relatively to the resolution of
topographic data, decreases with river slope that is equal to the ratio
between fault slip rate and knickpoint retreat rate. Vertical detectability
is only defined by the precision of the topographic data that sets the lower
magnitude leading to a discernible offset. Considering a retreat rate with a
dependency on knickpoint height leads to the merging of small knickpoints
into larger ones and larger than the maximum offset produced by individual
earthquakes. Moreover, considering simple scenarios of fault burial by
intermittent sediment cover, driven by climatic changes or linked to
earthquake occurrence, leads to knickpoint distributions and river profiles
markedly different from the case with no sediment cover. This highlights the
potential role of sediments in modulating and potentially altering the
expression of tectonic activity in river profiles and surface topography.
The correlation between the topographic profiles of successive parallel
rivers cutting the fault remains positive for distance along the fault of
less than half the maximum earthquake rupture length. This suggests that
river topography can be used for paleo-seismological analysis and to assess
fault slip partitioning between aseismic and seismic slip. Lastly, the
developed model can be coupled to more sophisticated landscape evolution
models to investigate the role of earthquakes on landscape dynamics.

The interactions among tectonics, climate and surface processes govern the evolution of the Earth's topography (e.g. Willet, 1999; Whipple, 2009). Among the potential link and feedbacks between tectonics and surface processes, the building of topographic slopes by tectonic deformation is critical. Erosion rates and most geomorphological processes are strongly sensitive to local slope, including river incision (e.g. Whipple and Tucker, 1999), glacial carving (e.g. Herman and Braun, 2008), soil creep (e.g. McKean et al., 1993) and hillslope mass wasting (e.g. Keefer, 1994). The dependency on slope can be linear or non-linear, mainly due to threshold effects or to a power-law behaviour. For instance, a theoretical model combined with a data compilation suggests that river incision rate is linearly dependent on slope at knickpoints and more than linearly dependent on slope for more gentle stream profiles (Lague, 2014). This is pivotal, as temporal variations in tectonic displacement and in slope building cannot be averaged out when considering river profile evolution using an erosion law with a non-linear dependency on slope. In addition to slope, the height of knickpoints (i.e. with a slope above average local slope) and waterfalls (i.e. with a slope close to infinity) appears as a fundamental ingredient of their survival, retreat rate and river incision (Hayakawa and Matsukura, 2003; Baynes et al., 2015; Scheingross and Lamb, 2017). Similar issues arise for hillslope dynamics impacted by fault scarp development (Arrowsmith et al., 1996) and possibly for faults in glaciated landscapes. Despite this, most landscape evolution models of topographic growth consider slope building as a continuous process resulting from a constant (or smoothly varying) uplift rate (e.g. Braun and Willett, 2013; Thieulot et al., 2014; Campforts et al., 2017). There is therefore a clear need to define how tectonic deformation builds topographic slopes in numerical models.

The expression of tectonic deformation on topographic slope is diverse, and its spatial and temporal scales range from metres to continents and from instantaneous to geological times, respectively. Tectonic deformation can (1) instantaneously generate steep-to-infinite slopes when earthquakes rupture the Earth's surface (e.g. Wells and Coppersmith, 1994); (2) induce progressive slope building at the orogen scale and over a seismic cycle by aseismic deformation (i.e. deformation not associated with earthquakes) and interseismic deformation (i.e. deformation occurring in between large-magnitude earthquakes) (e.g. Cattin and Avouac, 2000) or by the deformation associated with earthquakes with no surface rupture; and (3) lead to longer-term topographic tilting at the orogen-to-continental scale by isostatic readjustment (e.g. Watts, 2001) or viscous mantle flow (e.g. Braun, 2010). In this paper, we focus on the building of topographic slopes by fault slip at the intersection between a fault trace and a river. This is motivated first by the fact that the greatest slopes are expected to occur by faulting and second by the already-well-understood role of isostasy and viscous deformation on topography (e.g. Watts, 2001; Braun, 2010). In active mountain belts, displacement along frontal thrust faults can lead to the development of co-seismic waterfalls, knickpoints and knickzones than can reach several metres in elevation (e.g. Boulton and Whittaker, 2009; Yanites et al., 2010; Cook et al., 2013). These differential topographies, associated with high slopes, are referred to as slope patches in the following work. These slope patches have long been recognized as potential markers of the dynamic response of rivers (e.g. Gilbert, 1896) to transient conditions, not limited to changes in tectonic activity, and including base level fall and lithological contrasts, among others. Yet, in active tectonic areas, knickpoints are frequently associated with fault activity and transience in uplift rate (e.g. van der Beek et al., 2001; Quigley et al., 2006; Dorsey and Roering, 2006; Yildirim et al., 2011). These slope patches generated by frontal thrusts along a river migrate upstream by erosion and are expected to set the erosion rate of the entire landscape (Rosenbloom and Anderson, 1994; Royden and Perron, 2013; Yanites et al., 2010; Cook et al., 2013).

Fault slip and surface rupture classically occur by seismic slip during earthquakes. However, associating individual earthquakes with knickpoints or associating series of knickpoints with series of earthquakes remains challenging from field data. We therefore use in this paper a statistical model of earthquakes to simulate the expected slope and height distributions of the slope patches generated by earthquakes (i.e. fault seismic slip) and fault aseismic slip at the intersection between a thrust fault and a river. This model uses the branching aftershock sequence (BASS) model (Turcotte et al., 2007) to simulate temporal and spatial series of earthquakes based on the main statistical and scaling laws of earthquakes. The rupture extent and displacement of earthquakes are inferred using classical scaling laws (Leonard, 2010). We focus on the response of rivers and analyse the resulting knickpoint height distribution and their migration distance along a single river in near-fault conditions. We also infer the correlation between the topography of successive parallel rivers distributed along the strike of a single fault. The obtained results are then discussed with regards to the potential of knickpoints and waterfalls to offer paleo-seismological constraints and to the necessity of considering time-variable uplift accounting for earthquake sequences in landscape evolution models. It is important to stress that this study does not aim to investigate specific geomorphological settings but to give general theoretical and modelling arguments to the interpretation of river profiles upstream of active faults.

## 2.1 From fault slip and earthquakes to surface ruptures and knickpoints

In near-fault conditions, too few data characterizing fault rupture geometry at one location (e.g. along a river) exist to assess the distribution of the slope and height of surface ruptures resulting from earthquakes by local fault activity (e.g. Ewiak et al., 2015; Wei et al., 2015; Sun et al., 2016). Regional or global compilation of fault rupture by earthquakes (e.g. Wells and Coppersmith, 1994; Leonard, 2010; Boncio et al., 2018) offers another approach, that yet suffers from inescapable statistical biases mainly due to the use of faults with different slip rates, dimensions, seismogenic properties and records of paleo-earthquakes. In addition, small earthquakes associated with small rupture extents and co-seismic displacement are less likely to be identified in the field. For instance, using seismological scaling laws (Leonard, 2010), an earthquake of magnitude 3 on a thrust fault has a rupture length of 188 m and an average displacement of 1.2 cm. This displacement is clearly below the precision of current digital elevation models or in any case hidden by the inherent topographic roughness.

Statistical or theoretical inferences offer another means to associate fault activity and earthquakes with surface ruptures and knickpoints. Earthquakes tend to universally follow the Gutenberg–Richter frequency–magnitude distribution in Eq. (1):

where *M*_{w} is the magnitude, *N*(≥*M*_{w}) is the number of earthquakes with
magnitudes greater than or equal to *M*_{w}, *b* is the exponent of the tail
(referred to as the *b* value), generally observed to be close to 1
($\mathrm{0.5}<b<\mathrm{1.5}$), and *a* characterizes earthquake
productivity (Gutenberg and Richter, 1944). The definitions of all variables
used in this paper are summarized in Table C1. Assuming self-similarity, a
*b* value of 1 can be interpreted as the result of the successive segmentation
of larger earthquakes into smaller earthquakes (Aki, 1981; King, 1983) so
that any point along a 2-D fault plane, including the intersection between
the fault trace and a river, displays a uniform probability to be ruptured
by earthquakes of any magnitude. This inference only stands if the
distribution of earthquakes along the fault plane is uniform. However, fault
slip can occur by seismic slip but also by aseismic deformation, including
interseismic creep, post-seismic deformation and slow slip events (e.g.
Scholz, 1998; Peng and Gomberg, 2010; Avouac, 2015). The relative spatial
and temporal distribution of aseismic and seismic slip along a fault plane
is variable and still poorly understood. Yet, experimental results and the
depth distribution of earthquakes along subduction or intraplate thrust
faults suggest that shallow depths (<5 km) are favourable to
frictional stability and in turn to aseismic slip (Scholz, 1998). This
probably censors the magnitude range of earthquakes rupturing the surface
towards large magnitudes associated with rupture extent greater than this
minimum seismogenic depth.

## 2.2 Knickpoint formation

The transformation of surface ruptures into knickpoints remains a relatively enigmatic issue. Linking knickpoints to individual earthquakes is challenging, although some recently formed knickpoints have been clearly identified as the result of the surface rupture of a single large earthquake (e.g. Yanites et al., 2010; Cook et al., 2013). The transformation of individual surface ruptures into individual co-seismic knickpoints is not necessarily a bijective function and is more likely to be a surjective function. In other words, a knickpoint can be made of several surface ruptures. Indeed, if the time interval between two (or more) successive ruptures at the same location is less than a characteristic migration time required to segregate their topographic expressions, then the formed knickpoint will result from this succession of surface ruptures and earthquakes. An end-member setting favouring this behaviour is the case of fault scarps developing on hillslopes, whose degradation is generally assumed to follow a diffusion law (e.g. Nash, 1980; Avouac, 1993; Arrowsmith et al., 1996; Roering et al., 1999; Tucker and Bradley, 2010). Moreover, in the downstream part of rivers, fault scarps can remain buried under a sediment cover due to, for instance, the development of an alluvial fan (Finnegan and Balco, 2013; Malatesta and Lamb, 2018). Development of the fault scarp height by successive ruptures or the thinning of the alluvial cover can then expose the scarp, in turn potentially forming a knickpoint that can erode and migrate. This intermittent fault-burial mechanism can therefore produce knickpoints formed by the surface rupture of several earthquakes.

The burial of the fault during successions of aggradation–incision phases of an alluvial fan located immediately downstream of the fault (e.g. Carretier and Lucazeau, 2005) has not been considered in previous landscape evolution models. This mechanism is suggested to be a primary control of knickpoint and waterfall formation by allowing the merging of several small co-seismic scarps formed during burial phases into single high-elevation waterfalls that migrate during latter incision phases (Finnegan and Balco, 2013; Malatesta and Lamb, 2018).

## 2.3 Knickpoint migration and preservation

Once formed, knickpoints can migrate upstream due to river erosion. Over
geological timescales (>10^{3} years), rates of knickpoint
retreat for bedrock rivers typically range between $\sim {\mathrm{10}}^{-\mathrm{3}}$
and $\sim {\mathrm{10}}^{-\mathrm{1}}$ m year^{−1} (e.g. Van Heijst and Postma,
2001). This range is also consistent with the order of magnitude of
documented knickpoint retreat rates in eastern Scotland (Bishop et al.,
2005; Jansen et al., 2011), around $\sim {\mathrm{10}}^{-\mathrm{1}}$ m year^{−1}, in
the central Apennines, Italy, and in the Hatay Graben, southern Turkey
(Whittaker and Boulton, 2012), between $\sim {\mathrm{10}}^{-\mathrm{3}}$ and
$\sim {\mathrm{10}}^{-\mathrm{2}}$ m year^{−1}. However, on shorter timescales,
significantly higher rates can be found with values potentially reaching
∼10^{0} or even ∼10^{1} m year^{−1}. For
instance, Niagara Falls retreated at a rate of a few metres per year
over tens of years (Gilbert, 1907) and some knickpoints formed by the 1999
Chi-Chi earthquake in Taiwan even retreated by a few hundreds of metres over
about 10 years (Yanites et al., 2010; Cook et al., 2013). A more extensive
analysis of the range of knickpoint retreat rates in relation to the
observation timescale can be found in Van Heijst and Postma (2001) and in
Loget and Van Den Driessche (2009).

In detachment-limited conditions, the stream power incision model predicts that knickpoint horizontal migration or retreat follows a linear or non-linear kinematic wave in the upstream direction, depending on the slope exponent (e.g. Rosenbloom and Anderson, 1994; Tucker and Whipple, 2002; Whittaker and Boulton, 2012; Royden and Perron, 2013). This prediction is supported by the apparent correlation between retreat rate and drainage area or water discharge, deduced from field observation and experimental studies (Parker, 1977; Schumm et al., 1987; Rosenbloom and Anderson, 1994; Bishop et al., 2005; Crosby and Whipple, 2006; Loget et al., 2006; Berlin and Anderson, 2007). However, some experimental results show no dependency of retreat rate on water discharge (Holland and Pickup, 1976), possibly due to the self-regulatory response of river geometry to water discharge through change in river channel width (Baynes et al., 2018). Other factors influencing retreat rate include, among others, sediment discharge (e.g. Jansen et al., 2011; Cook et al., 2013), flood events (e.g. Baynes et al., 2015), rock strength (e.g. Stock and Montgomery, 1999; Hayakawa and Matsukura, 2003; Baynes et al., 2018), fracture density and orientation (Antón et al., 2015; Brocard et al., 2016) and the spacing and height of the waterfalls (Scheingross and Lamb, 2017).

Preservation of knickpoint shape during retreat is poorly understood as very little data exist on the temporal evolution of their shape. For instance, knickpoints along the Atacama Fault System are systematically reduced in height compared to the height of ruptures directly on the fault scarp (Ewiak et al., 2015). On the contrary, 10 years after Chi-Chi earthquake, the height of co-seismic knickpoints ranged from 1 to 18 m (Yanites et al., 2010), while the initial surface rupture was limited to 0.5 to 8 m in height (Chen et al., 2001). Theoretically, only the stream power model with a linear dependency on slope predicts the preservation of knickpoint shape, favoured by a parallel retreat (e.g. Rosenbloom and Anderson, 1994; Tucker and Whipple, 2002; Royden and Perron, 2013). A less-than-linear dependency on slope leads to concave knickpoints, while a more-than-linear dependency on slope leads to convex knickpoints. Transport-limited models that reduce to advection–diffusion laws lead to a diffusion of the differential topography associated with knickpoints. However, transport-limited models are likely more pertinent to predict the evolution of fault scarps along hillslopes (e.g. Rosenbloom and Anderson, 1994; Arrowsmith et al., 1996, 1998; Tucker and Whipple, 2002), and evidence points toward a linear dependency on slope for knickpoint erosion (Lague, 2014). Yet, the transformation of fault activity and slip during earthquakes to knickpoints and hillslope scarps and their preservation throughout their subsequent erosion and retreat remains a challenging issue.

## 3.1 Fault setting

The tectonic setting considered here is the one of a typical active
intracontinental thrust fault, able to generate earthquakes up to magnitude
7.3. The thrust fault has a length *L*=200 km, a width *W*=30 km and a dip
angle $\mathit{\theta}=\mathrm{30}{}^{\circ}$ so that the fault tip is located at a 15 km
depth. The duration of the simulation *T* is set to 10 kyr to cover many
seismic cycles and for earthquakes to be well distributed along the finite
fault plane. A schematic sketch illustrates the model setup (Fig. 1).

## 3.2 Mainshocks

Mainshocks are generated along the fault plane. The potential magnitude
range of mainshocks is bounded by fault width, which sets the maximum
earthquake rupture width and by a minimum rupture width, here chosen as 500 m. Based on Leonard (2010), the modelled thrust fault allows magnitudes
ranging from ${M}_{{\mathrm{w}}_{min}}=\mathrm{3.7}$ to ${M}_{{\mathrm{w}}_{max}}=\mathrm{7.3}$. Inside these bounds, the
magnitude of each mainshock is determined by randomly sampling the
Gutenberg–Richter distribution, with a *b* value of 1 (Fig. 2). The earthquake
productivity of the distribution is inferred based on the arbitrarily chosen
rate of mainshock (*R*=0.1 d^{−1}), leading to $a={\mathrm{log}}_{\mathrm{10}}\left(R\phantom{\rule{0.125em}{0ex}}T\right)+b\phantom{\rule{0.125em}{0ex}}{M}_{{\mathrm{w}}_{min}}=\mathrm{8.975}$. The time occurrence of each mainshock is randomly
sampled over the duration of the simulation. Each mainshock is therefore
considered independent, and the only relationship between mainshocks is that
their population statistically respects the Gutenberg–Richter distribution
(Gutenberg and Richter, 1944).

The spatial location of mainshocks inside the fault plane is sampled using a
2-D distribution that corresponds to a truncated normal distribution
across-strike and to a uniform distribution along-strike (Fig. 3). A normal
distribution with depth roughly mimics the depth distribution of natural
earthquakes in the upper crust, which tends to show a maximum number of
earthquakes at intermediate depth and less towards the top and the tip of
the fault (e.g. Sibson, 1982; Scholz, 1998). Therefore, we set the mean of
the normal distribution equal to a 7.5 km depth as the fault tip has a 15 km
depth, so earthquakes are more numerous at this intermediate depth.
We define two end-member models, referred to as (1) the “seismic and
aseismic slip” model using a variance of the normal distribution $\mathit{\sigma}=W/\mathrm{10}$, corresponding to a narrow depth distribution, and (2) the “only
seismic slip” model with *σ*=3.3 *W*, corresponding to an almost uniform
depth distribution. We impose that the maximum earthquake frequency, at
depth 7.5 km, is equal in between all the models.

## 3.3 Aftershocks

Each mainshock triggers a series of aftershocks that is determined based on
the BASS model (Turcotte et al., 2007). It
represents an alternative to the more classical epidemic-type aftershock
sequence (ETAS) models (Ogata, 1988), with the advantage of being fully
self-similar. We here only briefly describe the BASS model, as more details
can be found in Turcotte et al. (2007). Based on a mainshock, the BASS model
produces a sequence of aftershocks which respect four statistical laws: (1) the Gutenberg–Richter frequency–magnitude distribution (Gutenberg and
Richter, 1944; Fig. 1); (2) a modified Båth's law (Shcherbakov and
Turcotte, 2004), which controls the difference in the magnitude of a
mainshock and its largest aftershock; (3) a generalized form of Omori's law
describing the temporal decay of the rate of aftershocks (Shcherbakov et
al., 2004); and (4) a spatial form of Omori's law that controls the
spatial distribution of aftershocks (Helmstetter and Sornette, 2003). The
BASS model relies on six parameters: the *b* value that we set equal to *b*=1,
the magnitude difference Δ*M*_{w}=1.25 of Båth's law, the
exponent *p*=1.25 and offset *c*=0.1 d of the temporal Omori law, and
the exponent *q*=1.35 and offset *d*=4.0 m of the spatial Omori law.
The values of these aftershock parameters are based on Turcotte et al. (2007) and are constant for all the simulations performed in this paper.
Seismicity along the fault is therefore made of mainshocks and their
aftershocks. This aftershock model is also similar to the one developed by
Croissant et al. (2019).

## 3.4 Earthquake rupture

The length *L*_{rup}, width *W*_{rup} and average co-seismic displacement
*D* of each earthquake rupture, including mainshocks and aftershocks, are
determined using scaling laws with seismic moment *M*_{O}, empirically
determined from a set of intraplate dip-slip earthquakes (Leonard, 2010)
following Eqs. (2)–(4):

where *C*_{1}=17.5, ${C}_{\mathrm{2}}=\mathrm{3.8}\times {\mathrm{10}}^{-\mathrm{5}}$ and $\mathit{\beta}=\mathrm{2}/\mathrm{3}$ are constants, and
*μ*=30 GPa is the shear modulus (Fig. 1). The locations of the rupture
patches around each earthquake are positioned randomly to prevent
hypocentres being centred inside their rupture patches. The fault has some
periodic boundary conditions, in the sense that if the rupture patch of an
earthquake exceeds one of the fault limits, the rupture area in excess is
continued on the opposite side of this limit. This choice maintains a
statistically homogeneous pattern of fault slip rate on the fault plane in
the case of the “only seismic slip” model (which displays an almost
homogeneous distribution of mainshocks on the fault plane). Another
strategy, consisting in relocating each rupture in excess inside the fault
limits, was dismissed as it was leading to gradients of fault slip rates
close to fault tips.

## 3.5 Seismic and aseismic slip

Slip along the fault plane is partitioned between seismic and aseismic slip.
The average slip rate *V*_{F} of the fault over the duration of the
simulation is given by ${V}_{\mathrm{F}}={V}_{\mathrm{S}}+\phantom{\rule{0.125em}{0ex}}{V}_{\mathrm{A}}$, where ${V}_{\mathrm{S}}=\sum {M}_{O}/\left(\mathit{\mu}\mathrm{TWL}\right)$ is the seismic slip,
due to all the earthquakes rupturing the fault, and *V*_{A} is aseismic slip. The average degree of seismic coupling on the fault plane is $\mathit{\chi}={V}_{\mathrm{S}}/{V}_{\mathrm{F}}$ (Scholz, 1998) and represents the
proportion of fault slip that occurs by earthquakes and seismic slip. We
define the reference fault slip rate as equal to the seismic slip rate of
the “only seismic slip” model so that ${V}_{\mathrm{F}}={V}_{\mathrm{S}}\simeq \mathrm{15}$ mm year^{−1}.
This velocity is only given approximatively, as the model developed here is
stochastic and leads to intrinsic variability in the number and magnitude of
earthquakes for the same parametrization. We follow the paradigm of
statistically homogeneous long-term fault slip over the fault. The “only
seismic slip” model, with an almost uniform spatial distribution of
mainshocks, is therefore on average fully coupled, with *χ*=1, while the
“seismic and aseismic slip” model, displaying a large change with depth of
the distribution of mainshocks, is dominated by aseismic slip with *χ*≃0.25. In the modelling framework developed here, even a fully coupled
fault can display significant spatial variations of fault slip rate. Slip
rate on the fault plane of the “only seismic slip” model varies between
11.4 and 18.2 mm year^{−1} for an average value of ∼15 mm year^{−1}. However, these spatial variations are randomly distributed
and do not follow any specific pattern (Fig. 4).

## 3.6 River uplift

A virtual river, orientated orthogonally to the fault trace, crosses the
fault trace at its centre, at $x=L/\mathrm{2}$. This river witnesses the distribution
of co-seismic and aseismic displacement modifying its topography and slope.
For the sake of simplicity, we assume that (1) any surface rupture generates
displacement only in the vertical direction and (2) that co-seismic and
aseismic deformation lead to a block uplift of the hanging wall, homogeneous
along the river profile. These two assumptions clearly neglect the influence
of the fault dipping angle and of the spatial distribution of uplift in
surface during an earthquake, which depends mainly on earthquake magnitude,
depth, geometry and the crustal rheology. In turn, earthquakes that do
not rupture the surface at the location of the river have no effect on river
topography and slope in this simple model. The rate of uplift is equal to
*V*_{F} at the intersection between the fault trace and the river.

## 3.7 River topographic evolution

To model river erosion, we consider a simple model considering that
knickpoints migrate upstream at constant velocity *V*_{R} along the
*y* axis, which is perpendicular to the fault trace orientated along the
*x* axis. We show in Appendix A that this constant migration velocity
model corresponds to a prediction of the stream power law (Howard and Kerby,
1983; Howard et al., 1994; Whipple and Tucker, 1999; Lague, 2014) which holds if
drainage area is about constant over the region of interest. Our model is
therefore appropriate to model knickpoint migration in near-fault conditions
and for large drainage areas. In the following, we only consider the
migration of slope patches over short distances upstream, during the
*T*=10 kyr of the simulation. We set the horizontal retreat rate to
*V*_{R}=0.1 m year^{−1}, which corresponds to a high rate of knickpoint
retreat over geological timescales (>10^{3} years) but a
moderate one over shorter timescales (e.g. Van Heijst and Postma, 2001).

## 3.8 Numerical implementation

Numerically, we solve in 2-D the evolution of a river profile crossing a
fault, subjected to slip during earthquakes and to aseismic slip. After
having set the parameters, the model (1) generates mainshocks and
aftershocks, including their magnitude, location and timing, and (2) computes
the time evolution of the river profile subjected to uplift and erosion.
Time stepping combines a regular time step, to account for uplift by
aseismic slip, with the time of occurrence of each earthquake rupturing the
surface at the location of the river, to account for co-seismic slip. During
each aseismic time step, one node of coordinates (*h*=0, *y*=0) is added to
the river profile at the downstream end of the river (i.e. the location of
the fault trace). During each co-seismic time step, two nodes of coordinates
(*h*=0, *y*=0) and (*h*=*D*, *y*=0) are added to the river at the downstream
end of the river, to represent the vertical step associated with the
co-seismic knickpoint. The remaining nodes, located upstream, are uplifted
following the aseismic uplift rate *V*_{A} and potential co-seismic
displacement. River erosion is accounted for by horizontal advection of
river nodes following a constant velocity *V*_{R} along the *y* axis. As we
neglect the contribution of horizontal displacement due to fault slip, we do
not consider any horizontal advection induced by tectonics, contrary to some
previous studies (Miller et al., 2007; Castelltort et al., 2012; Thieulot et
al., 2014; Goren et al., 2015).

## 4.1 Magnitude distributions

We first use this model to investigate the distribution of earthquake
magnitudes that rupture (1) the fault, (2) the surface and (3) the surface at
the location of the river (Fig. 4). For clarity, the frequency–magnitude
distributions are shown as incremental distributions *N*(m) and
not as cumulative distributions *N*(≥*m*). Unsurprisingly, the
frequency–magnitude distribution of earthquakes on the fault follows a
negative power-law distribution with an exponent of *b*=1, following the
imposed Gutenberg–Richter distribution. Increasing the degree of seismic
coupling *χ* only shifts the distribution vertically by increasing the
total number of earthquakes.

The distribution of earthquakes rupturing the surface follows a negative power law with an exponent of −0.5 for the “only seismic slip” model with a high degree of seismic coupling. In the case of the “seismic and aseismic slip” model, characterized by a lower degree of seismic coupling, the distribution follows a more complex pattern. Below a threshold magnitude, here around 6, the distribution follows a negative power law with an exponent of −0.5. Above this threshold magnitude, the distribution rises to reach the Gutenberg–Richter distribution and then decreases following the trend of the Gutenberg–Richter distribution. This results from the non-uniformity of the distribution of earthquakes with depth. In this model, large-magnitude earthquakes can rupture the surface, without requiring their hypocentres to be at shallow depth, whereas small-magnitude earthquakes can only rupture the surface if their hypocentres are located close to the surface, which is unlikely due to the shape of the depth distribution of mainshocks (Fig. 3). The threshold magnitude depends on the depth distribution of mainshocks, and particularly on its upper limit but also on the aftershock depth distribution that extends the range of possible depths due to Omori's law in space.

The distribution of earthquake magnitude rupturing the river follows a uniform distribution for the “only seismic slip” model. This novel result has potentially large implications as it means that a river has an equal probability of being ruptured by large or small earthquakes. This homogeneous distribution results from considering earthquake ruptures at one location and is yet consistent with a Gutenberg–Richter distribution of magnitudes along the modelled 2-D fault plane. However, for the “seismic and aseismic slip” model, mostly large-magnitude earthquakes manage to have ruptures cutting the river profile. Low-magnitude earthquakes, except for a few events, do not rupture the river. The magnitude threshold for river rupture is close to 6, similar to the one observed for surface ruptures. To date, there is no universal model of the depth distribution of earthquakes and of the partitioning between aseismic and seismic slip at shallow depth for intra- or inter-plate faults (e.g. Marone and Scholz, 1988; Scholz, 1998, Schmittbuhl et al., 2015; Jolivet et al., 2015). Yet, our results, i.e. a uniform distribution of earthquake magnitude cutting the river in the fully seismic case or only large-magnitude earthquakes rupturing the river for the model dominated by aseismic slip at shallow depth, clearly offer a guide to analyse river profiles in terms of fault properties.

## 4.2 Displacement distributions

Fault displacement *D* during an earthquake scales linearly with seismic
moment *M*_{O} (Wells and Coppersmith, 1994; Leonard, 2010), which is related
to magnitude by a logarithmic function, ${M}_{\mathrm{w}}=\mathrm{2}/\mathrm{3}{\mathrm{log}}_{\mathrm{10}}\left({M}_{O}\right)-\mathrm{6.07}$ (Kanamori, 1977). It results that a uniform distribution of
earthquake magnitude, that is observed for earthquakes cutting the river in
the case of the “only seismic slip” model, should lead to a negative
exponential distribution of earthquake displacements. The same finding exists with
the numerical model (Fig. 4e). In the case of the “seismic and aseismic
slip” model, it is more difficult to quantitatively characterize the
resulting distributions due to the lack of events, but we observe a
relatively uniform distribution of surface displacements.

## 4.3 Temporal distributions

We now investigate the time distribution of earthquakes rupturing the surface at the location of the river and their associated displacement. The “seismic and aseismic slip” and the “only seismic slip” models have 20 and 299 earthquakes cutting the river, respectively. Their average co-seismic displacement is 1 and 0.5 m, respectively. This illustrates that models dominated by aseismic slip have less frequent earthquakes cutting the river but that their average displacement is greater, due to the censoring of surface ruptures associated with low-magnitude earthquakes (Fig. 4).

Consistent with this last result, the inter-event time Δ*t* in between
successive earthquakes cutting the river increases significantly from the
“only seismic slip” model to the most “seismic and aseismic slip” model.
In other words, the frequency of surface rupture is higher in the most
seismic models and decreases with aseismic slip. This inter-event time
distribution follows for each model an exponential decay (Fig. 5), which is
consistent with a Poisson process. For the “seismic and aseismic slip”
model, the low number of events, 20 earthquakes, precludes characterizing a
negative exponential distribution. This exponential decay implies that fault
properties have no major effect on the temporal structure of earthquakes
cutting a river, only on their frequency.

## 5.1 Constant knickpoint velocity

If the slope patches generated by differential motion across the fault do
not migrate horizontally, due to, for instance, a lack of erosion, the
succession of earthquakes would progressively build a vertical fault scarp
in this model. Here, we rather consider the case of a migrating topography
due to river backward erosion following a kinematic model with *V*_{R}=0.1 m year^{−1}. It results in an averaged river slope just upstream the fault
trace of $\phantom{\rule{0.125em}{0ex}}\mathit{\phi}={V}_{\mathrm{F}}/{V}_{\mathrm{R}}=\mathrm{0.15}$ or 8.5^{∘}, with
*V*_{F}=15 mm year^{−1} (see Appendix A and Fig. A1). River profiles are
obtained for the two different models (Fig. 6). We first only consider
seismic slip, so that only earthquakes rupturing the river contribute to
topographic building. After *T*=10 kyr of model duration, the models have
resulted in about 20 to 150 m of topographic building for the “seismic and
aseismic slip” and “only seismic slip” models, respectively. The local
ratio between *V*_{S} and *V*_{F} can depart from their fault-averaged
values *χ*, due (1) to the non-homogeneous distribution of co-seismic
slip on the fault for models with significant aseismic slip and (2) to the
stochasticity of each model. For instance, the “seismic and aseismic slip”
model shows an apparent ratio of $\mathrm{20}/\mathrm{150}\simeq \mathrm{0.13}$ compared to its average
value of *χ*=0.25. Each successive co-seismic knickpoint is separated by a
flat river section, due to the absence of slope building by aseismic slip.
As expected, the “only seismic model” displays a larger number of
co-seismic knickpoint than the aseismic model. Adding aseismic slip leads to
sloped reaches between each knickpoint (Fig. 6), with slopes equal to *V*_{A}∕*V*_{F}. There is obviously a larger slope variability in the models
dominated by seismic slip due to a larger number of knickpoints.

## 5.2 Knickpoint velocity that depends on knickpoint height

Even if most simulations of this paper are done with a simple kinematic
model using a constant knickpoint velocity, we now consider a model with a
knickpoint velocity that depends on knickpoint height with ${V}_{\mathrm{R}}=r{\left(\mathrm{1}+\mathrm{\Delta}h/\mathrm{\Delta}{h}_{\mathrm{0}}\right)}^{q}$, where *d* is a constant, set to the
previously used constant knickpoint retreat rate of 0.1 m year^{−1}, Δ*h* the knickpoint height, Δ*h*_{0}=1 m a reference knickpoint height
and *r*=0.1 an exponent representing the sensitivity of knickpoint velocity
to knickpoint height. This model is motivated by mechanical arguments
suggesting a dependency of knickpoint velocity to their height (Scheingross
and Lamb, 2017). We allow a quicker knickpoint of height Δ*h*_{i}
that encounters slower knickpoints of height Δ*h*_{j} to merge,
forming in turn a single knickpoint of height Δ*h*_{i}+Δ*h*_{j} and of greater speed than the former knickpoints. The resulting
river profile can be compared to the one obtained with the “only seismic
slip” model (Fig. 7a). The dependency of knickpoint speed to height leads
to a river profile with high but fewer knickpoints. The interdistance
between successive knickpoints increases with total retreat. Small
knickpoints only survive close to the fault before being “eaten” by
quicker and higher knickpoints during their retreat. Only the highest
knickpoints, reaching tens of metres in height, survive after a significant
distance of retreat. This behaviour is also evidenced when comparing the
distributions of knickpoint heights for these two models (Fig. 7b). The
dependency of knickpoint velocity to height leads to very few knickpoints,
with however a large proportion of them having a metric or decametric scale.
This highlights that even limited non-linearities in the knickpoint retreat
model can lead to river profiles with significant differences.

## 5.3 Sediment cover, fault burial and knickpoint formation

We have neglected up to now the role of sediments and their impact on
knickpoint formation. More specifically, fault scarps can remain buried
during the aggradation phase of an alluvial fan located immediately
downstream of the fault (e.g. Carretier and Lucazeau, 2005). This mechanism
is suggested to be a primary control of knickpoint and waterfall formation
by allowing the merging of several small co-seismic scarps formed during
burial phases into single high-elevation waterfalls that migrate during
latter incision phases (Finnegan and Balco, 2013; Malatesta and Lamb, 2018).
We test this mechanism and its impact on river profiles using a simple
description of fault burial by sediment cover (Fig. 8). At each time step,
the formation of a knickpoint can only occur if the fault scarp height,
*h*(*y*=0), is greater than the sediment thickness of the alluvial fan, *h*_{s}.
In this case, the formed knickpoint height is simply $h\left(y=\mathrm{0}\right)-{h}_{\mathrm{s}}$.

Temporal variations of sediment thickness are prescribed using four scenarios:
(1) no sediment cover, *h*_{s} =0, corresponding to the reference model (Fig. 8e); (2) a square-wave (or step-like) function with a periodicity of 2000 years
and a maximum amplitude of 10 m (Fig. 8d); (3) a sinusoidal function with a
periodicity of 2000 years and a maximum amplitude of 10 m (Fig. 8c); and (4) an
earthquake-driven sediment cover, where sediment increases instantaneously
after each earthquake that ruptures the river with an amplitude arbitrarily
defined proportional to (*M*_{w}−5)^{2}, followed by a linear decrease over
100 years, following results by Croissant et al. (2017) (Fig. 8b). This last
scenario mimics, in a very simplified manner, the potential transient
response of an alluvial fan to the observed increase of river sediment load
induced by earthquake-triggered landslides (Hovius et al., 2011; Howarth et
al., 2012; Croissant et al., 2017). Alternatively, the periodic scenario
mimics the potential response of sediment thickness to some climatic cycles.
These scenarios are purely illustrative and do not aim at offering an
accurate description of the impact of tectonic or climatic changes on
sediment cover dynamics. For each scenario, except the one with no sediment
cover, the mean sediment thickness is 5 m. For the sake of simplicity, we
only consider the “only seismic slip” model with the same temporal sequence of
earthquakes in each of the four scenarios. Knickpoint velocity is kept constant
and equal to *V*_{R}=0.1 m year^{−1}.

The square-wave model is useful to assess the impact of abrupt changes in sediment thickness. During the phase of a high sediment cover thickness that lasts 2000 years, the scarp progressively builds its height until it reaches 10 m during successive fault ruptures. Over this period, there is no knickpoint formation, while previously formed knickpoints continue to migrate upstream, leading to elongated flat river reaches upstream of the fault. Once the scarp is re-exposed, the following earthquakes generate knickpoints (yellow dots in Fig. 8a), with their individual height corresponding to each associated earthquake displacement. Then, the abrupt transition from 10 m of sediment thickness to no sediment thickness suddenly exposes 10 m or more of fault scarp that forms a migrating knickpoint of elevation much higher than the largest earthquake displacement, i.e. 1.8 m. Then, during the 2000 years that follow, with no sediment cover, each earthquake rupture generates a new knickpoint (blue dots in Fig. 8), as in the reference model.

The sinusoidal model, mimicking climatic oscillations (Fig. 8c), displays a relatively similar behaviour, except that it does not form 10 m high knickpoints during the phase of degradation of the sediment cover. Instead, this phase leads to the formation of “climatic knickpoints” as the rate of decrease in sediment thickness is greater than the rate of scarp building by fault slip. For the exact same reason, the phase of sediment aggradation is characterized by no knickpoint formation and by flat river reaches. Knickpoint formation and the signature of the river profile are therefore dominated by the climatic signal controlling sediment aggradation–degradation phases rather than by fault slip.

In these scenarios, the fault-burial mechanism by sediment cover does not necessarily lead to knickpoints with elevation greater than earthquake ruptures, except for abrupt removals of sediment cover such as in the square-wave model. Yet, in all these models, the fault-burial mechanism limits the periods of differential topography building, leading in turn to succession of steepened river reaches or knickzones, corresponding to periods of sediment removal, alternating with low slope river reaches, corresponding to periods of sediment aggradation. Figure 9 illustrates the role of sediment cover in modulating the surface expression of tectonics and co-seismic displacement. For the highest rates of sediment aggradation and removal, river profiles are dominated by the temporal evolution of the sediment cover and not by the activity of the fault, whereas for limited sediment aggradation and removal rates, the river profiles and the succession of knickpoints are dominated by the temporal occurrence of earthquakes and not by the temporal evolution of the sediment cover. These results are consistent with the ideas developed by Malatesta and Lamb (2018).

## 6.1 From single to several parallel rivers

We now explore the degree of spatial correlation in between the topographic profiles of several parallel rivers flowing in an across-strike manner along the fault trace. For the sake of simplicity, we ignore the role of sediment cover on knickpoint formation and use a constant knickpoint velocity. Paleo-seismological studies using knickpoints to infer fault activity generally consider the distributions of knickpoints along several subparallel rivers to lead to statistically robust analyses and to assess the spatial extent of each past earthquake (e.g. Ewiak et al., 2015; Wei et al., 2015; Sun et al., 2016). Correlating topography and knickpoints along the strike of a fault, using parallel rivers, also offers independent means to assess the rupture length and the magnitude of a past earthquake. Using multiple rivers is also less likely to be biased by potential heterogeneities occurring along single rivers.

We therefore consider a set of rivers separated by Δ*x*=1 km along the
strike of the fault, i.e. the *x* direction. Because (1) the drainage area of
each of these rivers can vary by orders of magnitude and (2) because
knickpoint retreat rates show a high variability, their knickpoint migration
rate *V*_{R} is randomly sampled in the range 0.001 to 0.1 m year^{−1}. Each profile of the 200 rivers shares some common topographic
characteristic, including their average number of knickpoints and total
elevation (Fig. 10). However, their average slopes and the horizontal
position *y* of the knickpoints largely differ due to the variability of
*V*_{R}. Knowing a priori *V*_{R} and the duration *T* of the simulation (i.e. the
age of the knickpoints) enables to define a normalized horizontal position
relative to the fault, *y*∕(*T**V*_{R}). Practically, several studies normalized
distance by the square root of drainage area, as drainage area is generally
used as a proxy for retreat rate (e.g. Crosby and Whipple, 2006).
Knickpoints generated at the same time, along different rivers with
different retreat rates, share the same normalized distance relative to the
fault. This representation is convenient to assess the spatial extent of an
earthquake rupturing several rivers along-strike. Non-normalized river
profiles are shown in Fig. B1.

## 6.2 Knickpoint correlation in between several parallel rivers crossing the fault

This representation is convenient to assess the degree of correlation of the
profiles of the successive rivers. Obviously, there is no significant
topographic correlation when considering rivers with such a high variability
in retreat rates, e.g. 0.001 to 0.1 m year^{−1}. We therefore compute
the matrix of correlation between each river elevation profile using the
river normalized horizontal distance (Fig. B1). River elevation is corrected
or “detrended” from its average slope to remove an obvious source of
topographic correlation. We then compute the average coefficient of
correlation for a given river interdistance Δ*x* ranging from 0 to
100 km (Fig. 11). The two models, the “only seismic slip” and the
“seismic and aseismic slip” models, show a similar pattern, with a
significant positive correlation (>0.5) for rivers separated by
less than 14 to 23 km (10 to 45 km if accounting for the standard
deviation). The maximum distance over which a correlation is significant
corresponds to about 35 km, half the maximum co-seismic rupture length of
∼70 km along the considered fault. This illustrates that
knickpoints should not be correlated for rivers separated by more than this
distance, considering the tectonic setting of this model, and fault
dimensions. This correlation distance could increase using a wider fault
generating larger-magnitude earthquakes with longer surface rupture. We also
find that the correlation is better for the model dominated by aseismic slip
and showing less knickpoints (Fig. B1). Positive correlations were obtained
using horizontal distance normalized by retreat rate. However, using only
catchments with similar retreat rates would also lead to positive and
significant correlation even when using non-normalized distance.

## 7.1 Knickpoint detectability for the reference model

River profiles are used in many studies to extract co-seismic knickpoints and to assess fault activity and local-to-regional seismic hazard (e.g. Ewiak et al., 2015; Wei et al., 2015; Sun et al., 2016). It is therefore required to investigate whether modelled knickpoints are statistically detectable. Knickpoint detection often relies on the use of digital elevation models and topographic data (e.g. Neely et al., 2017; Gailleton et al., 2019), which are obtained at a certain scale or resolution. The detectability of each individual knickpoint depends not only on its distance to its adjacent knickpoints but also on the horizontal resolution and vertical precision of the topographic data and on the roughness of the riverbed. In the following, we consider that a knickpoint is detectable if its height is greater than the vertical precision of topographic data and if its distance to adjacent knickpoints is greater than the horizontal resolution of topographic data.

Resolutions of topographic data available at the global scale (e.g. SRTM or ASTER) are between 10 and 100 m, with precision not better than a few metres. Local-to-regional topographic datasets obtained from current airborne lidar or photogrammetric data or derived from aerial or satellite imagery (e.g. Pléiades) display a resolution between 0.5 and about 1–5 m and a typical vertical precision of 10 cm above water. Moreover, in the vertical direction, knickpoint detectability depends also on the inherent bed roughness, mean alluvial deposit thickness and the local distribution of sediment grain size. Sediment grains of dimension greater than 0.1 m are commonly found in rivers located in mountain ranges (e.g. Attal and Lavé, 2006), especially at low drainage areas, and there is often a thin layer of sediment covering the channel bed, potentially hiding bedrock features. If we fully acknowledge the role of river roughness, we focus here on the issue of detectability relative to topographic resolution and precision, for the sake of simplicity, and using knickpoints formed by the “only seismic slip” model.

In terms of vertical precision, a precision of 0.1 m (e.g. lidar) enables the detection of knickpoints produced by an earthquake as low as magnitude ∼4.8 (Fig. 12a). For rivers permanently under water, traditional airborne lidar using near-infrared laser or photogrammetric data cannot measure river bathymetry imposing a detectability level and an uncertainty of knickpoint height of the order of the water depth. Topographic data with a precision of about 1 m would only enable to detect knickpoints for earthquakes of magnitude above 6.8. It results that about 18 % of the knickpoints are detectable using 1 m precision, while 72 % are detectable with lidar data and a precision of 0.1 m (Fig. 12b). SRTM or ASTER data have precisions of a few metres, at best, that would only enable the potential detection of earthquakes of magnitude ∼8 or more.

In terms of horizontal resolution, we assess knickpoint “detectability” by comparing knickpoint interdistance with the resolution of topographic data. The distribution of horizontal distance between successive knickpoints (that scales with the distribution of inter-event times) shows that knickpoint interdistance ranges between less than a millimetre to up to few tens of metres (Fig. 12a). Using a resolution of 10 m, only 7 % of the knickpoints can be detected, while a resolution of 1 m increases this percentage to 65 %. Combining horizontal and vertical detectability reduces even more the detectability of knickpoints, as only 2 or 45 % of the knickpoints are detectable using lidar (1 m resolution, 0.1 m precision) or digital elevation model (DEM) (10 m resolution, 1 m precision) characteristics, respectively.

## 7.2 Knickpoint detectability along rivers, colluvial channels and hillslopes

We now consider the issue of knickpoint detectability for a broader range of
model parameters, in particular fault dimensions, fault slip rate and
knickpoint retreat rate. Vertical detectability depends only on the range of
considered earthquake magnitude and displacement. As the maximum modelled
magnitude is directly limited by the dimension of the modelled fault,
considering greater knickpoints requires extending the dimension of the fault
and more specifically its width. In our model, horizontal detectability is
directly dependent on the river slope, $\mathit{\phi}={V}_{\mathrm{F}}/{V}_{\mathrm{R}}=\mathrm{0.15}$.
Indeed, the horizontal distance between successive knickpoints increases
linearly with knickpoint migration velocity *V*_{R}, while its decreases
linearly with the rate of fault slip *V*_{F} that sets the rate of
earthquake and knickpoint formation. Here, we investigate how knickpoint
detectability varies with slope *φ* (Fig. 12b). We consider that
river channels have a slope below 0.2, colluvial channels between 0.1 and
0.5, and hillslopes above 0.2, following classical slope–area relationships.
In terms of horizontal detectability, rivers with a slope $\mathit{\phi}<{\mathrm{10}}^{-\mathrm{2}}$ have a good detectability, more than 80 % using lidar or even
DEM resolution. Rivers with ${\mathrm{10}}^{-\mathrm{2}}<\mathit{\phi}<{\mathrm{2.10}}^{-\mathrm{1}}$ have a good
horizontal detectability using lidar data and a moderate one using DEMs (10 %
to 80 %). Colluvial channels have a moderate horizontal detectability
using lidar data and a poor one using DEMs (<10 %). Hillslopes
with *φ*>10^{0} have a poor detectability even with lidar data.
However, the overall detectability of DEM data is below 20 % due to the
issue of the vertical detectability; that is low for DEM data. The overall
detectability of lidar data can reach 70 % for low values of *φ*,
in the river domain, while it is moderate or even poor for colluvial channels
or hillslopes, respectively. This highlights the need for lidar data to
detect in a more systematic manner knickpoints along river or colluvial
channels.

## 8.1 Model limitations

To approach the problem of co-seismic knickpoint formation and their impact on river profile, we have made several simplifying assumptions. The spatial and temporal distribution of earthquakes, including mainshocks and aftershocks, only follow classical statistical and scaling laws. Fault stress state or friction properties, which are first-order controlling factors of earthquake triggering (e.g. Scholz, 1998), are not explicitly accounted for. Earthquake ruptures are assumed to be rectangular, to have dimensions scaling with seismic moment and to display a homogeneous displacement (Leonard, 2010), while natural ruptures display more variable behaviours. The relative contributions of seismic and aseismic processes to fault slip, and their spatial distributions are defined in a relatively ad hoc manner. Moreover, co-seismic displacement follows a block uplift mechanism, which contradicts observations and neglects the elasticity of the lithosphere. Yet, it is to be emphasized that block uplift in near-fault conditions for large-magnitude earthquakes corresponds to an asymptotic behaviour. A more realistic approach is to compute the surface displacement induced by each earthquake using, for instance, dislocations embedded into an elastic half-space (e.g. Okada, 1985). This alternative approach would also have the benefit of accounting for the surface displacement of earthquakes that do not rupture the surface. Moreover, surface rupture only occurs along a single fault and does not account for off-fault damage (e.g. Zinke et al., 2014), that could also generate knickpoints, or for more complex rupture geometry (e.g. Romanet et al., 2018). Knickpoint retreat along the river profile was modelled using a constant velocity, which corresponds to an asymptotic behaviour of the stream power incision model for small migration distance respective to the square root of river drainage area. If the migration of knickpoints or slope patches is classically modelled using the stream power incision model (Rosenbloom and Anderson, 1994; Whittaker and Boulton, 2012; Royden and Perron, 2013), this approach was recently questioned by experimental (Baynes et al., 2018) and field (Brocard et al., 2016) results, suggesting no obvious dependency of the migration rate on river discharge. Mechanistic models of waterfall erosion and retreat offer another more accurate but more complex approach (Scheingross and Lamb, 2017).

## 8.2 Model and results applicability to normal and strike-slip faults

The developed model, that was applied in this study to a continental thrust
fault, can also be directly applied to a normal fault. Indeed, the adopted
scaling relationships between earthquake rupture dimensions or displacements
and seismic moment (Leonard, 2010) apply to dip-slip earthquakes and
therefore to both normal and thrust faults. The main differences are the
polarity of motion between hanging walls and footwalls, and the dipping angle of
the fault. The latter difference vanishes in the developed approach as we
assume that rupture displacement occurs only in the vertical direction.
Under these limitations and simplifications, all the obtained results in
this paper can be therefore directly transcribed to normal faults. Because
normal faults tend to have a larger dipping angle, close to 60^{∘} on
average, than thrust faults, the approximation of purely vertical co-seismic
displacement is less incorrect for normal faults. Moreover, strike-slip
faults or dip-slip faults can also be accounted for by this model, by simply
tuning the parameters of the rupture scaling laws, i.e. *C*_{1}, *C*_{2} and
*β*. Assuming the depth distribution of seismicity along different
types of faults is identical (a likely incorrect hypothesis), changing the
type of fault would not have a major impact on the results presented in this
paper.

## 8.3 Knickpoints and horizontal tectonic displacement

Surface ruptures and displacements were only considered in the vertical
direction, clearly simplifying the variability in the orientation of natural
surface ruptures. If this paper is focused on the vertical expression of
faults along river profiles, future work should account for the influence of
horizontal tectonic displacement on river profiles (e.g. Miller et al.,
2007). Accounting for the dip angle of the fault and the associated
horizontal tectonic displacement can have two main effects: (1) move
knickpoints in the direction of tectonic motion and increase or decrease the
apparent retreat rate of knickpoints, in the case of normal or thrust
faults, respectively; and (2) move the position of the fault trace through
time, as, for instance, in the case of a thrust sheet when the hanging wall
moves over the footwall. We ignore the latter effect and focus on the
influence of tectonic motion on knickpoint retreat rate and on river slope.
Accounting for the dipping angle of the fault changes the expression of
river slope just upstream of the fault. Indeed, fault slip builds topography
in the vertical direction at a rate *V*_{F}sin(*θ*), while
knickpoints retreat by the cumulative effect of erosion and horizontal
tectonic displacement at a rate *V*_{R}±*V*_{F}cos(*θ*). The
sign ± is positive for normal faults but negative for thrust faults, as
knickpoints are displaced by tectonics towards the fault trace for the
latter. It results that the river slope becomes $\mathit{\phi}={V}_{\mathrm{F}}\mathrm{sin}\left(\mathit{\theta}\right)/({V}_{\mathrm{R}}+{V}_{\mathrm{F}}\mathrm{cos}(\mathit{\theta}\left)\right)$ for normal faults and $\mathit{\phi}={V}_{\mathrm{F}}\mathrm{sin}\left(\mathit{\theta}\right)/({V}_{\mathrm{R}}+{V}_{\mathrm{F}}\mathrm{cos}(\mathit{\theta}\left)\right)$ for thrust faults. For rivers which
generally have slopes lower than about 0.1, the horizontal tectonic
displacement *V*_{F}cos(*θ*) is likely to be
negligible compared to the retreat rate by erosion *V*_{R}, and the slope
can be approximated by $\mathit{\phi}\approx {V}_{\mathrm{F}}\mathrm{sin}\left(\mathit{\theta}\right)/{V}_{\mathrm{R}}$. However,
this approximation does not hold anymore for colluvial channels or for
hillslopes as the slope becomes closer to 1. Accounting for tectonic
displacement obviously changes the threshold of vertical detectability of
knickpoints as their height decreases when decreasing the fault dip angle.

## 8.4 Do mainshocks or aftershocks matter for knickpoints and river profiles?

Aftershocks play a secondary role in the seismicity model considered in this
paper. Indeed, for the “only seismic slip” models, aftershocks only
represent 18 % of the 442 188 earthquakes simulated on the fault.
Seismicity is therefore dominated by mainshocks. This is not surprising, as
about 95 % of earthquakes, that follow the Gutenberg–Richter
frequency–magnitude distribution, have a magnitude lower than 5 and have in
turn a very low probability to generate aftershocks because (1) the
aftershock model uses Båth's law with a magnitude difference between any
mainshock and their aftershocks, Δ*M*_{w}=1.25, and (2) the
minimum magnitude modelled is 3.7. Therefore, aftershocks will only be
triggered after intermediate- to large-magnitude earthquakes, *M*_{w}>5, which
only represent 5 % of the total number of earthquakes. It also results that
river profiles in our models are mostly built by mainshocks and not by
aftershocks that only represent 18 % of the cumulated uplift. Therefore,
developing an aftershock model to include earthquakes and their effects on
landscape evolution models represents an additional step in terms of model
complexity that is not mandatory. This means that simply accounting for
mainshocks by (1) sampling the Gutenberg–Richter distribution to determine
earthquake magnitude and (2) randomly sampling their location already
represents a consistent modelling approach towards including earthquakes in
landscape evolution models. Despite that, aftershocks can have significant
effects, punctually in time and space, for knickpoint formation, river
uplift or even landslide triggering (e.g. Croissant et al., 2019) that
justify for some studies the additional complexity of modelling them.

## 8.5 Knickpoint height distribution as a paleoseismological tool?

Co-seismic knickpoints are common geomorphological markers found in seismic areas (Boulton and Whittaker, 2009; Yanites et al., 2010; Cook et al., 2013). Several studies have offered constraints on fault seismogenic activity from the study of river profiles and knickpoint height (Boulton and Whittaker, 2009; Ewiak et al., 2015; Wei et al., 2015; He and Ma, 2015; Sun et al., 2016). Natural distributions of knickpoint height are systematically dominated by large heights, corresponding to earthquake magnitudes greater than 5. For instance, the magnitude of earthquakes deduced from knickpoints extracted along rivers crossing the Atacama Fault System follows a bell-shaped distribution favouring large-magnitude (5.8–6.9) earthquakes (Ewiak et al., 2015). Because the distributions of knickpoints were found to share similarities with the distribution of ruptures directly along the fault scarp, this rules out the hypothesis of fully eroded co-seismic knickpoints generated by small-magnitude earthquakes (Ewiak et al., 2015). This observation, of knickpoints dominated by large earthquakes and the censoring of small-magnitude earthquakes, is similar to the results obtained in this paper with the model dominated by aseismic slip at shallow depth (Fig. 3e). Alternative explanations for the apparent lack of small knickpoints or scarp ruptures in most natural datasets (Ewiak et al., 2015; Wei et al., 2015; He and Ma, 2015; Sun et al., 2016) include at least (1) the difficulty to detect the limited displacement induced by earthquakes of magnitude 5 or less and (2) the fault-burial mechanism (Finnegan and Balco, 2013; Malatesta and Lamb, 2018) that filters out small co-seismic surface ruptures. In any case, the depth distribution of earthquakes and of their rupture extent exerts fundamental control on the resulting height distribution of co-seismic knickpoint.

In turn, our results suggest that knickpoint datasets, that will become more and more accessible thanks to high-resolution topographic data, can be used to assess fault activity. Obviously, the height of knickpoints provides some form of evidence for the earthquakes that have generated them. A negative exponential distribution of knickpoint height points toward a purely seismic fault, while deviations from this trend can suggest aseismic slip or even a slip deficit at shallow depths. The main limitation is yet the poorly known impact of geomorphological processes on evolution of the shape of knickpoints. Some knickpoints along the Atacama Fault System have a reduced height compare to their initial rupture (Ewiak et al., 2015), while some knickpoints produced during Chi-Chi earthquake in 1999 were higher 10 years later (Yanites et al., 2010). These contrasting cases illustrate some potential, and poorly understood, pitfalls in using knickpoints to infer fault and seismic activity.

## 8.6 River dynamics: constant uplift or time-variable uplift with earthquakes?

Most numerical efforts attempting at modelling the long-term (>10–100 kyr) topographic building of mountainous or rift settings have used a
constant or smoothly varying uplift rate (e.g. Braun and Willett, 2013;
Thieulot et al., 2014; Campforts et al., 2017), not including the
variability of uplift rate during the seismic cycle. If using a stream power
incision model with a linear dependency on slope *n*=1, this choice is
acceptable as the variability of uplift rate and the associated variability
of slope patches shaped throughout the seismic cycle can be averaged out.
Moreover, knickpoint retreat rate is in this case independent of slope, as
this model corresponds to a linear kinematic wave equation (Rosenbloom and
Anderson, 1994; Tucker and Whipple, 2002; Whittaker and Boulton, 2012;
Royden and Perron, 2013). However, if using a non-linear dependency of
erosion rates to slope, with *n*≠1, and only considering a long-term
averaged uplift rate, and not its variability, is an approximation that
becomes more incorrect with the degree of non-linearity of the model. In
other words, the erosion rate of a river profile made of co-seismic
knickpoints separated by low-slope river sections built during aseismic
periods is not equivalent to the erosion rate of a smooth river profile with
the same average slope and built under a constant uplift rate. In a
non-linear stream power incision model, the retreat rate is sensitive to
slope at a power *n*−1. For *n*>1, greater slope patches will migrate
quicker than lower slope patches, and vice versa for *n*<1. While a large
proportion of the literature considers the linear stream power incision
model (or the unit stream power model) as the reference model, the
parametrization of the stream power incision and in particular of the slope
exponent *n* is still an open debate, as is its actual applicability to
model knickpoint migration (e.g. Lague, 2014). Moreover, the physics of
knickpoint or waterfall retreat likely depends on other variables such as
knickpoint height (Holland and Pickup, 1976; Hayakawa and Matsukura, 2003;
Haviv et al., 2010; Scheingross and Lamb, 2017), sediment supply (Jansen et
al., 2011), lithological structure (Lamb and Dietrich, 2009) and
lithological strength (Baynes et al., 2018). Even if this debate is clearly
out of the scope of this paper, the implication of this study for the
understanding of river erosion and dynamics should not be ignored. Indeed,
we have shown that even a slight sensitivity of knickpoint retreat rate to
knickpoint height leads to large differences in terms of river profile or
knickpoint height distribution, by rapidly merging out all small knickpoints
into larger ones associated with greater retreat rates (Fig. 7). Moreover, the
modelling results of this study show that the frequency–magnitude
distribution of earthquakes rupturing a river is uniform for purely seismic
faults and follows a bell shape, favouring large-magnitude earthquakes, for
faults with significant shallow aseismic slip. This result offers a
complementary – not an alternative – explanation to the fault-burial
mechanism (Malatesta and Lamb, 2018) for the apparent larger proportion of
high waterfalls.

## 8.7 Co-seismic displacements and knickpoints inside landscape evolution models?

Further implications on the impact of considering earthquakes in landscape dynamics can only be cast by using landscape evolution models (LEMs) (Croissant et al., 2017; Davy et al., 2017; Braun and Willett, 2013; Campforts et al., 2017; Egholm et al., 2011). The developed model in this paper can be implemented in most LEMs to investigate river and landscape response to earthquakes and their successions. However, the main foreseen difficulty is the large variability of inter-event times, which put strong constraints on the time stepping strategy. To overcome this difficulty, a minimum earthquake magnitude can be defined as a threshold: earthquakes with lower magnitudes are modelled as continuous fault slip, while earthquakes with greater magnitudes are modelled as discrete uplift events during a specific time step. A second difficulty is the spatial discretization of knickpoints that migrate inside the model domain. Most current LEMs use regular grids to discretize surface topography with a uniform resolution. To be consistent with their boundary conditions, such numerical schemes must adapt their spatial resolution to the typical modelled distance between successive knickpoints that can easily go below 1 m (Fig. 7b). This is problematic, as the efficiency of most LEMs scales at best with the number of model nodes (e.g. Braun and Willett, 2013). Using too-coarse resolutions would smooth out knickpoints and slope variability, leading to similar landscape evolution and dynamics as using a constant uplift, even with non-linear slope dependency. Another more adapted strategy is to use irregular grids, for instance, based on Delaunay triangulation, to discretize topography in LEMs (e.g. Braun and Sambridge, 1997; Steer et al., 2011). Despite being less commonly used in LEMs, irregular grids enable to properly account for co-seismic knickpoints and variable uplift rates by using fine resolutions close to knickpoints and coarser ones in other model domains. This in turn would lead to tractable model durations. Another benefit of irregular grids is their ability to be deformed in the horizontal directions. This is required to account for the horizontal components of co- or inter-seismic displacement that is systematically ignored in LEMs while being of greater amplitude than vertical displacement in convergent or strike-slip settings (e.g. Cattin and Avouac, 2000). Coupling inter- and co-seismic displacement with LEMs represents a future direction to further investigate the impact of earthquakes and tectonic deformation during the seismic cycle on landscape dynamics. The main remaining limitation is the development of mechanistic models for knickpoint retreat and evolution, a subject that has received recent attention (e.g. Scheingross and Lamb, 2017).

The accurate modelling of landscape evolution requires accounting for the temporal and spatial variability of surface uplift and displacement. We propose a statistical model of earthquakes, based on the BASS model (Turcotte et al., 2007), to simulate the slope and height distributions generated by earthquakes and aseismic slip at the intersection between a thrust fault and a river. The rupture extent and displacement of each earthquake are inferred using classical scaling laws (Leonard, 2010) that can be applied to strike-slip, normal or thrust faults. Slip along the fault plane is partitioned between seismic and aseismic slip using an ad hoc spatial distribution of mainshocks along the fault plane. Co-seismic uplift events, with rupture cutting rivers, generate knickpoints that migrate along the river profile following a constant retreat rate.

First, the developed model produces a uniform distribution of earthquake magnitude cutting the river that is obtained while imposing a Gutenberg–Richter frequency–magnitude distribution of earthquakes along the fault plane. In turn, the produced knickpoint heights follow a negative exponential height distribution. The interevent time distribution between successive knickpoints follows an exponential decay.

Second, partitioning shallow slip between seismic and aseismic slip censors the magnitude range of earthquakes rupturing the surface and cutting the river towards large magnitudes. Poorly coupled faults, dominated by shallow aseismic slip, generate mostly rare and on average high knickpoints, while fully coupled faults generate frequent knickpoints of moderate height, on average. Assuming no impact of geomorphological processes on the evolution of the shape of knickpoints, an unlikely hypothesis, these differences in the height distribution of knickpoints offer a guide to assess fault coupling and the shallow partitioning of fault slip over longer timescales than modern seismology.

Third, our results demonstrate the influence of earthquakes and of fault properties on river profiles. Using a constant knickpoint retreat rate, our simple model produces river profiles made of a succession of flat sections and knickpoints for fully coupled faults and straight river profiles with a constant slope and few knickpoints for poorly coupled faults. Accounting for a dependency of knickpoint retreat rate to knickpoint height leads to the progressive merging of small knickpoints into larger ones, with a height significantly greater than the vertical offset produced by the largest magnitude earthquakes. Moreover, fault burial by intermittent sediment cover can alter the surface expression of fault slip and earthquake activity, when the rate of sediment aggradation/degradation is greater than the rate of fault slip.

Fourth, knickpoint detectability, regarding the horizontal resolution and vertical precision of modern topographic datasets such as lidar or DEMs, directly depends on the river slope that is equal to the ratio between fault slip rate and knickpoint retreat rate. Decreasing the slope increases the horizontal distance between successive knickpoints and enhances knickpoint detectability. On the contrary, the vertical detectability is only limited to the vertical precision of topographic data relatively to the topographic offsets produced during earthquakes.

Fifth, when considering several parallel rivers distributed along the strike of the fault, a positive correlation between river profiles is obtained if the rivers are separated by less than half of the maximum rupture length occurring on the fault. This correlation is obtained using a horizontal distance normalized by knickpoint migration rates, or when considering rivers with similar migration rates. The coefficient of correlation becomes significantly positive (>0.5) when the river interdistance is less than about a quarter of the maximum rupture length. For a maximum earthquake magnitude of 7.3, this interdistance corresponds to 14 to 23 km and does not vary significantly with fault coupling.

Last, the developed model offers insights on the building of slopes and knickpoints by fault activity and earthquakes. This model could also be implemented in landscape evolution models to better infer the role of tectonics and earthquakes on landscape dynamics. This is pivotal to understanding how and why earthquakes build or destroy topography (Parker et al., 2011; Marc et al., 2016), investigating the feedbacks of erosion on fault dynamics over a seismic cycle (Vernant et al., 2013; Steer et al., 2014) or during orogenesis (Willet et al., 1999; Thieulot et al., 2014), isolating the feedbacks between river and hillslope dynamics (Valla et al., 2010; Jansen et al., 2011) or unravelling the source-to-sink relationships in seismically active landscapes (Howarth et al., 2012).

A simple MATLAB version of the model can be accessed through a GitHub and/or a Zenodo repository: https://github.com/philippesteer/RiverFault (last access: 3 May 2019) and https://doi.org/10.5281/zenodo.2654819 (Steer and Croissant, 2019).

A classical detachment-limited approach to describe the rate of river
erosion *E* is the stream power incision model (Howard and Kerby, 1983;
Howard et al., 1994; Whipple and Tucker, 1999; Lague, 2014) described in Eq. (5):

where *h* is the elevation of the bedrock bed of the river, *t* is the time,
*y* is the distance along the river (i.e. in an across-strike manner along the fault trace), $S=\mathrm{d}h/\mathrm{d}y$ is the local river slope, *A* is the upstream
drainage area, *K* is the erodibility, and *m* and *n* are two exponents.
Considering a linear dependency of erosion rates to slope, with *n*=1, the
stream power incision model is equivalent to a linear kinematic wave
equation. Under this condition, it can be demonstrated (Rosenbloom and
Anderson, 1994; Tucker and Whipple, 2002; Whittaker and Boulton, 2012;
Royden and Perron, 2013) that knickpoints or slope patches along the river
migrate upstream at a rate determined by Eq. (6):

Moreover, recent empirical results suggest that using *n*=1 and *m*=0.5 is
suited to describe knickpoint migration (Lague, 2014). If the total
migration distance is small compared to the entire river length, from its
source to the modelled frontal thrust fault, the migration velocity *V*_{R}
can be approximated as a constant. This condition holds only if *K**T*≪1,
considering that river length generally scales with about the square root of
drainage area (Hack, 1957). The horizontal knickpoint retreat rate
(*V*_{R}=0.1 m year^{−1}) can therefore be obtained for an infinite number of couples
of the *A* and *K* parameters, following the relationship *V*_{R}=*K**A*^{m}
that yet must at least satisfy the condition *K**T*≪1 (Fig. A1). Other
conditions exist, including the domain of validity in the *A* space of the
stream power incision model or that the slope generated for a given value of
*A* makes sense in terms of river steepness. However, they are not further
considered, as the scope of this paper is to develop a general quantitative
framework to investigate slope and topographic building by a fault.

Variable |
Definition |

M_{0} |
Earthquake moment magnitude |

M_{w} |
Earthquake magnitude |

N(≥M_{w}) |
Cumulative number of earthquakes of magnitude greater than M_{w} |

N(M_{w}) |
Incremental number of earthquakes of magnitude M_{w} |

a |
Gutenberg–Richter earthquake productivity |

b |
Gutenberg–Richter b value |

R |
Rate of mainshock |

ΔM_{w} |
Magnitude difference of Båth's law |

p |
Exponent of the temporal Omori law |

c |
Offset of the temporal Omori law |

q |
Exponent of the spatial Omori law |

d |
Offset of the spatial Omori law |

L_{rup} |
Earthquake rupture length |

W_{rup} |
Earthquake rupture width |

C_{1} |
Earthquake scaling law constant |

C_{2} |
Earthquake scaling law constant |

β |
Earthquake scaling law exponent |

μ |
Elastic shear modulus |

L |
Fault length |

W |
Fault width |

D |
Earthquake rupture mean displacement |

θ |
Fault dip angle |

V_{F} |
Average fault slip rate |

V_{S} |
Seismic fault slip rate |

V_{A} |
Aseismic fault slip rate |

χ |
Average degree of seismic coupling |

σ |
Variance of the normal depth distribution of mainshocks |

T |
Simulation duration |

t |
Simulation time |

x |
Along-fault coordinate |

y |
Along-river coordinate |

h |
Bedrock riverbed elevation |

h_{s} |
Sediment cover thickness |

V_{R} |
Knickpoint retreat rate |

r |
Retreat rate constant |

Δh |
Knickpoint height |

Δh_{0} |
Reference knickpoint height |

q |
Retreat rate exponent |

K |
Erodibility |

A |
Drainage area |

mn | Exponents of the stream power law |

φ |
Average river slope just upstream the fault |

PS wrote the paper and designed this study. PS and TC developed the accompanying numerical model. EB and DL motivated the paper through insightful discussions around river profiles and co-seismic knickpoints. All authors checked and revised the text and the figures of the paper, contributed to ideas developed in this study and discussed the implications for geomorphology and river profile analysis.

The authors declare that they have no conflict of interest.

We thank Louise Jeandet, Maxime Mouyen, Michaël Pons, Rodolphe Cattin and Philippe Davy for their helpful comments and for discussions about this work. We also thank Wolfgang Schwanghart, Robert Sare, George Hilley and an anonymous reviewer as well as the editor for their many positive and constructive comments that have contributed to enhance the quality of the manuscript.

This research has been supported by the Agence Nationale de la Recherche (grant no. ANR-14-CE33-0005) and the Conseil Régional de Bretagne (grant no. Boost'Europe EROQUAKE (9737)). This project has received funding from the European Research Council (ERC) under the European Union's Horizon Horizon 2020 research and innovation programme (grant agreement no. 803721). We also acknowledge support by Université Rennes 1, the CNRS and the French-Taiwanese International Laboratory D3E.

This paper was edited by Richard Gloaguen and reviewed by Wolfgang Schwanghart, Robert Sare, George Hilley, and one anonymous referee.

Aki, K.: A probabilistic synthesis of precursory phenomena, in: Earthquake prediction: an international review, edited by: Simpson, D. W. and Richards, P. G., AGU, Washington, D.C., USA, 4, 566–574, 1981.

Antón, L., Mather, A. E., Stokes, M., Munoz-Martin, A., and De Vicente, G.: Exceptional river gorge formation from unexceptional floods, Nat. Commun., 6, 7963, https://doi.org/10.1038/ncomms8963, 2015.

Arrowsmith, J. R., Pollard, D. D., and Rhodes, D. D.: Hillslope development in areas of active tectonics, J. Geophys. Res.-Sol. Ea., 101, 6255–6275, 1996.

Arrowsmith, J. R., Rhodes, D. D., and Pollard, D. D.: Morphologic dating of scarps formed by repeated slip events along the San Andreas Fault, Carrizo Plain, California, J. Geophys. Res.-Sol. Ea., 103, 10141–10160, 1998.

Attal, M. and Lavé, J.: Changes of bedload characteristics along the Marsyandi River (central Nepal): Implications for understanding hillslope sediment supply, sediment load evolution along fluvial networks, and denudation in active orogenic belts, in: Tectonics, climate, and landscape evolution, Geological Society of America Special paper, 398, 143–171, https://doi.org/10.1130/2006.2398(09), 2006.

Avouac, J. P.: Analysis of scarp profiles: evaluation of errors in morphologic dating, J. Geophys. Res.-Sol. Ea., 98, 6745–6754, 1993.

Avouac, J. P.: From geodetic imaging of seismic and aseismic fault slip to dynamic modeling of the seismic cycle, Annu. Rev. Earth Pl. Sc., 43, 233–271, 2015.

Baynes, E. R., Attal, M., Niedermann, S., Kirstein, L. A., Dugmore, A. J., and Naylor, M.: Erosion during extreme flood events dominates Holocene canyon evolution in northeast Iceland, P. Natl. Acad. Sci. USA, 112, 2355–2360, https://doi.org/10.1073/pnas.1415443112, 2015.

Baynes, E. R., Lague, D., Attal, M., Gangloff, A., Kirstein, L. A., and Dugmore, A. J.: River self-organisation inhibits discharge control on waterfall migration, Sci. Rep., 8, 2444, https://doi.org/10.1038/s41598-018-20767-6, 2018.

Berlin, M. M. and Anderson, R. S.: Modeling of knickpoint retreat on the Roan Plateau, western Colorado, J. Geophys. Res.-Earth, 112, F03S06, https://doi.org/10.1029/2006JF000553, 2007.

Bishop, P., Hoey, T. B., Jansen, J. D., and Artza, I. L.: Knickpoint recession rate and catchment area: the case of uplifted rivers in Eastern Scotland, Earth Surf. Proc. Land., 30, 767–778, 2005.

Boncio, P., Liberi, F., Caldarella, M., and Nurminen, F.-C.: Width of surface rupture zone for thrust earthquakes: implications for earthquake fault zoning, Nat. Hazards Earth Syst. Sci., 18, 241–256, https://doi.org/10.5194/nhess-18-241-2018, 2018.

Boulton, S. J. and Whittaker, A. C.: Quantifying the slip rates, spatial distribution and evolution of active normal faults from geomorphic analysis: Field examples from an oblique-extensional graben, southern Turkey, Geomorphology, 104, 299–316, 2009.

Braun, J.: The many surface expressions of mantle dynamics, Nat. Geosci., 3, https://doi.org/10.1038/ngeo1020, 2010.

Braun, J. and Sambridge, M.: Modelling landscape evolution on geological time scales: a new method based on irregular spatial discretization, Basin Res., 9, 27–52, 1997.

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, 170–179, 2013.

Brocard, G. Y., Willenbring, J. K., Miller, T. E., and Scatena, F. N.: Relict landscape resistance to dissection by upstream migrating knickpoints, J. Geophys. Res.-Earth, 121, 1182–1203, 2016.

Campforts, B., Schwanghart, W., and Govers, G.: Accurate simulation of transient landscape evolution by eliminating numerical diffusion: the TTLEM 1.0 model, Earth Surf. Dynam., 5, 4–66, https://doi.org/10.5194/esurf-5-47-2017, 2017.

Carretier, S. and Lucazeau, F.: How does alluvial sedimentation at range fronts modify the erosional dynamics of mountain catchments?, Basin Res., 17, 361–381, 2005.

Castelltort, S., Goren, L., Willett, S. D., Champagnac, J. D., Herman, F., and Braun, J.: River drainage patterns in the New Zealand Alps primarily controlled by plate tectonic strain, Nat. Geosci., 5, https://doi.org/10.1038/NGEO1582, 2012.

Cattin, R. and Avouac, J. P.: Modeling mountain building and the seismic cycle in the Himalaya of Nepal, J. Geophys. Res.-Sol. Ea., 105, 13389–13407, 2000.

Chen, Y. G., Chen, W. S., Lee, J. C., Lee, Y. H., Lee, C. T., Chang, H. C., and Lo, C. H.: Surface rupture of 1999 Chi-Chi earthquake yields insights on active tectonics of central Taiwan, B. Seismol. Soc. Am., 91, 977–985, 2001.

Cook, K. L., Turowski, J. M., and Hovius, N.: A demonstration of the importance of bedload transport for fluvial bedrock erosion and knickpoint propagation, Earth Surf. Proc. Land., 38, 683–695, 2013.

Croissant, T., Lague, D., Steer, P., and Davy, P.: Rapid post-seismic landslide evacuation boosted by dynamic river width, Nat. Geosci., 10, https://doi.org/10.1038/NGEO3005, 2017.

Croissant, T., Steer, P., Lague, D., Davy, P., Jeandet, L., and Hilton, R., G.: Seismic cycles, earthquakes, landslides and sediment fluxes: Linking tectonics to surface processes using a reduced-complexity model, Geomorphology, 339, 87–103 2019.

Crosby, B. T. and Whipple, K. X.: Knickpoint initiation and distribution within fluvial networks: 236 waterfalls in the Waipaoa River, North Island, New Zealand, Geomorphology, 82, 16–38, 2006.

Davy, P., Croissant, T., and Lague, D.: A precipiton method to calculate river hydrodynamics, with applications to flood prediction, landscape evolution models, and braiding instabilities, J. Geophys. Res.-Earth, 122, 1491–1512, 2017.

Dorsey, R. J. and Roering, J. J.: Quaternary landscape evolution in the San Jacinto fault zone, Peninsular Ranges of Southern California: Transient response to strike-slip fault initiation, Geomorphology, 73, 16–32, 2006.

Egholm, D. L., Knudsen, M. F., Clark, C. D., and Lesemann, J. E.: Modeling the flow of glaciers in steep terrains: The integrated second-order shallow ice approximation (iSOSIA), J. Geophys. Res.-Earth, 116, F02012, https://doi.org/10.1029/2010JF001900, 2011.

Ewiak, O., Victor, P., and Oncken, O.: Investigating multiple fault rupture at the Salar del Carmen segment of the Atacama Fault System (northern Chile): Fault scarp morphology and knickpoint analysis, Tectonics, 34, 187–212, 2015.

Finnegan, N. J. and Balco, G.: Sediment supply, base level, braiding, and bedrock river terrace formation: Arroyo Seco, California, USA, GSA Bulletin, 125, 1114–1124, 2013.

Gailleton, B., Mudd, S. M., Clubb, F. J., Peifer, D., and Hurst, M. D.: A segmentation approach for the reproducible extraction and quantification of knickpoints from river long profiles, Earth Surf. Dynam., 7, 211–230, https://doi.org/10.5194/esurf-7-211-2019, 2019.

Gilbert, G. K.: Niagara Falls and their history, in: National Geographic Society, The Physiography of the United States, The American Book Co., New York, USA, 203–236, 1896.

Gilbert, G. K.: Rate of recession of Niagara Falls, in: U.S. Geological Survey Bulletin, 306, 1907.

Goren, L., Castelltort, S., and Klinger, Y.: Modes and rates of horizontal deformation from rotated river basins: Application to the Dead Sea fault system in Lebanon, Geology, 43, 843–846, 2015.

Gutenberg, B. and Richter, C. F.: Frequency of earthquakes in California, B. Seismol. Soc. Am., 34, 185–188, 1944.

Hack, J. T.: Studies of longitudinal stream profiles in Virginia and Maryland, U.S. Geological Survey Bulletin, 294, 1957.

Haviv, I., Enzel, Y., Whipple, K. X., Zilberman, E., Matmon, A., Stone, J., and Fifield, K. L.: Evolution of vertical knickpoints (waterfalls) with resistant caprock: Insights from numerical modelling, J. Geophys. Res.-Earth, 115, F03028, https://doi.org/10.1029/2008JF001187, 2010.

Hayakawa, Y. and Matsukura, Y.: Recession rates of waterfalls in Boso Peninsula, Japan, and a predictive equation, Earth Surf. Proc. Land., 28, 675–684, 2003.

He, Z. and Ma, B.: Holocene paleoearthquakes of the Daqingshan fault detected from knickpoint identification and alluvial soil profile, J. Asian Earth Sci., 98, 261–271, 2015.

Helmstetter, A. and Sornette, D.: Foreshocks explained by cascades of triggered seismicity, J. Geophys. Res.-Sol. Ea., 108, 2457, https://doi.org/10.1029/2003JB002409, 2003.

Herman, F. and Braun, J.: Evolution of the glacial landscape of the Southern Alps of New Zealand: Insights from a glacial erosion model, J. Geophys. Res.-Earth, 113, F02009, https://doi.org/10.1029/2007JF000807, 2008.

Holland, W. N. and Pickup, G., Flume study of knickpoint development in stratified sediment, Geol. Soc. Am. Bull., 87, 76–82, 1976.

Hovius, N., Meunier, P., Lin, C. W., Chen, H., Chen, Y. G., Dadson, S., Ming-Jame, H., and Lines, M.: Prolonged seismically induced erosion and the mass balance of a large earthquake, Earth Planet. Sc. Let., 304, 347–355, 2011.

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

Howard, A. D., Dietrich, W. E., and Seidl, M. A.: Modeling fluvial erosion on regional to continental scales, J. Geophys. Res.-Sol. Ea., 99, 13971–13986, 1994.

Howarth, J. D., Fitzsimons, S. J., Norris, R. J., and Jacobsen, G. E.: Lake sediments record cycles of sediment flux driven by large earthquakes on the Alpine fault, New Zealand, Geology, 40, 1091–1094, 2012.

Jansen, J. D., Fabel, D., Bishop, P., Xu, S., Schnabel, C., and Codilean, A. T.: Does decreasing paraglacial sediment supply slow knickpoint retreat?, Geology, 39, 543–546, 2011.

Jolivet, R., Simons, M., Agram, P. S., Duputel, Z., and Shen, Z. K.: Aseismic slip and seismogenic coupling along the central San Andreas Fault, Geophys. Res. Lett., 42, 297–306, 2015.

Kanamori, H.: The energy release in great earthquakes, J. Geophys. Res., 82, 2981–2987, 1977.

Keefer, D. K.: The importance of earthquake-induced landslides to long-term slope erosion and slope-failure hazards in seismically active regions, in: Geomorphology and Natural Hazards, 265–284, 1994.

King, G.: The accommodation of large strains in the upper lithosphere of the earth and other solids by self-similar fault systems: the geometrical origin of b-value, Pure Appl. Geophys., 121, 761–815, 1983.

Lague, D.: The stream power river incision model: evidence, theory and beyond, Earth Surf. Proc. Land, 39, 38–61, 2014.

Lamb, M. P. and Dietrich, W. E.: The persistence of waterfalls in fractured rock, Geol. Soc. Am. Bull., 121, 1123–1134, 2009.

Leonard, M.: Earthquake fault scaling: Self-consistent relating of rupture length, width, average displacement, and moment release, B. Seismol. Soc. Am., 100, 1971–1988, 2010.

Loget, N. and Van Den Driessche, J.: Wave train model for knickpoint migration, Geomorphology, 106, 376–382, 2009.

Loget, N., Davy, P., and Van den Driessche, J.: Mesoscale fluvial erosion parameters deduced from modeling the Mediterranean sea level drop during the Messinian (late Miocene), J. Geophys. Res.-Earth, 111, F03005, https://doi.org/10.1029/2005JF000387, 2006.

Malatesta, L. C. and Lamb, M. P.: Formation of waterfalls by intermittent burial of active faults, GSA Bulletin, 130, 522–536, 2018.

Marc, O., Hovius, N., and Meunier, P.: The mass balance of earthquakes and earthquake sequences, Geophys. Res. Lett., 43, 3708–3716, 2016.

Marone, C. and Scholz, C. H.: The depth of seismic faulting and the upper transition from stable to unstable slip regimes, Geophys. Res. Lett., 15, 621–624, 1988.

McKean, J. A., Dietrich, W. E., Finkel, R. C., Southon, J. R., and Caffee, M. W.: Quantification of soil production and downslope creep rates from cosmogenic 10Be accumulations on a hillslope profile, Geology, 21, 343–346, 1993.

Miller, S. R., Slingerland, R. L., and Kirby, E.: Characteristics of steady state fluvial topography above fault-bend folds, J. Geophys. Res.-Earth, 112, F04004, https://doi.org/10.1029/2007JF000772, 2007.

Nash, D. B.: Morphologic dating of degraded normal fault scarps, J. Geol., 88, 353–360, 1980.

Neely, A. B., Bookhagen, B., and Burbank, D. W.: An automated knickzone selection algorithm (KZ-Picker) to analyze transient landscapes: Calibration and validation, J. Geophys. Res.-Earth, 122, 1236–1261, 2017.

Ogata, Y.: Statistical models for earthquake occurrences and residual analysis for point processes, J. Am. Stat. Assoc., 83, 9–27, 1988.

Okada, Y.: Surface deformation due to shear and tensile faults in a half-space, B. Seismol. Soc. Am., 75, 1135–1154, 1985.

Parker, R. N., Densmore, A. L., Rosser, N. J., De Michele, M., Li, Y., Huang, R., Whadcoat, S., and Petley, D. N.: Mass wasting triggered by the 2008 Wenchuan earthquake is greater than orogenic growth, Nat. Geosci., 4, https://doi.org/10.1038/NGEO1154, 2011.

Parker, R. S.: Experimental study of drainage basin evolution and its hydrologic implications, PhD dissertation, Colorado State University, Fort Collins, 1977.

Peng, Z. and Gomberg, J.: An integrated perspective of the continuum between earthquakes and slow-slip phenomena, Nat. Geosci., 3, https://doi.org/10.1038/ngeo940, 2010.

Quigley, M. C., Cupper, M. L., and Sandiford, M.: Quaternary faults of south-central Australia: palaeoseismicity, slip rates and origin, Australian J. Earth Sci., 53, 285–301, 2006.

Roering, J. J., Kirchner, J. W., and Dietrich, W. E.: Evidence for nonlinear, diffusive sediment transport on hillslopes and implications for landscape morphology, Water Resour. Res., 35, 853–870, 1999.

Romanet, P., Bhat, H. S., Jolivet, R., and Madariaga, R.: Fast and slow slip events emerge due to fault geometrical complexity, Geophys. Res. Lett., 45, 4809–4819, 2018.

Rosenbloom, N. A. and Anderson, R. S.: Hillslope and channel evolution in a marine terraced landscape, Santa Cruz, California, J. Geophys. Res.-Sol. Ea., 99, 14013–14029, 1994.

Royden, L. and Perron, J. T., Solutions of the stream power equation and application to the evolution of river longitudinal profiles, J. Geophys. Res.-Earth, 118, 497–518, 2013.

Scheingross, J. S. and Lamb, M. P.: A mechanistic model of waterfall plunge pool erosion into bedrock, J. Geophys. Res.-Earth, 122, 2079–2104, 2017.

Schmittbuhl, J., Karabulut, H., Lengliné, O., and Bouchon, M.: Seismicity distribution and locking depth along the Main Marmara Fault, Turkey, Geochem. Geophy. Geosy., 17, 954–965, 2015.

Scholz, C. H.: Earthquakes and friction laws, Nature, 391, 6662, https://doi.org/10.1038/34097, 1998.

Schumm, S. A., Mosley, M. P., and Weaver, W. E.: Experimental Fluvial Geomorphology, in: John Wiley and Sons, New York, 413 pp., 1987.

Shcherbakov, R. and Turcotte, D. L.: A modified form of Bath's law, B. Seismol. Soc. Am., 94, 1968–1975, 2004.

Shcherbakov, R., Turcotte, D. L., and Rundle, J. B.: A generalized Omori's law for earthquake aftershock decay, Geophys. Res. Lett., 31, L11613, https://doi.org/10.1029/2004GL019808, 2004.

Sibson, R. H.: Fault zone models, heat flow, and the depth distribution of earthquakes in the continental crust of the United States, B. Seismol. Soc. Am., 72, 151–163, 1982.

Steer, P. and Croissant, T.: philippesteer/RiverFault: Initial version, https://doi.org/10.5281/zenodo.2654819, 2019.

Steer, P., Cattin, R., Lavé, J., and Godard, V.: Surface Lagrangian Remeshing: A new tool for studying long term evolution of continental lithosphere from 2D numerical modelling, Comput. Geosci., 37, 1067–1074, 2011.

Steer, P., Simoes, M., Cattin, R., and Shyu, J. B. H.: Erosion influences the seismicity of active thrust faults, Nat. Commun., 5, 5564, https://doi.org/10.1038/ncomms6564, 2014.

Stock, J. D. and Montgomery, D. R.: Geologic constraints on bedrock river incision using the stream power law, J. Geophys. Res.-Sol. Ea., 104, 4983–4993, 1999.

Sun, C., Wan, T., Xie, X., Shen, X., and Liang, K.: Knickpoint series of gullies along the Luoyunshan Piedmont and its relation with fault activity since late Pleistocene, Geomorphology, 268, 266–274, 2016.

Thieulot, C., Steer, P., and Huismans, R. S.: Three-dimensional numerical simulations of crustal systems undergoing orogeny and subjected to surface processes, Geochem. Geophys. Geosy., 15, 4936–4957, 2014.

Tucker, G. E. and Bradley, D. N.: Trouble with diffusion: Reassessing hillslope erosion laws with a particle-based model, J. Geophys. Res.-Earth, 115, F00A10, https://doi.org/10.1029/2009JF001264, 2010.

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

Turcotte, D. L., Holliday, J. R., and Rundle, J. B.: BASS, an alternative to ETAS, Geophys. Res. Lett., 34, L12303, https://doi.org/10.1029/2007GL029696, 2007.

van der Beek, P., Pulford, A., and Braun, J.: Cenozoic landscape development in the Blue Mountains (SE Australia): lithological and tectonic controls on rifted margin morphology, J. Geol., 109, 35–56, 2001.

Van Heijst, M. W. I. M. and Postma, G.: Fluvial response to sea-level changes: a quantitative analogue, experimental approach, Basin Res., 13, 269–292, 2001.

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.-Earth, 115, F02010, https://doi.org/10.1029/2008JF001079, 2010.

Vernant, P., Hivert, F., Chery, J., Steer, P., Cattin, R., and Rigo, A.: Erosion-induced isostatic rebound triggers extension in low convergent mountain ranges, Geology, 41, 467–470, 2013.

Watts, A. B.: Isostasy and Flexure of the Lithosphere, in: Cambridge University Press, Cambridge, UK, 2001.

Wei, Z., Bi, L., Xu, Y., and He, H.: Evaluating knickpoint recession along an active fault for paleoseismological analysis: The Huoshan Piedmont, Eastern China, Geomorphology, 235, 63–76, 2015.

Wells, D. L. and Coppersmith, K. J.: New empirical relationships among magnitude, rupture length, rupture width, rupture area, and surface displacement, B. Seismol. Soc. Am., 84, 974–1002, 1994.

Whipple, K. X.: The influence of climate on the tectonic evolution of mountain belts, Nat. Geosci., 2, https://doi.org/10.1038/ngeo413, 2009.

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

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

Willett, S. D.: Orogeny and orography: The effects of erosion on the structure of mountain belts, J. Geophys. Res.-Sol. Ea., 104, 28957–28981, 1999.

Yanites, B. J., Tucker, G. E., Mueller, K. J., Chen, Y. G., Wilcox, T., Huang, S. Y., and Shi, K. W.: Incision and channel morphology across active structures along the Peikang River, central Taiwan: Implications for the importance of channel width, GSA Bulletin, 122, 1192–1208, 2010.

Yildirim, C., Schildgen, T. F., Echtler, H., Melnick, D., and Strecker, M. R.: Late Neogene and active orogenic uplift in the Central Pontides associated with the North Anatolian Fault: Implications for the northern margin of the Central Anatolian Plateau, Turkey, Tectonics, 30, TC5005, https://doi.org/10.1029/2010TC002756, 2011.

Zinke, R., Hollingsworth, J., and Dolan, J. F.: Surface slip and off-fault deformation patterns in the 2013 Mw 7.7 Balochistan, Pakistan earthquake: Implications for controls on the distribution of near-surface coseismic slip, Geochem. Geophy. Geosy., 15, 5034–5050, 2014.

- Abstract
- Introduction
- State of the art: linking fault slip to knickpoint formation and migration
- Methods
- Magnitude, displacement and temporal distributions of earthquakes and co-seismic knickpoints
- Knickpoints along single river profiles
- Knickpoints along successive parallel rivers
- Knickpoint detectability
- Discussion
- Conclusions
- Code availability
- Appendix A: Constant knickpoint retreat rate and the stream power law
- Appendix B: Parallel rivers' topographic correlation
- Appendix C: Table of variable definitions and notations
- Author contributions
- Competing interests
- Acknowledgements
- Financial support
- Review statement
- References

- Abstract
- Introduction
- State of the art: linking fault slip to knickpoint formation and migration
- Methods
- Magnitude, displacement and temporal distributions of earthquakes and co-seismic knickpoints
- Knickpoints along single river profiles
- Knickpoints along successive parallel rivers
- Knickpoint detectability
- Discussion
- Conclusions
- Code availability
- Appendix A: Constant knickpoint retreat rate and the stream power law
- Appendix B: Parallel rivers' topographic correlation
- Appendix C: Table of variable definitions and notations
- Author contributions
- Competing interests
- Acknowledgements
- Financial support
- Review statement
- References