Luke A Mcguire
email: lmcguire@arizona.edu

Scott W Mccoy

Odin Marc

William Struble

Katherine R Barnhart

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

Debris flows regularly traverse bedrock channels that dissect steep landscapes, but our understanding of the dominant controls that set rates and spatial patterns of landscape evolution by debris flows is still rudimentary. As a result, much of our understanding of steep bedrock channel networks is derived from application of geomorphic transport laws designed to represent erosion by water-dominated flows. To quantify the link between debris-flow mechanics and steep channel network form, establishment of a geomorphic transport law that quantifies bedrock erosion by debris flows and development of methods capable of handling the transient, non-local behavior of debris flows throughout the steep channel network are needed. Here, we propose a landscape evolution model to simulate longitudinal channel profiles that evolve in response to both debris-flow and fluvial processes. The model framework, which includes a methodology to explicitly track spatial variations in bulk debrisflow properties (e.g. flow depth, velocity) along the length of the channel profile, is designed to allow for the exploration of a range of potential debris-flow incision laws. We propose a relationship in which the debris-flow erosion rate is a function of debris-flow depth and channel slope. By comparing the morphology of modeled channel profiles to that typically seen in natural channels, we place constraints on a debris-flow incision law with this general form. Modeled channel profiles are consistent with observations of channel morphology in debris-flow-dominated terrain when the debris-flow incision rate is related to local channel slope raised to a power greater than one and approximately linearly related to debris-flow depth. Model results indicate that erosion by debris flows can explain the occurrence of a scaling break in the slope-area curve at low drainage areas and that upper-network channel morphology may be useful for inferring catchment-averaged erosion rates in quasi-steady landscapes. These results improve our ability to interpret topographic signals within steep channel networks and provide a general framework for exploring the role of debris-flow erosion within landscape evolution models.

Introduction

Debris flows are effective at transporting coarse sediment [START_REF] May | Debris flows through different forest age classes in the central Oregon Coast Range[END_REF][START_REF] May | Processes and rates of sediment and wood accumulation in headwater streams of the Oregon Coast Range, USA, Earth Surface Processes and Landforms[END_REF] and eroding bedrock [START_REF] Mccoy | Field measurement of basal forces generated by erosive debris flows[END_REF][START_REF] Mccoy | Infrequent, large-magnitude debris flows are important agents of landscape change[END_REF][START_REF] Stock | Erosion of steepland valleys by debris flows[END_REF] in steep, low-order channels where fluvial sediment transport may 1 https://doi.org/10.5194/esurf-2022-47 Preprint. Discussion started: 25 August 2022 c Author(s) 2022. CC BY 4.0 License. be inhibited by low runoff magnitudes and increases in thresholds for incipient sediment motion (e.g. [START_REF] Lamb | Is the critical Shields stress for incipient sediment motion dependent on channel-bed slope?[END_REF][START_REF] Prancevic | Incipient sediment motion across the river to debris-flow transition[END_REF]. The influence of debris-flow erosion on bedrock channel morphology, including longitudinal channel profiles [START_REF] Montgomery | Channel network source representation using digital elevation models[END_REF], has been recognized for decades. Although debris flows traverse channel networks in many steep landscapes and are capable of eroding bedrock via abrasion and plucking [START_REF] Hsu | Experimental study of bedrock erosion by granular flows[END_REF][START_REF] Stock | Erosion of steepland valleys by debris flows[END_REF][START_REF] Mccoy | Field measurement of basal forces generated by erosive debris flows[END_REF], their relative importance over geologic time compared to other geomorphic processes and the extent to which they affect landscape form at larger scales remains unclear. The spatial extent and magnitude of debris-flow erosion, for example, may be limited in terms of downstream extent due to a lack of mobility or erosive power on modest slopes. Additional work is needed to determine the controls on the magnitude of debris-flow erosion within different parts of the drainage network and the ensuing implications for landscape form [START_REF] Whipple | Bedrock rivers[END_REF].

Topographic signatures of geomorphic processes, which we define as quantitative connections between processes and the morphology of a landform, can be used to infer the presence and rates of geomorphic processes from topographic data. Bedrock channels dominated by fluvial erosion develop longitudinal profiles described by a power law relationship between slope and drainage area. A deviation from this power-law scaling relationship at small drainage areas (Fig. 1) has been interpreted as a topographic signature of debris-flow incision [START_REF] Montgomery | Channel network source representation using digital elevation models[END_REF][START_REF] Seidl | The problem of channel erosion into bedrock[END_REF][START_REF] Sklar | River longitudinal profiles and bedrock incision models: Stream power and the influence of sediment supply[END_REF]Stock andDietrich, 2006, 2003). More recently, the length of the channel network upstream of this transition zone, where debris-flow incision likely plays an important role, was suggested to increase with erosion rate in the San Gabriel Mountains (DiBiase et al., 2012). [START_REF] Penserini | A morphologic proxy for debris flow erosion with application to the earthquake deformation cycle, Cascadia Subduction Zone, USA[END_REF], working in the Oregon Coast Range, further found that there was a positive relationship between catchment-averaged erosion rate and the drainage area at which slope begins to decay as a power law function of drainage area. [START_REF] Stock | Valley incision by debris flows: Evidence of a topographic signature[END_REF] similarly proposed a conceptual model where the transition from a nearly linear debris-flow dominated long-profile to a concave-up fluvial-dominated long-profile migrates out to larger drainage areas as the rock uplift rate increases. Therefore, there is support for the idea that channels eroded by debris flows not only have a unique morphology (or topographic signature) that distinguishes them from purely fluvial channels, but they also record tectonic information.

These findings underscore the need to identify robust topographic signatures of debris-flow erosion, develop mechanistic explanations for the emergence of these signatures, and assess their sensitivity to climatic and tectonic forcing so that they can be exploited to gain process-based insights about the evolution of steep landscapes. In particular, there is a need to understand the relative importance of fluvial and debris-flow processes in setting the location and form of the morphologic transition zone (i.e., the zone where the relatively linear debris-flow-dominated long-profile transitions to a generally concave-up fluvialdominated long-profile) (Fig. 1). For example, the location of this transition may change with uplift rate due to the dynamics of fluvial erosion alone because channel steepness can vary nonlinearly with uplift rate as a result of relationships between runoff variability and fluvial erosion thresholds (e.g. DiBiase and [START_REF] Dibiase | The influence of erosion thresholds and runoff variability on the relationships among topography, climate, and erosion rate[END_REF]Lague, 2014). Debris-flow processes, in contrast, may exert a strong control on the location of the transition by setting the near-uniform slope that the channel approaches at small drainage areas. Critical first steps in understanding these topographic signatures of debris-flow erosion include formulating a debris-flow incision law and developing the methodology needed to incorporate that incision law into landscape evolution models. Previous studies have offered methods for modeling the effects of debris-flow erosion over geologic timescales. [START_REF] Stock | Erosion of steepland valleys by debris flows[END_REF] developed a one-dimensional (1d) landscape evolution model capable of reproducing channel profiles with the characteristic change in slope-area scaling observed in natural environments [START_REF] Montgomery | Channel network source representation using digital elevation models[END_REF][START_REF] Stock | Valley incision by debris flows: Evidence of a topographic signature[END_REF]. [START_REF] Stock | Erosion of steepland valleys by debris flows[END_REF] coupled empirical relationships for debris-flow properties with a debris-60 flow incision law based on inertial stress, but emphasized the need for improved methods to calculate spatially varying bulk debris-flow properties and additional studies to constrain debris-flow incision laws. Several subsequent studies have improved our understanding of the grain-scale processes that control debris-flow incision rates. [START_REF] Hsu | Experimental study of bedrock erosion by granular flows[END_REF][START_REF] Hsu | Mean and fluctuating basal forces generated by granular flows: Laboratory observations in a large vertically rotating drum[END_REF] used physical experiments in a rotating drum to suggest that debris-flow erosion rates scale with inertial stress, which can be cast as a function of flow shear rate. Following the erosion equation proposed by [START_REF] Sklar | A mechanistic model for river incision into bedrock by saltating bed load[END_REF] to account for bedrock 65 erosion due to discrete particle impacts, [START_REF] Mccoy | Controls on erosion and transport of mass by debris flows[END_REF] proposed that the debris-flow erosion rate could be a function of the product of impact frequency, defined as the rate of particle impacts per unit area per unit time, and impact force. Based on a series of discrete element simulations, McCoy (2012) estimated that the product of impact frequency and impact force scaled approximately linearly with flow depth and in a strongly nonlinear way with bed slope, namely with slope to approximately the sixth power. The utility of these relationships in a landscape evolution model, however, has been limited by our ability to simulate the spatial and temporal variability in debris-flow properties (e.g. depth, velocity, shear rate) throughout the channel network in landscape evolution models.

