Steady-state forms of channel profiles shaped by debris flow and fluvial processes

Debris flows regularly traverse bedrock channels that dissect steep landscapes, but our understanding of bedrock erosion by debris flows and their impact on steepland morphology is still rudimentary. Quantitative models of steep bedrock channel networks are based on geomorphic transport laws designed to represent erosion by water-dominated flows. To quantify the impact of debris flow erosion on steep channel network form, it is first necessary to develop methods to estimate spatial variations in bulk debris flow properties (e.g., flow depth, velocity) throughout the channel network that can be integrated into landscape evolution models. Here, we propose and evaluate two methods to estimate spatial variations in bulk debris flow properties along the length of a channel profile. We incorporate both methods into a model designed to simulate the evolution of longitudinal channel profiles that evolve in response to debris flow and fluvial processes. To explore this model framework, we propose a general family of debris flow erosion laws where erosion rate is a function of debris flow depth and channel slope. Model results indicate that erosion by debris flows can explain the occurrence of a scaling break in the slope–area curve at low-drainage areas and that upper-network channel morphology may be useful for inferring catchment-averaged erosion rates in quasi-steady landscapes. Validating specific forms of a debris flow incision law, however, would require more detailed model–data comparisons in specific landscapes where input parameters and channel morphometry can be better constrained. Results improve our ability to interpret topographic signals within steep channel networks and identify observational targets critical for constraining a debris flow incision law.

1 Introduction

Debris flows are effective at transporting coarse sediment (May, 2002; May and Gresswell, 2003) and eroding bedrock (McCoy et al., 2013; McCoy, 2015; Stock and Dietrich, 2006) in steep valleys and low-order channels where fluvial sediment transport may be inhibited by low runoff magnitudes and increases in thresholds for incipient sediment motion (e.g., Lamb et al., 2008; Prancevic et al., 2014). The influence of debris flow erosion on bedrock channel morphology, including longitudinal channel profiles (Montgomery and Foufoula‐Georgiou, 1993), has been recognized for decades. Although debris flows traverse channel networks in many steep landscapes and are capable of eroding bedrock via abrasion and plucking (Hsu et al., 2008; Stock and Dietrich, 2006; McCoy et al., 2013), their relative importance over geologic time compared to other geomorphic processes and the extent to which they affect landscape form at larger scales remains unclear. The spatial extent and magnitude of debris flow erosion, for example, may be limited in terms of downstream extent due to a lack of mobility or erosive power on modest slopes. Additional work is needed to determine the controls on the magnitude of debris flow erosion within different parts of the drainage network and the ensuing implications for landscape form (Whipple et al., 2013).

Topographic signatures of geomorphic processes, which we define as quantitative connections between processes and the morphology of a landform, can be used to infer the presence and rates of geomorphic processes from topographic data (e.g., Tucker and Bras, 1998). Bedrock channels dominated by fluvial erosion develop longitudinal profiles described by a power law relationship between slope and drainage area (Flint, 1974; Morisawa, 1962; Hack, 1957). A deviation from this power law scaling relationship at small drainage areas (Fig. 1) has been interpreted as a topographic signature of debris flow incision (Montgomery and Foufoula‐Georgiou, 1993; Seidl et al., 1992; Sklar and Dietrich, 1998; Stock and Dietrich, 2006, 2003). Stock and Dietrich (2003) found that the shape of longitudinal channel profiles that experience both debris flow and fluvial erosion can be described by a family of curves with the general form

$$\begin{array}{}\text{(1)}& S={\displaystyle \frac{S_{\mathrm{df}}}{\mathrm{1}+a_{\mathrm{1}}{A}^{a_{\mathrm{2}}}}},\end{array}$$

where A denotes upstream drainage area and S_df, a₁, and a₂ are empirical coefficients (Fig. 1). Here, we use the term channel in a general sense to refer to an axis of concentrated erosion along valley bottoms, but which may or may not reside within banks made of deposited sediment. The coefficient S_df is related to the slope that the channel approaches at low-drainage areas, a₂ controls the power law relationship between slope and drainage area when A is large, and a₁ (which has units of $\mathrm{1}/({\mathrm{length}}^{\mathrm{2}}{)}^{a_{\mathrm{2}}}$) controls the sharpness of the transition from fluvial power law scaling at large drainage areas to relatively constant slopes in debris-flow-dominated reaches at smaller drainage areas. The above expression for channel slope can be rewritten as

$$\begin{array}{}\text{(2)}& S={\displaystyle \frac{S_{\mathrm{df}}}{\mathrm{1}+(A/A_{\mathrm{df}}{)}^{a_{\mathrm{2}}}}},\end{array}$$

which is advantageous because A_df has units of length² and can be interpreted as the drainage area at which the slope–area relationship transitions from a near-constant slope with decreasing drainage area to the standard power law relationship expected in the fluvial network. In this sense, A_df provides one metric for identifying the transition between the debris flow domain, with a characteristic gradient tending towards S_df, and the fluvial process domain.

Figure 1 Examples of channel profiles from the San Gabriel Mountains, California, USA, along with best fit curves of the form $S=S_{\mathrm{df}}/(\mathrm{1}+(A/A_{\mathrm{df}}{)}^{a_{\mathrm{2}}})$. The locations of red circles coincide with A_df. Channel profiles were extracted from catchments with different erosion rates, E. Watershed outlets are located (UTM 11S) at 396887 m E 3799338 m N (black line), 384971 m E 3799277 m N (blue line), and 417896 m E 3792642 m N (orange line). Catchment-averaged erosion rates, E, are from DiBiase et al. (2010). Slope and drainage area values were extracted from the channel network and separated into 100 logarithmically spaced bins. The binned slope was aggregated using the mean slope within each bin. The R² value for each fit is 0.95 (black), 0.97 (blue), and 0.97 (orange).

Past work demonstrates that the length of the channel network upstream of the debris flow to fluvial transition zone, which we roughly associate with A_df, increases with erosion rate in two landscapes where debris flows are known to regularly traverse steep channels, namely the San Gabriel Mountains (DiBiase et al., 2012) and the Oregon Coast Range (Penserini et al., 2017). These results from DiBiase et al. (2012) and Penserini et al. (2017) are consistent with the conceptual model proposed by Stock and Dietrich (2003) where the transition from a nearly linear debris-flow-dominated long profile to a concave-up fluvial-dominated long profile migrates out to larger drainage areas as the rock uplift rate increases. Therefore, there is support for the idea that steady-state bedrock channels eroded by debris flows not only have a unique morphology (or topographic signature) that distinguishes them from purely fluvial channels, but they may also record tectonic information.

These findings underscore the need to develop a quantitative framework that can be used to explore topographic signatures generated by debris flow erosion, assess the sensitivity of topographic signatures to climatic and tectonic forcing, and ultimately interpret these signatures to gain process-based insights about the evolution of steep landscapes. In particular, there is a need to understand the relative importance of fluvial and debris flow processes in setting the location and form of the morphologic transition associated with A_df (Fig. 1). For example, the location of this transition may change with rock uplift rate due to the dynamics of fluvial erosion alone because channel steepness can vary nonlinearly with rock uplift rate as a result of relationships between runoff variability and fluvial erosion thresholds (e.g., DiBiase and Whipple, 2011; Lague, 2014). Debris flow processes, in contrast, may exert a strong control on the location of the transition by setting the near-uniform slope that the channel approaches at small drainage areas and the sharpness of the transition between near-linear and concave-up channels. Quantifying the roles of debris flow and fluvial processes on controlling the location of this transition, as approximated by A_df, would aid in determining the benefits and limitations of using the morphology of debris-flow-dominated channels as a proxy for erosion rate in quasi-steady landscapes (Penserini et al., 2017).

In landscape evolution models, fluvial erosion is modeled based on empirical relationships between local terrain attributes such as slope and drainage area, readily computed from a digital elevation model (DEM) (Tucker and Bras, 1998). However, incorporation of a debris flow incision law into this type of local framework is challenging owing to the nonlocal controls on debris flow erosion including initiation conditions, non-steady flow velocity and the finite and variable runout distance of discrete flows. For example, the frequency at which discrete flows traverse different parts of the landscape has been shown to be a key factor in their ability to sculpt topography (Shelef and Hilley, 2016). We are aware of only one attempt to model the effects of debris flow erosion over geologic timescales (Stock and Dietrich, 2006) aimed at reproducing channel profiles with the characteristic change in slope-area scaling observed in natural environments (Montgomery and Foufoula‐Georgiou, 1993; Stock and Dietrich, 2003). Stock and Dietrich (2006) coupled empirical relationships for debris flow properties with a debris flow incision law based on inertial stress but emphasized the need for improved methods to calculate spatially varying bulk debris flow properties and additional studies to constrain a debris flow incision law.

Several subsequent studies have improved our understanding of the grain-scale processes that control debris flow incision rates and their relationship to bulk flow properties that are more amenable to measurement and model simulation. Hsu et al. (2008, 2014) used physical experiments in a rotating drum to suggest that debris flow erosion rates scale with bulk inertial stress, which can be cast as a function of flow shear rate, commonly assumed to be proportional to the ratio of depth-averaged flow velocity and flow depth. Similarly, McCoy et al. (2013) showed that while bed impact forces beneath erosive debris flows were broadly distributed and the result of discrete particle–bed impacts, these force distributions scaled with bulk flow properties such as flow depth. Following the erosion equation proposed by Sklar and Dietrich (2004) to account for bedrock erosion due to discrete particle impacts, McCoy (2012) used a series of discrete element simulations to show that erosion by steady granular flows likely scales approximately linearly with flow depth and in a strongly nonlinear way with bed slope. The utility of these relationships in a landscape evolution model, however, requires tractable simulation of the spatial and temporal variability in the properties of individual debris flows (e.g., depth, velocity, shear rate) throughout the channel network and integration of the effects of numerous debris flow events on channel evolution over geologic timescales.

