Numerical Simulation of Non-Homogeneous Viscous Debris-Flows Based on the Smoothed Particle Hydrodynamics (SPH) Method

: Non-homogeneous viscous debris ﬂows are characterized by high density, impact force and destructiveness, and the complexity of the materials they are made of. This has always made these ﬂows challenging to simulate numerically, and to reproduce experimentally debris ﬂow processes. In this study, the formation-movement process of non-homogeneous debris ﬂow under three di ﬀ erent soil conﬁgurations was simulated numerically by modifying the formulation of collision, friction, and yield stresses for the existing Smoothed Particle Hydrodynamics (SPH) method. The results obtained by applying this modiﬁcation to the SPH model clearly demonstrated that the conﬁguration where ﬁne and coarse particles are fully mixed, with no speciﬁc layering, produces more ﬂuctuations and instability of the debris ﬂow. The kinetic and potential energies of the ﬂuctuating particles calculated for each scenario have been shown to be a ﬀ ected by the water content by focusing on small local areas. Therefore, this study provides a better understanding and new insights regarding intermittent debris ﬂows, and explains the impact of the water content on their formation and movement processes.


Introduction
As a frequent natural disaster, debris flow poses a serious threat to the lives and properties of people living in mountain areas [1][2][3][4]. This natural phenomenon needs to be taken into consideration when planning new developments in mountain catchments because its sudden, devastating, and extensive impacts could have strong consequences on the local economy. Debris flows can be divided into multiple categories such as viscous flow, dilatant flow, dilute flow, water-rock flow, etc. [5] based on the variety of materials and their combinations. The first category mentioned, the viscous debris flow, is particularly widely distributed in the mountainous areas of the southwest of China. This particular debris flow is characterized by high viscosity, wide particle-size distributions and uneven velocity distributions [6]. The difficulty in describing its motion process is due to the fact that the viscous debris flow is classified as a heterogeneous, non-constant, and non-Newtonian flow [7]. In previous studies, Johnson [8] attempted to describe the motion of viscous debris flow by using the Bingham Model, which is applicable to laminar viscous flows because in such fluids, the concentration of coarse particles is quite low. Bagnold [9] described it in terms of dilatant flow, and emphasized the discrete force between particles caused by collision. Based on the experiment conducted by Bagnold, Takahashi [10] proposed the idea of the presence of a vertical collision stress between particles, which supports the coarse particles not sinking into the fluid. Chen [11,12] comprehensively focused on the elastic shear force, the plastic shear force, and the laminar shear force in muddy fluids, implementing a modification of the motion resistance which could be suitable for viscoplastic fluids containing certain coarse particles. O'Brien [13,14] took into account the effects of plasticity, viscosity, collision, and turbulence on resistance forces and identified a general equation which includes the components due to cohesive yield stress, Mohr-Coulomb shear force, viscous shear stress, turbulent shear stress, and dispersive shear force, establishing a debris fluid model that could fully describe the effects due to different particle compositions and different densities.
Furthermore, the mesh-free numerical modeling techniques are being gradually introduced into this field. Dai [15] developed a three-dimensional model to simulate rapid landslide motions across 3D terrains. The artificial viscosities linearly related to the linear and quadratic terms of shear deformation were incorporated into the pressure terms in the momentum equation to dissipate energy for avoiding numerical oscillations. However, Dai [15] did not consider the effect of yield stress or the interaction between solid and liquid phases. Hosseini [16] adopted an innovative treatment similar to the one applied in Eulerian formulations to viscous terms, and hence facilitated the implementation of various inelastic non-Newtonian models. This approach utilized the exact forms of the shear strain rate tensor and its second principal invariant to calculate the shear stress tensor. Rodriguez-Paz [17] introduced a new frictional approach for the boundary conditions and modified constitutive equations in the SPH (Smoothed Particle Hydrodynamics) method to focus on the interaction between fluid particles and boundary conditions. The resulting technique was then applied for the numerical simulation of debris flows and the results were compared with those experimentally obtained and available in literature [18,19], providing good agreements.
In this paper, the SPH method was applied to estimate the movement of a non-homogeneous viscous debris-flow: the fluid under investigation was divided into solid-liquid phases. The solid phase was characterized by particles with larger size than the critical particle size, while the liquid phase was the mixture of water and particles with smaller size than the critical one. The constitutive relation for the liquid phase was characterized by yield stress, laminar viscous force, and turbulent viscous force, while the constitutive relation of the solid phase was related to support force, friction, and collision stresses. It is well known that the magnitude of shear deformation determines which force plays a dominant role in the process of the fluid movement [20]. Therefore, quantifying all the above forces, it could be possible to quantify and estimate the role of each component and further investigate how the shear sharp deformation could be reduced. However, due to the complexity of the materials composing debris flows, even the fluid with same rheological coefficients could generate effects attributable to different flow structures and characteristics. Therefore, the effect on the debris flow process of the initial vertical distribution of the two-phase solid-liquid is also considered in this study.

