Articles | Volume 13, issue 9
Model description paper
 | Highlight paper
31 Aug 2020
Model description paper | Highlight paper |  | 31 Aug 2020

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

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 interactions among these 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 TopoToolbox 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 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 are calculated using a Mohr–Coulomb stability analysis, while landslide-derived sediment is routed and deposited using a multiple-flow-direction, nonlinear 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 conditions. 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 provide 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 Tucker1999; 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 Meade1983). 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 Montgomery2012; Ouimet et al.2007; Korup et al.2010; Shobe et al.2016; Glade et al.2019).

Unraveling 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 interdependencies over various temporal and spatial scales (Tucker and Hancock2010). The past 20 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 Anderson2002; Turowski et al.2007). Most existing LEMs simulate river incision using one of two commonly used end-member models (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 riverbed and quantified using a scaling law between fluid stress and river incision rate (Seidl and Dietrich1992; Howard and Kerby1983; Campforts and Govers2015). 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 Voller2005). 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 Tucker1999).

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 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 Lague2009) and can simultaneously simulate fluvial sediment transport and river incision into bedrock (e.g., Lague2010). 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.

To understand how landslides influence landscape evolution, Densmore et al. (1998) proposed an approach, adapted by others (e.g., Champel2002; 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 (Taylor1948); (ii) landslides initialize near the toes of hillslopes; and (iii) landslide-derived sediment is spread under a constant slope, following the steepest downslope path. 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.

Other researchers have developed mechanistic models to simulate shallow landslide activity at the landscape scale (e.g., Montgomery and Dietrich1994; 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; Iverson2000; 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; see Iverson and George2016), they typically involve a range of geophysical processes and require input parameters which may not be adequately known at large spatial scales. This can make the more detailed models sensitive to equifinality (Beven and Freer2001). Moreover, deep-seated bedrock landslides, rather than shallow landslides, mobilize the largest volumes of sediment and therefore have the largest impact on landscape evolution (Burbank2002; Dussauge et al.2003; Jeandet et al.2019; Korup et al.2007; Broeckx et al.2020).

Notwithstanding the prominent role of landslides in shaping Earth's surface and controlling sediment supply and transport, few efforts have been made to actively simulate the impact of stochastic landsliding on landscape evolution and sediment dynamics over large spatial and temporal scales. 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 modeled using the SPACE fluvial incision framework (Shobe et al.2017), and hillslope-derived sediment discharge is 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 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 TopoToolbox 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:

(1) η t = R t + H t = U - E r fluv + D s fluv - E s fluv 1 - ϕ sed - E r hill + D s hill - E s hill 1 - ϕ sed ,

where η (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−1) is the rock uplift rate, and ϕsed is the sediment porosity. Erfluv (L T−1) 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. Esfluv (L T−1) is the fluvial volumetric entrainment flux of sediment per unit bed area, and Dsfluv (L T−1) is the fluvial volumetric deposition flux of sediment per unit bed area. Erhill (L T−1) is the volumetric flux of hillslope bedrock erosion due to landsliding per unit bed area, representing the amount of bedrock that is detached. Eshill (L T−1) is the volumetric entrainment flux of sediment erosion (produced by landsliding or creep) per unit bed area and Dshill (L T−1) is the volumetric deposition flux of hillslope-derived sediment per unit bed area. Note that HyLands does not explicitly distinguish between river or hillslope cells: all equations are applied to the entire landscape. Processes affecting sediment thickness and bedrock elevation in each cell can be either fluvial dynamics (SPACE) or landslides – or a combination of both – hence the hybrid nature of HyLands.

2.1 River sediment transport and bedrock erosion

HyLands uses the 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 Lague2009), 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 riverbed compositions from bare bedrock channels to mixed bedrock–alluvial channels to fully alluvial channels. This is accomplished by combining mass conservation of riverbed sediment with a bedrock erosion law to simultaneously solve for the time evolution of the bedrock and sediment surfaces. Note that this approach implies that all river cells in the landscape are assumed to occupy one grid cell of width 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. 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 discharge Qsfluv (L3 T−1) per unit width w (L) is written as

(2) Q s fluv / w x = E s fluv + 1 - F f fluv E r fluv - D s fluv ,

where Fffluv is a unitless fraction of fine fluvial sediment. The factor 1-Fffluv (-) 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 the 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 (L2 T−1) and local channel slope S (e.g., Howard and Kerby1983; Whipple and Tucker1999; Davy and Lague2009). The sediment entrainment rate Esfluv and the bedrock erosion rate Erfluv 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 increasing H/H* and a concomitant exponential decrease in bedrock erosion rate with increasing 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)

(3) E s fluv = K s q S n 1 - e - H / H *

for sediment and

(4) E r fluv = K r q S n e - H / H *

