A Quasi-Single-Phase Model for Debris Flows Incorporating Non-Newtonian Fluid Behavior

: Debris-ﬂow modeling is a great challenge due to its complex physical mechanism that remains poorly understood. The present research incorporates the effect of rheological features of the non-Newtonian ﬂuid into the complete quasi-single-phase mixture model, which explicitly accommodates the interactions between ﬂow, non-uniform sediment transport, and bed evolution. The effect of rheological features is estimated by Hersch–Bulkley–Papanastasiou model that can be simpliﬁed to Bingham or Newtonian models with speciﬁc coefﬁcients. The governing equations are solved by a fully conservative numerical algorithm, using an explicit ﬁnite volume discretization with well-balanced slope-limited centered scheme combined with an implicit discretization method. One set of large-scaled U.S. Geological Survey debris-ﬂow experiments is applied to investigate the inﬂuence of the non-Newtonian ﬂuid on debris-ﬂow modeling. It is found that the present model closed by Hersch–Bulkley–Papanastasiou model performs better than the models without considering effect of rheological features, which may facilitate the development of mixture modeling for debris ﬂows. x = 32 and 66 m. It is noted that thickest ﬂow depth lasts about 2.5 s for measurement while only 0.5 s for computation results at x = 66 m. This may be due to the effect of interphase (water-to-sediment) and particle–particle interactions of debris ﬂow cannot be fully expressed by the mixture model. The ﬂow front arrival points and the steady state after about t = 10 s are closer to the measured data.


Introduction
Debris flow, composed of highly concentrated mixtures of sediments and water [1], is extremely destructive and represents an enormous threat to people's life and property, public infrastructure, and ecological environment throughout their trajectories [2]. It is regarded as a physical process between landslide and flood flow. Hence, besides the characteristics of multiphase, fragmentation, flow, and yield, debris flow is also multiscale. The shape of the solid particles that constitute debris flow is extremely irregular, and the scale spans several orders of magnitude, ranging from clay, silt, fine sand, and gravel to boulders with diameter of several meters [3]. The dynamic process of debris flow presents a series of complex physical phenomena, composed of triggering, movement, and accumulation. It is difficult to observe the whole process in the field, and therefore, numerical simulation are always carried out to study the evolution of debris flow.
The numerical dynamic model of debris flow can be divided into continuous medium model, discrete medium model, and mixed medium model [4]. The continuous medium model is generally established by fluid researchers, while the discrete medium model is usually used by geotechnical researchers. Among these two models is the mixed medium model, where the continuity of fluid and the dispersion of particles are considered. The continuous medium can be divided into quasi-single-phase model [5][6][7] and multiphase model [8][9][10][11][12]. The former is characterized by a single momentum equation for watersediment mixture flow based on the assumption of the same velocity of fluid and solid phases, while the latter features the separate momentum equation for each phase, considering the velocity differences among phases. The two-phase model is promising for its more complete physical process to reveal the relative motions and interactions between the fluid and solid phases [11,12]. However, the high computing costs and the demand for extra closure relationships make it hard to be applied to large-scale field. In contrast, the quasisingle-phase model is more practical and has seen wide application for its much higher computation speed. However, the existing quasi-single-phase models used in engineering are often simplified. For example, some models are decoupled [13,14]; some are based on the sediment capacity assumption [15,16]; and some are restricted to uniform sediment transport cases and generally ignore the effects of debris flows fluctuations [5,17], which may not be generally justified from physical perspectives.
Meanwhile, rarely do models involve the rheological features of debris flows into a relatively complete single-phase model. However, it is known that debris flows in field cases often involve poorly sorted, heterogeneous particles in a non-Newtonian way rather than the ideal Newtonian liquid [18]. Hence, a rheological relation is warranted for numerical simulation to describe the non-Newtonian behavior of debris flow. The last several decades have witnessed the great efforts made by researchers to conduct their investigations on this field, for example, the Bingham model [19,20], Cross model [21,22], Herschel-Bulkley (HB) model [23], Herschel-Bulkley-Papanastasiou (HBP) model [24,25] and other non-Newtonian models. Among them, the Bingham model is widely used in debris-flow simulation, but the effective viscosity in the Bingham model will be infinite when the shearing rate approaches zero, leading to the divergence in the numerical simulation. Meanwhile, the Bingham model or the general Cross model in debris-flow simulations is the linear constitutive law between the shear stress tensor and the strain rate tensor under high strain rate [25]. Since many studies [26][27][28] have observed nonlinear behavior, particularly the shear thinning and shear thickening phenomenon in debris flow, the HB model is proposed and considered suitable to describe the nonlinear features of debrisflow [29]. However, the problem of numerical divergence still remains in the HB model, and Papanastasiou (1987) [24] introduced an improved version of the HB model (i.e., HBP model) to solve this problem. Here, a quasi-single-phase mixture model combined with HBP model is proposed to reveal the effect of rheological features of the non-Newtonian fluid for debris flows. The model is solved by a fully conservative numerical algorithm and then applied to the simulation of a U.S. Geological Survey debris-flow experiment.

