DebrisInterMixing-2 . 3 : a finite volume solver for three-dimensional debris-flow simulations with two calibration parameters – Part 1 : Model description

Here, we present a three-dimensional fluid dynamic solver that simulates debris flows as a mixture of two fluids (a Coulomb viscoplastic model of the gravel mixed with a Herschel–Bulkley representation of the fine material suspension) in combination with an additional unmixed phase representing the air and the free surface. We link all rheological parameters to the material composition, i.e., to water content, clay content, and mineral composition, content of sand and gravel, and the gravel’s friction angle; the user must specify only two free model parameters. The volume-of-fluid (VoF) approach is used to combine the mixed phase and the air phase into a single cellaveraged Navier–Stokes equation for incompressible flow, based on code adapted from standard solvers of the opensource CFD software OpenFOAM. This effectively singlephase mixture VoF method saves computational costs compared to the more sophisticated drag-force-based multiphase models. Thus, complex three-dimensional flow structures can be simulated while accounting for the pressureand shear-rate-dependent rheology.


Introduction
Debris flows typically occur in steep mountain channels. They are characterized by unsteady flows of water together with different contents of clay, silt, sand, gravel, and larger particles, resulting in a dense and often rapidly moving mixture mass. They are often triggered by heavy rainfall and can cause massive damage (Petley et al., 2007;Hilker et al., 2009). Their importance has increased due to extensive settlement in mountainous regions 20 and also due to their sensitivity to climate change (Guthrie et al., 2010). Their damage potential is not limited to direct impact; severe damage can also be caused by debris flows blocking channels and thus inducing over-topping of the banks by subsequent flows (Tang et al., 2011).
Modeling debris flows is a central part of debris-flow research, because measuring the 25 detailed processes in debris-flow experiments or in the field is challenging. It is still uncertain 2 how laboratory tests can be scaled to represent real flow events, and the inhomogeneous mixture of gravel and fine material brings about interactions of granular flow and viscous forces such as drag and pore-pressure that are difficult to track with the present measurement techniques at reasonable cost. As a consequence, the rheological behavior of debris flow material remains incompletely understood. 5 Typically, current numerical modeling approaches cannot predict run-out distances or impact pressures of debris flows without a parameter calibration that is based on simulations of previous well-documented events that happened at the same site. This clearly represents a challenge in practical applications, because reliable calibration data are rarely available. Due to the complex physics of debris flows, real flows can only be accurately described 10 by dynamical models that include strong phase interactions between granular and viscous fluid phases with several physical parameters (Pudasaini, 2012). From a practical application point of view and guided by considerations of computational efficiency, here we neglect such phase interactions and restrict ourselves to an effectively single-phase mixture flow simulation. All currently applied debris flow models that use a two-phase description of the 15 debris-flow material are depth-averaged or 2D. Three dimensional debris-flow models with momentum exchange between phases have, up to now, been limited to academic cases due to their high numerical costs. Depth-averaged approximations are most applicable to flows over smooth basal surfaces with gradual changes in slope. These approximations are less valid when topography changes are abrupt, such as close to flow obstacles. They are 20 also less applicable when the flow exhibits strong gradients in accelerations, such as during flow initiation and deposition, or in strongly converging and diverging flows. For such cases, we need a physically complete description of the flow dynamics without depth-averaging (Domnik and Pudasaini, 2012). The essential physics of debris flows can be better retained by developing a full-dimensional flow model and then directly solving the model equations 25 without reducing their dimension. The currently available models also contain many parameters that must be estimated based on measurements, or fitted to site-specific field data, limiting their applicability to real-world problems.
Here, we provide a largely simplified but effective solution linking the rheological model of the debris-flow material to field values such as grain size distribution and water content. The approach aims to link the knowledge of field experts for estimating the release volume and material composition with recent advances that account for complex flow phenomena using three-dimensional computational fluid dynamics with reasonable computational costs. 5 The parameters of the two resulting rheology models for the two mixing fluids are linked to material properties such that the model setup can be based on material samples collected from the field, yielding a model that has two free parameters for calibration. One mixture component represents the suspension of finer particles with water (also simply called slurry in this paper) and a second component accounts for the pressure-dependent flow behavior 10 of gravel. The two components result in a debris bulk mixture with contributions of the two different rheology models, weighted by the corresponding component concentrations. A third gas component is kept unmixed to model the free surface. The focus is on the flow and deposition process and the release body needs to be user-defined. The debris flow material can be considered as a combination of a granular component and an interstitial fluid composed of the fine material suspension. The latter was successfully modeled in the past as a shear-rate dependent Herschel-Bulkley fluid (Coussot et al., 1998). Because pressure and shear drive the energy 5 dissipation of particle-to-particle contacts, the shear rate substantially influences the energy dissipation within the granular phase. While the two-phase models of Iverson and Denlinger (2001) and Pitman and Le (2005) treated the granular phase as a shear-rate independent Mohr-Coulomb plastic material, dry granular material has been successfully modeled as a viscoplastic fluid by Ancey (2007), Forterre and Pouliquen (2008), 10 Balmforth and Frigaard (2007) and Jop et al. (2006). We follow the suggestions given by Pudasaini (2012) to account for the non-Newtonian behavior of the fluid and the shear-and pressure-dependent Coulomb-viscoplastic behavior of the granular phase, as proposed by Domnik et al. (2013). Several modeling approaches to account for the two-phase nature of debris flows used depth-averaged mass and momentum balance equation for each 15 phase coupled by drag models (eg. Bozhinskiy and Nazarov (2000), Pitman and Le (2005) and Bouchut et al. (2015)). Pudasaini (2012) proposed a more comprehensive two-phase mass flow model that includes general drag, buoyancy, virtual mass, and enhanced non-Newtonian viscous stress in which the solid volume fraction evolves dynamically. We apply a numerically more efficient method and treat the debris flow material as one mix-20 ture (Iverson and Denlinger (2001), Pudasaini (2005a)), with phase-averaged properties described by a single set of Navier-Stokes-type equations. The resulting reduction in numerical costs allows us to model the three-dimensional momentum transfer in the fluid as well as the free-surface flow over complex terrain and obstacles.
We assume that the velocity of the gravel is the same as the velocity of the fluid. This 25 assumption is motivated from application and would not be valid for general debris flows with interstitial fluid of small viscosity (i.e., a slurry with low concentration of fine material and large water content). The assumption of equal velocity of gravel and interstitial fluid in 6 one cell leads to a constant composition of the mixture by means of phase concentrations over the entire flow process. This basic model can be seen as a counterpart to the mixture model of Iverson and Denlinger (2001), extended by resolving the three-dimensional flow structure in combination with a pressure-and shear-rate-dependent rheology linked to the material composition. The basic model presented here focuses on the role of pressure- 5 and shear-rate-dependent flow behavior of the gravel, in combination with the shear-ratedependent rheology of the slurry. We base our model concept on the well-established finite volume solver interFoam, which is designed for incompressible two-phase flow simulations of immiscible fluids (Deshpande et al., 2012). A standard extension named interMixingFoam introduces two 10 mixing phases without momentum exchange, coupled to a third unmixed phase by surface tension. The present method limits the physics of flow in terms of the interactions between phases such as drag, buoyancy, virtual mass, non-Newtonian viscous stress and evolving volume fraction of the solid phase (Pitman and Le (2005); Pudasaini (2012)). Numerical costs are kept reasonable due to the Volume-of-Fluid (VoF) 15 method (Hirt and Nichols, 1981), which solves only one Navier-Stokes equation system for all phases. The viscosity and density of each grid cell is calculated as a concentrationweighted average of the viscosities and densities of the phases that are present in the cell. Between the two mixing phases of gravel and slurry, the interaction reduces to this averaging of density and viscosity. In this way, the coupling between driving forces, topography 20 and three-dimensional flow-dependent internal friction can be addressed for each phase separately, accounting for the fundamental differences in flow mechanisms of granular and viscoplastic fluid flow that arise from the presence or absence of Coulomb friction (Fig. 1). We apply linear concentration-weighted averaging of viscosities for estimating the bulk viscosity of a mixture for simplicity. Nevertheless, interacting forces such as drag between the 25 grains and the fluid do not appear in our model because we set solid velocity equal to fluid velocity, and thus from the physical point of view this is a debris bulk model. 7