Models designed to represent fluvial erosion within landscape evolution models rely on empirical relationships between local terrain attributes that are readily computed from a digital elevation model (DEM), such as slope and drainage area [START_REF] Tucker | Hillslope processes, drainage density, and landscape morphology[END_REF]. However, incorporation of a debris-flow incision law into this type of local framework is challenging owing to the nonlocal controls on debris-flow erosion including initiation conditions, non-steady flow velocity and the finite and variable runout distance of discrete flows. An additional challenge is the relative paucity of data available to develop empirical relationships between bulk debris-flow properties and terrain attributes. For example, quantifying the frequency at which discrete flows traverse different parts of the landscape has been shown to be a key factor in their ability to sculpt topography [START_REF] Shelef | A unified framework for modeling landscape evolution by discrete flows[END_REF].

Here, we develop a nonlocal modeling framework to predict the evolution of a 1d channel profile by using a process-based debris-flow runout model to route debris flows down the channel and calculate debris-flow properties (e.g. velocity, depth) and their temporal and spatial variations throughout the channel network. We compare results from this process-based runout model with those obtained when using a reduced complexity approach motivated by [START_REF] Gorr | A progressive flow-routing model for rapid assessment of debris-flow inundation[END_REF] that relies on empirical relationships [START_REF] Rickenmann | Empirical Relationships for Debris Flows[END_REF] to route debris-flows down a channel as opposed to solving a set of partial differential equations representing the dynamics of the debris-flow mixture in the process-based model [START_REF] Iverson | Flow of variably fluidized granular masses across three-dimensional terrain: 1. Coulomb mixture theory[END_REF]. The goal of this comparison is to assess the extent to which the empirical approach, which may be more readily integrated into twodimensional (2d) landscape evolution models, yields results that are consistent with the process-based approach. Calculations of spatial variations in debris-flow properties from these two methodologies provide a robust foundation for utilizing relationships between bulk debris-flow properties and particle impact forces [START_REF] Hsu | Experimental study of bedrock erosion by granular flows[END_REF][START_REF] Mccoy | Controls on erosion and transport of mass by debris flows[END_REF] to estimate rates of bedrock incision by debris flows. The main objectives of this study are to (1) develop a methodology to estimate spatial variations in bulk debris-flow properties along a channel profile using information typically available in landscape evolution models; (2) use this methodology to simulate long-profile evolution with a family of debris-flow incision laws based on topographic slope and debris-flow depth, including that proposed by [START_REF] Mccoy | Controls on erosion and transport of mass by debris flows[END_REF]; (3) assess the ability of these incision laws to reproduce channel profiles that have steady-state forms consistent with those in debris-flow-dominated landscapes; and (4) using the bestperforming incision law found in objective 3 and a standard sensitivity analysis, explore controls on the location and form of the topographic signature of debris-flow incision that has been observed in slope-area data.

Methods

Model framework

In the proposed 1d model framework, which is designed to simulate longitudinal channel profiles, the rate of change of elevation (z) with time (t) in a bedrock channel is driven by the uplift rate (U ), fluvial erosion (E f ), and erosion by debris flows (E df ) according to:

∂z ∂t = U -E f -E df . (1)
A brief description of model parameters is included in Table B1. We solve equation 1 numerically on a one-dimensional grid with a uniform spacing of ∆x = 5 m and use the standard explicit forward Euler method for time stepping. We chose a grid spacing of 5 m since it is small enough to resolve changes in the longitudinal channel profile and also large enough to keep model run times, which increase with decreasing grid spacing, manageable.

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

E f = KA ms S ns (2)
where A denotes upstream contributing area as a function of distance, x, from the channel head, S denotes slope, and K, m s = 1.4, and n s = 2.33 are empirical parameters that depend on relationships between discharge and channel width, w, hydraulic geometry and discharge variability, and grain size. Here, we assume that contributing area varies with distance from the channel head according to A = A 0 + 25x 5/3 , with A 0 = 1000 m 2 and that channel width increases with contributing area as w = k w A b , where k w = 0.05 and b = 0.3. Motivated by the geomorphic importance of debris flows in the San Gabriel

Mountains [START_REF] Lavé | Denudation processes and rates in the Transverse Ranges, southern California: Erosional response of a transitional landscape to external and anthropogenic forcing[END_REF], parameters related to channel geometry and fluvial incision were selected based on DiBiase and [START_REF] Dibiase | The influence of erosion thresholds and runoff variability on the relationships among topography, climate, and erosion rate[END_REF] and typical ranges reported by Lague (2014). Additional details on parameter choices for the stream power model are given in Appendix A. Unless noted otherwise, parameters used for the fluvial erosion model are listed in Table B2.

We propose a general formulation that can be used to estimate the erosion rate attributable to debris flows over geologic timescales, E df , at a point on the landscape, given information about the bulk properties of a representative debris flow.

Motivated by the relationships suggested by [START_REF] Mccoy | Controls on erosion and transport of mass by debris flows[END_REF] between impact frequency, impact force, debris-flow depth, and slope, we define E df as a function of slope and debris-flow depth, h. Debris-flow depth varies with position along the channel profile and, at any one location, it also varies with time throughout the course of a debris flow event. Letting t 0 denote the time a debris flow begins moving over a given location along the channel profile and t f the time when it has completely passed that location, then