Here, we address this gap by developing a nonlocal modeling framework to predict the evolution of a 1D channel profile eroded by both fluvial and debris flow processes. In this framework, the routing of debris flows down the channel profile as well as estimates of spatial and temporal variations in bulk debris flow properties are represented either through a process-based model that relies on a set of partial differential equations (Iverson and Denlinger, 2001) or a reduced-complexity approach motivated by Gorr et al. (2022) based on empirical relationships (Rickenmann, 1999). Calculations of spatial variations in debris flow properties from these two methodologies provide a robust foundation for utilizing relationships between bulk debris flow properties and particle impact forces (Hsu et al., 2008; McCoy et al., 2013) to estimate rates of bedrock incision by debris flows. We propose a general family of debris flow erosion laws, with erosion rate being a function of debris flow depth and channel slope, to illustrate how the proposed framework may be used to help constrain a debris flow erosion law and to explore model sensitivity. We examine the extent to which different erosion laws are capable of reproducing the relationship between slope and drainage area, as captured by Eq. (2), that has been observed in steep, debris-flow-prone landscapes and interpreted as a topographic signature of debris flows. The process-based and empirical routing approaches differ in their assumptions and complexity, and comparison yields insight into the relative importance of these differences as compared with the proposed form of a debris flow erosion law. We finish with a sensitivity analysis to explore controls on A_df and S_df, two metrics that capture basic aspects of debris-flow-dominated channel morphology.

2 Methods

2.1 Model framework

In the proposed 1D model framework, which is designed to simulate longitudinal channel profiles, the rate of change in elevation, z, with time, t, in a bedrock channel is driven by the rock uplift rate, U, fluvial erosion, E_f, and erosion by debris flows, E_df, according to

$$\begin{array}{}\text{(3)}& {\displaystyle \frac{\partial z}{\partial t}}=U-E_{\mathrm{f}}-E_{\mathrm{df}}.\end{array}$$

We solve Eq. (3) numerically on a one-dimensional grid with a uniform spacing of Δx=5 m and use the standard explicit forward Euler method for time stepping. We chose a grid spacing of 5 m since it is small enough to resolve changes in the longitudinal channel profile and also large enough to keep model run times, which increase with decreasing grid spacing, manageable.

Fluvial erosion is computed using the threshold-stochastic stream power incision model presented by Lague (2014),

$$\begin{array}{}\text{(4)}& E_{\mathrm{f}}=K{A}^{m_{\mathrm{s}}}{S}^{n_{\mathrm{s}}},\end{array}$$

where A denotes upstream drainage area as a function of distance, x, from the channel head; S denotes slope; and K, m_s, and n_s are empirical parameters that depend on relationships between discharge and channel width, w, hydraulic geometry and discharge variability, and grain size. Here, we assume that drainage area varies with distance from the channel head according to $A=A_{\mathrm{0}}+\mathrm{25}{x}^{\mathrm{5}/\mathrm{3}}$, with A₀=1000 m² and the exponent of $\mathrm{5}/\mathrm{3}$ chosen to be consistent with a Hack's exponent of $\mathrm{3}/\mathrm{5}$. We assume channel width increases with drainage area as w=k_wA^b, where k_w=0.05 and b=0.3. Motivated by the geomorphic importance of debris flows in the San Gabriel Mountains (Lavé and Burbank, 2004), parameters related to channel geometry, including k_w and b that are related to width–area scaling (Fig. B1), and fluvial incision were selected based on DiBiase and Whipple (2011) and typical ranges reported by Lague (2014). These parameter choices result in m_s=1.4 and n_s=2.33. Complete details on parameter choices for the stream power model are given in Appendix A. Unless noted otherwise, parameters used for the fluvial erosion model are listed in Table G1.

We propose a general formulation that can be used to estimate the erosion rate attributable to debris flows, E_df, at a point on the landscape, given information about the bulk properties of the flow. In this work, we assume that bulk properties of a debris flow for a given landscape position do not change. In other words, debris flow erosion is driven over time by repeatedly routing the same debris flow over the landscape. Motivated by observations that debris flow erosion rates scale with bulk inertial stress (Hsu et al., 2008, 2014), a function of shear rate, and that grain-scale bed–impact force distributions scale with flow depth (McCoy et al., 2013), it is reasonable to postulate an erosion law that includes debris flow depth and velocity. Since steady granular flows down inclined planes of increasing angles show a monotonically increasing relationship between slope angle and velocity (Silbert et al., 2001), slope may serve as a proxy for velocity. Here, we define E_df as a function of channel slope and debris flow depth, h. Debris flow depth varies with position along the channel profile and with time throughout the course of a debris flow event. Letting t₀ denote the time a debris flow begins moving over a given location along the channel profile and t_f the time when it has completely passed that location, then

$$\begin{array}{}\text{(5)}& E_{\mathrm{df}}=k_{\mathrm{df}}\underset{t_{\mathrm{0}}}{\overset{t_{\mathrm{f}}}{\int }}{S}^{\mathit{\alpha }}{h}^{\mathit{\beta }}\mathrm{\Phi }\mathrm{d}t,\end{array}$$

where the debris flow erodibility coefficient is defined as k_df=κ_dfF_df, κ_df is an empirical coefficient related to bedrock and flow properties (e.g., grain size), F_df is a term quantifying the frequency of debris flow, Φ is a threshold factor that reflects a reduction in incision when the debris flow is close to rest, and α and β are empirical exponents. We use this family of erosion laws to begin exploring the model framework we propose here. The model is designed in such a way that it would be straightforward to insert alternative erosion laws in the future. This formulation could also be extended to account for a distribution of representative debris flows with different properties given information about their relative recurrence. With the process-based debris flow routing model (Sect. 2.2), we compute time-varying flow properties required to determine a debris flow incision rate at each point along the channel profile using Eq. (5).

We also present a reduced-complexity routing algorithm, which closely follows the methodology presented by Gorr et al. (2022) to rapidly simulate debris flow runout for hazard assessment purposes, to compute spatial variations in bulk debris flow properties along a longitudinal channel profile (Sect. 2.3). In this approach, we use a set of empirical relationships that relate debris flow properties to topographic slope (Rickenmann, 1999) in order to estimate representative values for debris flow depth, h, at each point along the debris flow path, as well as the time it takes for the debris flow to pass over that point, t_f−t₀, which we denote as t_p. In other words, h varies along the channel profile, but we neglect variations in h that occur within individual debris flow events (i.e., the rise and fall in flow depth as a debris flow passes over a point on the landscape). In this case, we employ a debris flow erosion equation analogous to Eq. (5) that is simplified because h is constant for a given channel location,

$$\begin{array}{}\text{(6)}& E_{\mathrm{df}}=k_{\mathrm{df}}{t}_{\mathrm{p}}{S}^{\mathit{\alpha }}{h}^{\mathit{\beta }}\mathrm{\Theta },\end{array}$$

where Θ denotes a threshold factor that reflects a reduction in incision when the debris flow is close to rest. To assess the simplifying assumptions of this approach within the context of modeling the evolution of longitudinal channel profiles, we compare the morphology of modeled profiles using this reduced complexity algorithm with profiles generated using the process-based debris flow routing model. This comparison is limited, as described in the following sections in more detail, to cases where downstream changes in debris flow volume are assumed to be negligible. While we acknowledge this is not likely to be true in many natural settings (Santi et al., 2008; Santi and Morandi, 2013; Schürch et al., 2011), examining this end-member case allows for the most direct comparison between channel profiles produced by the model when using these two different debris flow routing methods.

2.2 Estimating debris flow incision with a process-based routing model

The initial step in computing the erosion rate attributable to debris flows is to determine the runout path of the debris flow as well as its bulk properties at different points along that path. The process-based debris flow routing model is based on a set of conservation laws for mass and momentum within a depth-averaged framework. This particular model formulation was chosen because it provides sufficient complexity to enable exploration of the links between flow properties and the morphology of the resulting channel profile. The governing equations represent the flow of a two-component mixture, solids suspended in a Newtonian fluid (Iverson and Denlinger, 2001), in a rectangular channel with variable width (Vázquez-Cendón, 1999):

$$\begin{array}{}\text{(7)}& {\displaystyle \frac{\partial h}{\partial t}}+{\displaystyle \frac{\partial \left(hv\right)}{\partial x}}=-{\displaystyle \frac{vh}{w}}{\displaystyle \frac{\partial w}{\partial x}},\text{(8)}& {\displaystyle }\begin{array}{rl}& {\displaystyle \frac{\partial \left(hv\right)}{\partial t}}+{\displaystyle \frac{\partial }{\partial x}}\left(h{v}^{\mathrm{2}}+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{g}_z{h}^{\mathrm{2}}\right)=\\ & g_{x}h-sgn\left(v\right)(\mathrm{1}-\mathit{\lambda })g_{z}h\mathit{\varphi }\left(I\right)-{\displaystyle \frac{\mathrm{2}v\mathit{\eta }{v}_{\mathrm{f}}}{\mathit{\rho }h}}-{\displaystyle \frac{v^{\mathrm{2}}h}{w}}{\displaystyle \frac{\partial w}{\partial x}}.\end{array}\end{array}$$

Here, w is the channel width; h is flow depth; v is velocity; ϕ(I) is the friction coefficient that depends on the inertial number (Jop et al., 2006); I, g_x, and g_z denote components of gravity in the downslope and slope-normal directions, respectively; $\mathit{\lambda }=p_{\mathrm{bed}}/\mathit{\rho }{g}_{z}h$ is the ratio of pore fluid pressure to total basal normal stress; v_f=0.5 is the fluid volume fraction; and η is the viscosity of the pore fluid. The first, second, and third source terms on the right-hand side of the momentum conservation equation (Eq. 8) account for variations in bed topography, frictional resistance associated with the solid phase of the flow, and viscous resistance associated with the fluid phase, respectively. The remaining source terms in the mass and momentum equations account for variations in channel width.

