Journal topic
Nat. Hazards Earth Syst. Sci., 18, 303–319, 2018
https://doi.org/10.5194/nhess-18-303-2018
Nat. Hazards Earth Syst. Sci., 18, 303–319, 2018
https://doi.org/10.5194/nhess-18-303-2018

Research article 22 Jan 2018

Research article | 22 Jan 2018

# Comparing thixotropic and Herschel–Bulkley parameterizations for continuum models of avalanches and subaqueous debris flows

Comparing thixotropic and Herschel–Bulkley parameterizations for continuum models of avalanches and subaqueous debris flows
Chan-Hoo Jeon1,2 and Ben R. Hodges2 Chan-Hoo Jeon and Ben R. Hodges
• 1Division of Marine Science, The University of Southern Mississippi, 1020 Balch Blvd, Stennis Space Center, Mississippi 39529, USA
• 2Center for Water and the Environment, The University of Texas at Austin, Austin, Texas 78712, USA

Correspondence: Chan-Hoo Jeon (chanhoo.jeon@usm.edu)

Abstract

Avalanches and subaqueous debris flows are two cases of a wide range of natural hazards that have been previously modeled with non-Newtonian 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 non-Newtonian 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 non-Newtonian 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 multi-material level-set 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 (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 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.

1 Introduction

A wide range of natural hazards involve gravity-driven flows down a slope, for example, landslides (terrestrial or submarine), flood-driven 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 . 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 (Iverson1997). 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 large-scale 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 natural-hazard 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 fluid-solid structural interactions.

Large-scale natural-hazard flows have been widely investigated with field observations, small-scale 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 . This flow complexity is illustrated by the classification of subaqueous mass movements by 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 non-moving to a flowing regime can involve spatial heterogeneity and time-dependent behavior that is not well-represented by parameterizations of the flowing regime. Real-world 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 for further insight on granular flows and for heterogenous flows. The fundamentals physics of such flows is presented in . Herein, we do not seek to distinguish between the differing physics of these various complex flows but rather focus on advancing the use of non-Newtonian viscosity models as a proxy for their general behavior. For simplicity in exposition, we will use the term “debris flow” to refer to any real-world mixture modeled as a continuum fluid using a non-Newtonian model.

Following , 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 ; (ii) merging soil and fluid mechanics models; and (iii) representing the heterogeneous debris as a continuum fluid with behaviors similar to a non-Newtonian 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 non-Newtonian fluids but merely that some of their behaviors can be captured with an appropriately parameterized viscosity model . Indeed, 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 fluid-solid distributions and interactions. In using a non-Newtonian rheological model, the real-world interaction between solid particles and surrounding fluid in a heterogeneous mixture can be thought of as similar to the microstructural behavior of a homogeneous non-Newtonian fluid where the local fluid viscosity is a function of the local stress. The main advantage of this approach is that a non-Newtonian rheological model is simply a time/space-dependent viscosity term for the Navier–Stokes equations. It follows that the time/space-varying eddy viscosities in a wide range of existing hydrodynamic codes can be readily adapted to non-Newtonian behavior and used for parameterized modeling of debris flows.

Note that the terminology of non-Newtonian flows can be confusing as “time-independent” models have viscosities that can change with both space and time throughout a flow. The difference between a “time-independent” and a “time-dependent” non-Newtonian fluid is whether the relation between stress and viscosity (i.e., non-Newtonian equation itself) is allowed to change with time. Thixotropic (time-dependent) fluids are defined as non-Newtonian fluids where the process of “aging” during a flow changes the underlying fluid microstructure and the relationship between stress and viscosity . Herein, we examine how the use of a thixotropic model provides the ability to model behaviors that cannot be represented with a time-independent non-Newtonian model. Our goal is to provide insight into the research needs for further experiments and model development into the natural hazards of gravity-driven debris flows across the transitions from inception to stalling.

Gravity-driven 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 (Iverson1997). Parameterized non-Newtonian fluid models are an obvious approach to approximate these behaviors. Time-independent rheological models have been widely used to simulate debris flows ; however the real-world flow characteristics include time-dependent behaviors that could be categorized as “thixotropic” . Our focus in this paper is examining how a thixotropic model behavior compares to the more common time-independent (Herschel–Bulkley) non-Newtonian fluid model.

From a macroscale perspective, debris flows have similar behaviors to “yield-stress fluids” that have been studied as a class of non-Newtonian fluids . A yield-stress 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 time-independent non-Newtonian models cited above. 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 behaviors similar to a lower-viscosity fluid that deforms more easily under stress. This change from high viscosity to low viscosity under stress is readily simulated with a conventional time-independent non-Newtonian Herschel–Bulkley model. Arguably, what is missing from a time-independent 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 time-dependent element in the non-Newtonian model.

The simplest non-Newtonian yield-stress 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 time-independent. The approach proposed by is a common approach for representing the general case of time-independent non-Newtonian 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 $\stackrel{\mathrm{˙}}{\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} is shear thinning, n>1 is shear thickening, and n=1 corresponds to the Bingham plastic model (Bingham1916).

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 $\stackrel{\mathrm{˙}}{\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 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. 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 . A more successful approach was that of , 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 .

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 exist 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 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 steady-state conditions for debris flows should be reasonably represented with time-independent models. Indeed, 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 time-independent Herschel–Bulkley model at high shear rates. 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 . 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. 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 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 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 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 renewal of structure (consolidation) 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 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 . 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. 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 , although a few studies have numerically investigated the characteristics of thixotropic flow on a simple inclined plane . 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. Thixotropic flows modeled at the laboratory scale typically use clays (e.g., bentonite, kaolin) to create the microstructure controlling non-Newtonian behavior . Preparation of a homogenous clay suspension for such experiments is a demanding task, the details of which can be found in , , and . Unfortunately, we cannot expect the structure of a heterogeneous large-scale debris flow to mimic the flow scales, yield stresses, and parameters for a homogeneous thixotropic laboratory flow. However, lacking data from a large-scale debris flow that could be adequately used for model comparisons, herein we take a first step by analyzing how thixotropic models compare to time-independent models for laboratory-scale flows.

Validating the use of a non-Newtonian model to represent a real-world debris flow presents challenges on two levels: first, does the model correctly represent a non-Newtonian flow? Second, does the non-Newtonian flow (when parameterized) represent a real-world debris flow? To date, successful non-Newtonian models of real-world flows have been parameterized using a time-independent approach, which limits the ability of the model to represent the transition phases, i.e., flow initiation and stalling . Unfortunately, data on transition phases for real-world flows are lacking and are severely limited even for laboratory-scale flows.

In this paper we evaluate a time-independent Papanastasiou model and a time-dependent Coussot model for simulations of laboratory-scale avalanche and subaqueous debris flows, with comparisons to available experimental measurements. The governing equations are presented in Sect. 2, and the non-Newtonian Papanastasiou and Coussot viscosity models in Sect. 3. 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 Sect. 4. The numerical solver, using a multi-material level-set 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.

2 Governing equations

The governing equations in conservation form for unsteady and incompressible fluid flow can be written as

$\begin{array}{}\text{(2)}& & \mathrm{\nabla }\cdot \mathbit{u}=\mathrm{0},\text{(3)}& & \frac{\partial \mathbit{u}}{\partial t}+\mathrm{\nabla }\cdot \left(\mathbit{u}\otimes \mathbit{u}\right)=\frac{\mathrm{1}}{\mathit{\rho }}\left(-\mathrm{\nabla }p+\mathrm{\nabla }\cdot \mathbf{T}+\mathbf{f}\right),\end{array}$

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; uu is the dyadic product of the velocity vector u; and T is the viscous stress tensor:

$\begin{array}{}\text{(4)}& \mathbf{T}=\mathrm{2}\mathit{\eta }\mathbf{D},\end{array}$

where η denotes the plastic viscosity and D is the rate of strain (deformation) tensor:

$\begin{array}{}\text{(5)}& \mathbf{D}=\frac{\mathrm{1}}{\mathrm{2}}\left[\mathrm{\nabla }\mathbit{u}+\left(\mathrm{\nabla }\mathbit{u}{\right)}^{T}\right],\end{array}$

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.

The non-Newtonian 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/space-varying 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 finite-volume numerical discretization . For simplicity in the present work, we limit ourselves to a two-dimensional (2-D) 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 or the second-order TVD (total variation diminishing) Superbee scheme 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 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 . To evaluate the values at the cell surfaces, the Green–Gauss method is used and the momentum interpolation scheme is applied. The code is parallelized with MPI (Message Passing Interface), and PETSc (Portable, Extensible Toolkit for Scientific Computation) is used for standard solver functions (e.g., the stabilized version of the biconjugate gradient squared method with pre-conditioning by the block Jacobi method). The developed code has been verified by the method of manufactured solutions (further details provided in Jeon2015).

3 Non-Newtonian fluid models

The Herschel–Bulkley model, Eq. (1), was made more practical for modeling a fluid flow continuum by , 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 $\stackrel{\mathrm{˙}}{\mathit{\gamma }}=\mathrm{2}\sqrt{|{\mathbf{II}}_{D}|}$, where IID is the second invariant of the rate of strain as (Mei2007)

$\begin{array}{}\text{(7)}& {\mathbf{II}}_{D}=\frac{\mathrm{1}}{\mathrm{2}}\left[{\left(\text{tr}\left(\mathbf{D}\right)\right)}^{\mathrm{2}}-\text{tr}\left({\mathbf{D}}^{\mathrm{2}}\right)\right]={D}_{\mathrm{11}}{D}_{\mathrm{22}}-{D}_{\mathrm{12}}^{\mathrm{2}}\end{array}$

and Dij 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 (thixotropic) model of introduces dependency on a time-varying microstructure parameter (λ) in the general form

$\begin{array}{}\text{(8)}& \mathit{\eta }={\mathit{\eta }}_{\mathrm{0}}\left(\mathrm{1}+\mathit{\omega }{\mathit{\lambda }}^{n}\right),\end{array}$

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:

$\begin{array}{}\text{(9)}& \frac{\mathrm{d}\mathit{\lambda }}{\mathrm{d}t}=\frac{\mathrm{1}}{{T}_{\mathrm{0}}}-\mathit{\alpha }\stackrel{\mathrm{˙}}{\mathit{\gamma }}\mathit{\lambda },\end{array}$

where T0 is the characteristic time of the microstructure, α is a material-specific parameter, and $\stackrel{\mathrm{˙}}{\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 . That is, larger values of α require greater microstructure homogenization (smaller λ) to drive the system to steady-state conditions ($\mathrm{d}\mathit{\lambda }/\mathrm{d}t\to \mathrm{0}$).

4 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 (T0), and two material-specific 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 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 :

$\begin{array}{}\text{(10)}& {\mathit{\tau }}_{\mathrm{c}}=\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 well-defined 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 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 vs. strain ($\mathit{\tau }:\stackrel{\mathrm{˙}}{\mathit{\gamma }}$), 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 $\mathit{\tau }:\stackrel{\mathrm{˙}}{\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 T0, η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 debris-flow model (Sect. 8), 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 T0 of the Coussot model in Eq. (10) can also be troublesome to estimate. This characteristic time for aging, which described as “spontaneous evolution of the microstructure”, is not widely, used and the literature does not provide insight on how to evaluate T0 as a function of other rheological characteristics. Furthermore, T0 has slightly different definitions by authors of several papers. 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, defined T0 as “the characteristic time of build-up of the microstructure at rest”. Their characteristic time is close to the rest time of . Therefore, we make the assumption that the rest time measured in the Chanson's experiments is the same with the T0 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 T0, 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 time-dependent parameters from the time-independent parameters, which provides confidence that time-dependent and time-independent models are being compared in a reasonable manner.

Figure 1Concept of a graphical method B for estimating a consistent set of ${\mathit{\lambda }}_{\mathrm{0}},\mathit{\alpha },and\mathit{\omega }$ parameters for the Coussot model.

5 Multi-material level-set method

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 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 multi-phase fluids 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 .

The multi-material level-set method herein follows with the addition of high-order numerical schemes . The “level set” of the ith fluid is designated as ϕi:

where $i=\mathit{\left\{}\mathrm{1},\mathrm{2},\mathrm{\cdots },{N}_{\mathrm{m}}\mathit{\right\}}$, Nm 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:

$\begin{array}{}\text{(12)}& \mathit{\rho }\equiv \sum _{i=\mathrm{1}}^{{N}_{\mathrm{m}}}{\mathit{\rho }}_{i}{H}_{i},\phantom{\rule{1em}{0ex}}\mathit{\eta }\equiv \sum _{i=\mathrm{1}}^{{N}_{\mathrm{m}}}{\mathit{\eta }}_{i}{H}_{i},\end{array}$

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=di. 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 :

$\begin{array}{}\text{(14)}& \frac{\partial {\mathit{\varphi }}_{i}}{\partial t}+\mathbit{u}\cdot \mathrm{\nabla }{\mathit{\varphi }}_{i}=\mathrm{0}.\end{array}$

The above is coupled to a 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 ${\mathit{\varphi }}_{i}\left(\mathbit{x},t\right)=\mathrm{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 ϕ(tt), the values of the level set will no longer satisfy the eikonal condition of $|\mathrm{\nabla }{\mathit{\varphi }}_{i}|=\mathrm{1}$; 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 . This problem is addressed with “reinitialization”, which resets the ϕ(tt) 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 . Let $\stackrel{\mathrm{^}}{\mathit{\varphi }}$ be an estimate of the reinitialized value for ϕ(tt) in the equation

$\begin{array}{}\text{(15)}& \frac{\partial \stackrel{\mathrm{^}}{{\mathit{\varphi }}_{i}}}{\partial \mathcal{T}}+\mathbf{S}\left({\stackrel{\mathrm{^}}{\mathit{\varphi }}}_{i}\right)\left(|\mathrm{\nabla }{\stackrel{\mathrm{^}}{\mathit{\varphi }}}_{i}|-\mathrm{1}\right)=\mathrm{0},\end{array}$

where 𝒯 is the pseudo-time, and S is the signed function as

The time-advanced set of ϕ(tt) is the starting guess for $\stackrel{\mathrm{^}}{\mathit{\varphi }}$, and the steady-state solution of $\stackrel{\mathrm{^}}{\mathit{\varphi }}$ will satisfy $|\mathrm{\nabla }{\stackrel{\mathrm{^}}{\mathit{\varphi }}}_{i}|=\mathrm{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 . For the reinitialization step of Eq. (15), the second-order ENO (essentially non-oscillatory) scheme and a smoothing approach are used for the spatial discretization (further details are provided in Jeon2015).

6 Poiseuille flow of Bingham fluid

A two-dimensional Poiseuille flow in a channel driven by a steady pressure gradient of $\partial p/\partial 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

$\begin{array}{}\text{(17)}& u\left(y\right)=\left\{\begin{array}{l}\frac{\mathrm{1}}{\mathrm{2}\mathit{\eta }}\left(-\frac{\partial p}{\partial x}\right)\left({F}^{\mathrm{2}}-{y}^{\mathrm{2}}\right)-\left(\frac{{\mathit{\tau }}_{\mathrm{0}}}{\mathit{\eta }}\right)\left(F-y\right)\\ \phantom{\rule{1em}{0ex}}\text{for}\phantom{\rule{1em}{0ex}}{F}_{D}\le y\le F\\ \frac{\mathrm{1}}{\mathrm{2}\mathit{\eta }}\left(-\frac{\partial p}{\partial x}\right)\left({F}^{\mathrm{2}}-{F}_{D}^{\mathrm{2}}\right)-\left(\frac{{\mathit{\tau }}_{\mathrm{0}}}{\mathit{\eta }}\right)\left(F-{F}_{D}\right)\\ \phantom{\rule{1em}{0ex}}\text{for}\phantom{\rule{1em}{0ex}}\mathrm{0}\le y<{F}_{D}\end{array}\right\,\end{array}$

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 FD is a length scale representing the relationship between yield stress and the pressure gradient:

${F}_{D}=\frac{{\mathit{\tau }}_{\mathrm{0}}}{\left(-\frac{\partial p}{\partial x}\right)}.$

A convenient set of Bingham fluid parameters for the Poiseuille test cases can be extracted from the dip-coating study of 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 cross-stream direction. A Neumann boundary condition is applied along the lower boundary of the simulation domain, so the simulation includes only the upper half-channel of this symmetric flow.

Using the Papanastasiou model of Eq. (6) to approximate a Herschel–Bulkley model of a Bingham fluid requires time-scale parameter m to provide smooth behavior across the yield-stress 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.

Table 1Bingham fluid Herschel–Bulkley model parameters used in Poiseuille flow test cases, from .

7 Thixotropic avalanches

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. developed dam-break 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 θ, d0, and l0 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.

Figure 2Comparison of analytical and numerical solutions for steady-state fluid velocity for Poiseuille flow of a Bingham fluid.

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

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 ) 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 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 time-dependent Coussot model of Eqs. (8) and (9). For a time-independent 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 time-independent 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 . The time-dependent 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 time-dependent model uses the same value (n=1.1) as the time-independent 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.

Table 2Dimensions and data for thixotropic avalanche simulations corresponding with experiments by .

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

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 (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 time-dependent 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 fourth-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.

Figure 4Analytical stress–strain for thixotropic avalanche case 1: shear stress (Pa) and rate of strain (s−1) with τ0=31Pa and τc=90Pa. The τ axis is scaled for comparison with Figs. 5 and 6, while the $\stackrel{\mathrm{˙}}{\mathit{\gamma }}$ axis has an individual scale for clarity.

Figure 5Analytical stress–strain for thixotropic avalanche case 2: shear stress (Pa) and rate of strain (s−1) with τ0=21.1Pa and τc=165Pa. The τ axis is scaled for comparison with Figs. 4 and 6, while the $\stackrel{\mathrm{˙}}{\mathit{\gamma }}$ axis has an individual scale for clarity.

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 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 , repeated here as

$\begin{array}{}\text{(18)}& {x}_{s}^{\ast }=\frac{\mathrm{sin}\mathit{\theta }}{\mathrm{2}}{\left({t}^{\ast }\right)}^{\mathrm{2}}.\end{array}$

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).

Figure 6Analytical stress–strain for thixotropic avalanche case 3: shear stress (Pa) and rate of strain (s−1) with τ0=14Pa and τc=50Pa. The τ axis is scaled for comparison with Figs. 4 and 5, while the $\stackrel{\mathrm{˙}}{\mathit{\gamma }}$ axis has an individual scale for clarity.

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 $\mathit{\eta }=K{\stackrel{\mathrm{˙}}{\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 time-dependent models in case 1 to develop reasonable flow conditions despite the higher τ0. No doubt the time-independent 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.

Figure 7Thixotropic avalanche case 1: comparison of numerical results and experimental data for non-dimensional front displacement (x) as a function of non-dimensional simulation time (t).

Figure 8Thixotropic avalanche case 2: comparison of numerical results and experimental data for non-dimensional front displacement (x) as a function of non-dimensional simulation time (t).

Figure 9Thixotropic avalanche case 3: comparison of numerical results and experimental data for non-dimensional front displacement (x) as a function of non-dimensional simulation time (t).

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 2-D 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 2-D simulations have good agreement up until ${t}^{\ast }\sim \mathrm{3}$ and then show significant underprediction of the experiments. As noted by , 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 2-D 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 time-dependent 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 time-independent 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 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 three-dimensional controls on the front propagation in the experiment that cannot be represented in the present 2-D model.

8 Subaqueous debris flows

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 . We conducted simulations matching the 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 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 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.

Figure 10Submarine landslide.

Figure 11Subaqueous debris case 1: shear stress (Pa) and rate of strain (s−1).

Figure 12Subaqueous debris case 2: shear stress (Pa) and rate of strain (s−1). Axis scalings are identical to Fig. 11 for comparison purposes.

The parameters for the time-independent fluid model from  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 time-dependent 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 T0 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 time-independent and the time-dependent fluid models are shown in Fig. 11 for case 1 and Fig. 12 for case 2.

Table 4Dimensions for simulations to match experiments by .

Table 5Parameters of the time-independent fluid model.

Table 6Parameters of the time-dependent fluid model.

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. Figures 15 and 16 show the evolution of the run-out distance (L) for case 1 and case 2, respectively. It can be seen that both time-independent and time-dependent models are reasonable approximations of the limited experimental data. Both types of models appear to underestimate the initial run-out and slightly overestimate later times.

Figure 14Profiles of debris and water (ϕ2: water).

Figure 15Subaqueous debris case 1: run-out distance (L, m) as a function of simulation time (t, s).

Figure 16Subaqueous debris case 2: run-out distance (L, m) as a function of simulation time (t, s).

Figure 17Subaqueous debris case 1: height of head flow (H, m) and water depth (D, m) as a function of simulation time (t, s).

Figure 18Subaqueous debris case 2: height of head flow (H, m) and water depth (D, m) as a function of simulation time (t, s).

Figure 19Subaqueous debris case 1: flow-front velocity (U, m s−1) as a function of simulation time (t, s).

Figure 20Subaqueous debris case 2: flow-front velocity (U, m s−1) as a function of simulation time (t, s).

Figures 17 and 18 show a comparison of H and D for simulations and experiments. Again, within the limited availability of experimental data, both time-independent and time-dependent models provide reasonable agreement. Figures 19 and 20 show similar agreement for the front velocities, although the experimental data are insufficient to validate the wave-like oscillation of the velocity in the simulations.

These results indicate the multi-material level-set model is capable of representing the key features in a subaqueous debris flow. For this flow, the use of the simpler time-independent 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

9 Discussion and conclusions

This work shows that a multi-phase flow solver using a multi-material level-set method with yield-stress models of non-Newtonian viscosity provides a means for numerical approximation of avalanches and subaqueous debris flows. This simulation approach was tested with both time-independent (Herschel–Bulkley, Papanastasiou, Bingham plastic) and time-dependent (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 time-independent 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 time-dependent approach provides the hysteresis that is desirable in a debris-flow model, in comparisons with experimental data the time-dependent Coussot model provides a clear advantage for only for a single case – where the Herschel–Bulkley and Bingham plastic models erroneously predicted near-zero flow. Nevertheless, much of the complexity in real-world 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 time-dependent 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 time-independent viscosity–stress relationships that are often used for non-Newtonian flow models of natural hazards 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 the structure is recovered, are both poorly understood. For time-dependent thixotropic models to be useful in modeling real-world avalanches and debris flows, there is a need for a consistent approach to defining the initial jamming (λ0), the characteristic time of aging (T0), and the asymptotic shear viscosity (η0), along with the material parameters ω and α for real-world systems. As yet, these parameters are not well defined for either simple laboratory models or complex real-world 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 real-world 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 time-dependent (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 real-world natural hazards can be modeled in this manner.

Data availability
Data availability.

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 (Chan-Hoo Jeon).

Competing interests
Competing interests.

The authors declare that they have no conflict of interest.

Acknowledgements
Acknowledgements.

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: Jean-Philippe Malet
Reviewed by: two anonymous referees

References

Ancey, C.: Plasticity and geophysical flows: A review, J. Non-Newton. Fluid, 142, 4–35, 2007. a

Ancey, C. and Cochard, S.: The dam-break problem for Herschel-Bulkley viscoplastic fluids down steep flumes, J. Non-Newton. Fluid, 158, 18–35, 2009. a

Aziz, M., Towhata, I., Yamada, S., Qureshi, M. U., and Kawano, K.: Water-induced 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/nhess-10-1229-2010, 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. Non-Newton. Fluid, 142, 219–243, 2007. a

Beverly, C. R. and Tanner, R. I.: Numerical analysis of three-dimensional Bingham plastic flow, J. Non-Newton. 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 debris-flow material in South-Gargano watersheds, Nat. Hazards, 54, 209–223, 2010. a

Bovet, E., Chiaia, B., and Preziosi, L.: A new model for snow avalanche dynamics based on non-Newtonian 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, Earth-Sci. 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 slip-debris flow initiation processes in pyroclastic deposits: the Sarno 1998 event, Nat. Hazards Earth Syst. Sci., 3, 53–69, https://doi.org/10.5194/nhess-3-53-2003, 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 macro-viscous 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 particle-laden flows: from laboratory experiments to field observations, J. Phys. D-Appl. Phys., 50, 053001, https://doi.org/10.1088/1361-6463/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 flow-front 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 Gummi-Benzollö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 debris-flow rheology myth, in: The 3rd International Conference on Debris-Flow 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 three-dimension terrain: 1. Coulomb mixture theory, J. Geophys. Res., 106, 537–552, 2001. a

Jeon, C.-H.: Modeling of Debris Flows and Induced Phenomena with Non-Newtonian 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 Fine-Grained 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 start-up 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 (Alpes-de-Haute-Provence, France), in: The 3rd International Conference on Debris-Flow 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/nhess-12-3377-2012, 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 PDE-based 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 fine-grained mud from Eastern Canada, Eng. Geol., 43, 31–44, 1996. a

Pierson, T. and Costa, J.: A rheologic classification of subaerial sediment-water 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 calibration-based approach for debris-flow forward-analyses, Nat. Hazards Earth Syst. Sci., 10, 1009–1019, https://doi.org/10.5194/nhess-10-1009-2010, 2010. a, b

Pudasaini, S. P.: A general two-phase 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/nhess-10-2507-2010, 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 non-oscillatory shock capturing schemes II, J. Comput. Phys., 83, 32–78, 1989.  a, b

Smith, G. A. and Lowe, D. R.: Lahars: volcano-hydrologic 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, interface-preserving 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 two-phase 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 two-phase 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/nhess-11-3053-2011, 2011. a, b