Modeling the Slump-Type Landslide Tsunamis Part I: Developing a Three-Dimensional Bingham-Type Landslide Model

Featured Application: This paper aims to develop a landslide model to simulate slump-type landslide tsunamis. The modiﬁed Bi-viscous model can be used to simulate the landslide movements in both the in-land area and ocean area. The model is able to describe the developments of the landslide as well as the slip surface. In addition, the model can be applied to hazard asscessments. Abstract: This paper incorperates Bingham and bi-viscosity rheology models with the Navier–Stokes solver to simulate the dynamics and kinematics processes of slumps for tsunami generation. The rheology models are integrated into a computational ﬂuid dynamics code, Splash3D, to solve the incompressible Navier–Stokes equations with volume of ﬂuid surface tracking algorithm. The change between un-yield and yield phases of the slide material is controlled by the yield stress and yield strain rate in Bingham and bi-viscosity models, respectively. The integrated model is carefully validated by the theoretical results and laboratory data with good agreements. This validated model is then used to simulate the benchmark problem of the failure of the gypsum tailings dam in East Texas in 1966. The accuracy of predicted ﬂood distances simulated by both models is about 73% of the observation data. To improve the prediction, a ﬁxed large viscosity is introduced to describe the un-yield behavior of tailings material. The yield strain rate is obtained by comparing the simulated inundation boundary to the ﬁeld data. This modiﬁed bi-viscosity model improves not only the accuracy of the spreading distance to about 97% but also the accuracy of the spreading width. The un-yield region in the modiﬁed bi-viscosity model is sturdier than that described in the Bingham model. However, once the tailing material yields, the material returns to the Bingham property. This model can be used to simulate landslide tsunamis.


Introduction
Tsunamis are potentially deadly and destructive sea waves. Most of the tsunamis are formed as a result of submarine earthquakes and submarine landslides. These landslides, in turn, are often triggered by earthquakes or volcanic eruptions [1]. Over the past 20 years, catastrophic tsunamis in Papua New Guinea (1998), Indian Ocean (2004), Japan (2011), Palu Bay Indonesia (2018), and Anak Krakatau Indonesia (2018) have driven major advances in understanding of earthquakes and submarine landslides as tsunami sources [2]. In fact, submarine landslides have become suspects in the creation
The rheological properties of BM can be presented as [38,39]: where µ A is the viscosity of the un-yield region, µ B is the viscosity of the yield zone, τ y is the yield stress, and . γ y is the yield strain rate, . γ ij = ∂ . Even with the simplicity of the BM, the stress is still indeterminate in the un-yield region, which means the exact shape and location of the yield surface(s) cannot be determined [40]. To remedy this drawback in the present work, the conventional bi-viscosity model (CBM) [41] was adopted. This idea allows a small deformation to occur in the un-yield region(s) by treating it as an extremely high viscosity fluid. In the yield region, the material is considered a Bingham fluid. This method makes it possible for the stress to be computable in the whole domain, including the un-yield region so that the location of the yield surface can be easily determined [40].
The rheological properties of CBM can be presented as [42,43]: Mathematically speaking, when . γ y approache zero, the CBM approaches BM. If the chosen . γ y is sufficiently small, we can practically replace the un-yield region viscosity with a higher viscosity. This guarantees that a viscous solver can handle the determination of the shape and the location of the plug surface [40].
However, a mud material will become sturdy when experiencing compaction or tamping processes experiencing. To describe the sturdy behavior in the plug zone, a larger µ A and a larger . γ y are required.
The larger µ A plays a role of keeping the rigid shape. The larger . γ y indicates that the material can sustain a large shear stress without deformation. To achieve that effect, the modified bi-viscosity model (MBM) is born and written as: In this model, the yield stress τ y and yield viscosity µ B of the mud material are exponentially dependent on material concentration [44]. The detailed descriptions are added in Section 5.2. To present the un-yield behavior, µ A is chosen to be infinite based on the suggestions of Assier Rzadkiewicz [34], Taibi [45], and Yu [28]. In this paper, the infinite number of viscosity µ A = 10 10 Pa s is chosen by the sensitivity analysis. The values of yield strain rate . γ y are also discussed in Section 5.2.
. γ y = 0.2 s −1 is adopted by a sensitivity analysis to illustrate the deformation in MBM.
The viscoplastic models, BM, CBM, and MBM, were coupled with the Splash3D model. The Splash3D model was renovated from the open-source software, Truchas, which was originally developed by Los Alamos National Laboratory [46]. The original program can simulate the incompressible flows with multi-fluid interfaces. The code solves three-dimensional continuity and Navier-Stokes equations by adopting the projection method [36,37,47] and the finite volume discretization method [48]. The Splash3D model was enhanced with several hydrodynamic modules such as the large eddy simulation (LES) turbulence module [40,49], and the moving-solid module [50] to deal with breaking waves and wave-obstacle interaction problems. Readers are encouraged to read the reference Chu [50] for the detailed numerical algorithm. In this study, the Splash3D is developed with the rheological model to solve mudflow problems. The fundamental governing equations are continuity and momentum equations: Appl. Sci. 2020, 10, 6501 5 of 21 where the subscripts i, j = 1, 2, 3 represent the x, y, z directions respectively; t is the time, u is the velocity; P is the pressure; the over-bar represents the spatially filtered value [51]; g is the gravitational acceleration; ρ is the density; and µ e is the effective viscosity. Under the influence of gravity, tailings flow out of the break at a high speed and might be in a turbulent state [52]. In this study, the large eddy simulation (LES) [53] is adopted to address the turbulence effect. The effective viscosity µ e is defined as: where µ . γ is the rheological viscosity of mud and µ t is the viscosity of the sub-grid scale turbulence. The Smagorinsky model [53] relates the residual stress to the rate of filtered strain. Based on the dimensional analysis, the subgrid-scale eddy viscosity is modeled as: where l S is the Smagorinsky length scale, which is a product of the Smagorinsky coefficient C s and the filter width ∆; S is the characteristic filtered rate of strain: In general, C s varies from 0.1 to 0.2 in different flows. The present simulations use a value of 0.15. ∆ is the filter width. Infinite volume discretization ∆ is the grid size: where ∆x 1 , ∆x 2 , ∆x 3 are the three components of the grid lengths. In the present model, the tailings fluid is treated as a single homogeneous mud material. The tailings fluid and air are assumed to be two incompressible and non-immiscible fluids. The free-surface between the tailings fluid and air is tracked by the volume of fluid (VOF) method [54]. The volume fraction, f m , is used to describe the fraction of the mth material in each cell. The volume fraction f m varies in [0,1] and should sum to unity everywhere. f m = 1 if the cell is fully occupied by the mth material; 0 < f m < 1 if the cell contains the interfaces of the mth material; f m = 0 if the cell contains no mth material. The VOF equation is given by Equation (10): In this study, the interfaces between different materials are solved by the VOF method. Because the fundamental assumption of the VOF method is that each fluid is immiscible, there is no difference in terms of the VOF equations between subaerial or submerged slumps. For both the subaerial and submerged (underwater) slumps, the VOF method is briefly described here [50,55,56].
Two stability conditions of dt need to be satisfied while solving Navier-Stokes equations: Max µ e f f (12) where dt c is the time step restricted by the advection term, C r is Courant number, which is defined as C r = Max(|u|)dt/dl, dl is the measure of the cell size, dt µ is the time step restricted by the diffusion term, V µ is the viscous number, which is defined as V µ = Max µ e f f dt/(dl) 2 .
When the viscosity is large in the un-yield zone [50,57,58], the time-step dt µ is very small, which might lead to the numerical divergence or even a crash in the solution procedure [59]. However, this small time-step restriction can be relaxed by adopting the implicit scheme. The viscous implicitness θ is used to calculate the velocity u at the time level n, and u θ = (1 − θ)u n + θu n+1 . In this study, θ is given as unity, which implies a fully implicit treatment, and dt µ is no longer restricted by Equation (12).
The algorithm of Splash3D solves equations of conservation of mass and momentum for any number of immiscible, incompressible fluids, and tracks the interfaces between them. Except for adding a water body to the computational domain in Part II to simulate slump-type landslide tsunamis, no additional numerical adjustment is required from Part I to Part II. Therefore, the findings in Part I can be applied to the underwater modeling in Part II. To simplify the complexness from the slump, water, and air, the slump-type landslide in the dry-land area, FGT66, is adopted in Part I for a better understanding of the model characteristics. The same numerical model will be applied to study the FGT66 in Part I and slump-type landslide tsunamis in Part II.
To summarize the numerical method, this paper adopts Splash3D model to solve the conservation and Navier-Stokes equations for any number of immiscible, incompressible fluids. The LES turbulence model with Smagorinsky closure is used to add the effect from turbulence flow field. The VOF method is adopted to track the interfaces between each fluid including slumps. The new contributions from this study in terms of the numerical modeling are incorporating the BM, CBM, and MBM into the Splash3D model with the implicit scheme to solve the slump material with a large viscosity number.