We assume that the ratio of pore fluid pressure to total basal normal stress decays with time since the debris flow entered the model domain, t, according to

$$\begin{array}{}\text{(9)}& \mathit{\lambda }={\mathit{\lambda }}_{\mathrm{0}}\left[\mathrm{1}-\mathrm{erfc}\left({\displaystyle \frac{\mathrm{2}h}{\sqrt{\mathrm{4}Dt}}}\right)\right],\end{array}$$

where λ₀=0.9 and D is the pore fluid pressure diffusivity. This approximation is consistent with an initially high pore fluid pressure shortly following initiation and subsequent linear diffusion of pore fluid pressure (Iverson and Denlinger, 2001) over time.

The friction coefficient is a function of the inertial number (Jop et al., 2006),

$$\begin{array}{}\text{(10)}& \mathit{\varphi }\left(I\right)={\mathit{\mu }}_s+({\mathit{\mu }}_{\mathrm{2}}-{\mathit{\mu }}_s)/(I_{\mathrm{0}}/I+\mathrm{1}),\end{array}$$

where $I=\dot{\mathit{\gamma }}{D}_{\mathrm{eff}}/(P/{\mathit{\rho }}_{\mathrm{s}}{)}^{\mathrm{0.5}}$, with P denoting the basal normal stress,; $\dot{\mathit{\gamma }}=\mathrm{2}v/h$ is the shear rate; ρ_s=2600 kg m⁻³ is the density of sediment; I₀=0.279 is a constant; μ_s=0.382; μ₂=0.644; and D_eff is a characteristic particle diameter. In this formulation, the friction coefficient increases with the inertial number and approaches μ₂ when I is large.

The governing equations are solved numerically on a grid with uniform spacing. We use a first-order, shock-capturing finite-volume method with a Harten–Lax–van Leer contact (HLLC) approximate Riemann solver (Toro, 2009) to compute the fluxes across each grid cell boundary (McGuire et al., 2016, 2017). Source terms are treated separately with an explicit, first-order forward Euler method for time stepping.

Debris flows enter the domain through the upper boundary, which can be conceptualized as the channel head, and are routed down the channel profile. We define a series of 20 ghost cells above the uppermost grid cell that effectively extend the model domain for the purpose of initializing a debris flow. Elevations of each ghost cell are determined by assuming that the slopes of all ghost cells are equal to the slope at the uppermost grid cell. Debris flows are initiated from a static pile of debris defined on the ghost cells. This procedure provides some time for debris to begin to flow before it enters the model domain, similar to what might be expected for debris flows that initiate in a colluvial hollow or gully upstream of a channel head. In nature, we expect debris flow volume to vary with drainage area as sediment is entrained and deposited along the runout path (Santi and Morandi, 2013; Schürch et al., 2011; Santi et al., 2008), but incorporating this effect into the source terms of the process-based routing model is beyond the scope of this study. When using the process-based debris flow routing model, we assume that debris flow volume is fixed and does not change along the flow path, although we do explore the effects of spatial variations in debris flow volume with the empirical routing approach described later. In addition, we perform a set of numerical experiments with the process-based routing model where we scale debris flow frequency with drainage area to account for an increase in the total volume of sediment transported by debris flows as drainage area increases. Regardless of which routing approach is used, however, we do not explicitly account for rock mass incorporated into the debris flow originating from bedrock incision. This sediment volume would be negligible compared to the total debris flow volume.

At each grid cell in the model domain (i.e., excluding ghost cells), the debris flow incision rate is computed using Eq. (5) based on the time-varying values of debris flow depth. More specifically, for a debris flow simulated over k time steps,

$$\begin{array}{}\text{(11)}& E_{\mathrm{df}}=k_{\mathrm{df}}{S}^{\mathit{\alpha }}\sum _{k=\mathrm{1}}^{k=n}{h}_k^{\mathit{\beta }}\mathrm{\Phi }\mathrm{\Delta }t,\end{array}$$

where Δt denotes the time step used when solving the flow equations, and we define the threshold factor, Φ, as Φ=1 when uh>0.01 m² s⁻¹ and Φ=0 otherwise (Fig. E1). To reduce computation time, the term ${\sum }_{k=\mathrm{1}}^{k=n}{h}_k^{\mathit{\beta }}\mathrm{\Phi }$ is not updated with each time step in the landscape evolution model. Small changes in topography, such as may occur during a single time step of the landscape evolution model, will not substantially affect flow mobility or spatial variations in flow depth along the runout path. Instead, we only route a debris flow down the channel profile to update ${\sum }_{k=\mathrm{1}}^{k=n}{h}_k^{\mathit{\beta }}\mathrm{\Phi }$ in the calculation of E_df whenever the channel slope has changed by 0.05 or more at any grid cell since the last time a debris flow was routed. We do, however, update E_df with every time step of the landscape evolution model to reflect changes in slope, S, since this requires little computation time compared with debris flow routing. This is one benefit of using slope as a proxy for velocity in the debris flow erosion law.

2.3 Estimating debris flow incision with an empirical routing model

We use a series of empirical relationships defined by Rickenmann (1999) to estimate representative values for debris flow depth, h, and passage time, t_p, at each point along the channel profile based on spatially variable estimates of debris flow volume, M; debris flow velocity, v; channel width, w; and topographic slope, S. We assume that debris flows initiate at or above the uppermost grid cell within the computational domain (i.e., the channel head), although their overall volume may change along the channel profile. We determine the downstream extent of debris flow runout by treating the debris flow as an idealized fluid with a prescribed yield strength, τ_y, and assuming that debris flow motion stops when shear stress at the base of the flow, τ=ρ_bgR_hsin θ, falls below τ_y (Whipple and Dunne, 1992; Gorr et al., 2022). Here, $R_h=wh/(w+\mathrm{2}h)$ denotes the hydraulic radius of the rectangular channel, θ is the channel slope angle, g=9.81 m s⁻² denotes gravitational acceleration, and ρ_b=1800 kg m⁻² is the bulk density of the debris flow. In practice, we determine R_h, t_p, and τ everywhere in the model domain, determine the downstream extent of debris flow runout, and then apply Eq. (6) to compute a non-zero value for E_df only along the debris flow travel path.

To begin, we specify debris flow volume passing through each grid cell as a function of upstream drainage area (A) according to $M=M_{\mathrm{0}}({\mathrm{10}}^{-\mathrm{6}}\cdot A{)}^{\mathit{\gamma }}$ (Santi and Morandi, 2013). This formulation assumes that debris flow volume increases downstream, reflecting entrainment of bed material or lateral inflow, but these volume changes can be neglected by setting γ=0. Debris flow volume may not increase monotonically along the runout path (Schürch et al., 2011), but the formulation proposed by Santi and Morandi (2013) provides a useful starting point for a general parameterization of downstream variations in debris flow volume, especially since there are data from a range of geographic regions to fit such a relationship. For example, Santi and Morandi (2013) demonstrate that the empirical coefficient, M₀, and exponent, γ, may vary considerably among landscapes. Santi and Morandi (2013) estimated M₀=3358 and γ=0.73 using data throughout the western and southwestern United States, M₀=10470 and γ=0.62 based on data from the Italian Alps, and M₀=18770 and γ=0.28 using data from the northwestern United States and southwestern Canada. Peak debris flow discharge can then be computed according to

$$\begin{array}{}\text{(12)}& Q=c_{\mathrm{1}}{M}^{c_{\mathrm{2}}},\end{array}$$

where c₁=0.135 and c₂=0.78 are empirical coefficients (Rickenmann, 1999). Noting that Q=wvh and using the relationship (Rickenmann, 1999)

$$\begin{array}{}\text{(13)}& v={\displaystyle \frac{\mathrm{1}}{\mathrm{3}\mathit{\mu }}}{\mathit{\rho }}_{\mathrm{b}}g{h}^{\mathrm{2}}S,\end{array}$$

where μ denotes the dynamic viscosity of the flow, it is possible to solve for flow depth,

$$\begin{array}{}\text{(14)}& h={\left({\displaystyle \frac{\mathrm{3}\mathit{\mu }{c}_{\mathrm{1}}{M}^{c_{\mathrm{2}}}}{{\mathit{\rho }}_{\mathrm{b}}gSw}}\right)}^{\mathrm{1}/\mathrm{3}}.\end{array}$$

Using the relationships between channel width, w, and area, A, and debris flow volume, M, and area, A, a representative flow depth for a given location along the channel profile can be written in terms of area and slope,

$$\begin{array}{}\text{(15)}& h={\left({\displaystyle \frac{\mathrm{3}\mathit{\mu }{c}_{\mathrm{1}}\left(M_{\mathrm{0}}\right({\mathrm{10}}^{-\mathrm{6}}\cdot A{)}^{\mathit{\gamma }}{)}^{c_{\mathrm{2}}}}{{\mathit{\rho }}_{\mathrm{b}}gS{k}_{\mathrm{w}}{A}^b}}\right)}^{\mathrm{1}/\mathrm{3}}.\end{array}$$

We treat this flow depth as a representative value for each point in the drainage network but acknowledge that it may overestimate flow depth because Eq. (12) is used to estimate peak debris flow discharge. We further define the passage time of the debris flow as

