DebrisInterMixing-2 . 3 : a Finite Volume solver for three dimensional debris flow simulations based on a single calibration parameter – Part 1 : Model description

DebrisInterMixing-2.3: a Finite Volume solver for three dimensional debris flow simulations based on a single calibration parameter – Part 1: Model description A. von Boetticher, J. M. Turowski, B. W. McArdell, D. Rickenmann, and J. W. Kirchner Department of Environmental Systems Science, ETH Zentrum, CHN H41, 8092 Zürich, Switzerland Helmholtz-Centre Potsdam GFZ German Research Center for Geosciences, Telegrafenberg, 14473 Potsdam, Germany Swiss Federal Research Institute WSL, Zürcherstrasse 111, 8903 Birmensdorf, Switzerland


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 fluid 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 intense settlement in mountainous regions 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.Modeling debris flows is a central part of debris-flow research, because measuring the detailed processes in debris-flow experiments or in the field is challenging.It is still uncertain 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 that are difficult to track with the present measurement techniques at reasonable cost.As a consequence, the rheological behavior of debris flow material is incompletely understood.
Typically, existing numerical modeling approaches cannot predict run-out distances or impact pressures of debris flows in known terrain without prior parameter calibration, based on simulating previous well-documented events that happened at the same site.
Clearly, this represents a challenge in practical applications, because often reliable calibration data are unavailable.In addition, the interactions between the granular and viscous phases, and the dynamic change in granular and viscous concentrations during the flow process, limit simple models to the narrow range of simulations that they have been calibrated for, where the fitted parameters account for these interactions.Complex models such as depth-averaged fluid simulations coupled to three dimensional particle methods are associated not only with high computational costs but also with a large number of model parameters, making model calibration the key issue for application to specific cases.This limits the possibilities of using debris flow models as a valid standard application in practice, because the user's ability to estimate values of poorly constrained parameters influences the results.
Here, we present an improved three-phase modeling approach as an alternative.We provide a coarse 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.The parameters of the two resulting rheology models for the two mixing phases are linked to material properties such that Introduction

Conclusions References
Tables Figures

Back Close
Full the model setup can be based on material samples collected from the field, yielding a model that has only one free parameter for calibration.
2 Modeling approach 2.1 A two phase model with pressure-and shear-rate-dependent rheology Based on a Finite Volume solver for a mixed three-phase incompressible Navier-Stokes equation, we apply two rheology models for two phases that can mix: one for the suspension of fine material and water, and one for gravel.We allow interactions of both rheologies, while keeping a third phase unmixed, accounting for the air forming the free surface.In this way, the coupling between driving forces, topography and flow-dependent internal friction can be addressed for each phase separately, accounting for the fundamental differences in flow mechanisms of granular and viscous flow (Fig. 1).Numerical costs are kept reasonable by using the Volume of Fluid method (Hirt and Nichols, 1981) such that only one Navier-Stokes equation system is solved for all three phases.We calculate the viscosity and density of each grid cell as a concentration-weighted average between the viscosities of the phases that are present in the cell.Phase interaction is reduced to this averaging of density and viscosity with the aim to avoid the standard approach of phase coupling by the use of drag force models.The drag between non-spherical grains and non-Newtonian fluids is still not well understood and would lead in our case to a model with a high number of unknown parameters.However, although non-linear concentration-weighted averaging of viscosities for estimating the viscosity of a mixture is common in the oil industry (Gao and Li, 2012), we apply linear averaging for simplicity, although non-linear averaging of viscosity between phases may be introduced in the future.The central assumption for concentration-weighted averaging is that the local velocity of the gravel is about the same as the velocity of the surrounding fluid.This assumption would not be valid for debris flows with little interstitial fluid content, or with interstitial fluid of small viscosity (i.e., Introduction

Conclusions References
Tables Figures

Back Close
Full a slurry with low concentrations of fine material).The assumption of equal velocities of both phases in one cell leads to a constant distribution of phase concentrations over the entire flow process.Nevertheless, this assumption avoids the need to model the drag forces between gravel and interstitial fluid, while still accounting for the pressuredependent flow behavior of the gravel in combination with the shear-dependent rheology of the slurry.
Because pressure and shear drive the energy 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 (2001b) and Pitman and Le (2005) treated the granular phase as a shear-rate independent Mohr-Coulomb plastic material, dry granular material was successfully modeled as a viscoplastic fluid by Ancey (2007), Forterre and Pouliquen (2008), 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 pressure-dependent Coulomb-viscoplastic behavior of the granular phase, as applied by Domnik et al. (2013).But instead of solving Navier-Stokes equations for each phase coupled by drag models, we apply the numerically more efficient method of Iverson and Denlinger (2001b) and treat the debris flow material as one mixture with phase-averaged properties described by a single Navier-Stokes equation.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.