Mathematical Model
The effect of rheological features of non-Newtonian fluid is incorporated into the present debris-flow model, which is an extension of the quasi-single-phase model by Xia et al. (2018) [7]. For one-dimensional debris flows over an inclined erodible bed composed of sediment with N size classes (d k denotes the diameter of the kth size sediment, and k = 1, 2, . . . , N), the complete depth-averaged, one-dimensional governing equations can be written as follows: where t = time; x = streamwise coordinate parallel to bed slope with the angle θ; h = debrisflow depth; U = depth-averaged mixture flow velocity in x direction; z b = bed elevation; c k = depth-averaged size-specific volumetric sediment concentration and C = ∑ c k = total volumetric sediment concentration; g = gravitational acceleration and g z = g cos θ, g x = g sin θ; τ b = bed shear stress; p = bed sediment porosity; ρ f , ρ s = densities of fluid and solid phases, respectively; ρ m = density of fluid-solid mixture; ρ 0 = density of the bed; E k = size-specific sediment entrainment flux and E T = ∑ E k = total sediment entrainment flux; D k = size-specific sediment deposition flux and D T = ∑ D k = total sediment deposition flux; δ = thickness of the active layer [30], and δ = 2d 84 , where d 84 is the particle size at which 84% of the sediment are finer; f ak = fraction of the kth size sediment in the active layer; ξ = z − δ = elevation of the bottom surface of the active layer; f Ik = fraction of the kth size sediment in the interface between the active layer and substrate layer; T R = depth-averaged stress due to fluctuations of debris flows in the x direction [12]; and T µ = depth-averaged viscous stress of debris flows. Notably, for ease of description, Equation (2) is written in a format similar to Xia et al. (2018) [7], and τ add = − τ e f f − τ N is the additional shear stress compared to Newtonian governing equations [7], where ; µ e f f = effective viscosity of debris flows; µ N = dynamic viscosity for Newtonian fluid; . γ = ∂U/∂z = shear rate; and z = the coordinate perpendicular to bed.
Based on the shallow water hydro-sediment-morphodynamic theory [31], the proposed model is fully coupled, explicitly incorporating the mass exchange between the flow and the bed. It equips the model with wider applications to both fixed and erodible bed. Meanwhile, a non-capacity approach is used, which determines the size-specific sediment transport by accommodating the contributions of advection due to mean flow velocity and of the mass exchange with the bed. It is noted that rheological models to reflect the non-Newtonian flow features are involved in equation derivations, which is shown by the effective viscosity µ e f f [22].
To close the governing Equations (1)-(5), a set of relationships must be introduced to determine the bed shear stresses, sediment entrainment and deposition, stresses due to fluctuations, shear rate, as well as the effective viscosity due to different rheological models. According to Xia et al. (2018) [7], the bed shear stresses for the debris flow are comprised of the bed shear stresses exerted on the fluid phase estimated by Manning's equation and solid phase evaluated by the revised Lucas formula [32]. Sediment exchange with the bed is one of the key components of computational models for debris flow over erodible beds. However, most existing bed entrainment rate equations fail to consider the effects of particle size, which may be questionable from a physical perspective that fine grains are easier to be eroded than large blocks [33]. Therefore, the closure model [34][35][36] is adopted to calculate the sediment entrainment and deposition according to Li et al. (2018b) [12]. Inspired by the study of Cao et al. (2015) [37] on roll waves over steep slopes, the stresses due to fluctuations can be analogous to turbulent flows, which are determined by the depth-averaged k − ε turbulence model [38] and a modification component [37,39]. Based on the assumption of linear velocity distribution along the depth, shear rate is estimated approximately by The key point to reveal the behavior of a non-Newtonian fluid is to evaluate the effective viscosity varying with shear rate [22]. The stress-strain relationships during plastic yield process can be estimated by the Hersch-Bulkley-Papanastasiou model [24], which is an enhanced version of the Hersch-Bulkley model that can express the shearthinning and shear-thickening phenomena. This HBP model is used and discussed in this study, which reads as where m, n = constant coefficients; µ = dynamic viscosity; and τ y = yield stress. They are estimated as follows [42]:  4 ; and σ = (C − C v0 )/C vm . It can be seen from Equation (6) that if m → ∞ , the HBP model can be reduced to Hersch-Bulkley model. n < 1 and n > 1 represent the shear-thinning and shear-thickening phenomena, respectively. When n = 1, m → ∞ , the HBP model can be reduced to Bingham model, and when n = 1, m = 0, the HBP model is simplified to Newtonian model.
The governing equations are solved by a fully conservative numerical algorithm using well-balanced slope-limited centered scheme [43,44] combined with double-sweep method [7,37].

