the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
Comparing thixotropic and Herschel–Bulkley parameterizations for continuum models of avalanches and subaqueous debris flows
ChanHoo Jeon
Ben R. Hodges
Avalanches and subaqueous debris flows are two cases of a wide range of natural hazards that have been previously modeled with nonNewtonian fluid mechanics approximating the interplay of forces associated with gravity flows of granular and solid–liquid mixtures. The complex behaviors of such flows at unsteady flow initiation (i.e., destruction of structural jamming) and flow stalling (restructuralization) imply that the representative viscosity–stress relationships should include hysteresis: there is no reason to expect the timescale of microstructure destruction is the same as the timescale of restructuralization. The nonNewtonian Herschel–Bulkley relationship that has been previously used in such models implies complete reversibility of the stress–strain relationship and thus cannot correctly represent unsteady phases. In contrast, a thixotropic nonNewtonian model allows representation of initial structural jamming and aging effects that provide hysteresis in the stress–strain relationship. In this study, a thixotropic model and a Herschel–Bulkley model are compared to each other and to prior laboratory experiments that are representative of an avalanche and a subaqueous debris flow. A numerical solver using a multimaterial levelset method is applied to track multiple interfaces simultaneously in the simulations. 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 (timedependent) fluid model shows reasonable agreement with all the experimental data. For most of the experimental conditions, the Herschel–Bulkley (timeindependent) model results were similar to the thixotropic model, a critical exception being conditions with a high yield stress where the Herschel–Bulkley model did not initiate flow. These results indicate that the thixotropic relationship is promising for modeling unsteady phases of debris flows and avalanches, but there is a need for better understanding of the correct material parameters and parameters for the initial structural jamming and characteristic time of aging, which requires more detailed experimental data than presently available.
A wide range of natural hazards involve gravitydriven flows down a slope, for example, landslides (terrestrial or submarine), flooddriven debris flows, mudflows, lahars, avalanches, and volcanic lava flows. Such flows range from relatively homogeneous particles (e.g., snow avalanches) to extremely heterogeneous particles (terrestrial landslides) and generally can be classified by solid concentration, material type, and mean velocity (Pierson and Costa, 1987; Smith and Lowe, 1991; Coussot and Meunier, 1996; Locat and Lee, 2002). Avalanches (e.g., snow, rock) are typically considered dry granular flows, whereas debris flows are liquid–solid mixtures where the solids are a dominant forcing, which can be contrasted to flood flows where sediment solids play a secondary role (Iverson, 1997). In theory, avalanche flows at the homogeneous end of the spectrum should be amenable to direct modeling as particles (granular flows), although it remains to be seen whether sufficient computer power can ever be practically applied for largescale natural hazards. Flows with heterogeneous mixtures of liquids and solids provide further challenges as we simply do not have an adequate and proven theory for representing their behavior at naturalhazard scales. Indeed, even if we develop a complete and practical theory for the movement of a mixture of fluid, particles, and entrained large objects across several magnitudes of scales, it is unclear how we would effectively capture the uncertainty associated with size and space distribution of solid objects (e.g., boulders in a landslide) that affect the flow propagation in any model attempting to directly represent fluidsolid structural interactions.
Largescale naturalhazard flows have been widely investigated with field observations, smallscale laboratory experiments, and numerical models. A common observation is that the complexity of the material composition and the effective rheological characteristics play important roles in material movement (Malet et al., 2003; Bisantino et al., 2010; Jeong, 2014; de Haas et al., 2015). This flow complexity is illustrated by the classification of subaqueous mass movements by Locat and Lee (2002) into five types with different behaviors: slides, topples, spreads, falls, and flows. At the “flow” end of the spectrum the water content is high, the particle sizes are small, and the flowing conditions are reasonably considered a fluid continuum. As the water content decreases and/or the particle size distribution covers more orders of magnitude, the theoretical basis for the fluid continuum approach becomes weaker and requires more empirical parameterization to capture other behaviors. Furthermore, the transition from a nonmoving to a flowing regime can involve spatial heterogeneity and timedependent behavior that is not wellrepresented by parameterizations of the flowing regime. Realworld debris flows include additional complexity as they erode and entrain material along the bottom and sides of the slope with the downstream flow. We take these issues as motivational for the present work and refer the reader to the recent review of Delannay et al. (2017) for further insight on granular flows and Shanmugam (2015) for heterogenous flows. The fundamentals physics of such flows is presented in Iverson (1997). Herein, we do not seek to distinguish between the differing physics of these various complex flows but rather focus on advancing the use of nonNewtonian viscosity models as a proxy for their general behavior. For simplicity in exposition, we will use the term “debris flow” to refer to any realworld mixture modeled as a continuum fluid using a nonNewtonian model.
Following Ancey (2007), the existing approaches to simulating debris flows can be categorized in three groups: (i) applying soil mechanics concept of coulomb behavior, which provides reasonable solutions for heterogeneous granular mass flows (Iverson and Denlinger, 2001; Iverson, 2003); (ii) merging soil and fluid mechanics models; and (iii) representing the heterogeneous debris as a continuum fluid with behaviors similar to a nonNewtonian fluid (the approach herein) where the transition from a stable structure to a moving fluid is handled as a viscous effect. This is not to imply that such flows are actually nonNewtonian fluids but merely that some of their behaviors can be captured with an appropriately parameterized viscosity model (e.g., Davies, 1986; Pierson and Costa, 1987; Coussot and Meunier, 1996; Pudasaini, 2012). Indeed, Iverson (2003) has referred to the rheological approach to debris flows as a “myth” and argued for its replacement with mixture models using separate solid–fluid components. However, their argument remains contentious, and it is not clear that the present state of mixture models is substantially less mythical than application of a rheological model when considering heterogenous mixtures over a wide range of scales. Given that debris flow covers such diverse phenomena and complex physics, it seems likely the “correct” model for the foreseeable future will be the model that best fits a specific event, experiment, or flow type of interest. In the absence of research that definitively solves the conundrum of debris flow, we follow the long history of using rheological models as a proxy. Such models are parsimonious in the number of coefficients and are effectively agnostic to the inherent uncertainties of fluidsolid distributions and interactions. In using a nonNewtonian rheological model, the realworld interaction between solid particles and surrounding fluid in a heterogeneous mixture can be thought of as similar to the microstructural behavior of a homogeneous nonNewtonian fluid where the local fluid viscosity is a function of the local stress. The main advantage of this approach is that a nonNewtonian rheological model is simply a time/spacedependent viscosity term for the Navier–Stokes equations. It follows that the time/spacevarying eddy viscosities in a wide range of existing hydrodynamic codes can be readily adapted to nonNewtonian behavior and used for parameterized modeling of debris flows.
Note that the terminology of nonNewtonian flows can be confusing as “timeindependent” models have viscosities that can change with both space and time throughout a flow. The difference between a “timeindependent” and a “timedependent” nonNewtonian fluid is whether the relation between stress and viscosity (i.e., nonNewtonian equation itself) is allowed to change with time. Thixotropic (timedependent) fluids are defined as nonNewtonian fluids where the process of “aging” during a flow changes the underlying fluid microstructure and the relationship between stress and viscosity (Moller et al., 2009). Herein, we examine how the use of a thixotropic model provides the ability to model behaviors that cannot be represented with a timeindependent nonNewtonian model. Our goal is to provide insight into the research needs for further experiments and model development into the natural hazards of gravitydriven debris flows across the transitions from inception to stalling.
Gravitydriven debris flows have a range of triggering mechanisms, and their composition evolves from initiation through motion and deposition or stalling, which can include a variety of behaviors that make modeling a challenge (Iverson, 1997). Parameterized nonNewtonian fluid models are an obvious approach to approximate these behaviors. Timeindependent rheological models have been widely used to simulate debris flows (e.g., Bovet et al., 2010; Pirulli, 2010; Tsai et al., 2011; Manga and Bonini, 2012); however the realworld flow characteristics include timedependent behaviors that could be categorized as “thixotropic” (Perret et al., 1996; Crosta and Dal Negro, 2003; Bagdassarov and Pinkerton, 2004; Aziz et al., 2010). Our focus in this paper is examining how a thixotropic model behavior compares to the more common timeindependent (Herschel–Bulkley) nonNewtonian fluid model.
From a macroscale perspective, debris flows have similar behaviors to “yieldstress fluids” that have been studied as a class of nonNewtonian fluids (Møller et al., 2006; Scotto di Santolo et al., 2010). A yieldstress fluid is effectively a solid (i.e., infinite viscosity) below a critical stress value (yield stress). This behavior is similar to what might be expected from a debris mixture of liquid and solids that is initially at rest and is triggered into motion as the yield stress is exceeded, which is the basis for prior timeindependent nonNewtonian models cited above. At the microscale under lowstress (nearrest) 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 higherviscosity 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 behaviors similar to a lowerviscosity fluid that deforms more easily under stress. This change from high viscosity to low viscosity under stress is readily simulated with a conventional timeindependent nonNewtonian Herschel–Bulkley model. Arguably, what is missing from a timeindependent model is that the destruction of the initial microstructure of the debris can change the effective macroscale viscosity and response to stress. If the flow stalls either globally or locally, it may take some time to reestablish its microstructure, so the yield stress for a recently stalled flow should be different than the yield stress after aging (consolidation). We can think of the behavior of a debris flow as controlled, at least partly, by the evolution of the microstructure and requiring a timedependent element in the nonNewtonian model.
The simplest nonNewtonian yieldstress fluids are Bingham plastics. More complex behaviors are associated with “shear thinning” and “shear thickening” where the effective viscosity nonlinearly changes with the rate of strain. For these standard cases, the relationship between viscosity and rate of strain is repeatable and timeindependent. The approach proposed by Herschel and Bulkley (1926) is a common approach for representing the general case of timeindependent nonNewtonian fluids wherein 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 $\dot{\mathit{\gamma}}$ is the scalar value of the rate of strain. The Herschel–Bulkley fluid index n controls the overall modeled behavior, where $\mathrm{0}<n<\mathrm{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 using a Herschel–Bulkley model is the viscosity is effectively infinite below the yield stress; i.e., the condition $\dot{\mathit{\gamma}}=\mathrm{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 illconditioned 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 biviscous 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 Sect. 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 “timeindependent” as the relationship between viscosity and shear is fixed by the selection of K, n, m, and τ_{0}. Arguably, there exist a wide range of debris flows over which the Papanastasiou approach should be adequate, as the timedependent characteristics of debris flows are, at least theoretically, principally confined to the initiation and 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 steadystate conditions for debris flows should be reasonably represented with timeindependent models. Indeed, O'Brien and Julien (1988) concluded, by their experiments, that mud flows whose volumetric sand concentration was 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 timeindependent 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 timeindependent 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 are focused on the steady or quasisteady conditions after the debris structure has (relatively) homogenized.
Thixotropic (timedependent) behavior, which is not represented in the Herschel–Bulkley model, provides an interesting avenue for representing the expected macroscale 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 shearthinning 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 renewal of structure (consolidation) occurs. For a real debris flow we expect different viscosity–stress behaviors during initiation, steadystate, and slowing phases (consistent with evolving microstructure), but a timeindependent Herschel–Bulkley model is effectively an assumption that the processes of destruction of microstructure and renewal are exactly reversible. 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 spatialtemporal destruction of microstructure that leads to changes in the effective viscosity that cannot be represented in the standard timeindependent models. Coussot et al. (2002a) proposed an empirical viscosity model for thixotropic fluids (presented in detail in Sect. 3, below), which captures these fundamental behaviors.
Prior research on thixotropic flows has mainly focused on laboratory experiments (Mohrig et al., 1999; Chanson et al., 2006; Sawyer et al., 2012; Haza et al., 2013), although a few studies have numerically investigated the characteristics of thixotropic flow on a simple inclined plane (Huynh et al., 2005; Hewitt and Balmforth, 2013). In general, numerical simulation results have not been well validated by the experimental data, arguably due to limitations in both nonNewtonian viscosity models and the sparsity of available laboratory data. Thixotropic flows modeled at the laboratory scale typically use clays (e.g., bentonite, kaolin) to create the microstructure controlling nonNewtonian behavior (Balmforth and Craster, 2001). Preparation of a homogenous clay suspension for such experiments is a demanding task, the details of which can be found in Coussot et al. (2002b), Huynh et al. (2005), and Chanson et al. (2006). Unfortunately, we cannot expect the structure of a heterogeneous largescale debris flow to mimic the flow scales, yield stresses, and parameters for a homogeneous thixotropic laboratory flow. However, lacking data from a largescale debris flow that could be adequately used for model comparisons, herein we take a first step by analyzing how thixotropic models compare to timeindependent models for laboratoryscale flows.
Validating the use of a nonNewtonian model to represent a realworld debris flow presents challenges on two levels: first, does the model correctly represent a nonNewtonian flow? Second, does the nonNewtonian flow (when parameterized) represent a realworld debris flow? To date, successful nonNewtonian models of realworld flows have been parameterized using a timeindependent approach, which limits the ability of the model to represent the transition phases, i.e., flow initiation and stalling (e.g., Bovet et al., 2010; Pirulli, 2010; Tsai et al., 2011; Manga and Bonini, 2012). Unfortunately, data on transition phases for realworld flows are lacking and are severely limited even for laboratoryscale flows.
In this paper we evaluate a timeindependent Papanastasiou model and a timedependent Coussot model for simulations of laboratoryscale avalanche and subaqueous debris flows, with comparisons to available experimental measurements. The governing equations are presented in Sect. 2, and the nonNewtonian Papanastasiou and Coussot viscosity models in Sect. 3. A key confounding issue for model–experiment comparisons is the estimation of parameters for a nonNewtonian fluid model (in particular the initial degree of jamming), which we discuss in Sect. 4. The numerical solver, using a multimaterial levelset method, is presented in Sect. 5. The solver is validated in Sect. 6 with the analytical solutions for the Poiseuille flow of a Bingham fluid. In Sect. 7 the solver is used to model a laboratory flow that is a reasonable proxy of a thixotropic avalanche. In Sect. 8 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 discuss the results and summarize conclusions in Sect. 9.
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:
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 nonNewtonian fluid.
The nonNewtonian fluid models herein use the local velocity rate of strain to update the plastic viscosity, η, as shown in Sect. 3, which makes the approach compatible with a wide range of numerical solvers that include a time/spacevarying eddy viscosity.
Equations (2) and (3) can be integrated over a control volume; by applying the Gauss divergence theorem, we obtain the basis for the common finitevolume numerical discretization (Ferziger and Perić, 2002). For simplicity in the present work, we limit ourselves to a twodimensional (2D) flow field for a downslope flow and the orthogonal (nearvertical) axis, which effectively assumes uniform flow in the crossstream axis. The external force term f represents the gravitational force only, neglecting surface tension forces and Coriolis. The advection term is discretized with the fifthorder WENO (weighted essentially nonoscillatory) scheme (Shi et al., 2002) or the secondorder TVD (total variation diminishing) Superbee scheme (Darwish and Moukalled, 2003) in separate numerical tests. The diffusion term on the righthand side of Eq. (3) is discretized with the secondorder central differencing scheme. The timederivative term for the momentum equations is integrated by the secondorder Crank–Nicolson implicit scheme. The deferredcorrection 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 firstorder 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 the biconjugate gradient squared method with preconditioning by the block Jacobi method). The developed code has been verified by the method of manufactured solutions (further details provided in Jeon, 2015).
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
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 $\dot{\mathit{\gamma}}=\mathrm{2}\sqrt{\left{\mathbf{II}}_{D}\right}$, where II_{D} is the 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 timeindependent.
In contrast, the timedependent (thixotropic) model of Coussot et al. (2002a) introduces dependency on a timevarying microstructure parameter (λ) in the general form
where η_{0} is the asymptotic viscosity at high shear rate, ω is a materialspecific 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 materialspecific parameter, and $\dot{\mathit{\gamma}}$ 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 values of α require greater microstructure homogenization (smaller λ) to drive the system to steadystate conditions ($\mathrm{d}\mathit{\lambda}/\mathrm{d}t\to \mathrm{0}$).
The timedependent Coussot model requires parameters for the asymptotic viscosity (η_{0}), Herschel–Bulkley fluid index (n), characteristic time (T_{0}), and two materialspecific parameters (ω and α) that control the response (destruction) of the microstructure. Additionally, an initial condition for λ_{0} is required to solve the ordinary differential equation presented as Eq. (9). The parameters η_{0} and n are easily obtained from the timeindependent 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):
$$\begin{array}{}\text{(10)}& {\mathit{\tau}}_{\mathrm{c}}={\displaystyle \frac{{\mathit{\eta}}_{\mathrm{0}}\left(\mathrm{1}+\mathit{\omega}{\mathit{\lambda}}_{\mathrm{0}}^{n}\right)}{\mathit{\alpha}{T}_{\mathrm{0}}{\mathit{\lambda}}_{\mathrm{0}}}}.\end{array}$$Unfortunately parameter values for ω and α also do not have welldefined estimates in the literature, so herein we adjust these to ensure real solutions for λ_{0}. However, in some simulations (see Sect. 7) this method appears to overestimate shear stress. Furthermore, obtaining real solutions for λ_{0} by perturbing α and ω can be timeconsuming.