for bedrock. Ks and Kr (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 does not need to be 1 for the SPACE model in general. There are a variety of ways to compute q. We use the common approach of calculating discharge as a function of drainage area such that q=kqAm, where m is a scaling exponent and kq is a coefficient subsumed into the fluvial erosion coefficients Ks and Kr.

Sediment deposition is implemented similarly to Davy and Lague (2009) such that the deposition flux depends on volumetric sediment discharge Qsfluv (L3 T−1) divided by the volumetric water discharge Q (L3 T−1) and the effective sediment settling velocity V (L T−1) as follows:

(5) D s fluv = Q s fluv Q V .

V is the net effective settling velocity, which represents the still-water particle settling velocity corrected for the upward effects of turbulence and the vertical gradient in sediment concentration through the water column (Davy and Lague2009). HyLands enables spatially variable values for V to distinguish between settling velocities in flooded and non-flooded areas.

Figure 1Sketch of fluvial SPACE component. Model setup and variable definitions for the SPACE bedrock–alluvial river erosion model. 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 rates of sediment entrainment Esfluv and bedrock erosion Erfluv are set by the ratio of sediment thickness H to the bedrock roughness height H* (Fig. 2). Modified from Shobe et al. (2017).


Figure 2Relative efficiency of fluvial sediment entrainment and bedrock erosion, fH/H*, as a function of the ratio of sediment thickness to bedrock roughness H/H*. fH/H* depends only on the ratio H/H*, varies between 0 and 1, and indicates the proportion of total stream power used to entrain sediment – see Eq. 3; fH/H*=1-e-H/H* (dashed line) – or erode bedrock – see Eq. 4; fH/H*=e-H/H* (solid 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. Modified from Shobe et al. (2017).


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 Eshill and Erhill, 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 θc bisects the local topographic slope β and the material's angle of internal friction ϕ (Densmore et al.1998; Champel2002) as follows:

(6) θ c = β + ϕ 2 .

Note that Eq. (6) implies that failure occurs at angles higher than ϕ, being the consequence of the force balance controlled by material properties (Eq. 8) and the mass of the wedge above the failure plane (Culmann1875). 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 θc is extended until it daylights (i.e., intersects the topographic surface).

Figure 3Sketch of landslide algorithm in two dimensions. Landslide erosion (red shaded area) is calculated using the Culmann approach (Culmann1875). β is the topographic slope; ϕ represents the angle of internal friction, and θ is the inclination of the rupture plane. Deposition of landslide-derived sediment (green shaded area) is calculated using a nonlocal diffusion equation (Eq. 11; see Carretier et al.2016). δ is the minimal deposit surface angle at 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 (Ffhill=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 downslope gradient which is lower than δ, the slope of the deposited volume is adjusted so that the deposit surface angle equals δ. Probability for sliding is calculated as the ratio of the local hillslope height Hs to the maximum hillslope height Hc (Eq. 7; see Densmore et al.1998). The inset plot illustrates that Hc depends on the rock strength (cohesion C and internal friction angle ϕ) and the topographic slope β. The plotted lines are calculated using Eq. 8, where ρ is 2700 kgm−3 and g = 9.81 m s−2.


Modeling landslide frequency and location depends critically on the identification of points in the landscape where landslides initiate. A wide variety of events ranging from coseismic activity and peak ground acceleration (Meunier et al.2007) over intense storm events (Marc et al.2018) to human hillslope destabilization (Guns and Vanacker2014) 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 cells as points in the landscape where the topographic slope (β) exceeds the angle of internal friction (ϕ). For all unstable cells, the probability for sliding, pLS, is calculated as

(7) p LS = H s H c ,

where Hs is the local hillslope height calculated as difference between every cell in the landscape and its highest neighbor (Fig. 3), and Hc is the maximum stable hillslope height, which is calculated as (Densmore et al.1998; Champel2002)

(8) H c = 4 C ρ g sin β cos ϕ 1 - cos ( β - ϕ ) .

Here C is the cohesion (M L−1 T−2), ρ is the rock density (M L−3) set to 2700 kg m−3, and g is the gravitational acceleration (g=9.81 ms−2). To simulate the random nature of landslides, grid cells where landslides initiate are selected using a stochastic sampling scheme:

(9) r d t t LS > p LS No landsliding < p LS Landsliding, select as critical cell, 

where r is a random number between 0 and 1, and tLS is the return time for landslides with tLS>=dt, where dt is the model time step. Unstable cells where a landslide is induced will further be referred to as critical cells (Fig. 3). At every model iteration, Eq. (9) is updated for all cells 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 Ochiai2006). Jeandet et al. (2019) inverted several landslide inventories and found effective cohesion values ranging between 10 and 35 kPa. These values are lower than geomechanical values representing large-scale rock strength (e.g., Densmore et al.1998; Champel2002), 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 tLS controls the absolute number of critical cells where landslides are initiated. If tLS equals the time step dt, the number of landslides per time step is solely controlled by the ratio HsHc. When tLSHs/Hc, however, the number of landslides per time step is reduced (Eq. 9). While the HsHc ratio controls the topographic location of landsliding onset and thus the landslide characteristics (size and volume), the landslide return time tLS controls the absolute number of landslides and therefore overall landslide erosion rates.

HyLands enables the simulation of landslides at every location in the landscape. At every iteration, landslides are induced at critical cells 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 cell, we build a stack of unstable (i.e., sliding) cells. The stack is initialized by adding the critical cell. Next, a recursive procedure is applied until the stack is empty. This procedure consists of the following steps: (i) select the first cell from the stack (thus starting with the critical cell). (ii) Evaluate all up-slope neighboring cells. If the elevation of a neighboring cell exceeds the elevation of the sliding plane defined by θc, the cell is identified as a sliding cell and added the to stack of landslide cells. (iii) Remove the first cell from the stack. This procedure is repeated until the stack is empty. HyLands offers the possibility to set a maximum landslide area (ALSmax). Once this maximum is achieved – or when no more cells are added to the stack of landslide cells – the landslide area is defined. All cells inside this area are eroded to the elevation of the sliding plane, thereby adjusting both Eshill and Erhill of Eq. (1) for all cells involved.

2.2.2 Hillslope-derived sediment

The spatial change in hillslope-derived volumetric sediment discharge Qshill (L3 T−1) per unit width w (L) is written as

(10) Q s Hill / w x = E s hill + 1 - F f hill E r hill - D s hill ,

where Ffhill is a unitless fraction of fine hillslope-derived sediment. The factor 1-Ffhill (–) represents the idea that some fraction of the hillslope-derived sediment is rapidly evacuated in high suspension (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 Ffhill=0, the system is fully mass conservative, and all sediments produced by landslide activity contribute to the sediment discharge (Eq. 1). Note that both Eshill and Dshill are corrected for dimensionless sediment porosity (ϕsed in Eq. 1), thereby factoring in the increase in bedrock-derived landslide sediment volumes during deposition due to the addition of pore spaces.

2.2.3 Deposition of landslide material

In the following we describe how deposition of landslide-derived material is being calculated. 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 (Pelletier2010). Moreover, the use 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 much longer distances than sediment deposited on steep parts. 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 Hanks1985). However, such an approach is not suited to simulate the distribution of landslide-derived sediment. While diffusion equations distribute sediment only between the neighboring cells of a cell, landslide-derived sediment has run-out distances that can be significantly longer than a single grid cell (Claessens et al.2007). Therefore, we adopt a nonlinear, nonlocal deposition scheme for landslide-derived sediment outlined by Carretier et al. (2016):

(11) D s hill = Q s Hill / w L ,

where Dshill (L T−1) is the volumetric deposition flux of hillslope-derived sediment per unit bed area and L (L) represents a sediment transport distance. The larger the L, the bigger the distance over which sediments are transported and the lower the local deposition rate. L is calculated for every grid cell as follows:

(12) L = d x 1 - S S c 2 ,

where Sc is a critical slope which we further assume to be equal to the angle of internal friction (ϕ). When the hillslope gradient SSc, 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 Roering2013; Carretier et al.2016). When S approaches Sc, L goes to infinity implying that no deposition will occur at the considered cell. At steep slopes, sediment transport therefore shows nonlocal behavior in the sense that erosion of nonlocal, upstream cells is integrated when calculating the sediment discharge Qs (Carretier et al.2016). When S>Sc, L is set to , and there is no deposition at that cell. For negative values of S, which might occur for flooded cells, L is set to dx. A minimal angle (δ) under which landslide-derived sediments should be deposited can be set as a parameter value but is not required to run HyLands (Fig. 3).

Contrary to the fluvial component of the model (SPACE) where a single flow direction algorithm (steepest descent) is used, landslide-derived sediment is spread over the landscape using a multiple-flow-direction algorithm redistributing sediment over the downstream cells in proportion to the local slope (Fig. 4; see 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 cells during every model iteration. Lakes are identified by filling all sinks in a landscape to the brim (Schwanghart and Scherler2014). By default, flooded cells do not erode but do allow for sediment deposition (Eq. 5). HyLands enables spatially variable values for V (Eq. 5) to distinguish between settling velocities in non-flooded versus flooded cells by changing the values for V and VLake, respectively (see Table 1).

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).

Figure 4Sketch of landslide algorithm in three dimensions. Cells shaded in blue indicate the critical cells where landslides potentially 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 (see Eqs. 11 and 12).


3 Verification and evaluation

3.1 Comparison to analytical solutions for the fluvial component

In the first three test runs (detachment-limited, transport-limited and mixed), a steady-state artificial landscape is simulated using a square grid of 20 by 20 cells with a spatial resolution of 100 m. The landsliding component is turned off for these runs (i.e., only fluvial erosion according to the SPACE model occurs). Each run is initialized from a surface with randomly generated microtopography. The initial surface is a tilted plane which drains towards the southwestern corner, the only open boundary cell. 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 time step is set to 10 years. Under detachment-limited conditions (imposed by setting Ff=1), the sediment thickness H equals 0 everywhere 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 (qAm) it has been shown that fluvial erosion results in the following steady-state slope–area relationship (Shobe et al.2017):

(13) S = U K r A m 1 / n .

Figure 5b 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 conditions where HH*, the theoretical slope–area relationship for fluvial incision can be written as (Shobe et al.2017)

(14) S = V r + 1 1 / n U K s 1 / n A - m / n ,

where r represents the runoff rate (L T−1). 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 5d illustrates the slope–area plot for all cells of the simulated steady-state landscape showing a close match with the analytical prediction (14). Moreover, HyLands also reproduces the theoretical steady-state sediment discharge relationship for transport-limited conditions with Fffluv=0 (Shobe et al.2017):

(15) Q s fluv = U A .

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

(16) S = U V K s A m r + U K r A m 1 / n .

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

(17) η t = R t + H t = 0 + 0 = 0 ,

which leads to a constant soil thickness over the landscape derived by (Shobe et al.2017)

(18) H = - H * ln 1 - V K s r K r + V .

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. 5g) matches the theoretical relationship in Eq. (16) and so does the soil thickness H, which evolves toward a constant thickness (Eq. 18), and Qs, adhering to the sediment discharge relationship (Eq. 15).

