# HyLands 1.0: a Hybrid Landscape evolution model to simulate the impact of landslides and landslide-derived sediment on landscape evolution

Benjamin Campforts, Charles M. Shobe, Philippe Steer, Matthias

Vanmaercke, Dimitri Lague, Jean Braun

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

## To cite this version:

Benjamin Campforts, Charles M. Shobe, Philippe Steer, Matthias Vanmaercke, Dimitri Lague, et al. HyLands 1.0: a Hybrid Landscape evolution model to simulate the impact of landslides and landslide-derived sediment on landscape evolution. Geoscientific Model Development, 2020, 13 (9), pp.3863-3886. 10.5194/gmd-2020-74. insu-02539385

HAL Id: insu-02539385

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

Submitted on 10 Apr 2020

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

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

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

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

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

# HyLands 1.0: a Hybrid Landscape evolution model to simulate the impact of landslides and landslide-derived sediment on landscape evolution

Benjamin Campforts¹,², Charles M. Shobe¹, Philippe Steer³, Matthias Vanmaercke⁴, Dimitri Lague³, and Jean Braun¹,⁵

¹Helmholtz Centre Potsdam, GFZ German Research Centre for Geosciences, Potsdam, Germany

²Institute for Arctic and Alpine Research, University of Colorado at Boulder, Boulder, CO, USA

³Univ Rennes, CNRS, Géosciences Rennes – UMR 6118, 35000 Rennes, France

⁴Université de Liège, UR SPHERES, Département de Géographie, Clos Mercator 3, 4000 Liège, Belgium

⁵Institute of Earth and Environmental Science, Universität Potsdam, Potsdam, Germany

**Correspondence:** Benjamin Campforts (benjamin.campforts@gfz-potsdam.de)

**Abstract.** Landslides are the main source of sediment in most mountain ranges. Rivers then act as conveyor belts, evacuating landslide-derived sediment. Sediment dynamics are known to influence landscape evolution through interactions among landslide sediment delivery, fluvial transport, and river incision into bedrock. Sediment delivery and its interaction with river incision therefore control the pace of landscape evolution and mediate relationships among tectonics, climate, and erosion. Numerical landscape evolution models (LEMs) are well suited to study the interaction among these earth surface processes. They enable evaluation of a range of hypotheses at varying temporal and spatial scales. While many models have been used to study the dynamic interplay between tectonics, erosion and climate, the role of interactions between landslide-derived sediment and river incision has received much less attention. Here, we present HyLands, a hybrid landscape evolution model integrated within the Topo Toolbox Landscape Evolution Model (TTLEM) framework. The hybrid nature of the model lies in its capacity to simulate both erosion and deposition at any place in the landscape due to fluvial bedrock incision, sediment transport and rapid, stochastic mass wasting through landsliding. Fluvial sediment transport and bedrock incision are calculated using the recently developed Stream Power with Alluvium Conservation and Entrainment (SPACE) model. Therefore, rivers in HyLands can dynamically transition from detachment-limited to transport-limited, and from bedrock to bedrock-alluvial to fully alluviated states. Erosion and sediment production by landsliding is calculated using a Mohr-Coulomb stability analysis while landslide-derived sediment is routed and deposited using a multiple flow direction, non-linear deposition method. We describe and evaluate the HyLands 1.0 model using analytical solutions and observations. We first illustrate the functionality of HyLands to capture river dynamics ranging from detachment-limited to transport-limited configurations. Second, we apply the model to a portion of the Namche-Barwa massif in Eastern Tibet and compare simulated and observed landslide magnitude-frequency and area-volume scaling relationships. Finally, we illustrate the relevance of explicitly simulating landsliding and sediment dynamics over longer timescales for landscape evolution in general and river dynamics in particular. With HyLands we pro-
<figure><img src="image_5.png" /></figure>

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

vide a new tool to understand both the long and short-term coupling between stochastic hillslope processes, river incision, and source-to-sink sediment dynamics.

# 1 Introduction

Landsliding is a highly effective erosional mechanism that dominates sediment mobilization rates in moderate to steep topographic settings (Hovius et al., 1997; Ouimet et al., 2007; Broeckx et al., 2020). Nonetheless, long term landscape evolution in non-glaciated settings is mainly controlled by the interplay between tectonic uplift and fluvial dynamics (Whipple and Tucker, 1999; Wobus et al., 2006). Fluvial channels in mountainous catchments play a dual role: they simultaneously incise into the bedrock and act as conveyor belts to carry eroded sediment out of the mountain range towards the ocean (Milliman and Meade, 1983). Through sediment evacuation and bedrock incision, fluvial incision lowers the base level for surrounding hillslopes, triggering hillslope failures. In turn, hillslope failure through mass wasting chokes the rivers with sediment and prevents further bedrock incision until landslide derived sediment has been evacuated from the system (Larsen and Montgomery, 2012; Ouimet et al., 2007; Korup et al., 2010; Shobe et al., 2016; Glade et al., 2019).

Unravelling the dynamic interplay between landslides and fluvial processes is key to understanding long-term landscape evolution and the associated sediment dynamics in mountainous terrain (Egholm et al., 2013). Increased insight into the spatial distribution of landslides has resulted in improved landslide susceptibility assessments (Guzzetti et al., 2006), but processes regulating landslide rate assessments (Broeckx et al., 2020) and landslide-derived sediment dynamics remain less well understood (Hovius et al., 2011; Croissant et al., 2017, 2019; Zhang et al., 2019; Broeckx et al., 2020).

Numerical models are excellent tools to study relationships between processes regulating Earth surface dynamics and their inter-dependencies over various temporal and spatial scales (Tucker and Hancock, 2010). The past twenty years have seen the development of a plethora of Landscape Evolution Models (LEMs), enabling studies of the interactions among climate, tectonics, and erosion. A crucial ingredient for any LEM is a fluvial erosion component regulating the way in which rivers transport sediment and incise into bedrock. Fluvial incision is controlled by both water and sediment cascading through river channels (Whipple et al., 2000; Hancock and Anderson, 2002; Turowski et al., 2007). Most existing LEMs simulate river incision using one of two commonly used end member models to simulate fluvial dynamics (Armitage et al., 2018). In one approach, river incision is simulated assuming a detachment-limited configuration where erosion is constrained by the power to erode particles from the river bed and quantified using a scaling law between fluid stress and river incision rate (Seidl and Dietrich, 1992; Howard and Kerby, 1983; Campforts and Govers, 2015). In the other approach, river incision is simulated assuming a transport-limited configuration where erosion is constrained by the capacity of the river to carry sediment, where the carrying capacity is a function of the fluid stress (Willgoose et al., 1991; Paola and Voller, 2005). These two formulations lead to similar outcomes in steady-state channels (where the river erosion rate equals the rock uplift rate), but noticeably different outcomes during transient river response to tectonic and climatic perturbations (Whipple and Tucker, 1999).

In real settings, however, even steep mountain channels undergoing long-term bedrock incision may experience bed cover by alluvial sediment. Further, over geologic time as tectonic and climatic forcings change, it is likely that any given channel
<figure><img src="image_7.png" /></figure>

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

55

transitions between detachment-limited and transport-limited behavior. Such heterogeneous configurations require a model setup that can dynamically transition between detachment-limited and transport-limited regimes (e.g., Davy and Lague, 2009) and can simultaneously simulate fluvial sediment transport and river incision into bedrock (e.g., Lague, 2010). Recently, the SPACE (Stream Power with Alluvium Conservation and Entrainment) model approach was proposed to meet both of these needs (Shobe et al., 2017). Because SPACE is purely a river incision model, it does not simulate hillslope or mass wasting processes. Additional model components are therefore needed to simulate the impact of mass wasting on landscape evolution and sediment dynamics.

60

To understand how landslides influence landscape evolution, Densmore et al. (1998) proposed an approach, adapted by others (e.g. Champel, 2002; Egholm et al., 2013), to integrate stochastic landslide dynamics in a numerical landscape evolution model. Densmore et al. (1998) assume that (i) all hillslopes behave as Mohr-Coulomb materials (Taylor, 1948), (ii) landslides initialize in river channels (i.e. at the base of hillslopes) and (iii) landslide-derived sediment is spread under a constant slope, following the steepest downslope path. Assuming that landslides initialize only in fluvial channels makes it computationally easier to implement a landslide model because sediment pathways on hillslopes need not be calculated. However, it is not realistic to neglect landslides that might initiate away from river channels: it has been shown that a large portion of landslide-derived sediment is stored along flow paths on hillslopes, rather than being immediately delivered to river channels (Broeckx et al., 2020; Hoffmann, 2015). The approach of Densmore et al. (1998) does not allow for landslide-derived sediment to be deposited and spread over hillslopes. Rather, landslide debris is spread as tongues of sediment filling up the river channel.

65

70

Other researchers have developed mechanistic models to simulate shallow landslide activity at the landscape scale (e.g. Montgomery and Dietrich, 1994; Claessens et al., 2007). Such models typically involve the explicit simulation of a soil layer and a coupled hydrologic model to calculate how changing pore water pressures trigger landslides (Van Asch et al., 1999; Iverson, 2000; Baum et al., 2010). Although such mechanistic models are useful for assessing landslide hazards (e.g. to simulate landslide liquefaction associated with the Oso landslide, cfr. Iverson and George, 2016), they typically involve a range of geophysical processes and associated input parameters which are not always available at large spatial scales. This can make the more detailed models sensitive to equifinality (Beven and Freer, 2001). Moreover, deep-seated bedrock landslides, rather than shallow landslides, mobilize the largest volumes of sediment and therefore have the largest impact on landscape evolution (Burbank, 2002; Dussauge et al., 2003; Jeandet et al., 2019; Korup et al., 2007; Broeckx et al., 2020).

75

80

In this paper, we present HyLands, a new Hybrid Landscape evolution model for simulating the interaction of landslide dynamics and fluvial processes. The model is intended to simulate Earth surface evolution at large spatial scales with a special focus on landsliding and the long-term effects of landslide-derived sediment. HyLands is integrated in the TTLEM 1.0 landscape evolution model (Campforts et al., 2017). Unlike the existing implementation of TTLEM, HyLands is a fully mass conservative model where fluvial dynamics are modelled using the SPACE fluvial incision framework (Shobe et al., 2017) and hillslope-derived sediment fluxes are explicitly simulated. In this paper, we first describe the fluvial and landslide components of the HyLands model. We verify the fluvial model component by comparing model behavior against known analytical solutions. Subsequently, we evaluate the performance of the landslide module by applying HyLands to a selected region of the landslide-prone Namche-Barwa massif in Eastern Tibet. We show that HyLands reproduces observed landslide scaling

85
<figure><img src="image_9.png" /></figure>

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

relationships. Next, we apply the model to a synthetic case to illustrate the potential of HyLands for studying the dynamic interaction between landslide activity and fluvial dynamics. We do this by evaluating how a steady-state landscape responds to an imposed pulse of landsliding activity. Finally, we discuss the current model limitations, future perspectives and a range of potential applications.

## 2 HyLands model description

HyLands is a Matlab model code building on the existing Topo Toolbox Landscape Evolution Model (TTLEM, Campforts et al., 2017). It simulates changes in bedrock height and sediment thickness on a regular grid. The model is mass conservative; sediment produced by river incision and hillslope processes such as landsliding is explicitly simulated in the model. At every model iteration, the elevation of all grid cells is updated according to the following conservation statement for sediment and rock:

$$ \frac{\partial \eta}{\partial t} = \frac{\partial R}{\partial t} + \frac{\partial H}{\partial t} \\ = U - E_{r_{fluv}} + \left( \frac{D_{s_{fluv}} - E_{s_{fluv}}}{1 - \phi_{sed}} \right) \\ - E_{r_{hill}} + \left( \frac{D_{s_{hill}} - E_{s_{hill}}}{1 - \phi_{sed}} \right) \qquad (1) $$

where $\eta$ [L] is the topographic elevation given by the sum of the bedrock elevation R [L] and the bed sediment thickness H [L]. U [L/T] is the rock uplift rate and $\phi_{sed}$ is the bed sediment porosity. $E_{r_{fluv}}$ [L/T] is the fluvial volumetric erosion flux of bedrock per unit bed area, representing the amount of bedrock that is detached and entrained into the water column. $E_{s_{fluv}}$ [L/T] is the fluvial volumetric entrainment flux of sediment per unit bed area and $D_{s_{fluv}}$ [L/T] is the fluvial volumetric deposition flux of sediment per unit bed area. $E_{r_{hill}}$ [L/T] is the volumetric flux of hillslope bedrock erosion (landsiding) per unit bed area, representing the amount of bedrock that is detached. $E_{s_{hill}}$ [L/T] is the volumetric entrainment flux of sediment erosion (produced by landsliding or creep) per unit bed area and $D_{s_{hill}}$ [L/T] is the volumetric deposition flux of hillslope-derived sediment per unit bed area.

### 2.1 River sediment transport and bedrock erosion

HyLands uses the Stream Power with Alluvium Conservation and Entrainment (SPACE) river erosion model of Shobe et al. (2017). SPACE has two key advantages for the purposes of modeling river response to landslide sediment delivery. First, because of its derivation from the erosion-deposition family of models (e.g., Beaumont et al., 1992; Davy and Lague, 2009), it can dynamically shift between detachment-limited (erosion is limited by the rate of sediment or bedrock detachment from the bed) and transport-limited (erosion is limited by the capacity of the flow to move detached sediment) behavior. Second, it can simulate the full continuum of possible river bed compositions from bare bedrock channels to mixed bedrock-alluvial channels to fully alluvial channels. This is accomplished by combining mass conservation of river bed sediment with a bedrock erosion
<figure><img src="image_11.png" /></figure>

law to simultaneously solve for the time evolution of the bedrock and sediment surfaces. We implement the SPACE model equations in the TTLEM MATLAB modeling framework. For a full overview of the SPACE model and comparison with other models for coupled sediment and bedrock channel evolution, see Shobe et al. (2017).