E df = k df t f t0 S α h β Φdt (3)
where the debris flow erodibility coefficient is defined as k df = κ df F df , κ df is an empirical coefficient related to bedrock and flow properties (e.g. grain size), F df is a term quantifying the frequency of debris flow, Φ is a threshold factor that reflects We also present a reduced-complexity routing algorithm, which closely follows the methodology presented by [START_REF] Gorr | A progressive flow-routing model for rapid assessment of debris-flow inundation[END_REF] to rapidly simulate debris-flow runout for hazard assessment purposes, to compute spatial variations in bulk debrisflow properties along a longitudinal channel profile (Section 2.3). In this approach, we use a set of empirical relationships that relate debris-flow properties to topographic slope [START_REF] Rickenmann | Empirical Relationships for Debris Flows[END_REF] in order to estimate representative values for debris-flow depth, h, at each point along the debris-flow path as well as the time it takes for the debris flow to pass over that point, t ft 0 , which we denote as t p . In other words, h varies along the channel profile but we neglect variations in h that occur within individual debris-flow events (i.e. the rise and fall in flow depth as a debris flow passes over a point on the landscape). In this case, we employ a debris-flow erosion equation analogous to equation 3 that is simplified because h is constant for a given channel location,

E df = k df t p S α h β Θ (4)
where Θ denotes a threshold factor that reflects a reduction in incision when the debris flow is close to rest. To assess the simplifying assumptions of this approach within the context of modeling the evolution of longitudinal channel profiles, we compare the morphology of modeled profiles using this reduced complexity algorithm with profiles generated using the process-based debris-flow routing model.

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

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

∂h ∂t + ∂(hv) ∂x = - vh w ∂w ∂x (5) ∂(hv) ∂t + ∂ ∂x hv 2 + 1 2 g 2 z h = g x h -sgn(v)(1 -λ)g z hφ(I) - 2vηv f ρh - v 2 h w ∂w ∂x . (6)
Here, w is the channel width, h is flow depth, v is velocity, φ(I) is the friction coefficient that depends on the inertial number (I) [START_REF] Jop | A constitutive law for dense granular flows[END_REF], g x and g z denote components of gravity in the x and z directions, respectively, λ = p bed /ρg z h is the ratio of pore fluid pressure to total basal normal stress, v f = 0.5 is the fluid volume fraction, and η is the viscosity of the pore fluid. We assume that the ratio of pore fluid pressure to total basal normal stress decays with time since the debris flow entered the model domain, t, according to

λ = λ 0 1 -erfc 2h √ 4Dt , (7)
where λ 0 = 0.9 and D = 10 -6 m 2 s -1 is the pore fluid pressure diffusivity. This approximation is consistent with an initially high pore fluid pressure shortly following initiation and subsequent linear diffusion of pore fluid pressure [START_REF] Iverson | Flow of variably fluidized granular masses across three-dimensional terrain: 1. Coulomb mixture theory[END_REF] over time.

The friction coefficient is a function of the inertial number [START_REF] Jop | A constitutive law for dense granular flows[END_REF],

φ(I) = µ s + (µ 2 -µ s )/(I/I 0 + 1) (8)
where I = γD ef f /(P/ρ s ) 0.5 , with P denoting the basal normal stress, γ = 2v/h is the shear rate, ρ s = 2600 kg m -3 is the density of sediment, I 0 = 0.279 is a constant, µ s = 0.382, µ 2 = 0.644, and D ef f is a characteristic particle diameter. In this formulation, the friction coefficient increases with the inertial number and approaches µ 2 when I is large.

The governing equations are solved numerically on a grid with uniform spacing. We use a first-order, shock-capturing finite volume method with a Harten-Lax-van Leer-Contact (HLLC) approximate Riemann solver [START_REF] Toro | The HLL and HLLC Riemann solvers[END_REF] to compute the fluxes across each grid cell boundary [START_REF] Mcguire | Constraining the relative importance of raindrop-and flow-driven sediment transport mechanisms in postwildfire environments and implications for recovery time scales[END_REF][START_REF] Mcguire | Debris flow initiation by runoff in a recently burned basin: Is grain-by-grain sediment bulking or en masse failure to blame?[END_REF]. Source terms are treated separately with an explicit, first-order forward Euler method for time stepping.

Debris flows enter the domain through the upper boundary, which can be conceptualized as the channel head, and are routed down the channel profile. We define a series of 20 ghost cells above the uppermost grid cell that effectively extend the model domain for the purpose of initializing a debris flow. Elevations of each ghost cell are determined by assuming that the slopes of all ghost cells are equal to the slope at the uppermost grid cell. Debris flows are initiated from a static pile of debris defined on the ghost cells. This procedure provides some time for debris to begin to flow before it enters the model domain, similar to what might be expected for debris flows that initiate in a colluvial hollow or gully upstream from a channel head. Some studies have reported a growth in debris-flow volume with drainage area [START_REF] Santi | Comparison of debris-flow volumes from burned and unburned areas[END_REF], reflecting entertainment of sediment as the flow moves downstream, but incorporating this effect into the source terms of the process-based routing model is beyond the scope of this study. When using the process-based debris-flow routing model, we assume that debris-flow volume is fixed and does not change along the flow path, although we do explore the effects of spatial variations in debris-flow volume with the empirical routing approach described later.

E df = k df S α k=n k=1 h β k Φ∆t. ( 9
)
where ∆t denotes the time step used when solving the flow equations and we define the threshold factor, Φ, as Φ = 1 when uh > 0.01 m 2 s -1 and Φ = 0 otherwise. To reduce computation time, the debris-flow depth is not updated after each landscape evolution model time step. Instead, a debris flow is routed down the channel profile and the debris-flow incision rate is updated to reflect any changes in flow thickness whenever the channel slope has changed by 0.05 or more at any grid cell since the last time a debris flow was routed. The debris-flow incision rate, however, is updated with every time step of the landscape evolution model to reflect changes in slope.

Estimating debris-flow incision with an empirical routing model

We use a series of empirical relationships defined by [START_REF] Rickenmann | Empirical Relationships for Debris Flows[END_REF] to estimate representative values for debris-flow depth, h, and passage time, t p , at each point along the channel profile based on spatially variable estimates of debris-flow volume, M , debris-flow velocity, v, channel width, w, and topographic slope, S. We assume that debris flows initiate at or above the uppermost grid cell within the computational domain (i.e. the channel head), athough their overall volume may change along the channel profile. We determine the downstream extent of debris-flow runout by treating the debris flow as an idealized fluid with a prescribed yield strength, τ y , and assuming that debris-flow motion stops when shear stress at the base of the flow, τ = ρ b ghS, falls below τ y [START_REF] Whipple | The influence of debris-flow rheology on fan morphology, Owens Valley, California[END_REF][START_REF] Gorr | A progressive flow-routing model for rapid assessment of debris-flow inundation[END_REF]. Here, g = 9.81 m s -2 denotes gravitational acceleration and ρ b = 1800 kg m -2 is the bulk density of the debris flow. In practice, we determine h, t p , and τ everywhere in the model domain, determine the downstream extent of debris-flow runout, and then apply equation 4 to compute a non-zero value for E df only along the debris-flow travel path.