Validation
Two cases of mudflows are simulated for the model validation. The results are compared with both analytical solutions and laboratory experiment data.

Bingham Fluid Driven by Pressure Gradients
Byron-Bird [49] derived analytical solutions for the Bingham flow in a channel, driven by a pressure gradient P 0 − P L . The channel was depicted as the length L and the width 2B. The no-slip boundary condition was applied to the surfaces of the channel (Figure 1).
Appl. Sci. 2020, 10, x FOR PEER REVIEW  6 of 22 is given as unity, which implies a fully implicit treatment, and is no longer restricted by Equation (12).
The algorithm of Splash3D solves equations of conservation of mass and momentum for any number of immiscible, incompressible fluids, and tracks the interfaces between them. Except for adding a water body to the computational domain in Part II to simulate slump-type landslide tsunamis, no additional numerical adjustment is required from Part I to Part II. Therefore, the findings in Part I can be applied to the underwater modeling in Part II. To simplify the complexness from the slump, water, and air, the slump-type landslide in the dry-land area, FGT66, is adopted in Part I for a better understanding of the model characteristics. The same numerical model will be applied to study the FGT66 in Part I and slump-type landslide tsunamis in Part II.
To summarize the numerical method, this paper adopts Splash3D model to solve the conservation and Navier-Stokes equations for any number of immiscible, incompressible fluids. The LES turbulence model with Smagorinsky closure is used to add the effect from turbulence flow field. The VOF method is adopted to track the interfaces between each fluid including slumps. The new contributions from this study in terms of the numerical modeling are incorporating the BM, CBM, and MBM into the Splash3D model with the implicit scheme to solve the slump material with a large viscosity number.

Validation
Two cases of mudflows are simulated for the model validation. The results are compared with both analytical solutions and laboratory experiment data.