Governing equations
Assuming isothermal incompressible phases without mass transfer, we separate the modeled space into a gas region denoting the air and a region of two mixed liquid components. The VoF method used here determines the volume fractions of all components in an arbitrary control volume by using an indicator function which yields a fraction field for 5 each component. The fraction field represents the probability that a component is present at a certain point in space and time (Hill, 1998). The air fraction may be defined in relation to the fraction of the mixed fluid α m as (1) and the mixed fluid α m may be defined as the sum of the constant fractions of the mixing 10 phases α 2 and α 3 : The debris bulk motion is defined by the continuity equation together with the transport and momentum equations: and 8 where t denotes time and U represents the debris bulk velocity field, T is the deviatoric viscous stress tensor for the mixture, ρ is the phase-averaged bulk density, p denotes pressure and f stands for body forces per unit mass. We assume incompressible material and all fractions are convected with the same bulk velocity. So, differences in phase velocities, and thus the interaction forces such as drag between the grains and the fluid are neglected. 5 An efficient technique of the VoF method convects the fraction field α m as an invariant with the divergence-free flow field U that is known from previous time steps: where U c is an artificial inter-facial compression velocity acting perpendicular to the interface between the gas region and the mixed liquid phases. The method allows a 10 reconstruction of the free surface with high accuracy if the grid resolution is sufficient (Berberović et al., 2009;Hoang et al., 2012;Deshpande et al., 2012;Hänsch et al., 2013). The details about the interface compression technique, the related discretization and numerical schemes to solve Eq. 6 are given in Deshpande et al. (2012). However, to allow diffusion between the mixing constituents of the slurry α 2 and the gravel α 3 in case of 15 initially unequally distributed concentrations, our modified version of the interMixingFoam solver applies Eq. 6 separately to each mixing component including diffusion: where i = 2, 3 denote the slurry and gravel constituents and D dif f is the diffusion constant that is set to a negligible small value within the scope of this paper. 20 The discrete form of Eq. 7 is derived by integrating over the volume V of a finite cell of a grid-discretization of the simulated space, which is done in the finite volume method by applying the Gauss Theorem over the cell faces. The advective phase fluxes φ 1..3 are obtained by interpolating the cell values of α 1 , α 2 and α 3 to the cell surfaces and by multiplying them with the flux φ through the surface, which is known from the current velocity field. To Equation 8 is the implemented scalar transport equation solved for each constituent using first-order Euler schemes for the time derivative terms, as has been recommended for liquid column breakout simulations (Hänsch et al., 2013). 10 After solving the scalar transport equations, the complete mass flux φ ρ is constructed from the updated volumetric fraction concentrations: where ρ 1..3 denote the constant densities of the corresponding phases and φ 1..3 are the corresponding fluxes. The complete mass flux φ ρ is used in the implementation to describe 15 the second term of the momentum equation (Eq. 14) as described in appendix A. Fig. 2 illustrates how the phase volume distributions α 1 (air), α 2 (slurry) and α 3 (gravel) are used to derive cell-averaged properties of the continuum. The conservation of mass and momentum is averaged with respect to the phase fraction α of each phase. The density field is defined as where ρ i denotes density of phase i and the density is assumed to be constant.
The deviatoric viscous stress tensor T is defined based on the mean strain rate tensor D that denotes the symmetric part of the velocity gradient tensor derived from the phaseaveraged flow field: Equation 12 was derived accounting for Eq. 3 and µ is the phase-averaged dynamic viscosity, which is simplified in analogy to Eq. 10 as the concentration-weighted average of the corresponding phase viscosities: 10 With the continuity equation 3 the term ∇ · T in Eq. 5 is written as ∇ · (µ∇U ) + ∇U · ∇µ to ease discretization. The body forces f in the momentum equation account for gravity and for the effects of surface tension. The surface tension at the interface between the debris flow mixture and air is modeled as a force per unit volume by applying a surface tension coefficient σ. Although the surface tension can be considered to have a minor influence on 15 debris flow behavior, it allows an adequate reproduction of surface flow patterns observed in laboratory scale experiments used for validation (v. Boetticher et al., 2015). The momentum conservation including gravitational acceleration g and surface tension is defined in our model as: 20 11 where κ denotes the local interfacial curvature and x stands for position. The modified pressure p d is employed in the solver to overcome some difficulties with boundary conditions in flow simulations with density gradients. In case the free surface lies within an inclined wall forming a no-slip boundary condition, the normal component of the pressure gradient must be different for the gas phase and the mixture due to the hydrostatic component ρg. It is 5 common to introduce a modified pressure p d related to the pressure p by The gradient of the modified pressure includes the static pressure gradient and contributions that arise from the density gradient as well as a body force due to gravity (Berberović et al., 2009). 10 Together with the continuity equation 3, Eq. 14 allows us to calculate the pressure and gravity driven velocities. The corresponding discretization and solution procedure with the PISO (Pressure-Implicit with Splitting of Operators (Issa, 1986)) algorithm is provided in appendix A. The set of equations governing the flow in our model are Eq. 8 to Eq. 15 together with the continuity equation 3. In the following section we present the rheology 15 models that define the viscosity components for Eq. 13.