## 2.1.1 Fluvial sediment and rock mass conservation

Conservation of sediment closely follows the erosion-deposition approach of Davy and Lague (2009), with the addition of terms that represent the entrainment of detached bedrock in the water column (Fig. 1). The spatial change in volumetric sediment flux $Q_{sfluv}$ [$L^3/T$] per unit width $w$ [$L$] is written as:

$$ \frac{\partial(Q_{sfluv}/w)}{\partial x} = E_{sfluv} + (1 - F_{fluv}) E_{rfluv} - D_{sfluv}. \qquad (2) $$

where $F_{fluv}$ is a unitless fraction of fine fluvial sediment. The factor $1 - F_{fluv}$ [−] represents the idea that a fraction of the bedrock particles detached from the bed may be small enough to stay in permanent suspension, and therefore should not be tracked as bed sediment.

## 2.1.2 Fluvial sediment entrainment, bedrock erosion, and sediment deposition

To evaluate the impact of landslide-derived sediment on landscape evolution, it is critical to have a model that simulates simultaneous sediment entrainment and bedrock erosion, and considers the influence of sediment cover on river erosion dynamics. In the SPACE model, sediment entrainment and bedrock erosion may occur simultaneously. Further, the magnitude of each process is set by the relative availability of sediment on the channel bed. SPACE accomplishes this by including the influence of sediment layer on sediment and bedrock erosion rates.

Sediment entrainment and bedrock erosion are both governed by a unit stream power expression in which erosive power is a function of water discharge per unit width $q$ [$L^2T^{-1}$] and local channel slope $S$ (e.g., Howard and Kerby, 1983; Whipple and Tucker, 1999; Davy and Lague, 2009). The sediment erosion rate $E_{sfluv}$ and the bedrock erosion rate $E_{rfluv}$ are modified by a term $H/H_*$ [−] that encapsulates the ratio of bed sediment thickness $H$ [$L$] to bedrock bed roughness $H_*$ [$L$]. High bed sediment thickness or low bedrock surface roughness leads to a condition in which $H/H_*$ is large and little bedrock is exposed to erosive flows. If bed sediment thickness is low or bedrock roughness is high, $H/H_*$ is small and most of the in-channel bedrock is exposed to the flow.

SPACE assumes an exponential increase in sediment entrainment rate with $H/H_*$ and a concomitant exponential decrease in bedrock erosion rate with $H/H_*$ (Fig. 2). Rates of sediment entrainment and bedrock erosion can therefore be written as (assuming a negligible erosion threshold; see Shobe et al. (2017) for equations that relax this assumption):

$$ E_{sfluv} = K_s q S^n (1 - e^{-H/H_*}) \qquad (3) $$

for sediment and

$$ E_{rfluv} = K_r q S^n e^{-H/H_*} \qquad (4) $$
<figure><img src="image_12.png" /></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 1. Sketch of fluvial SPACE component.** Model setup and variable definitions for the SPACE bedrock-alluvial river erosion model. Reproduced from Shobe et al. (2017). Entrainment and deposition of sediment, as well as erosion of bedrock, can occur simultaneously. This approach allows channels to dynamically transition among bedrock, bedrock-alluvial, and fully alluviated states. At a given stream power, the relative rate of sediment entrainment $E_{sfluv}$ and bedrock erosion $E_{rfluv}$ is set by the ratio of sediment thickness $H$ to the bedrock roughness height $H_*$ (Fig. 2).</figcaption></figure>

for bedrock. $K_s$ and $K_r$ [$L^{-1}$] are erodibility constants for sediment and rock, respectively. $n$ is a constant set to 1 for all simulations in this paper, but that need not be 1 for the SPACE model in general. There are a variety of ways to compute water discharge $q$. We use the common approach of calculating discharge as a function of drainage area such that $q = k_q A^m$, where $m$ is a scaling exponent and $k_q$ is a coefficient subsumed into the fluvial erosion coefficients $K_s$ and $K_r$.

Sediment deposition is implemented similar to Davy and Lague (2009) such that the deposition flux depends on sediment flux $Q_{sfluv}$ divided by the volumetric water discharge $Q$ [$L^3/T$] and the effective sediment settling velocity $V$ [$L/T$]:

$$D_s = \frac{Q_{sfluv}}{Q} V. \qquad (5)$$
<figure><img src="image_14.png" /></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 2.** Relative efficiency of fluvial sediment entrainment and bedrock erosion ($f(H/H_*)$) as a function of the ratio of sediment thickness to bedrock roughness $H/H_*$. $f(H/H_*)$ depends only on the ratio $H/H_*$, and varies between 0 and 1, and indicates the proportion of total stream power used to erode bedrock (solid line) or sediment (dashed line). As sediment thickness $H$ increases relative to the bedrock roughness length scale $H_*$, the sediment entrainment rate factor approaches 1 and the bedrock erosion rate factor approaches 0 because the bed becomes composed entirely of sediment and no bedrock is exposed. As sediment thickness declines relative to the bedrock roughness length scale, the bedrock erosion rate increases exponentially because more bedrock is exposed and the sediment entrainment rate declines exponentially due to a lack of available sediment. This approach implements the “cover effect,” in which the presence of sediment reduces bedrock erosion rates, but does not incorporate the “tools effect,” in which mobile sediment enhances bedrock erosion. Reproduced from Shobe et al. (2017).</figcaption></figure>

## 2.2 Landsliding

HyLands treats landslide erosion and deposition deterministically, but uses a stochastic approach to calculating landslide occurrence. HyLands simulates deep-seated gravitational landslides eroding simultaneously the sediment layer and the bedrock (erosion terms $E_{s,hill}$ and $E_{r,hill}$ respectively in Eq. 1). We assume that both the rock and the sediment layer behave as Mohr-Coulomb materials. In its current form, HyLands does not simulate shallow landslides where failure geometry is imposed by the depth and angle of soil-rock transitions. Landslide initiation does not involve a preceding triggering event (e.g. an earthquake) but is simulated using a probabilistic approach.

### 2.2.1 Landslide erosion

Following Densmore et al. (1998), we simulate landslide erosion using the Culmann theory for slope stability. Culmann (1875) proposes that hillslope failure will occur on the plane where the shear stress is balanced by the sliding resistance. Assuming Mohr-Coulomb materials, it has been shown that the failure plane with a dip $\theta_c$ bisects the local topographic slope $\beta$, and the
<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /><figcaption>**Figure 3.** Sketch of landslide algorithm in two dimensions. Landslide erosion (red shaded area) is calculated using the Culmann approach (Culmann, 1875). $\beta$ is the topographic slope, $\phi$ represents the angle of internal friction and $\theta$ is the inclination of the rupture plane. Deposition of landslide-derived sediment (green shaded area) is calculated using a non-local diffusion equation (Eq. 11, cfr. Carretier et al., 2016). $\delta$ is the minimal angle of the spreading slope under which landslide-derived sediment is distributed on the hillslope. This sketch illustrates a case where none of the landslide-derived sediment is in permanent suspension ($F_{f,hill} = 0$) and the amount of eroded sediment (red shaded area) equals the amount of deposited sediment (green shaded area). If the deposited volume creates a down-slope gradient which is lower than the minimum spreading angle, $\delta$, the slope of deposited volume is adjusted so that the spreading slope equals $\delta$. Probability for sliding is calculated as the ratio of the local hillslope height $H_s$ to the maximum hillslope height $H_c$ (Eq. 7, cfr. Densmore et al., 1998). The inset plot illustrates that $H_c$ depends on the rock strength (cohesion $C$ and internal friction angle $\phi$) and the topographic slope $\beta$. The plotted lines are calculated using Eq. 8 with $\rho$ is 2700 kg m$^{-3}$, and $g = 9.81 m s^{-2}$.</figcaption></figure>

material's angle of internal friction $\phi$ (Densmore et al., 1998; Champel, 2002):

$$ \theta_c = \frac{\beta + \phi}{2}. \qquad (6) $$

The implementation of the Culmann theory in HyLands is illustrated in Fig. 3. For all points within the landscape where a landslide is initialized, the failure plane dipping at $\theta_c$ is extended until it daylights (i.e., intersects the topographic surface).

Modeling and slide frequency and location depends critically on the identification of points in the landscape where landslides initiate. A wide variety of events ranging from co-seismic activity and peak ground acceleration (Meunier et al., 2007) over intense storm events (Marc et al., 2018) to human hillslope destabilisation (Guns and Vanacker, 2014) may trigger mass wasting through landslide activity. Although these triggers could be added, HyLands mainly aims to simulate the impact of topographic landscape configuration on landslide activity. Therefore, we follow Densmore et al. (1998) in identifying unstable grid nodes
<figure><img src="image_18.png" /></figure>

as points in the landscape where the topographic slope (β) exceeds the angle of internal friction (φ). For all unstable nodes, the probability for sliding, $p_{LS}$, is calculated as:

$$p_{LS} = \frac{H_s}{H_c} \qquad (7)$$

where $H_s$ is the local hillslope height calculated as difference between every pixel in the landscape and its highest neighbour (Fig. 3) and $H_c$ is the maximum stable hillslope height which is calculated as (Densmore et al., 1998; Champel, 2002):

$$H_c = \frac{4C}{\rho g} \frac{\sin\beta \cos\phi}{1 - \cos(\beta - \phi)} \qquad (8)$$

Here C is the cohesion [M L⁻¹ T⁻²], ρ is the rock density [M L⁻³] set to 2700 kg m⁻³, and g the gravitational acceleration (g = 9.81 ms⁻²). To simulate the random nature of landslides, grid nodes where landslides initiate are selected using a stochastic sampling scheme:

$$rnd \frac{dt}{t_{LS}} \begin{cases} > p_{LS} & \text{No Landsliding} \\ < p_{LS} & \text{Landsliding, select as critical node} \end{cases} \qquad (9)$$

where rnd is a random number between 0 and 1 and $t_{LS}$ is the return time for landslides with $t_{LS} >= dt$. Unstable nodes where a landslide is induced will further be referred to as critical nodes (Fig. 3). Every model iteration, Eq. 9 is updated for all nodes where the topographic gradient (β) exceeds the critical material friction angle (φ). From Eqs. 8 and 9, it follows that the probability for sliding depends on the topographic slope, β, and inversely correlates with the angle of internal friction, φ, and the cohesion of the material, C (see inset of Fig. 3). Cohesion is a scale-dependent variable and parameter values covering several orders of magnitude have been reported (Sidle and Ochiai, 2006). Jeandet et al. (2019) inverted several landslide inventories and found effective cohesion values ranging between 10 - 35 kPa. These values are lower than geomechanical values representing large-scale rock strength (e.g. Densmore et al., 1998; Champel, 2002) which is attributed to the decrease in rock cohesion in the vicinity of faults following earthquakes (Gallen et al., 2015). In our experiments, we will use values for cohesion in the range reported by Jeandet et al. (2019), which represent effective cohesion following an earthquake or a storm.

The landslide return time $t_{LS}$ controls the absolute number of critical nodes where landslides are initiated. If $t_{LS}$ equals the timestep $dt$, the number of landslides per timestep is solely controlled by the ratio $H_s/H_c$. When $t_{LS} >> H_s/H_c$, however, the number of landslides per timestep is reduced (Eq. 9). While the $H_s/H_c$ ratio controls the topographic location of landsliding onset and thus the landslide characteristics (size and volume), the landslide return time, $t_{LS}$ controls the absolute number of landslides and therefore overall landslide erosion rates.

HyLands enables the simulation of landslides at every location in the landscape. Every iteration, landslides are induced at critical nodes sampled using the probabilistic approach outlined above (Eq. 9, Figs. 3 and 4). We propose a recursive approach to calculate the magnitude of a single landslide. For every critical node, we build a stack of unstable (= sliding) nodes. The stack is initialized by adding the critical node. Next, a recursive procedure is applied until the stack is empty. This procedure exists of the following steps: (i) Select the first node from the stack (thus starting with the critical node). (ii) Evaluate all up-slope neighbouring nodes. If the elevation of a neighbouring node exceeds the elevation of the sliding plain defined by $\theta_c$, the
<figure><img src="image_19.png" /></figure>

node is identified as a sliding pixel and added to the stack of landslide pixels. (iii) Remove the first node from the stack. This procedure is repeated until the stack is empty. HyLands offers the possibility to set a maximum landslide area ($A_{LSmax}$). Once this maximum is achieved, or when no more pixels are added to the stack of landslide pixels, the landslide area is defined. All pixels inside this area are eroded to the elevation of the sliding plane, thereby adjusting both $E_{s_hill}$ and $E_{r_hill}$ of Eq. 1 for all pixels involved.

### 2.2.2 Flux of landslide derived sediment

The spatial change in hillslope derived volumetric sediment flux $Q_{sHill}$ [$L^3/T$] per unit width $w$ [$L$] is written as:

$$ \frac{\partial(Q_{sHill}/w)}{\partial x} = E_{sHill} + (1 - F_{fHill}) E_{rHill} - D_{sHill}. \qquad (10) $$

where $F_{fHill}$ is a unitless fraction of fine hillslope derived sediment. The factor $1 - F_{fHill}$ [-] represents the idea that some fraction of the hillslope derived sediment is instantaneously evacuated as dissolved or suspended sediment (Page et al., 1999; Hovius et al., 2000; Lin et al., 2008; Tenorio et al., 2018), and therefore should not be tracked as sediment. When $F_{fHill} = 0$, the system is fully mass conservative and all sediments produced by landslide activity contribute to the sediment flux (Eq. 1).

### 2.2.3 Deposition of landslide material

In HyLands, landslides can initiate at any point in the landscape (Fig. 4). A steepest descent flow routing algorithm is known to be unrealistic for flow and sediment redistribution on hillslopes (Pelletier, 2010). Moreover, the use of a constant spreading slope, as suggested by Densmore et al. (1998) does not take into account the topographic relief when redistributing landslide derived sediment. When using a constant sediment spreading slope, sediment deposited on flat parts of the landscape is spread out over a much longer distances than sediment deposited on steep parts, where the large difference between topographic slope and spreading slope can accommodate for large sediment volumes. This is not realistic as sediment travel distance should depend on topographic gradient: material traveling over steeper slopes should go farther, all else being equal (Roering et al., 1999; Campforts et al., 2016).