Model Comparison
The large-scale U.S. Geological Survey debris-flow experiments provided a systematic set of observed data that are well-suited for validating mathematical models [45][46][47]. In order to evaluate the influences of the rheological relation, the subset of the fixed bed experiment equipped with sand-gravel-mud sediment over a rough bed [46] is investigated by the model of Xia et al. (2018) [7] and the present model.

Case Description
The 95 m long, 2 m wide, and 1.2 m deep flume with a uniform bed slope of 31 degrees throughout most of its length was applied, the schematic flume geometry is shown in Figure 1. A 2 m high vertical headgate used to retain static debris prior to its release was set to be the reference of longitudinal distances, i.e., x = 0 m. At x = 74 m, the flume bed begins to flatten [46]. The flow-front positions were tracked, and the flow thickness at three cross-sections of the uniformly sloping flume located at x = 2 m, 32 m, and 66 m were also measured. Initially, approximately 10 m 3 of sand-gravel-mud mixture debris were released by the sudden opening of the two-piece steel headgate. They moved across the smooth bed from x = 0 to 6 m with the basal friction angle of 28 • and then the rough bed after x = 6 m with the basal friction angle of 40 • .    The present study mainly uses a one-dimensional model to simulate the debris-flow evolution; therefore, the computational domain consists of the uniformly sloping flume and the adjacent runout pad that is considered to have the same width as the flume.

Results
The numerical simulations are performed within the time period before the forward and backward waves reach the downstream and upstream boundaries, respectively; thus, the boundary conditions can be simply set according to the initial status.
The detailed grain-size distribution of the debris can be obtained according to Iverson et al. (2010) [46]. The coefficients in closure models of bed shear stresses, sediment entrainment and deposition, and stresses due to fluctuations are all the same as Xia et al. (2018) [7]. For effective viscosity, the dynamic viscosity for Newtonian fluid is set to be µ N = 0.001 Pa · s. First, we used the simplified HBP model with n = 1, m → ∞ (i.e., Bingham model) to evaluate the effect of rheological relation.  [7]. Just downslope from the flume headgate at x = 2 m, the computed flow depths are almost the same ( Figure 2B), so the rheological feature may not show its impact at this time. As the propagation of the debris flow, the stress-strain relationship due to the interaction forces of internal debris flow, described by the HBP model, contributes to a better agreement with the measured flow thickness as compared to Xia et al. (2018) [7] at cross-sections x = 32 and 66 m. It is noted that thickest flow depth lasts about 2.5 s for measurement while only 0.5 s for computation results at x = 66 m. This may be due to the effect of interphase (water-to-sediment) and particle-particle interactions of debris flow cannot be fully expressed by the mixture model. The flow front arrival points and the steady state after about t = 10 s are closer to the measured data.

