the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
A landslide runout model for sediment transport, landscape evolution, and hazard assessment applications
Jeffrey Keck
Erkan Istanbulluoglu
Benjamin Campforts
Gregory Tucker
Alexander HornerDevine
We developed a new rulebased, cellularautomaton algorithm for predicting the hazard extent, sediment transport, and topographic change associated with the runout of a landslide. This algorithm, which we call MassWastingRunout (MWR), is coded in Python and implemented as a component for the package Landlab. MWR combines the functionality of simple runout algorithms used in landscape evolution and watershed sediment yield models with the predictive detail typical of runout models used for landslide inundation hazard mapping. An initial digital elevation model (DEM), a regolith depth map, and the location polygon of the landslide source area are the only inputs required to run MWR to model the entire runout process. Runout relies on the principle of mass conservation and a set of topographic rules and empirical formulas that govern erosion and deposition. For the purpose of facilitating rapid calibration to a site, MWR includes a calibration utility that uses an adaptive Bayesian Markov chain Monte Carlo algorithm to automatically calibrate the model to match observed runout extent, deposition, and erosion. Additionally, the calibration utility produces empirical probability density functions of each calibration parameter that can be used to inform probabilistic implementation of MWR. Here we use a series of synthetic terrains to demonstrate basic model response to topographic convergence and slope, test calibrated model performance relative to several observed landslides, and briefly demonstrate how MWR can be used to develop a probabilistic runout hazard map. A calibrated runout model may allow for regionspecific and more insightful predictions of landslide impact on landscape morphology and watershedscale sediment dynamics and should be further investigated in future modeling studies.
 Article
(15836 KB)  Fulltext XML

Supplement
(1025 KB)  BibTeX
 EndNote