Bingham Fluid Driven by Pressure Gradients
Byron-Bird [49] derived analytical solutions for the Bingham flow in a channel, driven by a pressure gradient 0 − . The channel was depicted as the length and the width 2 . The no-slip boundary condition was applied to the surfaces of the channel (Figure 1).
if 0 = 0.0, 0 = 0.0, and the material is a Newtonian fluid. Four cases were proposed to validate the models with Equations (13) and (14). They are one Newtonian case and three Bingham cases with different parameters such as channel length , channel width 2 , one end's pressure 0 , Bingham viscosity , and yield stress 0 . Figure 2 shows good agreements between the theoretical and numerical results of the all runs. One of the important features of a Bingham fluid is the plug zone (Figure 2b-d), which cannot be seen in the Newtonian fluid ( Figure 2a). Note that the velocity of a Bingham fluid is constant in the plug region. In this The yield surface is located at x = x 0 where x 0 = τ 0 L P 0 −P L . The velocities in the plug region v > z , and in the liquefied region v < z are: if τ 0 = 0.0, x 0 = 0.0, and the material is a Newtonian fluid. Four cases were proposed to validate the models with Equations (13) and (14). They are one Newtonian case and three Bingham cases with different parameters such as channel length L, channel Appl. Sci. 2020, 10, 6501 7 of 21 width 2B, one end's pressure P 0 , Bingham viscosity µ B , and yield stress τ 0 . Figure 2 shows good agreements between the theoretical and numerical results of the all runs. One of the important features of a Bingham fluid is the plug zone (Figure 2b-d), which cannot be seen in the Newtonian fluid ( Figure 2a). Note that the velocity of a Bingham fluid is constant in the plug region. In this region, the rate of change of velocity (strain rate) is equal to zero. In the liquefied region, the strain rate is greater than zero and the stress-strain relation of the fluid is dependent on the plastic viscosity µ B . These figures demonstrate that the present numerical model can accurately describe the rheological behavior of Bingham fluids.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 7 of 22 region, the rate of change of velocity (strain rate) is equal to zero. In the liquefied region, the strain rate is greater than zero and the stress-strain relation of the fluid is dependent on the plastic viscosity . These figures demonstrate that the present numerical model can accurately describe the rheological behavior of Bingham fluids.

Spreading of Bingham Fluid on an Inclined Plane
The validation of the spreading of Bingham fluid on an inclined plane is set up based on the experiment of Liu [12]. Kaolinite was mixed with tap water to represent the mud. The mud was put in a reservoir. When the adjustable gate was opened, the mud flowed down onto an inclined dry bed with the inclined angle θ = 0.9 • . In this validation case, the openness of the gate was H = 0.0051 m. The fluid density ρ = 1106 kg m −3 , the yield stress, τ y = 0.875 Pa, the viscosity of the plug zone, µ A = 10 10 Pa s, the viscosity of liquefied zone, µ B = 0.034 Pa s. A two-dimensional numerical domain was set up as 3.5 × 0.2 m ( Figure 3) and was discretized into a regular mesh with grid size dx = 2.3, dz = 2.0 mm. Figure 4a shows the spreading of mud on an inclined plane versus time. Figure 4b shows that the numerical result of Bingham model matches well with the theoretical solution as well as the experimental data from Liu [12]. The mudflow develops into a self-similar front when time t > 8.0 s. Because of the yield stress, the free surface needs not to be horizontal when the mud fluid is in static equilibrium, nor parallel to the plane bed when it reaches a steady-state. The mud front, like a steady gravity current, eventually advances at a constant speed with the same profile when there is a steady upstream discharge of mud [12]. The numerical results present a similar pattern of analytical solutions to that in Liu [12].

Spreading of Bingham Fluid on an Inclined Plane
The validation of the spreading of Bingham fluid on an inclined plane is set up based on the experiment of Liu [12]. Kaolinite was mixed with tap water to represent the mud. The mud was put in a reservoir. When the adjustable gate was opened, the mud flowed down onto an inclined dry bed with the inclined angle = 0.9°. In this validation case, the openness of the gate was = 0.0051 m. The fluid density = 1106 kg m −3 , the yield stress, = 0.875 Pa, the viscosity of the plug zone, = 10 10 Pa s, the viscosity of liquefied zone, = 0.034 Pa s. A two-dimensional numerical domain was set up as 3.5 × 0.2 m ( Figure 3) and was discretized into a regular mesh with grid size dx = 2.3, dz = 2.0 mm. Figure 4a shows the spreading of mud on an inclined plane versus time. Figure 4b shows that the numerical result of Bingham model matches well with the theoretical solution as well as the experimental data from Liu [12]. The mudflow develops into a self-similar front when time t > 8.0 s. Because of the yield stress, the free surface needs not to be horizontal when the mud fluid is in static equilibrium, nor parallel to the plane bed when it reaches a steady-state. The mud front, like a steady gravity current, eventually advances at a constant speed with the same profile when there is a steady upstream discharge of mud [12]. The numerical results present a similar pattern of analytical solutions to that in Liu [12].   [12]. The numerical results are simulated by the Bingham model.

Numerical Setup
The present numerical model was applied to simulating the failure of the gypsum tailings dam in East Texas in 1966 (FGT66). The reservoir was a rectangular shape and reached a height of 11 m when the failure took place. The slide was triggered by seepage at the toe of the embankment. An estimated 80,000-130,000 m 3 of gypsum was released in this flow failure. The released material traveled about 300 m before it came to a stoppage, with an average velocity of 2.5-5.0 m/s [60]. In this paper, the numerical setup was composed based on the geometry reported by Jeyapalan [60], shown in Figure 5. The size of the tailings reservoir was 280 × 110 × 11 m, and the breach was 120 m long. Appl. Sci. 2020, 10, x FOR PEER REVIEW 8 of 22