Rheology model for the fine sediment suspension
The introduction of two phases is necessary to distinguish between the pressuredependent flow behavior of gravel and the shear-thinning viscosity of the suspension of finer particles with water.The rheology of mixtures of water with clay and sand can be described by the Herschel-Bulkley rheology law (Coussot et al., 1998), which defines Introduction

Conclusions References
Tables Figures

Back Close
Full 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 factor for the viscosity of the sheared material, γ is the shear rate and n defines the shearthinning (n < 1) or shear-thickening (n > 1) behavior.With n = 1 the model simplifies to the Bingham rheology model that was widely used to describe debris-flow behavior in the past.It may be reasonable to 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 was rarely applied in debris-flow modeling due to the large number of rheology parameters.We avoid this problem by assuming the rheology parameters to be defined by measurable material properties as described below.

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 one.If the coarser grain fraction is confined to the gravel phase, the Herschel-Bulkley parameters for the finer material can be linked to material properties as measured using simple standard geotechnical tests.According to Coussot et al. (1998), the exponent n can be assumed constant as 1/3, and k can be roughly estimated as 0.3 • τ y for mixtures with maximum grain-sizes < 0.4 mm (Kaitna et al., 2007).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: (2) Introduction

Conclusions References
Tables Figures

Back Close
Full where C is the volumetric solid concentration of the mixture, P 1 = 0.7P 0 when P 0 > 0.27 and P 1 = P 0 if P 0 is smaller, and where the subscript of C refers to the volumetric concentration 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 where τ 00 is the remaining free parameter which we use to account for the grid size dependency of the shear rate.We recommend a value of τ 00 = 30 Pa as a starting point for calibration.Yu et al. (2013) compared this method of estimating the yield stress τ y to experimental results they gained 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 determined as: where ρ is the density of the mixture, g the acceleration due to gravity, h the maximum accumulation thickness of the deposit, and φ the slope inclination.In addition, they compared the calculated yield stress of Eq. ( 2 of sand-clay mixtures with water can be estimated using Eq. ( 2) based on the volumetric concentration of the debris in the water-solids mixture and based on the percentages of different clays in the fraction of fine material.Adjustments to the numbers for calculating P 0 may be necessary to account for the activity of other clays.
The remaining uncertainties concern our assumptions are that n is constant at a value of 1/3, and that k = 0.3 • τ 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, and Remaitre et al. (2005) found n to vary from 0.27 to 0.36.Based on the laboratory scale experiments that are presented in von Boetticher et al. ( 2015) we have chosen n = 0.34 to obtain the best fit for the simulation of large-scale debris-flow experiments.

Representation of gravel by a Coulomb-viscoplastic rheology
For a satisfactory prediction of run-out and impact, an adequate simulation of the deposition process is necessary.During acceleration and high-speed flow, the shearthinning behavior of both the viscid and the granular phase dominate the viscosity.
However, pressure-dependent friction becomes important as soon as the material experiences high pressures, accompanied by reduction in shear due to decelerations caused by channel slope reduction.Flows of granular material have been successfully 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: where D is the rate-of-deformation tensor, pI is the normalized pressure times the identity matrix and ν s is the corresponding kinematic viscosity, which was modeled by Domnik and Pudasaini (2012) as:

Conclusions References
Tables Figures

Back Close
Full where ν min is a minimal kinematic viscosity, τ 0s is a density-normalized yield stress, and ||D|| is the norm of the strain-rate tensor: In Eq. ( 7), m y is a model parameter which we will keep constant as reasoned in the following section.Domnik et al. (2013) suggested replacing the yield stress by a pressuredependent Coulomb friction, p • sin(δ) where δ is the internal friction angle and p is the density-normalized pressure: Here, this Coulomb-viscoplastic rheology model is used to describe the gravel phase, by calculating the pressure-and shear-dependent viscosity in every cell.

Gravel phase properties
The For small friction angles, the modeled viscosity of the gravel phase decreases rapidly with increasing shear.Larger friction angles increase the viscosity and extend the pressure dependency to larger shear rates (Fig. 3).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.
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)) (where ν is the local kinematic viscosity and p the local densitynormalized pressure).Following Renardy (1986), we locally limit the viscosity to fulfill this stability criterion.