Over geologic timescales, landslides and their runout shape the topographic expression of mountain ranges and channel networks (e.g., Campforts et al., 2022; Korup, 2006; Larsen and Montgomery, 2012; Montgomery and Dietrich, 1988). Over more pragmatic engineering and environmental risk management timescales, landslides and their runout can inundate and destroy infrastructure (e.g., Kean et al., 2019) but also support numerous ecosystem benefits, including carbon and nutrient transport from hillslopes to channels and the creation of riparian habitat (Benda et al., 2003; Bigelow et al., 2007; Goode et al., 2012). Therefore, explicit representation of landslide runout is a necessary component of (1) landslide inundation hazard assessments, with an emphasis on inundation extent and flow depth (e.g., Frank et al., 2015; Han et al., 2015); (2) watershed sediment yield models, with an emphasis on the mobilization, deposition, and type of sediment carried by the landslide (e.g., Bathurst and Burton, 1998; Istanbulluoglu, et al., 2005); and (3) landscape evolution models, with an emphasis on topographic change prediction (e.g., Tucker and Bras, 1998; Istanbulluoglu and Bras, 2005; Campforts et al., 2022).
Landslide runout processes can be generalized into three phases: initiation, erosion, and deposition. After a landslide initiates, it may break apart and flow as a relatively dry debris slide, or it may mix with surface runoff to become a debris flow. The mobility of the masswasting material and resulting erosion–deposition pattern often varies as a function of runout topography and the initial relief and size of the landslide (Iverson, 1997). Mobility may also be impacted by substrate liquefaction (Hungr and Evans, 2004) and landslide basal cataclasis (Shaller et al., 2020). As the runout material moves downslope, flow depth varies as a function of channel width (Kean et al., 2019), which in turn impacts erosion rates (Schürch et al., 2011). Theoretical, field, and laboratory observations indicate that erosion rates may also depend on the moisture content of the channel bed (Iverson, 2012; McCoy et al., 2012), flow grain size (Egashira et al., 2001), and granular stress within the flow (Capart et al., 2015). The slope at which deposition begins is controlled by the graintowater ratio and friction angle of the runout material (Takahashi, 2014; Major and Iverson, 1999; Zhou et al., 2019), but the friction angle of the runout material may vary as a function of the grains in the flow and fluidization (Hutter et al., 1996). Lateral levees often form along the edges of the flow (Major, 1997; Whipple and Dunne, 1992; Shaller et al., 2020) and deposition at the distal end of the flow may occur as layered accretion (Major, 1997) or as the emplacement of a single, massive deposit (Shaller et al., 2020). If the water content of the runout material is high enough, as the solid fraction of the distal end of the flow compresses, the water is squeezed out and may continue as an immature debris flow (sensu Takahashi, 2014) or intense bedload (sensu Capart and Fracarollo, 2011), extending the runout distance (e.g., Shaller et al., 2020).
Landslide inundation hazard models aim to accurately predict the runout extent and/or flow depths of a runout event and may include some or most of the above processes in the model. Example models include (1) sitespecific empirical and statistical models that use simple geometric rules and an estimate of the total mobilized volume (initial landslide body + eroded volume) or a growth factor (e.g., Reid et al., 2016); (2) detailed, continuumbased mechanistic models, which conceptualize the runout process as a singlephase or multiphase flow using the depthintegrated Navier–Stokes equations for an incompressible, freesurface flow (i.e., shallowwater equations; Frank et al., 2015; Han et al., 2015; Iverson and Denlinger, 2001; Medina et al., 2008) and often (though not always) require preknowledge of the total mobilized volume (e.g., Barnhart et al., 2021; Han et al., 2015); (3) reduced or appropriatecomplexity flowrouting models (e.g., Murray, 2007) that use rulebased abstractions of the key physical processes that control the flow (Clerici and Perego, 2000; Guthrie and Befus, 2021; Gorr et al., 2022; Han et al., 2017, 2021; Horton et al., 2013; Liu et al., 2022) and are typically implemented using just the initial landslide location and volume but often rely on heavy, sitespecific parameterization; and (4) hybrid modeling approaches that combine mechanistic models with empirical and reducedcomplexity approaches (D'Ambrosio et al., 2003; Iovine et al., 2005; Lancaster et al., 2003; McDougall and Hungr, 2004).
For landscape evolution and watershed sediment yield applications (herein collectively referred to as watershed sediment models, WSMs), the runout model must be scalable in both space and time and capable of modeling the entire runout process given an internally modeled initial landslide body (e.g., Tucker and Bras, 1998; Doten et al., 2006; Campforts et al., 2022). As such, computationally efficient and parsimonious reducedcomplexity runout models that evolve the terrain and transfer sediment are often preferred in WSMs but with simplifications that can restrict model ability to accurately replicate observed inundation extent or depositional patterns. Such simplifications include omitting debris flow erosion and bulking in runout channels, limiting flow to only a single cell in the steepest downstream direction, and assuming debris flows only occupy the width of a single cell (e.g., Tucker and Bras, 1998; Istanbulluoglu and Bras, 2005) or link of a channel network (Benda and Dunne, 1997).
We developed a new reducedcomplexity landslide runout model, called MassWastingRunout (MWR), that bridges the scalable functionality of WSMs with the predictive accuracy of landslide inundation hazard models, without the computational overhead of a detailed mechanistic representation of the runout process or difficult parameterization typical of other models. MWR models landslide runout starting from the source area of the landslide, making it easily compatible with WSMs that internally determine the initial landslide body size and location. MWR tracks sediment transport and topographic change downstream and evolves the attributes of the transport material, making it suitable for sediment yield studies. MWR can be calibrated by adjusting just two parameters (S_{c} and q_{c}, described in Sect. 2) and is augmented with a Bayesian Markov chain Monte Carlo (MCMC) calibration utility that automatically parameterizes model behavior to observed runout characteristics (e.g., erosion, deposition, extent). MWR also includes a builtin utility called MWR Probability, designed for running an ensemble of simulations to develop probabilistic landslide runout hazard maps, making MWR suitable for hazard assessment applications.
In this paper, we present the conceptualization and numerical implementation of the MWR model (Sect. 2), describe the calibration utility and probabilistic implementation of MWR (Sect. 3), and demonstrate basic model response to topographic convergence and slope on a series of synthetic terrains (Sect. 4). Eventscale applications to replicate observed runout extent, sediment transport, and topographic change at four topographically and geologically unique field sites (see Fig. 1) are discussed (Sect. 5). We test MWR's predictive ability using the parameterization of one site to predict runout hazard at a nearby site and show a brief example of Monte Carlo model runs to determine runout probability from initial landslide source areas defined by an expertdetermined potentially unstable slope or a hydrologically driven landslide hazard model (Sect. 6). We conclude with a short summary of MWR model performance and suggest how a calibrated MWR can be incorporated into WSMs.
2.1 Overview of the cellularautomaton modeling approach
MWR is coded as a discrete cellularautomaton (CA) model. CA models apply a set of equations or rules (deterministic or probabilistic) to individual cells of a grid to change the numerical or categorical value of a cell state (e.g., Codd, 1968). In earth sciences, CA models are widely used to model everything from vegetation dynamics (e.g., Nudurupati et al., 2023) and lava flows (e.g., Barca et al., 1993) to geomorphic transport, in which gravitationally directed erosion and depositional processes modify a digital elevation model (DEM) representation of a landscape (e.g., Chase, 1992; Crave and Davy, 2001; Murray and Paola, 1994; Tucker et al., 2018). Existing CAbased landslide runout models include models by Guthrie and Befus (2021), D'Ambrosio et al. (2003), and Han et al. (2021). In all of these models, runout behavior is controlled by topographic slope and rules for erosion and deposition, but conceptualization and implementation differ.
In MWR, mass continuity is central to model conceptualization. Of the wide range of landslide runout processes described in the Introduction, MWR explicitly represents erosion, deposition, and flow resistance as a function of the initial landslide body and downslope terrain. Material exchange between the runout material and underlying regolith as well as flow resistance determines runout extent and landscape evolution. Model rules are designed such that they can be parameterized from field measurements. Finally, in MWR, most computations occur only at the location of moving debris, in a manner analogous to the “mobile” cellularautomaton implementation of Chase (1992).
Chase (1992) modeled precipitationdriven surface erosion by randomly placing single packets of precipitation on a DEM, which then moved from higherelevation to lowerelevation grid cells, eroding and transporting sediment as a function of the slope between the cells. The individual packets of precipitation were referred to as precipitons. In MWR, since we route the downslope progression of debris from a specified masswasting source area, we refer to these packets of debris as “debritons”. The debritons represent debris flux, here defined as a volume of debris transferred per model iteration per grid cell area, [m^{3} m^{−2} per iteration], and are equivalent to the flow depth in the cell.
The present implementation of the MWR algorithm is coded in Python and developed as a component of the Landlab earth surface modeling toolkit (Barnhart et al., 2020; Hobley et al., 2017). MWR uses the Landlab raster model grid, which consists of a lattice of equally sized, rectangular cells. Topographic elevation, derived topographic properties like slope and curvature, and other spatially varying attributes such as regolith depth and grain size are recorded at nodes in the center of each cell (see Fig. 5 of Hobley et al., 2017). In the subsequent sections we describe the model theory. All parameters and variables used in the theory are listed in Table B1.
2.2 Mobilization of the initial masswasting source material (Algorithm 1)
To initiate MWR, the user provides maps of initial topography, regolith depth, and the location and depth of the masswasting source material (e.g., the initial landslide body). Each raster model grid node in the masswasting source material is designated as a debriton (Fig. 2, iteration t=0) with a magnitude equal to the masswasting source material depth at the node and basal elevation equal to the initial topography minus the masswasting source material depth at the node. The basal elevation can be thought to represent the rupture or slip surface of the source material, and the redistribution (flux) of each debriton to its downslope nodes (receiver nodes) is determined as a function of the slope of the slip surface. At the lowestelevation debriton of the source material, flux to its downslope nodes is determined using the surface slope of the initial DEM (see the flow direction of the lowest node in Fig. 3a). This implementation helps to ensure that the lowestelevation debriton in the masswasting source material moves downslope and movement of upslope debritons is impacted by the geometry of the masswasting source material. For example, the receiver nodes of the lowestelevation debriton in the initial landslide body illustrated in Fig. 2 (node 45) would be identified as those among the eight neighboring nodes whose initial topographic elevation is less than the initial topographic elevation of node 45, while for the debriton at node 51, the receiver nodes would be identified as those among the eight neighboring nodes whose topographic elevation is less than the basal elevation of the debriton (see Fig. 3a).
2.3 Flow routing and rules for erosion, deposition, and resistance (Algorithm 2)
Algorithm 2 is essentially the runout model. It determines how each debriton traverses and modifies the landscape. After receiver nodes from the first model iteration are determined in Algorithm 1 (iteration t=0), Algorithm 2 is repeatedly implemented until all material has deposited (i.e., there are no debritons). Each debriton moves one grid cell per model iteration; the larger the landslide size, the more iterations are necessary to evacuate the landslide slip surface. As each debriton moves, it may erode or aggrade the landscape, impacting the movement of any upslope debritons. As is common with other reducedcomplexity models (e.g., Guthrie and Befus, 2021), we assume that inertial effects have a negligible impact on flow behavior (i.e., the kinematic flow approximation) and the downslope redistribution of a debriton or flux to each of a node's ith receiver nodes (${q}_{{R}_{i}}$) is determined as a function of topographic slope (slope of the terrain under the debriton). We do this using the Freeman (1991) multipleflowdirection algorithm:
where q_{O} is the total outgoing flux from the node and has units of depth in meters per model iteration, Nr is the number of receiving nodes, i is the index for each receiver node (e.g., $i=\mathrm{1},\phantom{\rule{0.125em}{0ex}}\mathrm{2}\mathrm{\dots}\text{Nr}$), and S_{i} is the underlying topographic slope to the ith receiver node (Fig. 3b). The Freeman (1991) multipleflowdirection algorithm is a commonly used approximation for twodimensional flow, and in this implementation it is handled by a preexisting Landlab flowrouting component. The exponent a controls how material is distributed to downslope nodes, with higher values causing narrower flow (Holmgren 1994). In a braided river cellularautomaton model, Murray and Paola (1997) used an approximation for turbulent shallowwater flow to justify a=0.5 (which is the exponent on the slope factor in channel friction laws). For our application, we found that MWR provided a closer fit to observed masswasting runout if a=1, suggesting that the material behavior is more similar to linear–viscous shear flow than to wallbounded turbulent shear flow (e.g., as the runout debris flows downslope, it tends to spread less than shallow, turbulent water). The total incoming flux (again, in units of meters per model iteration) towards a given node (q_{I}) is determined by summing the flux from each of the node's donor nodes:
where Nd is the number of donor nodes, and ${q}_{{D}_{j}}$ is the flux from node D_{j} (the jth donor node, $j=\mathrm{1},\mathrm{2}\mathrm{\dots}\text{Nd}$; Fig. 3b).
As noted by Tucker and Hancock (2010), the flow depths calculated from twodimensional flow approximations as in Eq. (1) can be influenced by the choice of grid size used to represent the terrain. Additionally, as simplified multidirectional flow models as in Eq. (1) neglect the pressure and momentum forces in the movement of flow, they can result in inaccurate flow width and depth estimates, depending on terrain slope and convergence. Rengers et al. (2016) noted these limitations when using a kinematic wave approximation of the shallowwater equations, as this approximation lacks a pressure term that facilitates the spreading of the modeled water surface. While the topographic controls on mass conservation are adequately represented by Eq. (1), our model bears such limitations when calculating flow depth and width. Moreover, in our model, flow depth is used to determine a depthdependent erosion rate. As such, in order to avoid unrealistically high flow depths (and erosion rates), we constrain flow depth to an upper limit as
where h_{max} is an effective upper limit to flow depth that in practice can be approximated as the maximum observed flow depth, as inferred from field indicators or assigned based on expert judgment (see Sect. 5), and h is the corrected flow depth used to calculate flow shear stress. This correction allows erosion rates to vary with flux but prevents unreasonably large values. This flow depth correction does not violate the conservation of mass and runout mass balance, as h is only used to calculate flow shear stress.
To determine aggradation (A) at a node, we use a critical slope (S_{c}) constraint that permits computationally rapid distribution of q_{I} over multiple nodes. Critical slope constraints or rules are common to many reducedcomplexity and landscape evolution models. Chen et al. (2023) showed that when flow inertia can be ignored, S_{c} can be approximated from the surface slope of observed deposits. Several landscape evolution models use an S_{c}based nonlinear, nonlocal aggradation scheme (e.g., Campforts et al., 2020; Carretier et al., 2016), but when this rule is implemented with the debriton framework described above, unreasonably tall deposits result when q_{I} is large and the slope at the node (S)≪S_{c}. To resolve this problem, aggradation depth can be limited to A≤S_{c}Δx (where Δx is the grid cell length), but we found that this constraint results in long deposits that parallel the underlying slope when q_{I} is large. Instead, MWR computes the aggradation depth at a node assuming that the aggradation will spread over N_{a} nodes until all of q_{I} is deposited and that the surface slope of the overall deposit will be equal to S_{c}, as shown in Fig. 4 and described as follows.
Aggradation at a node is determined as
where S is the steepest slope to the node's eight neighboring nodes, and ${A}_{\mathrm{p},{N}_{\mathrm{a}}}$ is a potential aggradation depth necessary to form the beginning of the overall deposit that (1) begins at the node and spreads over N_{a} consecutive nodes, (2) has a total volume equal to q_{I}Δx^{2}, (3) has a surface slope equal to the critical slope S_{c}, and (4) has an underlying topographic slope equal to S that is assumed to be constant over the N_{a} consecutive nodes of deposition. From these assumptions, we can analytically define ${A}_{\mathrm{p},{N}_{\mathrm{a}}}$ and N_{a} as a function of q_{I}, S_{c}, and S as follows.
First, q_{I}, calculated from Eq. (2), can be used to calculate A_{p,i} by expressing q_{I} as the sum of the N_{a} deposits that make up the overall deposit as
where A_{p,i} is the ith deposition amount in the overall deposit and i=1 is the last node of deposition (A_{p,1}; see Fig. 4). Since we assume the deposit slope and underlying topographic slope are uniform, the deposition amount at any of the N_{a} nodes can be determined from A_{p,1} as
From Eq. (6) we can rewrite Eq. (5) as a function of A_{p,1} and rearrange to define A_{p,1} as a function of q_{I}:
Substituting Eq. (7) into Eq. (6) and solving for i=N_{a}, we get an expression for ${A}_{\mathrm{p},{N}_{\mathrm{a}}}$:
Equation (8) can be rearranged into a quadratic equation and solved for N_{a} as
We use Eq. (8) to solve for ${A}_{\mathrm{p},{N}_{\mathrm{a}}}$ and Eq. (9) to solve for N_{a} assuming ${A}_{\mathrm{p},\mathrm{1}}=\mathrm{1}/\mathrm{2}\mathrm{\Delta}x{S}_{\mathrm{c}}$ and rounding the positive solution to the nearest integer. When implemented using a single debriton released on a twodimensional hillslope as illustrated in Fig. 4, the debriton deposits over N_{a} nodes at a uniform slope equal to S_{c}. When implemented on an actual threedimensional terrain, the interaction between multiple debritons in multiple directions creates a complex deposit whose slope changes with S_{c}.
To determine erosion depth (E) [m per iteration], we constrain E to the lesser of a potential erosion depth, h_{e}, and local regolith depth, h_{r}:
where h_{e} is computed as a function of the basal shear stress of the flow, τ [Pa] (Eqs. 12 and 13), and the critical shear stress (τ_{c}) of the regolith at the node [Pa],
The coefficient k is an erodibility parameter [m Pa^{−f}]. Stock and Dietrich (2006) showed that k encapsulates substrate properties. If h_{e} is used to represent erosion over geomorphic timescales, with repeated debris flow occurrences in a single model iteration, k becomes associated with debris flow length and frequency (Perron, 2017). In our application since we are modeling the erosion associated with a single runout event, as represented by the downslope movement of the debritons, the coefficient k therefore needs to scale h_{e} on the order of the average erosion depth caused by a single debriton. Using this logic, k can be computed using the observed average erosion depth and an estimated length of the runout material that caused the erosion. Further details on how we determine k from observed runout are included in Appendix A. The exponent f controls the nonlinearity of h_{e} with shear stress. Many authors (e.g., Chen and Zhang, 2015; Frank et al., 2015; Shen et al., 2020) use a value of 1 for f but field measurements by Schürch et al. (2011) (see their Fig. 3) suggest that f may be less than 1 if τ is assumed to vary linearly with flow depth, particularly at flow depths greater than 3 m.
MWR includes two options for defining τ: (1) a quasistatic basal shear stress approximation or (2) a grainsizebased shear stress approximation. The quasistatic basal shear stress approximation (e.g., Takahashi, 2014) is defined as
where ρ is the density of masswasting material (grain and water mixture) [kg m^{−3}], g is gravity [m s^{−2}], h is the adjusted flow depth described in Eq. (3), and θ is the topographic slope (tan ^{−1}(S)) measured in degrees.
The grainsizebased shear stress approximation is defined using an empirical formula by Bagnold (1954):
where σ is normal stress [Pa], and φ is the collision angle between grains measured from the vertical axis (see Bagnold, 1954) with a value of tan φ typically equal to 0.32. Stock and Dietrich (2006) defined σ as
where υ_{s} is the volumetric solid concentration, ρ_{s} is density of the solids [kg m^{−3}], u is flow velocity [m s^{−1}], z is depth below the flow surface [m], $\mathrm{d}u/\mathrm{d}z$ is the shear strain rate [s^{−1}], and D_{s} is the representative grain size [m]. Stock and Dietrich (2006) suggested that D_{s} corresponds to a small percentile of the coarsest fraction of the runout material (D_{88} to D_{96}) and they approximated $\mathrm{d}u/\mathrm{d}z$ as
Solely for the purpose of computing $\mathrm{d}u/\mathrm{d}z$, we approximate velocity at a node using a grainsizedependent empirical formula for debris flow velocity by Julien and Paris (2010) as
where u^{∗} is shear velocity $\left(\sqrt{gh\mathrm{tan}\mathit{\theta}}\right)$. Substituting Eqs. (16), (15), (14), and (13) into Eq. (11) yields a grainsizedependent approximation for h_{e} that mimics the nonlinear erosion response to flow depth in Schürch et al. (2011). Additionally, this form of τ is advantageous because it permits landslidedriven erosion rates to scale with landslide grain size, which can vary by lithologic region (e.g., RodaBoluda et al., 2018). As will be shown in Sect. 5, we obtained reasonable model calibration at multiple sites by defining D_{s} from the coarser grain sizes observed in the field at existing runout deposits, road cuts, and treethrow pits.
Once A [m] and E [m] have been determined, total outgoing flux per iteration, q_{O} [m], is determined as (see Fig. 3c)
where q_{c} is a threshold flux for deposition. When q_{I}<q_{c}, q_{I} deposits and q_{O} becomes zero. The threshold flux q_{c} conceptually represents the flow depth below which flow resistance is large enough to cease the forward momentum of the flow, whether in the form of internal friction or friction due to vegetation and obstructions (e.g., large clasts or logs). The density and water content of q_{I}, A, and E are treated as uniform, and surface runoff, such as channelized stream flow or hillslope infiltration excess runoff, that might mix with q_{I}, A, or E is ignored. Once q_{I}, A, q_{O}, and E have been determined, change in elevation at a node (Δη) is calculated as
Attributes (e.g., grain size, organic content, or any other attribute that is transferred in the flow) of the debriton and regolith are updated using a volumetricweighted average approach. First, for each regolith attribute being tracked by the model, the attribute value delivered to a node from its donor nodes (ξ_{D}) is determined as
where q_{D} is a vector containing each ${q}_{{D}_{j}}$ sent to the node, ξ_{D} is a vector containing the incoming attribute values for each ${q}_{{D}_{j}}$, and q_{I} is the sum of incoming flux from donor nodes defined by Eq. (2).
Second, the attribute value sent from a node to its receiver nodes (ξ_{R}) is determined as
where ξ_{t−1} is the attribute value at the node before any aggradation (i.e., the previous iteration attribute value). Finally, the attribute value at the node, updated to account for erosion and aggradation (ξ), is
Regolith thickness (h_{r}) and topographic elevation (η) are updated at a node as
where η_{t−1} and h_{rt−1} are the topographic surface elevation and regolith thickness at the node from the previous model iteration. After regolith thickness and topographic elevation have been updated for each debriton, the multidirection slope of the DEM, which is needed for implementing Eq. (1) in the next model iteration, is recomputed.
As the DEM is updated following each model iteration, topographic pits or flat topography may form. These features have no slope or slope inwards and obstruct debriton movement. To allow a debriton to pass an obstruction, we rely on a simple workaround: upon encountering the obstruction, the debriton is directed to itself and some portion of the debris is deposited based on Eq. (4). At the end of the model iteration, the node elevation and slope are updated. During the next iteration, if the remaining mobile debris is no longer obstructed, it moves to its downslope node(s). If the node is still obstructed, it is again sent to itself until either all material has deposited or the elevation of the node exceeds that of its neighbor nodes, allowing the debriton to move downslope.
3.1 Calibration utility
MWR includes a calibration utility that uses an adaptive Bayesian Markov chain Monte Carlo (MCMC) algorithm described by Coz et al. (2014) and Renard et al. (2006). The calibration utility determines (1) a single set of parameters that best match MWR output to an observed landslide runout dataset and (2) posterior parameter probability distribution functions (PDFs). The observed runout dataset can consist of a single landslide or multiple landslides. Depending on user input, MWR simultaneously or sequentially models runout from each landslide source area in one model run. To use the calibration utility, the user provides an initial (prior) guess of the parameter values and their respective PDFs that calibrate the MWR to a specific site. Then, the calibration utility randomly selects a set of trial parameter values (Λ) from the prior PDFs and runs MWR using Λ. Once the model has completed the run, the algorithm evaluates the posterior likelihood of the parameter set (L(Λ)) as the product of model ability to replicate observed runout (described below) and the prior likelihood of the parameter set. After the first L(Λ) has been determined, the utility selects a new set of parameters (Λ_{t+1}) by jumping some distance (described below) from each parameter in Λ space. Depending on the value of L(Λ_{t+1}), the algorithm either stays at Λ or moves to Λ_{t+1}. This Markov process is repeated a userspecified number of times. The jump direction is random, but the algorithm is adaptive because the jump distance changes depending on whether $L\left({\mathrm{\Lambda}}_{t+\mathrm{1}}\right)>L\left(\mathrm{\Lambda}\right)$ occurs more than a userspecified threshold value. For a detailed description of the algorithm see Coz et al. (2014).
The L(Λ) index is estimated as the product of the prior probability of the selected parameter values, p(Λ), and three other performance metrics as
where Ω_{T} is an index for evaluating model planimetric fit described by Heiser et al. (2017), and Δη_{E} and ${Q}_{{\mathrm{s}}_{\mathrm{E}}}$ are new dimensionless indices proposed for this study (described below). The index Δη_{E} is the volumetric error of the modeled topographic change over the entire model domain normalized by the observed total mobilized volume (initial landslide body + erosion volume). The index ${Q}_{{\mathrm{s}}_{\mathrm{E}}}$ is the mean cumulative sediment transport error along the modeled runout path normalized by the observed mean cumulative flow. Larger values of Ω_{T} and smaller values of Δη_{E} and ${Q}_{{\mathrm{s}}_{\mathrm{E}}}$ indicate that modeled runout more closely fits observed. Note that we add a value of 1 to Ω_{T} and use the squared reciprocal values of Δη_{E} and ${Q}_{{\mathrm{s}}_{\mathrm{E}}}$ in Eq. (24) so that the magnitude of L(Λ) is always equal to or greater than zero and increases with improved fit. The metric Ω_{T} is written as
where α, β, and γ are the areas of matching, overestimated, and underestimated runout extent, respectively.
The index Δη_{E} is determined as
where V is observed total mobilized volume and p is the number of nodes in the area made up of the matching, overestimated, and underestimated areas of modeled runout extent. Δη_{Mi} and Δη_{Oi} are the modeled and observed topographic change [m] at the ith node within that extent.
To calculate ${Q}_{{\mathrm{s}}_{\mathrm{E}}}$, we first determine the cumulative sediment transport volume (Q_{s}) at each node, j along the runout profile, in a manner similar to the flow volume–mass balance curves in Fannin and Wise (2001) and Hungr and Evans (2004):
where Δη_{i,j} is the topographic change [m] at the ith node located above the elevation of node j, and u_{j} is the total number of nodes located above node j. Q_{s} is computed for both the observed and modeled runout (Q_{sO} and Q_{sM}, respectively) and ${Q}_{{\mathrm{s}}_{\mathrm{E}}}$ is determined as
where r is the number of nodes along the runout profile, and $\stackrel{\mathrm{\u203e}}{{Q}_{\text{sO}}}$ is the observed mean cumulative flow.
As will be detailed in Sect. 5, field estimates for S_{c} and q_{c} vary over the length of the runout path. To account for the heterogeneity of S_{c} and q_{c}, we estimate prior PDFs of potential S_{c} and q_{c} values from field and/or remote sensing measurements. Then, from model calibration to a DEM of difference (prerunout DEM subtracted from the postrunout DEM; DoD) using the calibration utility, we produce posterior PDFs of S_{c} and q_{c} and find a single S_{c} and q_{c} pair that allows the modeled DoD to best replicate the observed DoD (i.e., the S_{c} and q_{c} pair with the highest L(Λ) value).
We run the calibration utility using a single Markov chain of 2000 repetitions. At most sites, the model converged relatively quickly on a solution and we therefore did not account for burnin or evaluate convergence (e.g., Gelman et al., 2021) and considered 2000 repetitions adequate. Future implementations of the calibration utility may include multiple chains, burnin, and a check for convergence. As a final note, many debris flow runout models are evaluated using Ω_{T} or variations of Ω_{T} alone (e.g., Gorr et al., 2022; Han et al., 2017) and the MWR calibration utility can also be run solely as a function of Ω_{T} (i.e., runout extent). However, we found that calibration based on Ω_{T} alone results in high parameter equifinality (e.g., Beven 2006); multiple parameter sets result in an equally calibrated model as evaluated by Ω_{T}. As such, we recommend calibrating debris flow and landslide runout models to an observed DoD. If repeated lidar is available, a DoD can be obtained from before and after scans of the observed runout event. Alternatively, a DoD can be created by hiking the observed runout event and mapping fieldinterpreted erosion and deposition depths. Additional details on how we prepared DoDs for multiple sites are included in the Supplement.
3.2 Mapping landslide runout hazard
MWR includes an additional utility called MWR Probability that produces landslide runout probability maps. MWR Probability repeatedly runs MWR a userspecified number of times (Np), each repetition with a different, randomly sampled parameter set from the posterior parameter PDFs produced by the calibration utility. MWR Probability includes three options for specifying the initial masswasting source material: (1) a userprovided landslide source area polygon based on field and/or remote sensing observations; (2) a userdefined hillslope that is susceptible to landslides (e.g., potentially unstable slope), where landslide area and location are randomly selected within but no larger than the hillslope – this option is useful when the extent of a potential landslide is unknown; and (3) a series of mapped landslide source areas within a watershed, as determined by an externally run Monte Carlo landslide initiation model (e.g., Hammond et al., 1992; Strauch et al., 2018) – this option is useful for regional runout hazard applications. If using Option 1, modeled runout probability represents uncertainty in MWR parameterization. If using Option 2 or 3, modeled runout probability reflects uncertainty in both MWR parameterization and landslide location and size.
For all three run options, each model iteration begins with the same initial topography. After Np model simulations, Np different versions of the postrunout landscape are created, and the probability of runout at each node is determined as
where num(Δη>0) is the number of times topographic elevation at a node changes as a result of erosion or deposition from the Np model runs. Probability of erosion or aggradation can be determined by replacing the numerator in Eq. (29) with num(Δη<0) or num(Δη>0), respectively.
We evaluate basic model behavior using a series of virtual experiments. The virtual experiments consist of six synthetic terrains, including (a) a planar slope that intersects a gently sloped plane (S=0.001), (b) a planar slope with a constriction that intersects a gently sloped plane, (c) a planar slope that has a bench midslope and then intersects a gently sloped plane, (d) a concaveup uniformconvergence slope, (e) a concaveup variableconvergence slope that widens (convergence decreases) in the downslope direction, and (f) a convexup variableconvergence slope that widens (convergence decreases) in the downslope direction. On each terrain, a 30 m wide, 50 m long, and 3 m deep landslide is released from the top of the terrain. All six terrains are covered by a 1 m thick regolith and use the same parameter values (S_{c}=0.03, q_{c}=0.2 m, k=0.01, D_{s}=0.2 m). Each terrain is represented using a 10 m grid. Experiment results are shown in Fig. 5.
On Terrain A, the landslide spread as it moved downslope and formed levees along the edge of the runout path. The width of the spread was a function of the multipleflowdirection algorithm and resistance along lateral margins of the runout as represented by q_{c}. At the slope break at the base of the slope, the material deposited at an angle controlled by S_{c}. On Terrain B, the flow initially eroded and deposited material identical to Terrain A, but near the slope break, the topographic constriction forced flow depth to increase and exceed q_{c}, minimizing the formation of levees (because q_{O}>q_{c}), and resulted in a slightly larger deposit at the base of the slope. On Terrain C, landslide runout was again initially identical to the runout on Terrain A; however, upon intersecting the midslope bench, most of the runout material deposited. A small, thinner portion did continue past the bench but eroded at a lower rate than upslope of the bench. Upon intersecting the flat surface at the base of the hillslope, the runout material was deposited.
On Terrain D, the landslide and its runout were confined to the center of the convergent terrain and were only deposited once the slope was less than S_{c}. The slide never widened because the uniformly convergent channel shape prevented spreading and the narrower flow width maintained a higher flow depth, which prevented the formation of levees. On Terrain E, the landslide again deposited once the slope was less than S_{c}, but because topographic convergence of Terrain E decreases in the downslope direction, as the runout material moved downslope, the deposit spread more than on Terrain D, which caused thinner flow and deposition along margins of the runout path. On the final terrain, Terrain F, the slope is always greater than S_{c}, so deposition was limited to levees along the edge of the flow that formed as the runout spread in response to decreasing convergence.
MWR model behavior can be summarized as follows: the displacement and deposition of landslide material predicted by MWR respond to topography in a reasonable manner. Flow width increases as convergence decreases (e.g., Terrain F), which in turn reduces flow depth. Lower flow depths cause lower erosion rates and reduce aggradation extent. Conversely, modeled flow depth increases when convergence increases (e.g., Terrain B). Where the flow encounters broadly convergent or planar slopes, lateral levee deposits form, a common feature of landslides reported in the literature and at sites reported here (see Sect. 5) that detailed mechanistic models can struggle to reproduce (e.g., Barnhart et al., 2021).
We do not attempt to compare MWR modeled flow with the output of shallowwaterequationbased models or observed granular flows (e.g., Medina et al., 2008; McDougall and Hungr, 2004; Iverson and Denlinger, 2001; Han et al., 2015). The cellularautomaton representation in MWR does not model the timedependent evolution of debris flow velocity and depth, and it conceptually moves debris instantaneously at each iteration, as driven by changes in the evolving topographic elevation field. Because of that, only the final outcome (modeled runout extent, sediment transport, and topographic change) of MWR can be compared with other models or observed runout, which we do in the next section. Also, as described in Sect. 2.3, the behavior of the multiple flow direction algorithm does vary with grid size. Using a coarser or finer grid, without adjusting model parameterization, could potentially change the runout patterns shown in Fig. 5.
5.1 Overview
In this section, we demonstrate the ability of a calibrated MWR to replicate observed runout extent, sediment transport, and topographic change at field sites located in the western USA and summarize model calibration results with an evaluation of MWR calibration relative to terrain attributes of the observed runout paths. Note that simply calibrating a model to match field data does not constitute a satisfactory test of model predictive ability (Iverson, 2003). Strategic testing, which involves calibrating the model to one site or period of time and then running the calibrated model at a separate site or period of time (Murray, 2013), is a better indicator. Two of our validation sites, the Cascade Mountains and Olympic Mountains sites, include two separate landslides and subsequent runout and we test model predictive ability at these sites in Sect. 6.
Calibrated model performance is demonstrated at the following field sites (see Fig. 6a for locations and observed runout extent): (1) two runout events over the same hillslope in the Cascade Mountains (Washington state – WA, USA) – a large debris avalanche in 2009 (Cascade Mountains, 2009) and a moderately sized debris flow in 2022 (Cascade Mountains, 2022) that inundated and flowed within a first to secondorder channel until perpendicularly intersecting a narrow river valley several hundred meters below the landslide (Fig. 1a); (2) debris flows in the Black Hills (WA) sourced from a small failure along the toe of a deepseated landslide (Black Hills, south) and a moderately sized debris avalanche from a large road fill (Black Hills, north) that flowed several kilometers along a relatively wide, broadly convergent channel before stopping (Fig. 1b); (3) a single, moderately sized debris avalanche in the Rocky Mountains (Rocky Mountains), the majority of which flowed several hundred meters over a broadly convergent to divergent hillslope in Colorado (Fig. 1c); and (4) a 30year chronology of small landslides and subsequent debris flows in the Olympic Mountains (WA) in steep, highly convergent channels that flowed well over a kilometer and coalesced into a single runout deposit in a dendritic, channelized watershed (Olympic Mountains; Fig. 1d). All landslides were initiated during heavy rainfall or rainplussnowmelt storm events (WRCC, 2024; NRCS, 2022; Table 1), but their runout varied in terms of erosion rate, grain size (Fig. 6b), depositional behavior (Fig. 6c), and the topographic convergence of the underlying terrain.
5.2 Model setup and field parameterization
Each model was set up on a 10 m grid representation of the preevent DEM with either a uniform or spatially varying regolith thickness (detailed for each site in the Supplement). The length (l) and area of the initial masswasting source material (e.g., the initial landslide body) were interpreted from a combination of lidar DEM, air photos, and field observations. The average depth of the initial landslide body was measured in the field or from the DoD. The volume of the initial landslide body was determined as the area times the average depth. An average width was determined as the area divided by the length. At the Olympic Mountains site, where the observed runout pattern formed as a result of multiple landslides (see the Supplement), landslide depth and width values listed in Table 1 are average values, and landslide length, area, and volume values are the average cumulative value upstream of each runout path. At all locations, we use Eq. (13) to approximate shear stress. We fieldsurveyed each site, noting the maximum flow depth (inferred from initial landslide body volume and height of scour marks as well as width of the channel in the erosion zone), typical deposition and erosion depths, and the size of the largest grains in the runout deposits.
We estimated parameter values from these field and remote observations (see Table 1). A sitespecific value for k was determined as a function of the observed average erosion depth (determined as total erosion volume divided by the erosion area, $\stackrel{\mathrm{\u203e}}{E}$) relative to the length of the runout debris, which we approximated as the length of the initial landslide body (l). Further details are described in the Appendix.
The volume of the initial landslide body ranged from 400 to 110 000 m^{3} across sites. At all sites, erosion and subsequent entrainment added to the total mobilized volume (initial landslide body + erosion volume), but the contribution was highly variable. The erosion volume divided by the total mobilized volume was as low as 0.19 for the Cascade Mountains 2022 landslide and as high as 0.97 for the Olympic Mountains landslides (Table 1).
The average maximum grain size varied from 0.2 m at the Black Hills sites to nearly 1 m at the Rocky Mountains site (Fig. 6b, Table 1). Values of $\stackrel{\mathrm{\u203e}}{E}/l$ ranged from 0.007 to 0.041 [m m^{−1}], with the highest rate occurring for the Rocky Mountains landslide and the lowest at the Black Hills sites. Details on grain size samples and data collected in the field are described in the Supplement. In terms of growth factors (average volumetric erosion per unit length of the erosiondominated region of the runout path; Hungr et al., 1984; Reid et al., 2016) values ranged from 10 m^{3} m^{−1} at the Black Hills south site to 95 m^{3} m^{−1} at the Rocky Mountains site (Table 1).
The median values of topographic slopes at which observed deposition occurred (i.e., Δη>0, inferred from the DoD) ranged between 0.1 and 0.3 across sites, while deposition was also observed in much steeper (>0.4) slopes and much flatter slopes at some sites (Fig. 6c). The slope of channel reaches where net deposition (cumulative erosion and deposition; e.g., Guthrie et al., 2010, inferred from field observations) was positive tended to be lowest at the Black Hills site (<1 % to 10 %) and highest at the Rocky Mountains site (16 % to 25 %).
We defined uniform prior distributions of S_{c} and q_{c} and then used the calibration utility to find the bestfit parameter values (parameter values corresponding to the highest L(Λ)). Minimum and maximum values of S_{c} were initially estimated from the range of observed slopes of areas of positive net deposition (Table 1). Minimum and maximum values of q_{c} were set as 0.01 to 1.75, which roughly represents the range of minimum observed thickness of debris flow termini in the field at all of the validation sites. For the purpose of implementing the calibration utility, we prepared a DoD of each site. The DoD was determined from either repeated lidar or field observations as detailed in the Supplement.
5.3 Calibration and model performance
Markov chains, colored according to the likelihood index, L(Λ), are plotted in the S_{c}–q_{c} domain, along with histograms of sampled S_{c} and q_{c} values for each landslide in Fig. 7. Each Markov chain includes 2000 model iterations. The runtime for 2000 model iterations depended on model domain, landslide size, and number of landslides modeled but varied from roughly 1.5 h for the Cascade Mountains 2022 landslide to 6 h for the Olympic Mountains site on a 2016 2.1 GHz Intel Core Xeon, 32 GB memory desktop. The chains show a wide array of sampling patterns and parameter ranges, but broadly speaking, at all sites, the algorithm jumped within S_{c}–q_{c} space towards higher L(Λ) to form bellshaped posterior distributions for each parameter. Depending on the landslide type, the calibration algorithm converged on different S_{c}–q_{c} pairs. For example, at the Cascade Mountains site, the calibration utility converged to smaller q_{c} and S_{c} values for the 2009 event (Fig. 7a), which permitted thinner flows over lower slopes and effectively made the 2009 modeled runout more mobile relative to the 2022 modeled runout (Fig. 7b). At the Rocky Mountains site (Fig. 7e), a relatively high q_{c} value helps limit the lateral extent of the modeled runout that in the observed runout was controlled by standing trees (see the Supplement).
Profile plots of modeled Q_{s} and maps of the modeled planimetric runout extent, colored to indicate where the runout matched (α), overestimated (β), or underestimated (γ) the observed runout, are shown in Fig. 8. Values of Ω_{T} we obtained with MWR are comparable to or higher than reported values of Ω_{T} in the literature that used a variety of models (Gorr et al., 2022; Barnhart et al., 2021; note that to compare Ω_{T} values to those studies, subtract 1 from values reported in this study). Across the sites, the volumetric error of the model, Δη_{E}, ranges between 6 % and 15 % (median 9.1 %) of the total mobilized volume. An overall <10 % volumetric error is reasonable considering the low number of parameters required to calibrate MWR and that empirical estimates of total mobilized volume used to run other runout models can vary by as much as an order of magnitude (e.g., Gartner et al., 2014; Barnhart et al., 2021). Model performance in predicting cumulative sediment transport along the runout profile was within similar error ranges. Except for the Rocky Mountains site where MWR consistently modeled widerthanobserved flow, the mean cumulative sediment transport error along the runout profile (${Q}_{{\mathrm{s}}_{\mathrm{E}}}$) was limited to between 5 % and 19 % of the mean cumulative flow determined from the observed DoD.
MWR generally successfully replicates observed sediment transport along the runout path via model parameterizations that are unique to each landslide. For example, the profile plots of Q_{s} at the Cascade Mountains site (Fig. 8a and b) show that during the 2009 landslide, all of the runout material flowed past the first 750 m of the runout path. During the 2022 landslide, material began to deposit just downslope of the initial landslide source area, as both observed Q_{s} and modeled Q_{s} reverse slope, indicating a decrease in cumulative flow. Model comparisons in the Cascade Mountains site were limited to the upper 750 m of the hillslope because a large portion of the runout material was lost to fluvial erosion in the valley (see the Supplement).
MWR also successfully replicates the observed sediment transport patterns at the Olympic Mountains site (profile plot of Q_{s} in Fig. 8f) and to a lesser degree the Rocky Mountains site (Fig. 8e). This finding is notable because at the Olympic Mountains and Rocky Mountains sites, observed runout extent and sediment depositional patterns were heavily impacted by woody debris or standing trees (see the Supplement).
Using a fixed grid size of 10 m might have impacted model performance in some areas like the Olympic Mountains and Cascade Mountains 2022 sites where MWR tended to overestimate the runout width (yellow zones in Fig. 8f and b). While a 10 m DEM is generally accepted as a good balance between model detail and computational limitations (e.g., Horton et al., 2013), for small landslides, the 10 m grid was close to the width of the channels that controlled observed runout (see Fig. 1d) and may not have accurately represented the terrain. Thus, modeled flow was less topographically constrained and tended to flow over a wider area than observed.
Also, because MWR does not have an explicit representation of flow momentum, it may show poor performance in regions of the runout path where flow momentum is the primary control on runout extent. For example, for the Cascade Mountains 2009 slide, MWR underestimates a short section of slopeperpendicular flow over a bench. Review of model behavior for this slide (Fig. 9) shows how MWR successfully mimics diverging flow around a broad ridge upslope of the bench (iteration t=28 in Fig. 9) but afterwards continues to follow topographic slope and converges too rapidly into a narrow ravine along the west edge of the bench (iteration t=40 in Fig. 9; compare to runout scar in the air photo and underestimated region on the topographic bench in Fig. 8a).
Nonetheless, overall, calibration was best for the Cascade Mountains 2009 landslide (values of Ω_{T} are highest and values of Δη_{E} and ${Q}_{{\mathrm{s}}_{\mathrm{E}}}$ are lowest) and poorest at the Rocky Mountains and Olympic Mountains sites (values of Ω_{T} are lowest and ${Q}_{{\mathrm{s}}_{\mathrm{E}}}$ and Δη_{E} are highest). At both the Rocky Mountains and Olympic Mountains sites, because we lacked repeat lidar, we created the DoD from a map of fieldestimated erosion and deposition depths and estimated the preevent DEM. The lower calibration scores may indicate that fieldestimated DoDs were not as accurate as those determined via lidar differencing. Another source of uncertainty that we have not addressed in our study is regolith thickness. At most sites, we used a uniform thickness. As regolith thickness limits maximum erosion depth (i.e., Eq. 10), using a spatially accurate regolith thickness may improve model performance. Finally, except for the Rocky Mountains site, where topography was unusually planar and the model seemed to consistently overestimate flow width, at most sites, MWR does not appear to have a strong systematic bias in modeled output, which suggests that when applied to convergent terrain, MWR may not have any structural weaknesses. In the next section, we evaluate model performance relative to runout path topography.
5.4 Runout path topography and model performance
Model behavior at the Rocky Mountains site suggests that MWR performance may systematically vary with topography (e.g., it may not perform as well on planar hillslopes). To check for systematic model variation, we compared model performance with three topographic indices described by Chen and Yu (2011). The indices are computed from the terrain in the observed runout extent and include the relief ratio ($H/L$), mean total curvature (κ), and mean specific stream power index (SPI). The index $H/L$ equals the average slope of the runout path (or relative relief), determined as the total topographic relief of the runout (measured from the center of the initial landslide to the end of the runout path) divided by the horizontal length of the runout, and indicates the mobility of the runout. The index κ represents topographic convergence, which is the second derivative of the terrain surface, with increasingly positive values of the index κ reflecting growing topographic convergence and a concaveup channel profile (e.g., Istanbulluoglu et al., 2008). The index SPI is determined as the natural log of the product of the contributing area and slope. Indices κ and SPI are computed at each node in the runout extent and model performance is compared to the mean value.
Comparison of model performance with respect to the topographic indices in Fig. 10 shows slightly improved model performance over runout paths that are less convergent (lower SPI and κ values) and on steeper terrain (higher $H/L$), but neither trend is significant. The latter finding appears to be mostly a result of how well modeled sediment transport and topographic change (${Q}_{{\mathrm{s}}_{\mathrm{E}}}$ and Δη_{E}) replicated observed change, as there does not appear to be a trend in Ω_{T} with $H/L$ and the two bestperforming models (both Cascade Mountains landslides) had the lowest (best) ${Q}_{{\mathrm{s}}_{\mathrm{E}}}$ values and low Δη_{E} values. Both findings are likely impacted by the grid size we used to represent terrain. As noted above, at all sites we used a 10 m grid, but at some sites 10 m does not quite capture the relief of channelized topography that controlled observed runout, leading to modeled runout that was considerably wider than observed and causing low Ω_{T} values (this is especially true at the Olympic Mountains site; Fig. 10a, b, and c). Also, it is important to note that these indices were calculated for the extent of the observed debris flows and may not represent the topographic form that controlled the model.
In summary, using the calibration utility, we showed how MWR can be calibrated to a range of different landslide types and runout terrains. To a certain degree, though calibration, MWR can be parameterized to compensate for deficiencies in the DEM or processes not explicitly represented in the model (momentum, woody debris). While model performance at the Rocky Mountains site suggests that MWR may not perform as well on planar hillslopes, a relationship between model performance and topography was not eminent. This finding is likely a result of the contributions of numerous factors other than the terrain form, such as the DEM resolution, the quality of the DoD, the accuracy of the regolith map, and the importance of processes not explicitly included in the model that also impact performance.
6.1 Strategic testing of MWR for hazard mapping applications
Having demonstrated calibrated model performance for a variety of landslides and runout terrains, we now strategically test MWR using the Cascade Mountains and Black Hills sites. Since both of these sites include two separate landslides, we can thus test model performance by swapping bestfit model parameters at each site, rerunning the models, and comparing results with the original calibrated results. At the Cascade Mountains site, the 2009 and 2022 landslides originated on the same hillslope (Fig. 8a and b). At the Black Hills site, the two landslides occurred on different hillslopes but in adjacent east–westoriented watersheds (Fig. 8c and d).
As shown in Fig. 11, at three of the landslides (both the Cascade Mountains landslides and the Black Hills north landslide), when the bestfit parameters from the other landslide are used to predict runout, the accuracy of modeled runout planimetric extent drops but resultant Ω_{T} values can still be as high as or higher than values reported in other studies (compare to equivalent Ω_{T} values in Gorr et al., 2022, and Barnhart et al., 2021). In terms of modeled sediment transport and topographic change, swapping bestfit parameters has a more substantial effect. For the Cascade Mountains 2009 landslide, using the 2022 bestfit parameter values causes about half of the modeled runout material to be prematurely deposited on the hillslope, reducing the amount of sediment that reaches the valley floor (${Q}_{{\mathrm{s}}_{\mathrm{E}}}$ increases by a factor of 9; Fig. 11). Using the Cascade Mountains 2009 parameter values on the Cascade Mountains 2022 landslide (Fig. 11b) increases modeled runout extent and results in nearly 4 times the entrainment and transport of sediment to the valley floor, causing ${Q}_{{\mathrm{s}}_{\mathrm{E}}}$ to increase by a factor of 20 and Δη_{E} by 83 %. At the Black Hills site, using the south basin bestfit model parameters at the north basin causes ${Q}_{{\mathrm{s}}_{\mathrm{E}}}$ and Δη_{E} to increase by 83 % and 39 %, respectively (Fig. 11c). Unlike the other three landslides, swapping bestfit parameters at the Black Hills south landslide results in both large sediment transport and runout extent error because the north landslide bestfit parameters applied to the south landslide causes the model to entrain too little and stop prematurely (Fig. 11d).
As landslide hazard models often forecast hazard probabilistically, an alternative test to simply swapping the bestfit parameters is to swap parameter PDFs determined from the calibration utility and compare the probability of runout at each model node (Eq. 29). As shown in Fig. 12, similar to the first test, for three of the landslides, using the parameter distribution associated with the neighboring landslide results in relatively minor changes in model output, as indicated by where runout is more likely to occur versus not occur (probability of runout ≥50 %; Fig. 12a, b, and d), but at the Black Hills south landslide, swapping parameter PDFs causes a large change in runout probability (Fig. 12c).
The results of these two tests suggest that sitespecific or even landslidetypespecific calibration may be needed to accurately predict runout behavior using MWR, especially when the user aims to apply MWR to sediment yield analyses and the study site consists of numerous landslide types, like the Black Hills and Cascade Mountain sites. While the need for calibration of MWR may limit model transferability across sites, this limitation holds true for most physicsbased models. Barnhart et al. (2021) compared the ability of three different detailed mechanistic models to replicate an observed postwildfire debris flow runout event in California, USA. All three models used a shallowwaterequationbased approach that conserved both mass and momentum, representing the flow as either a singlephase or doublephase fluid. All models gave comparable results in simulating the event, suggesting that there may not be a “true” best model. Despite the high level of detail and processes explicitly included in each model, all models were sensitive to and required an estimate of the total mobilized volume, and the ability to replicate observed runout ultimately depended on calibration.
We suspect that in regions where landslide processes are relatively uniform (like the Olympic Mountains site), parameterization determined through calibration to one landslide might be transferable across sites. Additionally, as noted in Sect. 3.1, we found that numerous parameter combinations allowed MWR to match observed runout extent. This finding suggests that if the project aim is limited to an evaluation of runout extent, model calibration to the site may not be as critical and parameter values from calibration to nearby landslides or even globally available repeated DEMs and air photos that show the slope of past landslide deposits (for S_{c}) and deposit thickness (for q_{c}) might be sufficient.
6.2 MassWastingRunout probability applications
In this section we briefly demonstrate how to determine runout probability from a probabilistically determined landslide hazard map or a specific, potentially unstable slope using MWR. The first application may be appropriate for watershed to regionalscale runout hazard assessments. The second application is an example hazard assessment for a potentially unstable hillslope. Both applications are demonstrated at the Olympic Mountains site where landslide size and type tended to be relatively uniform and parameter PDFs determined through calibration may therefore represent typical runout processes in the basin.
6.2.1 Runout probability from a landslide hazard map
To determine runout probability from a landslide hazard map, we ran MWR Probability using Option 3, reading a series of mapped landslide source areas created by an externally run Monte Carlo landslide initiation model. For the landslide initiation model, we used LandslideProbability (Strauch et al., 2018), an existing component in Landlab that computes landslide probability by iteratively calculating the factor of safety (FS: ratio of the resisting to the driving forces) at each node on the raster model grid Np times from randomly selected soil (regolith) hydrology properties (e.g., soil depth, saturated hydraulic conductivity), soil strength (friction angle, cohesion), and recharge rates (precipitation input rate minus evapotranspiration and soil storage). Landslide probability at a node is defined as the number of times FS<1 divided by Np.
We first ran LandslideProbability using a 50year precipitation event (WRCC, 2024) to determine landslide probability (Fig. 13a) over the entire Olympic Mountains model domain and create the series of Np FS maps. Details on the LandslideProbability setup are included in the Supplement. We then read the series of Np FS maps into MWR Probability, treating all nodes with FS<1 as a landslide source, and ran MWR Np times. Each iteration, MWR read a new FS map and randomly selected a new set of parameter values from S_{c}–q_{c} parameter PDFs created by the calibration utility.
Runout probability, which reflects MWR parameter uncertainty (as illustrated in Fig. 7f) and uncertainty in the initial landslide size and location caused by a 50year precipitation event, is illustrated in Fig. 13b and shows that the probability of runout is high in many of the secondorder channels but low at the basin outlet. As discussed in Sect. 3, the probability of aggradation or erosion caused by the runout can be determined by adjusting the numerator of Eq. (29). As an example, the probability of deposition greater than 1 m is shown in Fig. 13c.
6.2.2 Runout probability for a specific, potentially unstable slope
When field evidence or other data indicate that a specific hillslope may be potentially unstable but the exact area of a potential landslide on that slope is unknown, MWR can be used to generate a hazard estimate that takes into account the uncertainty in the landslide area. For this application, MWR Probability is run using Option 2, which requires a polygon representing the extent of the potentially unstable slope. We designated a 0.6 ha, convergent hillslope in the headwaters of the Olympic Mountains site as a potentially unstable slope (Fig. 13d). For each model repetition, a landslide area can form anywhere within the potentially unstable slope and is at least as large as a userdefined minimum size but no larger than the potentially unstable slope. This example shows that, given uncertainty in the landslide size and location and uncertainty in MWR parameterization, if a landslide were to initiate on the potentially unstable slope, the probability of the runout reaching the basin outlet is less than 5 %.
In this study, we described, calibrated, and tested MassWastingRunout (MWR), a new cellularautomaton landslide runout model that combines the functionality of simple runout algorithms used in WSMs (landscape evolution and watershed sediment yield models) with the predictive detail typical of runout models used for landslide inundation hazard mapping. MWR is implemented in Python as a component for the Landlab earth surface modeling toolkit and is designed for WSM and probabilistic landslide hazard assessment applications. MWR includes a Markov chain Monte Carlo calibration utility that determines the bestfit parameter values for a site as well as empirical probability density functions (PDFs) of the parameter values. MWR also includes a utility called MWR Probability that takes the PDF output from the calibration utility to determine runout probability.
Given a DEM and map of approximate regolith depth, MWR needs only the location and geometry of an initial landslide source area to model the entire runout process. Results indicate that despite its simple conceptualization, MWR shows skill in modeling the final runout extent, sediment transport, and topographic change associated with a landslide. When compared to other models capable of replicating observed landslide inundation patterns, the strength of MWR lies in its use of fieldinferable parameters, its ability to internally estimate the total mobilized volume (initial landslide body + erosion volume), and its relatively parsimonious model design.
MWR can be calibrated to a site using just two parameters (critical slope, S_{c}, and a threshold flux for deposition, q_{c}) and the MWR calibration utility enables the user to calibrate the model for a watershed within several hours on a standard desktop (Sect. 5.3). Although the predictive power of MWR hinges on calibration – a common requirement for mechanistic models – its reliance on two calibration parameters serves to constrain model uncertainty. Sitespecific calibration may be needed when MWR is used for sediment yield analysis, but if the aim is limited to mapping runout extent, it may be possible to infer parameterization from nearby landslides or possibly from globally available repeated DEMs and air photos that show where past masswasting flows have stopped (for S_{c}) and how thick their frontal lobes are at the point of deposition (for q_{c}). Nonetheless, as a rulebased, cellularautomaton model, MWR is not designed to accurately simulate flow depth. For accurate flow depths or debris flow impact forces, a detailed mechanistic modeling approach should be used.
MWR shows a rich set of intuitive responses to topographic curvature and slope (Sect. 4). When calibrated to the runout of six different observed landslides, the volumetric error of MWR, Δη_{E}, ranged between 6 % and 15 % (median 9.1 %) of the observed total mobilized volume. Except for the Rocky Mountains site where MWR consistently modeled widerthanobserved flow, the cumulative flow error along the runout profile (${Q}_{{\mathrm{s}}_{\mathrm{E}}}$) was limited to 5 %–19 % of the mean cumulative flow determined from the observed DEM of difference (DoD). These are considered acceptable levels of performance given that the total mobilized volume of many debris flow models assumes an orderofmagnitude range of confidence. A notable finding of this paper is that at most sites, MWRmodeled runout did not have any strong systematic bias in predictions (toward unrealistically short or wide flows, for example), which suggests that MWR is structurally sound. However, MWR may underperform compared to mechanistic models when flow momentum is the primary driver of runout extent (e.g., in areas of slopeperpendicular flow).
Finally, we briefly showed how to couple MWR with the landslide initiation model LandslideProbability to map debris flow hazard when the initial landslide location is uncertain. As a component of the Landlab earth surface modeling toolkit, MWR is designed to be compatible with other models and thus relatively easy to integrate into a WSM. An example WSM that incorporates MWR might include models for landslide initiation, hillslope diffusion, and fluvial incision to investigate the role of landslides and their runout in longterm landscape evolution. Future studies will explore largescale application in landscape evolution or sediment yield models and characterize model parameters for different landslide types and hydroclimatic conditions. The use of a calibrated runout model in WSMs might allow for regionspecific and more insightful predictions of landslide impact on landscape morphology and watershedscale sediment dynamics.
The average erosion depth caused by the observed runout ($\stackrel{\mathrm{\u203e}}{E}$) can be determined from the DoD as the total erosion volume (∑EΔx^{2}) divided by the erosion area (A_{e}):
where ∑EΔx^{2} and A_{e} exclude the initial landslide body volume and area, areas of deposition (Δη>0), and areas with no change in elevation (Δη=0). In terms of the debriton conceptualization used in MWR, $\stackrel{\mathrm{\u203e}}{E}$ can also be written as a function of the mean number of times a debriton would need to pass over a grid cell ($\stackrel{\mathrm{\u203e}}{n}$) multiplied by an average erosion depth per debriton (${\stackrel{\mathrm{\u203e}}{h}}_{\mathrm{e}}$) to equal $\stackrel{\mathrm{\u203e}}{E}$ as
An estimate for $\stackrel{\mathrm{\u203e}}{n}$ can be determined from the average length of the runout material divided by the cell width:
At most sites, we approximate the average length of the runout material simply as the mapped landslide length (l). As the debritons move down slopes in excess of S_{c}, they entrain material, split, and spread, and the runout material tends to lengthen. Using the initial landslide length to represent the runout length thus represents a minimum value for $\stackrel{\mathrm{\u203e}}{n}$, and if needed, the numerator of Eq. (A3) can be multiplied by a coefficient to better match runout length. Combining Eqs. (A2) and (A3), ${\stackrel{\mathrm{\u203e}}{h}}_{\mathrm{e}}$ can be defined as the average erosion rate per unit length of runout debris ($\stackrel{\mathrm{\u203e}}{E}/l$) times the cell width:
Rewriting Eq. (11) as a function of the average shear stress in the erosiondominated reaches of the runout path ($\stackrel{\mathrm{\u203e}}{\mathit{\tau}}$) and assuming τ_{c}≅0, the debris flow erodibility parameter k can be estimated as
To solve for k, we estimated $\stackrel{\mathrm{\u203e}}{\mathit{\tau}}$ from fieldapproximated debris flow depth and channel slope measurements in the erosiondominated reaches of the runout path (Table 1). We used Eq. (13) to define $\stackrel{\mathrm{\u203e}}{\mathit{\tau}}$. For D_{s}, we used the average maximum grain size observed over the whole runout path. If τ is defined as a function of the graincollisiondependent shear stress approach in Eq. (13) and k is determined as a function of f, as in Eq. (A5), the impact of f on model behavior is relatively small.
MassWastingRunout and several tutorial notebooks area available at https://github.com/landlab/landlab (Landlab, 2024).
Validation field site data, including the landslide source area and DEM of difference, are available at http://www.hydroshare.org/resource/55813a5e01764546b76641a7385c2236 (Keck, 2024).
Jupyter Notebook tutorials that allow the user to recreate and modify the plots in Fig. 5 or model runout at the Cascade Mountains 2022 site can be run online or downloaded to a personal computer. The online tutorial and directions for downloading to a personal computer are available at the Landlab website: https://landlab.github.io/#/ (year of publication: 2024, Landlab, 2017).
The supplement related to this article is available online at: https://doi.org/10.5194/esurf1211652024supplement.
JK conceptualized and coded MWR and the calibration and probability utilities, hiked and surveyed each validation site, prepared the validation site model input data, calibrated and assessed model performance, wrote the manuscript, and prepared the figures. EI served as an advisor and helped with model conceptualization, writing, figure formatting, and funding. BC helped with field surveys, model conceptualization, and writing. GT helped with model conceptualization and writing. AHD helped with figure formatting and model conceptualization.
The contact author has declared that none of the authors has any competing interests.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.
This research was partially supported by the following programs: National Science Foundation (NSF) PREEVENTS program (ICER1663859), NSF OAC2103632, and NASA Disasters Program (grant number 80NSSC23K1103). Critical and helpful referee reviews as well as multiple detailed reviews by associate editor Wolfgang Schwanghart significantly improved the manuscript. Early conceptualization of the model as well as development of the ${Q}_{{\mathrm{s}}_{\mathrm{E}}}$ metric greatly benefited from discussion with Hervé Capart. Discussions with TzuYin Kasha Chen and ChiYao Hung led to model improvements. Stephen Slaughter fieldreviewed the Cascade Mountains 2009 landslide and the Black Hills landslides the year they occurred and provided photos and field observations that aided author interpretation. John Jenkins helped with field reconnaissance, and Miles Micheletti prepared the structurefrommotion DEM at the Cascade Mountains 2022 landslide. Eli Schwat helped with field reconnaissance at the Olympic Mountains site. This work also benefitted from coding guidance from Eric Hutton and support from the staff and researchers at CSDMS.
This research has been supported by the National Science Foundation (NSF) PREEVENTS program (grant nos. ICER1663859 and OAC2103632) and the NASA Disasters Program (grant no. 80NSSC23K1103).
This paper was edited by Wolfgang Schwanghart and reviewed by Saskia de Vilder and two anonymous referees.
Bagnold, R. A.: Experiments on a gravityfree dispersion of large solid spheres in a Newtonian fluid under shear, P. R. Soc. Lond., 225, 49–63, https://doi.org/10.1098/rspa.1954.0186, 1954.
Barca, D., Crisci, G., Di Gregorio, S., and Nicoletta, F.: Cellular automata method for modelling lava flows: Simulation of the 1986–1987 eruption, Mount Etna, Sicily, in Active lavas: Monitoring and modeling, edited by: Kilburn, C., and Luongo, G., University College of London Press, London, 291–309, https://doi.org/10.4324/978100332517815, 1993.
Barnhart, K. R., Hutton, E. W. H., Tucker, G. E., Gasparini, N. M., Istanbulluoglu, E., Hobley, D. E. J., Lyons, N. J., Mouchene, M., Nudurupati, S. S., Adams, J. M., and Bandaragoda, C.: Short communication: Landlab v2.0: a software package for Earth surface dynamics, Earth Surf. Dynam., 8, 379–397, https://doi.org/10.5194/esurf83792020, 2020.
Barnhart, K. R., Jones, R., George, D. J., McArdell, B. W., Rengers, F. K., Staley, D. M., and Kean, J. W.: MultiModel Comparison of Computed Debris Flow Runout for the 9 January 2018 Montecito, California PostWildfire Event, J. Geophys. Res.Earth, 126, e2021JF006245, https://doi.org/10.1029/2021jf006245, 2021.
Benda, L. and Dunne, T.: Stochastic forcing of sediment supply to channel networks from landsliding and debris flow, Water Resour. Res., 33, 2849–2863, https://doi.org/10.1029/97wr02388, 1997.
Benda, L., Veldhuisen, C. P., and Black, J.: Debris flows as agents of morphological heterogeneity at loworder confluences, Olympic Mountains, Washington, Geol. Soc. Am. Bull., 115, 1110, https://doi.org/10.1130/b25265.1, 2003.
Beven, K.: A manifesto for the equifinality thesis, J. Hydrol., 320, 18–36, https://doi.org/10.1016/j.jhydrol.2005.07.007, 2006.
Burton, A. and Bathurst, J. C.: Physically based modelling of shallow landslide sediment yield at a catchment scale, Environ. Geol., 35, 89–99, https://doi.org/10.1007/s002540050296, 1998.
Bigelow, P., Benda, L., Miller, D., and Burnett, K. M.: On Debris Flows, River Networks, and the Spatial Structure of Channel Morphology, Forest Sci., 53, 220–238, https://doi.org/10.1093/forestscience/53.2.220, 2007.
Campforts, B., Shobe, C. M., Overeem, I., and Tucker, G. E.: The Art of Landslides: How Stochastic Mass Wasting Shapes Topography and Influences Landscape Dynamics, J. Geophys. Res.Earth, 127, e2022JF006745, https://doi.org/10.1029/2022jf006745, 2022.
Campforts, B., Shobe, C. M., Steer, P., Vanmaercke, M., Lague, D., and Braun, J.: HyLands 1.0: a hybrid landscape evolution model to simulate the impact of landslides and landslidederived sediment on landscape evolution, Geosci. Model Dev., 13, 3863–3886, https://doi.org/10.5194/gmd1338632020, 2020.
Capart, H. and Fraccarollo, L.: Transport layer structure in intense bedload, Geophys. Res. Lett., 38, L20402, https://doi.org/10.1029/2011gl049408, 2011.
Capart, H., Hung, C., and Stark, C. R.: Depthintegrated equations for entraining granular flows in narrow channels, J. Fluid Mech., 765, R4, https://doi.org/10.1017/jfm.2014.713, 2015.
Carretier, S., Martinod, P., Reich, M., and Godderis, Y.: Modelling sediment clasts transport during landscape evolution, Earth Surf. Dynam., 4, 237–251, https://doi.org/10.5194/esurf42372016, 2016.
Chase, C. G.: Fluvial landsculpting and the fractal dimension of topography, Geomorphology, 5, 39–57, https://doi.org/10.1016/0169555x(92)90057u, 1992.
Chen, C. and Yu, F.: Morphometric analysis of debris flows and their source areas using GIS, Geomorphology, 129, 387–397, https://doi.org/10.1016/j.geomorph.2011.03.002, 2011.
Chen, H. X. and Zhang, L. M.: EDDA 1.0: integrated simulation of debris flow erosion, deposition and property changes, Geosci. Model Dev., 8, 829–844, https://doi.org/10.5194/gmd88292015, 2015.
Chen, T.Y. K., Wu, Y.C., Hung, C.Y., Capart, H., and Voller, V. R.: A control volume finiteelement model for predicting the morphology of cohesivefrictional debris flow deposits, Earth Surf. Dynam., 11, 325–342, https://doi.org/10.5194/esurf113252023, 2023.
Clerici, A. and Perego, S.: Simulation of the Parma River blockage by the Corniglio landslide (Northern Italy), Geomorphology, 33, 1–23, https://doi.org/10.1016/s0169555x(99)000951, 2000.
Codd, E. F.: Cellular Automata (1st ed.), Academic Press, New York, 122 pp., ISBN 0121788504, 1968.
Coz, J. L., Renard, B., Bonnifait, L., Branger, F., and Boursicaud, R. L.: Combining hydraulic knowledge and uncertain gaugings in the estimation of hydrometric rating curves: A Bayesian approach, J. Hydrol., 509, 573–587, https://doi.org/10.1016/j.jhydrol.2013.11.016, 2014.
Crave, A. and Davy, P.: A stochastic “precipiton” model for simulating erosion/sedimentation dynamics, Comput. Geosci., 27, 815–827, https://doi.org/10.1016/s00983004(00)001679, 2001.
D'Ambrosio, D., Di Gregorio, S., Iovine, G., Lupiano, V., Rongo, R., and Spataro, W.: First simulations of the Sarno debris flows through Cellular Automata modelling, Geomorphology, 54, 91–117, https://doi.org/10.1016/s0169555x(03)000588, 2003.
Doten, C. O., Bowling, L. C., Lanini, J. S., Maurer, E. P., and Lettenmaier, D. P.: A spatially distributed model for the dynamic prediction of sediment erosion and transport in mountainous forested watersheds, Water Resour. Res., 42, W04417, https://doi.org/10.1029/2004wr003829, 2006
Egashira, S., Honda, N., and Itoh, T.: Experimental study on the entrainment of bed material into debris flow, Phys. Chem. Earth, Parts A/B/C, 26, 645–650, https://doi.org/10.1016/s14641917(01)000629, 2001.
Fannin, R. J. and Wise, M. P.: An empiricalstatistical model for debris flow travel distance, Can. Geotech. J., 38, 982–994, https://doi.org/10.1139/t01030, 2001.
Frank, F., McArdell, B. W., Huggel, C., and Vieli, A.: The importance of entrainment and bulking on debris flow runout modeling: examples from the Swiss Alps, Nat. Hazards Earth Syst. Sci., 15, 2569–2583, https://doi.org/10.5194/nhess1525692015, 2015.
Freeman, T. G.: Calculating catchment area with divergent flow based on a regular grid, Comput. Geosci., 17, 413–422, https://doi.org/10.1016/00983004(91)90048i, 1991.
Gartner, J. E., Cannon, S. H., and Santi, P. M.: Empirical models for predicting volumes of sediment deposited by debris flows and sedimentladen floods in the transverse ranges of southern California, Eng. Geol., 176, 45–56, https://doi.org/10.1016/j.enggeo.2014.04.008, 2014.
Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari, A., and Rubin, D. B.: Bayesian Data Analysis (3rd ed.), Electronic Edition, http://www.stat.columbia.edu/~gelman/book/BDA3.pdf (last access: 17 April 2023), 2021.
Goode, J. R., Luce, C. H., and Buffington, J. M.: Enhanced sediment delivery in a changing climate in semiarid mountain basins: Implications for water resource management and aquatic habitat in the northern Rocky Mountains, Geomorphology, 139–140, 1–15, https://doi.org/10.1016/j.geomorph.2011.06.021, 2012.
Gorr, A., McGuire, L. A., Youberg, A., and Rengers, F. K.: A progressive flowrouting model for rapid assessment of debrisflow inundation, Landslides, 19, 2055–2073, https://doi.org/10.1007/s1034602201890y, 2022.
Guthrie, R., Hockin, A., Colquhoun, L., Nagy, T., Evans, S. G., and Ayles, C. P.: An examination of controls on debris flow mobility: Evidence from coastal British Columbia, Geomorphology, 114, 601–613, https://doi.org/10.1016/j.geomorph.2009.09.021, 2010.
Guthrie, R. and Befus, A.: DebrisFlow Predictor: an agentbased runout program for shallow landslides, Nat. Hazards Earth Syst. Sci., 21, 1029–1049, https://doi.org/10.5194/nhess2110292021, 2021.
Hammond, C. J., Prellwitz, R. W., and Miller, S. M.: Landslides hazard assessment using Monte Carlo simulation. in: Proceedings of 6th international symposium on landslides, edited by: Bell, D. H., Christchurch, New Zealand, Balkema, 2, 251–294, ISBN 9789054100324, 1992.
Han, Z., Chen, G., Li, Y., Tang, C., Xu, L., He, Y., Huang, X., and Wang, W.: Numerical simulation of debrisflow behavior incorporating a dynamic method for estimating the entrainment, Eng. Geol., 190, 52–64, https://doi.org/10.1016/j.enggeo.2015.02.009, 2015.
Han, Z., Li, Y., Huang, J., Chen, G., Xu, L., Tang, C. Y., Zhang, H., and Shang, Y.: Numerical simulation for runout extent of debris flows using an improved cellular automaton model, B. Eng. Geol. Environ., 76, 961–974, https://doi.org/10.1007/s1006401609026, 2017.
Han, Z., Ma, Y., Li, Y., Zhang, H., Chen, N., Hu, G., and Chen, G.: Hydrodynamic and topography based cellular automaton model for simulating debris flow runout extent and entrainment behavior, Water Res., 193, 116872, https://doi.org/10.1016/j.watres.2021.116872, 2021.
Heiser, M., Scheidl, C., and Kaitna, R.: Evaluation concepts to compare observed and simulated deposition areas of mass movements, Comput. Geosci., 21, 335–343, https://doi.org/10.1007/s1059601696099, 2017.
Hobley, D. E. J., Adams, J. M., Nudurupati, S. S., Hutton, E. W. H., Gasparini, N. M., Istanbulluoglu, E., and Tucker, G. E.: Creative computing with Landlab: an opensource toolkit for building, coupling, and exploring twodimensional numerical models of Earthsurface dynamics, Earth Surf. Dynam., 5, 21–46, https://doi.org/10.5194/esurf5212017, 2017.
Holmgren, P.: Multiple flow direction algorithms for runoff modelling in grid based elevation models: An empirical evaluation, Hydrol. Process., 8, 327–334, https://doi.org/10.1002/hyp.3360080405, 1994.
Horton, P., Jaboyedoff, M., Rudaz, B., and Zimmermann, M.: FlowR, a model for susceptibility mapping of debris flows and other gravitational hazards at a regional scale, Nat. Hazards Earth Syst. Sci., 13, 869–885, https://doi.org/10.5194/nhess138692013, 2013.
Hungr, O., Morgan, G. J., and Kellerhals, R.: Quantitative analysis of debris torrent hazards for design of remedial measures, Can. Geotech. J., 21, 663–677, https://doi.org/10.1139/t84073, 1984.
Hungr, O. and Evans, S. G.: Entrainment of debris in rock avalanches: An analysis of a long runout mechanism, Geol. Soc. Am. Bull., 116, 1240, https://doi.org/10.1130/b25362.1, 2004.
Hutter, K., Svendsen, B., and Rickenmann, D.: Debris flow modeling: A review, Continuum Mech. Therm., 8, 1–35, https://doi.org/10.1007/bf01175749, 1996.
Iovine, G., D'Ambrosio, D., and Di Gregorio, S.: Applying genetic algorithms for calibrating a hexagonal cellular automata model for the simulation of debris flows characterised by strong inertial effects, Geomorphology, 66, 287–303, https://doi.org/10.1016/j.geomorph.2004.09.017, 2005.
Istanbulluoglu, E. Bras R. L.: Vegetationmodulated landscape evolution: Effects of vegetation on landscape processes, drainage density, and topography, J. Geophys. Res., 110, F02012, https://doi.org/10.1029/2004jf000249, 2005.
Istanbulluoglu, E., Bras R. L., FloresCervantes, H., and Tucker, G. E.: Implications of bank failures and fluvial erosion for gully development: Field observations and modeling, J. Geophys. Res., 110, F01014, https://doi.org/10.1029/2004JF000145, 2005.
Istanbulluoglu, E., O. Yetemen, E. R. Vivoni, H. A. Gutie' rrezJurado, and R. L. Bras, Ecogeomorphic implications of hillslope aspect: Inferences from analysis of landscape morphology in central New Mexico, Geophys. Res. Lett., 35, L14403, https://doi.org/10.1029/2008GL034477, 2008.
Iverson, R. M.: The physics of debris flows, Rev. Geophys., 35, 245–296, https://doi.org/10.1029/97rg00426, 1997.
Iverson, R. M.: How should mathematical models of geomorphic processes be judged?, in: Prediction in Geomorphology, edited by: Wilcock, P. and Iverson, R., American Geophysical Union, ISBN 0875909930, 2003.
Iverson, R. M.: Elementary theory of bedsediment entrainment by debris flows and avalanches, J. Geophys. Res., 117, F03006, https://doi.org/10.1029/2011jf002189, 2012.
Iverson, R. M. and Denlinger, R. P.: Flow of variably fluidized granular masses across threedimensional terrain: 1. Coulomb mixture theory, J. Geophys. Res., 106, 537–552, https://doi.org/10.1029/2000jb900329, 2001.
Julien, P. Y. and Paris, A.: Mean Velocity of Mudflows and Debris Flows, J. Hydraul. Eng., 136, 676–679, https://doi.org/10.1061/(asce)hy.19437900.0000224, 2010.
Kean, J. W., Staley, D. M., Lancaster, J., Rengers, F., Swanson, B., and Coe, J.: Inundation, flow dynamics, and damage in the 9 January 2018 Montecito debrisflow event, California, USA: Opportunities and challenges for postwildfire risk assessment, Geosphere, 15, 1140–1163, https://doi.org/10.1130/GES02048.1, 2019.
Keck, J.: MassWastingRunout field site data, HydroShare [data set], http://www.hydroshare.org/resource/55813a5e01764546b76641a7385c2236 (last access: 17 September 2024), 2024.
Korup, O.: Effects of large deepseated landslides on hillslope morphology, western Southern Alps, New Zealand, J. Geophys. Res., 111, F01018, https://doi.org/10.1029/2004jf000242, 2006.
Lancaster, S. T., Hayes, S. K., and Grant, G. E.: Effects of wood on debris flow runout in small mountain watersheds, Water Resour. Res., 39, 1168, https://doi.org/10.1029/2001wr001227, 2003.
landlab: landlab, GitHub [code], https://github.com/landlab/landlab, 15 August 2024.
Landlab: Landlab, https://landlab.github.io/#/ (last access: 17 September 2024), 2017.
Larsen, I. J. and Montgomery, D. R.: Landslide erosion coupled to tectonics and river incision, Nat. Geosci., 5, 468–473, https://doi.org/10.1038/ngeo1479, 2012.
Liu, J., Wu, Y., Gao, X., and Zhang, X.: A Simple Method of Mapping Landslides Runout Zones Considering Kinematic Uncertainties, Remote Sens.Basel, 14, 668, https://doi.org/10.3390/rs14030668, 2022.
Major, J. J.: Depositional Processes in LargeScale DebrisFlow Experiments, J. Geol., 105, 345–366, https://doi.org/10.1086/515930, 1997.
Major, J. J. and Iverson, R. M.: Debrisflow deposition: Effects of porefluid pressure and friction concentrated at flow margins, Geol. Soc. Am. Bull., 111, 1424–1434, https://doi.org/10.1130/00167606(1999)111<1424:DFDEOP>2.3.CO;2, 1999.
McCoy, S. W., Kean, J. W., Coe, J. A., Tucker, G. S., Staley, D. M., and Wasklewicz, T. A.: Sediment entrainment by debris flows: In situ measurements from the headwaters of a steep catchment, J. Geophys. Res., 117, F03016, https://doi.org/10.1029/2011jf002278, 2012.
McDougall, S. and Hungr, O.: A model for the analysis of rapid landslide motion across threedimensional terrain, Can. Geotech. J., 41, 1084–1097, https://doi.org/10.1139/t04052, 2004.
Medina, V., Hürlimann, M., and Bateman, A.: Application of FLATModel, a 2D finite volume code, to debris flows in the northeastern part of the Iberian Peninsula, Landslides, 5, 127–142, https://doi.org/10.1007/s1034600701023, 2008.
Montgomery, D. R. and Dietrich, W. E.: Where do channels begin?, Nature, 336, 232–234, https://doi.org/10.1038/336232a0, 1988.
Murray, B. A. and Paola, C.: A cellular model of braided rivers, Nature, 371, 54–57, https://doi.org/10.1038/371054a0, 1994.
Murray, A. B. and Paola, C.: Properties of a cellular braidedstream model, Earth Surf. Proc. Land., 22, 1001–1025, https://doi.org/10.1002/(SICI)10969837(199711)22:11<1001::AIDESP798>3.0.CO;2O, 1997.
Murray, A. B.: Reducing model complexity for explanation and prediction, Geomorphology, 90, 178–191, https://doi.org/10.1016/j.geomorph.2006.10.020, 2007.
Murray A. B.: Which Models Are Good (Enough), and When?. in: Treatise on Geomorphology, edited by: John F. Shroder, Volume 2, 50–58, San Diego, Academic Press, ISBN 9780123983534, 2013.
Natural Resources Conservation Service (NRCS): Snow and Water Interactive Map (n.d.), Natural Resources Conservation Service [data set], https://www.nrcs.usda.gov/resources/dataandreports/snowandwaterinteractivemap (last access: 15 April 2022), 2022.
Nudurupati, S. S., Istanbulluoglu, E., Tucker, G. E., Gasparini, N. M., Hobley, D. E. J., Hutton, E. W. H., Barnhart, K. R., and Adams, J. M.: On transient semiarid ecosystem dynamics using Landlab: vegetation shifts, topographic refugia, and response to climate, Water Resour. Res., 59, e2021WR031179, https://doi.org/10.1029/2021wr031179, 2023.
Perron, J. T.: Climate and the Pace of Erosional Landscape Evolution, Annu. Rev. Earth Pl. Sc., 45, 561–591, https://doi.org/10.1146/annurevearth060614105405, 2017.
Reid, M. J., Coe, J. A., and Brien, D. L.: Forecasting inundation from debris flows that grow volumetrically during travel, with application to the Oregon Coast Range, USA, Geomorphology, 273, 396–411, https://doi.org/10.1016/j.geomorph.2016.07.039, 2016.
Renard, B., Garreta, V., and Lang, M. J.: An application of Bayesian analysis and Markov chain Monte Carlo methods to the estimation of a regional trend in annual maxima, Water Resour. Res., 42, W12422, https://doi.org/10.1029/2005wr004591, 2006.
Rengers, F. K., McGuire, L. A., Kean, J. W., Staley, D. M., and Hobley, D. E. J.: Model simulations of flood and debris flow timing in steep catchments after wildfire, Water Resour. Res., 52, 6041–6061, https://doi.org/10.1002/2015WR018176, 2016.
RodaBoluda, D. C., D'Arcy, M., McDonald, J., and Whittaker, A. C.: Lithological controls on hillslope sediment supply: insights from landslide activity and grain size distributions, Earth Surf. Proc. Land., 43, 956–977, https://doi.org/10.1002/esp.4281, 2018.
Shaller, P. J., Doroudian, M., and Hart, M. W.: The Eureka Valley Landslide: Evidence of a dual failure mechanism for a LongRunout Landslide, Lithosphere, 2020, 8860819, https://doi.org/10.2113/2020/8860819, 2020.
Schürch, P., Densmore, A. L., Rosser, N., and McArdell, B. W.: Dynamic controls on erosion and deposition on debrisflow fans, Geology, 39, 827–830, https://doi.org/10.1130/g32103.1, 2011.
Shen, P., Zhang, L. M., Wong, H., Peng, D., Zhou, S., Zhang, S., and Chen, C.: Debris flow enlargement from entrainment: A case study for comparison of three entrainment models, Eng. Geol., 270, 105581, https://doi.org/10.1016/j.enggeo.2020.105581, 2020.
Stock, J. P. J. and Dietrich, W. E.: Erosion of steepland valleys by debris flows, Geol. Soc. Am. Bull., 118, 1125–1148, https://doi.org/10.1130/b25902.1, 2006.
Strauch, R., Istanbulluoglu, E., Nudurupati, S. S., Bandaragoda, C., Gasparini, N. M., and Tucker, G. E.: A hydroclimatological approach to predicting regional landslide probability using Landlab, Earth Surf. Dynam., 6, 49–75, https://doi.org/10.5194/esurf6492018, 2018.
Takahashi, T.: Debris Flow (2nd ed.), CRC Press, Taylor and Francis Group, ISBN 9780203576748, 2014.
Tucker, G. E. and Bras, R. L.: Hillslope processes, drainage density, and landscape morphology, Water Resour. Res., 34, 2751–2764, https://doi.org/10.1029/98wr01474, 1998.
Tucker, G. E. and Hancock, G. J.: Modelling landscape evolution, Earth Surf. Proc. Land., 35, 28–50, https://doi.org/10.1002/esp.1952, 2010.
Tucker, G. E., McCoy, S. W., and Hobley, D. E. J.: A lattice grain model of hillslope evolution, Earth Surf. Dynam., 6, 563–582, https://doi.org/10.5194/esurf65632018, 2018.
Western Regional Climate Center (WRCC): RAWS USA Climate Archive, https://wrcc.dri.edu/, last access: 17 September 2024.
Whipple, K. X. and Dunne, T.: The influence of debrisflow rheology on fan morphology, Owens Valley, California, Geol. Soc. Am. Bull., 104, 887–900, https://doi.org/10.1130/00167606(1992)104<0887:TIODFR>2.3.CO;2, 1992.
Zhou, G. G. D., Li, S., Song, D., Choi, C. E., and Chen, X.: Depositional mechanisms and morphology of debris flow: physical modelling, Landslides, 16, 315–332, https://doi.org/10.1007/s1034601810959, 2019.
 Abstract
 Introduction
 Description of the MassWastingRunout model
 Calibration and MWR Probability
 Basic model behavior
 Model validation
 Discussion
 Concluding remarks
 Appendix A: Determination of k
 Appendix B
 Code availability
 Data availability
 Executable research compendium (ERC)
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References
 Supplement
 Abstract
 Introduction
 Description of the MassWastingRunout model
 Calibration and MWR Probability
 Basic model behavior
 Model validation
 Discussion
 Concluding remarks
 Appendix A: Determination of k
 Appendix B
 Code availability
 Data availability
 Executable research compendium (ERC)
 Author contributions
 Competing interests
 Disclaimer
 Acknowledgements
 Financial support
 Review statement
 References
 Supplement