Spreading of Bingham Fluid on an Inclined Plane
The validation of the spreading of Bingham fluid on an inclined plane is set up based on the experiment of Liu [12]. Kaolinite was mixed with tap water to represent the mud. The mud was put in a reservoir. When the adjustable gate was opened, the mud flowed down onto an inclined dry bed with the inclined angle = 0.9°. In this validation case, the openness of the gate was = 0.0051 m. The fluid density = 1106 kg m −3 , the yield stress, = 0.875 Pa, the viscosity of the plug zone, = 10 10 Pa s, the viscosity of liquefied zone, = 0.034 Pa s. A two-dimensional numerical domain was set up as 3.5 × 0.2 m ( Figure 3) and was discretized into a regular mesh with grid size dx = 2.3, dz = 2.0 mm. Figure 4a shows the spreading of mud on an inclined plane versus time. Figure 4b shows that the numerical result of Bingham model matches well with the theoretical solution as well as the experimental data from Liu [12]. The mudflow develops into a self-similar front when time t > 8.0 s. Because of the yield stress, the free surface needs not to be horizontal when the mud fluid is in static equilibrium, nor parallel to the plane bed when it reaches a steady-state. The mud front, like a steady gravity current, eventually advances at a constant speed with the same profile when there is a steady upstream discharge of mud [12]. The numerical results present a similar pattern of analytical solutions to that in Liu [12].   [12]. The numerical results are simulated by the Bingham model.

Numerical Setup
The present numerical model was applied to simulating the failure of the gypsum tailings dam in East Texas in 1966 (FGT66). The reservoir was a rectangular shape and reached a height of 11 m when the failure took place. The slide was triggered by seepage at the toe of the embankment. An estimated 80,000-130,000 m 3 of gypsum was released in this flow failure. The released material traveled about 300 m before it came to a stoppage, with an average velocity of 2.5-5.0 m/s [60]. In this paper, the numerical setup was composed based on the geometry reported by Jeyapalan [60], shown in Figure 5. The size of the tailings reservoir was 280 × 110 × 11 m, and the breach was 120 m long.