Rheology model for the fine sediment suspension
The viscosity of the gas phase, µ 1 is chosen constant. The introduction of two mixing phases is necessary to distinguish between the pressure-dependent flow behavior of gravel and the shear-thinning viscosity of the suspension of finer particles with water. The rheology 20 of mixtures of water with clay and sand can be described by the Herschel-Bulkley rheology law (Coussot et al., 1998), which defines the shear stress in the fluid as: where τ y is a yield stress below which the fluid acts like a solid, k is a consistency index for the viscosity of the sheared material,γ is the shear rate and n defines the shear-thinning (n < 1) or shear-thickening (n > 1) behavior. In OpenFOAM, the shear rate is derived in 3D from the strain rate tensor D: 5 The shear rate is based on the strain rate tensor to exclude the rotation velocity tensor that does not contribute to the deformation of the fluid body. The model can be rewritten as a generalized Newtonian fluid model to define the shear-rate-dependent effective kinematic viscosity of the slurry phase as: 10 if the viscosity is below an upper limit µ 0 and if the viscosity is higher, to ensure numerical stability.
With n = 1 the model simplifies to the Bingham rheology model that has been widely used to describe debris-flow behavior in the past (Ancey, 2007). It may be reasonable to 15 imagine the rheology parameters to be dependent on the state of the flow. However, even with the implicit assumption that the coefficients are a property of the material and not of the state of the flow, the Herschel-Bulkley rheology law has been rarely applied in debrisflow modeling due to three rheology parameters. We avoid this problem by assuming some rheology parameters to be defined by measurable material properties as described below. 20 13