Table 1Parameter values.

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 (Fffluv=1). Landslide parameters are only relevant for models where landslide activity is simulated. a The synthetic landscape evolution model consists of three 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 Time spent to complete one full model iteration on a Windows PC, with an Intel® Core i9-9880H CPU at 2.30 GHz and a RAM of 32 GB. c HyLands enables spatially variable values for V (Eq. 5) to distinguish between settling velocities in non-flooded versus flooded cells by changing the values for V and VLake, respectively.

Download Print Version | Download XLSX

Figure 5Verification of fluvial SPACE component. (a) Longitudinal profile of the trunk stream at steady state under detachment-limited conditions where no sediment is present because Fffluv=1, and all produced sediment is assumed to be evacuated instantaneously. At steady state, the concave-upward profile is in equilibrium with the imposed rock uplift pattern. The steady-state slope–area relationship (b) matches the analytical solution (Eq. 13). (c) Longitudinal profile of the trunk stream at steady state under transport-limited conditions. At steady state, the concave-upward profile is in equilibrium with the imposed rock uplift pattern. The steady-state slope–area relationship (d) matches the analytical solution (Eq. 14). Panel (e) illustrates the steady-state fluvial sediment discharge (Qsfluv) as a function of the drainage area, which matches the analytical area–sediment discharge relationship (Eq. 15). (f) Longitudinal profile of the trunk stream at steady state under mixed bedrock–alluvial conditions. At steady state, both the topographic and bedrock profiles are in equilibrium with the imposed rock uplift pattern. The steady-state slope–area relationship (g) matches the analytical solution (Eq. 16). Panel (h) illustrates the steady-state fluvial sediment discharge (Qsfluv) as a function of the drainage area, which matches the analytical area–sediment discharge relationship (Eq. 15).


3.2 Evaluation of the landslide component

Because landsliding is a stochastic process, it is 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 Turcotte1999; Stark and Hovius2001; Guzzetti et al.2002; Korup2005; Guns and Vanacker2014; Larsen and Montgomery2012). 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 Hovius2001; Malamud and Turcotte1999; Guzzetti et al.2002). Magnitude–frequency distributions are typically described using a three-parameter inverse gamma distribution as (Malamud et al.2004)

(19) p ( A L ; ρ l , a l , s l ) = 1 a l Γ ( ρ l ) a l A L - s l ρ l + 1 exp - a l A L - s l ,

where AL is the landslide area (L2), p(AL) is the probability density of a landslide area (AL), al, sl, and ρl are empirical parameters, and Γ(ρl) is the gamma function of ρl. A second empirical universal relationship is the volume–area scaling relationship where the volume VL of a given landslide is a function of its area AL as (Hovius et al.1997)

(20) V L = α l A L γ l ,

where αl is an intercept and γ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 predating 1974 and an inventory containing all landslide events between 1974 and 2007 (Fig. 8a). 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. We selected the Namche Barwa–Gyala Peri massif to evaluate the performance of HyLands given its unique geomorphologic configuration featuring amongst the highest globally documented river stream power in combination with very active hillslope processes (Larsen and Montgomery2012). With HyLands being designed to couple the role of fluvial and hillslope processes, this region makes for a good test environment. Note however that we do not intend to calibrate nor validate the model but run it using fixed, theoretical model parameters (Sect. 3.2.3). Applications of HyLands aiming to constrain the model through parameter calibration and validation (Sect. 4.3) would require additional data to ground-truth landslide inventories and to provide detailed records on landslide triggers such as earthquakes and storms.

Figure 6Namche Barwa–Gyala Peri massif used for model evaluation. The red rectangle on the inset figure indicates the location of the study area. The green dashed rectangle indicates the part of the DEM used to evaluate HyLands. The shaded colors show elevation derived from the 30 m SRTM v3 DEM (Farr et al.2007) and resampled to a higher resolution of 20 m using a bicubic interpolation method. Main map is produced with TopoToolbox (Schwanghart and Scherler2014). Inset map is made in QGIS 3.10.9, using Natural Earth vector and raster map data available at (last access: 25 August 2020).

3.2.3 Model parameterization