Numerical Setup
The present numerical model was applied to simulating the failure of the gypsum tailings dam in East Texas in 1966 (FGT66). The reservoir was a rectangular shape and reached a height of 11 m when the failure took place. The slide was triggered by seepage at the toe of the embankment. An estimated 80,000-130,000 m 3 of gypsum was released in this flow failure. The released material traveled about 300 m before it came to a stoppage, with an average velocity of 2.5-5.0 m/s [60]. In this paper, the numerical setup was composed based on the geometry reported by Jeyapalan [60], shown in Figure 5. The size of the tailings reservoir was 280 × 110 × 11 m, and the breach was 120 m long. The center cross-section of the breach was located at y = 220 m. The computational domain (510 in length, 400 in width, and 12 m in height) was discretized into a uniform mesh with a grid size dx = 2.0, dy = 2.0, dz = 1.0 m. The number of the grid was 612,000. The bottom boundary (at z = 0 m) was a no-slip boundary condition. The downstream (x = 400 m) and lateral boundaries (y = 0 and y = 400 m) were free-slip walls. The downstream and lateral boundaries would not affect the simulation results because the domain was set to be much larger than the predicted tailing pattern. The gypsum tailings material was expected as a Bingham material. Based on the parameters reported by Jeyapalan [60], Pastor [61], and Chen [20], the yield stress of the tailings was τ y = 10 3 Pa, the viscosity of the liquefied zone was µ B = 50 Pa s, and the density was ρ = 1400 kg m −3 . The viscosity of the plug zone was suggested to be infinite (e.g., µ A = 10 10 Pa s) by Assier Rzadkiewicz [34], Taibi [45], and Yu [28].
Appl. Sci. 2020, 10, x FOR PEER REVIEW 9 of 22 The center cross-section of the breach was located at y = 220 m. The computational domain (510 in length, 400 in width, and 12 m in height) was discretized into a uniform mesh with a grid size dx = 2.0, dy = 2.0, dz = 1.0 m. The number of the grid was 612,000. The bottom boundary (at z = 0 m) was a no-slip boundary condition. The downstream (x = 400 m) and lateral boundaries (y = 0 and y = 400 m) were free-slip walls. The downstream and lateral boundaries would not affect the simulation results because the domain was set to be much larger than the predicted tailing pattern. The gypsum tailings material was expected as a Bingham material. Based on the parameters reported by Jeyapalan [60], Pastor [61], and Chen [20], the yield stress of the tailings was = 10 3 Pa, the viscosity of the liquefied zone was = 50 Pa s, and the density was = 1400 kg m −3 . The viscosity of the plug zone was suggested to be infinite (e.g., = 10 10 Pa s) by Assier Rzadkiewicz [34], Taibi [45], and Yu [28].  Figure 6b,c is produced by CBM with ̇= 1 × 10 −4 s −1 and ̇= 2 × 10 −1 s −1 , respectively. In the CBM results, a high viscosity number is used to represent the un-yield phase, and a low viscosity number is used to represent the yield phase. The tailings shapes are similar between results of ̇= 1 × 10 −4 s −1 (Figure 6b) and ̇= 2 × 10 −1 s −1 (Figure 6c). They are also similar to those in the BM result (Figure 6a). However, the mud profiles at the center streamline are slightly different. The 'ridge' in the center of the breach (Figure 6b) is replaced by a relatively smooth hump (Figure 6c). In Figure 6, the results from CBM were nearly identical to those from BM. The flood distances predicted by BM and CBM were about 220 m. Compared to the field observation [60], the error of the predicted flood distance was about 27%. The field photo (Figure 7b) shows that the flood boundary was longer and narrower than the simulated results from BM and CBM. This might result from the sturdy behavior in the un-yield region. This sturdy behavior can be reached by increasing the un-yield viscosity and yield strain ̇. In this study, the un-yield viscosity = 10 10 Pa s was chosen by the suggestion of Assier Rzadkiewicz [34], Taibi [34], and Yu [28]. The yield strain rate was chosen to be ̇= 2 × 10 −1 s −1 by matching the flood distance. Figure 7a shows the  (Figure 6c). They are also similar to those in the BM result (Figure 6a). However, the mud profiles at the center streamline are slightly different. The 'ridge' in the center of the breach (Figure 6b) is replaced by a relatively smooth hump (Figure 6c). In Figure 6, the results from CBM were nearly identical to those from BM. The flood distances predicted by BM and CBM were about 220 m. Compared to the field observation [60], the error of the predicted flood distance was about 27%. The field photo (Figure 7b) shows that the flood boundary was longer and narrower than the simulated results from BM and CBM. This might result from the sturdy behavior in the un-yield region. This sturdy behavior can be reached by increasing the un-yield viscosity µ A and yield strain . γ y . In this study, the un-yield viscosity µ A = 10 10 Pa s was chosen by the suggestion of Assier Rzadkiewicz [34], Taibi [34], and Yu [28]. The yield strain rate was chosen to be . γ y = 2 × 10 −1 s −1 by matching the flood distance. Figure 7a shows the simulation of the deposited tailings from MBM. The flood distance at the freezing time t = 110 s was 310 m, which was 97% accurate to the filed data [60]. The result from MBM also showed a longer and narrower shape. However, it shall be noted that the white line segments in the aerial photo (Figure 7b) were not the elevation contour lines. The white line segments represent the horizontal displacement of the gypsum tailings. They also indicate that the velocities along the central streamline were faster than those in the other regions. An indirect validation can be seen in the free-surface velocity profile, shown in Figure 9 at t = 30 s. However, the free-surface velocity will gradually reduce to zero as the freezing time is approaching.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 10 of 22 simulation of the deposited tailings from MBM. The flood distance at the freezing time t = 110 s was 310 m, which was 97% accurate to the filed data [60]. The result from MBM also showed a longer and narrower shape. However, it shall be noted that the white line segments in the aerial photo ( Figure  7b) were not the elevation contour lines. The white line segments represent the horizontal displacement of the gypsum tailings. They also indicate that the velocities along the central streamline were faster than those in the other regions. An indirect validation can be seen in the freesurface velocity profile, shown in Figure 9 at t = 30 s. However, the free-surface velocity will gradually reduce to zero as the freezing time is approaching.   9 show the time evolution of the free-surface velocity of the mudflows from t = 0~110 s by using BM and MBM, respectively. The velocity at the early stage (about t = 0~10 s) of BM was higher than that of MBM. However, the simulated mudflow in BM stopped earlier (around 70-90 s). The surface velocity gradually approached zero from t = 70 s, and the flow came to a full stop at t = 90 s. The gypsum tailings distance was around 220 m, and the mean velocity was around 2.4-3.1 m s −1 as shown in Figure 8. On the other hand, the flood distance from MBM was 310 m which was much closer to the field observation of 300 m [60]. The larger un-yield viscosity and yield strain rate limited the flood velocity at the early stage of the event (Figure 9). The mud started to liquefy and collapse in a small region near the breach in the first 10 s. The spreading shape of the tailings was symmetric along the centerline of the breach during t = 0-20 s. Because the supply of the tailings from the impoundment was asymmetric, the spreading shape gradually became asymmetric when t > 20 s. The maximum velocity of the released tailings occurred at t = 30 s. Then, the flow velocity gradually decreased. After t = 90 s, the tailings slowed down and stopped moving at t = 110 s. The mean velocity of the tailings flow was estimated at around 2.8-3.4 m s −1 . Appl. Sci. 2020, 10, x FOR PEER REVIEW 10 of 22 simulation of the deposited tailings from MBM. The flood distance at the freezing time t = 110 s was 310 m, which was 97% accurate to the filed data [60]. The result from MBM also showed a longer and narrower shape. However, it shall be noted that the white line segments in the aerial photo ( Figure  7b) were not the elevation contour lines. The white line segments represent the horizontal displacement of the gypsum tailings. They also indicate that the velocities along the central streamline were faster than those in the other regions. An indirect validation can be seen in the freesurface velocity profile, shown in Figure 9 at t = 30 s. However, the free-surface velocity will gradually reduce to zero as the freezing time is approaching.   9 show the time evolution of the free-surface velocity of the mudflows from t = 0~110 s by using BM and MBM, respectively. The velocity at the early stage (about t = 0~10 s) of BM was higher than that of MBM. However, the simulated mudflow in BM stopped earlier (around 70-90 s). The surface velocity gradually approached zero from t = 70 s, and the flow came to a full stop at t = 90 s. The gypsum tailings distance was around 220 m, and the mean velocity was around 2.4-3.1 m s −1 as shown in Figure 8. On the other hand, the flood distance from MBM was 310 m which was much closer to the field observation of 300 m [60]. The larger un-yield viscosity and yield strain rate limited the flood velocity at the early stage of the event (Figure 9). The mud started to liquefy and collapse in a small region near the breach in the first 10 s. The spreading shape of the tailings was symmetric along the centerline of the breach during t = 0-20 s. Because the supply of the tailings from the impoundment was asymmetric, the spreading shape gradually became asymmetric when t > 20 s. The maximum velocity of the released tailings occurred at t = 30 s. Then, the flow velocity gradually decreased. After t = 90 s, the tailings slowed down and stopped moving at t = 110 s. The mean velocity of the tailings flow was estimated at around 2.8-3.4 m s −1 .  Figures 8 and 9 show the time evolution of the free-surface velocity of the mudflows from t = 0~110 s by using BM and MBM, respectively. The velocity at the early stage (about t = 0~10 s) of BM was higher than that of MBM. However, the simulated mudflow in BM stopped earlier (around 70-90 s). The surface velocity gradually approached zero from t = 70 s, and the flow came to a full stop at t = 90 s. The gypsum tailings distance was around 220 m, and the mean velocity was around 2.4-3.1 m s −1 as shown in Figure 8. On the other hand, the flood distance from MBM was 310 m which was much closer to the field observation of 300 m [60]. The larger un-yield viscosity and yield strain rate limited the flood velocity at the early stage of the event (Figure 9). The mud started to liquefy and collapse in a small region near the breach in the first 10 s. The spreading shape of the tailings was symmetric along the centerline of the breach during t = 0-20 s. Because the supply of the tailings from the impoundment was asymmetric, the spreading shape gradually became asymmetric when t > 20 s. The maximum velocity of the released tailings occurred at t = 30 s. Then, the flow velocity gradually decreased. After t = 90 s, the tailings slowed down and stopped moving at t = 110 s. The mean velocity of the tailings flow was estimated at around 2.8-3.4 m s −1 . Appl. Sci. 2020, 10, x FOR PEER REVIEW 11 of 22  The observation and simulation results of inundation distance, freezing time, and mean velocity provided in this study, as well as some historical studies, are listed in Table 1 for comparisons. The  The observation and simulation results of inundation distance, freezing time, and mean velocity provided in this study, as well as some historical studies, are listed in Table 1 for comparisons. The The observation and simulation results of inundation distance, freezing time, and mean velocity provided in this study, as well as some historical studies, are listed in Table 1 for comparisons. The simulation results of the MBM were not only accurate in the freezing time and the mean velocity, but also in the inundation distance.  [20] 360 120 3.0 Bingham model (Figure 8) 220 70-90 2.4-3.1 Modified bi-viscosity model (Figure 9) 310 90-110 2.8-3.4 Figure 10 shows the strain rate and the un-yield/yield zones of MBM results in the center cross-section of the breach (y = 220 m). The vertical axis was ten times exaggerated. The yield strain rate . γ y = 0.2 s −1 was chosen to identify the interface between the un-yield/yield zones. The liquefied zone was located at the breach's front when t = 10 s, due to the large strain rate caused by gravitational force. The liquefied zone gradually shifted to the ground region after t = 10 s, resulting from the large shear stress that occurs near the bottom. However, the neighboring areas remained un-yield. The interface between the un-yield and yield zone was presented during the period of 10-40 s. It was caused by the viscosity's discontinuity of un-yield and yield zones. At t = 90-110 s, the liquefied zone shrunk gradually and disappeared at t = 110 s, due to the zero velocity in the entire flow field. This zero-velocity phenomenon was very close to the real landslide situation in which the velocity ceases to zero eventually.
Appl. Sci. 2020, 10, x FOR PEER REVIEW 12 of 22 simulation results of the MBM were not only accurate in the freezing time and the mean velocity, but also in the inundation distance.  [20] 360 120 3.0 Bingham model (Figure 8) 220 70-90 2.4-3.1 Modified bi-viscosity model ( Figure 9) 310 90-110 2.8-3.4 Figure 10 shows the strain rate and the un-yield/yield zones of MBM results in the center crosssection of the breach (y = 220 m). The vertical axis was ten times exaggerated. The yield strain rate ̇= 0.2 s −1 was chosen to identify the interface between the un-yield/yield zones. The liquefied zone was located at the breach's front when t = 10 s, due to the large strain rate caused by gravitational force. The liquefied zone gradually shifted to the ground region after t = 10 s, resulting from the large shear stress that occurs near the bottom. However, the neighboring areas remained un-yield. The interface between the un-yield and yield zone was presented during the period of 10-40 s. It was caused by the viscosity's discontinuity of un-yield and yield zones. At t = 90-110 s, the liquefied zone shrunk gradually and disappeared at t = 110 s, due to the zero velocity in the entire flow field. This zero-velocity phenomenon was very close to the real landslide situation in which the velocity ceases to zero eventually.  Figure 11 shows the velocity magnitude profiles obtained from BM and MBM in the centerline cross-section (y = 220 m). In the BM results, the tailings moved faster than that in MBM. The maximum velocity of the tailings' front at t = 10 s was approximately 6.0-8.0 m s −1 and decreased sharply during t = 10-40 s. It made the inundation distance (around 220 m) at t = 110 s shorter than that in MBM. In the MBM results, the maximum moving velocity was about 5.0-7.0 m s −1 , which was slightly smaller than that in the results of BM. However, the MBM tailings flow took a longer time to reach the zero-velocity stage and the resulting inundation distance was longer than that in the BM results. The inundation distance (310 m) predicted by MBM was very close to the field observed (300 m). The yield strain rate played an important role in MBM. As long as the strain rate was smaller than the yield strain rate, a large viscosity was applied to slowing down the deformation. Due to the absence of measuring techniques, the yield strain rate in this study was obtained from the sensitivity analysis by taking the flood distance from the field observation as the criterion.

