# Spatial Evolution Analysis of Tailings Flow from Tailings Dam Failure Based on MacCormack-TVD
## 2 Tailings Pond Breach Modeling
### Rheological Equations of Tailing Sand Flow  
The depth-integrated mass conservation Equation (4) and the momentum conservation Equations (5) and (6) are reformulated in vector form according to the MacCormack finite difference format with total variation diminishing (TVD) [25]:  
\[\frac{\partial X}{\partial t}+\frac{\partial F}{\partial x}=S;\ \frac{\partial X}{ \partial t}+\frac{\partial G}{\partial y}=T \tag{4}\]  
Figure 1: Schematic of the model.  
\[\text{where}\ X=\begin{cases}h\\ q_{x}\\ q_{y}\end{cases};F=\begin{cases}\frac{q_{x}}{h}+\frac{k_{xy}g_{z}h^{2}}{2}\\ \frac{q_{x}q_{y}}{h}+\frac{k_{xy}g_{z}h^{2}}{2}\end{cases};G=\begin{cases}\frac{ q_{y}}{q_{z}q_{y}}\\ \frac{q_{y}^{2}}{h}+\frac{k_{xy}g_{z}h^{2}}{2}\end{cases};\\ S=\begin{cases}-g_{z}h\frac{\frac{q_{z}}{q_{y}}}{1+\left(\frac{q_{x}}{2h}\right) +\left(\frac{q_{z}}{q_{y}}\right)}+S_{fx}\\ 0\end{cases};T=\begin{cases}0\\ 0\\ -g_{z}h\frac{\frac{q_{z}}{q_{y}}}{1+\left(\frac{q_{z}}{2h}\right)+\left(\frac{ q_{z}}{q_{y}}\right)}+S_{fly}\end{cases}.\end{cases}\]  
In the formula, \(q_{x}=h\overline{v}_{x},q_{y}=h\overline{v}_{y}\) denote the flow rate per unit width in the x-direction and y-direction; \(S_{f}=\left(S_{fx},S_{gy}\right)\) is the frictional resistance.  
The solution at moment \((n+1)\Delta t\) is:  
\[X_{i,j}^{n+1}=L_{x}^{2}\left(\frac{\Delta t}{2}\right)L_{y}^{2}\left(\frac{ \Delta t}{2}\right)L_{y}^{1}\left(\frac{\Delta t}{2}\right)L_{x}^{1}\left( \frac{\Delta t}{2}\right)X_{i,j}^{n} \tag{5}\]  
In the x-direction, \(L_{x}^{1}\) denotes the first execution of a \"prediction-correction-averaging\" solution process, and \(L_{x}^{2}\) denotes the second execution of the same process. Similarly, \(L_{y}^{1}\) and \(L_{y}^{2}\) denote the first and second execution of a \"prediction-correction-averaging\" solution process in the y-direction, respectively. Since the \"prediction-correction-averaging\" solving process is identical in both directions, only the process in the x-direction will be illustrated here. The solving process is as follows:  
First, perform a prediction step:  
\[X_{i,j}^{p}=X_{i,j}^{n}-\left(F_{i,j}^{n}-F_{i-1,j}^{n}\right)\cdot\frac{ \Delta t}{2\Delta x}+S_{i,j}^{n}\cdot\frac{\Delta t}{2} \tag{6}\]  
A correction step is being performed:  
\[X_{i,j}^{c}=X_{i,j}^{n}-\left(F_{i+1,j}^{p}-F_{i,j}^{p}\right)\cdot\frac{ \Delta t}{2\Delta x}+S_{i,j}^{p}\cdot\frac{\Delta t}{2} \tag{7}\]  
Finally, the values obtained from the prediction and correction steps are averaged to obtain the numerical solution:  
\[X_{i,j}^{n+1/2}=\frac{\left(X_{i,j}^{n}+X_{i,j}^{n}\right)}{2}+\left[G\left( \tau_{i,j}^{c}\right)+G\left(\tau_{i,j}^{c}\right)\right]\Delta X_{i+1/2,j}^{ n}-\left[G\left(\tau_{i,-1,j}^{c}\right)+G\left(\tau_{i,j}^{c}\right) \right]\Delta X_{i-1/2,j}^{n} \tag{8}\]  
where \(G(x)=0.5\times C\times[1-\varphi(x)]\).  
In the formula,  
\[\varphi(x)=max(0,\min(2x,1))\] \[C=\begin{cases}C_{r}\times(1-C_{r})\ \ C_{r}\leq 0.5\\ \In the formula, \(s=\frac{\rho_{w}}{\rho}\); \(\rho_{voellmy}\) is the Voellmy friction coefficient of the fluid; \(C_{z}\) is the Chezy coefficient; \(\rho_{w}\) is the density of the water; \(h\) is the height of the fluid; \(g_{z}\) is the gravitational component acting perpendicularly to the slope; \(\sigma\) is the angle of basal friction.