## 2 Rheological Model and Numerical Algorithm
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:  
\[\frac{\partial\overline{u}_{i}}{\partial\overline{x}_{i}}=0 \tag{4}\]  
\[\frac{\partial}{\partial t}(\rho\overline{u}_{i})+\frac{\partial}{\partial x_{j }}\Big{(}\rho\overline{u_{i}u_{j}}\Big{)}=-\frac{\partial\overline{P}}{ \partial x_{i}}+\frac{\partial}{\partial x_{j}}\bigg{[}\mu_{e}\bigg{(}\frac{ \partial\overline{u_{i}}}{x_{j}}+\frac{\partial\overline{u_{j}}}{x_{i}}\bigg{)} \bigg{]}+\rho g_{i} \tag{5}\]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; \(\rho\) is the density; and \(\mu_{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 \(\mu_{e}\) is defined as:  
\[\mu_{e}=\mu\left(\dot{\gamma}\right)+\mu_{t} \tag{6}\]  
where \(\mu\left(\dot{\gamma}\right)\) is the rheological viscosity of mud and \(\mu_{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:  
\[\mu_{t}=\rho l_{S}^{2}\overline{S}=\rho(C_{s}\Delta)^{2}\overline{S} \tag{7}\]  
where \(l_{S}\) is the Smagorinsky length scale, which is a product of the Smagorinsky coefficient \(C_{s}\) and the filter width \(\Delta\); \(\overline{S}\) is the characteristic filtered rate of strain:  
\[\overline{S}=\sqrt{\frac{1}{2}\dot{\gamma_{ij}}\dot{\gamma_{ij}}} \tag{8}\]  
In general, \(C_{s}\) varies from 0.1 to 0.2 in different flows. The present simulations use a value of 0.15. \(\Delta\) is the filter width. Infinite volume discretization \(\Delta\) is the grid size:  
\[\Delta=(\Delta x_{1}\times\Delta x_{2}\times\Delta x_{3})^{1/3} \tag{9}\]  
where \(\Delta x_{1}\), \(\Delta x_{2}\), \(\Delta 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 \(m\)th 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 \(m\)th material; \(0<f_{m}<1\) if the cell contains the interfaces of the \(m\)th material; \(f_{m}=0\) if the cell contains no \(m\)th material. The VOF equation is given by Equation (10):  
\[\frac{\partial(f_{m})}{\partial t}+\
abla\cdot(u_{i}f_{m})=0 \tag{10}\]  
In this study, the interfaces between different materials are solved by the VOF method.