Difference between the BM and MBM
Appl. Sci. 2020, 10, x FOR PEER REVIEW 13 of 22 Figure 10. Snapshots on the centerline cross-section of the breach (y = 220 m) predicted by MBM. (a) The profiles of strain rate. The color bar is set from 0.0 to 1.0 s −1 to emphasize the interface between the plug zone and liquefied zone at ̇= 0.2 s −1 . (b) The discontinuity of the plug zone and liquefied zone. Figure 11 shows the velocity magnitude profiles obtained from BM and MBM in the centerline cross-section (y = 220 m). In the BM results, the tailings moved faster than that in MBM. The maximum velocity of the tailings' front at t = 10 s was approximately 6.0-8.0 m s −1 and decreased sharply during t = 10-40 s. It made the inundation distance (around 220 m) at t = 110 s shorter than that in MBM. In the MBM results, the maximum moving velocity was about 5.0-7.0 m s −1 , which was slightly smaller than that in the results of BM. However, the MBM tailings flow took a longer time to reach the zerovelocity stage and the resulting inundation distance was longer than that in the BM results. The inundation distance (310 m) predicted by MBM was very close to the field observed (300 m). The yield strain rate played an important role in MBM. As long as the strain rate was smaller than the yield strain rate, a large viscosity was applied to slowing down the deformation. Due to the absence of measuring techniques, the yield strain rate in this study was obtained from the sensitivity analysis by taking the flood distance from the field observation as the criterion.   (Figure 12a). On the other hand, in MBM results, the yield strain rate ̇= 0.2 s −1 was introduced as the indicator to identify the plug and sheared zone. Because the un-yield viscosity = 10 10 Pa s was much greater than /̇, a discontinuous pattern of the strain rate could be observed in Figure 12b. The yield strain rate ̇= 0.2 s −1 kept the plug zone rigid. The  Figure 12 illustrates the strain rate profile of the initiation process of the tailings flow. The strain rate profiles in BM results showed a smooth and continuous feature. A large amount of tailing material deformed and slid down (Figure 12a). On the other hand, in MBM results, the yield strain rate . γ y = 0.2 s −1 was introduced as the indicator to identify the plug and sheared zone. Because the un-yield viscosity µ A = 10 10 Pa s was much greater than τ y / . γ y , a discontinuous pattern of the strain rate could be observed in Figure 12b. The yield strain rate . γ y = 0.2 s −1 kept the plug zone rigid. The initiation process of the mudslide in MBM results was different from that in BM results. A high strain rate appeared not only near the toe of the breach but also in the gate area, which caused the sliding process and formed a slip surface. The slip surface was the interface between the un-yield and yield regions. In the bank of homogeneous mud, the slip surface of failure could be determined by the empirical method, which follows the arc of a circle that usually intersects the toe of the bank [52,53]. However, the slip surface was developed automatically by MBM. It is worth a more profound study in the future.

