Interactive comment on “Model development for simulating mudslide and the case study of the failure of the gypsum tailings dam in East Texas in 1966”

the we can’t reach the conclusion of stratiﬁcation. we conclude that the un-yield region of tailing material is sturdier based on the result with a larger un-yield viscosity and yield strain rate. been largely modiﬁed to satisfy the standard of NHESS. The abstract has been revised as below: This paper adopts Bingham and Bi-viscosity rheology models to simulate the dynamic and kinematic processes of mudslides. The rheology models are integrated into a computational ﬂuid dynamics code, Splash3D, to solve the incompressible Navier-Stokes equations with Volume of Fluid surface tracking algorithm. The changes of the un-yield and yield phases of the mudslide material are 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 study the failure of the gypsum tailings dam in East Texas in 1966. The results show that 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 tailing material. The yield strain rate is obtained by comparing the simulated ﬂood distance to the ﬁeld data. This modiﬁed Bi-viscosity model improves not only the accuracy of the ﬂood distance to about 97% but also the accuracy of the spread width. This result suggests that the un-yield region in the modiﬁed Bi-viscosity model was sturdier than that described in the Bingham model. However, once the tailing material yields, the tailing ﬂuid returns to the conventional Bingham liqueﬁed material.

these aspects clear. Answer: This paragraph has been modified to be more clear as below: The rheological properties are crucial issues in a tailing flow simulation. The tailings fluid formed by mixing tailings and water is a non-Newtonian fluid in nature (Henriquez and Simms, 2009) with complex rheological properties. The travel distance and the spreading of a tailings flow are affected by the rheological equation (Yu et al., 2020). As for concerning the determination of the rheological parameters, Henriquez and Simms, (2009) determined the yield stress and viscosity of tailings flow using rheometer and slump tests. The mixture of different materials leads to a complex, yet not well understood rheological behavior (Von Boetticher et al., 2017). Field observations of mudflow behavior and rheology are challenging and still rare. Numerical modeling is chosen when an assessment of mudflow behavior is needed for planning, zoning, and hazard assessment (Scheuner et al., 2011;Christen et al., 2012;Kattel et al., 2016;Mergili et al., 2017). Most models require direct calibration to capture site-specific behavior. However, reliable calibration data are scarce, and laboratory experiments are difficult to upscale to field situations (Von Boetticher et al., 2017).
Answer: These sentences have been modified as below: Mudflow can be simulated by Bingham model (Schamber and MacArthur, 1985, Liu and Mei, 1989, Herschel-Bulkley model (Bates and Ancey, 2017;Huang and García, 1998), and the Coulomb-viscoplastic model (Pudasaini and Mergili, 2019) by incorporating them with the depth-integrated equations models. Referee #1: L51: "with strong vertical fluid particle acceleration." -> "with strong vertical fluid particle acceleration (Domnik et al., 2013: https://doi.org/10.1016/j.jnnfm.2012.03.001, https://doi.org/10.1016/j.jnnfm.2013.07.005)." Answer: These sentences have been modified as below: The depth-integrated model is simplified from Navier-Stokes C3 equations by ignoring the vertical acceleration (Han et al., 2019). The depth-integrated model is suitable for predicting the flow without strong vertical acceleration or a sharp velocity shearing (Khattri and Pudasaini, 2019). The vertical acceleration and velocity shearing are important in the case of a slope with rugged topography or mudslide overtopping a structure. Inside a complex 3D flow structure, the tailings material will transfer from an un-yield/plug zone to a yield/liquefied/sheared zone if the shear stress is greater than the yield shear stress. Before the tailings material reaching this plug zone, the liquefied zone might dominate the entire flow field due to strong shear. Adopting a three-dimensional rheology model is an alternative option to fully and globally describe the strongly converging and diverging flows (Domnik et al., 2013;Yu et al., 2020). For more detailed results, the three-dimensional Navier-Stokes equations are recommended (Lee et al., 2010;. The Navier-Stokes equations modes were utilized to study mudflow in the 1990s. Many of these studies concentrated on solving 2D problems (O'Brien et al., 1993;Assier Rzadkiewicz et al., 1997;Huang and García, 1998). With the advance of computers in the early twentieth century, the 3D Navier-Stokes equations were able to study mudflow by Smoothed Particle Hydrodynamics (SPH) method (Dai et al., 2014;, and projection method (VonBoetticher et al., 2017;Abadie et al., 2019;Yu et al., 2020).
Referee #1: L55-60: Please improve writing. Answer: The writing of the entire article has been greatly improved.
Referee #1: L74-82: Very basic, could be removed. Answer: Agree. They are removed. Referee #1: L82: "However, Navier-Stokes equations derived under the assumption C4 of Eulerian fails in describing the un-yield motion.": What does it mean? Not clear. Answer: This sentence is not clear. It is removed.
Referee #1: L86-89: Please explain how you would obtain stratification, because this is the main aspect of this MS. And, it seems to be mostly parameter fit, increase mu_A as high as you get what you want. So, what is the mechanical aspect of this paper, and how can you justify using those mu_A values. This is a critical aspect here. Please clearly justify, discuss. Answer: Thanks for the comments. After a long discussion with the co-authors, we agree that we can't reach the conclusion of stratification. Instead, we can conclude that the un-yield region of tailing material is sturdier based on the result with a larger un-yield viscosity and yield strain rate. The mechanical aspect of this paper has been clarified as below: The rheological properties of BM can be presented as (Tanner, 1982;Byron-Bird et al., 1983): µ(γ ÌG )={âŰȃ(µ_A=∞ and γ ÌG=0, if τ < τ _y@µ_B+τ _y/ √ (1/2 γ ÌG_ij:γ ÌG_ij ) and γ ÌG>0, if τ ≥τ _y )âŤd' (1) where µ_A is the viscosity of the un-yield region, µ_B is the viscosity of the yield zone, τ _y is the yield stress, and γ ÌG_y is the yield strain rate, γ ÌG_ij=(∂u ÌG_i)/(∂x_j )+(∂u ÌG_j)/(∂x_i ). The symbol γ ÌG is the second invariant of the γ ÌG_ij, which is defined as γ ÌG= √ (1/2 γ ÌG_ij γ ÌG_ij ). Despite the simplicity of the BM, the fact that the stress is indeterminate in the un-yield region, which means the exact shape and location of the yield surface(s) cannot be determined (Khabazi et al., 2016). To remedy this drawback, in the present work, the conventional Bi-viscosity model (CBM) proposed by Beverly and Tanner (1992) is 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 (Khabazi et al., 2016). The rheological properties of CBM can be presented as (Beverly and Tanner, 1992;Khabazi et al., 2016): µ(γ ÌG )={âŰȃ(µ_A=τ _y/γ ÌG_y ,if γ ÌG<γ ÌG_y@µ_B+τ _y/ √ (1/2 γ ÌG_ij:γ ÌG_ij ) ,if γ ÌG≥γ ÌG_y )âŤd' (2) Mathematically speaking, when γ ÌG_y approaches to zero, the CBM will approach to BM. By choosing γ ÌG_y sufficiently small, we can practically replace the C5 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(s) (Khabazi et al., 2016). However, if a mud material experiences compaction or tamping processes, the material will become sturdy, and resulting in a larger µ_A. To describe the sturdy behavior in the plug zone, a larger µ_A and a larger γ ÌG_y are required. The larger µ_A, playing a role of keeping the rigid shape, and the larger γ ÌG_y, indicating the material can sustain a larger deformation. This modified Bi-viscosity model ( , (2020). In this paper, the infinite number of viscosity µ_A=ãȂŰ10ãȂŮˆ10 Pa s is chosen by a sensitivity analysis. The values of yield strain rate γ ÌG_y are also discussed in Section 5.2. By sensitivity analysis, γ ÌG_y=0.2 s-1 is adopted to illustrate the deformation in MBM.
Does (3) avoid discontinuity, or increase it? Consideration of (3) does not seem to constitute a new model as mentioned in the Title. This appears to be defined for simulation purpose to try to obtain pattern seen in the field event. If not, please provide arguments. Answer: The Eq.
(3) is used to introduce the discontinuity. Agree. Eq. (3) is not a new model. So, we have modified the title as: "Numerical study on the failure of the gypsum tailings dam in East Texas in 1966".
Referee #1: L108: It is a bit confusing. Sometime you talk about un-yield sometime about fluid, and now it appears that the material is dry. What type of material, and composition is considered? Further, how do you determine r? As the material density changes during motion, do you have a transport equation for r? The VOF method was also used by Yu et al. (2020) to track the interface between the tailing fluid and air. Also, improve English in the following lines. Answer: Thanks for the comments. This part has been modified as below: 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-immersible fluids. The free-surface between the tailings fluid and air is tracked by the Volume of Fluid (VOF) method (Hirt and Nichols, 1981).  (6) and (7), or only (8). (ii) If (6) and (7) are used what is the advantage? (iii) Can you provide some references? Or, do you develop all these equations (3)-(11)? The readers would wonder. (iv) Can you define what is un-yield and fluid? It is completely confusing. You mention sediment and fluid, but never mention what material/composition you are considering. Are you using un-yield + fluid two-phase mixture material? Does not seem so. Answer: (i)∼(iii) This part has been removed out of the manuscript because we have just published a new paper (Chu et al., 2020), which described the details of the projection method. We also present the detailed numerical algorithm here: The projection method (DeLong, 1997) was used to decouple the velocity and pressure in the Navier-Stokes equations and to solve the Poisson pressure equation (PPE). The Navier-Stokes equations were solved by the two-step projection method. The first step was the predictor step: (ˆ(n+1) uˆ*-ˆn uˆn)/∆t=-âĹG·(uu)ˆn-âĹGPˆn+âĹG·(µˆ(n+1) (âĹGu+âĹGˆT u)ˆn )+ˆn gˆn (1) where the superscripts n and n+1 represent the old and new time step, respectively; superscript * represents an intermediate time between time n and n+1. The density of the cellˆ(n+1) in the cell was calculated by the volume-weight method:ˆ(n+1)=ˆ(n+1) _t (2) where _t is the density of the tailing fluid. Different from the conventional two-step projection method, the old-time step pressure was included in the predictor step to increase the accuracy of uˆ*. Because uˆ* was calculated from the old-time solution explicitly, the fluid velocity uˆ* does not satisfy the Navier-Stokes equations and needs to be corrected. The corrector step is: (ˆ(n+1) uˆ(n+1)-ˆ(n+1) uˆ*)/∆t=-âĹGãȂŰδPãȂŮˆ(n+1)+ˆ(n+1) gˆ(n+1)-ˆn gˆn (3) where the change in pressure was calculated as ãȂŰδPãȂŮˆ(n+1)=Pˆ(n+1)-Pˆn. Taking divergence to equation (6) yields the Pressure Poisson Equation (PPE) for solving Pˆ(n+1): âĹG·(âĹGδPˆ(n+1))/ˆ(n+1) =âĹG·(uˆ*/∆t+F_bˆ(n+1)-F_bˆn ) (4) In this corrector step, the new velocity uˆ(n+1) can be solved when the pressure Pˆ(n+1) is known. In this study, the tailing fluid is assumed to be Non-Newtonian fluids. Viscous forces are incorporated into the predictor step to estimate the cell-centered velocity field by using the * time step velocity field. This implicit approximation eliminates the limitation of the C8 time step for the stability criteria.
Answer: (iv) Sorry if it is confusing you. In the present model, the tailing fluid has been treated as a single homogeneous mud material. The rheology of tailing fluid is defined into two zones: un-yield/plug zone and yield/liquefied/sheared zone. The tailings material will transfer from an un-yield/plug zone to a yield/liquefied/sheared zone if the shear stress is greater than the yield shear stress, or vice versa.
Referee #1: L138: "a large viscosity parameter is imposed": again, no physical explanation. Answer: The physical explanation is added as below: The rheological properties of BM can be presented as (Tanner, 1982;Byron-Bird et al., 1983): µ(γ ÌG )={âŰȃ(µ_A=∞ and γ ÌG=0, if τ < τ _y@µ_B+τ _y/ √ (1/2 γ ÌG_ij:γ ÌG_ij ) and γ ÌG>0, if τ ≥τ _y )âŤd' (1) where µ_A is the viscosity of the un-yield region, µ_B is the viscosity of the yield zone, τ _y is the yield stress, and γ ÌG_y is the yield strain rate, γ ÌG_ij=(∂u ÌG_i)/(∂x_j )+(∂u ÌG_j)/(∂x_i ). The symbol γ ÌG is the second invariant of the γ ÌG_ij, which is defined as γ ÌG= √ (1/2 γ ÌG_ij γ ÌG_ij ). Despite the simplicity of the BM, the fact that the stress is indeterminate in the un-yield region, which means the exact shape and location of the yield surface(s) cannot be determined (Khabazi et al., 2016). To remedy this drawback, in the present work, the conventional Bi-viscosity model (CBM) proposed by Beverly and Tanner (1992) is 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 (Khabazi et al., 2016). The rheological properties of CBM can be presented as (Beverly and Tanner, 1992;Khabazi et al., 2016): µ(γ ÌG )={âŰȃ(µ_A=τ _y/γ ÌG_y ,if γ ÌG<γ ÌG_y@µ_B+τ _y/ √ (1/2 γ ÌG_ij:γ ÌG_ij ) ,if γ ÌG≥γ ÌG_y )âŤd' (2) Mathematically speaking, when γ ÌG_y approaches to zero, the CBM will approach to BM. By choosing γ ÌG_y 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 C9 the shape and the location of the plug(s) (Khabazi et al., 2016). However, if a mud material experiences compaction or tamping processes, the material will become sturdy, and resulting in a larger µ_A. To describe the sturdy behavior in the plug zone, a larger µ_A and a larger γ ÌG_y are required. The larger µ_A, playing a role of keeping the rigid shape, and the larger γ ÌG_y, indicating the material can sustain a larger deformation. This modified Bi-viscosity model ( , (2020). In this paper, the infinite number of viscosity µ_A=ãȂŰ10ãȂŮˆ10 Pa s is chosen by a sensitivity analysis. The values of yield strain rate γ ÌG_y are also discussed in Section 5.2. By sensitivity analysis, γ ÌG_y=0.2 s-1 is adopted to illustrate the deformation in MBM.
Referee #1: L158-159: Could be removed, not necessary, very basic. Table 1 is not needed, all parameters are in Fig. Answer: Thank you. Those parts have been removed.
Referee #1: L176: "Because of the yield stress, the free surface parallels to the plane bed": But not in the front. L177: "The mud front, just like a steady gravity current, eventually advances at a constant speed": Can you show it analytically, or you are talking only about numerical results? Answer: Sorry for the confusion. This part has been improved as: 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 (Liu and Mei, 1989). The numerical results present a similar pattern of analytical solutions to that in Liu and Mei, (1989).

C10
Referee #1: L190: "lateral boundaries (y = 0 m and y = 400 m) are free-slip boundary": Strange! Why don't you explain why this condition is used while the base is no-slip? Answer: Thank you for the comments. This part is improved as: The downstream (x = 400 m) and lateral boundaries (y = 0 m and y = 400 m) are free-slip walls. The downstream and lateral boundaries will not affect the simulation results because the domain is set to be much larger than the predicted tailing pattern.
Referee #1: L193: The value of mu_A is huge! Can you justify it mechanically, or prove with evidence as the paper is mainly based on such values? Answer: The mechanism and the reference have been added into the manuscript as below: The rheological properties of BM can be presented as (Tanner, 1982;Byron-Bird et al., 1983): µ(γ ÌG )={âŰȃ(µ_A=∞ and γ ÌG=0, if τ < τ _y@µ_B+τ _y/ √ (1/2 γ ÌG_ij:γ ÌG_ij ) and γ ÌG>0, if τ ≥τ _y )âŤd' (1) where µ_A is the viscosity of the un-yield region, µ_B is the viscosity of the yield zone, τ _y is the yield stress, and γ ÌG_y is the yield strain rate, γ ÌG_ij=(∂u ÌG_i)/(∂x_j )+(∂u ÌG_j)/(∂x_i ). The symbol γ ÌG is the second invariant of the γ ÌG_ij, which is defined as γ ÌG= √ (1/2 γ ÌG_ij γ ÌG_ij ). Despite the simplicity of the BM, the fact that the stress is indeterminate in the un-yield region, which means the exact shape and location of the yield surface(s) cannot be determined (Khabazi et al., 2016). To remedy this drawback, in the present work, the conventional Bi-viscosity model (CBM) proposed by Beverly and Tanner (1992) is 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 (Khabazi et al., 2016). The rheological properties of CBM can be presented as (Beverly and Tanner, 1992;Khabazi et al., 2016): µ(γ ÌG )={âŰȃ(µ_A=τ _y/γ ÌG_y ,if γ ÌG<γ ÌG_y@µ_B+τ _y/ √ (1/2 γ ÌG_ij:γ ÌG_ij ) ,if γ ÌG≥γ ÌG_y )âŤd' (2) Mathematically speaking, when γ ÌG_y approaches to zero, the CBM will approach to BM. By choosing γ ÌG_y 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 C11 the shape and the location of the plug(s) (Khabazi et al., 2016). However, if a mud material experiences compaction or tamping processes, the material will become sturdy, and resulting in a larger µ_A. To describe the sturdy behavior in the plug zone, a larger µ_A and a larger γ ÌG_y are required. The larger µ_A, playing a role of keeping the rigid shape, and the larger γ ÌG_y, indicating the material can sustain a larger deformation. This modified Bi-viscosity model ( , (2020). In this paper, the infinite number of viscosity µ_A=ãȂŰ10ãȂŮˆ10 Pa s is chosen by a sensitivity analysis. The values of yield strain rate γ ÌG_y are also discussed in Section 5.2. By sensitivity analysis, γ ÌG_y=0.2 s-1 is adopted to illustrate the deformation in MBM.
Referee #1: L196: "Splash3D": Please explain it briefly, not all the readers might know it. Answer: Thank you for your comment. The brief explanation of Splash3D is added into the manuscript as below: In this paper, the viscoplastic models, BM, CBM, and MBM are coupled with the Splash3D model. The Splash3D model is renovated from the open-source software, Truchas, which was originally developed by Los Alamos National Laboratory (Lam et al., 2007). 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 (Chorin, 1968(Chorin, , 1997Chu et al., 2020) and the finite volume discretization method (Eymard et al., 2000). The Splash3D model was enhanced with several hydrodynamic modules such as the Large-Eddy Simulation (LES) turbulence module (Wu, 2004;Liu, 2009), moving-solid module (Chu et al., 2020) to deal with the breaking wave, and wave-obstacle iteration problems. Readers are encouraged to read the reference Chu et al., (2020) for the detailed numerical algorithm. In this study, the Splash3D is developed with the rheological C12 model to solve the mudflow problems.
Referee #1: L211-215: The most important aspect of this paper is concerned with the generation of the stratification pattern observed in the field event. However, the authors could not make it quite clear what exactly is the mechanism in their model that produces these patterns. This is perhaps the major drawback of the MS. How can just gravity do that job? Please make it clear. Answer: Referred to the early reply, we can't reach the conclusion of stratification. Instead, we can conclude that the un-yield region of tailing material is sturdier based on the result with a larger un-yield viscosity and yield strain rate. Besides, several conclusions are modified in section 6. The flood distance predicted by MBM is 310 m which is closer to the filed observation with accuracy at about 97%. Not only the flood distance but also the spread width is improved in the MBM result. This result suggests that the un-yield region of tailing material in MBM results is sturdier than one in the pure Bingham material. However, once tailing material yields, the rheology returns to the Bingham properties.
Referee #1: L232-233, 250-252: "Not only freezing time, and mean velocity but also the value of inundation distance receives good agreements with observed values.", "the MBM result takes a longer time to reach the zero-velocity stage and the resulting inundation distance is longer than the BM result.": This is very good! However, not yet clear which mechanism produces the real pattern. And the pattern seen in the simulation is very weak as compared with the real field event. Please discuss. Answer: Thanks for the comments. The pattern (the white line segments) in the photo represents the horizontal displacement of the gypsum tailing. They also indicate that the velocities along the flood central streamline are faster than the other region. This part has been improved as : Fig 7 (a) shows the simulated deposited tailings from MBM. The flood distance at the freezing time t = 110 s is 310 m, which is much closer to the filed data (Jeyapalan et al., 1983). The result from MBM also shows a longer and narrower shape. However, it shall be noted that the white line segments in the aerial photo (Error! Reference source not found. (b)) are not the elevation contour lines. The white line C13 segments represent the horizontal displacement of the gypsum tailings. They also indicate that the velocities along the central streamline are faster than those in the other regions. An indirect validation can be seen in the free-surface velocity profile, shown in Error! Reference source not found. at t = 30 s. However, the free-surface velocity will gradually reduce to zero as the freezing time is approaching.
Referee #1: L254-255. "The sliding slope of the mudslide is developed automatically by MBM. This is big progress in studying mudslide flows.": No! This slope is also seen in BM. So, not clear what you want to say.
Answer: Sorry if it is confusing you. "The sliding slope" shall be "The slip surface". This part has been improved as: Error! Reference source not found. illustrates the strain rate profile of the initiation process of the tailing flow. The strain rate profiles in BM results show a smooth and continuous feature (Error! Reference source not found. (a)). A large amount of tailing material deforms and slides down (Error! Reference source not found. a)). On the other hands, in MBM results, the yield strain rate γ ÌG_y=0.2 s-1 is introduced as the indicator to identify the plug and sheared zone. Because the unyield viscosity µ_A=ãȂŰ10ãȂŮˆ10 Pa s is much greater than τ _y/γ ÌG_y, a discontinuity pattern of the strain rate can be observed in Error! Reference source not found. (b). The yield strain rate γ ÌG_y=0.2 s-1 keeps the plug zone rigid. The initiation process of mudslide in MBM results is different from the ones in BM results. A high strain rate appears not only near the toe of the breach but also in the gate area, which causes the sliding process and forms a slip surface. The slip surface is the interface between the un-yield and yield parts. In the bank of homogeneous mud, the slip surface of failure can be determined from the empirical method, which follows the arc of a circle that usually intersects the toe of the bank (Sun et al., 2008;Fredlund et al., 2012). However, the slip surface is developed automatically by MBM. It is worth a more profound study in the future. Error! Reference source not found. shows the strain rate profiles of BM and MBM. The slip surface (Error! Reference source not found. (b) at t = 10 s), as well as the interface between the plug/sheared zones (Fig 10 (b) and Fig 13 (b) at t = 40 s),