Solver description
The assumption of negligible velocity differences between the gravel particles and the slurry within a finite-volume cell allows the solution of an averaged Navier-Stokes equation for the three phases air, gravel and fluid.Each phase is treated as a continuum.The phases for the slurry and for the gravel inter-penetrate each other, while the air phase is kept separate.The conservation of mass and momentum is averaged with respect to the phase fraction α of each phase.The phase fraction is the probability that a phase is present at a certain point in space and time (Hill, 1998).With the phase fractions as a representation of the phase volume distribution, a transport equation for each phase can be solved to obtain the change of phase distribution over the last time step.With the updated phase distribution, the pressure and velocity fields are calculated by solving the phase-averaged Navier-Stokes equations.The Volume-offluid (VOF) method (Hirt and Nichols, 1981) is used to reconstruct the free surface with Introduction

Conclusions References
Tables Figures

Back Close
Full convection schemes from the volume fraction distribution (Rusche, 2002).The method allows a surface reconstruction with high accuracy if the grid resolution is sufficient (Hänsch et al., 2013;Hoang et al., 2012).Figure 4 illustrates how the phase volume distributions α 1 (air), α 2 (slurry) and α 3 (gravel) are used to derive cell-averaged properties of the continuum.
Because most debris-flow models are depth-averaged and use shallow-water approximated equations, it could be questioned why a three-dimensional approach is necessary.Brodani-Minussi and deFreitas Maciel (2012) compared dam-break experiments of a Herschel-Bulkley fluid and its numerical simulations using the VOF approach with published shallow-water equation based models.Especially for the first instant after the material release, the application of shallow-water equations seems to introduce errors that are propagated throughout the process leading to erroneous runout estimates.In addition to the three-dimensional approach, we introduced an iterative step to determine the shear-dependent viscosity without delay for the model to be able to deal with the challenges of a dam-break release.
We describe the solver below, beginning with a brief introduction to the PISO (Pressure-Implicit with Splitting of Operators Issa, 1986) algorithm for solving the incompressible Navier-Stokes equations for phase-averaged mass and momentum conservation.We then address the numerical solution for each phase flow using advectiondiffusion equations with focus on the mixing, and finally we describe some aspects of the dependencies between grid resolution, rheology models and solver stability are described.The incompressible Navier-Stokes momentum equation takes the form:

Calculation of the velocity field
where u i and u j (i , j = 1, 2, 3) are the velocity components in the Cartesian directions 1, 2, 3 at a place with coordinates x i and x j .p stands for the local density-normalized pressure and ν denotes the kinematic viscosity of the fluid.Integrating Eq. ( 10) over the volume V of a finite cell of a grid-discretization of the simulated space leads to a cellsurface based conservation of momentum for the volume V by applying the Gauss Theorem.Integrals over the volume are thereby replaced by integrals over the cell surface: The Finite Volume Method replaces the integral over a cell surface by a discrete value, obtaining Eq. ( 12) from Eq. ( 11).The index k denotes values at cell face k, and field parameters without indices denote the corresponding value in the middle of the cell.
Together with the conservation of mass, this leads to an equation system defining the velocities depending on pressure for an incompressible fluid.To reduce the system of equations to the number of unknowns it is necessary to calculate the values at the cell surfaces from interpolation between the values at the cell centers of neighboring cells using interpolation schemes, Javurek (2006) gives a brief description and summary of implemented schemes.All simulations were carried out using first-order Euler schemes for the time derivative terms, as has been recommended for liquid column breakout simulations (Hänsch et al., 2013).Standard Gaussian finite volume integration with linear 6360 Introduction

