## 2 Materials and Methods
### GeoFlow-SPH Numerical Code  
The GeoFlow-SPH code was developed in Madrid by an expert research team for almost a decade. It has previously been applied to theoretical [25], experimental [2], and real case histories [10]. The latest theoretical framework is based on a two-phase, depth-integrated model recently developed by [PERSON] et al. (2021) [26]. This model is a general approach for two-phase modeling capable of considering many essential physical aspects, such as the phenomena of pore-water pressure evolution and reproducing the dynamic behavior of debris flows by taking into account their soil properties, including permeability and volume stiffness. The two-phase model consists of balance equations of mass and linear momentum expressed, respectively, in quasi-Lagrangian form as follows:  
\[\frac{\partial h_{\alpha}}{\partial t}+\frac{\partial}{\partial x_{ai}}(h_{ \alpha}\bar{\sigma}_{ai})=\bar{n}_{\alpha}e_{R} \tag{8}\]  
\[ \begin{split}\rho_{a}h_{a}\frac{\bar{d}^{(\alpha)}\bar{\sigma}_{ a}}{dt}&=\rho_{a}\text{grad}\bar{P}_{a}-\biggl{(}\frac{1}{2}\rho_{w}b_{3}h^{2}- \Delta\bar{p}_{w}h\biggr{)}\text{grad}\bar{n}_{\alpha}\\ &\quad+\tau_{B}^{(\alpha)}+\rho_{a}bh_{\alpha}+\bar{R}h-\rho_{a} \bar{\sigma}_{a}\bar{n}_{\alpha}e_{R}\end{split}  \tag{9}\]  
where sub-index \(\alpha\) denotes the solid (s) or fluid (\(w\)) phase, \(\bar{n}_{a}\) is the porosity of the solid or fluid phase (\(\bar{n}_{s}=1-\bar{n}\) and \(\bar{n}_{w}=\bar{n}\)), \(e_{R}\) is the erosion rate, and \(\bar{h}_{\alpha}=h\bar{n}_{\alpha}\). \(b_{3}(=-\text{g})\)is gravitational force and its axis is vertical and points upwards. In the diffusion term of Equation (9), the averaged pressures \((\bar{P}_{a})\) acting on solid or fluid phases are defined as:  
\[\bar{P}_{s}=\frac{1}{2}b_{3}hh_{s}+\frac{\Delta\bar{p}_{w}ht\bar{n}}{ \rho_{s}}\ \ \text{ForSolid} \tag{10}\] \[\bar{P}_{w}=\frac{1}{2}b_{3}hh_{w}-\frac{\Delta\bar{p}_{w}ht\bar{n }}{\rho_{w}}\ \ \ \text{ForFluid}\]  
which indicates that, in our special case of interest, total pore-water pressure \((p_{w})\) is composed of a hydrostatic part \(\left(p_{hyd}\right)\), which varies linearly from zero at the surface to \(pgh\) at the bottom, and an excess pore-water pressure \((\Delta p_{w})\) which should be computed for each time and portion of space. Next, the mathematical equations given in Equation (9) are completed using a rheological law to compute the basal shear stress \((\tau_{B})\). In this study, the numerical analysis was performed through voellmy's rheological law, which has the same features as the frictional rheological model where the cohesion and viscous terms are disregarded. It is capable of considering the evolution of pore-water pressure at the basal surface. In the case of a pure frictional mass, the basal shear stress is given by:  
\[\tau_{B}=\left(\rho_{d}^{{}^{\prime}}gh-\Delta p_{w}^{b}\right)\frac{\bar{ \sigma}_{i}}{|\bar{\sigma}|}\tan\phi_{B}+\rho\frac{|\bar{\sigma}|}{\bar{\sigma} }\bar{\sigma}_{i} \tag{11}\]  
where \(h\) is the propagation height, \(\phi_{B}\) the basal friction angle, \(\bar{\sigma}\) the depth averaged flow velocity, \(\bar{\zeta}\) the turbulence coefficient, and \(\Delta p_{w}^{b}\) the excess pore-water pressure at the basal surface, which is computed by using a consolidation equation which is given later in this section. It can be seen in the above equation that the basal shear stress \(\tau_{B}\) will depend on the basal excess pore pressure, and it is modified in accordance with pore-water pressure's evolution at each node and time step. Take into account that the higher the pore-water pressure, the lower the shear stress at the basal surface. Regarding entrainment, the empirical erosion law of Hungr has been implemented in the GeoFlow-SPH code.