The SPH Method
SPH is a kind of mesh-free method based on a pure Lagrangian description, which has been widely applied in multiple engineering and science fields [21][22][23][24]. Compared with the mesh method based on continuum theory, the SPH method avoids the problem of mesh distortion in dealing with the flow issue since there is no connectivity between the particles, since it is developed on a uniform smoothed particle hydrodynamic framework. By adopting this technique, the goal is to provide accurate and stable numerical solutions for integral equations or partial differential equations (PDEs) using a series of arbitrarily distributed particles carrying field variables, such as mass, density, energy, and stress tensors [25]. There are two main steps to transform PDEs equations into the SPH form, called particle approximation and kernel approximation, respectively [26]. The first step consists of representing a function in continuous form as an integral representation by using an interpolation function. In this step, the approximation of the function and its derivatives is based on the evaluation of the smooth kernel function and its derivatives. The second step involves representing the problem domain by using a set of discrete particles within the influence area of the particle at location x, and then estimating the field variables for those particles as follows: where x' denotes the position of continuum in the influence domain before the discretization; W denotes the smoothing function; h is a parameter that defines the size of the kernel support, known as the smoothing length; Ω represents the problem space whose radius is taken as several times of h according to different smoothing functions; N is the total number of neighboring particles; m is the mass; and ρ is the density. Kernel approximation is the technique of approximating the values of both the field function and the derivative of the field function. The kernels used in the SPH method approximate a δ function (the Dirac function). Monaghan [27] suggested that a suitable Kernel approximation must have a compact support in order to ensure zero interactions outside its computational range. The original calculations of Gingold and Monaghan [28] used a Gaussian Kernel. The Gaussian Kernel function is simple to use and has high accuracy. Especially for the case of disordered particle distribution, this technique generates stable and accurate approximation results. However, the Gaussian Kernel function does not have a strict compact support unless the size of the Kernel support approaches the infinity value. Additionally, further various Kernels forms with a compact support (such as spline [29], super-Gaussian [30], polynomial [31], and cosine [32]) were proposed in previous studies but the one of the most popular Kernels more commonly utilized is based on the spline functions [29] as defined by: where q = |r|/h and r is the separation distance between the particles. In this study, Equation (3) has been used to approximate the values of the field under investigation. This Kernel has compact support so that its interactions are exactly zero for r > 2h. The smoothing distance or so called "Kernel range" h determines the degree with which a particle interacts with adjacent particles.

Gradient and Divergence
As a standard procedure, the gradient and divergence operators need to be formulated in a SPH algorithm if the simulation of the Navier-Stokes equations is to be attempted. In this work, the following commonly used forms are employed for gradient of a scalar A and divergence of a vector A [33]: where ∇ a W ab is the gradient of the Kernel function W x − x j , h with respect to x j , subscripts a and b represent the target particles and the particles in the influence domain, respectively, affecting the position of particle i. This choice of discretization operators ensures that an exact projection algorithm is produced. To date, there are various options to represent these operators, but only certain specific ones [34,35] have proven to be more convenient in terms of the accuracy and robustness of the method.