Sensitivity Analysis
Flow front and shear stress are two important factors to judge whether a num model can successfully simulate the main characteristics of debris flow or not. Hen interesting to investigate the sensitivity of the simulation results by the present m

Sensitivity Analysis
Flow front and shear stress are two important factors to judge whether a numerical model can successfully simulate the main characteristics of debris flow or not. Hence, it is interesting to investigate the sensitivity of the simulation results by the present model to coefficients n and m. First, m is set to be infinity, and n is tuned by 30% to check its effect on the results. Then, n is set to be 1, and the results of the track of flow front are computed by m = 0, 0.01 and infinity, respectively. Figure 3A compares the computed flow front positions with the measured data. It seems that the results by the model with different n are almost overlapped despite the improvement compared to Xia et al. (2018) [7]. According to Equation (6), the ratio of shear stress τ e f f to yield stress τ y can be written as: γ (9) ater 2022, 14, x FOR PEER REVIEW with the results shown in Figure 3A.    n reflects the shear-thinning (n < 1) and shear-thickening (n > 1) behavior of a non-Newtonian fluid, as shown in Figure 3B. With the variation of shear rate . γ, the ratio of shear stress τ e f f to yield stress τ y shows a linear relationship when n = 1, while the ratio increases if n > 1 and decreases if n < 1 when shear rate . γ is over 1. Furthermore, from Figure 3B, the coefficient n mainly controls the linear or nonlinear behavior in the high shearing rate range, while the most values of shear rate in this study case are below 300 with limit differences among shear stress S τ e f f , which indicates the little effects on the flow front tracking ( Figure 3A) due to tuned n. When n = 1, m = 0 and infinity form the result interval of flow front position, and m = 0.01 is in between, shown in Figure 3C. As the m gets smaller, the computed results are closer to that of Newtonian model. Notably, when n = 1 and m = 0, stress-strain relationship is reduced to Newtonian model (red line in Figure 3D), and when n = 1 and m → ∞ , stress-strain relationship is reduced to Bingham model Figure 3D), and the coefficient m mainly controls the initial rapid growth of shearing stress.
The L 1 norm of flow front location is deployed to quantify the impacts of m and n on debris-flow modeling, as shown in Table 1. The L 1 norm is defined as , where x f = computed flow front location; x * f = measured flow front location. The L 1 norms of tuned n are almost the same, which is consistent with the results shown in Figure 3A. The L 1 value of flow front simulation with Newtonian model (n = 1, m = 0) is 2.56. It is bigger than the results using the non-Newtonian model, which are 2.03 and 1.702 for m = 0.01 and m → ∞ , respectively. Generally, the coefficient n has limited influence on flow front simulation, maybe due to the scale of the experiment and the composition of the sediment. However, considering the non-Newtonian model can improve the accuracy of the model in tracking the flow front.

Effects of Additional Shear Stress Due to Non-Newtonian Fluid
Compared to the model of Xia et al. (2018) [7], based on Newtonian fluid assumption, the present model includes the additional shear stress τ add to describe the rheological features of debris flows, which actually is the difference between viscosity of non-Newtonian fluid (i.e., µ e f f ) and Newtonian fluid (i.e., µ N ). To provide insight into the contributions of the rheological features of the non-Newtonian fluid to debris-flow modeling shown by the additional shear stress, the spatial distribution ratios of S τ e f f /S G (S τ e f f = τ add ; S G = ρ m g z h) at specific instants along with the flow thicknesses are present, as shown in Figure 4. It illustrates that the ratio is negative around the flow front. Although the value is little from the trough to the peak of the debris flow, it is considerable around the debris-flow front, and the additional shear stress S τ e f f is by no means negligible compared with the gravitational term S G . This effect is similar to that of stress terms due to fluctuations in Xia et al. (2018) [7], and their effects are almost at the same level as the ratios and have the same order of magnitude ( Figure 9 in Xia et al. (2018) [7] compared to Figure 4 in present paper).
in Xia et al. (2018) [7], and their effects are almost at the same level as the ratios a the same order of magnitude (Figure 9 in Xia et al. (2018) [7] compared to Fig  present paper).