Determination of rheology model parameters based on material properties
Results from recent publications allow the reduction of the number of free Herschel-Bulkley parameters to two. If the coarser grain fraction is assumed to be in the gravel phase, the Herschel-Bulkley parameters for the finer material can be linked to material properties that can be measured using simple standard geotechnical tests. According to Coussot et al. 5 (1998), the exponent n can be assumed constant at 1/3, and k can be roughly estimated as bτ y , where the constant b = 0.3s −n for mixtures with maximum grainsizes < 0.4 mm (Coussot et al., 1998). An approach for estimating the yield stress τ y based on water content, clay fraction and composition, and the solid concentration of the entire debris flow material was proposed by Yu et al. (2013) as: where C (a constant) is the volumetric solid concentration of the mixture (the volume of all solid particles including fine material relative to the volume of the debris flow material including water), P 1 = 0.7P 0 when P 0 > 0.27 and P 1 = P 0 if P 0 <= 0.27, and 15 where the subscript of C refers to the volumetric concentration (relative to the total volume of all solid particles and water) of the corresponding mineral. The discontinuity of P 1 at a modified clay concentration of P 0 = 0.27 is a coarse adjustment to a more-or-less sudden change observed in the experimental behavior.
For C < 0.47, τ 0 is equal to τ 00 and otherwise τ 0 can be calculated by 20 τ 0 = τ 00 e 5(C−0.47) 14 where τ 00 is the remaining free parameter which we use to account for the grid size dependency of the shear rate (Yu et al., 2013). The change in the definition of τ 0 at solid concentrations exceeding 0.47 accounts for a threshold between weak coherent and medium coherent debris flow material. Yu et al. (2013) compared this method of estimating the yield stress τ y to experimental 5 results they obtained from a set of 514 flume experiments with mixtures of water and clay with fine and coarse sand and less than 5 % gravel. They determined the yield stress by releasing the material mixture from a reservoir into an inclined channel of 0.2 m width and by increasing the inclination slightly until remobilization occurred after the material came to rest. The experimental yield stress τ y−exp was then simplified as: where ρ exp is the density of the applied mixture, h is the maximum accumulation thickness of the deposit, and β is the slope inclination. In addition, they compared the calculated yield stress of Eq. 20 with experimental yield stresses reported by a number of authors including Coussot et al. (1998) and Ancey and Jorrot (2001). Ancey and Jorrot (2001)  the numbers for calculating P 0 may be necessary to account for the activity of other clays. The remaining uncertainties concern the assumptions about the value of n, and that k can be defined in such simple dependency to τ y in the presence of coarser sand. Experiments seem to confirm that n increases in presence of coarser material (Imran et al., 2001), but further research is needed to quantify this effect. Remaitre et al. (2005) found n to vary 25 from 0.27 to 0.36. Schatzmann et al. (2003) used n = 0.33 to reproduce measured curves obtained with a mixture of 27.5 volumetric percent slurry with 30 % gravel where gravel grainsizes ranged from 3 to 10 mm. He applied n = 0.5 to fit the Herschel-Bulkley model to the experiment with 22.5 % slurry and 30 % gravel. Based on the laboratory-scale experiments that are presented in v. Boetticher et al. (2015) we have chosen n = 0.34 to obtain the best fit for the simulation of large-scale debris-flow experiments. However, other debris flow mixtures may demand a recalibration of n because natural debris flow mixtures cover 5 a wide spectrum from shear thinning to shear thickening behavior. Thus we consider n as the second model calibration parameter. The model calibration is constricted to τ 00 which accounts for the grid resolution sensitivity of the shear-rate and n, which allows to adjust the model to shear thinning or shear thickening mixtures. 10 With a novel model, Domnik and Pudasaini (2012) showed that introducing a Coulombviscoplastic sliding law to replace a simple no-slip boundary condition completely changes the flow dynamics and the depth variations in velocity and dynamic pressure. They show that with a Coulomb-viscoplastic sliding law, the observable shearing mainly occurs close to the sliding surface, in agreement with observations but in contrast to the no-slip bound- 15 ary condition. In our model, during acceleration and high-speed flow, the shear-thinning behavior of both the fluid and the granular phase dominate the viscosity. However, the representation of gravel by a pressure and rate-dependent Coulomb-viscoplastic rheology becomes important as soon as the material experiences high pressures, accompanied by reduction in shear due to decelerations. This may be caused by channel slope reductions, 20 diverging or converging flows, and interactions with obstacles (see, e.g., Pudasaini (2005b); Domnik et al. (2013)). Flows of granular material could be modeled as viscoplastic fluids (Ancey, 2007;Forterre and Pouliquen, 2008;Balmforth and Frigaard, 2007;Jop et al., 2006) as cited by Domnik and Pudasaini (2012). Based on Ishii (1975), the granular Cauchy stress tensor T s can be written as: 25 T s = −pI + 2µ s D,