To begin, we specify debris-flow volume at each grid cell as a function of upstream drainage area (A) according to M = M 0 (10 -6 •A) γ [START_REF] Santi | Comparison of debris-flow volumes from burned and unburned areas[END_REF]. This formulation assumes that debris-flow volume increases downstream, reflecting entrainment of bed material or lateral inflow, but these volume changes can be neglected by setting γ = 0. The empirical coefficient, M 0 , and exponent, γ, may vary considerably among landscapes. [START_REF] Santi | Comparison of debris-flow volumes from burned and unburned areas[END_REF] estimated M 0 = 3358 and γ = 0.73 using data throughout the western and southwestern United States, M 0 = 10470 and γ = 0.62 based on data from the Italian Alps, and M 0 = 18770 and γ = 0.28 using data from the northwestern United States and southwestern Canada. Peak debris flow discharge can then be computed according to

Q = c 1 M c2 , (10)
where c 1 = 0.135 and c 2 = 0.78 are empirical coefficients [START_REF] Rickenmann | Empirical Relationships for Debris Flows[END_REF]. Noting that Q = wvh and using the relationship [START_REF] Rickenmann | Empirical Relationships for Debris Flows[END_REF]) where µ denotes the dynamic viscosity of the flow, it is possible to solve for flow depth,

v = 1 3µ ρ b gh 2 S, (11
h = 3µc 1 M c2 ρ b gSw 1/3 . ( 12
)
Using the relationships between channel width, w, and area, A, and debris-flow volume, M , and area, A, a representative flow depth for a given location along the channel profile can be written in terms of area and slope,

h = 3µc 1 (M 0 A γ ) c2 ρ b gSk w A b 1/3 . ( 13
)
We treat this flow depth as a representative value for each point in the drainage network but acknowledge that it may overestimate flow depth because equation 10 is used to estimate peak debris flow discharge. We further define the passage time of the debris flow as

t p = M Q = M 0 A γ c 1 (M 0 A γ ) c2 . ( 14
)
Finally, we define the threshold factor, Θ, such that the debris-flow incision rate decreases as the flow approaches the end of its travel path and the shear stress at the base of the flow approaches the yield strength. Specifically,

Θ = 1 - τ y τ (15)
The debris-flow erosion rate can then be determined according to equation 4.

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

Numerical experiments

Our numerical experiments have two goals, which will be treated in turn. First, we assess which erosion laws, as defined by different values of α and β, can reproduce the first-order characteristics of observed channel longitudinal profiles as well as how this may be affected by the choice of debris-flow routing model (i.e. process-based or empirical). Second, we aim to understand the sensitivity of predicted channel profiles to model parameters.

A family of debris-flow incision laws

We explore model behavior for different values of α and β in the family of incision laws described by equations 3 and 4 by comparing modeled, steady-state longitudinal profiles with those typical of debris-flow-dominated terrain (Fig. 1). We do not try to recreate the channel morphology observed within specific watersheds or geographic regions. Instead, we aim to provide some constraints on α and β by identifying ranges for these two exponents that result in longitudinal channel profiles that are consistent with observed changes in the relationship between slope and contributing area in natural channels traversed by debris flows [START_REF] Stock | Valley incision by debris flows: Evidence of a topographic signature[END_REF]. In particular, there is a power law relationship between slope and contributing area in portions of the drainage network where erosion is dominated by fluvial processes. This power law scaling breaks down at smaller drainage areas in channels that are traversed by debris flows [START_REF] Stock | Valley incision by debris flows: Evidence of a topographic signature[END_REF], where slope becomes less sensitive to variations in contributing area (Fig. 1). Channel slope may increase slowly or remain constant as drainage area decreases toward the channel head. [START_REF] Stock | Valley incision by debris flows: Evidence of a topographic signature[END_REF] found that the shape of longitudinal channel profiles that experience both debris flow and fluvial erosion can be described by a family of curves with the general form

S = S df 1 + a 1 A a2 (16)
where S df , a 1 , and a 2 are empirical coefficients (Fig. 1). The coefficient S df is related to the slope that the channel network approaches at low contributing areas, a 2 controls the power law relationship between slope and contributing area when A is large, and a 1 (which has units of 1/(length 2 ) a2 ) controls the sharpness of the transition from fluvial power law scaling at large drainage areas to relatively constant slopes in debris-flow-dominated reaches at smaller drainage areas. Here, we work with a modified form of the above expression for channel slope, given by

S = S df 1 + (A/A df ) p (17)
which is advantageous for the present study because A df has units of length 2 and can be interpreted as the drainage area at which the slope-area relationship transitions from a constant slope with decreasing drainage area to the standard power law relationship expected in the fluvial network. In this sense, A df provides one metric for identifying the transition between the debris-flow domain, with a characteristic gradient tending towards S df , and the fluvial process domain.

A landscape evolution model designed to simulate the evolution of channel longitudinal profiles in response to both debris flow and fluvial erosion should produce steady-state channel profiles that are well described by equation 17. We assess model performance for different α and β in two ways. First, we compute the R 2 associated with the best fit to equation 17. We allow S df , A df , and p to vary freely when fitting modeled slope-area data to equation 17. Second, we examine the difference between the maximum slope along the channel profile, S max , and the slope at the channel head, S ch . The second criteria focuses on checking a basic morphologic property observed in natural channels, namely that channel slope generally increases or remains constant as drainage area decreases in quasi-steady-state landscapes.

We assessed performance of the landscape evolution model with different values of α and β when using the process-based routing model and when using the empirical routing model. Using the process-based model, we performed a numerical experiment where we varied α, β, pore pressure diffusivity (D), viscosity of the pore fluid (η), the friction parameter (µ 2 ), the debris flow erodibility coefficient (k df ), and instantaneous fluvial erodibility coefficient (k e ) within the ranges specified in Table B3.

We allowed some variation in model parameters other than α and β to ensure trends between α, β, and model performance metrics were not specific to a particular subset of the parameter space. We selected 500 parameter sets using a Latin Hypercube sampling strategy. We performed an analogous numerical experiment using the empirical routing model where we sampled B4). We were able to perform a greater number of simulations using the empirical model because it is less computationally demanding.

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

S = (U/K) 1/ns A -ms/ns . ( 18
)
Simulations ended once an approximate steady state had been reached, which typically took 10 6 -10 7 years.

Sensitivity analysis

We performed sensitivity analyses using both the process-based and empirical routing models to explore how the topographic signature of debris-flow incision is likely to be expressed in different settings. Motivated by [START_REF] Mccoy | Controls on erosion and transport of mass by debris flows[END_REF] and the results of our numerical experiments to constrain α and β, we set α = 6 and β = 1 in these experiments. We focus, in particular, on understanding relationships between model parameters, A df , and S df in steady-state longitudinal channel profiles.

To perform the sensitivity analysis with the process-based routing model, we used a Latin Hypercube sampling strategy to select 1500 sets of parameters where instantaneous fluvial erodibility (k e ), debris flow erodibility (k df ), viscosity of the pore fluid (η), pore fluid diffusivity (D), the friction factor (µ 2 ), and uplift rate (U ) varied within the ranges defined in Table B5.

The sensitivity analysis using the empirical routing approach was analogous, but we were able to perform a greater number of simulations. We used a Latin Hypercube sampling strategy to select 7000 sets of parameters where instantaneous fluvial erodibility (k e ), debris flow erodibility (k df ), viscosity (µ), yield strength (τ y ), uplift rate (U ), and debris-flow volume parameters M 0 and γ varied within the ranges defined in Table B6.

Results

Constraints on a debris-flow incision law

Process-based routing model