C14
can be identified in the results of MBM. From the comparisons of Fig 13 (a) and (b) at t = 10 s, and also Fig (a) and (b) at t = 10 s, we can see that the slip surface is relatively sharp in the MBM results than the ones in MB.
Referee #1: L314: References: Improve including suggested. Answer: Thank you very much. The references have been included.
Referee #1: Fig. 4: Interchange panel a and b. Fig. 6: panel d very basic, and e is in Fig. 7. So, you could remove these panels. Answer: Thank you very much. Those figures have been updated as suggested.
Referee #1: Fig. 7: "liquefied tailing" -> "deposited tailing". Local elevation changes in b (field photo) is very strong, however in a (simulation) it is very weak. Explain, and try to improve. Answer: Thanks for the comments. The word "liquefied tailings" was adopted directly based on the reference Jeyapalan et al., (1983). The aerial photo of the FGT66 has been updated to a clear one. In the photo, the pattern of white lines indicates the moving speed of the moving tailings before the freezing time. However, the pattern presented in the numerical simulation is the free-surface elevation. The result from MBM matches the observation well in terms of the flood distance, flood speed, freezing time, and the boundary of the flooded area.
Referee #1: Fig. 9: Panel with time = 30 s: Not the same time, and not the geometry, this is somehow similar to the real geometric pattern in Fig. 7 b. I would expect to see similar pattern of the final geometry. Answer: The white line segments shown in the aerial photo indicates the deformation of the tailing. It also recorded the velocity during the flooding process. Fig 9 at t = 30 s also shows the velocity profile of the free-surface from MBM. This result indirectly matches the field observation in terms of the fast-moving speed at the centerline area. However, the free-surface velocity will gradually reduce to zero as the freezing time is approaching. The white line segments in the aerial photo are not the elevation contours. It is difficult to extract the elevation profile directly from the aerial photo. In this paper, the flood distance and the flood area C15