A Modular, Non-Newtonian, Model, Library Framework (DebrisLib) for Post-Wildfire Flood Risk Management

Wildfires increase flow and sediment load through removal of vegetation, alteration of soils, decreasing infiltration, and production of ash commonly generating a wide variety of geophysical flows (i.e., hyperconcentrated flows, mudflows, debris 10 flows, etc.). Numerical modellers have developed a variety of Non-Newtonian algorithms to simulate each of these processes, and therefore, it can be difficult to understand the assumptions and limitations in any given model or replicate work. This diversity in the processes and approach to non-Newtonian simulations makes a modular computation library approach advantageous. A computational library consolidates the algorithms for each process and discriminates between these processes and algorithms with quantitative non-dimensional thresholds. This work presents a flexible numerical library framework (DebrisLib) to simulate large15 scale, post-wildfire, non-Newtonian geophysical flows using both kinematic wave and shallow-water models. DebrisLib is derived from a variety of non-Newtonian closure approaches that predict a range of non-Newtonian flow conditions. It is a modular code designed to operate with any Newtonian, shallow-water parent code architecture. This paper presents the non-Newtonian model framework and demonstrates its effectiveness by calling it from two very different modelling frameworks developed by the U.S. Army Corp of Engineers (USACE), specifically, within the one-dimensional and two-dimensional Hydrologic Engineering Centre 20 River Analysis System (HEC-RAS) and two-dimensional Adaptive Hydraulics (AdH) numerical models. The development and linkage-architecture were verified and validated using two non-Newtonian flume experiments selected to represent a range of nonNewtonian flow conditions (i.e., hyperconcentrated flow, mudflow, debris flow) commonly associated with post-wildfire flooding.