At large drainage areas, modeled profiles exhibit a power law scaling between slope and drainage area that is expected based on the fluvial incision law (Figs. 2,3). The R 2 value associated with a fit to equation 17 was greater than 0.95 for approximately 93% of the modeled profiles. At drainage areas of approximately 0.1 km 2 , however, some parameter combinations produced channel profiles where slope began to decrease as drainage area decreased. In other words, the difference between the maximum slope, S max , along the channel profile and the slope at the channel head, S ch , was positive and regularly exceeded 0.2 in cases where β > 2 (Fig. 2). This decrease in slope at low drainage areas is inconsistent with observations that indicate slope continues to increase or remain constant as drainage area decreases, as described by equation 17. Differences between S max and S ch decreased rapidly as α/β increased (Fig. 2).

Empirical routing model

Modeled profiles exhibited the expected power law scaling between slope and drainage area at large drainage areas where fluvial incision dominated debris-flow incision. The coefficient of determination (R 2 ) value associated with a fit to equation 17 was greater than 0.95 for approximately 99% of the modeled profiles. As with results obtained using the process-based routing model, some parameter combinations produced channel profiles where the maximum channel slope was not observed at the 315 channel head (Fig. 3). Differences between S max and S ch are minor when α = 6 and β = 1, but become more substantial as α decreases and/or as β increases (Fig. 3). More generally, the extent to which modeled channel profiles exhibit a decrease in slope at small drainage areas depends on α, β, and the exponent γ that controls the relationship between debris-flow volume and drainage area (Fig. 4). In cases where γ < 0.3, numerous combinations of α and β lead to decreases in slope at low drainage areas. For any choice of 2 ≤ α ≤ 8 and β ≈ 1.5 or less, S max -S ch was always less than 0.05 (Fig. 4). For cases where β > 2, there was no choice of α between two and eight that maintained S max -S ch < 0.05 for all values of γ. Differences between S max and S ch increase as β increases and/or as α decreases. We were unable to explore the effects of spatial variations in debris-flow volume using the process-based routing model, where we neglect changes in debris-flow volume along the flow path (i.e. γ = 0). The difference between the maximum slope (Smax) and the slope at the channel head (S ch ) generally decreases with α, the slope exponent in the debris-flow erosion law as seen when comparing across panels a-f, and increases with β , the depth exponent in the debris-flow erosion law (as seen in marker color). Model performance, as measured by Smax -S ch , is most sensitive to changes in α and β when γ is less than approximately 0.3.

Steady-state forms of channel profiles

Process-based routing model

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

The sensitivity of steady-state long-channel profiles to changes in uplift rate leads to power law relationships between A df and U and between S df and U (Fig. 7). By randomly sampling model parameters, including uplift rate, within prescribed ranges, we assume that none are correlated with each other. However, this is unlikely in natural landscapes and correlations are expected. For example, we may expect that debris-flow frequency, F df , increases with uplift. We explore the consequence of such a correlation on the emergent channel morphology by sub-sampling results to only include parameter sets that meet specific criteria. To force a linear increase in k df with U , which could result from an increase in debris-flow frequency with U , we selected simulations where 4 • 10 -5 + 6 • 10 -5 U ≤ k df ≤ 6 • 10 -5 + 6 • 10 -5 U . Forcing this linear relationship between k df and U leads to an increase in the power law exponent relating U to A df and a decrease in the exponent relating U to S df (Fig. 7).

Empirical routing model

Simulations indicate that S df is most sensitive to changes in k df (e.g. debris-flow frequency, erodibility) and γ, which controls the relationship between debris-flow volume and drainage area (Fig. 8) as well as uplift rate (Fig. 9). Typical values of S df decrease with k df and increase with γ but appear independent of k e and M 0 . The area at which there is a transition to fluvial power law scaling between slope and area, A df , appears most sensitive to γ, k df , k e , and U whereas it is relatively insensitive to M 0 (Figs. 8,9). Mean values of A df tend to decrease strongly with k e , increase with k df , and decrease with γ. Parameters more directly related to the physical properties of the debris flows, viscosity (µ) and yield strength (τ y ), had relatively minor control over S df and A df (Figs. 8,9).

There is a power-law relationship between A df and uplift (U ), specifically A df = 0.52U 0.42 although there is considerable scatter due to the wide range of parameter values included in the sensitivity analysis (Fig. 10). By randomly sampling the model parameters within prescribed ranges, we assume that none are correlated with each other. However, this is unlikely in natural landscapes and correlations are expected. For example we may expect that F df and/or γ increase with uplift. Again, we explore the consequence of such correlations by sub-sampling results to only include parameter sets that met specific criteria.

To force a linear increase in k df with U , which could result from an increase in debris-flow frequency with U , we selected simulations where 0.003 + 0.0037U ≤ k df ≤ 0.005 + 0.0037U . Fitting a power law relationship between A df and U based on these 258 sub-sampled data points results in a greater power law exponent (Fig. 10), specifically A df = 0.66U 0.67 . The exponent in the power law describing the relationship between A df and U therefore depends on the relationship between U and k df , and is greater when k df increases with uplift rate and lower when k df and uplift rate are independent.

A linear increase in k df with uplift rate results in a decrease in the best fit power law exponent describing the relationship between S df and U . In particular, S df = 0.51U 0.17 when k df and U are independent while S df = 0.45U 0.03 when k df increases linearly with U . We similarly observe shifts in the power law relationships between A df and U and between S df and U when data are sub-sampled to force a linear increase in the volume-area scaling exponent (γ) with uplift rate by requiring 0.2 + 0.5U ≤ k df ≤ 0.3 + 0.5U (Fig. 10). More specifically, the power-law exponent relating A df and U decreases from 0.42 to 0.14 whereas the exponent relating S df and U increases from 0.17 to 0.3.

Discussion

Constraints on a geomorphic transport law for debris-flow incision

Results indicate that the debris-flow incision law proposed by [START_REF] Mccoy | Controls on erosion and transport of mass by debris flows[END_REF] based on the changing frequency and magnitude of grain-scale bed impact forces observed beneath a variety of granular flows, namely α ≈ 6 and β ≈ 1, produces channel profiles that are consistent with observations from natural landscapes (Figs. 2,3,C2). This is true within a wide range of the parameter space explored here, including for a range of γ that covers the variability observed across several different geographic regions reported by [START_REF] Santi | Comparison of debris-flow volumes from burned and unburned areas[END_REF]. Although some combinations of α and β produce substantial decreases in slope at small drainage areas, indicating they are not reasonable combinations for a debris-flow incision law, values of α = 6 and β = 1 are not unique in their ability to lead to reasonable model performance. Other combinations of α and β also yielded profiles with high R 2 along with infrequent instances where slope decreases substantially with decreases in drainage area. Data and numerical experiments presented here are not capable of differentiating among these potential debris-flow incision laws, although cases where α < 3 and/or β > 2 generally performed poorly (Figs. 2,3).