A common approach to simulate sediment transport and deposition on hillslopes while considering the topographic gradient, is the use of linear or nonlinear diffusion equations (Roering et al., 1999; Andrews and Hanks, 1985). However, such approach is not suited to simulate the distribution of landslide-derived sediment. While diffusion equations distribute sediment only between the neighbouring cells of a pixel, landslide derived sediment has run-out distances which can be significantly longer than a single grid cell (Claessens et al., 2007). Therefore, we adopt a non-linear, non-local deposition scheme for landslide derived sediment outlined by Carretier et al. (2016):

$$ D_{sHill} = \frac{Q_{sHill}/w}{L} \qquad (11) $$

where $D_{sHill}$ [$L/T$] is the volumetric deposition flux of hillslope derived sediment per unit bed area and $L$ [$L$] represents a sediment transport distance. The larger $L$, the bigger the distance over which sediments are transported and the lower the local
<figure><img src="image_20.png" /></figure>

deposition rate. L is calculated for every grid cell as:

$$L = \frac{dx}{1 - \left(\frac{S}{Sc}\right)^2} \qquad (12)$$

where $S_c$ is a critical slope, which we further assume to be equal to the angle of internal friction ($\phi$). When the hillslope gradient $S \ll S_c$, most of the incoming sediment will be deposited and the resulting outcome is similar to the one obtained using a regular diffusion equation, also referred to as a local solution (Furbish and Roering, 2013; Carretier et al., 2016). When $S$ approaches $S_c$, $L$ goes to infinity implying that no deposition will occur at the considered cell. At steep slopes, sediment transport therefore shows non-local behaviour in the sense that erosion activity of non-local, upstream cells is integrated when calculating the sediment flux $Q_s$ (Carretier et al., 2016). When $S > S_c$, $L$ is set to $inf$. For negative values of $S$, which might occur for flooded nodes, $L$ is set to $dx$. In HyLands, a minimal spreading angle ($\delta$) can still be imposed under which landslide derived sediments are deposited but is not required (Fig. 3).

Contrary to fluvial dynamics (SPACE) where a single flow direction algorithm (steepest descent) is used, landslide-derived sediment is spread over the hillslopes using a multiple flow direction algorithm redistributing sediment over the downstream cells in proportion to the local slope (Fig. 4, cfr. Carretier et al., 2016). When a lot of sediment is debouched into fluvial channels, rivers can be blocked by landslide dams. HyLands uses a lake identification algorithm to identify flooded nodes during every model iteration. Lakes are identified by filling all sinks in a landscape to the brim. By default, flooded nodes do not erode but do allow for sediment deposition (Eq. 5).

In the remainder of this paper, we will first evaluate the performance of the fluvial and landslide components of HyLands through a set of verification and validation runs. Next, the coupling between landslide activity and long term landscape evolution will be evaluated using a synthetic model setup where a steady-state landscape is exposed to a pulse of landslides. All model experiments executed in the framework of this paper are available as executable Matlab scripts and as dynamic landscape evolution movies (Table 2).

# 3 Verification and evaluation

## 3.1 Comparison to analytical solutions for the fluvial dynamics

In the first three test runs (detachment, transport-limited and mixed), a steady-state artificial landscape is simulated using a square grid of 20 by 20 nodes with a coarse spatial resolution of 100 m. The run is initialized from a surface with randomly generated microtopography. The initial surface is a tilted plain which drains towards the southwestern corner, the only open boundary node. Therefore sediment and water can only leave the domain through this southwestern corner. The setup is identical to the one proposed by Shobe et al. (2017) in order to facilitate comparison. The timestep is set to 10 years. Under detachment-limited conditions (imposed by setting $F_f = 1$), the sediment thickness $H$ equals 0 everywhere and through the entire model run and sediment produced by river incision into bedrock is instantaneously evacuated from the simulated domain. When assuming that water discharge is proportional to the drainage area ($q \propto A^m$) it has been shown that under steady-state
<figure><img src="image_21.png" /></figure>

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

<figure><img src="image_23.png" /><figcaption>**Figure 4.** Sketch of landslide algorithm in three dimensions. Cells shaded in blue indicate the critical nodes where landslides initiate. Cells shaded in red represent the landslide source areas. After mass failing, sediment will be redistributed over the downslope cells using a multiple flow direction algorithm (indicated with green arrows). Sediment deposition rate depends on a transport distance L (cfr. Eqs. 11 and 12).</figcaption></figure>

conditions, fluvial erosion results in the following slope-area relationship (Shobe et al., 2017):

$$S = \left( \frac{U}{K_r A^m} \right)^{1/n} \qquad (13)$$

Fig. 5.b illustrates that when using parameter values listed in Table 1, HyLands reproduces the slope-area relationship given by Eq. 13. Similarly, it can be shown that under transport-limited configurations where $H \gg H_*$, the theoretical slope-area relationship for fluvial incision can be written as (Shobe et al., 2017):

$$S = \left[ \frac{V}{r} + 1 \right]^{1/n} \left[ \frac{U}{K_s} \right]^{1/n} A^{-m/n} \qquad (14)$$

where r represents the runoff rate [L/T]. To mimic a transport limited configuration, we run HyLands assigning an initial sediment thickness H of 100 m. Other parameters values are shown in Table 1. Figure 5.d illustrates the slope-area plot for all nodes of the simulated steady-state landscape showing a close match with the analytical prediction (14). Moreover, HyLands also reproduces the theoretical steady-state sediment flux relationship for transport-limited conditions (Shobe et al., 2017):

$$Q_{S_{fluv}} = UA. \qquad (15)$$

Finally, we evaluate the hybrid nature of the SPACE component in simultaneously simulating fluvial bedrock incision and sediment dynamics. Under such configuration the slope-area relationship can be written as (Shobe et al., 2017):

$$S = \left[ \frac{UV}{K_s A^m r} + \frac{U}{K_r A^m} \right]^{1/n} \qquad (16)$$
<figure><img src="image_24.png" /><figcaption>**Table 1.** Parameter values</figcaption></figure>

<table><thead><tr><th rowspan="2"></th><th colspan="3">No Landsliding</th><th colspan="4">Landsliding</th></tr><tr><th>Detachment-limited</th><th>Transport-limited</th><th>Mixed</th><th>Namche-Barwa</th><th>Pre</th><th>Synthetic(a)<br>LS-Event</th><th>Post</th></tr></thead><tbody><tr><td>Number of rows (–)</td><td>20</td><td>20</td><td>20</td><td>1918</td><td></td><td>75</td><td></td></tr><tr><td>Number of columns (–)</td><td>20</td><td>20</td><td>20</td><td>1149</td><td></td><td>75</td><td></td></tr><tr><td>Node spacing (m)</td><td>100</td><td>100</td><td>100</td><td>20</td><td></td><td>20</td><td></td></tr><tr><td>Time step (yr)</td><td>10</td><td>10</td><td>10</td><td>5</td><td></td><td>5</td><td></td></tr><tr><td>Run time (kyr)</td><td>100</td><td>100</td><td>200</td><td>500</td><td>5×10<sup>6</sup></td><td>100</td><td>5×10<sup>4</sup></td></tr><tr><td>Initial H (m)</td><td>0</td><td>100</td><td>0</td><td>0</td><td>0</td><td>varying</td><td>varying</td></tr><tr><td>U (myr<sup>-1</sup>)</td><td>1×10<sup>-4</sup></td><td>1×10<sup>-4</sup></td><td>1×10<sup>-4</sup></td><td>0</td><td></td><td>1×10<sup>-3</sup></td><td></td></tr><tr><td>K<sub>r</sub> (m<sup>-1</sup>)</td><td>1×10<sup>-3</sup></td><td>1×10<sup>-4</sup></td><td>5×10<sup>-3</sup></td><td>5×10<sup>-4</sup></td><td></td><td>5×10<sup>-5</sup></td><td></td></tr><tr><td>K<sub>s</sub> (m<sup>-1</sup>)</td><td>0.01</td><td>0.01</td><td>0.01</td><td>1×10<sup>-3</sup></td><td></td><td>7.5×10<sup>-5</sup></td><td></td></tr><tr><td>m (–)</td><td>0.5</td><td>0.5</td><td>0.5</td><td>0.5</td><td></td><td>0.5</td><td></td></tr><tr><td>n (–)</td><td>1</td><td>1</td><td>1</td><td>1</td><td></td><td>1</td><td></td></tr><tr><td>H<sub>*</sub> (m)</td><td>1</td><td>1</td><td>1</td><td>2</td><td></td><td>.5</td><td></td></tr><tr><td>φ<sub>sed</sub> (–)</td><td>0</td><td>0</td><td>0</td><td>0</td><td></td><td>0</td><td></td></tr><tr><td>F<sub>ffluv</sub> (–)</td><td>1</td><td>0</td><td>0</td><td>0</td><td></td><td>0</td><td></td></tr><tr><td>V<sup>(b)</sup> (myr<sup>-1</sup>)</td><td>1</td><td>5</td><td>5</td><td>2</td><td></td><td>2</td><td></td></tr><tr><td>V<sup>(b)</sup><sub>Lake</sub> (myr<sup>-1</sup>)</td><td>1</td><td>5</td><td>5</td><td>10</td><td></td><td>10</td><td></td></tr><tr><td>C (kPa)</td><td>-</td><td>-</td><td>-</td><td>15</td><td>-</td><td>15</td><td>-</td></tr><tr><td>φ (°)</td><td>-</td><td>-</td><td>-</td><td>38</td><td>-</td><td>35</td><td>-</td></tr><tr><td>t<sub>LS</sub> (yr)</td><td>-</td><td>-</td><td>-</td><td>2×10<sup>4</sup></td><td>-</td><td>2×10<sup>3</sup></td><td>-</td></tr><tr><td>δ (°)</td><td>-</td><td>-</td><td>-0.01</td><td>-</td><td>0.01</td><td>-</td><td>-</td></tr><tr><td>F<sub>fhill</sub> (–)</td><td>1</td><td>0</td><td>0</td><td>0.25</td><td>-</td><td>0.25</td><td>-</td></tr></tbody></table>

Not all parameters will influence the model outcome in all cases. For example, the value of V is irrelevant for the detachment-limited case when all eroded bedrock passes out of the model domain as permanently suspended fine sediment ($F_{ffluv} = 1$). Landslide parameters are only relevant for models where landslide activity is simulated.

(a): The synthetic landscape evolution model consists of three consequent stages, a pre-landslide stage, a landslide stage and a post-landslide stage. Only parameter values which differ for these stages are listed in the table

(b): HyLands enables spatially variable values for V (Eq. 5) to distinguish between settling velocities over flooded versus non-flooded nodes by changing the values for respectively V and $V_{Lake}$

Under steady-state, both the height of the bedrock and the sediment layer should remain unchanged so that:

$$ \frac{\partial \eta}{\partial t} = \frac{\partial R}{\partial t} + \frac{\partial H}{\partial t} = 0 + 0 = 0 \qquad (17) $$

285 which leads to a constant soil thickness over the landscape given by (Shobe et al., 2017):

$$ H = -H_* \ln \left[ 1 - \frac{V}{\frac{K_s r}{K_r} + V} \right]. \qquad (18) $$
<figure><img src="image_25.png" /></figure>

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

<figure><img src="image_27.png" /><figcaption>**Figure 5. Verification of fluvial SPACE component.** (a) Longitudinal profile of the trunk stream at steady state simulated under detachment-limited conditions where no sediment is present because $F_{ffluv}$ = 1 and all produced sediment is assumed to be evacuated instantaneously. At steady state, the profile evolves towards a concave-upward profile, being in equilibrium with the imposed uplift which is reflected in the steady-state slope-area relationship (b) matching the predicted analytical solution (Eq. 13). (c) Longitudinal profile of the trunk stream at steady state simulated under transport-limited conditions. At steady state, the profile evolves towards a concave-upward profile, being in equilibrium with the imposed uplift which is reflected in the steady-state slope-area relationship (d) matching the predicted analytical solution (Eq. 14). (e) illustrates the steady-state fluvial sediment flux ($Q_{sfluv}$) as a function of the drainage area and matches with the predicted analytical flux-area relationship (Eq. 15). (f) Longitudinal profile of the trunk stream at steady state simulated under mixed alluvial-bedrock conditions. At steady state, both the topographic and bedrock profiles are in equilibrium with the imposed uplift which is reflected in the steady-state slope-area relationship (g) matching the predicted analytical solution (Eq. 16). (h) illustrates the steady-state fluvial sediment flux ($Q_{sfluv}$) as a function of the drainage area and matches with the predicted analytical flux-area relationship under hybrid sediment-bedrock fluvial incision dynamics (Eq. 15).</figcaption></figure>

To evaluate the performance of the hybrid fluvial dynamics, we run HyLands to a steady-state, starting from an initial surface without any sediment cover and using parameter values listed in Table 1. The obtained slope-area relationship (Fig. 2.g) matches with the theoretical relationship Eq. 16, so does the soil thickness *H* which evolves toward a constant thickness (Eq. 18) and the sediment flux relationship honoring (Eq. 15).

290
<figure><img src="image_28.png" /></figure>

## 3.2 Evaluation of the landslide component

Because landslide activity is a stochastic process, it not possible to derive an exact, analytical, solution to evaluate the performance of the landslide component in HyLands. However, it has been shown that most landslide inventories obey consistent magnitude-frequency and magnitude-volume relationships (Malamud and Turcotte, 1999; Stark and Hovius, 2001; Guzzetti et al., 2002; Korup, 2005; Guns and Vanacker, 2014; Larsen and Montgomery, 2012). To evaluate the performance of HyLands, we run the model over a limited amount of time for an area where both relationships are well constrained. The performance of the landslide module is evaluated based on its capacity to reproduce those calibrated relationships.