Method B: our second approach (which is preferred) is to approximate the critical shear stress (τ_{c}) of a timedependent fluid model using the maximum shear stress (τ_{max}) of a timeindependent fluid model. This implies that, on a graph of stress vs. strain ($\mathit{\tau}:\dot{\mathit{\gamma}}$), the critical stress–strain point of the timeindependent model should match the maximum stress point of the timedependent model (i.e., the point where hysteresis causes the timedependent model to operate along a different $\mathit{\tau}:\dot{\mathit{\gamma}}$ 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 (Sect. 7), the point Q is known because the critical shear stress is given in the experimental paper. However, for our debrisflow model (Sect. 8), only timeindependent 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 they 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 (Sect. 7). For simulations of subaqueous debris flow (Sect. 8), the experiments did not report any timescales that could be used to estimate T_{0}, so we included it as an unknown in the method B described above. In general, the graphical method B provides a simple means to estimate a consistent set of timedependent parameters from the timeindependent parameters, which provides confidence that timedependent and timeindependent models are being compared in a reasonable manner.
Some types of debris flow, such as 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 nonNewtonian 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 multimaterial levelset method so that any number of fluids with differing Newtonian and nonNewtonian properties can be considered. When only two fluids are considered, the multimaterial levelset method corresponds to the general levelset method for twophase flow. The levelset 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 multimaterial levelset method herein follows Merriman et al. (1994) with the addition of highorder numerical schemes (Shu and Osher, 1989; Shi et al., 2002). The “level set” of the ith fluid is designated as ϕ_{i}:
where $i=\mathit{\{}\mathrm{1},\mathrm{2},\mathrm{\cdots},{N}_{\mathrm{m}}\mathit{\}}$, 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 multiplefluid 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 levelset 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 nondiffusive transport equation (Osher and Fedkiw, 2001):
The above is coupled to a solution of momentum and continuity, Eqs. (2) and (3), to form a complete levelset solution for fluid flow. The continuous interface i at time t is located where ${\mathit{\varphi}}_{i}(\mathit{x},t)=\mathrm{0}$. In general, the i interface will be between the discrete grid points of the numerical solution, so it is found by multidimensional 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 $\left\mathrm{\nabla}{\mathit{\varphi}}_{i}\right=\mathrm{1}$; that is, the levelset 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 pseudotime to steady state such that the steadystate equation satisfies the eikonal condition (Sussman et al., 1998). Let $\widehat{\mathit{\varphi}}$ be an estimate of the reinitialized value for ϕ(t+Δt) in the equation
where 𝒯 is the pseudotime, and S is the signed function as (Sussman et al., 1998)
The timeadvanced set of ϕ(t+Δt) is the starting guess for $\widehat{\mathit{\varphi}}$, and the steadystate solution of $\widehat{\mathit{\varphi}}$ will satisfy $\left\mathrm{\nabla}{\widehat{\mathit{\varphi}}}_{i}\right=\mathrm{1}$ to numerical precision.
For the present work, the advection term in Eq. (14) is discretized with the fifthorder WENO scheme, and the timederivative term is integrated by the thirdorder TVD Runge–Kutta method (Shu and Osher, 1989). For the reinitialization step of Eq. (15), the secondorder ENO (essentially nonoscillatory) scheme (Sussman et al., 1998) and a smoothing approach (Peng et al., 1999) are used for the spatial discretization (further details are provided in Jeon, 2015).
A twodimensional Poiseuille flow in a channel driven by a steady pressure gradient of $\partial p/\partial x$ provides a validation case for the nonNewtonian 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 dipcoating 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. Our model grid uses 320 cells in the flow direction and 32 cells in the crossstream direction. A Neumann boundary condition is applied along the lower boundary of the simulation domain, so the simulation includes only the upper halfchannel of this symmetric flow.
Using the Papanastasiou model of Eq. (6) to approximate a Herschel–Bulkley model of a Bingham fluid requires timescale parameter m to provide smooth behavior across the yieldstress threshold. We tested values of $m=\left\{\mathrm{100},\mathrm{400}\right\}$ s. As shown in 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.
An avalanche is a granular flow of an initially solid field that is triggered from rest into a downslope flow. A thixotropic model of an avalanche as a fluid continuum can represent a rapid progression from local to global release of the initial structural jamming, λ_{0}. Chanson et al. (2006) developed dambreak experiments with a thixotropic fluid that provide reasonable facsimiles of avalanche flows if the timescale to remove the dam is smaller than the timescale for release of structural jamming. The initial conditions of the Chanson experiments are shown in Fig. 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 multimaterial levelset 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”, which 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 were provided for only five experiments (Fig. 6 in 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 timeindependent model.
We simulate the three cases of Table 2 with the timeindependent Papanastasiou model of Eq. (6) and the timedependent Coussot model of Eqs. (8) and (9). For a timeindependent Bingham 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 Sect. 6, above. For a timeindependent Herschel–Bulkley model, we use the same K and m as the Bingham plastic model, but with n=1.1 as was used in the detailed technical report on the same experiments by Chanson et al. (2004). The timedependent model requires specification of parameters $\left\{n,\phantom{\rule{0.125em}{0ex}}{T}_{\mathrm{0}},\phantom{\rule{0.125em}{0ex}}\mathit{\alpha},\phantom{\rule{0.125em}{0ex}}{\mathit{\eta}}_{\mathrm{0}},\phantom{\rule{0.125em}{0ex}}{\mathit{\lambda}}_{\mathrm{0}},\phantom{\rule{0.125em}{0ex}}\mathit{\omega}\right\}$ as discussed in Sect. 4. The Herschel–Bulkley index in the timedependent model uses the same value (n=1.1) as the timeindependent model. Two sets of values for $\left\{\mathit{\alpha},\phantom{\rule{0.125em}{0ex}}{\mathit{\lambda}}_{\mathrm{0}},\phantom{\rule{0.125em}{0ex}}\mathit{\omega}\right\}$ are determined by the two methods (A and B) outlined in Sect. 4, 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 noslip 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 (gsin θ, −gcos θ). The density and viscosity of air are 1.0 kg m^{−3} and 1.0 × 10^{−5} Pa s, respectively.
The analytical relationships between shear stress and rate of strain for the different viscosity models are presented in Figs. 4 through 6. In these 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 Sect. 4 for determining timedependent parameters with Eqs. (8) and (9). The estimated parameters of λ_{0}, ω, and α by two methods that are used in these figures are shown in Table 3. These results illustrate the challenge of using method A (the critical stress relationship) for estimating λ_{0}. The numerical solutions of the Coussot model ordinary differential equation, Eq. (9), are obtained by the Runge–Kutta fourthorder method. The resulting timedependent stress–strain relationship can be far from the timeindependent 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 nondimensionalized front location and simulation time after gate opening are ${x}^{\ast}=x/{d}_{\mathrm{0}}$ and ${t}^{\ast}=t\sqrt{g/{d}_{\mathrm{0}}}$, respectively. A simple theoretical estimate for the wave front propagation suitable for short timescales 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 dashed line represents the theoretical solution for 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 timeindependent models (Bingham and Herschel–Bulkley) completely fail for case 1 (larger τ_{0}) even though the timedependent models continue to perform reasonably well. The failure appears to be due to an inability of the timeindependent models in case 1 to develop sufficient strain to move out of the $\mathit{\eta}=K{\dot{\mathit{\gamma}}}^{n\mathrm{1}}+m{\mathit{\tau}}_{\mathrm{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) allows the timedependent 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 approach was to use coefficients that could be set a priori based on data from the experiments and a plausible m value from Sect. 6.
We observe that the simplified theoretical front prediction from Eq. (18), the dashed line in the figures, is a good representation of Chanson's type II flows (case 1 and case 2) up until ${t}^{\ast}\sim \mathrm{3}$, but it diverges rapidly thereafter. Our 2D simulations consistently overpredict 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}^{\ast}>\mathrm{4}$. However, for case 3 (type I flow), the simplified theory is relatively poor, while the 2D simulations have good agreement up until ${t}^{\ast}\sim \mathrm{3}$ and then show 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 the 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 (cases 1 and 2) and at longer times is too slow for type I flows (case 3). The challenge, of course, is that 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 timedependent model appears to be selecting the appropriate values of $\left\{{\mathit{\lambda}}_{\mathrm{0}},\mathit{\alpha},\mathit{\omega}\right\}$ that are consistent with experimentally determined values of $\left\{{\mathit{\eta}}_{\mathrm{0}},{\mathit{\tau}}_{\mathrm{0}},{\mathit{\tau}}_{\mathrm{c}},n,{T}_{\mathrm{0}}\right\}$. Although the more parsimonious timeindependent model (with fewer parameters) performs reasonably well for our cases 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}^{\ast}\sim \mathrm{4}$. In contrast, the models initially show a more rapid acceleration and an inflection point to deceleration at ${t}^{\ast}\sim \mathrm{1}$. Interestingly, the simulated front locations in cases 1 and 2 are not unreasonable predictions for ${t}^{\ast}>\mathrm{4}$, but they get there along slightly different paths than the experiments. The case 3 (type I) models show different behaviors: they perform quite well for ${t}^{\ast}\le \mathrm{3}$ and then show deceleration at the same time as the experiment appears to be accelerating. Unfortunately, the experiments of Chanson et al. (2006) did not extend beyond ${t}^{\ast}\sim \mathrm{6.5}$, so it is impossible to know whether the experiments are showing an inflection point to deceleration at ${t}^{\ast}\sim \mathrm{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 threedimensional controls on the front propagation in the experiment that cannot be represented in the present 2D model.
In general, subaqueous debris flows are heterogenous gravity flows where the interaction of the overlying water with the downslope flow of the debris has a significant effect on momentum. Such flows are expected to be qualitatively similar to the subaqueous mud flow examined in the laboratory by Haza et al. (2013). We conducted simulations matching the Haza et al. (2013) experimental cases with the largest density difference between the water and mud. These conditions provide the largest effective negative buoyancy for the debris and minimize effect of turbidity. 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 noslip 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 normal to the 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 Table 6, and the densities of water and air are 1000.0 and 1.0 kg m^{−3}, respectively.
The parameters for the timeindependent 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 Sect. 6, above. The parameters for the timedependent fluid model are estimated from method B in Sect. 4 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 to the experiments. The analytical relationships between the shear stress and the rate of strain for the timeindependent and the timedependent 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 runout distance from the initial position (L), and the flowfront 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. Figures 15 and 16 show the evolution of the runout distance (L) for case 1 and case 2, respectively. It can be seen that both timeindependent and timedependent models are reasonable approximations of the limited experimental data. Both types of models appear to underestimate the initial runout and slightly overestimate later times.
Figures 17 and 18 show a comparison of H and D for simulations and experiments. Again, within the limited availability of experimental data, both timeindependent and timedependent models provide reasonable agreement. Figures 19 and 20 show similar agreement for the front velocities, although the experimental data are insufficient to validate the wavelike oscillation of the velocity in the simulations.
These results indicate the multimaterial levelset model is capable of representing the key features in a subaqueous debris flow. For this flow, the use of the simpler timeindependent viscosity model seems justified, although this is likely a function of the experimental conditions. An important limitation of the tested subaqueous debris flows is that they do not have the restructuralization in the downstream flow or the strongly jammed initial structure seen in the experiments of Chanson et al. (2006)
This work shows that a multiphase flow solver using a multimaterial levelset method with yieldstress models of nonNewtonian viscosity provides a means for numerical approximation of avalanches and subaqueous debris flows. This simulation approach was tested with both timeindependent (Herschel–Bulkley, Papanastasiou, Bingham plastic) and timedependent (thixotropic Coussot) models of viscosity, which are implemented using continuum mechanics solutions for multiple fluids. A key problem is that the Coussot model requires more parameters than the timeindependent fluid models, but available experimental data are insufficient to definitively set parameter values. To resolve this issue, two different approaches were used to evaluating the Coussot parameters. Overall, the numerical results showed reasonable agreement with prior experimental data.
Although stress–strain relationships indicate the timedependent approach provides the hysteresis that is desirable in a debrisflow model, in comparisons with experimental data the timedependent Coussot model provides a clear advantage for only for a single case – where the Herschel–Bulkley and Bingham plastic models erroneously predicted nearzero flow. Nevertheless, much of the complexity in realworld behavior for debris mixtures is due to interactions across spatial scales for heterogeneous mixtures, which leads to significantly different stress–strain relationships during structural breakdown and restructuralization that should require a timedependent model. Unfortunately, for experimental simplicity most researchers expend significant effort to create a homogeneous mixture as an initial condition for a debris flow, and the extent to which the structural breakdown results in temporary heterogeneous scales is unknown. Existing laboratory data do not provide sufficiently detailed insight into the processes controlling destruction of jamming or the restructuralization of the flow, which leaves significant uncertainty in specification of the correct parameters.
The timeindependent viscosity–stress relationships that are often used for nonNewtonian flow models of natural hazards are a subset of possible viscosity–stress models. We believe that more complex models may be necessary for realworld 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 timeindependent 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 the structure is recovered, are both poorly understood. For timedependent thixotropic models to be useful in modeling realworld avalanches and debris flows, there is a need for a consistent approach to defining the initial jamming (λ_{0}), the characteristic time of aging (T_{0}), and the asymptotic shear viscosity (η_{0}), along with the material parameters ω and α for realworld systems. As yet, these parameters are not well defined for either simple laboratory models or complex realworld flows. To improve our understanding of the thixotropic model, there is a need for a comprehensive sensitivity analysis of these five driving parameters for the expected scales of realworld systems (which are as yet unknown). Furthermore, with or without the thixotropic model, 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 stress–strain relationship. The present work shows that a timedependent (thixotropic) viscosity model may be an effective proxy for representing the inception and stalling of an avalanche or debris flow, but much work remains to be done before realworld natural hazards can be modeled in this manner.
The experimental data in the figures can be found in Chanson et al. (2006) and Haza et al. (2013). The simulation data can be obtained from the corresponding author (ChanHoo Jeon).
The authors declare that they have no conflict of interest.
The authors acknowledge the Texas Advanced Computing Center (TACC) at The
University of Texas (UT) at Austin for providing HPC, visualization, database, or
grid resources that have contributed to the research results reported within
this paper (http://www.tacc.utexas.edu). Publication support was provided
by the Center for Water and the Environment at UT Austin and the Carl Ernest
and Mattie Ann Muldrow Reistle Jr. Centennial Fellowship in Engineering.
Edited by: JeanPhilippe Malet
Reviewed by: two anonymous referees
Ancey, C.: Plasticity and geophysical flows: A review, J. NonNewton. Fluid, 142, 4–35, 2007. a
Ancey, C. and Cochard, S.: The dambreak problem for HerschelBulkley viscoplastic fluids down steep flumes, J. NonNewton. Fluid, 158, 18–35, 2009. a
Aziz, M., Towhata, I., Yamada, S., Qureshi, M. U., and Kawano, K.: Waterinduced granular decomposition and its effects on geotechnical properties of crushed soft rocks, Nat. Hazards Earth Syst. Sci., 10, 1229–1238, https://doi.org/10.5194/nhess1012292010, 2010. a
Bagdassarov, N. and Pinkerton, H.: Transient phenomena in vesicular lava flows based on laboratory experiments with analogue materials, J. Volcanol. Geoth. Res., 132, 115–136, 2004. a
Balay, S., Abhyankar, S., Adams, M. F., Brown, J., Brune, P., Buschelman, K., Dalcin, L., Eijkhout, V., Gropp, W. D., Kaushik, D., Knepley, M. G., McInnes, L. C., Rupp, K., Smith, B. F., Zampini, S., Zhang, H., and Zhang, H.: PETSc Web page, available at: http://www.mcs.anl.gov/petsc, last access: 5 December 2016. a
Balmforth, N. J. and Craster, R. V.: Geomorphological Fluid Mechanics, in: Lecture Notes in Physics, vol. 582, Springer, Berlin, Heidelberg, Germany, 2001. a
Balmforth, N. J., Craster, R. V., Rust, A. C., and Sassi, R.: Viscoplastic flow over an inclined surface, J. NonNewton. Fluid, 142, 219–243, 2007. a
Beverly, C. R. and Tanner, R. I.: Numerical analysis of threedimensional Bingham plastic flow, J. NonNewton. Fluid, 42, 85–115, 1992. a, b
Bingham, E. C.: An investigation of the laws of plastic flow, U.S. Bur. Stand. Bullet., 13, 309–353, 1916. a
Bisantino, T., Fischer, P., and Gentile, F.: Rheological characteristics of debrisflow material in SouthGargano watersheds, Nat. Hazards, 54, 209–223, 2010. a
Bovet, E., Chiaia, B., and Preziosi, L.: A new model for snow avalanche dynamics based on nonNewtonian fluids, Meccanica, 45, 753–765, 2010. a, b, c, d
Chang, Y. C., Hou, T. Y., Merriman, B., and Osher, S.: A Level Set Formulation of Eulerian Interface Capturing Methods for Incompressible Fluid Flows, J. Comput. Phys., 124, 449–464, 1996. a
Chanson, H., Coussot, P., Jarny, S., and Toquer, L.: A Study of Dam Break Wave of Thixotropic Fluid: Bentonite Surges down an Inclined Plane, Tech. Rep. CH54/04, Department of Civil Engineering, The University of Queensland, St. Lucia QLD, Australia, 2004. a
Chanson, H., Jarny, S., and Coussot, P.: Dam break wave of thixotropic fluid, J. Hydraul. Eng., 132, 280–293, 2006. a, b, c, d, e, f, g, h, i, j, k, l, m
Coussot, P. and Meunier, M.: Recognition, classification and mechanical description of debris flows, EarthSci. Rev., 40, 209–227, 1996. a, b
Coussot, P., Nguyen, Q. D., Huynh, H. T., and Bonn, D.: Avalanche behavior in yield stress fluids, Phys. Rev. Lett., 88, 175501, https://doi.org/10.1103/PhysRevLett.88.175501, 2002a. a, b, c
Coussot, P., Nguyen, Q. D., Huynh, H. T., and Bonn, D.: Viscosity bifurcation in thixotropic, yielding fluids, J. Rheol., 46, 573–589, 2002b. a, b
Crosta, G. B. and Dal Negro, P.: Observations and modelling of soil slipdebris flow initiation processes in pyroclastic deposits: the Sarno 1998 event, Nat. Hazards Earth Syst. Sci., 3, 53–69, https://doi.org/10.5194/nhess3532003, 2003. a
Darwish, M. S. and Moukalled, F.: TVD schemes for unstructured grids, Int. J. Heat Mass Tran., 46, 599–611, 2003. a
Davies, T. R. H.: Large debris flows: a macroviscous phenomenon, Acta Mech., 63, 161–178, 1986. a
De Blasio, F. V., Engvik, L., Harbitz, C. B., and Elverhøi, A.: Hydroplaning and submarine debris flows, J. Geophys. Res., 109, C01002, https://doi.org/10.1029/2002JC001714, 2004. a
de Haas, T., Braat, L., Leuven, J. R. F. W., Lokhorst, I. R., and Kleinhans, M. G.: Effects of debris flow composition on runout, depositional mechanisms, and deposit morphology in laboratory experiments, J. Geophys. Res.Earth, 120, 1949–1972, 2015. a
Delannay, R., Valance, A., Mangeney, A., Roche, O., and Richard, P.: Granular and particleladen flows: from laboratory experiments to field observations, J. Phys. DAppl. Phys., 50, 053001, https://doi.org/10.1088/13616463/50/5/053001, 2017. a
Dent, J. D. and Lang, T. E.: A viscous modified Bingham model of snow avalanche motion, Ann. Glaciol., 4, 42–46, 1983. a
Ferziger, J. H. and Perić, M.: Computational methods for fluid dynamics, 3rd edn., Springer, Berlin, Heidelberg, Germany, 2002. a, b, c
Filali, A., Khezzar, L., and Mitsoulis, E.: Some experiences with the numerical simulation of Newtonian and Bingham fluids in dip coating, Comput. Fluids, 82, 110–121, 2013. a, b
Guermond, J. L., Minev, P., and Shen, J.: An overview of projection methods for incompressible flows, Comput. Meth. Appl. Mech. Engrg., 195, 6011–6045, 2006. a
Haza, Z. F., Harahap, I. S. H., and Dakssa, L.: Experimental studies of the flowfront and drag forces exerted by subaqueous mudflow on inclined base, Nat. Hazards, 68, 587–611, 2013. a, b, c, d, e
Herschel, W. H. and Bulkley, R.: Konsistenzmessungen von GummiBenzollösungen, Kolloid Zeitschrift, 39, 291–300, 1926. a
Hewitt, D. R. and Balmforth, N. J.: Thixotropic gravity currents, J. Fluid Mech., 727, 56–82, 2013. a
Huynh, H. T., Roussel, N., and Coussot, P.: Ageing and free surface flow of a thixotropic fluid, Phys. Fluids, 17, 033101, https://doi.org/10.1063/1.1844911, 2005. a, b
Iverson, R. M.: The physics of debris flows, Rev. Geophys., 35, 245–296, 1997. a, b, c
Iverson, R. M.: The debrisflow rheology myth, in: The 3rd International Conference on DebrisFlow Hazard Mitigation: Mechanics, Prediction, and Assessment, edited by: Rickenmann, D. and Chen, C.L., 303–314, Mills Press, Davos, Switzerland, 2003. a, b
Iverson, R. M. and Denlinger, R. P.: Flow of variably fluidized granular masses across threedimension terrain: 1. Coulomb mixture theory, J. Geophys. Res., 106, 537–552, 2001. a
Jeon, C.H.: Modeling of Debris Flows and Induced Phenomena with NonNewtonian Fluid Models, PhD dissertation, The University of Texas at Austin, Austin, Texas, USA, 2015. a, b
Jeong, S. W.: The Effect of Grain Size on the Viscosity and Yield Stress of FineGrained Seiments, J. Mt. Sci., 11, 31–40, 2014. a
Liu, K. F. and Mei, C. C.: Slow spreading of a sheet of Bingham fluid on an inclined plane, J. Fluid Mech., 207, 505–529, 1989. a
Liu, W. and Zhu, K.Q.: A study of startup flow of thixotropic fluids including inertia effects on an inclined plane, Phys. Fluids, 23, 013103, https://doi.org/10.1063/1.3536654, 2011. a
Locat, J. and Lee, H. J.: Submarine landslides: advances and challenges, Can. Geotech. J., 39, 193–212, 2002. a, b
Malet, J. P., Remaître, A., Maquaire, O., Ancey, C., and Locat, J.: Flow susceptibility of heterogeneous marly formations: implications for torrent hazard control in the Barcelonnette Basin (AlpesdeHauteProvence, France), in: The 3rd International Conference on DebrisFlow Hazard Mitigation: Mechanics, Prediction, and Assessment, edited by: Rickenmann, D. and Chen, C.L., 351–362, Mills Press, Davos, Switzerland, 2003. a
Manga, M. and Bonini, M.: Large historical eruptions at subaerial mud volcanoes, Italy, Nat. Hazards Earth Syst. Sci., 12, 3377–3386, https://doi.org/10.5194/nhess1233772012, 2012. a, b
Mei, C. C.: Lecture notes on fluid dynamics, Massachusetts Institute of Technology (MIT), Cambridge, Massachusetts, USA, 2007. a
Merriman, B., Bence, J. K., and Osher, S. J.: Motion of multiple junctions: a level set approach, J. Comput. Phys., 112, 334–363, 1994. a
Mohrig, D., Elverhøi, A., and Parker, G.: Experiments on the relative mobility of muddy subaqueous and subaerial debris flows and their capacity to remobilize antecedent deposits, Mar. Geol., 154, 117–129, 1999. a
Moller, P., Fall, A., Chikkadi, V., Derks, D., and Bonn, D.: An attempt to categorize yield stress fluid behaviour, Philos. T. Roy. Soc. A, 367, 5139–5155, 2009. a
Møller, P. C. F., Mewis, J., and Bonn, D.: Yield stress and thixotropy: on the difficulty of measuring yield stresses in practice, Soft Matter, 2, 274–283, 2006. a, b
Murthy, J. Y. and Mathur, S.: Periodic flow and heat transfer using unstructured meshes, Int. J. Numer. Meth. Fl., 25, 659–677, 1997. a
O'Brien, J. S. and Julien, P. Y.: Laboratory analysis of mudflow properties, J. Hydraul. Eng., 114, 877–887, 1988. a
Osher, S. and Fedkiw, R. P.: Level set methods: An Overview and Some Recent Results, J. Comput. Phys., 169, 463–502, 2001. a, b
Papanastasiou, T. C.: Flow of materials with yield, J. Rheol., 31, 385–404, 1987. a, b, c
Peng, D., Merriman, B., Osher, S., Zhao, H., and Kang, M.: A PDEbased Fast Local Level Set Method, J. Comput. Phys., 155, 410–438, 1999. a, b
Perret, D., Locat, J., and Martignoni, P.: Thixotropic behavior during shear of a finegrained mud from Eastern Canada, Eng. Geol., 43, 31–44, 1996. a
Pierson, T. and Costa, J.: A rheologic classification of subaerial sedimentwater flows, Rev. Eng. Geol., VII, 1–12, 1987. a, b
Pignon, F., Magnin, A., and Piau, J.: Thixotropic colloidal suspensions and flow curves with minimum: identification of flow regimes and rheometric consequences, J. Rheol., 40, 573–587, 1996. a
Pirulli, M.: On the use of the calibrationbased approach for debrisflow forwardanalyses, Nat. Hazards Earth Syst. Sci., 10, 1009–1019, https://doi.org/10.5194/nhess1010092010, 2010. a, b
Pudasaini, S. P.: A general twophase debris flow model, J. Geophy. Res., 117, F03010, https://doi.org/10.1029/2011JF002186, 2012. a
Sawyer, D. E., Flemings, P. B., Buttles, J., and Mohrig, D.: Mudflow transport behavior and deposit morphology: Role of shear stress to yield strength ratio in subaqueous experiments, Mar. Geol., 307–310, 28–39, 2012. a
Scotto di Santolo, A., Pellegrino, A. M., and Evangelista, A.: Experimental study on the rheological behaviour of debris flow, Nat. Hazards Earth Syst. Sci., 10, 2507–2514, https://doi.org/10.5194/nhess1025072010, 2010. a
Shanmugam, G.: The landslide problem, J. Palaeogeogr., 4, 109–166, 2015. a
Shi, J., Hu, C., and Shu, C.W.: A technique of treating negative weights in WENO schemes, J. Comput. Phys., 175, 108–127, 2002. a, b
Shu, C. W. and Osher, S.: Efficient implementation of essentially nonoscillatory shock capturing schemes II, J. Comput. Phys., 83, 32–78, 1989. a, b
Smith, G. A. and Lowe, D. R.: Lahars: volcanohydrologic events and deposition in the debris flow – hyperconcentrated flow continuum, in: Sedimentation in Volcanic Settings, edited by: Fisher, R. V. and Smith, G. A., vol. 45, 59–70, SEPM Special Publication, Tulsa, Oklahoma, USA, 1991. a
Sussman, M. and Fatemi, E.: An efficient, interfacepreserving level set redistancing algorithm and its application to interfacial incompressible fluid flow, SIAM J. Sci. Comput., 20, 1165–1191, 1999. a
Sussman, M., Smereka, P., and Osher, S.: A level set approach for computing solutions to incompressible twophase flow, J. Comput. Phys., 114, 146–159, 1994. a, b
Sussman, M., Fatemi, E., Smereka, P., and Osher, S.: An improved level set method for incompressible twophase flows, Comput. Fluids, 27, 663–680, 1998. a, b, c, d
Tsai, M. P., Hsu, Y. C., Li, H. C., Shu, H. M., and Liu, K. F.: Application of simulation technique on debris flow hazard zone delineation: a case study in the Daniao tribe, Eastern Taiwan, Nat. Hazards Earth Syst. Sci., 11, 3053–3062, https://doi.org/10.5194/nhess1130532011, 2011. a, b
 Abstract
 Introduction
 Governing equations
 NonNewtonian fluid models
 Estimation of parameters for timedependent Coussot model
 Multimaterial levelset method
 Poiseuille flow of Bingham fluid
 Thixotropic avalanches
 Subaqueous debris flows
 Discussion and conclusions
 Data availability
 Competing interests
 Acknowledgements
 References
 Abstract
 Introduction
 Governing equations
 NonNewtonian fluid models
 Estimation of parameters for timedependent Coussot model
 Multimaterial levelset method
 Poiseuille flow of Bingham fluid
 Thixotropic avalanches
 Subaqueous debris flows
 Discussion and conclusions
 Data availability
 Competing interests
 Acknowledgements
 References