$$\begin{array}{}\text{(16)}& t_{\mathrm{p}}={\displaystyle \frac{M}{Q}}={\displaystyle \frac{M_{\mathrm{0}}({\mathrm{10}}^{-\mathrm{6}}\cdot A{)}^{\mathit{\gamma }}}{c_{\mathrm{1}}\left(M_{\mathrm{0}}\right({\mathrm{10}}^{-\mathrm{6}}\cdot A{)}^{\mathit{\gamma }}{)}^{c_{\mathrm{2}}}}}.\end{array}$$

Finally, we define the threshold factor, Θ, such that the debris flow incision rate decreases as the flow approaches the end of its travel path and the shear stress at the base of the flow approaches the yield strength. Specifically,

$$\begin{array}{}\text{(17)}& \mathrm{\Theta }=\mathrm{1}-{\displaystyle \frac{{\mathit{\tau }}_{\mathrm{y}}}{\mathit{\tau }}}.\end{array}$$

The debris flow erosion rate can then be determined according to Eq. (6).

In this study, we fix all model parameters within a given simulation. As such, the channel profiles that develop can be thought of as reflecting the morphology of a channel shaped by the repeated impacts of a characteristic debris flow. Future studies could explore the effects of debris flows characterized by a distribution of parameters to better reflect natural variations in flow properties.

2.4 Simplified analytical solution

When using empirical relationships to determine flow properties along the debris flow runout path, we can derive an approximate analytical solution for the slope of the upper, debris-flow-dominated reach of the channel at steady state. We begin by considering debris flow erosion as quantified by Eqs. (6), (15), and (16). To arrive at an analytical solution, we then make several simplifying assumptions. First, we assume that fluvial erosion is negligible. Second, we assume that the channel is sufficiently steep so that shear stress at the base of the debris flow greatly exceeds the yield strength (i.e., τ>>τ_y). This implies that debris flows always traverse the entire channel reach that we are modeling and that it is reasonable to neglect the entertainment threshold (i.e., Θ=1). Enforcing the condition that the channel profile has reached a steady state yields

$$\begin{array}{}\text{(18)}& \begin{array}{rl}{\displaystyle \frac{\partial z}{\partial t}}& =\mathrm{0}=U-k_{\mathrm{df}}{S}^{\mathit{\alpha }}\left({\displaystyle \frac{M_{\mathrm{0}}({\mathrm{10}}^{-\mathrm{6}}\cdot A{)}^{\mathit{\gamma }}}{c_{\mathrm{1}}\left(M_{\mathrm{0}}\right({\mathrm{10}}^{-\mathrm{6}}\cdot A{)}^{\mathit{\gamma }}{)}^{c_{\mathrm{2}}}}}\right)\\ & {\left({\displaystyle \frac{\mathrm{3}\mathit{\mu }{c}_{\mathrm{1}}\left(M_{\mathrm{0}}\right({\mathrm{10}}^{-\mathrm{6}}\cdot A{)}^{\mathit{\gamma }}{)}^{c_{\mathrm{2}}}}{\mathit{\rho }gS{k}_{\mathrm{w}}{A}^b}}\right)}^{\mathit{\beta }/\mathrm{3}}\end{array}.\end{array}$$

Solving for slope as a function of drainage area, we obtain

$$\begin{array}{}\text{(19)}& S={\left({\displaystyle \frac{U}{k_{\mathrm{df}}{\mathit{\lambda }}_{\mathrm{1}}{\mathit{\lambda }}_{\mathrm{2}}{\mathit{\lambda }}_{\mathrm{3}}{\mathit{\lambda }}_{\mathrm{4}}}}\right)}^{\frac{\mathrm{1}}{\mathit{\alpha }-\mathit{\beta }/\mathrm{3}}}{A}^N,\end{array}$$

where

$$\begin{array}{}\text{(20)}& N={\displaystyle \frac{\mathit{\beta }/\mathrm{3}(b-\mathit{\gamma }{c}_{\mathrm{2}})+\mathit{\gamma }(c_{\mathrm{2}}-\mathrm{1})}{\mathit{\alpha }-\mathit{\beta }/\mathrm{3}}},\end{array}$$

and λ₁, λ₂, λ₃, λ₄ are given by

$$\begin{array}{}\text{(21)}& {\displaystyle }{\mathit{\lambda }}_{\mathrm{1}}=M_{\mathrm{0}}^{c_{\mathrm{2}}(\mathit{\beta }/\mathrm{3}-\mathrm{1})+\mathrm{1}},\text{(22)}& {\displaystyle }{\mathit{\lambda }}_{\mathrm{2}}={\mathrm{10}}^{-\mathrm{6}\mathit{\gamma }\left(c_{\mathrm{2}}\right(\mathit{\beta }/\mathrm{3}-\mathrm{1})+\mathrm{1})},\text{(23)}& {\displaystyle }{\mathit{\lambda }}_{\mathrm{3}}=c_{\mathrm{1}}^{\mathit{\beta }/\mathrm{3}-\mathrm{1}},\text{(24)}& {\displaystyle }{\mathit{\lambda }}_{\mathrm{4}}={\left({\displaystyle \frac{\mathrm{3}\mathit{\mu }}{\mathit{\rho }g{k}_{\mathrm{w}}}}\right)}^{\mathit{\beta }/\mathrm{3}}.\end{array}$$

From this analytical solution, we can see that slope increases with rock uplift rate and decreases with k_df, which would increase with bedrock erodibility and/or debris flow frequency. The sensitivity of slope to changes in U and k_df is strongest when $\mathit{\alpha }-\mathit{\beta }/\mathrm{3}$ is small and gets weaker as $\mathit{\alpha }-\mathit{\beta }/\mathrm{3}$ increases. In addition, the relation between S and A will depend on b, β, and γ. The sign of the exponent, N, controls whether slope decreases or increases with drainage area, A. Without downstream increases in debris flow volume (i.e., γ=0), channel widening leads to debris flow thinning and therefore to steepening of the steady-state channel slope with increasing drainage area. Slope decreases with drainage area when N<0 or equivalently when

$$\begin{array}{}\text{(25)}& \mathit{\beta }/\mathrm{3}(b-\mathit{\gamma }{c}_{\mathrm{2}})+\mathit{\gamma }(c_{\mathrm{2}}-\mathrm{1})<\mathrm{0}.\end{array}$$

Assuming c₂, the exponent in the power law relating peak debris flow discharge to debris flow volume is less than 1 (Rickenmann, 1999), then γ(c₂−1) will always be negative. Therefore, N will also always be negative when $b-\mathit{\gamma }{c}_{\mathrm{2}}<\mathrm{0}$, or γ>0.38 given values of b and c₂ used here, though this is a more restrictive condition than is necessary to ensure N<0. Plots of N as a function of α and β for different values of γ demonstrate that N is less than zero in cases where γ=0.25 and β<2 (Fig. C1). In any case, we see that high α and low β values promote limited variations in S with A by reducing the magnitude of N. However, a caveat is that the exact dependency of S on A may also be highly influenced by the relationship between A and k_df, which accounts for debris flow frequency (Stock and Dietrich, 2006). These dependencies will not be formally explored in this first study but remain to be explored through future field or modeling studies. Furthermore, the dependence of N on b motivates field observation of the debris flow channel width as a function of drainage area. In the following numerical experiments we aim to confirm the validity of this analytical solution and its consistency to a more process-based model. We further aim to determine the range of α and β values that produce channel profiles consistent with the observational constraint on the relationship between slope and drainage area as summarized by Eq. (2).

2.5 Numerical experiments

Our numerical experiments have two goals, which are treated in turn. First, we assessed which erosion laws, as defined by different values of α and β, can reproduce the first-order characteristics of observed channel longitudinal profiles, as well as how this may be affected by the choice of debris flow routing model (i.e., process-based or empirical). Second, we performed a series of simulations aimed at understanding the sensitivity of A_df and S_df to model parameters.

2.5.1 A family of debris flow incision laws

We explored model behavior for different values of α and β in the family of incision laws described by Eqs. (5) and (6) by comparing modeled, steady-state longitudinal profiles with those typical of debris-flow-dominated terrain. We did not try to recreate the channel morphology observed within specific watersheds or geographic regions. Instead, we aimed to provide some constraints on α and β by identifying ranges for these two exponents that resulted in longitudinal channel profiles that are consistent with observed changes in the relationship between slope and contributing area in natural channels traversed by debris flows (Stock and Dietrich, 2003).

A landscape evolution model designed to simulate the evolution of channel longitudinal profiles in response to both debris flow and fluvial erosion should produce steady-state channel profiles that are well described by Eqs. (1) and (2). Equation (1), which was formulated by Stock and Dietrich (2003) as part of an analysis of channel morphology across a range of geographic areas, suggests that channel slope increases or remains approximately constant as drainage area decreases. We performed an analysis of 31 channel longitudinal profiles in the San Gabriel Mountains, USA, to determine the frequency with which channel slope decreased as drainage area decreased. The San Gabriel Mountains were chosen for this analysis because some of our model parameter choices are based on previous studies in this mountain range and topography is in an approximate steady state (DiBiase et al., 2012). We extracted channel profiles for a subset of catchments with ¹⁰Be catchment-averaged erosion rates (DiBiase et al., 2010), where we eliminated catchments with signs of disequilibrium such as knickpoints. In 30 of the 31 catchments, in which erosion rate varied widely from <0.1 mm yr⁻¹ to more than 1 mm yr⁻¹, slope increased or remained approximately constant as drainage area decreased. In one catchment, there was a difference of 0.03 between the maximum slope along the channel profile and the top of channel profile (Fig. B2).

We therefore assessed model performance for different α and β in two ways. First, we computed the R² associated with the best fit to Eq. (2). We allowed S_df, A_df, and a₂ to vary freely when fitting to Eq. (2). Second, we examined the difference between the maximum slope along the channel profile, S_max, and the slope at the channel head, S_ch. The second criteria focused on checking a basic morphologic property observed in natural channels, namely that channel slope generally increases or remains constant as drainage area decreases in quasi-steady-state landscapes.

