Interactive comment on “ Comparing Thixotropic and Herschel-Bulkley Models for Avalanches and Subaqueous Debris Flows

Debris flows such as avalanches and landslides are heterogeneous mixtures of solids and liquids but are often simulated as homogeneous non-Newtonian fluids using a Herschel-Bulkley model. By representing the heterogeneous debris as a homogeneous non-Newtonian fluid, it is possible to use standard numerical approaches for the Navier-Stokes equations where viscosity is allowed to vary in time and space (e.g. eddy-viscosity turbulence models). Common non-Newtonian models are time-independent so that the relationship between the time-space-varying effective viscosity and flow stress is unchanging. However, the complex behaviors of debris flows at flow initiation (jamming) and cessation (restructuralization) imply that the viscosity-stress relationships should have time-dependent behaviors, which is a feature of thixotropic non-Newtonian fluids. In this paper, both Herschel-Bulkley and thixotropic non-Newtonian fluid models are evaluated for simulating avalanches along a slope and subaqueous debris flows. A numerical solver using a multi-material level set method is applied to track multiple interfaces simultaneously. The numerical results are validated with analytical solutions and available experimental data using parameters selected based on the experimental setup and without post-hoc calibration. The thixotropic (time-dependent) fluid model shows reasonable agreement with all the experimental data. For most of the experimental conditions, the Herschel-Bulkley (time-independent) model results were similar to the thixotropic model, a critical exception being conditions with a high yield stress. Where the flow initiation is strongly dominated by the structural jamming and the initial yield behavior the time-independent model performed poorly.