Governing Equations
The governing equations for transient compressible fluid flow include the conservation of mass and momentum equations. In a Lagrangian framework, these can be written as follows: Dv where t is time, v is the particle velocity vector, g is the gravitational acceleration, P stands for pressure, and D/Dt refers to the material derivative. The density ρ has been intentionally kept in the equations to be able to enforce the incompressibility of the fluid. Using an appropriate constitutive equation to model the shear stress tensor τ, Equations (6) and (7) can be used to solve both Newtonian and non-Newtonian flows.
The momentum equations include three driving force terms, i.e., body force, forces due to divergence of the stress tensor, and the pressure gradient, and these always have to be considered together with the incompressibility constraint. In a SPH formulation, the above system of governing equations must be solved for each particle at each time-step, and the order with which the force terms are incorporated into the momentum equations can be different from one algorithm to another.

Equation of State
SPH method can be formulated for incompressible and compressible flows. The equation of state, giving the relationship between particle density and fluid pressure, can be written as follows [28]: where P 0 represents a constant value of pressure, usually expressed in terms of initial pressure and ρ 0 is the reference density.

Viscous Terms
In the context of the SPH method, several forms of viscosity terms were introduced by Lucy [25], Gingold and Monaghan [29], Wood [34], Loewenstein and Mathews [36], and Shao and Lo [37]. As the purpose of this work was to solve non-Newtonian fluids, a new description of viscosity was developed to facilitate the modeling of such flow characteristics. Viscosity of incompressible Newtonian fluids depends only on the second principal invariant of the shear strain rate: Water 2019, 11, 2314 5 of 16 In solving Equations (8) and (9), the finite difference method should firstly be used to solve the total derivative between two particles, and then decompose the results into x, y directions, following Lo and Shao [38].
The viscous debris flow studied in this paper is a non-Newtonian fluid whose constitutive equation has the following form: Bingham fluid : viscoplastic fluid : τ = 1 2 where τ B refers to the yield stress; µ B , µ 1 are the coefficients of friction and µ D , µ 2 are the coefficients of collision stresses. According to the different flows considered, the symbols in Equation (13) can represent different meanings. When the object considered corresponds to solid particles, µ 0 indicates the static support force between the solid particles, µ 1 is the particle friction coefficient, and µ 2 is the particle collision coefficient. When considering a liquid phase slurry, µ 0 is the yield stress, µ 1 is the laminar viscosity coefficient, and µ 2 is the turbulent viscosity coefficient. By comparing Equations (10)-(12), Equations (10) and (11) can be regarded as a special form of Equation (12), because slurry flows consists of water and fine particles (liquid phase) and coarse particles (solid phase). In the Bingham model, since the fluid turbulence is not considered, the turbulent viscosity coefficient is µ 2 = 0. In the expansion model, the inter-particle frictional force is negligible relative to the particles' collision, so the second term in Equation (11) is zero. By substituting single components of Equations (9) and (10), the second term of right hand side in Equation (7) can be written as follows: By substituting Equation (13) into Equation (14), the viscous term in the x, y direction can be represented by the following discretization: It can be seen from Equation (15) that when deformation |D| is relatively small, yield stress (particle support force) has a greater impact on the fluid's acceleration, but there should be an upper limit to this effect in order to prevent excessive acceleration which could cause the local instability of the fluid investigated. In existing studies [16], a lower limit is usually set for |D|. When |D| is less than this lower limit, the relationship between stress and |D| satisfies linearity: where σ is the limiting factor.
In [16], Hosseini considers that the viscosity of "solid zone" fluid is much greater than that of the main fluid (100 times). In this study, although the turbulent stress term 1 2 µ 2 |D|D is introduced, it can be ignored in the region of low velocity. By calibrating this value against experimental results obtained for this study, 200µ 1 was the limited factor selected. Figure 1 shows the facility where the experiments for the numerical validation of the method previously described in Section 2 were conducted. This facility is located in the Debris Flow Observation and Research Station at Jiangjia Gully, the largest field research center in China-also known as "the debris flow museum". For this study, water is stored upstream to the material prepared under multiple configurations and water is always released by opening a steel gate at the velocity of 3 m/s. The flume can then be divided into two parts: (i) a steep upstream reach (6.0 m long and the chute slope can vary from 15 • to 40 • (always set up at 25 • for these experiments); and (ii) a flat-bottomed downstream section (3.0 m long). The slope of the flume's bed can be manually adjusted.

Experimental Setup and Boundary Conditions
Details regarding the experimental procedure conducted can be found in [6]. For this study, three slurries were testes with densities ρ = 1400 kg/m 3 ; ρ = 1500 kg/m 3 and ρ = 1600 kg/m 3 . Different layers patterns were selected according to the different configurations displayed in Figure 2. By adding water to the flume, when the water level reached the height of the mixing fluid, the front-end steel gate of the mixing area was released at a speed of 3 m/s. In order to maintain the driving force of the mixtures, the water level behind the mixtures was kept at h = 0.2 m during the experiment.    There are three configurations: Figure 2a shows the configuration a, with coarse particles as a top layer and fine sediment positioned below them; Figure 2b shows the configuration b, a realistic configuration with the original sediments collected from the Jiangjiagou Valley, fully mixed and distributed along the entire measurement area; Figure 2c shows the configuration c, with coarse particles distributed at the bottom and fine grains on the top of them. For this study, three types of liquid slurries were considered: (i) ρ = 1400 kg/m 3 , (ii) ρ = 1500 kg/m 3 , and (iii) ρ = 1600 kg/m 3 . The slurry rheology coefficient was measured by the MCR301 advanced rotary rheometer manufactured by Anton Paar, Austria [39]. Values of viscosity µ for the fluids simulated are displayed in Table 1. The particle size distribution curves are shown in Figure 3, Ψ stands for the percentage of accumulated mass of particles and d denotes particle size. The particle size distribution curves are shown in Figure 3, Ψ stands for the percentage of accumulated mass of particles and d D denotes particle size.

Simulation and Analysis
By simulating the movement processes of different viscous debris flows, a series of experiments has been completed where free surface heights, fluid velocities, pressures, and shear deformations associated with the movement of the fluid were measured. The numerical simulation was carried out to generate the experimental conditions previously described (Figure 4). The numerical simulation

Simulation and Analysis
By simulating the movement processes of different viscous debris flows, a series of experiments has been completed where free surface heights, fluid velocities, pressures, and shear deformations associated with the movement of the fluid were measured. The numerical simulation was carried out to generate the experimental conditions previously described (Figure 4). The numerical simulation replicated N = 113,054 particles, particle spacing d p = 0.0025 m, solid particle density ρ s = 2200 kg/m 3 , thickness of solid phase h s = 0.1 m and thickness of liquid phase h l = 0.1 m for configuration a and c in Figure 2. Similar to the tests conducted on the experimental facility, three different viscosity coefficients for the liquid phases (as shown in Table 1) were selected in the numerical simulation. The inflow conditions were the same as those applied experimentally, and the water level as driving force of the debris flow was kept at 0.2 m for each entire simulation. replicated N = 113,054 particles, particle spacing dp = 0.0025 m, solid particle density ρs = 2200 kg/m 3 , thickness of solid phase hs = 0.1 m and thickness of liquid phase hl = 0.1 m for configuration a and c in Figure 2. Similar to the tests conducted on the experimental facility, three different viscosity coefficients for the liquid phases (as shown in Table 1) were selected in the numerical simulation. The inflow conditions were the same as those applied experimentally, and the water level as driving force of the debris flow was kept at 0.2 m for each entire simulation.  Figure 5 represents the free surface values recorded at different times for each of the tests conducted in Table 1, after the debris flow initiation generated by the release of water upstream. The x-axis represents the length L of the debris flow, and the y-axis represents the height H of the debris flow. After 0.7 s, the distance reached by the debris flow (considering all the configurations) is within the range 2.48-2.66 m, and the maximum velocity range is 4.65-5.12 m/s. Authors noticed that when   Table 1, after the debris flow initiation generated by the release of water upstream. The x-axis represents the length L of the debris flow, and the y-axis represents the height H of the debris flow. After 0.7 s, the distance reached by the debris flow (considering all the configurations) is within the range 2.48-2.66 m, and the maximum velocity range is 4.65-5.12 m/s. Authors noticed that when the debris flow has similar viscous properties but for the vertical distribution of different particles, the differences in velocities are not so significant and can considered almost negligible in most of the cases. However, the maximum velocity recorded for configuration b (Shown in Figure 2b correspond to tests 2, 5 and 8 in Figure 5) is similar to the one recorded for configuration a (Shown in Figure 2a and corresponding to tests 1, 4, and 7 in Figure 5), while configuration c (Shown in Figure 2c and corresponding to tests 3, 6, and 9 in Figure 5) was characterized by higher values of velocities and elevations measured. By comparing the shapes of head under different vertical distributions, it was found that for the tests conducted in Table 1, free surface values measured for configuration b fluctuate more than in configuration a and c, demonstrating that this scenario is typical of intermittent debris flows.

Characterization of Intermittent Debris Flow
In order to study the causes of this phenomenon, the characteristics of the fluid movement processes associated to configuration b were analyzed. Figure 6 shows the evolution of the solidliquid phase at different locations simulated numerically. Figure 6a shows that in the horizontal region, most of the particles still retain under the laminar form. When moving into the upstream part of the sloped section, the fluid height decreases and the liquid phase group is stretched, as shown in Figure 6b. Then, due to the slope, velocity increases while the fluid height decreases, and different layers of liquid and solid particles will appear almost as parallel mixing within the entire width of the debris flow, as shown in Figure 6c. At this stage, the altering layers interact changing continuously positions demonstrating that mixing processes are taking place and when the mixing is finally completed, the fluctuation amplitude reduce becoming more stable, while the influenced range of the fluctuations can spread over a longer length, as shown in Figure 6d. Figure 6e shows the effect of the liquid phase on the height of the debris flow. It is clear that when liquid particles accumulate due to the mixing phenomena (highlighted as circles in Figure 6e), there is a correspondent decrease of the height of the debris flow (pointed out using arrows in Figure 6e). This inverse relationship is very interesting especially because it demonstrates how the gathering and accumulation of liquid particles

Characterization of Intermittent Debris Flow
In order to study the causes of this phenomenon, the characteristics of the fluid movement processes associated to configuration b were analyzed. Figure 6 shows the evolution of the solid-liquid phase at different locations simulated numerically. Figure 6a shows that in the horizontal region, most of the particles still retain under the laminar form. When moving into the upstream part of the sloped section, the fluid height decreases and the liquid phase group is stretched, as shown in Figure 6b. Then, due to the slope, velocity increases while the fluid height decreases, and different layers of liquid and solid particles will appear almost as parallel mixing within the entire width of the debris flow, as shown in Figure 6c. At this stage, the altering layers interact changing continuously positions demonstrating that mixing processes are taking place and when the mixing is finally completed, the fluctuation amplitude reduce becoming more stable, while the influenced range of the fluctuations can spread over a longer length, as shown in Figure 6d. Figure 6e shows the effect of the liquid phase on the height of the debris flow. It is clear that when liquid particles accumulate due to the mixing phenomena (highlighted as circles in Figure 6e), there is a correspondent decrease of the height of the debris flow (pointed out using arrows in Figure 6e). This inverse relationship is very interesting especially because it demonstrates how the gathering and accumulation of liquid particles tends to appear towards the bottom side of the debris flow layer.  On this basis, the relationship between the moisture content φ (the amount of liquid particles divided by the amount of solid particles, N l /N s ), the kinetic energy of particles E k and the height of the free surface H (related to the potential energy of particles E p ) were calculated for the tests conducted for configuration b. As shown in Figure 7a, the height of the free surface H decreases as the debris flow develops. There is a noticeable correlation between the fluctuation of the free surface associated with the fluctuation of the moisture content. In the regions of L = 1.00-1.38 m and L = 1.82−2.04 m, the height of free surface decreases linearly, and in these two regions, the water content remains in the range 0.2-0.65. The points that obviously exceed this threshold are L = 0.90, L = 1.52, L = 1.6, L = 1.74, and L = 1.80−1.84, and the height of free surface is different from that of linear decline in these areas or vicinity. When the moisture content is within the range 0.20-0.65, the free surface of the debris flow is characterized by a linear change, but when the moisture content exceeds this range, it generates an impact on the free surface. From Figure 7b, it can be seen that there is a more obvious negative correlation between the particles kinetic energy Ek and the moisture content ϕ. To almost every peak of the kinetic energy Ek (highlighted as green circles in Figure 7b) calculated corresponds a peak of the moisture content ϕ (highlighted as blue circles in Figure 7b), which indicates that kinetic energy Ek and moisture content ϕ interact directly. However, this effect can only be assigned to small-scale portions of the particle kinetic energy fluctuations. Looking at Figure 7b, at the location of L = 1.74 m, the moisture content value corresponds to ϕ = 0.7619 and it is the maximum value measured in this region, and the corresponding kinetic energy Ek records its minimum value. But because of the large kinetic energy of the particles recorded in this region, the corresponding kinetic energy Ek = 5.4848 J is still higher than that recorded at the position of L = 1.64 m in the adjacent one Ek = 0.30263 J. From Figure 7b, it can be seen that there is a more obvious negative correlation between the particles kinetic energy E k and the moisture content φ. To almost every peak of the kinetic energy E k (highlighted as green circles in Figure 7b) calculated corresponds a peak of the moisture content φ (highlighted as blue circles in Figure 7b), which indicates that kinetic energy E k and moisture content φ interact directly. However, this effect can only be assigned to small-scale portions of the particle kinetic energy fluctuations. Looking at Figure 7b, at the location of L = 1.74 m, the moisture content value corresponds to φ = 0.7619 and it is the maximum value measured in this region, and the corresponding kinetic energy E k records its minimum value. But because of the large kinetic energy of the particles recorded in this region, the corresponding kinetic energy E k = 5.4848 J is still higher than that recorded at the position of L = 1.64 m in the adjacent one E k = 0.30263 J. Finally, the energy conversion curves of fluids with different viscous coefficients were inspected and confronted, as shown in Figure 9. It was found that the gravitational potential energy (E p = mgH) and the total energy (E 0 = E k + E p ) of fluids decreases at a similar rate. The difference between three fluids is mainly reflected on kinetic energies. When comparing the set of fluids with the smallest density (ρ = 1400 kg/m 3 , µ 0 = 0.00004, µ 1 = 0.0048, and µ 2 = 0.0197), results shows that velocity values increase from t = 0.0 s up to t = 0.6 s, reaching almost the highest values, and then the kinetic energy of the three fluids tends to be equal. As the time progresses, the same order appears again in the kinetic energy magnitude arrangement, which is E k,ρ = 1400 kg/m 3 > E k,ρ = 1500 kg/m 3 > E k,ρ = 1600 kg/m 3 . This phenomenon is also due to the stronger fluctuation of the less dense fluids and these effects caused by different viscous fluids on debris flow array and collision, and friction forces on debris flow movement, will require further investigation in the future. Finally, the energy conversion curves of fluids with different viscous coefficients were inspected and confronted, as shown in Figure 9. It was found that the gravitational potential energy (Ep = mgH) and the total energy (E0 = Ek + Ep) of fluids decreases at a similar rate. The difference between three fluids is mainly reflected on kinetic energies. When comparing the set of fluids with the smallest density (ρ = 1400 kg/m 3 , 0 μ = 0.00004, 1 μ = 0.0048, and 2 μ = 0.0197), results shows that velocity values increase from t = 0.0 s up to t = 0.6 s, reaching almost the highest values, and then the kinetic energy of the three fluids tends to be equal. As the time progresses, the same order appears again in the kinetic energy magnitude arrangement, which is Ek,ρ = 1400 kg/m 3 > Ek,ρ = 1500 kg/m 3 > Ek,ρ = 1600 kg/m 3 . This phenomenon is also due to the stronger fluctuation of the less dense fluids and these effects caused by different viscous fluids on debris flow array and collision, and friction forces on debris flow movement, will require further investigation in the future.

Summary and Conclusions
Debris flows are a natural phenomenon causing a lot of economical and human losses worldwide. Due to their nature, debris flows travel long distances at high speeds, and the time-space evolution of the relationship between soil and water content strongly affects the propagation stage. Thus, a quantitative modeling of this phenomenon is crucial to design strategies to be adopted to reduce the negative impacts. This paper contributes to this topic through the use of a SPH model investigating multiple combinations of fine and coarse particles with water content.
The SPH model was used to simulate tests conducted in the experimental facility that is located in the Debris Flow Observation and Research Station at Jiangjia Gully, the largest field research center in China-also known as "the debris flow museum". As previously demonstrated, the SPH model is capable of properly reproducing the main characteristics of debris flows (propagation height and velocity, and more importantly to correctly simulate the time-space evolution of solid and liquid particles during the whole process from initiation to propagation over an impervious/permeable bottom boundary).

Summary and Conclusions
Debris flows are a natural phenomenon causing a lot of economical and human losses worldwide. Due to their nature, debris flows travel long distances at high speeds, and the time-space evolution of the relationship between soil and water content strongly affects the propagation stage. Thus, a quantitative modeling of this phenomenon is crucial to design strategies to be adopted to reduce the negative impacts. This paper contributes to this topic through the use of a SPH model investigating multiple combinations of fine and coarse particles with water content.
The SPH model was used to simulate tests conducted in the experimental facility that is located in the Debris Flow Observation and Research Station at Jiangjia Gully, the largest field research center in China-also known as "the debris flow museum". As previously demonstrated, the SPH model is capable of properly reproducing the main characteristics of debris flows (propagation height and velocity, and more importantly to correctly simulate the time-space evolution of solid and liquid particles during the whole process from initiation to propagation over an impervious/permeable bottom boundary).
Based on the theory of solid-liquid two phase flows, the viscous term in the SPH model was modified to make it suitable for nonhomogeneous viscous debris flows. It was found that the denser the fluid is, the greater are the yield stress and the turbulent viscous coefficient. However for laminar viscous coefficients, the fluid with the density of ρ = 1500 kg/m 3 has the largest values. The results obtained can be summarized as follows:

1.
By comparing the shape and velocity of debris flows under different configurations, it was found that the vertical distribution of particles played a very important role in debris flow fluctuation, with a greater influence than on the viscous coefficient. The third configuration with mixed fine and coarse particles showed to fluctuate more violently, and this outcome confirmed one of the main assumptions for intermittent debris flows.

2.
By analyzing the characteristics of the fluid movement processes, it was found that when the two layers (fine and coarse particles) are mixed with the water, liquid particles tended to gather towards the bottom side of the debris flow causing a correspondent decrease of height. However, this effect could only be observed at small-scale areas. The potential energy was greatly affected by the magnitude of the moisture content, while the kinetic energy was significantly affected by the derivative of moisture content in the L direction. 3.
The differences of the energy conversion curves associated to different viscous coefficients were mainly noticed in kinetic energies. Fluids with smaller densities exhibited higher initiation velocities and higher fluctuations values.