Representation of gravel by a Coulomb-viscoplastic rheology
16 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | where pI is the pressure times the identity matrix and µ s is the corresponding dynamic viscosity, which was modeled by Domnik and Pudasaini (2012) as: where µ min is a minimal dynamic viscosity, τ 0s is a yield stress, and ||D|| is the norm of the strain-rate tensor defined by the authors as: In Eq. 25, m y is a numerical parameter with units of s which we will keep constant. Domnik et al. (2013) derived the yield stress as a pressure-dependent Coulomb friction, p · sin(δ) where δ is the internal friction angle: 10 Here, this Coulomb-viscoplastic rheology model is used to describe the gravel phase. The pressure-and shear-dependent viscosity is calculated in every cell with the corresponding local pressure p and strain-rate tensor D derived from the phase-averaged flow field.

Gravel phase properties
The Coulomb-viscoplastic rheology law developed by Domnik et al. (2013) includes two 15 parameters: one material parameter that can be measured (the friction angle δ) and a numerical parameter m y introduced for numerical reasons to influence the smooth transition between yielded and unyielded regions. For smaller values of m y , the transition is smoother. In the absence of shear, to achieve a viscosity representing a Coulomb friction equal to p · sin(δ) where p is the local pressure, m y needs to be equal to 1 s. However, the development of µ s under large pressure or strong shear is the same for both m y = 1 s and m y = 0.2 s. So, we choose m y to be constant and equal to 0.2 s for all simulations. However, we mention that parts of the nearly immobile material that experience little pressure (in general, immobile material close to the surface) show a significant reduction in viscosity 5 when m y = 0.2 s (Fig. 3). As a consequence, m y minimally affects debris flow release and flow at large scales, but material with a shallow flow depth in a run-out plane close to deposition may develop front fingering (which is dependent on, and sensitive to, the value of m y ) by allowing sudden local solidification.
For small friction angles, the modeled viscosity of the gravel phase decreases rapidly with 10 increasing shear. Larger friction angles increase the viscosity and extend the pressure dependency to larger shear rates (Fig. 4). We estimated the friction angle δ based on the maximum angle of repose in tilt-table tests of the gravel. In our laboratory experiments, we determined the friction angle in a simple adaptation of the method of Deganutti et al. (2011) by tilting a large box with loose material until a second failure of the material body occurred. 15 In analogy to the Herschel-Bulkley implementation, an upper limit for the viscosity is implemented to maintain numerical stability. Pressure-dependent viscosity in the incompressible Navier-Stokes equations causes numerical instability as soon as the eigenvalues of the symmetric part of the local velocity gradient become larger than 1/(2(δµ/δp)). Following Renardy (1986), we locally limit the viscosity to keep it below a corresponding local stability 20 limit.