Introduction
Avalanches, landslides, mudflows, and volcanic lava form hazardous gravity-driven debris flows that are typically multiphase heterogeneous mixtures of solids and liquids (Pudasaini, 2012;Davies, 1986).Debris flows in the real world show timedependent (thixotropic) characteristics (Aziz et al., 2010;Bagdassarov and Pinkerton, 2004;Crosta and Dal Negro, 2003;Perret et al., 1996).However, time-independent rheological models have been widely used to simulate debris flows (Bovet et al., 2010;Manga and Bonini, 2012;Pirulli, 2010;Tsai et al., 2011).Following Ancey (2007), approaches to simulating debris flows can be categorized in three groups: (i) applying soil mechanics concept of Coulomb behavior, (ii) merging soil and fluid mechanics models, and (iii) representing the heterogeneous debris as a solid/liquid mixture with behaviors similar to a non-Newtonian fluid.This research uses the third approach.The interaction between solid particles and surrounding fluid is 1 Nat.Hazards Earth Syst.Sci.Discuss., https://doi.org/10.5194/nhess-2017-258Manuscript under review for journal Nat.Hazards Earth Syst.Sci. Discussion started: 17 July 2017 c Author(s) 2017.CC BY 4.0 License.modeled as creating the microstructural behavior of a yield stress fluid that is approximated using non-Newtonian fluid models.
The main advantage of this approach is that rheological models for debris flows are easily included in the viscosity term of the Navier-Stokes equations.It follows that a wide range of existing hydrodynamic codes can be readily adapted to non-Newtonian behavior and used as models of debris flows.
From a macroscale perspective, debris flows have similar behaviors to "yield stress fluids" that have been studied as a class of non-Newtonian fluids (Moller et al., 2009;Scotto di Santolo et al., 2010).A yield stress fluid is a solid (or flows with extremely high viscosity) below the critical stress value (the yield stress).This behavior is similar to what might be expected under low stress conditions from a debris mixture of liquid and solids that is initially at rest.At the microscale under low stress (near rest) conditions the fluid flow around the solids in a debris mixture is inhibited by viscous boundary layers and inertia of the solids, which provides effects similar to a higher-viscosity fluid at the macroscale (i.e.low deformation under stress).
Once the solids in the debris have accelerated the effects of particle lift, drag, and rotation induced by the surrounding turbulent fluid flow, as well as solid-solid impacts and particle disintegration, will provide to behaviors similar to a lower-viscosity fluid that deforms more easily under stress.The destruction of the initial microstructure of the debris will change the effective macroscale viscosity and response to stress.Thus, we can think of the behavior of a debris flow as controlled, at least partly, by the evolution of the microstructure.
Non-Newtonian yield stress fluids can be classified as either time-independent or time-dependent.The former, such as Bingham plastics, have a constant and repeatable relationship between effective viscosity, shear stress, and yield stress.In contrast the viscosity relationship in a time-dependent, or "thixotropic" fluid evolves as the flow accelerates from rest (Moller et al., 2009).The terminology can be confusing because both approaches allow the viscosity to vary over time, but only the thixotropic fluid allows the relationship between viscosity and stress to change with time.The approach proposed by Herschel and Bulkley (1926), commonly known as the Herschel-Bulkley model, is the standard approach for representing the general case of time-independent Bingham plastics as well as shear-thinning and shear-thickening fluids.In this approach, the plastic viscosity, η, is conditional on the yield stress, τ 0 , as where K is the consistency parameter, n is the Herschel-Bulkley fluid index, and γ is the scalar value of the rate of strain.
The Herschel-Bulkley fluid index n controls the overall modeled behavior, where 0 < n < 1 is shear thinning, n > 1 is shear thickening and n = 1 corresponds to the Bingham plastic model (Bingham, 1916).
A recognized problem with numerical simulation of a fluid using a Herschel-Bulkley model is the viscosity is effectively infinite below the yield stress, i.e. the condition γ = 0 in Eq. ( 1) is identical to η = ∞ for modeling a fluid continuum that becomes solid below the yield stress.An infinite (or even very large) viscosity creates an ill-conditioned matrix in a discrete solution of the partial differential equations for fluid flow.Furthermore, the instantaneous transition from infinite to finite viscosity as the yield stress is crossed provides a sharp change that can lead to unstable numerical oscillations.Dent and Lang (1983) attempted to resolve this issue with a bi-viscous Bingham fluid model for computing motion of snow avalanches.Their approach was shown to be reasonable using comparisons with experimental data, but was later determined to be invalid for conditions where the shear stresses are much lower than the yield stress (Beverly and Tanner, 1992).A more successful approach was that of Papanastasiou (1987), who proposed modifying the Herschel-Bulkley model with an exponential parameter, m.
The Papanastasiou model (presented in detail in Section 3, below), with appropriate values for m, shows good approximations at low shear rates for Bingham plastics (Beverly and Tanner, 1992).
Although a flow simulated with the Papanastasiou model will have changes in the viscosity with time (as the shear changes with time), the model is still deemed "time-independent" as the relationship between viscosity and shear is fixed by the selection of K, n, m, and τ 0 .Arguably, there exists a wide range of debris flows over which the Papanastasiou approach should be adequate, as the time-dependent characteristics of debris flows are, at least theoretically, principally confined to the initiation add cessation of the flow, i.e. when the microstructure of the debris is evolving and changing the relationship between shear and viscosity.It follows that steady-state conditions for debris flows should be reasonably represented with time-independent models.Indeed, O'Brien and Julien (1988) concluded, by their experiments, that mud flows whose volumetric sand concentration were less than 20% showed the behavior of a silt-clay mixture, which can be described reasonably well by the Bingham plastic model at low shear rates and a time-independent Herschel-Bulkley model at high shear rates.Liu and Mei (1989) reported good agreement for theory and experiment with a Bingham plastic model and a homogeneous mud flow that provides a steady front propagation speed (necessarily long after the initiation phase).The Herschel-Bulkley model has also been used to simulate debris flow along a slope, but reported results have discrepancies with experimental data, especially in the early stages (Ancey and Cochard, 2009;Balmforth et al., 2007).Bovet et al. (2010) applied the time-independent Papanastasiou model to simulate snow avalanches with some success, but again their results showed more significant discrepancies with experiments during flow initiation.De Blasio et al. ( 2004) simulated both subaerial and subaqueous debris flows with a Bingham fluid model.Their results for the subaerial debris flows were in a reasonable agreement with laboratory data, but their subaqueous simulations showed a significant discrepancy with measurements.A clear challenge in validating models of debris flows beyond steady conditions is that the most commonly available experimental data is focused on the steady or quasi-steady conditions after the debris structure has (relatively) homogenized.
Thixotropic (time-dependent) behavior, which is not represented in the Herschel-Bulkley model, provides an interesting avenue for representing the expected macro-scale behavior of a debris flow near initiation.At rest, debris solids provide structural resistance to flow (for denser solids), and a greater inertial resistance to motion than the fluid.Thus, it is reasonable to expect initial behavior similar to a Bingham plastic; i.e. initially-infinite viscosity with a high yield stress.However, the onset of motion for the debris flow begins the destruction of the microstructure, homogenization of the debris, and a change in the relationship between stress and viscosity, which might be thought of as shear-thinning behavior.A key difference between a Herschel-Bulkley model and the real world is that the former requires a return to structure whenever the internal stress drops below the yield stress, however, in a debris flow we expect the destruction of microstructure to significantly reduce the stress at which re-structuralization occurs.For a real debris flow we expect different viscosity-stress behaviors during initiation, steady-state, and slowing phases (consistent with evolving microstructure), but a time-independent Herschel-Bulkley model is effectively an assumption that the fluid is unaffected by the phase of the flow (i.e. the scale of the microstructure is constant).
For a thixotropic fluid the time-dependency can occur as part of spatial gradients that evolve over time; e.g.high shear stress is localized in a small region by heterogeneity of particles, and in this region the fluid begins to yield (Pignon et al., 1996).Thus, in a thixotropic fluid there is spatial-temporal destruction of microstructure that leads to changes in the effective viscosity that cannot be represented in the standard time-independent models.Coussot et al. (2002a) proposed an empirical viscosity model for thixotropic fluids (presented in detail in Section 3.1, below), which captures these fundamental behaviors.
Prior research on thixotropic avalanches and subaqueous debris flows has mainly focused on laboratory experiments (Chanson et al., 2006;Haza et al., 2013;Mohrig et al., 1999;Sawyer et al., 2012), although a few studies have numerically investigated the characteristics of thixotropic flow on a simple inclined plane (Hewitt and Balmforth, 2013;Huynh et al., 2005).In general, numerical simulation results have not been well validated by the experimental data, arguably due to limitations in both non-Newtonian viscosity models and the sparsity of available laboratory data.
In this paper we evaluate a time-independent Papanastasiou model and a time-dependent Coussot model for simulations of a thixotropic avalanche and subaqueous debris flow, with comparisons to available experimental measurements.The governing equations are presented in Section 2, and the non-Newtonian Papanastasiou and Coussot viscosity models in Section 3.1.A key confounding issue for model/experiment comparisons is the estimation of parameters for a non-Newtonian fluid model (in particular the initial degree of jamming), which we discuss in Section 3.2.The numerical solver, using a multi-material level set method, is presented in Section 4. The solver is validated in Section 5 with the analytical solutions for the Poiseuille flow of a Bingham fluid.In Section 6 the solver is used to model a laboratory flow that is a reasonable proxy of a thixotropic avalanche.
In Section 7 we present the numerical simulations of subaqueous debris flows with three interfaces: debris-water, debris-air, and water-air, and compare our results to prior experimental data.We draw conclusions from our work in Section 8.
2 Governing Equations: The governing equations in conservation form for unsteady and incompressible fluid flow can be written as (Ferziger and Perić, 2002) where u is the velocity vector, ρ is the density, p is the pressure, f includes additional forces such as gravitational force, surface tension force, and Coriolis force, u ⊗ u is the dyadic product of the velocity vector u, and T is the viscous stress tensor: where η denotes the plastic viscosity and D is the rate of strain (deformation) tensor: makes the approach compatible with a wide range of numerical solvers.
Governing equations ( 2) and ( 3) can be integrated over a control volume and, by applying the Gauss divergence theorem, we obtain the basis for the common finite-volume numerical discretization (Ferziger and Perić, 2002).For simplicity in the present work, we limit ourselves to a two-dimensional flow field for a downslope flow and the orthogonal (near-vertical) axis, which effectively assumes uniform flow in the cross-stream axis.The external force term f represents the gravitational force only, neglecting surface tension forces and Coriolis.The advection term is discretized with the fifth-order WENO (Weighted Essentially Non-Oscillatory) scheme (Shi et al., 2002) or the second-order TVD (Total Variation Diminishing) Superbee scheme (Darwish and Moukalled, 2003) in separate numerical tests.The diffusion term on the right hand side of Eq. ( 3) is discretized with the second-order central differencing scheme.The time derivative term for the momentum equations is integrated by the second-order Crank-Nicolson implicit scheme.The deferred correction scheme (Ferziger and Perić, 2002) is applied and ghost nodes are evaluated by the Richardson extrapolation method for high accuracy at the boundaries.The pressure gradient term is calculated explicitly and then corrected by the first-order incremental projection method (Guermond et al., 2006).To evaluate the values at the cell surfaces, the Green-Gauss method is used and the momentum interpolation scheme (Murthy and Mathur, 1997) is applied.The code is parallelized with MPI (Message Passing Interface), and PETSc (Portable, Extensible Toolkit for Scientific Computation) (Balay et al., 2016) is used for standard solver functions (e.g. the stabilized version of Biconjugate Gradient Squared method with pre-conditioning by the block Jacobi method).The developed code has been verified by the method of manufactured solutions (Jeon, 2015).