We run HyLands using the NASA Shuttle Radar Topography Mission (SRTM) v3.0 elevation dataset as an initial surface (Farr et al.2007), resampled to a higher resolution of 20 m using a bicubic interpolation method. We resampled the DEM to a resolution of 20 m in order to evaluate the capacity of HyLands to reproduce the rollover in the magnitude frequency distribution, often reported to occur for landslide areas <900 m2, which would be the minimum landslide area when using the original SRTM data. 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 and 4000 W m2; Finnegan et al.2008. The simulated grid is composed of 1918×1149 cells, covering a total area of ca. 960 km2. 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 Whipple2012).

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 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 (ϕsed=0 and Fffluv=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 (Ffhill=0.25). We assume that the angle of internal friction (ϕ) is comparable to the mode of the topographical slope distribution (Burbank et al.1996; Korup2008; Montgomery and Gran2001), reported to range between 37 and 39 and here set to 38 (Larsen and Montgomery2012). Cohesion C is known to vary over a wide range and strongly depends on rock mechanical properties (Wyllie and Mah2017). 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 (12–20 Pa in Jeandet et al.2019). Cohesion and the angle of internal friction influence the size distribution of landslides in several ways (Jeandet et al.2019). The angle of internal friction ϕ controls the angle of the potential rupture plane such that lower values of ϕ 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 Hc in Eq. (6). Larger values for C will result in larger values for Hc, thus decreasing the probability of landslides on gentler slopes and resulting in fewer small landslides. The minimum value of the minimal deposit surface angle under which landslide-derived sediment is redistributed on hillslopes is set to 0.01. We set the return time for landsliding tLS to 2×104 years. tLS regulates the probability that unstable cells evolve into a landslide (Eq. 9) and therefore controls the number of landslide events per time step Δt. When applying HyLands to reconstruct or predict landslide activity, tLS should be a function of the frequency (or return time) of triggering events (large earthquakes or rainfall events). Evaluation of model sensitivity to changing values for tLS is a natural avenue for further work. A full overview of the model parameters is given in Table 1.


Table 2Simulation movies and script names.

Download Print Version | Download XLSX

Figure 7Time slices of HyLands model run for the Namche Barwa–Gyala Peri massif after 5, 165, 330 and 500 model years. Panels (a), (d), (g) and (j) indicate the location of the landslides at the given time step (black diamonds). Panels (b), (e), (h) and (k) zoom into the red squares in (a), (d), (g) and (j) showing simulated landslide erosion–deposition. The colors represent the square root (for ease of visualization) of landslide erosion () and deposition (+) during the presented time step. Panels (c), (f), (i) and (l) show the sediment thickness, H, generated through river incision (SPACE module) and landsliding. The underlying hillshade was derived from the 30 m SRTM v3 DEM (Farr et al.2007).

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 cells) are well spread over the landscape. The number of landslide events (the number of black diamonds in Fig. 7) is mainly controlled by the return time for landsliding tLS. 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 cells become susceptible to erosion and deposition as landscape evolution continues (e.g., the deposition pattern in Fig. 7h reflects the shape of erosion patterns resulting from previous landslide activity).

To quantitatively evaluate the performance of the HyLands landslide module, we compare modeled landslide properties against observed scaling relationships (Eqs. 19 and 20). Figure 8a compares the modeled 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 (see Sect. 2.2.1) and the topography of the studied region (Jeandet et al.2019). Figure 8b shows that HyLands is capable of approaching the shape of the universal area–volume relationships found by Larsen et al. (2010). While HyLands seems to overestimate simulated landslide volumes for very small landslides, the area–volume relationship simulated with HyLands approaches a linear relationship in a log–log space for larger landslides, similar to the shape of the observed area–volume relationship. Note that overall, landslide volumes simulated with HyLands are overpredicted compared to observations. Study-area-specific model calibration would improve this fit but is beyond the scope of this model evaluation in which we evaluate the capacity of HyLands to reproduce the shape of the universal area–volume relationship. We attribute the positive deviation from the linear area–volume relationship in a log–log space for smaller landslides to the nature of the landslide algorithm. HyLands simulates deep-seated landslides; several of the smaller landslides are likely to be shallow landslides, which are currently not simulated. Moreover, HyLands does not allow any sediment to be deposited within the landslide scar, while this typically does occur in nature. Future developments of the algorithm could allow for shallow landsliding and in-scar deposition for more realistic simulations. Furthermore, 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.

Figure 8Comparison between modeled and observed characteristic landslide scaling relationships. (a) Magnitude–frequency relationship. The grey and black lines represent the best-fitting inverse gamma distribution (Eq. 19) of the landslide activity mapped before 1974 and between 1974 and 2012, respectively (Larsen and Montgomery2012). The fitting parameters are, respectively, (al=768×10-6, sl=-32.6×10-6, ρl=1.27) and (al=6100×10-6, sl=-311×10-6, ρ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 γ=1.32±0.02 and 10logα=-0.49±0.06. The red dots represent the geometry of the landslides simulated with HyLands.


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 riverbed 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 that 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).

Figure 9Synthetic model run showing landscape evolution to steady state followed by an intense landsliding period of 100 years. (a–c) Time slices showing evolution of the landscape to steady state, before the landslide period. The upper-left subplots show the evolution of topography through time. The upper-right subplots show the evolution of sediment thickness (H) through time. On both subplots, the blue line represents the location of the river plotted in the lower subplots. These lower subplots show the topographic and bedrock elevation (red and black line, respectively). The difference between the topographic elevation and the elevation of the bedrock represents the sediment thickness. With respect to total elevation, sediment thickness is small, so sediment thickness (orange line) is also plotted against the righthand y axis. The grey shaded area represents bedrock underlying the river profile. (d–f) Time slices showing the landslide period where intense landsliding is occurring over a period of 100 years. The upper-left subplots show the landslide activity. The location of landslides is indicated with black diamonds. The colors represent the square root of the landslide erosion () and deposition (+) during the presented time step. The upper-right subplots show the evolution of sediment thickness (H) through time. On both subplots, the blue line represents the location of the river plotted in the lower subplots. These lower subplots show the topographic and bedrock elevation (red and black line, respectively) as well as the volume occupied by sediments and water (orange and blue shaded area, respectively). 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).


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 (Sect. 3), i.e., with a single open boundary cell in the southwestern corner. To simulate more realistic landscape scales, we use a domain of 75 by 75 cells 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. 9a–c). Model behavior 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 time steps, the drainage network establishes, and the landscape gradually approaches a steady state with uniform sediment thickness across the entire landscape.

In the second stage (Fig. 9d–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 (tLS=2×103 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 (Ffhill=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 riverbed can change due to landslides and landslide-derived sediment rearranging the valley-bottom topography.

Immediately after the intense landsliding period, the trunk stream of the drainage network is still choked with sediment, and landslide dams are abundant (Fig. 10a–f). In the first few thousand years following the intense landsliding period, the lakes gradually fill 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 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 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, 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 base-level signal.

There are two distinct mechanisms for the generation of irregularities in the channel profiles: drainage rerouting 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 Figs. 9 and 10). The latter creates variability in the magnitude but not the direction of the bedrock slope. The drainage rerouting mechanism dominates in the simulations presented here and results in the formation of epigenetic river gorges (Fig. 10). Epigenetic gorges are characterized by rivers incising into the bedrock of former valley walls due to the blockage of the formal channel by landslide-derived sediment (Ouimet et al.2008).

Figure 10Synthetic model run showing recovery from intense landsliding illustrated in Fig. 9. (a–f) Time slices showing reestablishment of the landscape to steady state. Bedrock bumps created by landslide-induced drainage redirection are eroded, and the channel reattains its smoothly concave-up, steady-state configuration. For a detailed description of subplots and labels, see Fig. 9.


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 Anderson2011; Glade et al.2019). Much effort has been devoted to understanding the relationship between fluvial erosion efficiency and climate variability both through theoretical developments (Tucker2004; Lague2014) and observations (DiBiase and Whipple2011; Ferrier et al.2013). A main finding 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.2020; Glade et al.2019). However, the role of landslides in long-term landscape evolution, 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 (Figs. 9 and 10) 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 but also it results in the formation of bedrock knickpoints and associated retreating incision waves. The HyLands output thereby corroborates 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; Lague2010). However, like SPACE, HyLands does not simulate the effect of increased bedrock incision efficiency due to mobile sediment acting as eroding tools (Sklar and Dietrich2004). 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 modeling 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 discharge Qs (see e.g., Gasparini et al.2007; Hobley et al.2011; Zhang et al.2018).