Modeled interaction between granular and fluid constituents
As stated before, debris flows are multiphase flow processes where solid particles and interstitial fluid and probably gas phases interact. The model presented applies cell-averaged 25 values of the mixtures and does not allow for consideration of interactions between a gran-18 ular and a fluid phase. However, a dynamic evolution of the gravel and slurry concentration is the necessary next step. A two-way coupling to a Lagrangian particle simulation could deliver an alternative solution.

Advantages of full dimensional mass flow
Laboratory and field observations show that the rapid flow regime is characterized by rel- 5 atively uniform velocity profiles with depth, and dominant sliding at the base. By contrast, in the depositional regime and particularly during transitions from rapid flow to the deposition regime, shock-like structures form and a predominantly basal sliding flow changes into a surface boundary layer flow, which, further downstream, quickly slows and eventually settles. This transition has been observed in the laboratory with granular PIV (parti-10 cle image velocimetry) measurements for rapidly moving granular material impinging on a rigid wall perpendicular to the bed, leading to strong shearing through the flow depth (Pudasaini, 2007). The scaling of impact pressures on rigid obstacles with obstacle size is complex (Bugnion et al., 2012); the spatial distribution depends on the local flow structure (Scheidl et al., 2013) and the distribution of forces over the flow depth is non-trivial 15 and highly dynamic (Berger et al., 2011).Thus, although the rapid flow regime can be reasonably approximated by a depth-integrated dynamical model, adequately simulating the deposition regime and interactions with obstacles requires resolving the three-dimensional flow structure. Our attempts to consider the three dimensional flow structure of debris flows were motivated by the need to predict the dynamic loading of shallow landslide impacts 20 to flexible protection barriers (v. Boetticher, 2013) and the corresponding research in full dimensional model development is a topic of ongoing interest (Leonardi et al., 2016). In the model presented here we aim to account for both the three-dimensional flow structure and the corresponding influence of pressure. Domnik et al. (2013) developed a two-dimensional Coulomb-viscoplastic model and applied it for inclined channel flows of granular materials 25 from initiation to deposition. They proposed a pressure-dependent yield strength to account for the frictional nature of granular materials. The interaction of the flow with the solid boundary was modelled by a pressure-and rate-dependent Coulomb-viscoplastic sliding law. In regions where depth-averaging becomes inaccurate, like in the initiation and deposition regions, and particularly in interactions with obstacles, three-dimensional models are essential because in these regions the momentum transfer must be considered in all directions. Thus, prediction of the velocity variations with depth, and the full dynamical and internal pressure, can only be obtained through three-dimensional modeling (Domnik et al., 2013). 5 These dynamical quantities can be adequately described by a Coulomb-viscoplastic material with a pressure-dependent yield stress expressed in terms of the internal friction angle, the only material parameter in their model. The pressure dependency introduces an additional need to account for the three-dimensional flow structure. 10 Because most debris-flow models are depth-averaged and use shallow-water approximations, one could ask why a three-dimensional approach is necessary. Brodani-Minussi and deFreitas Maciel (2012) simulated dam-break experiments of a Herschel-Bulkley fluid using the VoF approach, and compared the results with published shallow-water-equation-based models. Especially for the first instant after the material release, the application of shallow-events at the same location and similar conditions. But changes in release volume or position can lead to different accelerations and corresponding changes in the automatic time step control can alter the development of rheology over time. As long as a flow stage is reached where the flow stops accelerating, the influence on the final front velocity should be negligible. Other debris flow models, which do not iteratively adjust viscosity, cannot ac-5 curately simulate impacts. Here, our model constitutes a significant improvement, since in the three-dimensional solver we presented, the viscosity bias was reduced by implementing a corrector step: taking the average between the viscosity at the beginning of the time step and the viscosity that corresponds to the velocity field at the end of the time step, the time step is solved again, leading to a better calculation of the velocity. This step can be 10 repeated, according to user specifications, to correct the viscosity several times. Although this procedure increases numerical calculation time, it clearly reduces the time-step dependency of the simulation. Some dependency on the time step is still present when modeling the collapse of material columns, but the origin of this problem is different because it occurs also for Newtonian fluids.