Time-independent and time-dependent model
The Herschel-Bulkley model, Eq. ( 1), was made more practical for modeling a fluid flow continuum by Papanastasiou (1987), whose approach can be represented as : second invariant of the rate of strain as (Mei, 2007) and D ij denotes the (i, j) component of the strain tensor D in Eq. ( 5).As with the Herschel-Bulkley model on which it is based, the Papanastasiou model is time-independent.
In contrast, the time-dependent Coussot model (Coussot et al., 2002a) introduces dependency on a time-varying microstructure parameter (λ) in the general form: where η 0 is the asymptotic viscosity at high shear rate, ω is a material-specific parameter, and n is the Herschel-Bulkley fluid index.The microstructural parameter of the fluid, λ, is evaluated using a temporal differential equation: where T 0 is the characteristic time of the microstructure, α is a material-specific parameter, and γ is the rate of strain (as in the Herschel-Bulkley and Papanastasiou models, above).Here α represents the strength of the shear effect associated with inhomogeneous microstructure (Liu and Zhu, 2011).That is, larger of values of α require greater microstructure homogenization (smaller λ) to drive the system to steady-state conditions (dλ/dt → 0).

Estimation of parameters for time-dependent Coussot model
The time-dependent Coussot model requires parameters for the asymptotic viscosity (η 0 ), Herschel-Bulkley fluid index (n), characteristic time (T 0 ), and two material-specific parameters (ω and α) that control the response (destruction) of the microstructure.Additionally, an initial condition for λ 0 is required to solve ODE (9).The parameters η 0 and n are easily obtained from the time-independent Herschel-Bulkley model, which are typically available in experimental studies.However, the other parameters of the Coussot model are more troublesome.
As λ represents the microstructure in the Coussot model, λ 0 can be thought of as the initial degree of jamming caused by the microstructure (i.e. the structure that must be broken down to create fluid flow).As yet, there does not appear to be an accepted method to estimate λ 0 .We propose two methods evaluating λ 0 , and test these in the accompanying simulations.As discussed below, Method A is a simple analytical approach based on the critical stress, whereas Method B uses a graphical approach.
-Method A: Assuming all other parameters of the fluid are known, including the critical stress τ c , the initial condition, λ 0 , can be evaluated using the Coussot equation for the critical stress as Coussot et al. (2002a): Unfortunately parameter values for ω and α also do not have well-defined estimates in the literature, so herein we adjust these to ensure real solutions for λ 0 .However, in some simulations (see Section 6) this method appears to over-estimate shear stress.Furthermore, obtaining real solutions for λ 0 by perturbing α and ω can be time-consuming.-Method B: Our second approach (which is preferred) is to approximate the critical shear stress (τ c ) of a time-dependent fluid model using the maximum shear stress (τ max ) of a time-independent fluid model.This implies that on a graph of stress v. strain (τ : γ), the critical stress-strain point of the time-independent model should match the maximum stress point of the time-dependent model (i.e. the point where hysteresis causes the time-dependent model to operate along a different τ : γ curve).This point is labeled Q in Fig. 1.It is a relatively simple graphical trial and error exercise to adjust λ 0 , ω, and α to obtain the correct Q for a given T 0 , η 0 , and n.In this approach, the most important question is how to set the matching point, Q.In our avalanche model ( §6), the point Q is known because the critical shear stress is given in the experimental paper.However, for our debris flow model ( §7), only time-independent parameters are given in the corresponding experimental report.Thus, the matching point Q for this case was set where the maximum rate of strain of the thixotropic model was the same as the maximum rate of strain of the Herschel-Bulkley model.
The T 0 of the Coussot model in Eq. ( 10) can also be troublesome to estimate.This characteristic time for aging, which Coussot et al. (2002b) described as "spontaneous evolution of the microstructure," is not widely used and the literature does not provide insight on how to evaluate T 0 as a function of other rheological characteristics.Furthermore, T 0 has slightly different definitions by authors of several papers.Chanson et al. (2006) defined it as the characteristic time without any further measurement in their experiments, but provided another parameter, "rest time", used to set up the Bentonite suspensions in laboratory experiments in the result tables.However, Møller et al. (2006) defined T 0 as "the characteristic time of buildup of the microstructure at rest".Their characteristic time is close to the rest time of Chanson et al. (2006).Therefore, we make the assumption that the "rest time" measured in the Chanson's experiments is the same with the T 0 of Coussot for the thixotropic avalanche simulations (Section 6).For simulations of subaqueous debris flow (Section 7), the experiments did not report any time scales that could be used to estimate T 0 , so we included it as an unknown in the Method B described above.In independent parameters, which provides confidence that time-dependent and time-independent models are being compared in a reasonable manner.
4 Multi-material Level Set Method: Some types of debris flow, e.g.avalanches, can be reasonably modeled as a single fluid with a free surface where dynamics of the overlying fluid (in this example, air) are neglected.In contrast, subaqueous debris flows are more likely to require coupled modeling between lighter overlying water (Newtonian fluid) and heavier non-Newtonian debris.It is also possible to imagine more complex configurations where simultaneous solution of multiple debris layers or perhaps debris, water, and air might be necessary.For general purposes, it is convenient to apply a multi-material level set method so that any number of fluids with differing Newtonian and non-Newtonian properties can be considered.When only two fluids are considered, the multi-material level set method corresponds to the general level set method for two-phase flow.The level set method has a long history in multiphase fluids (Sussman et al., 1994;Chang et al., 1996;Sussman et al., 1998;Peng et al., 1999;Sussman and Fatemi, 1999;Bovet et al., 2010), and is based on using a φ i distance (level set) function to represent the distance of the i material (or material phase) from an interface with another material (Osher and Fedkiw, 2001).
The multi-material level set method herein follows Merriman et al. (1994) with the addition of high-order numerical schemes (Shi et al., 2002;Shu and Osher, 1989).The "level set", of the i-th fluid is designated as φ i , where where i = {1, 2, • • • , N m }, N m is the number of materials, Γ i is the interface of fluid i, and d is the distance from the interface.
The density and viscosity at a computational node for the multiple fluid system are evaluated from a combination of the individual fluid characteristics based on an approximate Heaviside function that provides a continuous transition over some distance on either side of an interface: where the Heaviside function for fluid i is where 2 is therefore the finite thickness of the numerical interface between fluids.
The level set initial condition is simply the distance from any grid point in the model to an initial set of interfaces, i.e. φ i = d i .
Note that each point has a distance to each i interface.The level set is treated as a conservatively-advected variable that evolves according to a simple non-diffusive transport equation (Osher and Fedkiw, 2001): The above is coupled to solution of momentum and continuity, Eqs. ( 2) and (3) to form a complete level-set solution for fluid flow.The continuous interface i at time t is located where φ i (x, t) = 0.In general, the i interface will be between the discrete grid points of the numerical solution, so it is found by multi-dimensional interpolation from the discrete φ i values.
After advancing the level set from φ(t) to φ(t + ∆t), the values of the level set will no longer satisfy the Eikonal condition of that is, the level set values on the grid cells obtained by solving Eq. ( 14) are no longer equidistant from the interface (i.e. the zero level set).It is known that if the level sets are naively evolved through time without satisfying the Eikonal condition the Heaviside functions will become increasingly inaccurate (Sussman et al., 1994).This problem is addressed with "reinitialization," which resets the φ(t + ∆t) to satisfy the Eikonal condition.The simplest approach to reinitialization is iterating an unsteady equation in pseudo-time to steady state such that the steady-state equation satisfies the Eikonal condition (Sussman et al., 1998).Let φ be an estimate of the reinitialized value for φ(t + ∆t) in the equation where T is the pseudo time, and S is the signed function as (Sussman et al., 1998) The time-advanced set of φ(t + ∆t) is the starting guess for φ, and the steady-state solution of φ will satisfy |∇ φi | = 1 to numerical precision.
For the present work, the advection term in Eq. ( 14) is discretized with the fifth-order WENO scheme, and the time derivative term is integrated by the third-order TVD Runge-Kutta method (Shu and Osher, 1989).For the reinitialization step of Eq. ( 15), the second-order ENO (Essentially Non-Oscillatory) scheme (Sussman et al., 1998) and a smoothing approach (Peng et al., 1999) are used for the spatial discretization (Jeon, 2015).