Like SPACE, HyLands does not include process-based approaches (Kean and Smith2004; Wobus et al.2006; Davy and Lague2009; Coulthard et al.2013), simplified adjustment rules (Lague2010; Yanites2018) 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). It has however been shown that river width might vary as a function of sediment discharge or under varying tectonic configurations (Amos and Burbank2007; Turowski et al.2009; Pfeiffer et al.2017). Recent work using a two-dimensional hydro-sedimentary numerical model (Davy et al.2017) based on the Saint-Venant equations has shown that river reorganization and narrowing after landslide events might strongly increase sediment transport capacity and alter sediment evacuation time after large landslide events (Croissant et al.2017). While simulating dynamic river width reorganization at the landscape scale is currently not possible over longer timescales due to computational limitations, generic approximations for 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 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 surfaces (Jeandet et al.2019), the use of more advanced rupture plane geometries has been proposed. Gallen et al. (2015) for example propose the use of concave-upward rupture planes to simulate coseismic landsliding. However, their approach is based on the statistical aggregation of one-dimensional slope-stability solutions and therefore simplifies the three-dimensional complexity of the topographic surface. 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 Dietrich1994; 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 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 over-parameterization of the model (Van Rompaey and Govers2002).

Our probabilistic sliding mechanism neglects seismic or hydrological landslide triggers (e.g., Keefer1984, 2002; Keefer and Larsen2007; 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 (see Croissant et al.2019). The probability of coseismic 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 nonlinear, 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 simulated with HyLands with runout distances simulated with higher-complexity models (Iverson2000; George2011; Zhou et al.2020) or medium-complexity approaches such as RAMMS or Flow-R (Horton et al.2013; Fan et al.2017). 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 scales and timescales.

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 scales. Therefore, a range of datasets can be used to constrain model parameters. We identify three main categories of potential calibration data.

  1. First are 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 (see Fig. 7).

  2. Second are 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 decadal timescales (Hovius et al.1997; Larsen and Montgomery2012). Data on sediment redistribution following landslide events are however more difficult to collect (Fan et al.2019). 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 mobilization rates and predict landslide sediment production rates in regions which are otherwise difficult to access.

  3. Third are catchment-averaged cosmogenic radionuclide (CRN)-derived erosion rates. CRN data have been used to calibrate river incision models that explicitly integrate the stochastic nature of fluvial incision over time (DiBiase and Whipple2011; Scherler et al.2017; Campforts et al.2020). However, CRN data are sensitive to landslide activity (Niemi et al.2005; Yanites et al.2009; Wang et al.2017; Tofelde et al.2018), and the calibration of stochastic river incision models has been shown to be sensitive to landslide activity (Campforts et al.2020). HyLands directly simulates stochastic landsliding and hence enables explicit simulation of the impact of landsliding on CRN-derived erosion rates.

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.

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 (Korup2005; Korup et al.2007). By altering sediment discharge, landslides 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). A particular question which remains open for debate is the way in which landslides influence the evolution of a landscape to steady state. Although the stochastic nature of landslides prevents landscapes from evolving towards perfectly time- and space-invariant topographies, landslide-influenced landscapes will evolve towards a quasi-steady state if external drivers such as climate and tectonics remain constant. Though the study of steady-state landscapes improves our mechanistic understanding of landscape evolution, an even more interesting and challenging question would be to study the impact of landslides on the dynamic evolution of a landscape towards such a steady state. The latter problem is more relevant for most real-world landscapes, which are thought to be transient rather than in steady state (Mudd2017).

Second, HyLands can be used to evaluate the response of a landscape to a major landslide triggering event and to understand the timescales over which landslide-derived sediments are exported from the landscape (Wang et al.2015; Li et al.2016; Schwanghart et al.2016; Robinson et al.2017; Roback et al.2018). 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 (Sect. 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-discharge-dependent river 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 Lague2009; 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 Scherler2014), 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 evolution and response to 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: (last access: 25 August 2020). The exact version of the software code used to produce the results presented in this paper is archived on Zenodo (, Campforts2020a). Documentation, installation instructions and software dependencies for the HyLands project can be found at (last access: 25 August 2020). Detailed scripts and user manuals for the simulations illustrated in this paper can be found at (last access: 25 August 2020), which is also archived on Zenodo: (Campforts2020i). HyLands is platform independent and requires MATLAB 2018a or higher and the Image Processing Toolbox. The HyLands modeling framework is distributed under a MIT open-source license.

Data availability

Digital elevation models used in this paper are derived from the 30 m SRTM v3 DEM (Farr et al.2007). The resampled DEM is available through

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.


Benjamin Campforts received funding from the Research Foundation Flanders, FWO grant agreement no. 12Z6518N. Charles M. Shobe received funding from the European Union's Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement no. 833132. Philippe Steer and Dimitri Lague have received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement no. 803721).

Financial support

This research has been supported by the FWO – Research Foundation Flanders (grant no. FWO – Postdoctoral fellowship/12Z6518N), the H2020 Marie Skłodowska-Curie Actions (grant no. STRATASCAPE (833132)) and the H2020 European Research Council (grant no. FEASIBLe (803721)).

The article processing charges for this open-access
publication were covered by a Research
Centre of the Helmholtz Association.

Review statement

This paper was edited by Andrew Wickert and reviewed by Alexander Densmore and one anonymous referee.


Amos, C. B. and Burbank, D. W.: Channel width response to differential uplift, J. Geophys. Res., 112, F02010,, 2007. a

Andrews, D. J. and Hanks, T. C.: Scarp degraded by linear diffusion: Inverse solution for age, J. Geophys. Res., 90, 10193,, 1985. a

Armitage, J. J., Whittaker, A. C., Zakari, M., and Campforts, B.: Numerical modelling of landscape and sediment flux response to precipitation rate change, Earth Surf. Dynam., 6, 77–99,, 2018. a

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, J. Geophys. Res.-Earth, 113, 1–16,, 2008. a

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, J. Geophys. Res., 115, F03013,, 2010. a

Beaumont, C., Fullsack, P., and Hamilton, J.: Erosional control of active compressional orogens, in: Thrust Tectonics, Springer Netherlands, Dordrecht, 1–18,, 1992. a, b

Beer, A. R., Turowski, J. M., and Kirchner, J. W.: Spatial patterns of erosion in a bedrock gorge, J. Geophys. Res.-Earth, 122, 191–214,, 2017. a