### 3.2.1 Landslide scaling relationships

A first empirical universal relationship is the landslide magnitude-frequency distribution that describes the number of landslide events of a given size. This relationship is characterized by a negative power law for landslides having an area greater than a given threshold value and a characteristic rollover for smaller landslides (Stark and Hovius, 2001; Malamud and Turcotte, 1999; Guzzetti et al., 2002). Magnitude-frequency distributions are typically described using a three parameter inverse gamma distribution as (Malamud et al., 2004):

$$p(A_L; \rho_l, a_l, s_l) = \frac{1}{a_l * \Gamma(\rho_l)} \left[ \frac{a_l}{A_L - s_l} \right]^{\rho_l + 1} \exp \left[ \frac{a_l}{A_L - s_l} \right] \quad (19)$$

where $A_L$ is the landslide area [$L^2$], $p(A_L)$ is the probability density of a landslide area ($A_L$), $a_l$, $s_l$ and $\rho_l$ are empirical parameters, and $\Gamma(\rho_l)$ is the gamma function of $\rho_l$. A second empirical universal relationship, is the volume-area scaling relationship where the volume $V_L$ of a given landslide is a function of its area $A_L$ as (Hovius et al., 1997):

$$V_L = \alpha_l A_L^{\gamma_l} \quad (20)$$

where $\alpha_l$ is an intercept and $\gamma_l$ a scaling exponent.

### 3.2.2 Applying HyLands to the Namche Barwa-Gyala Peri massif

To evaluate to the performance of HyLands, the model is applied to a digital elevation model (DEM) of the Eastern Himalaya where the Yarlung Tsangpo river cuts through the Namche Barwa-Gyala Peri massif (Fig. 6). The area is characterized by rapid exhumation (King et al., 2016), steep topography, and steep river gradients causing high stream power (Finnegan et al., 2008). To quantify erosion rates in the area, Larsen and Montgomery (2012) mapped more than 15,000 landslides and constructed an inventory of landslides pre-dating 1974 and an inventory containing all landslide events between 1974 and 2007 (Fig 8.a). We use this area solely to demonstrate and evaluate the performance of HyLands and do not aim to reproduce exact features of landscape exhumation in this region.
<figure><img src="image_29.png" /></figure>

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

<figure><img src="image_31.png" /><figcaption>**Figure 6.** Namche Barwa-Gyala Peri massif used for model evaluation. The red rectangle on the inset figure indicates the geographical location of the study area. The green dashed rectangle indicates the part of the DEM used to evaluate HyLands. The shaded colours indicate elevations, which were derived from the 30 m SRTM v3 DEM (NASA JPL, 2013) and resampled to a higher resolution of 20 m using a bicubic interpolation method. Main map is produced with TopoToolbox (Schwanghart and Scherler, 2014). Inset map is made in QGis 3©, using Natural Earth vector and raster map data available at www.naturalearthdata.com.</figcaption></figure>

### 3.2.3 Model parameterization

We run HyLands using the NASA Shuttle Radar Topography Mission (SRTM) v3.0 elevation data as an initial surface (Farr et al., 2007), resampled to a higher resolution of 20 m using a bicubic interpolation method. As shown in Fig. 6, we only simulate part of the larger Namche Barwa-Gyala Peri massif studied by (Larsen and Montgomery, 2012). We selected a region mostly free of glaciers and surrounding the section of the Yarlung Tsangpo river where unit stream power is very high (ranging between 500 - 4000 W.m⁻², (Finnegan et al., 2008). The simulated grid is composed of 1918 × 1149 nodes, covering a total area of ca. 960 km². We simulate landscape evolution over 500 years, using time steps of 5 years. For this experimental run, we assume that there is no uplift. We acknowledge that this condition is not met in the area, but imposing an uplift field is not necessary for evaluating the performance of the HyLands landsliding algorithm. Inserting realistic uplift patterns to simulate the dynamic evolution of the area would require (i) implementation of the complex tectonic configuration of the area (King et al., 2016) and (ii) simulation of a bigger area to capture the dynamic interplay between uplift and river dynamics. This is beyond the scope of the application and we therefore assume that the tectonic configuration controlling landscape evolution over the limited timescale simulated in this experiment (500 years) is captured by the topography of the area (Kirby and Whipple, 2012).

We run HyLands assuming open boundary conditions: all sediment produced within the domain through river incision and landsliding can be exported from the domain across any of the four boundaries. For simplicity, we use a simple stream power
<figure><img src="image_32.png" /></figure>

<figure><img src="image_33.png" /><figcaption>**Table 2.** Simulation movies and script names</figcaption></figure>

<table><thead><tr><th></th><th></th><th>Script name on GitHub</th><th>Link to movie</th></tr></thead><tbody><tr><td rowspan="3">No Landsliding</td><td>Detachment-limited</td><td>HyLands-NoLS-DL.m</td><td>https://doi.org/10.5446/45969</td></tr><tr><td>Transport-limited</td><td>HyLands-NoLS-TL.m</td><td>https://doi.org/10.5446/45967</td></tr><tr><td>Mixed</td><td>HyLands-NoLS-Mixed.m</td><td>https://doi.org/10.5446/45968</td></tr><tr><td rowspan="4">Landsliding</td><td>Real DEM, Namche-Barwa</td><td>HyLands-LS-NB.m</td><td>https://doi.org/10.5446/45973</td></tr><tr><td>Before intense LS period</td><td>HyLands-LS-B-LS.m</td><td>https://doi.org/10.5446/45970</td></tr><tr><td>Synthetic intense LS period</td><td>HyLands-LS-LS.m</td><td>https://doi.org/10.5446/45971</td></tr><tr><td>After LS period</td><td>HyLands-LS-A-LS.m</td><td>https://doi.org/10.5446/45972</td></tr></tbody></table>

formulation for river incision where thresholds for both sediment entrainment and bedrock erosion are negligible. Standard scaling exponents are used ($m = 0.5$ and $n = 1$ in Eq. 4) and the bed sediment porosity and the fraction of fine river sediments are assumed to be zero ($\phi_{sed} = 0$ and $F_{f_{ftuv}} = 0$ in Eq. 1). We calculate landslide activity using the landslide module of HyLands and assume that 25% of the landslide-derived sediment is evacuated out of the system as fine material ($F_{f_{hill}} = 0.25$). We assume that the angle of internal friction ($\phi$) is comparable to the mode of the topographical slope distribution (Burbank et al., 1996; Korup, 2008; Montgomery and Gran, 2001), reported to range between 37°- 39° and here set to 38° (Larsen and Montgomery, 2012). Cohesion $C$ is known to vary over a wide range and strongly depends on rock mechanical properties (Wyllie and Mah, 2017). Site-specific calibration would require detailed mapping of lithological units and we therefore set $C$ to 15 kPa, a value in the range of previously optimized cohesion for the Himalaya (Jeandet et al., 2019, 12-20Pa in ). Cohesion and the angle of internal friction influences the size-distribution of landslides in several ways (Jeandet et al., 2019). The angle of internal friction $\phi$ controls the angle of the potential rupture plane such that lower values of $\phi$ will result in lower rupture dipping angles, which, for the same topographical configuration, results in thicker and larger landslides. The effective rock cohesion value $C$ influences the critical hillslope height $H_c$ in Eq. 6. Larger values for $C$ will result in larger values for $H_c$, thus decreasing the probability of landslides on less steep slope sections and resulting in fewer small landslides. The minimum value of the spreading slope under which landslide-derived sediment is redistributed on hillslopes is set to 0.01°. We set the return time for landsliding $t_{LS}$ to $2 \times 10^4$ years. $t_{LS}$ regulates the probability that unstable cells evolve into a landslide (Eq. 9) and therefore controls the number of landslide events per timestep $\Delta t$. When applying HyLands to reconstruct or predict landslide activity, $t_{LS}$ should be a function of the frequency (or return time) of triggering events (large earthquakes or rainfall events). A full overview of the model parameters is given in Table 1.

### 3.2.4 Model evaluation results

Figure 7 shows time slices of the model run after 5 (initial iteration), 165, 330 and 500 (final iteration) model years. Locations for landslide initiation (critical nodes) are well spread over the landscape. The number of landslide events (the number of
<figure><img src="image_34.png" /></figure>

<figure><img src="image_35.png" /><figcaption>**Figure 7.** Timeslices of HyLands model run for the Namche Barwa region after 5, 165, 330 and 500 model years. (a, d, g, j) Indicate the location of the landslides at the given timestep (black dots). The colors represent the square root (used for ease of visualization) of the landslide erosion (-) and deposition (+) during the given time step. (b, e, h, k) Zoom of simulated landslide erosion/deposition indicated with the red square in respectively (a, d, g and j). (c, f, i, l) Represent the square root of the total amount of sediment, generated through river incision (SPACE module) and landsliding after respectively 5, 500, 1500 and 2000 model years. The grey shaded colours indicate elevations, which were derived from the 30 m SRTM v3 DEM (NASA JPL, 2013).</figcaption></figure>
<figure><img src="image_36.png" /></figure>

black dots in Fig. 7) is mainly controlled by the return time for landsliding $t_{LS}$. The fan-shaped deposition zones of landslide-derived sediments reflect the use of a multiple flow sediment routing algorithm. Landslide-derived sediment predominantly accumulates at hillslope toes, as well as in or near river channels. Some accumulation also occurs on hillslopes. Overall, landslides strongly influence the thickness of the alluvial bed sediment layer. Note that the shape of the erosion and deposition zones adjust through the course of the model run. The presence of previous landslide activity alters the topographic relief and hence determines which pixels become susceptible to erosion and deposition as landscape evolution continues (Fig. 7).

360

To quantitatively evaluate the performance of the HyLands landslide module, we compare modelled landslide properties against observed scaling relationships (Eq. 19 and 20). Figure 8.a compares the modelled and observed magnitude-frequency distribution. We observe good correspondence between the model and the data, with the power-law tail of the distribution falling within the envelope defined by the two inventories of Larsen and Montgomery (2012). Similarly to observed magnitude-frequency distributions, HyLands simulates the rollover or the transition from an increasing magnitude-frequency relationship to a decreasing one. This observation confirms that the shape of landslide magnitude-frequency distributions can be explained using mechanical landslide processes (cfr. 2.2.1) and the geometry of the studied region (Jeandet et al., 2019). Figure 8.b shows that HyLands is capable of approaching the universal Area-Volume relationships found by Larsen et al. (2010). While HyLands seems to overestimate simulated landslide volumes for very small landslides, the fit between HyLands and the observed relationship improves for large landslides. We attribute overestimation of landslide volumes for smaller landslides to the nature of the landslide algorithm: while HyLands does not allow any sediment to be deposited within the landslide scar, this typically does occur in nature. Future developments of the algorithm could allow in-scar deposition for more realistic simulations. Moreover, there is a resolution effect: due to DEM noise or heterogeneity, the algorithm might select small landslides of one or two cells on very steep hillslope patches thus resulting in high landslide volumes. However, such steep hillslope patches might represent noise in the DEM rather than actual steep slopes. The use of high resolution DEMs could partly resolve this issue.

365

370

375

## 3.3 Model application

The explicit coupling of landslides and landslide-derived sediment to long-term landscape evolution enables the study of a wide range of interactions which otherwise can only be inferred or partially simulated. A common application is to evaluate the coupling between landslides and river-bed morphology. To evaluate the impact of a landslide event on long-term channel profile evolution, we run a synthetic landscape evolution model to steady state. After 5 million years, we simulate a period of 100 years with intense landslide activity, analogous to a period of elevated landslide activity triggered by a series of seismic events. After 100 years of landslide activity, we assume landslides are no longer triggered and let the landscape evolve back to its original steady state. Such an experiment not only allows evaluation of the extent to which landslides perturb the topography of river profiles, but also enables estimates of the time required for a landscape to respond to a major perturbation (e.g. a series of earthquake-triggered landslides).

380

385

The model run consists of three stages. In the first stage, the model is run to a steady state. For reasons of comparability, we simulate landscape evolution on a grid similar to the one used for the verification runs (section 3), i.e. with a single open boundary node in the southwestern corner. To simulate more realistic landscape scales, we use a domain of 75 by 75 nodes
<figure><img src="image_37.png" /></figure>

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

<figure><img src="image_39.png" /><figcaption>**Figure 8.** Comparison between modelled and observed characteristic landslide scaling relationships. (a) Magnitude frequency relationship. The grey and black line represent the best fitting inverse gamma distribution (Eq. 19) of the landslide activity mapped before 1974 and between 1974 - 2012, respectively (Larsen and Montgomery, 2012). The fitting parameters are respectively ($a_l = 768 \times 10^{-6}$, $s_l = -32.6 \times 10^{-6}$ and $\rho_l = 1.27$) and ($a_l = 6100 \times 10^{-6}$, $s_l = -311 \times 10^{-6}$ and $\rho_l = 0.96$)). The red dots represent the magnitude-frequency distribution simulated using HyLands. (b) Area-Volume relationship. The grey dashed zone represents the expected area-volume scaling relationship, observed for bedrock landslides in the Himalaya (Larsen et al., 2010). Fit is calculated using Eq. 20 with fitting parameters $\gamma = 1.32 \pm 0.02$ and $^{10}\log\alpha = -0.49 \pm 0.06$. The red dots represent the geometry of the data simulated with HyLands.</figcaption></figure>

with a higher resolution of 20 m. A complete overview of model parameter values is given in Table 1. The evolution of the landscape over time is shown in a series of time slices (Fig. 9.a-c). Model behaviour at this stage is as expected for the SPACE river erosion model when hillslope processes are not explicitly simulated (Shobe et al., 2017). During the first timesteps, the
<figure><img src="image_40.png" /></figure>

To steady state

<figure><img src="image_41.png" /><figcaption>(b) HYLANDS: 50000 model years</figcaption></figure>

<figure><img src="image_42.png" /><figcaption>(c) HYLANDS: 5000000 model years</figcaption></figure>

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