Difference between the BM and MBM
Appl. Sci. 2020, 10, x FOR PEER REVIEW 14 of 22 initiation process of the mudslide in MBM results was different from that in BM results. A high strain rate appeared not only near the toe of the breach but also in the gate area, which caused the sliding process and formed a slip surface. The slip surface was the interface between the un-yield and yield regions. In the bank of homogeneous mud, the slip surface of failure could be determined by the empirical method, which follows the arc of a circle that usually intersects the toe of the bank [52,53]. However, the slip surface was developed automatically by MBM. It is worth a more profound study in the future.    Because of adopting the larger un-yield viscosity µ A , larger yield strain rate . γ y , same yield stress τ y , and same yield viscosity µ B to simulate the FGT66, the un-yield zone of tailings in MBM was sturdier than that in the BM. However, once tailings yielded, the rheology returned to the conventional Bingham properties. Because of adopting the larger un-yield viscosity , larger yield strain rate ̇, same yield stress , and same yield viscosity to simulate the FGT66, the un-yield zone of tailings in MBM was sturdier than that in the BM. However, once tailings yielded, the rheology returned to the conventional Bingham properties.

The Role of the Grid Resolution
The grid resolution was a key factor in the strain rate calculation. To understand the sensitivity of the grid resolution of the flood distance, three cases of BM and nine cases of MBM were performed with the resolutions varying from dx: 1.8-2.2, dy: 1.8-2.2, and dz: 0.8-1.2 m. Figure 14 shows the deposited boundary of the tailings predicted by BM and MBM. The results show that BM was less sensitive to the resolution than MBM. The results from MBM show that the deposited boundary was more sensitive with dz than dx and dy. The inundation distance was shorter with dz = 1.2 m (blue lines) and longer with dz = 0.8 m (red lines). To give an overview, the inundation width and inundation distance were convergent as dz < 1.2 m.

The Role of the Grid Resolution
The grid resolution was a key factor in the strain rate calculation. To understand the sensitivity of the grid resolution of the flood distance, three cases of BM and nine cases of MBM were performed with the resolutions varying from dx: 1.8-2.2, dy: 1.8-2.2, and dz: 0.8-1.2 m. Figure 14 shows the deposited boundary of the tailings predicted by BM and MBM. The results show that BM was less sensitive to the resolution than MBM. The results from MBM show that the deposited boundary was more sensitive with dz than dx and dy. The inundation distance was shorter with dz = 1.2 m (blue lines) and longer with dz = 0.8 m (red lines). To give an overview, the inundation width and inundation distance were convergent as dz < 1.

The Role of the Yield Strain Rate
In MBM, the yield strain rate defined the fluid behavior in the regime of un-yield and yield zones. If the yield strain rate was zero, the material returned to Bingham fluid. As the yield strain rate became higher, the plug zone became wider. The fluid was harder to transfer from un-yield to yield phases. The behavior of the tailings was not monotone when increased. There were two kinds of behavior of the tailings, presented in Figure 15. In the regime of 0.0 0.3 s −1 , the result

The Role of the Yield Strain Rate
In MBM, the yield strain rate . γ y defined the fluid behavior in the regime of un-yield and yield zones. If the yield strain rate was zero, the material returned to Bingham fluid. As the yield strain rate became higher, the plug zone became wider. The fluid was harder to transfer from un-yield to yield phases. The behavior of the tailings was not monotone when . γ y increased. There were two kinds of behavior of the tailings, presented in Figure 15. In the regime of 0.0 ≤ . γ y ≤ 0.3 s −1 , the result shows that when . γ y increased, the inundation widths were narrower and the inundation lengths were longer.
In the regime of 0.3 ≤ . γ y ≤ 0.6 s −1 , when . γ y increased, not only were the inundation widths narrower, the inundation lengths were also shorter as well.