Conclusions References
Tables Figures

Back Close
Full interpolation, which is of second order with unbounded numerical behavior, was chosen for the gradient derivative, convection and Laplacian terms.A limited surface-normal gradient scheme was applied that blends corrected and uncorrected treatment of cell orthogonality.An unconventional compressive scheme could be used to sharpen the interface, but here a conventional upwind scheme was used, because more research is necessary on the performance of the compressive scheme (Raees et al., 2011).
The PISO algorithm uses an optional implicit predictor for the velocity field followed by explicit corrector steps.The predictor uses the pressure field of the old time step to estimate a velocity field for the current time step which is in general not divergence-free.With the idea that a correct velocity field should be divergence-free due to continuity, a Poisson equation for the first corrected pressure field at the current time step is formed by taking the divergence of the equation that defines velocity as a function of pressure.With the corrected pressure, a corrected velocity field can be calculated and the corrector step can be repeated until divergence-free velocities are found.For a detailed description with the corresponding matrix equations, implementation and further literature see OpenFOAMWiki (2015).

Solving the multi-phase flows
By knowing the velocity field from the previous time step, the current time step starts with calculating the current phase concentrations, accounting for the changes by advection and dispersion.
where D diff is the diffusion constant.Because diffusion is neglected in our model, Eq. ( 13 known from the current velocity field.To keep the air phase unmixed, it is necessary to determine the flux φ r through the surfaces between air and the debris flow mixture, and to subtract it from the calculated phase fluxes φ 1..3 to achieve an impermeable surface.Inherited from the original interMixingFoam solver (OpenFOAM-Foundation, 2014), limiters are applied during this step to bound the fluxes to keep phase concentrations between 0 and 1.With known fluxes φ 1..3 for each phase, the scalar transport equation for each phase takes the form where i = 1, 2, 3 denote the phases air, slurry and gravel.
It is necessary to limit the fluxes or the phase concentrations such that the sum of volumetric concentrations in each cell adds up to one; otherwise the concept of incompressible flow is violated.However, this limiting constraint can lead to changes of total phase volumes conflicting with the conservation of mass.The standard implementation of the interMixingFOAM solver encounters such difficulties.However, we achieved a good solution in our modified code by first solving Eq. ( 14) for the slurry phase, then limiting the updated slurry concentrations to values greater than or equal to zero, then solving Eq. ( 14) for the gravel concentrations and limiting the updated gravel concentrations to a range between 0 and (1.0 -slurry concentration).
Because the influence of the air phase on the debris flow simulation is of small importance, a simple solution is obtained when deriving the final air phase concentration as In this way, for each phase of the debris flow material, negative phase concentrations or values larger than one can be avoided, while still ensuring that the sum of all three phase concentrations adds up to one.The error is concentrated mainly in the air phase and the gravel phase and the chosen bounding results in a stable solver with sufficient conservation of mass for the debris flow material.Introduction

Conclusions References
Tables Figures

Back Close
Full After solving the transport equations, the complete mass flux φ ρ from the updated volumetric phase concentrations is constructed: where ρ 1..3 denotes the densities of the corresponding phases and φ 1..3 the corresponding fluxes.

Effect of grid resolution and time step on rheology
Since the shear rate influences both viscosity models, a strong influence of grid resolution on viscosity results, 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 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 viscous rheology model, and can be used to adjust the simulation to the grid resolution.
As long as the gravel phase and grid resolution do not change, it should be possible to model different water and clay contents based on one calibration test.However, as the composition changes, both τ y and τ 00 must change commensurately, since the change in yield stress affects the shear rate.Our procedure for 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 τ 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 6363 Introduction

Conclusions References
Tables Figures

Back Close
Full 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 granular 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 no front instability is caused by a large aspect ratio of the cell dimensions.
Another major problem in many models that simulate the release or impact of material with velocity-dependent rheology is that the viscosity is kept constant over each time step although it actually changes with the changing flow.In our model, during release of immobile material that accelerates, the viscosity is overestimated over each time step.As a consequence, the velocity at the end of the time step is underestimated, which again amplifies the overestimation of viscosity in the next time step.Conversely, at an impact, the sudden deceleration causes an underestimation of viscosity over the time step length, leading to an overestimated velocity that again amplifies the underestimation of the viscosity in the next time step.In both situations, the error sums up from time step to time step.As a result, flow velocities change with changing time step size.Avalanche codes such as RAMMS (Christen et al., 2007) deal with this problem by calibrating the model to data from previous 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 change 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 accurately simulate impacts.Here, our model constitutes a significant improvement, since in the three-dimensional solver presented here, 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 Introduction

Conclusions References
Tables Figures

Back Close
Full step can be 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.

Discussion
Field observations and experiments indicate that the debris flow rheology varies from nearly rigid to highly fluid due to local and temporal variations of pore fluid pressure, shear rate and particle dynamics.In the past, debris flow models have often treated the flowing mixture as a single homogeneous phase with either viscous or granular flow characteristics.The traditional approach to the fluid dynamics of debris flows accounted for the vertical momentum exchange in the flow process in a simple manner by assuming a velocity profile over the flow height (concept of depth-averaged shallowwater equations).A good presentation of depth-averaged single phase models, sorted by complexity, is given in Luca et al. (2009).Granular debris flow models addressed the energy dissipation through dispersive shear stress, kinetic stress and collision stress (Savage, 1984).Dispersive shear stress is caused by the friction between particles in contact as they move past one another during the macroscopic shearing motion.Kinetic stress arises when particles at one level in the velocity profile move up or down to another level, and collision stress should account for the sum of energy dissipation due to particle collision (Takahashi and Tsujimoto, 2000).The model of Savage and Hutter (1989) is an example of such a granular perspective linking the energy dissipation in a depth-averaged approach to the basal friction.However, the obvious dependency of the bed friction angle on local topography and velocity led to models with friction angles varying during simulation.One approach was to change the basal friction angle as a function of the Savage number (Iverson et al., 2010) leading to a shear rate and Introduction

Conclusions References
Tables Figures

Back Close
Full grain size dependent energy dissipation providing a first attempt to link the rheology of the flow to known material properties.
In contrast to the approaches outlined above, which addressed debris flows from the granular avalanche perspective, viscoplastic (Coulomb-viscous) theories postulated debris flows as a homogeneous viscoplastic continuum (Iverson and Denlinger, 1987).In viscoplastic approaches, the mechanical behavior of debris flow material is seen as dominated by the influence of a muddy matrix that fills the space between coarser grains.Initially related to Bingham fluids, this matrix behaves like a solid if shear stresses do not exceed a yield stress, and like Newtonian fluids with constant viscosity, if the yield stress is exceeded.If the yield stress is modeled as dependent on the normal stress acting on planes of shearing, one obtains the Coulomb-viscous model (Johnson, 1965(Johnson, , 1970(Johnson, , 1984)).This model can reproduce the ability of debris flows to move steadily over different slopes.The concept of a yield stress could also explain the observed concentration of deformation in thin bands of sheared layers close to the flow boundaries in channelized debris flows.To account for shear thinning, the Bingham rheology was replaced by a Herschel-Bulkley rheology that can reduce the viscosity at higher shear rates (shear thinning) or increase it (shear thickening), depending on parameter settings.While sand and clay mixtures with water show a shear thinning effect, more granular mixtures exhibit shear thickening (Lashgari et al., 2014).Iverson (1997) generalized the granular flow model of Savage and Hutter (1989) to account for the presence of pore fluid.With a mixing theory framework, this model was further developed (Iverson and Denlinger, 2001a) into the Coulomb mixture model, which could simulate a wide spectrum of grain-fluid flows from initiation to deposition, with no redefinition of parameters.It opened the transition from single-phase models to flows composed of solid-fluid mixtures.Iverson (2003) pointed out that models that treat the debris flow material as a single phase with one rheology are not capable of representing the real behavior, and that a two-phase approach is necessary where one phase accounts for the viscid forces while the second phase models the granular forces between the grains.Pudasaini fying the single-phase models of Savage and Hutter (1989), the debris-mixture model of Iverson and Denlinger (2001a) and the two-fluid debris flow model of Pitman and Le (2005).Our model can be considered as based on the concept of the Coulomb mixture model (Iverson and Denlinger, 2001a) but with a state-of-the-art Herschel-Bulkley representation of the fluid and a pressure-dependent Coulomb-viscoplastic representation of the gravel (Domnik et al., 2013) in a three-dimensional approach without depthaveraging.
In comparison to drag-force-based Eulerian multiphase models, the Volume of Fluid approach applied here provides significant reduction in calculation time.For an estimate we compared our model with the OpenFOAM standard solver multiphaseEuler-Foam.We selected the official tutorial case damBreak4phaseFine but turned the water phase into mercury to gain a 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 comparable collapse of water columns (Fig. 5).
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.
The model was also applied to an open clear water channel experiment with about 50.6 L s −1 discharge in a 40 m long and 1.1 m wide rectangular smooth channel with 0.026 % inclination (Fischer, 1966).The slurry phase was initialized as water together with a zero gravel phase concentration.as a comparison between a measured and simulated vertical velocity profile suggests (Fig. 6).

Conclusions
The new debris-flow solver has two main strengths.First, it can model threedimensional flows and their impact against complexly shaped objects, representing the processes at a high level of detail.Second, its design allows simulating different debris flow material compositions without recalibrating the one free parameter, as long as the simulation grid does not change.Due to the solver's 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 accompanying paper.By systematically excluding unknown parameters from the model architecture and by accounting for known flow phenomena, we have developed a debris flow model whose parameters can be estimated directly from site geometry and material composition rather than from extensive calibration.The concept is promising, but due to high numerical costs and long calculation times the model is still limited to small simulations of several hundred square meters in surface area.Introduction

Conclusions References
Tables Figures
Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | ) with experimental yield stresses reported by a number of authors includingCoussot et al. (1998) andAncey and Jorrot (2001).Ancey and Jorrot (2001) used 2 and 3 mm glass beads in a kaolinite dispersion as well as fine sand-kaolinite-water mixtures.Up to yield stresses of about 200 Pa the yield stresses estimated by Eq. (2) fit the observed ones well.Thus, the yield stresses Introduction Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Coulomb-viscoplastic rheology law developed byDomnik et al. (2013) includes two parameters: the friction angle δ, and the parameter m y influencing the 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 density-normalized pressure, m y needs to be equal to one.However, the development of ν s under large pressure or strong shear is the same for both m y = 1 and m y = 0.2, but parts of the nearly immobile material that face little pressure (in general, immobile material close to the surface) show a significant reduction in viscosity when m y = 0.2 (Fig.2).As a consequence, m y minimally affects debris flow release and flow at large scales, but material with a low 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.We choose m y to be constant and equal to 0.2 for all simulations.Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | The solver is based on an adaption and extension of the interMixingFoam solver, which is one of the standard solvers of the open source Finite Volume Code OpenFOAM (OpenFOAM-Foundation, 2014).The Finite Volume method used here is based on a discretization of the incompressible Navier-Stokes equation to describe the fluid dynamicsDiscussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | ) reduces to the advection equation, which can be solved based on the advective phase fluxes φ 1..3 for each phase.The phase fluxes 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 Introduction Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | Discussion Paper | (2012) developed a general two-phase flow model, uni-Discussion Paper | Discussion Paper | Discussion Paper | A Hybrid URANS-LES model was applied to account for the turbulent flow.Instead of an inlet discharge the model applied periodic inlet and outlet boundary conditions and the flow was driven by gravity.The debrisInter-MixingFoam solver predicted the discharge of the turbulent channel flow with an underestimation of 15 % and underestimated the corresponding surface elevation by 2.5 %.However, the deviations in predicted and measured average flow velocities are probably related to shortcomings of the URANS turbulence model at the bottom boundary, Discussion Paper | Discussion Paper | Discussion Paper |

Figure 1 .Figure 2 .Figure 3 .Figure 4 .Figure 5 .
Figure 1.Viscosity distribution (indicated by color scale) along a 28 cm long section through the modeled 0.01 m 3 release block 0.2 s after release, corresponding to the experimental setup of Hürlimann et al. (2015).The starting motion (black velocity arrows) with corresponding viscosity distribution of the mixture (left) is a consequence of blending pure shear-rate dependent mud phase rheology (center) with the pressure-and shear-rate-dependent gravel phase rheology (right).Because the gravel concentration in this example is low, its effect on the overall viscosity is small.