We assessed performance of the landscape evolution model with different values of α and β when using the process-based routing model and when using the empirical routing model. Using the process-based model, we performed a numerical experiment where we varied α; β; pore pressure diffusivity, D; viscosity of the pore fluid, η; the friction parameter, μ₂; the debris flow erodibility coefficient, k_df; and instantaneous fluvial erodibility coefficient, k_e (Appendix A) within the ranges specified in Table G2. We allowed some variation in model parameters other than α and β to ensure trends between α, β, and model performance metrics were not specific to a particular subset of the parameter space. We selected 500 parameter sets using a Latin hypercube sampling strategy. We performed an analogous numerical experiment using the empirical routing model where we sampled 4000 different parameter sets with varying values of α, β, instantaneous fluvial erodibility (k_e), debris flow erodibility (k_df), viscosity (μ), yield strength (τ_y), and debris flow volume parameters M₀ and γ within prescribed ranges (Table G3). We were able to perform a greater number of simulations using the empirical model because it is less computationally demanding.

All simulations began with an initial condition determined by the analytical solution for a steady-state fluvial channel, specifically

$$\begin{array}{}\text{(26)}& S=(U/K{)}^{\mathrm{1}/n_{\mathrm{s}}}{A}^{-m_{\mathrm{s}}/n_{\mathrm{s}}}.\end{array}$$

Simulations ended once an approximate steady state had been reached, which typically took 10⁶–10⁷ years.

2.5.2 Sensitivity analysis

We performed sensitivity analyses using both the process-based and empirical routing models to explore how the topographic signature of debris flow incision is likely to be expressed in different settings. Motivated by the results of our numerical experiments to constrain α and β and by insights from the simplified analytical solution, we set α=6 and β=1 for the sensitivity analysis. The analytical solution for slope suggests that the exponent, N, that controls whether slope increases or decreases with drainage area will be negative for γ>0.25 and relatively low in absolute value for γ=0 when α=6 and β=1. Therefore, we focused on this combination of α and β as it is likely to yield results that are consistent with observations as summarized by Eq. (2). We focused, in particular, on understanding relationships between model parameters and resulting longitudinal profile form as quantified by A_df and S_df in steady-state longitudinal profiles.

To perform the sensitivity analysis with the process-based routing model, we used a Latin hypercube sampling strategy to select 1500 sets of parameters where instantaneous fluvial erodibility, k_e; debris flow erodibility, k_df; viscosity of the pore fluid, η; pore fluid diffusivity, D; the friction factor, μ₂; and rock uplift rate, U varied within the ranges defined in Table G4. The sensitivity analysis using the empirical routing approach was analogous, but we were able to perform a greater number of simulations. We used a Latin hypercube sampling strategy to select 4000 sets of parameters where instantaneous fluvial erodibility, k_e; debris flow erodibility, k_df; viscosity, μ; yield strength, τ_y; rock uplift rate, U; and debris flow volume parameters M₀ and γ varied within the ranges defined in Table G5.

We performed a qualitative sensitivity analysis by visually examining model output using colored scatter plots and also performed a quantitative global sensitivity analysis using the PAWN method (Pianosi and Wagener, 2015) as implemented with the SALib Python package (Herman and Usher, 2017). The PAWN method is a density-based global sensitivity analysis method. The output of a PAWN sensitivity analysis consists of a sensitivity index for each input variable that summarizes its relative contribution to uncertainty in the output. The sensitivity index varies from 0 to 1, with greater values indicating a greater relative importance of the parameter. By comparing the magnitudes of the sensitivity indices for different input parameters, we were able to rank them in terms of relative importance. We separately assessed sensitivity to each of two model outputs, A_df and S_df, since these two metrics summarize basic morphologic information about steady-state channel profiles. We performed PAWN sensitivity analyses separately for models that employ the process-based and empirical routing approaches. The PAWN sensitivity analysis allowed us to rank input variables in terms of their relative importance for determining A_df and S_df.

3 Results

3.1 Constraints on a debris flow incision law

3.1.1 Process-based routing model

At large drainage areas, modeled profiles exhibit a power law scaling between slope and drainage area that is expected based on the fluvial incision law (Figs. 2, 3). The R² value associated with a fit to Eq. (2) was greater than 0.95 for approximately 93 % of the modeled profiles. At low-drainage areas, however, all parameter combinations produced channel profiles where slope began to decrease as drainage area decreased. In other words, the difference between the maximum slope, S_max, along the channel profile and the slope at the channel head, S_ch, was positive and regularly exceeded 0.2 in cases where β>2 (Fig. 2). This decrease in slope at low-drainage areas, particularly a decrease in magnitude of more than 0.05, is inconsistent with Eq. (2) and observations (Fig. B2) that indicate slope continues to increase or remain approximately constant as drainage area decreases. Differences between S_max and S_ch decreased rapidly as $\mathit{\alpha }/\mathit{\beta }$ increased (Fig. 2). An increase in slope with drainage area near the channel head, however, is not an inevitable consequence of using the process-based routing model (Appendix D).

Figure 2 Numerical experiments using the process-based routing model to determine which slope (α) and depth (β) exponents in the debris flow erosion law create profiles consistent with those seen in nature. (a) The difference between the maximum slope (S_max) and the slope at the channel head (S_ch), which is the first point plotted on the profiles, generally decreases with α. Color indicates the value of β, highlighting an increase in S_max−S_ch when β increases, and hence poor model performance in many cases with β>1. (b) The difference between the maximum slope and the slope at the channel head generally increases with β. Color indicates the value of α, with greater values of α generally leading to smaller S_max−S_ch and better model performance. (c) The difference between the maximum slope and the slope at the channel head decreases rapidly with $\mathit{\alpha }/\mathit{\beta }$. (d–i) Representative profiles for different α and β. Substantial reductions in channel slope at small drainage areas are inconsistent with observations (e.g., Eq. 2).

Figure 3 Numerical experiments using the empirical routing model to determine which slope (α) and depth (β) exponents in the debris flow erosion law create profiles consistent with those seen in nature. (a) The difference between the maximum slope (S_max) and the slope at the channel head (S_ch) generally decreases with α. Color indicates the value of β, highlighting poor model performance, as measured by S_max−S_ch, in many cases when β increases. (b) The difference between the maximum slope and the slope at the channel head generally increases with β. Color indicates the value of α, with greater values of α generally leading to better model performance. (c) The difference between the maximum slope and the slope at the channel head decreases rapidly with $\mathit{\alpha }/\mathit{\beta }$. (d–i) Representative profiles for different α and β for cases where γ=0. Substantial reductions in channel slope at small drainage areas are inconsistent with observations (e.g., Eq. 2).

3.1.2 Empirical routing model

Modeled profiles exhibited the expected power law scaling between slope and drainage area at high-drainage areas where fluvial incision dominated debris flow incision. The coefficient of determination (R²) value associated with a fit to Eq. (2) was greater than 0.95 for all modeled profiles. As with results obtained using the process-based routing model, some parameter combinations produced channel profiles where the maximum channel slope was not observed at the channel head (Fig. 3). Differences between S_max and S_ch are minor when α=6 and β=1 but become more substantial as α decreases and/or as β increases (Fig. 3).

Figure 4 Numerical experiments using the empirical routing model that highlight the importance of the volume–area scaling exponent γ. The difference between the maximum slope (S_max) and the slope at the channel head (S_ch) generally decreases with α, the slope exponent in the debris flow erosion law as seen when comparing across panels (a)–(f), and increases with β, the depth exponent in the debris flow erosion law (as seen in marker color). Model performance, as measured by S_max−S_ch, is most sensitive to changes in α and β when γ is less than approximately 0.3.

More generally, the extent to which modeled channel profiles exhibit a decrease in slope at small drainage areas depends on α, β, and the exponent γ that controls the relationship between debris flow volume and drainage area (Fig. 4). In cases where γ<0.25, numerous combinations of α and β lead to decreases in slope at low-drainage areas. For any choice of $\mathrm{2}\le \mathit{\alpha }\le \mathrm{8}$ and β≈1.5 or less, S_max−S_ch was always less than 0.05 (Fig. 4). For cases where β>2, α needed to be approximately 5 or greater to maintain $S_{\max }-S_{\mathrm{ch}}<\mathrm{0.05}$ for all values of γ. Differences between S_max and S_ch increase as β increases and/or as α decreases. We were unable to directly explore the effects of spatial variations in debris flow volume using the process-based routing model, where we neglect changes in debris flow volume along the flow path (i.e., γ=0).

3.2 Steady-state forms of channel profiles

3.2.1 Process-based routing model

Two defining characteristics of the simulated steady-state channel profiles, the near-constant slope that they approach near the channel head and the minimum drainage area at which there is a power law scaling between slope and drainage area, can be summarized by the following two metrics: S_df and A_df. Results of the sensitivity analysis demonstrate that neither A_df nor S_df are particularly sensitive to parameters that primarily affect flow mobility, including viscosity of the pore fluid (η), friction parameters (μ₂), and pore fluid pressure diffusivity (D) (Figs. 5, 6, Table 2). Rather, A_df is most sensitive to the instantaneous fluvial erodibility (k_e), debris flow erodibility (k_df), and rock uplift rate (U), whereas S_df is controlled predominantly by k_df and U (Table 2).

Table 1 Model parameters.

Table 2 PAWN sensitivity indices of the process-based model.