Conclusions
Under the framework of shallow water hydro-sediment-morphodynamics plete numerical model for debris flow is proposed. Unlike some previous simplifi els, the present model not only incorporates the effect of rheological features of t Newtonian fluid but also explicitly accommodates the interactions between flo uniform sediment transport, and bed evolution. This promotes the accuracy of flow modeling from a theoretical perspective. A fully conservative numerical al is used to solve the governing equations. The hyperbolic system that can be readil by many existing schemes is disrupted due to the fluctuation stress. The explic volume discretization combined with the implicit discretization method is ap solve this mathematical problem by efficient in-house program code. It lays the tion for application to both laboratory and field-scale debris-flow modeling.
One set of large-scaled U.S. Geological Survey debris-flow experiments was to investigate the non-Newtonian fluid on debris-flow modeling. Hersch-B Papanastasiou model is used to estimate the rheological feature of debris flow be is more generalized and considered suitable for describing the nonlinear features o flow. It is found that the present model closed by HBP model performs better t existing quasi-single-phase mixture model. The coefficients n and m in the HBP m discussed. The results illustrate that the coefficient n plays a role in the high shear range, while the coefficient m has effects on the initial growth of shearing stress. the additional shear stress due to non-Newtonian fluid is also investigated, and i be negligible as compared with the gravitational term.
Further studies are needed to be carried out because the rheological features varied due to the composition of sediments, soil moisture, and even the scale o flows. More suitable stress-strain relationships are desiderated to be found. Be dimensional model can be extended in the near future.

Conclusions
Under the framework of shallow water hydro-sediment-morphodynamics, a complete numerical model for debris flow is proposed. Unlike some previous simplified models, the present model not only incorporates the effect of rheological features of the non-Newtonian fluid but also explicitly accommodates the interactions between flow, non-uniform sediment transport, and bed evolution. This promotes the accuracy of debris-flow modeling from a theoretical perspective. A fully conservative numerical algorithm is used to solve the governing equations. The hyperbolic system that can be readily solved by many existing schemes is disrupted due to the fluctuation stress. The explicit finite volume discretization combined with the implicit discretization method is applied to solve this mathematical problem by efficient in-house program code. It lays the foundation for application to both laboratory and field-scale debris-flow modeling.
One set of large-scaled U.S. Geological Survey debris-flow experiments was applied to investigate the non-Newtonian fluid on debris-flow modeling. Hersch-Bulkley-Papanastasiou model is used to estimate the rheological feature of debris flow because it is more generalized and considered suitable for describing the nonlinear features of debris flow. It is found that the present model closed by HBP model performs better than the existing quasi-single-phase mixture model. The coefficients n and m in the HBP model are discussed. The results illustrate that the coefficient n plays a role in the high shearing rate range, while the coefficient m has effects on the initial growth of shearing stress. Besides, the additional shear stress due to non-Newtonian fluid is also investigated, and it cannot be negligible as compared with the gravitational term.
Further studies are needed to be carried out because the rheological features may be varied due to the composition of sediments, soil moisture, and even the scale of debris flows. More suitable stress-strain relationships are desiderated to be found. Besides, 2-dimensional model can be extended in the near future.
Author Contributions: Conceptualization, C.X.; methodology, C.X.; validation, H.T.; writingoriginal draft preparation, C.X. and H.T.; writing-review and editing, C.X. and H.T. All authors have read and agreed to the published version of the manuscript.
Funding: This research was funded by the National Natural Science Foundation of China, grant number 12002310.