The Role of the Yield Strain Rate
In MBM, the yield strain rate ̇ defined the fluid behavior in the regime of un-yield and yield zones. If the yield strain rate was zero, the material returned to Bingham fluid. As the yield strain rate became higher, the plug zone became wider. The fluid was harder to transfer from un-yield to yield phases. The behavior of the tailings was not monotone when ̇ increased. There were two kinds of behavior of the tailings, presented in Figure 15. In the regime of 0.0 ≤̇≤ 0.3 s −1 , the result shows that when ̇ increased, the inundation widths were narrower and the inundation lengths were longer. In the regime of 0.3 ≤̇≤ 0.6 s −1 , when ̇ increased, not only were the inundation widths narrower, the inundation lengths were also shorter as well.  Rheological properties of hyper-concentration are generally formulated as a function of the fluid concentration. Julien [44] recommended empirical formulas with the exponential relationships for yield stress and viscosity at large concentrations of fines. The typical values of coefficients for different types of mud, clay, and lahar are presented in Table 2. Kaolinite and Typical soils are utilized to describe the features of BM and MBM in this section. Eighteen numerical cases, including nine Bingham cases and nine modified Bi-viscosity cases with different concentration C v , are performed. The yield shear stress and the viscosity are presented in Table 3. The yield strain rate is specified as . γ y = 0.0 s −1 for the Bingham cases and . γ y = 0.2 s −1 for the modified Bi-viscosity cases. The goal of this analysis is to find a material, which has a similar property to the tailings in FGT66. A similar inundation profile is the key. Figure 16 shows the simulation results with different concentration C v by BM and MBM. The deposited boundary of Kaolinite at C v = 0.5 (red line in Figure 16b) has the best fit to the MBM result of FGT66 (red line in Figure 15). It is due to the similarity of τ y between two materials: Kaolinite τ y = 1580 Pa and the gypsum tailings τ y = 1500 Pa.

Conclusions and Future Work
This study applied the rheology models, BM, CBM, and MBM, to the mudslide simulation. The rheology models were integrated into the 3D Navier-Stokes equations coupled with the LES turbulent model to give a detailed description of the vertical acceleration. The free-surface kinematic was described by the VOF method with a PLIC surface tracking scheme. The BM results were validated by the channel flow and the spreading slow data, which received good agreements in both cases.
The Bingham model (BM) and conventional bi-viscosity model (CBM) were then used to simulate FGT66. The predicted inundation distance was 220 m with accuracy at about 73% of the observed data. To improve the result, a modification to CBM was raised. A large viscosity number of un-yield region, = 10 Pa s, was applied to representing a more rigid behavior of the material as suggested by Assier Rzadkiewicz [34], Taibi [45], and Yu [28]. A series of sensitivity analyses on the yield strain rate was performed by matching the simulated tailings' boundary to the field observation. A larger yield strain rate indicated the tailings with sturdier behavior. The yield strain rate was suggested to be = 0.2 s −1 in the MBM to simulate the FGT66. By the sensitivity analysis, the material of FGT66 was close to the kaolinite with concentration = 0.5. The results show that the modified bi-viscosity model (MBM) could provide a better prediction than BM and CBM in terms

Conclusions and Future Work
This study applied the rheology models, BM, CBM, and MBM, to the mudslide simulation. The rheology models were integrated into the 3D Navier-Stokes equations coupled with the LES turbulent model to give a detailed description of the vertical acceleration. The free-surface kinematic was described by the VOF method with a PLIC surface tracking scheme. The BM results were validated by the channel flow and the spreading slow data, which received good agreements in both cases.
The Bingham model (BM) and conventional bi-viscosity model (CBM) were then used to simulate FGT66. The predicted inundation distance was 220 m with accuracy at about 73% of the observed data. To improve the result, a modification to CBM was raised. A large viscosity number of un-yield region, µ A = 10 10 Pa s, was applied to representing a more rigid behavior of the material as suggested by Assier Rzadkiewicz [34], Taibi [45], and Yu [28]. A series of sensitivity analyses on the yield strain rate . γ y was performed by matching the simulated tailings' boundary to the field observation. A larger yield strain rate indicated the tailings with sturdier behavior. The yield strain rate was suggested to be . γ y = 0.2 s −1 in the MBM to simulate the FGT66. By the sensitivity analysis, the material of FGT66 was close to the kaolinite with concentration C v = 0.5. The results show that the modified bi-viscosity model (MBM) could provide a better prediction than BM and CBM in terms of the flood distance and the spreading width. The development of the flood free-surface, velocity, strain rate, and un-yield/yield zone were presented and discussed.
The important results and conclusions are listed as follows: • The Bingham model was successfully integrated into the 3D CFD model to simulate a mudslide by adopting the implicit scheme for the viscosity term; • If the yield strain was small enough, the conventional bi-viscosity model would converge into the Bingham model numerically;

•
The flood distance predicted by BM and CBM was 220 m, with the accuracy at about 73% of the field observation; • To improve the result, MBM was introduced by giving a large viscosity to the un-yield phase and adjusting the yield strain rate with a sensitivity analysis; • The flood distance predicted by MBM was 310 m, which was closer to the filed observation with 97% of the accuracy; • Not only did the flood distance improve, the spreading width improved in the result of MBM; • The free-surface, velocity, strain rate, and the un-yield/yield profiles were presented and discussed; • The sensitivity analysis of the grid resolution was performed. The grid resolution in the BM was less sensitive than that in the MBM; • The analysis of the yield strain rate was carried out. A larger yield strain rate indicated the tailings with sturdier behavior. The yield strain rate was suggested to be . γ y = 0.2 s −1 in the MBM in the case of FGT66; • The un-yield zone in MBM was sturdier than that in the pure BM. However, once tailings material yielded, the rheology returned to the conventional Bingham liquefied properties; • By the sensitivity analysis, the material of FGT66 was close to the kaolinite with concentration C v = 0.5.
This work is the Part I of modeling the slump-type landslide tsunamis. The Bingham-type rheology model was developed to simulate slumps and tailings flows. Model validation and sensitivity analysis were performed. Tsunamis excited by slumps will be studied in Part II.