Figure 5 Scatter plots summarizing results of the sensitivity analysis with the process based debris flow routing model. Sensitivity of A_df (left column) and S_df (right column) to particular model parameters is indicated when there is a gradient in color, whereas plots with no spatial pattern in color indicate a lack of sensitivity. The relationship between A_df, k_e, and different parameters related to debris flow processes illustrate sensitivity to the debris flow erodibility coefficient, k_df, and the instantaneous fluvial erodibility coefficient (k_e). The morphologic parameter S_df (right column) is most sensitive to k_df. Parameter definitions and units can be found in Table G2.

Figure 6 Steady-state longitudinal profiles produced by the process-based model for various parameter combinations. Variations in (a) the debris flow erodibility coefficient, k_df; (b) the instantaneous fluvial erodibility, k_e; and (c) rock uplift rate, U, drive greater changes in the channel morphology, as summarized by the relationship between channel slope and drainage area, relative to variations in (d) the viscosity of the pore fluid, η; (e) friction factor, μ₂; and (f) pore pressure diffusivity, D. Default parameter values are as follows: k_df=0.0001 m^1−β s⁻²; $k_{\mathrm{e}}=\mathrm{5}\times {\mathrm{10}}^{-\mathrm{14}}$ ${\mathrm{m}}^{\mathrm{5}/\mathrm{2}}$ s² ${\mathrm{kg}}^{-\mathrm{3}/\mathrm{2}}$; U=0.5 mm yr⁻¹; η=60 Pa s; μ₂=0.62; $D=\mathrm{5}\times {\mathrm{10}}^{-\mathrm{6}}$ m² s⁻¹.

Figure 7 Results from the 1500 process-based model simulations in the sensitivity analysis with α=6 and β=1 that show a power law relationship between (a) rock uplift rate (U) and A_df and (b) U and S_df. The dashed line indicates a best-fit curve. Patterns in symbol color, which represents the magnitude of the debris flow erodibility coefficient, k_df, can be used to infer that systematic variations in k_df with rock uplift rate, U, would influence the relationship between U and A_df and U and S_df.

Figure 8 Scatter plots summarizing results of the sensitivity analysis with the empirical debris flow routing model. The relationship between A_df, k_e, and different parameters related to debris flow erosion (left panels) illustrate sensitivity to the volume–area scaling exponent (γ), debris flow erodibility coefficient, (k_df), and the instantaneous fluvial erodibility coefficient (k_e). The morphologic parameter S_df (right panels) is most sensitive to the volume exponent (γ) and k_df. Parameter definitions and units can be found in Table 1.

Figure 9 Steady-state longitudinal profiles produced by the empirical routing model for various parameter combinations. Variations in (a) the debris flow erodibility coefficient, k_df; (b) the instantaneous fluvial erodibility, k_e; (c) rock uplift rate, U; and (d) the volume–area scaling exponent, γ, drive greater changes in the channel morphology, as summarized by the relationship between channel slope and drainage area, relative to variations in (e) the volume–area scaling coefficient, M₀, and (f) viscosity, η. Default parameter values are as follows: k_df=0.0001 m^1−β s⁻²; $k_{\mathrm{e}}=\mathrm{5}\times {\mathrm{10}}^{-\mathrm{14}}$ ${\mathrm{m}}^{\mathrm{5}/\mathrm{2}}$ s² ${\mathrm{kg}}^{-\mathrm{3}/\mathrm{2}}$; U=0.5 mm yr⁻¹; γ=0.5; M₀=3000 m^3−2γ; μ=500 Pa s.

The sensitivity of steady-state long-channel profiles to changes in rock uplift rate leads to power law relationships between A_df and U and between S_df and U (Fig. 7). By randomly sampling model parameters, including rock uplift rate, within prescribed ranges, we assume that none are correlated with each other. However, this is unlikely in natural landscapes, and correlations are expected. For example, we may expect that debris flow frequency, F_df, increases with rock uplift rate. The consequences of such a correlation can be seen by examining the effects of k_df on A_df and S_df for a given rock uplift rate (Fig. 7). Increases in k_df, for a given rock uplift rate, lead to increases in A_df and decreases in S_df. A correlation between rock uplift rate and k_df would therefore be likely to influence the fit between S_df and A_df.

Figure 10 (a, b) Results from the 4000 empirical routing model simulations in the sensitivity analysis with α=6 and β=1 show a power law relationship between (a, c) uplift (U) and A_df and (b, d) between U and S_df. Dashed lines show best-fit curves. For a given rock uplift rate, changes in symbol color also highlight how A_df and S_df vary with k_df such that a systematic variation in k_df with U would influence the power law fit.

3.2.2 Empirical routing model

Simulations indicate that S_df is most sensitive to changes in γ, which controls the relationship between debris flow volume and drainage area, and k_df, which is related to debris flow frequency and bedrock erodibility followed by rock uplift rate (Figs. 8, 9, Table 3). Typical values of S_df decrease with k_df and increase with γ and U but are not strongly controlled by k_e, τ_y, μ, and M₀ (Table 3). The area at which there is a transition to fluvial power law scaling between slope and area, A_df, is most sensitive to γ, k_df, k_e, U, and k_df, whereas it is relatively insensitive to M₀, μ, and τ_y (Figs. 8, 9, Table 3). Mean values of A_df tend to decrease strongly with k_e, increase with k_df, and decrease with γ. Parameters more directly related to the physical properties of the debris flows; viscosity, μ; and yield strength, τ_y, had relatively minor control over S_df and A_df (Figs. 8, 9).

Table 3 PAWN sensitivity indices of the empirical model.

There is a power law relationship between A_df and rock uplift rate, U, although there is considerable scatter due to the wide range of parameter values included in the sensitivity analysis (Fig. 10). By randomly sampling the model parameters within prescribed ranges, we assume that none are correlated with each other. However, this is unlikely in natural landscapes, and correlations are expected. For example, we may expect that F_df and/or γ increase with rock uplift rate. Again, we explore the consequence of such correlations by examining patterns in colored scatter plots that can help visualize the impact of k_df on the relationship between U and either A_df or S_df (Fig. 10). A considerable amount of scatter in the relationship between U and A_df appears to be attributable to variations in γ and k_df, as expected from results of the PAWN sensitivity analysis (Table 3).

4 Discussion

4.1 Constraints on a geomorphic transport law for debris flow incision

Results indicate that many members within the proposed family of debris flow incision laws, as formulated by Eqs. (5) and (6), produce channel profiles that are consistent with observations from natural landscapes (Figs. 2, 3). This is true within a wide range of the parameter space explored here, including for a range of γ that covers the variability observed across several different geographic regions reported by Santi and Morandi (2013). Data and numerical experiments presented here are not capable of differentiating among these potential debris flow incision laws, although there are cases where α<3 and/or β>2 generally performed poorly (Figs. 2, 3). Additional work is needed to formulate and test a debris flow incision law, including incision laws not restricted to the form of Eq. (5) (e.g., Stock and Dietrich, 2006). Stock and Dietrich (2006), for example, present a debris flow incision law based on inertial stress. Analogously, there are a range of exponents used in the generalized stream power incision law for fluvial erosion and work continues in an effort to constrain those exponents (e.g., Clubb et al., 2016; Turowski, 2018, 2021). Based on the extent to which key characteristics of long-channel profile morphology is affected by k_df, k_e, γ, and U, ideal landscapes for testing a debris flow incision law would be ones where there are constraints on these parameters (Figs. 5, 7, 8, 10).

Here, we assess different debris flow incision laws based on their ability to reproduce a general pattern in slope–area data (i.e., Eq. 2) observed in debris-flow-dominated landscapes (e.g., Fig. B2). To produce the observed steady-state morphology of debris-flow-dominated long profiles with a slope that is approximately constant or slowly decreasing with A, examination of the analytical solution for slope in the upper channel network given by Eq. (19) demonstrates that $\mathit{\gamma }>b/c_{\mathrm{2}}$ is a sufficient, though more restrictive than necessary, condition to ensure that steady-state slopes decrease as drainage area increases. For simulations presented here, c₂=0.78 and b=0.3, which implies γ>0.38. Examination of Eq. (20) shows that when γ=0.25, N is generally small in magnitude, particularly for greater values of α (Fig. C1). Results of numerical simulations are consistent with the analytical solution, with γ>0.25 being sufficient to ensure that slope is approximately constant or decreases as drainage area increases (Fig. 4). Equation (19) also demonstrates how the magnitude of the exponent, N, is modulated by α and β (Fig. C1). The analytical solution is consistent with numerical simulations using both the empirical and process-based routing models that show a general trend toward greater differences between maximum channel slope and slope at the channel head as β increases and as α decreases (Figs. 2, 3). Since both numerical and analytical solutions to the model equations demonstrate how b and γ exert a strong control on determining the basic morphology of channel profiles in debris-flow-dominated reaches near the channel head (i.e., the sign of N in Eq. 19), additional model evaluation criteria beyond those proposed here would be needed to test a debris flow incision law. One possibility would be to assess the ability of a debris flow incision law to reproduce observed trends between erosion rate and A_df or a₁, such as that observed by Penserini et al. (2017) in the Oregon Coast Range.

In general, the proposed empirical and process-based approaches for estimating bulk debris flow properties along the channel profile do not appear to result in different model behavior (Figs. 2, 3, 6, 9). It is not possible to directly compare the longitudinal profiles produced by the two different routing models with the same values of α and β because the parameters that determine debris flow mobility are different among the two models. For example, the process-based model has no yield strength parameter, yet this parameter plays a key role in determining debris flow runout in the empirical model. Although the flow depth, velocity, and passage time predicted by the two routing models will undoubtedly vary, these variations are not sufficient to alter the extent to which different values of α and β produce modeled profiles that are consistent or inconsistent with observations of channel morphology in debris-flow-dominated terrain (Figs. 2, 3). Furthermore, results using the process-based and empirical routing approaches both highlight the sensitivity of A_df and S_df to rock uplift rate and k_e, k_df, and γ relative to model parameters related to flow mobility (Figs. 6, 9). These similarities in model behavior are encouraging because the empirical routing approach provides a framework to estimate bulk debris flow properties using quantities that can be computed from a digital elevation model, specifically upstream contributing area and slope, that could be used in future efforts to more efficiently explore alternative debris flow incision laws.