100 years of LS events

<figure><img src="image_44.png" /><figcaption>(e) HYLANDS: 50 model years</figcaption></figure>

<figure><img src="image_45.png" /><figcaption>(f) HYLANDS: 100 model years</figcaption></figure>

<figure><img src="image_46.png" /><figcaption>**Figure 9.** Synthetic model run showing landscape evolution to steady state followed by an intense landsliding period of 100 years.</figcaption></figure>

(a-d) Time slices showing evolution of the landscape to steady state, before the landslide period. (e-h) Time slices showing the landslide period where intense landsliding is occurring over a period of 100 years. Note that, during landsliding, both pure landslide dams arise as well as irregularities in the bedrock profile (the grey bumps). The latter originate from the river being redirected after landsliding forming epigenetic gorges (see text).
<figure><img src="image_47.png" /></figure>

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

drainage network establishes and the landscape gradually approaches a steady state with uniform sediment thickness across the entire landscape.

In the second stage (Fig. 9.d-f), we simulate a period of intense landslide activity by triggering a large number of landslides. Landslides are initiated based on their probability of sliding (Eq. 7) assuming an internal friction angle of 35° and a low landslide return time ($t_{LS} = 2 \times 10^3$ years). Under this configuration, many of the steep portions of the landscape become prone to landslide erosion and transform into a landslide source area. We assume that 25% of the landslide-derived sediment is instantaneously evacuated out of the system as fine material ($F_{hill} = 0.25$). Landslides trigger the formation of landslide dams, resulting in flooded river sections. Landslide dams not only alter the topographic elevation of the simulated domain but also change the drainage network. The location of the river bed can change due to landslides and landslide-derived sediment rearranging the valley-bottom topography. This is why the bedrock profile of the plots shown in Fig. 9.e-h has a bumpy shape at several locations along the profile.

Immediately after the intense landsliding period, the trunk stream of the drainage network is choked with sediment and landslide dams are abundant (Fig. 10.a-f). In the first few thousand years following the intense landsliding period, the lakes gradually fill in with sediment. After 1500 years, most of the landslide-dammed lakes are filled with sediment. The fluvial profile is now characterized by a chain of knickpoints characteristic for fluvial profiles experiencing the delivery of immobile debris by landslides (Ouimet et al., 2007) or other hillslope processes (Shobe et al., 2016, 2018). Where the in-channel bedrock is not covered with sediment, river incision into bedrock continues. However, upstream of landslide dams, the bed is choked with sediment and the alluvial cover is too thick for bedrock incision to continue. As bedrock uplift continues and sediment is slowly being evacuated from filled lakes, the bedrock profile adjusts and small knickpoints are created along the river profile. While the specific cause is different (landslide dams ponding sediment vs. delivery of large-grained colluvium), the mechanism of knickpoint generation is similar to the numerical simulations of Shobe et al. (2016) and Shobe et al. (2018) in that a bare-bedrock reach downstream of a sediment-mantled reach can undergo faster erosion, thereby generating knickpoints that are decoupled from the baselevel signal.

There are two distinct mechanisms for the generation of irregularities in the channel profiles: drainage re-routing due to landslide dams and knickpoint generation due to spatially varying sediment cover triggering differential erosion. The former mechanism can result in reaches where the bedrock slope is adverse relative to the water surface slope (the bumps in the bedrock profile in Fig. s 9 and 10). The latter creates variability in the magnitude, but not the direction, of bedrock slope. The drainage re-routing mechanism dominates in the simulations presented here (Fig. 10).

# 4 Discussion

Landscapes are the outcome of external perturbations, such as climate or tectonic variability, and internal dynamics originating from the coupling between fluvial incision and hillslope response (Burbank and Anderson, 2011; Glade et al., 2019). Much effort has been devoted towards understanding the relationship between fluvial erosion efficiency and climate variability both through theoretical developments (Tucker, 2004; Lague, 2014) and observations (DiBiase and Whipple, 2011; Ferrier et al.,
<figure><img src="image_49.png" /></figure>

Post landslide events

<figure><img src="image_50.png" /><figcaption>(d) HYLANDS: 3500 model years</figcaption></figure>

<figure><img src="image_51.png" /><figcaption>(b) HYLANDS: 1000 model years</figcaption></figure>

<figure><img src="image_52.png" /><figcaption>(e) HYLANDS: 6500 model years</figcaption></figure>

<figure><img src="image_53.png" /><figcaption>(c) HYLANDS: 1500 model years</figcaption></figure>

<figure><img src="image_54.png" /><figcaption>(f) HYLANDS: 500000 model years</figcaption></figure>

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

**Figure 10.** Synthetic model run showing recovery from intense landsliding illustrated in Fig. 9. (a-h) Time slices showing reestablishment of the landscape steady state. Bedrock bumps created by landslide-induced drainage redirection are eroded and the channel re-attains its smoothly concave-up, steady-state configuration.
<figure><img src="image_56.png" /></figure>

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

2013). A main finding of those authors is that the role of allogenic fluvial response (i.e. transient adjustment to an external perturbation) can only be understood when considering autogenic fluvial dynamics such as the existence of incision thresholds (Snyder et al., 2003; Lague et al., 2005) and the internal lithological heterogeneity in a landscape (Campforts et al., 2019; Glade et al., 2019). However, the role of landslides in long-term landscape evolution, and especially the dynamic interaction between river incision and landslide activity, is only poorly understood. HyLands offers a tool to study dynamic feedbacks between landslides and river incision. The role of sediment dynamics in altering fluvial erosion and sediment transport is clearly illustrated in the numerical experiment (Fig. 9) where 5-10 kyr are required for the landscape to evolve back to a steady state after a pulse of landsliding. Not only does the delivery of landslide-derived sediment to the channel bed alter the topography of the channel profile, it also results in the formation of bedrock knickpoints and associated retreating incision waves. HyLands thereby corroborate earlier observations that landslides, and tight channel-hillslope couplings in general, are autogenic mechanisms altering the way in which landscapes respond to external (allogenic) perturbations (Ouimet et al., 2007; Shobe et al., 2016; Glade et al., 2019). HyLands is designed to study the dynamic feedbacks between landslides and river erosion at large spatial and temporal scales. To do so, the model integrates an algorithm for deep-seated landsliding with a recently proposed model for fluvial incision (Shobe et al., 2017). HyLands enables simulations over several millions of years and reproduces analytical predictions for fluvial dynamics and observed landslide scaling relationships.

## 4.1 Fluvial component

The SPACE river erosion model, which governs river evolution in HyLands, advances on existing river incision models in that it explicitly simulates the role of sediment in reducing the efficiency of bedrock incision (Beaumont et al., 1992; Lague, 2010). However, like SPACE, HyLands does not simulate the effect of increased bedrock incision efficiency due to mobile sediment acting as eroding tools (Sklar and Dietrich, 2004). Field observations warrant consideration of the tool effect (Cook et al., 2013; Beer et al., 2017), and theoretical predictions have shown that the interaction between sediment and bedrock incision is adjusted when the tool effect is considered (Gasparini et al., 2007). The impact of explicitly simulating the tool effect due to landslide-derived sediment has been evaluated in a numerical modelling study (Egholm et al., 2013). Egholm et al. (2013) concluded that landslide activity and its delivery of abrasive agents to the channel accelerate fluvial incision in actively uplifting mountain regions, whereas the lack of landslides in tectonically inactive mountain ranges strongly decreases erosion efficiency and enables topographic preservation. Adding the tool effect and evaluating its potential importance is therefore a primary goal for further model development in HyLands. Simulating the tool effect of sediments can be achieved by making the bedrock erosion function dependent (Eq. 4) on the sediment flux $Q_s$ (cfr. e.g. Gasparini et al., 2007; Hobley et al., 2011).

Like SPACE, HyLands does not include process-based approaches (Kean and Smith, 2004; Wobus et al., 2006; Davy and Lague, 2009; Coulthard et al., 2013), simplified width adjustment rules (Lague, 2010; Yanites, 2018), or empirical closures (Attal et al., 2008) to dynamically calculate river width adjustments though time. HyLands assumes a relationship between drainage area and river width depending on the scaling exponent $m$ (fixed in our simulations to 0.5; see also Table 1). This approach implies that all river cells in the landscape are assumed to occupy 1 grid cell with distance $dx$, that channel width may be less than, equal to, or greater than $dx$, and that river width is only a function of contributing drainage area. It has
<figure><img src="image_58.png" /></figure>

however been shown that river width might vary as a function of sediment flux or under varying tectonic configurations (Amos and Burbank, 2007; Turowski et al., 2009). Recent work using a 2D hydro-sedimentary numerical model (Davy et al., 2017) based on the Saint-Venant Equations has shown that river re-organisation and narrowing after landslide events might strongly increase sediment transport capacity and alter sediment evacuation time after big landslide events (Croissant et al., 2017). While simulating dynamic river width reorganisation at the landscape scale is currently not possible over longer timescales due to computational limitations, generic approximations for the landslide triggered channel narrowing (Croissant et al., 2019) could be integrated in future versions of HyLands.

## 4.2 Landslides

The landslide algorithm in HyLands is based on finite slope mechanics and assumes a planar rupture plane geometry. Although our approach reproduces observed magnitude-frequency and area-volume scaling relationships and is supported by previous work where landslides have been simulated using planar rupture planes (Jeandet et al., 2019), the use of more advanced rupture plane geometries has been proposed. Gallen and Wegmann (2017) for example propose the use of concave-upward rupture planes to simulate co-seismic landsliding. However, their approach is based on the statistical aggregation of one-dimensional slope-stability solutions and therefore does not fully honor the three-dimensional topographic complexity of the topographic surface like HyLands does. Evaluating the role of varying rupture plane geometries in three dimensions is one of the potential future developments of HyLands.

At this stage HyLands does not explicitly simulate shallow landsliding which typically occurs at the interface between the bedrock and the overlying sediment/regolith cover. Given the existing ability of HyLands to simultaneously simulate bedrock evolution and sediment thickness, adding a shallow landslide algorithm is feasible (Montgomery and Dietrich, 1994; Claessens et al., 2007) and would further our understanding of the coupling between climate variability and landscape stability (Parker et al., 2016). For shallow landslides to be added as a component in HyLands, the implementation and calibration of a regolith formation and soil flux model will however be required (Campforts et al., 2016). Using additional model components to simulate soil formation and transport requires calibration of several additional processes components; care is needed to prevent over-calibration and -parameterization of the model (Van Rompaey and Govers, 2002).

Our probabilistic sliding mechanism neglects seismic or hydrological landslide triggers (e.g. Keefer, 1984, 2002; Keefer and Larsen, 2007; Marc et al., 2015, 2018, 2019). Rather we simulate landslides as a stochastic process based on the mechanical stability of slope patches. Future developments could however adjust the spatial probability of landsliding by coupling an explicit earthquake model to HyLands (cfr., Croissant et al., 2019). The probability for co-seismic landslide activity can then be directly obtained by using constrained relationships between Peak Ground Acceleration (PGA) and landslide initiation probabilities (Meunier et al., 2007).

HyLands uses a non-linear, multiple-flow sediment redistribution scheme depending on topographic slope. Accurately simulating landslide sediment run-out distances is however a challenging process which is difficult to constrain and often simulated using empirical approximations (e.g. based on the absolute height difference within a landslide (Claessens et al., 2007)). A potential way to validate and calibrate landslide runout distance would be to compare landslide-derived sediment distributions
<figure><img src="image_59.png" /></figure>

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

simulated with HyLands with runout distances simulated with higher complexity models (Iverson, 2000; Iverson and George, 2016; Zhou et al., 2020). Nevertheless, in its current form, HyLands reproduces characteristics of landscapes and channel profiles dominated by deep-seated landslides (Ouimet et al., 2007), and is therefore a useful tool to study the interaction between river incision and landslide dynamics at landscape evolution space and time scales.

495

## 4.3 Calibration of HyLands

A main challenge when applying HyLands in real settings is the calibration of both the river incision and landslide parameters. The power of any LEM lies in its capacity to integrate data over multiple spatial and temporal timescales. Therefore, a range of datasets can be used to constrain model parameters. We identify three main categories of potential calibration data.

500

1. Topographic parameters that can be derived from DEMs. These include a range of metrics describing river characteristics (drainage density, river steepness, river stream power), hillslope properties (slope distribution, mean and median slope angles, aspect), and landslide scaling relationships (magnitude-frequency and area-volume distributions). All these metrics can be derived from topographic data and subsequently used to constrain HyLands erosion parameters (cfr. Fig. 7).

505

2. Data directly constraining the integrated effects of river incision or landslide erosion. Ongoing efforts to map mass movements in landslide prone areas now enable estimates of erosion rates over decal timescales (Hovius et al., 1997; Larsen and Montgomery, 2012). Data on sediment redistribution following landslide events is however more difficult to collect. Although initial compilations now exist of global landslide sediment mobilization rates (Broeckx et al., 2020), such inventories remain incomplete. Data on landslide mobilization rates can be used to train HyLands while in turn, HyLands can be used to further extend datasets on landslide mobilisation rates and predict landslide sediment production rates in regions which are otherwise difficult to access.

510

3. Catchment-averaged cosmogenic radionuclide (CRN) derived erosion rates. CRN data has been used to calibrate river incision models that explicitly integrate the stochastic nature of fluvial incision over time (DiBiase and Whipple, 2011; Scherler et al., 2017; Campforts et al., 2019). However, CRN data is sensitive to landslide activity (Niemi et al., 2005; Yanites et al., 2009; Wang et al., 2017) and the calibration of stochastic river incision models has been shown to be sensitive to landslide activity (Campforts et al., 2019). HyLands directly simulates stochastic landsliding and hence enables explicit simulation of the impact of landsliding on CRN-derived erosion rates, making them a promising tool to constrain HyLands.

515

520

## 4.4 Potential applications

