## 2 Rheological Model and Numerical Algorithm  
For flow rheology, the Bingham model (BM) has been widely used to simulate mudflows [4; 29; 35], lava flows [36], landslides [37], and multi-phase mass flows [16]. The rheological properties of BM can be presented as [38; 39]:  
\[\mu\Big{(}\dot{\gamma}\Big{)}=\left\{ \begin{array}{cc}\mu_{A}=\infty\;and\; \dot{\gamma}=0,&\mbox{if }\tau<\tau_{y}\\ \mu_{B}+\frac{\tau_{y}}{\sqrt[2]{\dot{\gamma}}_{ij}\dot{\gamma}_{ij}}\;and\; \dot{\gamma}>0,&\mbox{if }\tau\geq\tau_{y}\end{array} \right. \tag{1}\]  
where \(\mu_{A}\) is the viscosity of the un-yield region, \(\mu_{B}\) is the viscosity of the yield zone, \(\tau_{y}\) is the yield stress, and \(\dot{\gamma}_{y}\) is the yield strain rate, \(\dot{\gamma}_{ij}=\frac{\partial\dot{u}_{i}}{\partial x_{j}}+\frac{\partial\dot {u}_{j}}{\partial x_{i}}\). The symbol \(\dot{\gamma}\) is the second invariant of the \(\dot{\gamma}_{ij}\), which is defined as \(\dot{\gamma}=\sqrt{\frac{1}{2}\dot{\gamma}_{ij}\dot{\gamma}_{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]:  
\[\mu\Big{(}\dot{\gamma}\Big{)}=\left\{ \begin{array}{cc}\mu_{A}=\frac{\tau_{y}} {\dot{\gamma}_{y}},&\quad\text{if }\dot{\gamma}<\dot{\gamma}_{y}\\ \mu_{B}+\frac{\tau_{y}}{\sqrt{\frac{1}{2}\dot{\gamma}_{ij}\dot{\gamma}_{ij}}}, &\quad\text{if }\dot{\gamma}\geq\dot{\gamma}_{y}\end{array} \right. \tag{2}\]  
Mathematically speaking, when \(\dot{\gamma}_{y}\) approache zero, the CBM approaches BM. If the chosen \(\dot{\gamma}_{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 \(\mu_{A}\) and a larger \(\dot{\gamma}_{y}\) are required. The larger \(\mu_{A}\) plays a role of keeping the rigid shape. The larger \(\dot{\gamma}_{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:  
\[\mu\Big{(}\dot{\gamma}\Big{)}=\left\{ \begin{array}{cc}\mu_{A}>\frac{\tau_{y} }{\dot{\gamma}},&\quad\text{if }\dot{\gamma}<\dot{\gamma}_{y}\\ \mu_{B}+\frac{\tau_{y}}{\sqrt{\frac{1}{2}\dot{\gamma}_{ij}\dot{\gamma}_{ij}}}, &\quad\text{if }\dot{\gamma}\geq\dot{\gamma}_{y}\end{array} \right. \tag{3}\]  
In this model, the yield stress \(\tau_{y}\) and yield viscosity \(\mu_{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, \(\mu_{A}\) is chosen to be infinite based on the suggestions of [PERSON] [34], [PERSON] [45], and [PERSON] [28]. In this paper, the infinite number of viscosity \(\mu_{A}=10^{10}\) Pa s is chosen by the sensitivity analysis. The values of yield strain rate \(\dot{\gamma}_{y}\) are also discussed in Section 5.2. \(\dot{\gamma}_{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.