4.2 Steady-state forms of longitudinal channel profiles

Model results help clarify the roles played by debris flow and fluvial erosion processes in setting longitudinal profile form in the upper channel network. Changes in parameters related solely to fluvial erosion do not have a strong influence on S_df, which simulations demonstrate is primarily controlled by changes in debris flow processes (Figs. 5, 8, 9). Specifically, increases in rock uplift rate, U, and γ or decreases in k_df promote increases in S_df, assuming all else is fixed. In contrast, A_df is controlled by a combination of fluvial and debris flow processes (Figs. 5, 6, 8, 9). On average, increases in the instantaneous fluvial erodibility lead to decreases in A_df, whereas increases in the debris flow erodibility coefficient promote increases in A_df. If S_df was a constant set, for example, by soil geotechnical properties related to slope stability, then A_df could be estimated by the area at which the steady-state fluvial channel gradient reaches S_df. However, this type of threshold behavior of S_df is not what we observe. Rather, simulations demonstrate that S_df varies with U, γ, and k_df and independent changes to either debris flow incision processes or fluvial processes are sufficient to influence A_df. Thus, accounting for both debris flow and fluvial processes is important to understand how the steep channel network will respond to changes in tectonic or climatic forcing.

Simulations, assuming α=6 and β=1, indicate that debris flows may frequently traverse channel reaches at larger drainage areas without influencing longitudinal channel profile form in a substantial way. Debris flows routed with the empirical model traversed the entire model domain (8 km²) in greater than 99 % of simulations, but the median A_df was approximately 0.4 km². A primary reason for this is the sensitivity of the debris flow incision law to slope when α=6, which means that portions of the channel profile could regularly be traversed by debris flows but they would do little work to erode bedrock when the channel slope is modest. This result indicates that the presence of debris flow activity in natural channels, as indicated by debris flow deposits, may not be a reliable indicator of the importance of debris flow incision. However, additional work is needed to constrain the relationship between slope and debris flow incision rates as well as to explore the influence of debris flows mobilizing large-caliber sediment below A_df that would otherwise shield the bed from subsequent fluvial erosion. An additional consequence of the nonlinear (α>>1) relationship between slope and debris flow incision and the observation that the majority of debris flows remain mobile throughout the entire model domain is that A_df and S_df are less sensitive to parameters related to debris flow properties, namely viscosity and yield strength in the case of the empirical routing model (Figs. 9, 10) and viscosity of the pore fluid, friction parameters, and pore pressure diffusivity in the case of the process-based routing model (Figs. 5, 6). As long as debris flows are sufficiently mobile to traverse moderate slopes, these parameters primarily affect the debris flow incision rate by changing debris flow depth and/or passage time. Because both debris flow depth and passage time are linearly related to the debris flow incision rate when β=1, a factor of 2 increase or decrease in flow depth or passage time would only require a relatively small adjustment in slope to compensate for the ability of debris flows to balance the imposed rock uplift rate at steady state. PAWN sensitivity analyses conducted using results of 4000 simulations for cases where α=6 and β=2 (Table F1) and 4000 simulations where α=4 and β=1 (Table F2) demonstrate this conclusion, namely A_df and S_df being less sensitive to parameters related to debris flow properties, also holds for other combinations of α and β.

4.3 Tectonics from debris flow processes and topography

Model results support previous observations indicating that the morphology of channel profiles in debris-flow-dominated landscapes may provide constraints on erosion rates in steady-state landscapes (Figs. 7, 10). Penserini et al. (2017) document an inverse relationship between a₁ and E in the central Oregon Coast Range based on analyses of channel profiles in six watersheds with catchment-averaged erosion rates determined from cosmogenic radionuclide analysis. Cast in terms of the morphologic variables used here to describe channel profiles, specifically using A_df in place of a₁, the results from Penserini et al. (2017) indicate that A_df increases with E. Simulations confirm this pattern of increasing A_df with E but also highlight the importance of constraining relationships between E and k_df, as well as E and γ, in order for A_df to serve as a proxy for erosion rates in an absolute sense (Figs. 7, 10). This indicates a need to prioritize constraining relationships between E and debris flow frequency and bed erodibility (which together control k_df) and E and the rate at which debris flow volume increases with drainage area (e.g., γ) rather than potential relationships between E and debris flow mobility parameters, such as viscosity and yield strength. Although not explored here, differences in m_s and n_s, which vary among landscapes, are also likely to influence relationships between E and A_df. Similarly, several landscape parameters assumed fixed in this study, such as the area at the channel head, A₀, or the channel width scaling, b, may be influenced by debris flow erosion and exert control over the long-term channel profile morphology. Analytical (i.e., Eq. 25) and numerical (Figs. 5, 8) model results indicate that relationships between drainage area and debris flow volume, as well as drainage area and channel width, in particular, play key roles in determining the morphology of the debris-flow-dominated channels.

Simulations indicate increases in erosion rate, or equivalently rock uplift rate in a steady-state landscape, lead to an increase in S_df (Figs. 7, 10). Interestingly, Penserini et al. (2017) found no systematic variation in S_df with erosion rate in the Oregon Coast Range. However, the lack of any relationship between S_df and erosion rate (E) in the Oregon Coast Range may result from a correlation between E and debris flow frequency and/or E and debris flow volume. Simulations of steady-state channel profiles indicate that S_df increases with rock uplift rate in cases where there is no relationship between rock uplift rate (U) and k_df but imposing an increase in k_df with U would alter that relationship (Figs. 7, 10). Thus, two avenues are essential for improving our ability to use the morphology of debris-flow-dominated channels as a proxy for erosion rate, namely (1) extracting upland channel morphology in a larger number of catchments with constrained erosion rates and (2) gathering evidence on the interconnections of key parameters (Stock and Dietrich, 2006).

4.4 Model applications and limitations

We describe a model that provides a framework for exploring the effect of episodic debris flows on channel longitudinal profiles. The model also comes with several limitations. We assume that all debris flows initiate at the channel head, whereas debris flow initiation locations in natural landscapes will be more varied. Past work highlights the role that network structure, specifically the number of debris flow initiation locations upstream of a given channel reach, may play in controlling channel form (Stock and Dietrich, 2006). Variations in the spatial distribution of debris flow initiation locations within a watershed could be explored within the model presented here by prescribing debris flow frequency as a function of drainage area (e.g., Fig. D1) or explicitly modeling multiple debris flows from different initiation locations. Additionally, the scaling between channel width and drainage area may differ in the upper portion of the channel network from previously reported relationships that are derived from data at larger drainage areas where fluvial processes are dominant (DiBiase and Whipple, 2011). However, as additional data become available to better parameterize channel morphometry in steep landscapes (Neely and DiBiase, 2023), particularly at small drainage areas, new parameterizations can be exchanged with those presented here. Similarly, when using the empirical debris flow routing model, we rely on relationships between debris flow volume and drainage area that were derived using data collected primarily at drainage areas greater than 0.1 km² (Santi and Morandi, 2013). Lastly, the process-based debris flow routing model presented here assumes that debris flow volume is constant and does not change along the travel path, although quantifying controls on sediment entrainment by debris flows and incorporating entertainment into process-based debris flow routing models are areas of active research (Iverson, 2012; Iverson and Ouyang, 2015; McCoy et al., 2012; Haas and Woerkom, 2016). Advances in our understanding of how debris flows entrain sediment could allow for more detailed comparisons between empirical and process-based approaches to sediment bulking in the proposed landform evolution model.

The landform evolution model presented here may serve as a basis for future studies that aim to test or validate potential debris flow incision laws, incorporate debris flow incision into 2D landscape evolution models, or explore how the upper channel network responds to tectonic or climatic perturbations. The process-based routing model may be best suited for modeling 1D channel profiles where changes in flow volume can be neglected and debris flow constituents are sufficiently well known to allow for estimates of the model parameters, thereby minimizing the number of numerical experiments needed to characterize model behavior. The empirical debris flow routing algorithm provides an efficient framework for investigating the effects of different debris flow bulking relationships and exploring large parameter spaces. It is also particularly promising for application in 2D landscape evolution models given its simplicity relative to process-based debris flow routing models and its ability to connect slope and drainage area, which are readily available in nearly all landscape evolution models, with bulk debris flow properties relevant to debris flow incision.

5 Conclusions

We present a novel framework for incorporating erosion by debris flows into a model for channel profile evolution. We propose two methods to estimate debris flow runout and bulk debris flow properties (e.g., depth, velocity) throughout the channel network, one based on a process-based debris flow routing model and the other based on an empirical routing approach. Combined with a geomorphic transport law describing the relationship between debris flow depth, channel slope, and a debris flow incision rate, we are able to quantify spatial variations in the debris flow incision rate throughout the channel network. We explore the performance of a family of potential debris flow incision laws by comparing the form of modeled longitudinal channel profiles with those typically observed in debris-flow-dominated landscapes. Results demonstrate that a debris flow incision law based on flow depth, slope, and debris flow passage time can reproduce the relationship between slope and drainage area that has been interpreted as a topographic signature of debris flows, given general constraints on empirical exponents that relate flow depth and local channel slope to the incision rate. Since a large subset of the proposed family of erosion laws is capable of reproducing this topographic signature of debris flows, additional criteria and more precise bounds on poorly constrained model parameters are needed to test and validate debris flow erosion laws. Simulations indicate that both A_df and S_df have potential to serve as a morphologic proxy for the catchment-averaged erosion rate. However, both A_df and S_df are sensitive to debris flow frequency and debris flow erodibility (k_df), and the empirical exponent characterizing how debris flow volume increases with drainage area (γ), indicating that the utility of such a proxy would depend on the extent to which relationships between erosion rate, k_df, and γ could be constrained. Results provide a general framework that can be used to test debris flow incision laws and explore the relative importance of debris flow versus fluvial processes in shaping channel profiles in steep landscapes. Results take initial steps toward the broader inclusion of bedrock erosion by debris flows into landscape evolution models and provide insights into the relationship between debris flow processes and channel profile morphology in steep landscapes.