Effect of grid resolution on rheology
Since the shear rate influences both viscosity models, viscosity is sensitive to grid resolution, because the shear rate is averaged over the cell size. For flows over rough topography this may be less critical, but for laboratory flume experiments with thin shear bands the results may depend on grid resolution. When simulating laboratory flume experiments where 20 debris-flow material accelerated in a relatively narrow and short channel (Scheidl et al., 2013), a cell height of 1.5 mm, which is of the order of the laboratory rheometer gap, was still not fine enough to reach the limit of grid sensitivity. The free model parameter τ 00 influences the shear-rate-dependent term of the viscoplastic rheology model and can be used to adjust the simulation to the grid resolution. As long as the gravel phase and grid resolu-adjusting to different mixtures based on one calibrated test is to perform one iteration step for the yield stress of the new mixture; by calculating τ y based on the original τ 00 value from the calibration test but with the new material composition, an updated yield stress of the new mixture is determined. Raising or lowering τ 00 by the same ratio as the change from the original yield stress of the calibration test to the updated yield stress generates the final 5 τ y as it is applied to the simulation of the new mixture.
The viscosity of the granular phase is averaged over the cell faces to avoid discontinuous viscosity jumps between cells, which may cause instability due to the sensitivity of incompressible solvers to pressure-dependent viscosity. However, thin cells that allow a precise calculation of the shear gradient lead to a preferred direction of the smoothing of the granu-10 lar phase's viscosity which may introduce physically unrealistic behavior. Cell length (in the flow direction), cell width and cell height should at least be of the same order. Especially when front fingering is of interest, a grid resolution test should be carried out, ensuring that front instability is not caused by a large aspect ratio of the cell dimensions. 15 Because the purpose of this paper is to illustrate the solver structure and model basis, we defer a detailed discussion of model performance to our next paper, in which the model is validated against laboratory tests, large scale experiments and natural hill-slope debris flow events. Here, we discuss only the efficiency of the solver itself, together with tests of the model accuracy in a gravity-driven open channel flow.

Comparison to a drag-force-based Eulerian multiphase model
In comparison to drag-force-based Eulerian multiphase models, the Volume of Fluid approach applied here significantly reduces in calculation time. For an estimate we compared our model with the OpenFOAM standard solver multiphaseEulerFoam. We selected the official tutorial case damBreak4phaseFine but turned the water phase into mercury to gain a 25 three-phase test case, and applied the standard solver settings from the case to our model. On a CentOS 6.3 Linux machine with 31 GiB memory and sixteen Intel Xeon CPU E5-2665 @ 2.40 GHz processors, our model resulted in a 5.5 times faster calculation with a com-23 parable collapse of material columns (Fig. 8). For the sake of completeness our calculation included one iterative viscosity correction step, thus the model efficiency can be estimated to be about ten times higher than a drag-force-based phase coupling approach.

Conclusions
The new debris-flow solver has two main strengths. First, it can model three-dimensional 5 flows and their impact against complexly shaped objects, representing the processes at a high level of detail and reasonable computational costs. Second, its design allows simulating different debris flow material compositions with a reduced calibration process as long as the simulation grid does not change, because the calibration parameters τ 00 and n are largely insensitive to changes in water content, channel roughness or release volume. Due 10 to the used pressure-and shear-dependent rheology, realistic deposit geometries and release dynamics can be achieved, as presented and discussed on the basis of test cases in the next paper. By systematically excluding unknown parameters from the model architecture and by accounting for known flow phenomena in a simplified way, we have developed a debris flow model whose parameters can be roughly estimated from material composition, 15 leaving only two calibration parameters. The concept is promising, however important parts of phase interactions are neglected in favor of lower numerical costs.