Given the capacity of HyLands to explicitly simulate the interaction between fluvial dynamics and landslide triggering, it provides a unique toolbox to advance the field of geomorphology on several fronts. In the following we outline two example applications.
<figure><img src="image_61.png" /></figure>

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

First, HyLands can be used as an experimental environment to test how landslides influence landscape response to external perturbations. Landslides are known to mediate long term landscape evolution (Korup, 2005; Korup et al., 2007). By altering sediment fluxes, they fundamentally alter the dynamic equilibrium between hillslopes and rivers, resulting in long-term implications for landscape evolution (Egholm et al., 2013). Moreover, large landslides are reported to critically alter drainage networks by causing major river captures (Korup et al., 2007; Dahlquist et al., 2018).

Second, HyLands can be used to evaluate the response time of a landscape to a landslide event and to understand the timescales over which landslide-derived sediments are exported from the landscape (Wang et al., 2015; Li et al., 2016; Yanites et al., 2010). As illustrated in this paper, a landscape requires a certain response time to recover from a landslide event—or a series of landslide events—and to evolve back to a steady-state configuration (Fig. 10). Landslide activity is typically manifested in downstream sediment dynamics and LEMs are the right tool to simulate large-scale landscape response to landslide activity in upstream mountain regions. HyLands will enable prediction of downstream sediment response to landsliding, provided that the model can be calibrated effectively (section 4.3).

# 5 Conclusions

We presented a new, fully coupled model for river incision into bedrock, sediment transport, and bedrock landsliding. HyLands couples a mass conservative sediment-flux-dependent incision model (SPACE, Shobe et al., 2017) with a deep-seated landslide algorithm (Densmore et al., 1998) and a multiple-flow sediment redistribution algorithm (Carretier et al., 2018). HyLands is designed to simulate landscape evolution at large temporal and spatial scales. The fluvial component of the model matches known, steady-state analytical solutions developed in earlier work (Davy and Lague, 2009; Shobe et al., 2017). Landslides produced by HyLands replicate observed scaling relationships indicating the realism of the simulations. HyLands is implemented in the TopoToolbox GIS interface (Schwanghart and Kuhn, 2010), thereby facilitating the use of rasterized field data for calibration and providing direct access to a wide range of GIS analysis tools. In an example application, we illustrated how HyLands can be used to evaluate the impact of landslide activity on fluvial and hillslope characteristics. We showed how landslide activity triggers the formation of landslide-dammed lakes and how HyLands is capable of simulating subsequent lake infilling and knickpoint formation, similar to reported landscape changes following landslide activity (Ouimet et al., 2007). The foremost advantage of HyLands is its capacity to explicitly simulate the role of landslides, landslide-derived sediment and fluvial dynamics at the landscape scale. The model is well-suited to address a range of new questions related to how channel-hillslope coupling modulates landscape response to external perturbations.

**Code availability.** The HyLands 1.0 TTLEM component as well as all other TTLEM components used in this paper are part of TopoToolbox version 2. The source code and future updates are available in the GIT repository: https://github.com/BCampforts/topobox. The exact version of the software code used to produce the results presented in this paper is archived on Zenodo (https://zenodo.org/badge/latestdoi/78645261). Upon publication we will submit a pull request to the master repository where TopoToolbox is housed: (wschwanghart/topobox). Documentation, installation instructions, and software dependencies for the HyLands project can be found at https://github.
<figure><img src="image_63.png" /></figure>

com/BCampforts/topobox. Detailed scripts and user manuals for the simulations illustrated in this paper can be found at https://github.com/BCampforts/pub_hylands_campforts_etal_GMD which is also archived on Zenodo: https://zenodo.org/badge/latestdoi/247779084. HyLands is platform independent and requires MATLAB 2014b or higher and the Image Processing Toolbox. The HyLands modeling framework is distributed under a MIT open-source license.

560

*Data availability.* Digital elevation models used in this paper are derived from the 30 m SRTM v3 DEM (NASA JPL, 2013). The resampled DEM is available through https://github.com/BCampforts/topobox

*Video supplement.* Videos are described in Table 2 of the main text which contains hyperlinks to the following movies: https://doi.org/10.5446/45969; https://doi.org/10.5446/45967; https://doi.org/10.5446/45968; https://doi.org/10.5446/45973; https://doi.org/10.5446/45970; https://doi.org/10.5446/45971; https://doi.org/10.5446/45972

565

*Author contributions.* BC and CMS conceived the conceptual model on landslide-derived sediment dynamics and designed the research project. BC developed the algorithm with help from CMS. BC and CMS took the lead in writing the paper. Concepts to verify and evaluate model components were conceptualized by PS, DL, MVM and JB. All authors contributed to shaping the research, analyses, and paper.

*Competing interests.* The authors declare that they have no conflict of interest.

*Acknowledgements.* BC received funding from the Research Foundation Flanders, FWO grant agreement no. 12Z6518N. CMS received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement no. 833132. PS and DL have received funding from the European Research Council (ERC) under the European Union's Horizon Horizon 2020 research and innovation programme (grant agreement no. 803721).

570
<figure><img src="image_64.png" /></figure>

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

References

575 Amos, C. B. and Burbank, D. W.: Channel width response to differential uplift, Journal of Geophysical Research, 112, F02010, https://doi.org/10.1029/2006JF000672, http://doi.wiley.com/10.1029/2006JF000672, 2007.

Andrews, D. J. and Hanks, T. C.: Scarp degraded by linear diffusion: Inverse solution for age, Journal of Geophysical Research, 90, 10193, https://doi.org/10.1029/JB090iB12p10193, http://doi.wiley.com/10.1029/JB090iB12p10193, 1985.

580 Armitage, J. J., Whittaker, A. C., Zakari, M., and Campforts, B.: Numerical modelling of landscape and sediment flux response to precipitation rate change, Earth Surface Dynamics, 6, 77–99, https://doi.org/10.5194/esurf-6-77-2018, https://www.earth-surf-dynam-discuss.net/esurf-2017-34/https://www.earth-surf-dynam.net/6/77/2018/, 2018.

Attal, M., Tucker, G. E., Whittaker, A. C., Cowie, P. A., and Roberts, G. P.: Modelling fluvial incision and transient landscape evolution: Influence of dynamic Channel adjustment, Journal of Geophysical Research: Earth Surface, 113, 1–16, https://doi.org/10.1029/2007JF000893, 2008.

585 Baum, R. L., Godt, J. W., and Savage, W. Z.: Estimating the timing and location of shallow rainfall-induced landslides using a model for transient, unsaturated infiltration, Journal of Geophysical Research, 115, F03013, https://doi.org/10.1029/2009JF001321, http://doi.wiley.com/10.1029/2009JF001321, 2010.

Beaumont, C., Fullsack, P., and Hamilton, J.: Erosional control of active compressional orogens, in: Thrust Tectonics, pp. 1–18, Springer Netherlands, Dordrecht, https://doi.org/10.1007/978-94-011-3066-0_1, http://link.springer.com/10.1007/978-94-011-3066-0{\_}1, 1992.

590 Beer, A. R., Turowski, J. M., and Kirchner, J. W.: Spatial patterns of erosion in a bedrock gorge, Journal of Geophysical Research: Earth Surface, 122, 191–214, https://doi.org/10.1002/2016JF003850, http://doi.wiley.com/10.1002/2016JF003850, 2017.

Beven, K. and Freer, J.: Equifinality, data assimilation, and uncertainty estimation in mechanistic modelling of complex environmental systems using the GLUE methodology, Journal of Hydrology, 249, 11–29, https://doi.org/10.1016/S0022-1694(01)00421-8, http://dx.doi.org/10.1016/S0022-1694(01)00421-8http://linkinghub.elsevier.com/retrieve/pii/S0022169401004218, 2001.

595 Broeckx, J., Rossi, M., Lijnen, K., Campforts, B., Poesen, J., and Vanmaercke, M.: Landslide mobilization rates: A global analysis and model, Earth-Science Reviews, 201, 102972, https://doi.org/10.1016/j.earscirev.2019.102972, https://linkinghub.elsevier.com/retrieve/pii/S0012825219302557, 2020.

Burbank, D., Meigs, A., and Brozović, N.: Interactions of growing folds and coeval depositional systems, Basin Research, 8, 199–223, https://doi.org/10.1046/j.1365-2117.1996.00181.x, 1996.

600 Burbank, D. W.: Rates of erosion and their implications for exhumation, Mineralogical Magazine, 66, 25–52, https://doi.org/10.1180/0026461026610014, https://www.cambridge.org/core/product/identifier/S0026461X0004500X/type/journal/{\_}article, 2002.

Burbank, D. W. and Anderson, R. S.: Tectonic Geomorphology, John Wiley & Sons, Ltd, Chichester, UK, https://doi.org/10.1002/9781444345063, http://books.google.com/books?id=83FuAvtSwE4C{\}&pg=PT270http://doi.wiley.com/10.1002/9781444345063, 2011.

605

Campforts, B. and Govers, G.: Keeping the edge: A numerical method that avoids knickpoint smearing when solving the stream power law, Journal of Geophysical Research: Earth Surface, 120, 1189–1205, https://doi.org/10.1002/2014JF003376, http://doi.wiley.com/10.1002/2014JF003376, 2015.

Campforts, B., Vanacker, V., Vanderborght, J., Baken, S., Smolders, E., and Govers, G.: Simulating the mobility of meteoric 10 Be in the landscape through a coupled soil-hillslope model (Be2D), Earth and Planetary Science Letters, 439,

610
<figure><img src="image_66.png" /></figure>

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

143–157, https://doi.org/10.1016/j.epsl.2016.01.017, http://dx.doi.org/10.1016/j.epsl.2016.01.017 http://linkinghub.elsevier.com/retrieve/pii/S0012821X16000339, 2016.

Campforts, B., Schwanghart, W., and Govers, G.: Accurate simulation of transient landscape evolution by eliminating numerical diffusion: the TTLEM 1.0 model, Earth Surface Dynamics, 5, 47–66, https://doi.org/10.5194/esurf-5-47-2017, http://www.earth-surf-dynam.net/5/47/2017/, 2017.

Campforts, B., Vanacker, V., Herman, F., Vanmaercke, M., Schwanghart, W., Tenorio, G. E., Willems, P., and Govers, G.: Lithology and orographic precipitation control river incision in the tropical Andes, Earth Surface Dynamics Discussions, 2019, 1–43, https://doi.org/10.5194/esurf-2019-48, 2019.

Carretier, S., Martinod, P., Reich, M., and Godderis, Y.: Modelling sediment clasts transport during landscape evolution, Earth Surface Dynamics, 4, 237–251, https://doi.org/10.5194/esurf-4-237-2016, http://www.earth-surf-dynam.net/4/237/2016/, 2016.

Carretier, S., Tolorza, V., Regard, V., Aguilar, G., Bermúdez, M. A., Martinod, J., Guyot, J. L., Hérail, G., and Riquelme, R.: Review of erosion dynamics along the major N-S climatic gradient in Chile and perspectives, Geomorphology, 300, 45–68, https://doi.org/10.1016/j geomorph.2017.10.016, https://doi.org/10.1016/j.geomorph.2017.10.016, 2018.

Champel, B.: Growth and lateral propagation of fault-related folds in the Siwaliks of western Nepal: Rates, mechanisms, and geomorphic sig-nature, Journal of Geophysical Research, 107, 2111, https://doi.org/10.1029/2001JB000578, http://doi.wiley.com/10.1029/2001JB000578, 2002.

Claessens, L., Schoorl, J., and Veldkamp, A.: Modelling the location of shallow landslides and their effects on landscape dynamics in large watersheds: An application for Northern New Zealand, Geomorphology, 87, 16–27, https://doi.org/10.1016/j.geomorph.2006.06.039, http://linkinghub.elsevier.com/retrieve/pii/S0169555X06003382, 2007.

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 Surface Processes and Landforms, 38, 683–695, https://doi.org/10.1002/esp.3313, http://doi.wiley.com/10.1002/esp.3313, 2013.

Coulthard, T. J., Neal, J. C., Bates, P. D., Ramirez, J., de Almeida, G. A. M., and Hancock, G. R.: Integrating the LISFLOOD-FP 2D hydrodynamic model with the CAESAR model: implications for modelling landscape evolution, Earth Surface Processes and Landforms, 38, 1897–1906, https://doi.org/10.1002/esp.3478, http://doi.wiley.com/10.1002/esp.3478, 2013.

Croissant, T., Lague, D., Steer, P., and Davy, P.: Rapid post-seismic landslide evacuation boosted by dynamic river width, Nature Geoscience, 10, 680–684, https://doi.org/10.1038/ngeo3005, http://www.nature.com/articles/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, https://doi.org/10.1016/j.geomorph.2019.04.017, https://doi.org/10.1016/j.geomorph.2019.04.017, 2019.

Culmann, K.: Die Graphische Statike, Verlag von Meyer & Zeller, Zurich, https://archive.org/details/diegraphischest01culmgoog/page/n8/mode/2up, 1875.

Dahlquist, M. P., West, A. J., and Li, G.: Landslide-driven drainage divide migration, Geology, 46, 403–406, https://doi.org/10.1130/G39916.1, https://pubs.geoscienceworld.org/gsa/geology/article/46/5/403/528890/Landslidedriven-drainage-divide-migration, 2018.

Davy, P. and Lague, D.: Fluvial erosion/transport equation of landscape evolution models revisited, Journal of Geophysical Research: Solid Earth, 114, 1–16, https://doi.org/10.1029/2008JF001146, 2009.
<figure><img src="image_68.png" /></figure>

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

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, Journal of Geophysical Research: Earth Surface, 122, 1491–1512, https://doi.org/10.1002/2016JF004156, http://doi.wiley.com/10.1002/2016JF004156, 2017.

650

Densmore, A. L., Ellis, M. A., and Anderson, R. S.: Landsliding and the evolution of normal-fault-bounded mountains, Journal of Geophysical Research: Solid Earth, 103, 15 203–15 219, https://doi.org/10.1029/98JB00510, http://doi.wiley.com/10.1029/98JB00510, 1998.