Beven, K. and Freer, J.: Equifinality, data assimilation, and uncertainty estimation in mechanistic modelling of complex environmental systems using the GLUE methodology, J. Hydrol., 249, 11–29,, 2001. a

Broeckx, J., Rossi, M., Lijnen, K., Campforts, B., Poesen, J., and Vanmaercke, M.: Landslide mobilization rates: A global analysis and model, Earth-Sci. Rev., 201, 102972,, 2020. a, b, c, d, e

Burbank, D., Meigs, A., and Brozović, N.: Interactions of growing folds and coeval depositional systems, Basin Res., 8, 199–223,, 1996. a

Burbank, D. W.: Rates of erosion and their implications for exhumation, Mineral. Mag., 66, 25–52,, 2002. a

Burbank, D. W. and Anderson, R. S.: Tectonic Geomorphology, John Wiley & Sons, Ltd, Chichester, UK,, 2011. a

Campforts, B.: BCampforts/pub_hylands_campforts_etal_GMD: pub_hylands_campforts_etal_GMD (Version v1.03), Zenodo,, 2020a. a

Campforts, B.: HyLands: No Landsliding – Mixed fluvial incision, HyLands 1.0: a Hybrid Landscape evolution model to simulate the impact of landslides and landslide-derived sediment on landscape evolution., 2020b. a, b

Campforts, B.: HyLands: No Landsliding – Transport-limited, HyLands 1.0: a Hybrid Landscape evolution model to simulate the impact of landslides and landslide-derived sediment on landscape evolution., 2020c. a, b

Campforts, B.: HyLands: No Landsliding – Mixed fluvial incision, HyLands 1.0: a Hybrid Landscape evolution model to simulate the impact of landslides and landslide-derived sediment on landscape evolution., 2020d. a, b

Campforts, B.: HyLands: Landsliding – Synthetic: Before intense LS period, HyLands 1.0: a Hybrid Landscape evolution model to simulate the impact of landslides and landslide-derived sediment on landscape evolution., 2020e. a, b

Campforts, B.: HyLands: Landsliding – Synthetic: During intense LS period, HyLands 1.0: a Hybrid Landscape evolution model to simulate the impact of landslides and landslide-derived sediment on landscape evolution., 2020f. a, b

Campforts, B.: HyLands: Landsliding – Synthetic: After intense LS period, HyLands 1.0: a Hybrid Landscape evolution model to simulate the impact of landslides and landslide-derived sediment on landscape evolution., 2020g. a, b

Campforts, B.: HyLands: Landsliding – Real DEM, Namche-Barwa, HyLands 1.0: a Hybrid Landscape evolution model to simulate the impact of landslides and landslide-derived sediment on landscape evolution., 2020h. a, b

Campforts, B.: BCampforts/topotoolbox: topotoolbox-v2.4-HyLands-v1.0 (Version v2.4-HyLands-v1.0), Zenodo,, 2020i. a

Campforts, B. and Govers, G.: Keeping the edge: A numerical method that avoids knickpoint smearing when solving the stream power law, J. Geophys. Res.-Earth, 120, 1189–1205,, 2015. a

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 Planet. Sc. Lett., 439, 143–157,, 2016. a, b

Campforts, B., Schwanghart, W., and Govers, G.: Accurate simulation of transient landscape evolution by eliminating numerical diffusion: the TTLEM 1.0 model, Earth Surf. Dynam., 5, 47–66,, 2017. a, b

Campforts, B., Vanacker, V., Herman, F., Vanmaercke, M., Schwanghart, W., Tenorio, G. E., Willems, P., and Govers, G.: Parameterization of river incision models requires accounting for environmental heterogeneity: insights from the tropical Andes, Earth Surf. Dynam., 8, 447–470,, 2020. a, b, c

Carretier, S., Martinod, P., Reich, M., and Godderis, Y.: Modelling sediment clasts transport during landscape evolution, Earth Surf. Dynam., 4, 237–251,, 2016. a, b, c, d, e

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,, 2018. a

Champel, B.: Growth and lateral propagation of fault-related folds in the Siwaliks of western Nepal: Rates, mechanisms, and geomorphic signature, J. Geophys. Res., 107, 2111,, 2002. a, b, c, d

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,, 2007. a, b, c, d

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

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 Surf. Proc. Land., 38, 1897–1906,, 2013. a

Croissant, T., Lague, D., Steer, P., and Davy, P.: Rapid post-seismic landslide evacuation boosted by dynamic river width, Nat. Geosci., 10, 680–684,, 2017. a, b

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

Culmann, K.: Die Graphische Statike, Verlag von Meyer & Zeller, Zurich, 1875. a, b, c

Dahlquist, M. P., West, A. J., and Li, G.: Landslide-driven drainage divide migration, Geology, 46, 403–406,, 2018. a

Davy, P. and Lague, D.: Fluvial erosion/transport equation of landscape evolution models revisited, J. Geophys. Res.-Sol. Ea., 114, 1–16,, 2009. a, b, c, d, e, f, g, h

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

Densmore, A. L., Ellis, M. A., and Anderson, R. S.: Landsliding and the evolution of normal-fault-bounded mountains, J. Geophys. Res.-Sol. Ea., 103, 15203–15219,, 1998. a, b, c, d, e, f, g, h, i, j, k

DiBiase, R. A. and Whipple, K. X.: The influence of erosion thresholds and runoff variability on the relationships among topography, climate, and erosion rate, J. Geophys. Res., 116, F04036,, 2011. a, b

Dussauge, C., Grasso, J.-R., and Helmstetter, A.: Statistical analysis of rockfall volume distributions: Implications for rockfall dynamics, J. Geophys. Res.-Sol. Ea., 108,, 2003. a

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,, 2013. a, b, c, d, e

Fan, L., Lehmann, P., McArdell, B., and Or, D.: Linking rainfall-induced landslides with debris flows runout patterns towards catchment scale hazard assessment, Geomorphology, 280, 1–15,, 2017. a

Fan, X., Scaringi, G., Korup, O., West, A. J., Westen, C. J., Tanyas, H., Hovius, N., Hales, T. C., Jibson, R. W., Allstadt, K. E., Zhang, L., Evans, S. G., Xu, C., Li, G., Pei, X., Xu, Q., and Huang, R.: Earthquake‐Induced Chains of Geologic Hazards: Patterns, Mechanisms, and Impacts, Rev. Geophys., 57, 421–503,, 2019. a

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, Rev. Geophys., 45, RG2004,, 2007. a, b, c, d

Ferrier, K. L., Huppert, K. L., and Perron, J. T.: Climatic control of bedrock river incision, Nature, 496, 206–209,, 2013. a

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, Geol. Soc. Am. Bull., 120, 142–155,, 2008. a, b

Furbish, D. J. and Roering, J. J.: Sediment disentrainment and the concept of local versus nonlocal transport on hillslopes, J. Geophys. Res.-Earth, 118, 937–952,, 2013. a

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,, 2015. a, b

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, J. Geophys. Res., 112, F03S09,, 2007. a, b