Poiseuille Flow of Bingham Fluid:
A two-dimensional Poiseuille flow in a channel driven by a steady pressure gradient of ∂p/∂x provides a validation case for the non-Newtonian fluid solver.If gravity is considered negligible and the flow is approximated as symmetric about a centerline between two walls, then the analytical solution for the flow on one side of the centerline is (Papanastasiou, 1987): where F is the distance from the center to a channel wall, y is the Cartesian axis normal to the flow direction with y = 0 at the centerline of the flow between the two walls, τ 0 is the yield stress, and F D is a length scale representing the relationship between yield stress and the pressure gradient: A convenient set of Bingham fluid parameters for the Poiseuille test cases can be extracted from the dip coating study of Filali et al. (2013) as shown in Table 1.In the simulations, the distance from the centerline to a side wall is 0.05 m.  3 where θ, d 0 , and l 0 represent the angle of a slope, the height of the initial avalanche that is normal to the slope, and the length of the avalanche along a slope, respectively.We modeled this same setup with our multi-material level-set solver.
The Chanson experiments identified four thixotropic flow types that were functions of the relative effect of initial structural jamming.Weak jamming (i.e.small λ 0 ) characterizes Type I, such that inertial effects dominate the downstream flow (highest Re) and the flow only ceases when it reaches the experiment outfall.It follows that Type I is effectively a model of an avalanche that propagates until it is stopped by an obstacle or change in slope.Type II flows had intermediate initial jamming, which showed rapid initial flow followed by deceleration until "restructuralization" that effectively stops the downstream progression.
Type II is a model of an avalanche that dissipates itself on the slope.The Type III flows, with the highest λ 0 , have complicated behavior with separation into identifiable packets of mass (typically two, but sometimes more) with different velocities.Type IV behavior was the extremum of zero flow.Chanson reported 28 experiments in total, but data on wave front propagation was provided for only five experiments (Fig. 6 in Ref. Chanson et al. (2006)) of Type I and II behavior.We simulated three of these experiments that covered a wide range of characteristics and behaviors, as shown in Table 2.Note that Chanson et al. (2006) used τ c2 to designate the critical shear stress during unloading (restructuralization), which we consider an approximation for the yield stress, τ 0 , for a time-independent model.We simulate the three cases of Table 2 with the time-independent Papanastasiou model of Eq. ( 6) and the Coussot timedependent model of Eqs. ( 8) and (9).For a Bingham time-independent model, we use n = 1 with K = η 0 from the Chanson experiments.The smoothing value of m = 100 was selected based on the Poiseuille flow modeled in Section 5, above.For a Herschel-Bulkley time-independent model, we use the same K and m as the Bingham plastic model, but with n = 1.1 that was used in the detailed technical report on the same experiments by Chanson et al. (2004).The time-dependent model requires specification of parameters {n, T 0 , α, η 0 , λ 0 , ω} as discussed in Section 3.2.The Herschel-Bulkley index in the time-dependent model uses the same value (n = 1.1) as the time-independent model.Two sets of values for {α, λ 0 , ω} are determined by the two methods (A and B) outlined in section 3.2, above.Method A uses Eq. ( 10), which requires a value for τ c ; herein this  is taken as Chanson's critical loading stress (τ c1 in Chanson et al., 2006) during the initial structural breakdown.Similarly, Method B requires a τ max for point Q in Fig. 1, which is also set to the critical loading stress.
For all simulations, the no-slip wall condition is applied to the bottom wall, and the number of computational cells is 512 × 80.The computational domain is rotated so the x axis is along the sloping bed, which means that computational cell faces are either orthogonal or parallel to the slope.The gravitational constant (g = 9.81 m • s −2 ) is divided into two components of (g sin θ, −g cos θ).The density and viscosity of air are 1.0 kg • m −3 and 1.0E-5 Pa • s, respectively.
The analytical relationships between shear stress and rate of strain for the different viscosity models are presented in Figs. 4 -6.In the figures, "Herschel-Bulkley" and "Bingham" lines are the results of Eq. ( 6) with n = 1.1 and n = 1.0, respectively.The "case A" and "case B" lines denote results of Methods A and B from Section 3.2 for determining time-dependent parameters with Eqs. ( 8) and ( 9).The estimated parameters of λ 0 , ω, and α by two methods that are used in Figs. 4 -6 are shown in   3.These figures illustrate the challenge of using Method A (the critical stress relationship) for estimating λ 0 .The numerical solutions of the Coussot model ordinary differential, Eq. ( 9), are obtained by the Runge-Kutta 4th-order method.The resulting time-dependent stress-strain relationship can be far from the time-independent relationship that is otherwise thought to be a reasonable model.
Propagation of the fluid wave front provides a simple means of directly comparing the temporal and spatial evolution of the model and experiments.To facilitate comparisons across experimental scales, the non-dimensionalized front location and time are x * = x/d 0 and t * = t g/d 0 , respectively.A simple theoretical estimate for the wave front propagation suitable for short time scales was derived from equations of motion as Eq. ( 26) in Chanson et al. (2006), repeated here as: The simulation, experiment, and theory results are shown in Figs. 7, 8, and 9 for Cases 1, 2, and 3 of Table 2, respectively.The 10 dashed line represents the theoretical solution for the propagating the front of Eq. ( 18).
The most striking feature in the results is that the simulations for Cases 2 and 3 (smaller τ 0 ) are relatively similar for all the models, whereas the time-independent models (Bingham and Herschel-Bulkley) completely fail for Case 1 (larger τ 0 ) even though the time-dependent models continue to perform reasonably well.The failure appears to be due to an inability of the time-independent models in Case 1 to develop sufficient strain to move out of the η = K γn−1 + mτ 0 regime that governs viscosity below the yield stress in Eq. ( 6).In contrast, the microstructural aging process that is inherent in Eqs. ( 8) and ( 9) allow the time-dependent models in Case 1 to develop reasonable flow conditions despite the higher τ 0 .No doubt the timeindependent models could be made to perform better in Case 1 by further manipulation of the model coefficients; however, our 5 approach was to use coefficients that could be set a priori based on data from the experiments and a plausible m value from Section 5.
We observe that the simplified theoretical front prediction from Eq. ( 18 consistently overpredicts the experimental front propagation in the early stages for Cases 1 and 2, but show better agreement with experiments than the simplified theory for t * > 4.However, for Case 3 (a Type I flow), the simplified theory is relatively poor, while the 2D simulations have good agreement up until t * ∼ 3, and then have significant underprediction of the experiments.As noted by Chanson et al. (2006), the Case 3 (Type I) experiments are at higher Reynolds numbers that, although theoretically laminar, may be transitioning to weakly turbulent.Because simplified theory of Eq. ( 18) is derived by neglecting inertia, it is not surprising that its performance degrades with increasing Reynolds number.
Although the simulations results have reasonable global agreement with experiments, on closer examination it can be seen that the 2D simulations predict a front movement that is initially too rapid in Type II flows (Case 1 and 2), and at longer times is too slow for Type I flows (Case 3).The challenge, of course, is the model error is integrative: if λ is wrong at a given time, then the dλ/dt will be wrong as well and the frontal position error will accumulate.Thus, an important issue for the time-dependent model appears to be selecting the appropriate values of {λ 0 , α, ω} that are consistent with experimentally-determined values of {η 0 , τ 0 , τ c , n, T 0 }.Although the more parsimonious time-independent model (with fewer parameters) performs reasonably well for our Case 2 and 3, it performs poorly in Case 1 and so should only be used with caution and careful calibration.
The above observations lead to a conclusion that the accelerative behaviors in the simulations and experiments are not well matched.The problem is shown most clearly in Fig. 7 for Case 1, where the experiments initially follow the acceleration implied by Eq. ( 18), but diverge with an inflection point and deceleration occurring somewhere near t * ∼ 4. In contrast, the models initially show a more rapid acceleration and an inflection point to deceleration at t   (2006) did not extend beyond t * ∼ 6.5, so it is impossible to know whether the experiments are showing an inflection point to deceleration at t * ∼ 5, but it seems likely given the results of the Case 1 and 2 studies.If there is an inflection point for Case 3, then it would appear that the consistent problem with the models is getting the correct transition from frontal acceleration to deceleration.To date, our experiments have not shown that we can significantly alter the model acceleration inflection points by altering parameters, which may indicate that there is a need to further consider the fundamental forms of the Coussot and Papanastasiou models when used for thixotropic flows.An alternative explanation may be that there are three-dimensional controls on the front propagation in the experiment that cannot be represented in the present 2D model.
7 Subaqueous Debris Flows: We have simulated subaqueous debris flows with our non-Newtonian fluid model and compared to experiments by Haza et al. (2013).We matched experimental cases with the largest density difference for the subaqueous debris to provide the largest effective negative buoyancy for the debris.The selected cases are 35% and 30% KCC (Kaolin Clay Content).The schematic design is shown in Fig. 10, with dimensions provided in Table 4.The gravitational constant for all simulations is g = 9.81 m • s −2 .
The simulation uses 340 × 100 rectangular cells.The no-slip wall boundary condition is applied to the bottom boundary.The computational domain is rotated so the x-axis is parallel to the slope, which allows the bottom to be represented as a straight surface without using cut grid cells or unstructured grids.This rotation also provides convenience in measuring the variables The parameters for the time-independent fluid model from Haza et al. (2013) are shown in Table 5.For all simulations, m = 100 for the exponential smoothing parameter is used based on results from Section 5, above.The parameters for the time-dependent fluid model are estimated from Method B in Section 3.2 and are shown in Table 6.The experiments did not report a rest time, so T 0 was set at a small positive value that provided a reasonable match the experiments.The analytical relationships between the shear stress and the rate of strain for the time-independent and the time-dependent fluid models are shown in Fig. 11 for Case 1 and Fig. 12 for Case 2.
Figure 13 provides a reference for measurements used to compare the model and experiments.These include the height of head-flow (H), the water depth at the front of head-flow (D), the run-out distance from the initial position (L), and the flow-front velocity (U ). Figure 14 shows evolution of the zero level sets for water (φ 2 ), which provides the continuous line separating the water from both the debris and the air.The time-independent viscosity-stress relationships that are often used for non-Newtonian flows are a subset of possible viscosity-stress models.We believe that more complex models may be necessary for real-world heterogeneous mixtures that include hysteresis in the stress/strain relationship as microstructure evolves with time.In particular, where a fluid at rest has a strongly jammed structure or undergoes restructuralization as the flow slows, the time-independent Bingham plastic and Herschel-Bulkley models will likely be inadequate.Unfortunately, the processes by which the initial jamming is locally overcome, and the processes through which a structure is recovered, are both poorly understood.For improved modeling of these flows, there is clearly a need for (1) more detailed experimental measurements during flow initiation and restructuralization, and (2) a better understanding of the relationship between measurable microstructure parameters and the effective viscosity-stress relationship. Nat ) where the superscript 'T ' indicates a matrix transpose.The η in the above is constant in time and uniform in space for a Newtonian fluid, but is potentially some nonlinear function of other flow variables for a non-Newtonian fluid.Note that common 4 Nat.Hazards Earth Syst.Sci.Discuss., https://doi.org/10.5194/nhess-2017-258Manuscript under review for journal Nat.Hazards Earth Syst.Sci. Discussion started: 17 July 2017 c Author(s) 2017.CC BY 4.0 License.incompressible Navier-Stokes solvers often use η as an eddy viscosity to model small-scale nonlinear advective effects as effective momentum diffusion, so the inclusion of nonlinearity and time-space dependency of η = f (u) is already a relatively common feature in hydrodynamic models that are ostensibly only for Newtonian fluids.That is, common Newtonian fluid solvers for turbulent flow can be used as a non-Newtonian solver as long as the viscosity is strictly a function of the velocity field and not other properties (e.g.density or pressure).For the non-Newtonian fluid models herein, both the time-independent and time-dependent methods use the local velocity rate-of-strain to update the plastic viscosity, η, as shown in Section 3, which 6) here m has dimension of time such that as m → ∞ we recover the original Herschel-Bulkley model with η → ∞, whereas m = 0 is a simple Newtonian fluid.The scalar value of the rate of strain is obtained from γ = 2 |II D | where II D is the Nat.Hazards Earth Syst.Sci.Discuss., https://doi.org/10.5194/nhess-2017-258Manuscript under review for journal Nat.Hazards Earth Syst.Sci. Discussion started: 17 July 2017 c Author(s) 2017.CC BY 4.0 License.