To produce the observed steady-state morphology of debris-flow dominated long-profiles with a slope that is approximately constant or slowly decreasing with A, the erosion rate driven by debris flows (E df ) must remain constant or decrease slightly with A. Assuming that the role of the threshold term in E df is negligible over the majority of the debris-flow dominated reach, show a power law relationship between uplift (U ) and A df (dashed best fit curve) and between U and S df . The exponent in the power law relationship changes when sub-sampling from the 7000 trials to only include data points within the parameter space such that there is a linear increase in k df as U increases (orange dots, solid black line). (c,d) The exponent in the power law best fit is also sensitive to changes in the relationship between U and the volume exponent (γ), as illustrated when sub-sampling to include only points within the parameter space such that there is a linear increase in γ as U increases (orange dots, solid black line).

0.38 in order for h, as given by equation 13, to be constant or increase with increasing A. This is consistent with modeled channel profiles where the maximum slope is achieved at the channel head for γ > 0.3 (Fig. 4). Additional work, particularly targeted at landscapes where there are tighter controls on assumptions that influence model outcomes such as relationships between debris-flow volume and drainage area, is needed to further constrain a debris-flow incision law. Analogously, there are a range of exponents used in the generalized stream power incision law for fluvial erosion and work continues in an effort to constrain those exponents (e.g. [START_REF] Clubb | The relationship between drainage density, erosion rate, and hilltop curvature: Implications for sediment transport processes[END_REF][START_REF] Turowski | Alluvial cover controlling the width, slope and sinuosity of bedrock channels[END_REF]Turowski, , 2021)).

In general, the proposed empirical and process-based approaches for estimating bulk debris-flow properties along the channel ,9). These similarities in model behavior are encouraging because the empirical routing approach provides a framework to estimate bulk debris-flow properties using quantities that can be computed from a DEM, specifically upstream contributing area and slope, that could be used in future efforts to more efficiently explore alternative debris-flow incision laws. [START_REF] Stock | Erosion of steepland valleys by debris flows[END_REF], for example, present a debris-flow incision law based on inertial stress that would require an estimate for flow velocity as well as flow depth. Although we do not directly use velocity estimates in this work, it would be straightforward to obtain estimates of velocity along the debris-flow path using Eq. 11.

Steady-state forms of longitudinal channel profiles

Model results support a conceptual model in which both debris-flow and fluvial erosion play critical roles in setting longitudinal profile form in the upper channel network. Changes in parameters related solely to fluvial erosion do not have a strong influence on S df , which simulations demonstrate is primarily controlled by changes in debris-flow processes (Figs. 5,6,8,9).

Specifically, increases in uplift rate U and γ, or decreases in k df promote increases in S df , assuming all else is fixed. In contrast, A df is controlled by a combination of fluvial and debris-flow processes (Figs. 5,6,8,9) ,10b). Thus, two avenues are essential to allow a more detailed comparison between natural data and the modeling framework presented here, namely (1) extracting upland channel morphology in a larger number of catchments with constrained erosion rates, and (2) gathering evidence on the interconnections of key parameters [START_REF] Stock | Erosion of steepland valleys by debris flows[END_REF]. watersheds with catchment-averaged erosion rates determined from cosmogenic radionuclide analysis. Cast in terms of the morphologic variables used here to describe channel profiles, specifically using A df in place of a 1 , the results from [START_REF] Penserini | A morphologic proxy for debris flow erosion with application to the earthquake deformation cycle, Cascadia Subduction Zone, USA[END_REF] indicate that A df increases with E. Simulations confirm this pattern of increasing A df with E, but also highlight the importance of constraining relationships between E and k df as well as E and γ in order for A df to serve as a proxy for erosion rates in an absolute sense (Figs. 7,10). This indicates a need to constrain relationships between E and debris-flow frequency and bed erodibility (which together control k df ) as well as E and the rate at which debris-flow volume increases with contributing area (e.g. γ) rather than potential relationships between E and debris-flow mobility parameters, such as viscosity and yield strength. Although not explored here, differences in m s and n s , which vary among landscapes (Lague, 2014), are also likely to influence relationships between E and A df . Similarly, several landscape parameters assumed fixed in this study, such as the area at the channel head, A 0 , or the channel width scaling, b, may be influenced by debris-flow erosion and exert control over the long-term channel morphology.

Model applications and limitations

We describe a model that provides a framework for exploring the effect of episodic debris flows on channel longitudinal profiles.

The model also comes with several limitations. We assume that all debris flows initiate at the channel head, whereas debris-flow initiation locations in natural landscapes will be more varied. Past work highlights the role that network structure, specifically the number of debris-flow initiation locations upstream from a given channel reach, may play in controlling channel form [START_REF] Stock | Erosion of steepland valleys by debris flows[END_REF]. Variations in the spatial distribution of debris-flow initiation locations within a watershed could be explored within the model presented here by prescribing debris-flow frequency as a function of drainage area or explicitly modeling multiple debris flows from different initiation locations. Additionally, the scaling between channel width and drainage area may differ in the upper portion of the channel network from previously reported relationships that are derived from data at larger drainage areas where fluvial processes are dominant (DiBiase and [START_REF] Dibiase | The influence of erosion thresholds and runoff variability on the relationships among topography, climate, and erosion rate[END_REF]. Similarly, when using the empirical debris-flow routing model, we rely on relationships between debris-flow volume and drainage area that were derived using data collected primarily at drainage areas greater than 0.1 km 2 [START_REF] Santi | Comparison of debris-flow volumes from burned and unburned areas[END_REF]. Lastly, the process-based debris-flow routing model presented here assumes that debris-flow volume is constant and does not change along the travel path, although quantifying controls on sediment entrainment by debris flows and incorporating entertainment into process-based debris-flow routing models are areas of active research [START_REF] Iverson | Elementary theory of bed-sediment entrainment by debris flows and avalanches[END_REF][START_REF] Iverson | Entrainment of bed material by Earth-surface mass flows: Review and reformulation of depth-integrated theory[END_REF][START_REF] Mccoy | Sediment entrainment by debris flows: In situ measurements from the headwaters of a steep catchment[END_REF][START_REF] Haas | Bed scour by debris flows: experimental investigation of effects of debris-flow composition[END_REF]. Advances in our understanding of how debris flows entrain sediment could allow for more detailed comparisons between empirical and process-based approaches to sediment bulking in the proposed landscape evolution model.

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

Conclusions

We present a novel framework for incorporating erosion by debris flows into a 1d landscape evolution model. We demonstrate two methods to estimate debris-flow runout and bulk debris-flow properties (e.g. depth, velocity) throughout the channel network, one based on a process-based debris-flow routing model and the other based on an empirical routing approach.

Combined with a geomorphic transport law describing the relationship between debris-flow depth, channel slope, and a debrisflow incision rate, we are able to quantify spatial variations in the debris-flow incision rate throughout the channel network.

We explore the performance of a family of potential debris-flow incision laws by comparing the form of modeled longitudinal channel profiles with those typically observed in debris-flow-dominated landscapes. Results demonstrate that a debris-flow incision law based on flow depth, slope, and debris-flow passage time is reasonable given general constraints on empirical

Figure 1 .

 1 Figure 1. Examples of channel profiles from the San Gabriel Mountains, California, USA, along with best fit curves of the form S = S df /(1 + (A/A df ) p ). The locations of red circles coincide with A df . Channel profiles were extracted from catchments with different erosion rates, E. Watershed outlets are located at 396964 m E 3799133 m N (black line), 385052 m E 3799080 m N (blue line), and 417980 m E 3792440 m N (orange line). Catchment-averaged erosion rates, E, are from DiBiase et al. (2010). Slope and drainage area values were extracted from the channel network and separated into 100 logarithmically spaced bins. Binned slope was aggregated using the mean slope within each bin.

5

  https://doi.org/10.5194/esurf-2022-47 Preprint. Discussion started: 25 August 2022 c Author(s) 2022. CC BY 4.0 License. a reduction in incision when the debris flow is close to rest, and α and β are empirical exponents. This formulation can be extended to account for a distribution of representative debris flows with different properties given information about their relative recurrence. With the process-based debris-flow routing model (Section 2.2), we compute time-varying flow properties required to determine a debris-flow incision rate at each point along the channel profile using equation 3.

6

  https://doi.org/10.5194/esurf-2022-47 Preprint. Discussion started: 25 August 2022 c Author(s) 2022. CC BY 4.0 License. The first, second, and third source terms on the right hand side of the momentum conservation equation (equation 6) account for variations in bed topography, frictional resistance associated with the solid phase of the flow, and viscous resistance associated with the fluid phase, respectively. The remaining source terms in the mass and momentum equations account for variations in channel width.

7

  https://doi.org/10.5194/esurf-2022-47 Preprint. Discussion started: 25 August 2022 c Author(s) 2022. CC BY 4.0 License. At each grid cell in the model domain (i.e. excluding ghost cells), the debris-flow incision rate is computed using equation 3 based on the time-varying values of debris-flow depth. More specifically, for a debris flow simulated over k timesteps,

  ) 8 https://doi.org/10.5194/esurf-2022-47 Preprint. Discussion started: 25 August 2022 c Author(s) 2022. CC BY 4.0 License.

  https://doi.org/10.5194/esurf-2022-47 Preprint. Discussion started: 25 August 2022 c Author(s) 2022. CC BY 4.0 License.