around the world (Rowe et al., 1954;Lane et al., 2006;Shin, 2010;Shakesby, 2011;Moody et al., 2013). In the years following a wildfire, ecotone shifts, gully formation, and channel incision alter the hydrologic system response, resulting in dramatic, and sometimes prolonged changes to downstream systems. Wildfires impact hydrology by removing the rainfall interception canopy, generating ash and charred material, reducing organic binding materials in soils, increasing hydrophobic soils, and modifying the intrinsic and hydraulic properties of soils and sediments (Certini, 2005;Moody et al., 2009;Ebel et al., 2012). Lahars and mine 35 tailing dams also threaten communities and infrastructure around the world and, like post-wildfire flows, these events diverge from classic, Newtonian, shallow-water flow assumptions.
Researchers have developed a range of numerical models to simulate non-Newtonian flows for hazard assessment, flood risk evaluation, and mitigation design (O'Brien et al., 1993;Iverson and Denlinger, 2001;Hungr et al., 2005). These physics-based numerical models simulate advection with constitutive laws of fluid mechanics in one-and two-dimensions. Savage and Hutter 40 (1989) provided the foundation for predicting non-Newtonian flow using Saint-Venant based shallow-water equations which have been commonly applied for numerical modelling of rapid mass movements over irregular geometry (e.g., Chen and Lee, 2000;Iverson and Denlinger, 2001;Pudasaini et al., 2005;Naef et al., 2006;Martinez et al., 2007;Pastor et al., 2009;Luna et al., 2012).
There is an increasing need to develop a post-wildfire modelling framework and enhance non-Newtonian numerical modelling 45 capabilities within hydrologic and hydraulic engineering models, for decision making in emergency management and flood risk management operations. However, with the diversity of post-wildfire processes (floods, mudflows, debris flows, etc.) and multiple algorithms available to simulate each process, a formal, modular, library framework that consolidates these algorithms and provides them to multiple hydrodynamic codes will make non-Newtonian implementation more transparent, repeatable, and robust as demonstrated here. The non-Newtonian algorithm library (DebrisLib) provides a flexible, modular, non-Newtonian, numerical-50 modelling framework that developers can call from any one-dimensional or two-dimensional shallow-water-based hydraulic and hydrologic model. Consolidating multiple non-Newtonian closures associated with the continuum of geophysical flow processes in a modular library and sharing it between hydrodynamic codes has three main advantages: 1. The library consolidates diverse literature, making a wide suite of algorithms available to each linked model, avoiding duplication of effort. 55 2. The library makes non-Newtonian modelling more transparent by standardizing algorithm implementation.
3. The library leverages the validation and verification activities of multiple development communities, accelerating the quality control of the code. https://doi.org/10.5194/hess-2020-509 Preprint. Discussion started: 9 November 2020 c Author(s) 2020. CC BY 4.0 License. This paper documents the numerical library development, enhancements, and linkage architecture necessary for predicting postwildfire non-Newtonian flood events across different numerical modelling libraries, specifically one-and two-dimensional USACE 60 numerical models. Subsequent papers will look at application of this work to real-world post-wildfire conditions.

Background
Wildfires can significantly shock or perturb semi-arid and arid environments (Pierce et al., 2004;Orem and Pelletier, 2016). Most western U.S. watersheds require decades to recover. These events increase flow and sediment loads, creating immediate and longterm management concerns for federal, state, and local agencies. These management challenges include predicting the post-wildfire 65 flood response of a particular wildfire-precipitation mix. System responses can range from minimal impact, to increased clear water flows, to more destructive ash flows, mud floods, and debris flows. Mudflows, debris flows, and ash flows are unsteady gravity-driven events that involve complex mixtures of sediment, water, and entrained material (i.e., organics, woody debris, unconsolidated substrate). These non-Newtonian flows have high sediment concentrations and can carry large boulders, trees, and uprooted infrastructure (e.g., upstream structures or bridge debris). These flows are commonly modelled using shallow-water 70 equations as either 1) single phase, 2) two-phase, or 3) mixture theory, using non-Newtonian closure approaches and approximations (e.g., Iverson, 1997;Jin and Fread, 1997;. The grain-size distribution, sediment concentration, and flow stress state (as a function of slope) determine which geophysical flow particular wildfire-precipitation events generate.
In addition to the process complexity, modellers face algorithmic diversity challenges. Each of these processes have multiple modeling approaches. Researchers have developed a wide range of algorithms to help apply these rheological and geotechnical 75 modeling libraries to different classifications of geophysical flows. The methods that simulate these post-wildfire flood events are commonly grouped into four main categories: 1) Linear (Bingham 1922) and non-linear viscoplastic models (Jin and Fread 1997;. These approaches use rheological modes as heuristics for cohesive and/or viscous stresses in the fluid and are commonly used to describe the rheology of laminar mud flows. 80 2) Dispersive-turbulent stress models (O'Brien et al. 1993). This approach adds a turbulent stress term to the dispersive model to describe the inter-particle mechanics within a clay, silt, and organic matrix. This model is commonly applied to sediment mixtures containing cohesive sediment in quantities greater than 10% by volume.

100
The research team developed a non-Newtonian library (DebrisLib) to simulate the continuum of post-wildfire flood responses along these concentration and grain size gradients, based on rheology and Coulomb theory, to operate independently of the shallowwater parent codes. The non-Newtonian library assembles closure approaches that different shallow-water codes can leverage and https://doi.org/10.5194/hess-2020-509 Preprint. Discussion started: 9 November 2020 c Author(s) 2020. CC BY 4.0 License. determines the flow conditions using non-dimensional threshold parameters. In addition, the library estimates a non-Newtonian drag coefficient (Happel and Brenner, 1965;Julien 2010), addresses hindered settling (Richardson and Zaki, 1954;Baldock et al., 105 2004;Cuthbertson et al., 2008), computes relative buoyancy terms, and accounts for increased viscosities and mass density (Einstein and Chien, 1955;Dabak and Yucel, 1986). The diversity of the processes and sub-diversity of algorithms for each process makes a shared library approach particularly useful for non-Newtonian model development. Future versions of DebrisLib will address both sediment erosion, transport, and deposition.

110
Applying non-Newtonian transport in a shallow-water model requires computes internal losses by adding a slope term to the friction slope (Sf) in the momentum equation and bulking (or increasing) the flow to account for the volume of the solids. The shallow-water friction slope (Sf) accounts for the forces acting against the flow at the fluid boundary (e.g., the channel) while the "mud and debris slope" (SMD) accounts for the internal losses due to the viscosity, turbulence, and/or dispersion within the fluid.
The library computes a non-Newtonian shear stress based on the flow classification (e.g., mudflow, debris flow, etc.) and the 115 appropriate rheological approach (i.e., stress-strain model). Then the library integrates the viscous, turbulent, and dispersion shear from the stress-strain model into the momentum equation by converting the shear to a slope and adding this slope to the friction slope. Additionally, because these flows can include between 5% and 70% solids by volume, a fixed bed model must augment the flow volume to account for the impact of the sediment on the mass and depth of the flow.

Shallow-Water Numerical Models
The shallow-water flow equations solve the continuity and momentum equations simultaneously to compute water stage and velocity. The frictional forces between the fluid and the solid boundary are the primary resisting forces in the standard Newtonian clear-water hydraulic equations. In one dimension, the conservation of mass and momentum are the Saint-Venant equations: 125 In two dimensions the shallow-water equations vertically integrate the mass and momentum equations under the assumptions of incompressible flow and hydrostatic pressure. Assuming negligible free surface shear and pressure variations at the free surface, the two-dimensional shallow-water equations are: The Reynolds stresses are determined using the Boussinesq approach to the gradient in the mean current: The depth-averaged Shallow Water Equations model in HEC-RAS solves volume and momentum conservation equations and 160 includes temporal and spatial accelerations as well as horizontal mixing while the diffusive wave equation model ignores these processes but is therefore simpler and more computationally efficient. The 2D volume conservation of the water-solid mixture is given by: where is the flow surface elevation, is time, ℎ is the water depth, is the velocity vector, and is a source or sink term, to 165 account for external and internal fluxes. The depth-averaged momentum conservation equations may be written as (Hergarten and Robl, 2015): in which is the gravitational acceleration, is a turbulent eddy viscosity, = ( , ) is the total basal stress, is the water-170 solid mixture density, is the hydraulic radius, , is the water surface slope, and is the inclination angle of the current velocity direction. In the above equations, the second term on the right-hand-side represents the horizontal mixing due to turbulence and also in the case of a debris flow, horizontal mixing due to particle collisions. Utilizing the conservative form of the mixing terms is essential for accurate momentum conservation. The bottom friction coefficient is computed utilizing the Manning's roughness coefficient as = + , where is the turbulent stress and is the mud and debris 175 stress which includes all non-Newtonian stresses. The turbulence bottom shear stress is computed as a function of the Manning's roughness coefficient. The mud and debris stress are described in detail in the section "Rheological Models". It is noted that when the non-Newtonian stress is equal to zero and the cosine functions (slope corrections) are removed, the above 2D shallowwater equations reduce to the clear-water equations utilized in HEC-RAS.

Non-Newtonian Shallow-Water Closure 180
Mud and debris flows generate additional resisting forces. Increasing the solid content increases the viscosity of non-Newtonian flows, generating internal resisting forces within the fluid. At higher concentrations, particularly with coarse particles, particle https://doi.org/10.5194/hess-2020-509 Preprint. Discussion started: 9 November 2020 c Author(s) 2020. CC BY 4.0 License. (1) cohesion between particles, (2) internal friction between fluid and sediment particles, (3) turbulence, and (4) inertial impact between particles. The quadratic model separates the stress-strain relationships into these four, additive components, such that the shear stress is:   Takahashi (1980) identified experimentally that the Bagnold impact coefficient ( ) ranges between 0.35 and 0.5 and that is significantly larger than the value recommended from Bagnold (1954. Iverson (1997) where = the Von Karmen coefficient (≅ 0.41) = the proportional distance from the boundary (bed). 230 This quadratic model combines linear and non-linear rheological modes to compute internal shear. Rheological models do not perform as well as the mixtures get more clastic (i.e., high concentrations of coarse particles). DebrisLib simulates clastic debris flows with a Coulomb approximation based on the Johnson and Rodine (1984) Coulomb viscous model (Naef et al., 2006). This approach replaces the Bingham yield strength ( ) a geotechnical, Coulomb yield stress defined as, = Coulomb friction angle (°) with typically ranges between 30° and 40° (Iverson, 1997;McArdell et al., 2007). https://doi.org/10.5194/hess-2020-509 Preprint. Discussion started: 9 November 2020 c Author(s) 2020. CC BY 4.0 License.
The Coulomb friction angle is a function of the individual grain friction angle and the packing geometry of the particles along the 240 failure plane.

Numerical Model Discretization
The

Verification and Validation Datasets
The verification and validation process involved simulation of multiple flume experiment selected to represent the continuum of non-Newtonian flow behaviour under both steady and unsteady conditions. This includes the Jeyepalan (1981) and Hungr (1995) Iverson et al. (1992 and. Initial verification of DebrisLib was conducted using the 1D analytical solutions from an idealized 30-meter-high tailings dam subjected https://doi.org/10.5194/hess-2020-509 Preprint. Discussion started: 9 November 2020 c Author(s) 2020. CC BY 4.0 License.
to instant liquefaction failure from Jeyapalan (1981) and Hungr (1995). The primary verification variables included failure runout distance, velocity, depth, and cessation of flow. The initial validation flume experiments are from Haldenwang et al. (2006) who sand, and 7% mud particles (Iverson et al., 2010). The HEC-RAS grid for the USGS flume experiment has a constant grid resolution of 0.15 and 0.2 in the downslope and transverse directions respectively with a total of approximately 13,350 cells.

Results 285
The model library and the linkage-architecture were evaluated by simulating two non-Newtonian verification data sets by calling the library from two different hydrodynamic codes. To evaluate the model library a selection of flume experiments was used to demonstrate the library's flexibility using both small-and large-scale flows of various grain size distributions and sediment concentrations. These flume experiments were selected to demonstrate the flexibility of the model library and its linkage architecture with USACE shallow-water models. The data sets simulated include: equilibrium, hyperconcentrated flow (Cv=10%) 290 from Haldenwang et al. (2006) experiments and an analytical, unsteady non-Newtonian dam break with a higher concentration (Cv=55%), from Jeyapalan 1981 andHungr (1995) running out until flow stops or yields. Input conditions and calibration data for Hungr 1995, Haldenwang et al., 2006, and Iverson et al., 2010 experiments are provided in Table 1. It is important to note that this paper only presents our general findings regarding USACE modelling capabilities. Subsequent publications will focus on the additional details for the provided flume experiments. The 2D HEC-RAS and AdH results for the Jeyapalan 1981 andHungr (1995) dam breach simulations are included in Figure 2 (left and right respectively) representing conditions of dynamic unsteady conditions following an impoundment failure. The initial profile represents the 30 m high pre-dam break conditions and the final profile is the analytical solution for the failure runout profile when the flow cesses. This analytical 300 dam removal computes "run out length," which is the distance required for this high-energy mass flow to come to rest (i.e., the shear stress drops below the yield strength). Both numerical simulations compute run out lengths very similar to Hungr's (1995) result, and have similar final depth profiles. Both HEC-RAS and AdH sufficiently replicated the impoundment failure with to include unsteady runout, flow depth and velocity. These results were consistent with other publications based on this analytical solution (e.g. Naef et al., 2005). Simulations were conducted using Bingham, Herschel-Bulkley, and O'Brien quadratic model with 305 only Bingham rheology conditions presented with all cases exhibiting similar flow behaviour and results. representation of the flow and depths for the 5-degree flume conditions. This is consistent for both HEC-RAS and AdH simulation.
The differences between HEC-RAS and AdH results are likely due to a range of factors, to include how each model accounts for initial conditions, model discretization, grid size, and time steps size. All rheology closure models simulated under predict depthflow curve at flow ranges smaller than approximately 7.5 L/s likely a result of several factors, to include the strain rate 320 approximation of 3 ℎ � approximation, change in velocity profile, potential sediment deposition, or a combination. This will be investigated more thoroughly in subsequent research. Interestingly all models did better at predicting transitional and turbulent conditions versus laminar ((as shown with the red line in Figure 3). Only on previous study was found that used the Haldenwang data set for comparison with a shallow-water model (Perez, 2017) but our results are consistent with those findings. Additional model research is ongoing using this dataset aimed at betters quantifying turbulent behaviour at lower sediment concentrations 325 analogous to hyperconcentrated flows. This will also help better understand the deviations between model and flume results for flow conditions less that about 5 L/s).

330
AdH and HEC-RAS are validated with a representative experiment with a mixture of sand, gravel, and mud (silts and clays) and a section of bumpy concrete tiles. The comparison of AdH and HEC-RAS water levels compared to the measurements from one of the USGS experiments is presented in Figure 4. Both AdH and HEC-RAS water surface time series compared reasonably well with pressure measurements. Differences in AdH and HEC-RAS results due to differences in the model formulation and numerical 335 assumptions and solution methods. This demonstrates that the optimal rheological parameters will likely be different for different models. As shown in Figure 4 AdH more closely replicated the USGS flume results compared to the HEC-RAS results. Both models tend to do better at the upstream locations specifically at x = 2 m and x = 32 m. Model results continue to deviate from flume conditions slightly for x = 66 and x = 90 with more pronounced deviation seen in the HEC-RAS results. These differences are the result of the flume slope being approximately 31 degrees which violates the shallow-water assumptions of most depth 340 averaged models, including these. Surprisingly they still perform reasonable considering this limitation. https://doi.org/10.5194/hess-2020-509 Preprint. Discussion started: 9 November 2020 c Author(s) 2020. CC BY 4.0 License.

Discussion 345
The post-wildfire non-Newtonian model library DebrisLib expands upon the quadratic model (O'Brien et al., 1993) to predict a range of post-wildfire sedimentation behaviour across multiple numerical model platforms. The flexible non-Newtonian library allows for prediction of a range of post-wildfire flood conditions demonstrated using a two shallow-water numerical models. We introduce a modification to existing non-Newtonian taxonomy and classification describing common non-Newtonian post-wildfire flows and demonstrated using both 1D and 2D shallow-water numerical models for a range of non-Newtonian flow conditions. To 350 date, existing non-Newtonian models predict a limited range of flow conditions within a single modelling platform thus multiple models must be stitched together to address the full range of flow conditions. The process and approach diversity of the non-https://doi.org/10.5194/hess-2020-509 Preprint. Discussion started: 9 November 2020 c Author(s) 2020. CC BY 4.0 License. Newtonian literature in particular make a modular computation library-based framework advantageous by consolidating the algorithms for each process. With a few common variables like depth, velocity, mixture density, and slope the model library can be linked with most existing shallow-water Newtonian numerical models, regardless of dimension, library, or solution method, to 355 predict non-Newtonian flows. Currently the model library is under development with planned released scheduled soon. Simulating a range of flow conditions across multiple non-Newtonian transport models will increase access to post-wildfire modelling capabilities and reduce uncertainty associated with comparisons from multiple independent model approaches. Furthermore, this research constitutes the first non-Newtonian simulations using USACE shallow-water models.
The flexibility will allow broader access to the non-Newtonian modelling capabilities for practitioners, managers, and researchers 360 to improve understanding of the model limitations and operational modelling for post-wildfire emergency management and flood risk management. During development of the model library the team identified some key disadvantages and advantages associate with library-based frameworks in research and development. An overview of the advantages and challenges associated with development of a flexible library-based model framework are provided below in Table 3. Table 2. Advantages and challenges associated with development of Library-based frameworks.

Conclusion 375
High-intensity wildfires remove vegetation, alter soils, modify subsurface root structures, purge organic soil, and create widespread hydrophobic soils. This increases runoff, erosion, and sediment transport which results in destructive flooding. Wildfire effects can generate gravity-driven surface runoff and erosion events that involve complex mixtures of water, ash, sediment, and entrained debris (i.e., destroyed upstream infrastructure, woody debris, and very large sediment clasts). Several other geophysical flows, including lahars and mine tailing dam failures have similar non-Newtonian dynamics and require these closures to model them 380 effectively. Hydrologic and hydraulic models that assess wildfire impacts on flood risk management need to develop and improve non-Newtonian numerical modelling capabilities. This work consolidates the state-of-the-practice in a modular non-Newtonian and sediment transport library. Then demonstrates that it can easily connect to multiple modelling platforms with different dimensionality, interfaces, and mesh requirements, and accurately predict commonly post-wildfire non-Newtonian flow conditions. This modular library improves and enhances prediction capabilities to assist with planning, management, and mitigation 385 in post-wildfire environments using practical science-based approaches and smart integrated numerical approaches.