DiBiase, R. A. and Whipple, K. X.: The influence of erosion thresholds and runoff variability on the relationships among topography, climate, and erosion rate, Journal of Geophysical Research, 116, F04036, https://doi.org/10.1029/2011JF002095, http://doi.wiley.com/10.1029/2011JF002095, 2011.

655

Dussauge, C., Grasso, J.-R., and Helmstetter, A.: Statistical analysis of rockfall volume distributions: Implications for rockfall dynamics, Journal of Geophysical Research: Solid Earth, 108, https://doi.org/10.1029/2001JB000650, http://doi.wiley.com/10.1029/2001JB000650, 2003.

Egholm, D. L., Knudsen, M. F., and Sandiford, M.: Lifespan of mountain ranges scaled by feedbacks between landsliding and erosion by rivers, Nature, 498, 475–478, https://doi.org/10.1038/nature12218, http://www.nature.com/articles/nature12218, 2013.

660

Farr, T. G., Rosen, P. A., Caro, E., Crippen, R., Duren, R., Hensley, S., Kobrick, M., Paller, M., Rodriguez, E., Roth, L., Seal, D., Shaffer, S., Shimada, J., Umland, J., Werner, M., Oskin, M., Burbank, D., and Alsdorf, D.: The Shuttle Radar Topography Mission, Reviews of Geophysics, 45, RG2004, https://doi.org/10.1029/2005RG000183, http://doi.wiley.com/10.1029/2005RG000183, 2007.

Ferrier, K. L., Huppert, K. L., and Perron, J. T.: Climatic control of bedrock river incision, Nature, 496, 206–209, https://doi.org/10.1038/nature11982, http://dx.doi.org/10.1038/nature11982, 2013.

665

Finnegan, N. J., Hallet, B., Montgomery, D. R., Zeitler, P. K., Stone, J. O., Anders, A. M., and Yuping, L.: Coupling of rock uplift and river incision in the Namche Barwa-Gyala Peri massif, Tibet, Geological Society of America Bulletin, 120, 142–155, https://doi.org/10.1130/B26224.1, https://pubs.geoscienceworld.org/gsabulletin/article/120/1-2/142-155/125457, 2008.

Furbish, D. J. and Roering, J. J.: Sediment disentainment and the concept of local versus nonlocal transport on hillslopes, Journal of Geophysical Research: Earth Surface, 118, 937–952, https://doi.org/10.1002/jgrf.20071, http://doi.wiley.com/10.1002/jgrf.20071, 2013.

670

Gallen, S. F. and Wegmann, K. W.: River profile response to normal fault growth and linkage: An example from the Hellenic forearc of south-central Crete, Greece, Earth Surface Dynamics, 5, 161–186, https://doi.org/10.5194/esurf-5-161-2017, 2017.

Gallen, S. F., Clark, M. K., and Godt, J. W.: Coseismic landslides reveal near-surface rock strength in a high-relief, tectonically active setting, Geology, 43, 11–14, https://doi.org/10.1130/G36080.1, https://pubs.geoscienceworld.org/geology/article/43/1/11-14/131646, 2015.

675

Gasparini, N. M., Whipple, K. X., and Bras, R. L.: Predictions of steady state and transient landscape morphology using sediment-flux-dependent river incision models, Journal of Geophysical Research, 112, F03S09, https://doi.org/10.1029/2006JF000567, http://doi.wiley.com/10.1029/2006JF000567, 2007.

Glade, R. C., Shobe, C. M., Anderson, R. S., and Tucker, G. E.: Canyon shape and erosion dynamics governed by channel-hillslope feedbacks, Geology, 47, 650–654, https://doi.org/10.1130/G46219.1, https://pubs.geoscienceworld.org/gsa/geology/article/47/7/650/570313/canyon-shape-and-erosion-dynamics-governed-by, 2019.

680

Guns, M. and Vanacker, V.: Shifts in landslide frequency-area distribution after forest conversion in the tropical Andes, Anthropocene, 6, 75–85, https://doi.org/10.1016/j.ancene.2014.08.001, http://dx.doi.org/10.1016/j.ancene.2014.08.001, 2014.

Guzzetti, F., Malamud, B. D., Turcotte, D. L., and Reichenbach, P.: Power-law correlations of landslide areas in central Italy, Earth and Planetary Science Letters, 195, 169–183, https://doi.org/10.1016/S0012-821X(01)00589-1, https://linkinghub.elsevier.com/retrieve/pii/S0012821X01005891, 2002.

685
<figure><img src="image_70.png" /></figure>

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

Guzzetti, F., Reichenbach, P., Ardizzone, F., Cardinali, M., and Galli, M.: Estimating the quality of landslide susceptibility models, Geomorphology, 81, 166–184, https://doi.org/10.1016/j geomorph.2006.04.007, https://linkinghub.elsevier.com/retrieve/pii/S0169555X06001371, 2006.

Hancock, G. S. and Anderson, R. S.: Numerical modeling of fluvial strath-terrace formation in response to oscillating climate, Bulletin of the Geological Society of America, https://doi.org/10.1130/0016-7606(2002)114<1131:NMOFST>2.0.CO, 2002.

Hobley, D. E., Sinclair, H. D., Mudd, S. M., and Cowie, P. A.: Field calibration of sediment flux dependent river incision, Journal of Geophysical Research: Earth Surface, 116, F04017, https://doi.org/10.1029/2010JF001935, http://doi.wiley.com/10.1029/2010JF001935, 2011.

Hoffmann, T.: Sediment residence time and connectivity in non-equilibrium and transient geomorphic systems, Earth-Science Reviews, 150, 609–627, https://doi.org/10.1016/j.earscirev.2015.07.008, http://dx.doi.org/10.1016/j.earscirev.2015.07.008, 2015.

Hovius, N., Stark, C. P., and Allen, P. A.: Sediment flux from a mountain belt derived by landslide mapping, Geology, 25, 231–234, https://doi.org/10.1130/0091-7613(1997)025<0231:SFFAMB>2.3.CO;2, 1997.

Hovius, N., Stark, C. P., Hao-Tsu, C., and Jiun-Chuan, L.: Supply and Removal of Sediment in a Landslide-Dominated Mountain Belt: Central Range, Taiwan, The Journal of Geology, 108, 73–89, https://doi.org/10.1086/314387, http://www.journals.uchicago.edu/doi/10.1086/314387, 2000.

Hovius, N., Meunier, P., Lin, C. W., Chen, H., Chen, Y. G., Dadson, S., Horng, M. J., and Lines, M.: Prolonged seismically induced erosion and the mass balance of a large earthquake, Earth and Planetary Science Letters, 304, 347–355, https://doi.org/10.1016/j.epsl.2011.02.005, http://dx.doi.org/10.1016/j.epsl.2011.02.005, 2011.

Howard, A. D. and Kerby, G.: Channel changes in badlands, Geological Society of America Bulletin, 94, 739, https://doi.org/10.1130/0016-7606(1983)94<739:CCIB>2.0.CO;2, https://pubs.geoscienceworld.org/gsabulletin/article/94/6/739-752/202874, 1983.

Iverson, R. M.: Landslide triggering by rain infiltration, Water Resources Research, 36, 1897–1910, https://doi.org/10.1029/2000WR900090, http://doi.wiley.com/10.1029/2000WR900090, 2000.

Iverson, R. M. and George, D. L.: Modelling landslide liquefaction, mobility bifurcation and the dynamics of the 2014 Oso disaster, Géotechnique, 66, 175–187, https://doi.org/10.1680/jgeot.15.LM.004, http://www.icevirtuallibrary.com/doi/10.1680/jgeot.15.LM.004, 2016.

Jeandet, L., Steer, P., Lague, D., and Davy, P.: Coulomb Mechanics and Relief Constraints Explain Landslide Size Distribution, Geophysical Research Letters, 46, 4258–4266, https://doi.org/10.1029/2019GL082351, https://onlinelibrary.wiley.com/doi/abs/10.1029/2019GL082351, 2019.

Kean, J. W. and Smith, J. D.: Flow and boundary shear stress in channels with woody bank vegetation, Riparian Vegetation and Fluvial Geomorphology, 8, 237–252, https://doi.org/10.1029/008WSA17, http://www.agu.org/books/ws/v008/008WSA17/008WSA17.shtml, 2004.

Keefer, D. K.: Landslides caused by earthquakes, Geological Society of America Bulletin, 95, 406, https://doi.org/10.1130/0016-7606(1984)95<406:LCBE>2.0.CO;2, https://pubs.geoscienceworld.org/gsabulletin/article/96/8/1091-1092/203067, https://pubs.geoscienceworld.org/gsabulletin/article/95/4/406-421/10164, 1984.

Keefer, D. K.: Investigating landslides caused by earthquakes - A historical review, Surveys in Geophysics, 23, 473–510, https://doi.org/10.1023/A:1021274710840, 2002.

Keefer, D. K. and Larsen, M. C.: Assessing Landslide Hazards, Science, 316, 1136–1138, https://doi.org/10.1126/science.1143308, https://www.sciencemag.org/lookup/doi/10.1126/science.1143308, 2007.

King, G. E., Herman, F., Lambert, R., Valla, P. G., and Guralnik, B.: Multi-OSL-thermochronometry of feldspar, Quaternary Geochronology, 33, 76–87, https://doi.org/10.1016/j.quageo.2016.01.004, http://dx.doi.org/10.1016/j.quageo.2016.01.004, 2016.
<figure><img src="image_72.png" /></figure>

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

Kirby, E. and Whipple, K. X.: Expression of active tectonics in erosional landscapes, Journal of Structural Geology, 44, 54–75, https://doi.org/10.1016/j.jsg.2012.07.009, http://dx.doi.org/10.1016/j.jsg.2012.07.009, 2012.

725

Korup, O.: Large landslides and their effect on sediment flux in South Westland, New Zealand, Earth Surface Processes and Landforms, 30, 305–323, https://doi.org/10.1002/esp.1143, 2005.

Korup, O.: Rock type leaves topographic signature in landslide-dominated mountain ranges, Geophysical Research Letters, 35, L11402, https://doi.org/10.1029/2008GL034157, http://doi.wiley.com/10.1029/2008GL034157, 2008.

730

Korup, O., Clague, J. J., Hermanns, R. L., Hewitt, K., Strom, A. L., and Weidinger, J. T.: Giant landslides, topography, and erosion, Earth and Planetary Science Letters, 261, 578–589, https://doi.org/10.1016/j.epsl.2007.07.025, https://www.sciencedirect.com/science/article/pii/S0012821X07004785https://linkinghub.elsevier.com/retrieve/pii/S0012821X07004785, 2007.

Korup, O., Densmore, A. L., and Schlunegger, F.: The role of landslides in mountain range evolution, Geomorphology, https://doi.org/10.1016/j geomorph.2009.09.017, 2010.

735

Lague, D.: Reduction of long-term bedrock incision efficiency by short-term alluvial cover intermittency, Journal of Geophysical Research: Earth Surface, 115, n/a–n/a, https://doi.org/10.1029/2008JF001210, http://doi.wiley.com/10.1029/2008JF001210, 2010.

Lague, D.: The stream power river incision model: evidence, theory and beyond, Earth Surface Processes and Landforms, 39, 38–61, https://doi.org/10.1002/esp.3462, http://doi.wiley.com/10.1002/esp.3462, 2014.

740

Lague, D., Hovius, N., and Davy, P.: Discharge, discharge variability, and the bedrock channel profile, Journal of Geophysical Research: Earth Surface, 110, n/a–n/a, https://doi.org/10.1029/2004JF000259, http://doi.wiley.com/10.1029/2004JF000259, 2005.

Larsen, I. J. and Montgomery, D. R.: Landslide erosion coupled to tectonics and river incision, Nature Geoscience, 5, 468–473, https://doi.org/10.1038/ngeo1479, http://dx.doi.org/10.1038/ngeo1479http://www.nature.com/articles/ngeo1479, 2012.

Larsen, I. J., Montgomery, D. R., and Korup, O.: Landslide erosion controlled by hillslope material, Nature Geoscience, 3, 247–251, https://doi.org/10.1038/ngeo776, http://www.nature.com/articles/ngeo776, 2010.

745

Li, G., West, A. J., Densmore, A. L., Hammond, D. E., Jin, Z., Zhang, F., Wang, J., and Hilton, R. G.: Connectivity of earthquake-triggered landslides with the fluvial network: Implications for landslide sediment transport after the 2008 Wenchuan earthquake, Journal of Geophysical Research: Earth Surface, 121, 703–724, https://doi.org/10.1002/2015JF003718, http://doi.wiley.com/10.1002/2015JF003718, 2016.

Lin, G.-W., Chen, H., Chen, Y.-H., and Horng, M.-J.: Influence of typhoons and earthquakes on rainfall-induced landslides and suspended sediments discharge, Engineering Geology, 97, 32–41, https://doi.org/10.1016/j.enggeo.2007.12.001, https://linkinghub.elsevier.com/retrieve/pii/S0013795207002426, 2008.

750

Malamud, B. D. and Turcotte, D. L.: Self-organized criticality applied to natural hazards, in: Natural Hazards, vol. 20, pp. 93–116, https://doi.org/10.1023/A:1008014000515, 1999.

Malamud, B. D., Turcotte, D. L., Guzzetti, F., and Reichenbach, P.: Landslide inventories and their statistical properties, Earth Surface Processes and Landforms, 29, 687–711, https://doi.org/10.1002/esp.1064, http://doi.wiley.com/10.1002/esp.1064, 2004.

755

Marc, O., Hovius, N., Meunier, P., Uchida, T., and Hayashi, S.: Transient changes of landslide rates after earthquakes, Geology, 43, 883–886, https://doi.org/10.1130/G36961.1, https://pubs.geoscienceworld.org/geology/article/43/10/883-886/131712, 2015.