7000

  different parameter sets with varying values of α, β, instantaneous fluvial erodibility (k e ), debris flow erodibility (k df ), 10 https://doi.org/10.5194/esurf-2022-47 Preprint. Discussion started: 25 August 2022 c Author(s) 2022. CC BY 4.0 License. viscosity (µ), yield strength (τ y ), and debris-flow volume parameters M 0 and γ within prescribed ranges (Table

Figure 2 .

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

Figure 3 .

 3 Figure 3. Numerical experiments using the empirical routing model to determine which slope (α) and depth (β) exponents in the debris-flow erosion law create profiles consistent with those seen in nature. (a) The difference between the maximum slope (Smax) and the slope at the channel head (S ch ) generally decreases with α. Color indicates the value of β, highlighting poor model performance, as measured by Smax -S ch , in many cases when β increases. (b) The difference between the maximum slope and the slope at the channel head generally increases with β. Color indicates the value of α, with greater values of α generally leading to better model performance. (c) The difference between the maximum slope and the slope at the channel head decreases rapidly with α/β. (d-i) Representative profiles for different α and β.

Figure 4 .

 4 Figure 4. Numerical experiments using the empirical routing model that highlight the importance of the volume-area scaling exponent γ.

  14 https://doi.org/10.5194/esurf-2022-47 Preprint. Discussion started: 25 August 2022 c Author(s) 2022. CC BY 4.0 License.

Figure 5 .

 5 Figure 5. Scatter plots summarizing results of the sensitivity analysis with the process based debris-flow routing model. Sensitivity of A df (left column) and S df (right column) to particular model parameters is indicated when there is a gradient in color, whereas plots with no spatial pattern in color indicate a lack of sensitivity. The relationship between A df , ke, and different parameters related to debris-flow processes illustrate sensitivity to the debris-flow erodibility coefficient, k df , and the instantaneous fluvial erodibility coefficient (ke). The morphologic parameter S df (right column) is most sensitive to k df .

Figure 6 .

 6 Figure 6. Steady-state longitudinal profiles produced by the process-based model for various parameter combinations. Variations in (a) debris flow erodibility coefficient, k df , (b) instantaneous fluvial erodibility, ke, and (c) uplift, U , drive greater changes in channel morphology relative to variations in (d) viscosity of the pore fluid, η, (e) friction factor, µ2, and (f) pore pressure diffusivity, D. Default parameter values are: k df = 7.5 • 10 -5 ; ke = 5 • 10 -14 ; U = 0.0006; η = 60; µ2 = 0.62; µ = 5 • 10 -6 .

Figure 7 .

 7 Figure 7. Results from the 1500 process-based model simulations in the sensitivity analysis with α = 6 and β = 1 (blue dots) show a power law relationship between (a) uplift (U ) and A df (dashed best fit curve) as well as (b) U and S df . The exponent in the power law relationship changes when sub-sampling from the 322 trials to only include data points within the parameter space such that there is a linear increase in k df as U increases (orange dots, solid black line).

Figure 8 .

 8 Figure 8. Scatter plots summarizing results of the sensitivity analysis with the empirical debris-flow routing model. The relationship between A df , ke, and different parameters related to debris-flow erosion (left panels) illustrate sensitivity to the volume-area scaling exponent (γ), debris flow erodibility coefficient, k df , and the instantaneous fluvial erodibility coefficient (ke). The morphologic parameter S df (right panels) is most sensitive to the volume exponent (γ) and k df .

Figure 9 .

 9 Figure 9. Steady-state longitudinal profiles produced by the empirical routing model for various parameter combinations. Variations in (a) the debris flow erodibility coefficient, k df , (b) the instantaneous fluvial erodibility, ke, and (c) uplift, U , and (d) the volume-area scaling exponent, γ, drive greater changes in the channel morphology, as summarized by the relationship between channel slope and drainage area, relative to variations in (e) the volume-area scaling coefficient, M0, and (f) viscosity, η. Default parameter values are: k df = 0.0025; ke = 5 • 10 -14 ; U = 0.0006; γ = 0.5; M0 = 3000; µ = 1100.

380

  this requires that h remains constant or increases with A. Examination of equation 13 implies that γ = b/c 2 for h to stay constant as A increases. This equivalence would imply that the increase in channel width with increasing A is matched by an increase in debris-flow discharge. The chosen values of c 2 = 0.78 and b = 0.3 imply that γ must be greater than approximately

Figure 10 .

 10 Figure 10. (a,b) Results from the 7000 empirical routing model simulations in the sensitivity analysis with α = 6 and β = 1 (blue dots)