Appendix A: Stochastic stream power incision model

The parameterization for the stochastic stream power model is not tuned to any particular landscape or geographic region, but relies on values and relationships that are based on typical values reported by Lague (2014) for high discharge settings and by DiBiase and Whipple (2011) for the San Gabriel Mountains. Given the parameters listed in Table G1, we follow Lague (2014) and compute the critical shear stress for bed load transport, τ_c according to ${\mathit{\tau }}_{\mathrm{c}}=\mathrm{0.045}g({\mathit{\rho }}_{\mathrm{s}}-{\mathit{\rho }}_{\mathrm{w}})D_{\mathrm{eff}}$, with D_eff=0.09 m the effective grain size. The stochastic-threshold prediction for the slope exponent in the fluvial incision law, n_s, is given by $n_{\mathrm{s}}={\mathit{\beta }}_{\mathrm{f}}/{\mathit{\alpha }}_{\mathrm{f}}(k+\mathrm{1})/(\mathrm{1}-{\mathit{\omega }}_{\mathrm{s}})$, where β_f=0.7 is the slope exponent in the hydraulic friction law, α_f=0.6 is the discharge exponent in the hydraulic friction law, ω_s is the at-a-station width scaling exponent, and k=0.5 is the discharge variability coefficient. The prediction for the area exponent is given by $m_{\mathrm{s}}=(c-b)(k+\mathrm{1})/(\mathrm{1}-{\mathit{\omega }}_{\mathrm{s}})$, where c=1 is the mean discharge–area scaling exponent and b=0.3 is the width–area scaling coefficient (Lague, 2014). Lastly, the fluvial erodibility coefficient is calculated as (Lague, 2014)

$$\begin{array}{}\text{(A1)}& K=k_s{k}_{wq}^{-\mathit{\alpha }{n}_{\mathrm{s}}/\mathit{\beta }}{R}_c^{m_{\mathrm{s}}/c},\end{array}$$

where R_c=0.28 m denotes the mean annual runoff, $k_{wq}=R_c^{b/c}/k_{\mathrm{w}}$ denotes a width factor, and

$$\begin{array}{}\text{(A2)}& k_s=\left({\displaystyle \frac{a{\mathit{\alpha }}_{\mathrm{f}}(\mathrm{1}-{\mathit{\omega }}_{\mathrm{s}})k_t^{n_{\mathrm{s}}/{\mathit{\beta }}_{\mathrm{f}}}{k}^{k+\mathrm{1}}\mathrm{\Gamma }(k+\mathrm{1}{)}^{-\mathrm{1}}}{(k+\mathrm{1})(k+\mathrm{1}-a{\mathit{\alpha }}_{\mathrm{f}}(\mathrm{1}-{\mathit{\omega }}_{\mathrm{s}}\left)\right)}}\right){\mathit{\tau }}_{\mathrm{c}}^{a-n_{\mathrm{s}}/{\mathit{\beta }}_{\mathrm{f}}}{k}_{\mathrm{e}}.\end{array}$$

Here, $k_t=g{\mathit{\rho }}_{\mathrm{w}}{N}^{\mathrm{3}/\mathrm{5}}$, N=0.05 denotes the Manning friction coefficient, and a=1.5 is a shear stress exponent. The rate of fluvial erosion is then computed according to $E_{\mathrm{f}}=K{A}^{m_{\mathrm{s}}}{S}^{n_{\mathrm{s}}}$.

Appendix B: San Gabriel Mountains channel morphology

Scaling relationships that relate drainage area and channel width are often derived from data that include drainage areas greater than those modeled in this study. Recall that in this study we use the term channel to broadly refer to a concentrated axis of erosion along valley bottoms. Since debris flows initiate and traverse channels at low-drainage areas, we quantified channel width at drainage areas less than 3 km² in the San Gabriel Mountains. More specifically, we focused on a region in the San Gabriel Mountains burned by the 2016 Fish Fire. A series of rainstorms in the first year following the fire led to runoff-generated debris flows that scoured many low-order channels to bedrock (Rengers et al., 2021; Tang et al., 2019). A post-event DEM derived from airborne lidar provided an opportunity to quantify channel width in this location. We estimated channel width by visually examining cross-channel profiles and identifying channel banks or distinct breaks in slope that indicated a shift from the hillslope to channel. Data indicate that channel width increases as a power law function of drainage area with an exponent of approximately 0.28 (Fig B1). These data provide support for the width–area scaling used here, but we acknowledge that better characterization of the morphometry of debris-flow-dominated channels would improve the utility of landform evolution models that represent steep, low-order channels.

Figure B1 Estimates of channel width, estimated from a high-resolution lidar-derived digital elevation model, as a function of drainage area for a portion of the San Gabriel Mountains, USA. The area burned in the 2016 Fish Fire and experienced a series of debris flows during the first rainy season following the fire that scoured valleys and channels to bedrock in many places. For comparison, relationships between channel width and drainage area as determined by DiBiase et al. (2012) are also shown.

Figure B2 We fit a line, S_fit, to binned slope–area data for areas between A₀ and A_df. (a) Based on analysis of 31 channel profiles in the San Gabriel Mountains, differences between S_fit(A_df) and S_fit(A₀) range from approximately 0.05 to 0.2. (b) Example of binned slope–area data for a channel profile where slope increases with decreasing drainage area from A_df to A₀. (c) Example of binned slope area data where slope decreases slightly as drainage area decreases from A_df to A₀.

Appendix C: Analytical solution

Figure C1 The analytical solution for steady-state channel slope indicates that slope is a power law function of drainage area with an exponent, N, given by Eq. (20). Slope increases with drainage area when N>0 and decreases with drainage area when N<0. The above plots show how the magnitude and sign of N vary with the two exponents in the debris flow incision law, α and β, for increasing values of γ, the debris flow volume–area scaling exponent, from (a) γ=0, (b) γ=0.25, γ=0.5, and γ=0.75.

Appendix D: Spatially variable debris flow frequency

The process-based routing model does not directly account for downstream changes in debris flow volume. When using the empirical routing model, for example, we prescribe debris flow volume as a function of drainage area according to $M=M_{\mathrm{0}}({\mathrm{10}}^{-\mathrm{6}}\cdot A{)}^{\mathit{\gamma }}$. However, we can scale debris flow frequency with drainage area in a way parameterizes an overall increase in the debris flow volume as drainage area increases. For a basic illustration of this parameterization and its effects on the model solution, we performed a series of simulations where we scale k_df by a factor of $\mathrm{1000}({\mathrm{10}}^{-\mathrm{6}}\cdot A{)}^{\mathit{\gamma }}/M_{\mathrm{0}}$. This parameterization leads to an increase in the total sediment transported by debris flows as a function of drainage area that is consistent with the way in which debris flow volume increases downstream when using the empirical routing model. Results are summarized in Fig. D1 and lead to patterns that are qualitatively consistent with those obtained when parameterizing downstream increases in debris flow volume with the empirical routing model (Fig. 9d).

Figure D1 Steady-state channel profiles using the process-based routing model where we parameterize an increase in debris flow frequency with drainage area, A, by scaling k_df by $\mathrm{1000}({\mathrm{10}}^{-\mathrm{6}}\cdot A{)}^{\mathit{\gamma }}/M_{\mathrm{0}}$. This parameterization leads to an increase in the total sediment transported by debris flows as a function of drainage area that is consistent with the way in which debris flow volume increases downstream when using the empirical routing model. All other parameters are fixed: $k_{\mathrm{df}}=\mathrm{5}\times {\mathrm{10}}^{-\mathrm{4}}$; $k_{\mathrm{e}}=\mathrm{5}\times {\mathrm{10}}^{-\mathrm{14}}$; η=60; μ₂=0.67; $\mathit{\mu }=\mathrm{5}\times {\mathrm{10}}^{-\mathrm{6}}$.

Appendix E: Process-based debris flow routing model

Figure E1 (a) A debris flow being routed down a channel profile using the process-based flow routing model. Flow depth is multiplied by a factor of 50 for display purposes. (b) The debris flow erosion rate varies spatially due to differences in slope and debris flow depth, including differences in temporal changes in flow depth as it passes over each point. The fluvial erosion rate is computed based on slope and drainage area at each point.

Appendix F: PAWN sensitivity analysis

Here, we report results for the PAWN sensitivity analysis when using erosion laws with α=6 and β=2 (Table F1) and as α=4 and β=1 (Table F2).

Table F1 PAWN sensitivity indices for the empirical model: α=6, β=2.

Table F2 PAWN sensitivity indices for the empirical model: α=4, β=1.

Appendix G: Model Parameters

Tables below provide details on the value or range of model parameters used in different numerical experiments.

Table G1 Stochastic stream power model parameters.

Table G2 Parameters used when running numerical experiments with the process-based routing model to constrain α and β.

Table G3 Parameters used when running numerical experiments with the empirical routing model to constrain α and β.

Table G4 Parameters used in the sensitivity analysis with the process-based routing model.

Table G5 Parameters used in the sensitivity analysis with the empirical routing model.