Marc, O., Stumpf, A., Malet, J.-P., Gosset, M., Uchida, T., and Chiang, S.-H.: Initial insights from a global database of rainfall-induced landslide inventories: the weak influence of slope and strong influence of total storm rainfall, Earth Surface Dynamics, 6, 903–922, https://doi.org/10.5194/esurf-6-903-2018, https://www.earth-surf-dynam.net/6/903/2018/, 2018.
<figure><img src="image_74.png" /></figure>

760 Marc, O., Behling, R., Andermann, C., Turowski, J. M., Illien, L., Roessner, S., and Hovius, N.: Long-term erosion of the Nepal Himalayas by bedrock landsliding: The role of monsoons, earthquakes and giant landslides, Earth Surface Dynamics, 7, 107–128, https://doi.org/10.5194/esurf-7-107-2019, 2019.

Meunier, P., Hovius, N., and Haines, A. J.: Regional patterns of earthquake-triggered landslides and their relation to ground motion, Geophysical Research Letters, 34, L20408, https://doi.org/10.1029/2007GL031337, http://doi.wiley.com/10.1029/2007GL031337, 2007.

765 Milliman, J. D. and Meade, R. H.: World-Wide Delivery of River Sediment to the Oceans, The Journal of Geology, https://doi.org/10.1086/628741, 1983.

Montgomery, D. R. and Dietrich, W. E.: A physically based model for the topographic control on shallow landsliding, Water Resources Research, 30, 1153–1171, https://doi.org/10.1029/93WR02979, http://doi.wiley.com/10.1029/93WR02979, 1994.

770 Montgomery, D. R. and Gran, K. B.: Downstream variations in the width of bedrock channels, Water Resources Research, 37, 1841–1846, https://doi.org/10.1029/2000WR900393, http://doi.wiley.com/10.1029/2000WR900393, 2001.

NASA JPL: NASA Shuttle Radar Topography Mission Global 1 arc second DEM, Version 3, Tech. rep., NASA LP DAAC, https://doi.org/http://doi.org/10.5067/MEaSURES/SRTM/SRTMGL1.003, 2013.

775 Niemi, N. A., Oskin, M., Burbank, D. W., Heimsath, A. M., and Gabet, E. J.: Effects of bedrock landslides on cosmogenically determined erosion rates, Earth and Planetary Science Letters, 237, 480–498, https://doi.org/10.1016/j.epsl.2005.07.009, https://www.sciencedirect.com/science/article/pii/S0012821X0500467Xwww.elsevier.com/locate/epslhttp://linkinghub.elsevier.com/retrieve/pii/S0012821X0500467X, 2005.

Ouimet, W. B., Whipple, K. X., Royden, L. H., Sun, Z., and Chen, Z.: The influence of large landslides on river incision in a transient landscape: Eastern margin of the Tibetan Plateau (Sichuan, China), Geological Society of America Bulletin, 119, 1462–1476, https://doi.org/10.1130/B26136.1, https://pubs.geoscienceworld.org/gsabulletin/article/119/11-12/1462-1476/125374, 2007.

780 Page, M. J., Reid, L. M., and Lynn, I. H.: New Zealand Hydrological Society Sediment production from Cyclone Bola landslides, Waipaoa catchment, Source Journal of Hydrology (New Zealand), 38, 289–308, http://www.jstor.org/stable/43944823, 1999.

Paola, C. and Voller, V. R.: A generalized Exner equation for sediment mass balance, Journal of Geophysical Research: Earth Surface, 110, n/a–n/a, https://doi.org/10.1029/2004JF000274, http://doi.wiley.com/10.1029/2004JF000274, 2005.

785 Parker, R. N., Hales, T. C., Mudd, S. M., Grieve, S. W. D., and Constantine, J. A.: Colluvium supply in humid regions limits the frequency of storm-triggered landslides, Scientific Reports, 6, 34438, https://doi.org/10.1038/srep34438, http://dx.doi.org/10.1038/srep34438http://www.nature.com/articles/srep34438, 2016.

Pelletier, J. D.: Minimizing the grid-resolution dependence of flow-routing algorithms for geomorphic applications, Geomorphology, 122, 91–98, https://doi.org/10.1016/j geomorph.2010.06.001, https://linkinghub.elsevier.com/retrieve/pii/S0169555X10002606http://dx.doi.org/10.1016/j.geomorph.2010.06.001, 2010.

790 Roering, J. J., Kirchner, J. W., and Dietrich, W. E.: Evidence for nonlinear, diffusive sediment transport on hillslopes and implications for landscape morphology, Water Resources Research, 35, 853–870, https://doi.org/10.1029/1998WR900090, http://doi.wiley.com/10.1029/1998WR900090, 1999.

Scherler, D., DiBiase, R. A., Fisher, G. B., and Avouac, J.-P.: Testing monsoonal controls on bedrock river incision in the Himalaya and Eastern Tibet with a stochastic-threshold stream power model, Journal of Geophysical Research: Earth Surface, 122, 1389–1429, https://doi.org/10.1002/2016JF004011, http://doi.wiley.com/10.1002/2016JF004011, 2017.

795 Schwanghart, W. and Kuhn, N. J.: TopoToolbox: A set of Matlab functions for topographic analysis, Environmental Modelling and Software, 25, 770–781, https://doi.org/10.1016/j.envsoft.2009.12.002, http://dx.doi.org/10.1016/j.envsoft.2009.12.002, 2010.
<figure><img src="image_75.png" /></figure>

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

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

800 Seidl, M. A. and Dietrich, W. E.: The problem of channel erosion into bedrock, Catena Supplement, 23, 101–124, 1992.

Shobe, C. M., Tucker, G. E., and Anderson, R. S.: Hillslope-derived blocks retard river incision, Geophysical Research Letters, 43, 5070–5078, https://doi.org/10.1002/2016GL069262, http://doi.wiley.com/10.1002/2016GL069262, 2016.

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

Shobe, C. M., Tucker, G. E., and Rossi, M. W.: Variable-Threshold Behavior in Rivers Arising From Hillslope-Derived Blocks, Journal of Geophysical Research: Earth Surface, 123, 1931–1957, https://doi.org/10.1029/2017JF004575, https://onlinelibrary.wiley.com/doi/abs/10.1029/2017JF004575, 2018.

Sidle, R. C. and Ochiai, H.: Landslides: Processes, Prediction, and Land Use, Water Resources Monograph, American Geophysical Union, Washington, D. C., https://doi.org/10.1029/WM018, https://pubs.geoscienceworld.org/eeg/article/13/1/90-91/136800, http://doi.wiley.com/10.1029/WM018, 2006.

Sklar, L. S. and Dietrich, W. E.: A mechanistic model for river incision into bedrock by saltating bed load, Water Resources Research, 40, W06 301, https://doi.org/10.1029/2003WR002496, http://doi.wiley.com/10.1029/2003WR002496, 2004.

Snyder, N. P., Whipple, K. X., Tucker, G. E., and Merritts, D. J.: Importance of a stochastic distribution of floods and erosion thresholds in the bedrock river incision problem, Journal of Geophysical Research: Solid Earth, 108, 2388, https://doi.org/10.1029/2001JB001655, http://doi.wiley.com/10.1029/2003JB002649, http://doi.wiley.com/10.1029/2001JB001655, 2003.

Stark, C. P. and Hovius, N.: The characterization of landslide size distributions, Geophysical Research Letters, 28, 1091–1094, https://doi.org/10.1029/2000GL008527, http://doi.wiley.com/10.1029/2000GL008527, 2001.

Taylor, D.: Fundamentals of Soil Mechanics, Wiley, New York, 1948.

Tenorio, G. E., Vanacker, V., Campforts, B., Álvarez, L., Zhiminaicela, S., Vercruysse, K., Molina, A., and Govers, G.: Tracking spatial variation in river load from Andean highlands to inter-Andean valleys, Geomorphology, 308, 175–189, https://doi.org/10.1016/j.geomorph.2018.02.009, http://linkinghub.elsevier.com/retrieve/pii/S0169555X18300576, 2018.

Tucker, G. E.: Drainage basin sensitivity to tectonic and climatic forcing: Implications of a stochastic model for the role of entrainment and erosion thresholds, Earth Surface Processes and Landforms, 29, 185–205, https://doi.org/10.1002/esp.1020, 2004.

Tucker, G. E. and Hancock, G. R.: Modelling landscape evolution, Earth Surface Processes and Landforms, 35, 28–50, https://doi.org/10.1002/esp.1952, http://www3.interscience.wiley.com/journal/121517813/abstract, http://doi.wiley.com/10.1002/esp.1952, 2010.

Turowski, J. M., Lague, D., and Hovius, N.: Cover effect in bedrock abrasion: A new derivation and its implications for the modeling of bedrock channel morphology, Journal of Geophysical Research, 112, F04 006, https://doi.org/10.1029/2006JF000697, http://doi.wiley.com/10.1029/2006JF000697, 2007.

Turowski, J. M., Lague, D., and Hovius, N.: Response of bedrock channel width to tectonic forcing: Insights from a numerical model, theoretical considerations, and comparison with field data, Journal of Geophysical Research, 114, F03 016, https://doi.org/10.1029/2008JF001133, http://doi.wiley.com/10.1029/2008JF001133, 2009.

Van Asch, T., Buma, J., and Van Beek, L.: A view on some hydrological triggering systems in landslides, Geomorphology, 30, 25–32, https://doi.org/10.1016/S0169-555X(99)00042-2, https://linkinghub.elsevier.com/retrieve/pii/S0169555X99000422, 1999.
<figure><img src="image_77.png" /></figure>

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

Van Rompaey, A. J. J. and Govers, G.: Data quality and model complexity for regional scale soil erosion prediction, International Journal of Geographical Information Science, 16, 663–680, https://doi.org/10.1080/13658810210148561, http://www.informaworld.com/openurl? genre=article{&}doi=10.1080/13658810210148561{&}magic=crossref, 2002.

Wang, J., Jin, Z., Hilton, R. G., Zhang, F., Densmore, A. L., Li, G., and West, A. J.: Controls on fluvial evacuation of sediment from earthquake-triggered landslides, Geology, 43, 115–118, https://doi.org/10.1130/G36157.1, http://pubs.geoscienceworld.org/geology/ article/43/2/115/131758/Controls-on-fluvial-evacuation-of-sediment-from, 2015.

Wang, W., Godard, V., Liu-Zeng, J., Scherler, D., Xu, C., Zhang, J., Xie, K., Bellier, O., Ansberque, C., and de Sigoyer, J.: Perturbation of fluvial sediment fluxes following the 2008 Wenchuan earthquake, Earth Surface Processes and Landforms, 42, 2611–2622, https://doi.org/10.1002/esp.4210, http://doi.wiley.com/10.1002/esp.4210, 2017.

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, Journal of Geophysical Research: Solid Earth, 104, 17661–17674, https://doi.org/10.1029/1999JB900120, http://doi.wiley.com/10.1029/1999JB900120, 1999.

Whipple, K. X., Hancock, G. S., and Anderson, R. S.: River incision into bedrock: Mechanics and relative efficacy of plucking, abrasion, and cavitation, Geological Society of America Bulletin, 112, 490–503, https://doi.org/10.1130/0016-7606(2000)112<490:RIIBMA>2.0.CO;2, https://pubs.geoscienceworld.org/gsabulletin/article/112/3/490-503/183620, 2000.

Willgoose, G., Bras, R. L., and Rodrigueziturbe, I.: Results from a New Model of River Basin Evolution, Earth Surface Processes and Landforms, 16, 237–254, https://doi.org/10.1002/esp.3290160305, 1991.

Wobus, C., Whipple, K. X., Kirby, E., Snyder, N., Johnson, J., Spyropolou, K., Crosby, B., and Sheehan, D.: Tectonics from topography: Procedures, promise, and pitfalls, in: Tectonics, Climate, and Landscape Evolution, vol. 398, pp. 55–74, Geological Society of America, https://doi.org/10.1130/2006.2398(04), https://pubs.geoscienceworld.org/books/book/569/chapter/3802975/https://pubs.geoscienceworld.org/books/book/569/chapter/3802975, 2006.

Wyllie, D. C. and Mah, C. W.: Rock Slope Engineering, CRC Press, https://doi.org/10.1201/9781315274980, https://www.taylorfrancis.com/books/9781482265125, 2017.

Yanites, B. J.: The Dynamics of Channel Slope, Width, and Sediment in Actively Eroding Bedrock River Systems, Journal of Geophysical Research: Earth Surface, 123, 1504–1527, https://doi.org/10.1029/2017JF004405, https://onlinelibrary.wiley.com/doi/abs/10.1029/2017JF004405, 2018.

Yanites, B. J., Tucker, G. E., and Anderson, R. S.: Numerical and analytical models of cosmogenic radionuclide dynamics in landslide-dominated drainage basins, Journal of Geophysical Research, 114, F01007, https://doi.org/10.1029/2008JF001088, http://doi.wiley.com/10.1029/2008JF001088, 2009.

Yanites, B. J., Tucker, G. E., Mueller, K. J., and Chen, Y.-G.: How rivers react to large earthquakes: Evidence from central Taiwan, Geology, 38, 639–642, https://doi.org/10.1130/G30883.1, http://pubs.geoscienceworld.org/geology/article/38/7/639/130315/How-rivers-react-to-large-earthquakes-Evidence, 2010.

Zhang, J., J Van Westen, C., Tanyas, H., Mavrouli, O., Ge, Y., Bajrachary, S., Raj Gurung, D., Raj Dhital, M., and Raj Khanal, N.: How size and trigger matter: Analyzing rainfall- And earthquake-triggered landslide inventories and their causal relation in the Koshi River basin, central Himalaya, Natural Hazards and Earth System Sciences, 19, 1789–1805, https://doi.org/10.5194/nhess-19-1789-2019, 2019.

Zhou, S., Ouyang, C., An, H., Jiang, T., and Xu, Q.: Comprehensive study of the Beijing Daanshan rockslide based on real-time videos, field investigations, and numerical modeling, Landslides, https://doi.org/10.1007/s10346-020-01345-2, http://link.springer.com/10.1007/s10346-020-01345-2, 2020.