390

  profile do not appear to result in different model behavior(Figs. 2, 3, 6, 9). It is not possible to directly compare the longitudinal profiles produced by the two different routing models with the same values of α and β because the parameters that affect bulk debris-flow properties are different among the two models. For example the process-based model has no yield strength parameter, yet this parameter plays a key role in determining debris-flow runout in the empirical model. Although 20 https://doi.org/10.5194/esurf-2022-47 Preprint. Discussion started: 25 August 2022 c Author(s) 2022. CC BY 4.0 License. the flow depth, velocity, and passage time predicted by the two routing models will undoubtedly vary, these variations do not appear to be sufficient to alter the extent to which different values of α and β produce modeled profiles that are consistent or inconsistent with observations of channel morphology in debris-flow-dominated terrain (Figs. 2, 3). Furthermore, results using the process-based and empirical routing approaches both highlight the sensitivity of A df and S df to uplift, k e , and k df relative to model parameters related to flow mobility (Figs. 6

  Simulations, assuming α = 6 and β = 1, indicate that debris flows may frequently traverse channel reaches at larger drainage areas without influencing longitudinal channel profile form in a substantial way. Debris flows regularly traversed the entire 21 https://doi.org/10.5194/esurf-2022-47 Preprint. Discussion started: 25 August 2022 c Author(s) 2022. CC BY 4.0 License.model domain (8 km 2 ), but A df was generally less than 0.4 km 2 (Fig.C1). A primary reason for this is the sensitivity of the debris-flow incision law to slope, which means that portions of the channel profile could regularly be traversed by debris flows but they would do little work to erode bedrock when the channel slope is modest. This result indicates that the presence of debris-flow activity in natural channels, as indicated by debris-flow deposits, may not be a reliable indicator of the importance of debris-flow incision. However, additional work is needed to further constrain the slope exponent, α, in the debris-flow incision law, as well as to explore the influence of debris flows mobilizing large-caliber sediment below A df that would otherwise shield the bed from subsequent fluvial erosion. An additional consequence of the nonlinear (α >> 1) relationship between slope and debris-flow incision and the observation that the majority of debris flows remain mobile throughout the entire model domain, is that A df and S df are less sensitive to parameters related to debris-flow properties, namely viscosity and yield strength in the case of the empirical routing model (Figs.8, 9) and viscosity of the pore fluid, friction parameters, and pore pressure diffusivity in the case of the process-based routing model(Figs. 5, 6). As long as debris flows are sufficiently mobile to traverse moderate slopes, these parameters primarily affect the debris-flow incision rate by changing debris-flow depth and/or passage time. Because both debris-flow depth and passage time are linearly related to the debris-flow incision rate (when β = 1), a factor-of-two increase or decrease in flow depth or passage time would only require a relatively small adjustment in slope to compensate for the ability of debris flows to balance the imposed uplift rate at steady state.4.3 Tectonics from debris-flow processes and topographyModel results support previous observations indicating that the morphology of channel profiles in debris-flow-dominated landscapes may provide constraints on erosion rates in steady-state landscapes(Figs. 7, 10).[START_REF] Penserini | A morphologic proxy for debris flow erosion with application to the earthquake deformation cycle, Cascadia Subduction Zone, USA[END_REF] document an inverse relationship between a 1 and E in the central Oregon Coast Range based on analyses of channel profiles in six

  22 https://doi.org/10.5194/esurf-2022-47 Preprint. Discussion started: 25 August 2022 c Author(s) 2022. CC BY 4.0 License.

  . On average, increases in the instantaneous fluvial erodibility lead to decreases in A df whereas increases in the debris-flow erodibility coefficient promote increases in A df . If S df was a constant set, for example, by soil geotechnical properties related to slope stability, then A df could be estimated by the area at which the steady-state fluvial channel gradient reaches S df . However, such simple threshold behavior of S df is not what we observe. Rather, simulations demonstrate that S df varies with U , γ and k df and independent changes to either debris-flow incision processes or fluvial processes are sufficient to influence A df . Thus, accounting for both debris-flow and fluvial processes is important to understand how the steep channel network will respond to changes in tectonic or climatic forcing.It is noteworthy that[START_REF] Penserini | A morphologic proxy for debris flow erosion with application to the earthquake deformation cycle, Cascadia Subduction Zone, USA[END_REF] found no systematic variation in S df with erosion rate in the Oregon Coast Range, whereas our model results indicate increases in erosion rate, or equivalently uplift rate in a steady-state landscape, lead to an increase in S df . However, the lack of any relationship between S df and erosion rate (E) in the Oregon Coast Range may result from a correlation between E and debris-flow frequency and/or E and debris-flow volume. Simulations of steady-state channel profiles indicate that S df increases with uplift rate in cases where there is no relationship between uplift rate (U ) and k df but that S df is roughly independent of U when k df increases linearly with U (Figs. 7b

Acknowledgements. This material is based upon work supported by the National Science Foundation under Grant No. 1951274. Jason Kean, Francis Rengers, Leslie Hsu, Ryan Gold, and Janet Carter provided comments as part of U.S. Geological Survey review. 30 https://doi.org/10.5194/esurf-2022-47 Preprint. Discussion started: 25 August 2022 c Author(s) 2022. CC BY 4.0 License.

exponents that relate flow depth and local channel slope to the incision rate. Simulations indicate that both A df and S df have potential to serve as a morphologic proxy for the catchment-averaged erosion rate. However, both A df and S df appear sensitive to debris-flow frequency and debris-flow erodibility (k df ), and the empirical exponent characterizing how debrisflow volume increases with contributing area (γ), indicating that the utility of such a proxy would depend on the extent to which relationships between erosion rate, k df , and γ could be constrained. Results provide a general form for a debris-flow incision law, a framework that can be used to explore alternative debris-flow incision laws, as well as insight into the relative importance of debris-flow versus fluvial processes in shaping channel profiles in steep landscapes. The results pave the way to future inclusion of debris-flow incision into 2d landscape evolution models and improved interpretations of high-resolution topographic signatures in rapidly eroding landscapes.

Appendix A: Stochastic stream power incision model

The parameterization for the stochastic stream power model is not tuned to any particular landscape or geographic region, but relies on values and relationships that are based on typical values reported by Lague (2014) for high discharge settings and by DiBiase and [START_REF] Dibiase | The influence of erosion thresholds and runoff variability on the relationships among topography, climate, and erosion rate[END_REF] for the San Gabriel Mountains. Given the parameters listed in Table B2, we follow Lague (2014) and compute the critical shear stress for bed load transport, τ c according to τ c = 0.045g(ρ sρ w )D ef f , with D ef f = 0.09 the effective grain size. The stochastic-threshold prediction for the slope exponent in the fluvial incision law, n s , is given by n s = β f /α f (k + 1)/(1ω s ), where β f = 0.7 is the slope exponent in the hydraulic friction law, α f = 0.6 is the discharge exponent in the hydraulic friction law, ω s = 0.25 is the at-a-station width scaling exponent, and k = 0.5 is the discharge variability coefficient. The prediction for the area exponent is given by m s = (cb)(k + 1)/(1ω s ), where c = 1 is the mean discharge-area scaling exponent and b = 0.3 is the width-area scaling coefficient (Lague, 2014). Lastly, the fluvial erodibility coefficient is calculated as (Lague, 2014)

c /k w denotes a width factor, and

Here, k t = gρ w N 3/5 , N = 0.05 denotes the Manning friction coefficient, and a = 1.5 is a shear stress exponent. The rate of fluvial erosion is computed then computed according to E f = KA ms S ns .

Appendix B: Model parameters

This appendix contains tables with information about model parameter for different numerical experiments. A df = 0.38 km 2 , at which there is a transition from a near-constant slope to a the power-law relationship between slope and drainage area.

Once the channel profile has reached a steady state, debris flows traverse the entire channel profile in approximately 88% of simulations and stop at drainage areas less than 1 km 2 in fewer than 1% of simulations. Competing interests. The authors declare that they have no competing interests.