George, D. L.: Adaptive finite volume methods with well-balanced Riemann solvers for modeling floods in rugged terrain: Application to the Malpasset dam-break flood (France, 1959), Int. J. Numer. Meth. Fl., 66, 1000–1018,, 2011. a

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,, 2019. a, b, c, d

Guns, M. and Vanacker, V.: Shifts in landslide frequency-area distribution after forest conversion in the tropical Andes, Anthropocene, 6, 75–85,, 2014. a, b

Guzzetti, F., Malamud, B. D., Turcotte, D. L., and Reichenbach, P.: Power-law correlations of landslide areas in central Italy, Earth Planet. Sc. Lett., 195, 169–183,, 2002. a, b

Guzzetti, F., Reichenbach, P., Ardizzone, F., Cardinali, M., and Galli, M.: Estimating the quality of landslide susceptibility models, Geomorphology, 81, 166–184,, 2006. a

Hancock, G. S. and Anderson, R. S.: Numerical modeling of fluvial strath-terrace formation in response to oscillating climate, GSA Bulletin, 114, 1131–1142,<1131:NMOFST>2.0.CO;2, 2002. a

Hobley, D. E., Sinclair, H. D., Mudd, S. M., and Cowie, P. A.: Field calibration of sediment flux dependent river incision, J. Geophys. Res.-Earth, 116, F04017,, 2011. a

Horton, P., Jaboyedoff, M., Rudaz, B., and Zimmermann, M.: Flow-R, a model for susceptibility mapping of debris flows and other gravitational hazards at a regional scale, Nat. Hazards Earth Syst. Sci., 13, 869–885,, 2013. a

Hovius, N., Stark, C. P., and Allen, P. A.: Sediment flux from a mountain belt derived by landslide mapping, Geology, 25, 231–234,<0231:SFFAMB>2.3.CO;2, 1997. a, b, c

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, J. Geol., 108, 73–89,, 2000. a

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 Planet. Sc. Lett., 304, 347–355,, 2011. a

Howard, A. D. and Kerby, G.: Channel changes in badlands, GSA Bulletin, 94, 739–752,<739:CCIB>2.0.CO;2, 1983. a, b

Iverson, R. M.: Landslide triggering by rain infiltration, Water Resour. Res., 36, 1897–1910,, 2000. a, b

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,, 2016. a

Jeandet, L., Steer, P., Lague, D., and Davy, P.: Coulomb Mechanics and Relief Constraints Explain Landslide Size Distribution, Geophys. Res. Lett., 46, 4258–4266,, 2019. a, b, c, d, e, f, g

Kean, J. W. and Smith, J. D.: Flow and boundary shear stress in channels with woody bank vegetation, Water Sci. Appl., 8, 237–252, 2004. a

Keefer, D. K.: Landslides caused by earthquakes, GSA Bulletin, 95, 406–421, 1984. a

Keefer, D. K.: Investigating landslides caused by earthquakes – A historical review, Surv. Geophys., 23, 473–510,, 2002. a

Keefer, D. K. and Larsen, M. C.: Assessing Landslide Hazards, Science, 316, 1136–1138,, 2007. a

King, G. E., Herman, F., Lambert, R., Valla, P. G., and Guralnik, B.: Multi-OSL-thermochronometry of feldspar, Quatern. Geochronol., 33, 76–87,, 2016. a, b

Kirby, E. and Whipple, K. X.: Expression of active tectonics in erosional landscapes, J. Struct. Geol., 44, 54–75,, 2012. a

Korup, O.: Large landslides and their effect on sediment flux in South Westland, New Zealand, Earth Surf. Proc. Land., 30, 305–323,, 2005. a, b

Korup, O.: Rock type leaves topographic signature in landslide-dominated mountain ranges, Geophys. Res. Lett., 35, L11402,, 2008. a

Korup, O., Clague, J. J., Hermanns, R. L., Hewitt, K., Strom, A. L., and Weidinger, J. T.: Giant landslides, topography, and erosion, Earth Planet. Sc. Lett., 261, 578–589,, 2007. a, b, c

Korup, O., Densmore, A. L., and Schlunegger, F.: The role of landslides in mountain range evolution, Geomorphology, 120, 77–90,, 2010. a

Lague, D.: Reduction of long-term bedrock incision efficiency by short-term alluvial cover intermittency, J. Geophys. Res.-Earth, 115, F02011,, 2010. a, b, c

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

Lague, D., Hovius, N., and Davy, P.: Discharge, discharge variability, and the bedrock channel profile, J. Geophys. Res.-Earth, 110, F04006,, 2005. a

Larsen, I. J. and Montgomery, D. R.: Landslide erosion coupled to tectonics and river incision, Nat. Geosci., 5, 468–473,, 2012. a, b, c, d, e, f, g, h, i

Larsen, I. J., Montgomery, D. R., and Korup, O.: Landslide erosion controlled by hillslope material, Nat. Geosci., 3, 247–251,, 2010. a, b

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, J. Geophys. Res.-Earth, 121, 703–724,, 2016. a

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, Eng. Geol., 97, 32–41,, 2008. a

Malamud, B. D. and Turcotte, D. L.: Self-organized criticality applied to natural hazards, Nat. Hazards, 20, 93–116,, 1999. a, b

Malamud, B. D., Turcotte, D. L., Guzzetti, F., and Reichenbach, P.: Landslide inventories and their statistical properties, Earth Surf. Proc. Land., 29, 687–711,, 2004. a

Marc, O., Hovius, N., Meunier, P., Uchida, T., and Hayashi, S.: Transient changes of landslide rates after earthquakes, Geology, 43, 883–886,, 2015. a

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 Surf. Dynam., 6, 903–922,, 2018. a, b

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 Surf. Dynam., 7, 107–128,, 2019. a

Meunier, P., Hovius, N., and Haines, A. J.: Regional patterns of earthquake-triggered landslides and their relation to ground motion, Geophys. Res. Lett., 34, L20408,, 2007. a, b

Milliman, J. D. and Meade, R. H.: World-Wide Delivery of River Sediment to the Oceans, J. Geol., 91, 1–21,, 1983. a

Montgomery, D. R. and Dietrich, W. E.: A physically based model for the topographic control on shallow landsliding, Water Resour. Res., 30, 1153–1171,, 1994. a, b

Montgomery, D. R. and Gran, K. B.: Downstream variations in the width of bedrock channels, Water Resources Research, 37, 1841–1846,, 2001. a

Mudd, S. M.: Detection of transience in eroding landscapes, Earth Surf. Proc. Land., 42, 24–41,, 2017. a

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 Planet. Sc. Lett., 237, 480–498,, 2005. a

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), Geol. Soc. Ame. Bull., 119, 1462–1476,, 2007. a, b, c, d, e, f

Ouimet, W. B., Whipple, K. X., Crosby, B. T., Johnson, J. P., and Schildgen, T. F.: Epigenetic gorges in fluvial landscapes, Earth Surf. Proc. Land., 33, 1993–2009,, 2008. a