Appendix A
The following section describes the detailed implemetation of the PISO iteration procedure as described in Deshpande et al. (2012). By applying the continuum surface force model of 20 Brackbill et al. (1992), the volume integral of Eq. 14 is given as The computational domain is discretized into finite-volume cells. Each cell is considered as the owner of exactly one face that it shares with an adjacent neighbor cell, thus each face has a defined owner cell. A surface normal vector S f with magnitude equal to the surface area of the face is defined on the face pointing outward from the owner cell (Fig. 9). The 5 value at face f of any variable χ calculated in the cell centers as χ P and χ N (Fig. 9) can be derived by interpolation using a mixture of central and upwind schemes: with a weighting factor γ that can account for the flow direction based on the chosen interpolation scheme and flux limiter. In case of a linear interpolation scheme and a flux limiter 10 ψ, γ can be defined as where d is the distance between the cell centers P and N and f N is the distance from the face center to the cell center N . The face flux denoted as φ f serves as a switch to account for the flow direction since it turns negative when 15 the flow is from N to P (Berberović et al., 2009). Several limiters are implemented (OpenFOAM-Foundation, 2016b); we chose the vanLeer scheme and assumed uniform grid spacing to simplify the following explanations with f N/d = 0.5. Variables that are evaluated at the cell faces are subscripted by f . Due to stability problems that arise from the pressure-velocity coupling in collocated meshes 20 (Ferziger and Peric, 2002), the pressure is solved for the cell centers whereas the velocity is interpolated to the cell faces within the PISO loop. 25 With the switch function the velocity U f at face f can be written based on Eq. A2 and A3 as and the corresponding face-perpendicular velocity gradient is given by 5 Deshpande et al. (2012) as At the present time step t n the phase averaged density of the next time step ρ n+1 is known from solving the transport equations. In a first approximation, the corresponding viscosity field µ n+1 can be derived accordingly. A simplified formulation of the momentum 10 equation A1 without pressure, surface tension and gravity terms discretized for cell P could then be formulated as The tilde stands for the velocity at cell P predicted in the current iterative step, for which Eq. 7 yields an explicit expression. The sum over the face density fluxes on the left hand 15 side of equation 7 are known from the mass flux φ ρ derived from equation 9.

26
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | For that purpose, Eq. 5 and 6 are inserted into Eq. 7 using the velocity of the prior iteration step, U m , in all neighbor cells (Deshpande et al., 2012). The explicit expression for the estimated velocity is and by including surface tension and gravity this leads to The detailed composition of H(U m ) and A P formulated with respect to the splitting between neighbor and owner cells can be found in Deshpande et al. (2012); here it is sufficient to keep in mind that H(U m ) contains all off-diagonal contributions of the linear system. The next step is to assemble the approximated face flux where the subscript f indicates that the variable values at the faces are used. The final flux is found by adding the pressure contribution The sum of the flux over the cell faces needs to be zero due to mass conservation for the 15 incompressible flow 27 Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Thus the pressure is defined by the linear equation system for the updated pressure p d m+1 and can be solved with the preconditioned conjugate gradient (PCG) algorithm, to mention one of several options implemented in OpenFOAM. With the updated pressure p d m+1 , 5 the face fluxes φ m+1 f are derived from Eq. 11 and the updated velocity field U m+1 is obtained from the explicit velocity correction which is the end of the PISO loop. After updating the index m to m + 1, the iteration restarts by recalculating H with the updated velocity from equation 8, repeating the loop until a 10 divergence-free velocity field is found.

Code availability
The source code can be downloaded from the supplement application.zip. Please follow the instructions given in the README.pdf file for installation.  . The starting motion (black velocity arrows) with corresponding viscosity distribution of the mixture (left) is a consequence of blending pure shear-rate dependent slurryphase rheology (center) with the pressure-and shear-rate-dependent gravel phase rheology that accounts for Coulomb friction (right). Because the gravel concentration in this example is low, its effect on the overall viscosity is small.