Figure 1 .
Figure 1.Concept of a graphical Method B for estimating a consistent set of λ0, α, ω parameters for Coussot model.
general, graphical Method B provides a simple means to estimate a consistent set of time-dependent parameters from the time-Nat.Hazards Earth Syst.Sci.Discuss., https://doi.org/10.5194/nhess-2017-258Manuscript under review for journal Nat.Hazards Earth Syst.Sci. Discussion started: 17 July 2017 c Author(s) 2017.CC BY 4.0 License.

Table 1 .
Bingham fluid Herschel-Bulkley model parameters used in Poiseulille flow test cases, from Filali et al.

Figure 2 .
Fig.2, the numerical results are in very good agreement with the analytical solutions for both values.For this simulation, the lower value of m = 100 s is reasonable for a Papanastasiou model.

Figure 3 .
Figure 3. Definition sketch for initial conditions of an avalanche along a slope.

Figure 7 .Figure 8 .Figure 9 .
Figure 7. Thixotropic avalanche Case 1: comparison of numerical results and experimental data for non-dimensional front displacement (x * ) as a function of non-dimensional time (t * )
slope (e.g.front distance, and water/mud thicknesses at the front.)These simulations include three fluids: mud, water, and air.The density of mud for each case is shown in Table6, and the densities of water and air are 1000.0kg • m −3 and 1.0 kg • m −3 , respectively.

Table 3 .
Parameters of the time-dependent fluid model for thixotropic avalanche simulations using Method A and Method B for setting values.