Page, M. J., Reid, L. M., and Lynn, I. H.: New Zealand Hydrological Society Sediment production from Cyclone Bola landslides, Waipaoa catchment, J. Hydrol., 38, 289–308, 1999. a

Paola, C. and Voller, V. R.: A generalized Exner equation for sediment mass balance, J. Geophys. Res.-Earth, 110, 1–8,, 2005. a

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, Sci. Rep., 6, 34438,, 2016. a

Pelletier, J. D.: Minimizing the grid-resolution dependence of flow-routing algorithms for geomorphic applications, Geomorphology, 122, 91–98,, 2010. a

Pfeiffer, A. M., Finnegan, N. J., and Willenbring, J. K.: Sediment supply controls equilibrium channel geometry in gravel rivers, P. Natl. Acad. Sci. USA, 114, 3346–3351,, 2017. a

Roback, K., Clark, M. K., West, A. J., Zekkos, D., Li, G., Gallen, S. F., Chamlagain, D., and Godt, J. W.: The size, distribution, and mobility of landslides caused by the 2015 Mw7.8 Gorkha earthquake, Nepal, Geomorphology, 301, 121–138,, 2018. a

Robinson, T. R., Rosser, N. J., Densmore, A. L., Williams, J. G., Kincey, M. E., Benjamin, J., and Bell, H. J. A.: Rapid post-earthquake modelling of coseismic landslide intensity and distribution for emergency response decision support, Nat. Hazards Earth Syst. Sci., 17, 1521–1540,, 2017. a

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,, 1999. a, b

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, J. Geophys. Res.-Earth, 122, 1389–1429,, 2017. a

Schwanghart, W. and Scherler, D.: Short Communication: TopoToolbox 2 – MATLAB-based software for topographic analysis and modeling in Earth surface sciences, Earth Surf. Dynam., 2, 1–7,, 2014. a, b, c

Schwanghart, W., Bernhardt, A., Stolle, A., Hoelzmann, P., Adhikari, B. R., Andermann, C., Tofelde, S., Merchel, S., Rugel, G., Fort, M., and Korup, O.: Repeated catastrophic valley infill following medieval earthquakes in the Nepal Himalaya, Science, 351, 147–150,, 2016. a

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

Shobe, C. M., Tucker, G. E., and Anderson, R. S.: Hillslope-derived blocks retard river incision, Geophys. Res. Lett., 43, 5070–5078,, 2016. a, b, c, d

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, Geosci. Model Dev., 10, 4577–4604,, 2017. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o, p, q

Shobe, C. M., Tucker, G. E., and Rossi, M. W.: Variable‐Threshold Behavior in Rivers Arising From Hillslope‐Derived Blocks, J. Geophys. Res.-Earth, 123, 1931–1957,, 2018. a, b

Sidle, R. C. and Ochiai, H.: Landslides: Processes, Prediction, and Land Use, Water Resources Monograph, American Geophysical Union, Washington, D. C.,, 2006. a

Sklar, L. S. and Dietrich, W. E.: A mechanistic model for river incision into bedrock by saltating bed load, Water Resour. Res., 40, W06301,, 2004. a

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, J. Geophys. Res.-Sol. Ea., 108, 2388,, 2003. a

Stark, C. P. and Hovius, N.: The characterization of landslide size distributions, Geophys. Res. Lett., 28, 1091–1094,, 2001. a, b

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

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,, 2018. a

Tofelde, S., Duesing, W., Schildgen, T. F., Wickert, A. D., Wittmann, H., Alonso, R. N., and Strecker, M.: Effects of deep-seated versus shallow hillslope processes on cosmogenic 10Be concentrations in fluvial sand and gravel, Earth Surf. Proc. Land., 43, 3086–3098,, 2018. a

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 Surf. Proc. Land., 29, 185–205,, 2004. a

Tucker, G. E. and Hancock, G. R.: Modelling landscape evolution, Earth Surf. Proc. Land., 35, 28–50,, 2010. a

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, J. Geophys. Res., 112, F04006,, 2007. a

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, J. Geophys. Res., 114, F03016,, 2009. a

Van Asch, T., Buma, J., and Van Beek, L.: A view on some hydrological triggering systems in landslides, Geomorphology, 30, 25–32,, 1999. a

Van Rompaey, A. J. J. and Govers, G.: Data quality and model complexity for regional scale soil erosion prediction, Int. J. Geogr. Inf. Sci., 16, 663–680,, 2002. a

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,, 2015. a

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 Surf. Proc. Land., 42, 2611–2622,, 2017. a

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

Whipple, K. X., Hancock, G. S., and Anderson, R. S.: River incision into bedrock: Mechanics and relative efficacy of plucking, abrasion, and cavitation, Geol. Soc. Am. Bull., 112, 490–503,<490:RIIBMA>2.0.CO;2, 2000. a

Willgoose, G., Bras, R. L., and Rodrigueziturbe, I.: Results from a New Model of River Basin Evolution, Earth Surf. Proc. Land., 16, 237–254,, 1991. a

Wobus, C., Whipple, K. X., Kirby, E., Snyder, N., Johnson, J., Spyropolou, K., Crosby, B., and Sheehan, D.: Tectonics from topography: Procedures, promise, and pitfalls, in: Tectonics, Climate, and Landscape Evolution, Geological Society of America, 398, 55–74,, 2006. a, b

Wyllie, D. C. and Mah, C. W.: Rock Slope Engineering, CRC Press,, 2017. a

Yanites, B. J.: The Dynamics of Channel Slope, Width, and Sediment in Actively Eroding Bedrock River Systems, J. Geophys. Res.-Earth, 123, 1504–1527,, 2018. a

Yanites, B. J., Tucker, G. E., and Anderson, R. S.: Numerical and analytical models of cosmogenic radionuclide dynamics in landslide-dominated drainage basins, J. Geophys. Res., 114, F01007,, 2009. a

Zhang, J., van Westen, C. J., Tanyas, H., Mavrouli, O., Ge, Y., Bajrachary, S., Gurung, D. R., Dhital, M. R., and Khanal, N. R.: How size and trigger matter: analyzing rainfall- and earthquake-triggered landslide inventories and their causal relation in the Koshi River basin, central Himalaya, Nat. Hazards Earth Syst. Sci., 19, 1789–1805,, 2019.  a

Zhang, L., Stark, C., Schumer, R., Kwang, J., Li, T., Fu, X., Wang, G., and Parker, G.: The Advective‐Diffusive Morphodynamics of Mixed Bedrock‐Alluvial Rivers Subjected to Spatiotemporally Varying Sediment Supply, J. Geophys. Res.-Earth, 123, 1731–1755,, 2018. a

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, 17, 1217–1231,, 2020. a

Short summary
Landslides shape the Earth’s surface and are a dominant source of terrestrial sediment. Rivers, then, act as conveyor belts evacuating landslide-produced sediment. Understanding the interaction among rivers and landslides is important to predict the Earth’s surface response to past and future environmental changes and for mitigating natural hazards. We develop HyLands, a new numerical model that provides a toolbox to explore how landslides and rivers interact over several timescales.