Depth-Integrated Two-Phase Modeling of Two Real Cases: A Comparison between r.avaﬂow and GeoFlow-SPH Codes

: Due to the growing populations in areas at high risk of natural disasters, hazard and risk assessments of landslides have attracted signiﬁcant attention from researchers worldwide. In order to assess potential risks and design possible countermeasures, it is necessary to have a better understanding of this phenomenon and its mechanism. As a result, the prediction of landslide evolution using continuum dynamic modeling implemented in advanced simulation tools is becoming more important. We analyzed a depth-integrated, two-phase model implemented in two different sets of code to stimulate rapid landslides, such as debris ﬂows and rock avalanches. The ﬁrst set of code, r.avaﬂow, represents a GIS-based computational framework and employs the NOC-TVD numerical scheme. The second set of code, GeoFlow-SPH, is based on the mesh-free numerical method of smoothed particle hydrodynamics (SPH) with the capability of describing pore pressure’s evolution along the vertical distribution of ﬂowing mass. Two real cases of an Acheron rock avalanche and Sham Tseng San Tsuen debris ﬂow were used with the best ﬁt values of geotechnical parameters obtained in the prior modeling to investigate the capabilities of the sets of code. Comparison of the results evidenced that both sets of code were capable of properly reproducing the run-out distance, deposition thickness, and deposition shape in the benchmark exercises. However, the values of maximum propagation velocities and thickness were considerably different, suggesting that using more than one set of simulation code allows us to predict more accurately the possible scenarios and design more effective countermeasures.


Introduction
Every year, hundreds of landslide-related phenomena, including rock avalanches and debris flows, cause loss of human life and the destruction of much else [1]. The risk analysis of this kind of natural hazard requires the evaluation of the stabilization of the source area and the prediction of propagation behavior (runout distance, velocities, thickness, and shape of deposit). In this study, we dealt with the latter issue by using reliable tools for investigating the dynamics of these flows. Once the flow characteristics are predicted, we can develop more precise hazard maps and design more effective countermeasures. There exist many different debris-resisting structures to cope with the threat related to debris avalanches and debris flows, such as bottom drainage screens [2] and artificial barriers [3].
Once unstabilized areas where mass flows are likely to release have been identified, existing programs based on dynamic numerical models are employed to predict flow behaviors and identify possible impact areas. Nowadays, there exists a large number of numerical models [4][5][6] developed to be useful for risk assessment and runout analysis, and the choice of the most suitable simulation code, given the advances in computational capabilities, has raised another important issue. It is always recommended to use more than one simulation to predict specific events or a detailed back-analysis. This allows us to compare numerical results, understand the main differences, and assess potential risks more precisely. Regardless of the selected simulation code's capabilities, many other factors, including the choice of rheological law [7], the pore-water pressure evolution [8], the presence of erosion [9], and accurate definitions of parameters, affect the results significantly.
This research's motivation was the significant challenges researchers and users face when attempting to obtain the same results from two different sets of code. We also aimed to evaluate the performance of a depth-integrated, two-phase model implemented in two different sets of code, and investigate the capabilities of these two widely validated sets of code, GeoFlow-SPH [10] and r.avaflow [11], which are based on the same mathematical model have with fundamentally different numerical solvers. The most important thing was to comprehend the methodologies behind the two sets of code, in order to choose an appropriate rheological law and appropriate geotechnical properties, rather than using the same rheological and empirical laws for both sets of code. This is because engineers will use the numerical results obtained for mapping dangerous areas and designing structural countermeasures.
The paper is structured as follows. The first section is dedicated to introducing the two different continuum-based programs and briefly describing the implemented depthintegrated, two-phase mathematical model consisting of the balance equations of mass and linear momentum, and suitable rheological and empirical equations. In the next section, we report the simulations of two different real events. The simulations used the selected sets of code while aiming to evaluate the ability of the depth-integrated two-phase model to reproduce rapid, large-scale landslides. Finally, the obtained numerical results are compared and discussed.

Introduction
Nowadays, several advanced simulation tools [4][5][6] are available to define areas in danger of landslides, foresee landslides' propagation paths, and obtain parameters to design structural and non-structural countermeasures in order to prevent the damage they cause. Here, we have selected GeoFlow-SPH and r.avaflow, which were developed for the propagation modeling of rapid landslides, such as debris flows and rock avalanches, and have been widely validated through the back analysis of real events [11][12][13][14]. GeoFlow-SPH and r.avaflow both employ a two-phase model which uses an interacting solid and fluid mixture; it is similar to the model proposed by Pitman and Le (2005) [15], and based on the depth-integrated mathematical model of Zienkiewicz and Shiomi (1984) [16].
In propagation modeling, the mathematical formulations implemented in sets of code are usually simplified and reduced to two dimensions through depth integration approximations. Since many flow-like landslides can be long and wide compared to their depths, the governing equations can be integrated along the vertical axis. The resulting 2D depth-integrated model provides a good balance between accuracy and time consumption. GeoFlow-SPH and r.avaflow are also based on integrated solutions of the balance equations of linear momentum and mass using the shallow water flow assumption.
We start by explaining how the mentioned simulation programs can be represented. Since the equations are depth-integrated, the particles cannot be defined along the vertical axis. Thus, in depth-integrated models, the flowing mass is divided into a finite number of columns represented by particles. Therefore, two heights corresponding to the solid (h s ) and fluid (h w ) phases will appear in each model by considering the assumption that the fluidized soil is fully saturated, and the sum of these partial heights is the total height (h) of the mixture. Therefore, the phases are simulated in different layers and interact with each other by interaction force. In Figure 1, the evolution of both phases' heights with respect to the total height is depicted. Letn w =n andn s = (1 −n) be the fractions of volume (V) occupied by the water and solid particles, respectively. The heights of the solid phase (h s ) and fluid phase (h w ) are given by:  [11,17], is an open-source computational tool designed for simulating the dynamics of rapid mass flows, including avalanches and debris flows. It represents a GIS-based computational framework, and the sourcing mass and topography can be defined through raster maps. The tool employs the NOC-TVD numerical scheme [18] along with the Pudasaini two-phase model [19], or a simple singlephase, Voellmy-type model. Besides, it considers the interactions between the phases and the entrainment of material along a path. The interested reader will find in the article by Mergili et al. (2017) [11], full details about the r.avaflow numerical code.
In this tool, the flow's evolution is computed through depth-averaged, two-phase balance equations developed by Pudasaini (2012) [19]. The mass balance equations for the solid and fluid phases, respectively, are given by: The momentum equations for the solid and the fluid phases are written in conservative form, respectively, as: where u vm s = γC u f − u s and u vm f = (α s /α s )C u f − u s are the virtual masses for the solid and fluid phases, respectively. γ is the fluid to solid density ratio. C = (N vm − 1)/(n s /n w + γ) is the virtual mass force coefficient. N vm = (m w /h w − m s /h s )/(u w − u s ) is the virtual mass number. β s = K s g z (1 − γ) and β w = g z are the hydraulic pressure coefficients for the solid and fluid phases, respectively. K is the earth pressure coefficient and g z the gravity.
The source terms of Equation (3) are as follows: where g is the non-dimensional gravitational acceleration. δ the internal friction angle. p b s = g z and p b w = p b s /(1 − γ) are effective basal pressures for the solid and fluid phases. ε is the typical height (H) to the typical extent of debris flow (L). X is the viscous shearing coefficient for fluid and ξ the vertical distribution of the particle concentration. N R = gLρ w H /(n w ρ w ) and N R A = gLρ w H /(Aρ w ) are the quasi-Reynolds numbers and the mobility Reynolds numbers, respectively. A is the mobility of fluid at the interface. In the special case in which A = n w , these quasi-Reynolds numbers coincide. C DG is the generalized drag coefficient.
Next, a set of additional functions and laws are needed to complete the conservation of mass and momentum equations presented in the current section. Here, we provide the following items: (1) the basal friction laws, (2) the interaction force laws, and (3) the erosion laws.
Regarding the rheological models implemented in the r.avaflow code, the solid stress is computed based on Mohr-Coulomb plasticity, whereas the fluid stress is modeled as a solid volume-fraction-gradient-enhanced non-Newtonian viscous stress. These models rely on the empirical coefficients C AD for ambient drag and C FF for fluid friction. C AD is a coefficient to be multiplied with the frontal surface and the velocity of the flow to derive air resistance. C FF is Manning's coefficient used as the fluid friction coefficient, to consider the effects of the roughness of basal surface.
Entrainment is computed through an empirical law in which a user-defined entrainment coefficient (E s ) is multiplied with the total momentum (M s + M w ) at each raster cell. Then, the solid and fluid entrainment rates q E,s and q E,w can be obtained as follows: wheren s,E is the solid fraction of the entrainable material. In the r.avaflow, the changes of the basal topography are considered by computing the solid and fluid erosion heights, H E,s and H E,w , at each time step, as follows: where H Emax,s and H Emax,w are the maximum entrainable depths at the given cell, t is the time passed at the end of the time step, ∆t is the time step length, and β is the local slope of the basal surface. The solid-fluid interaction is considered by using the two-phase generalized drag model proposed by Pudasaini (2012) [19], a modified Pitman and Le model [15]. It is capable of taking into account the effect of fluid viscosity on drag. The generalized drag coefficient (C DG ) is given by: where V T = gd/γ is the terminal velocity. d the particle diameter. P ∈ (0, 1) is the parameter for a combination of solid and fluid-like contributions to drag resistance. F = γ(n w /n s ) 3 Re p /180. Re p = (ρ w du T )/η w is the particle Reynolds number. It is important to note that the excess pore pressure terms are disregarded in this model due to the assumption that dilation can be neglected or the permeability is large. Therefore, only hydrostatic fluid pressure is considered: ∂p/∂z = g z .
Numerical schemes are a fundamental tool for solving differential equations and obtaining approximations of physical problems. In the r.avaflow code, the system of equations is discretized on a staggered grid, in which the system is moved half of the cell size with every time step. Then, the total variation diminishing non-oscillatory central differencing (TVD-NOC) scheme [20][21][22] is employed to solve the model equations. The TVD-NOC has successfully been applied to a large number of mass flow problems [18,21,23,24].

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, depthintegrated model recently developed by Pastor 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: where sub-index α denotes the solid (s) or fluid (w) phase,n α is the porosity of the solid or fluid phase (n s = 1 −n andn w =n), e R is the erosion rate, and is gravitational force and its axis is vertical and points upwards. In the diffusion term of Equation (9), the averaged pressures (P α ) acting on solid or fluid phases are defined as: which indicates that, in our special case of interest, total pore-water pressure (p w ) is composed of a hydrostatic part p hyd , which varies linearly from zero at the surface to ρgh at the bottom, and an excess pore-water pressure (∆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 (τ 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: where h is the propagation height, φ B the basal friction angle,v the depth averaged flow velocity, ξ the turbulence coefficient, and ∆p b w 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 τ 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. It is based on an algorithm in which the total volume of debris increases in accordance with a specified rate and is given by: which indicates a direct proportionality between the entrainment rate and the product of depth average velocity (v) and mobilized soil depth (h). E s can be obtained directly from the initial and final volumes of the material and the distance traveled as E s ≈ ln((V final /V 0 )/(distance)). The interaction between solid and fluid phases is given by the depth-integrated value ofR. It is provided by the Anderson and Jackson law [27], which is capable of considering large relative velocities. It is given by: where V T is the terminal velocity, g is the acceleration of gravity, and m is a constant. In depth-integrated models, the vertical structure of the magnitudes is lost as the only available information is their depth-integrated values. To overcome this limitation, an additional equation has been implemented in GeoFlow-SPH's model to describe how pore pressure evolves over time and depth [28]. Here, we recall the consolidation equation describing the evolution of pore pressure along the vertical axis, as follows: where ρ d (= (1 − n)(ρ s − ρ w )) is the effective density, c v the consolidation coefficient, and E m the odometric modulus. Based on the above equation, the pore-water pressures are also influenced by the propagation height and its velocity, and the porosity variations.
The consolidation parameter of c v = k v k w plays an important role, as a high value of the coefficient allows the rapid dissipation of excess pore-water pressure through the consolidation process. In the article by Pastor (2021) [26], the interested reader will find a detailed explanation of the consolidation equation. Next, we briefly describe a Lagrangian meshless numerical method that is used to discretize the presented depth-integrated equations. Smoothed particle hydrodynamics (SPH) uses it to transform the problems that are basically in the form of partial differential equations (PDEs) to forms suitable for particle-based simulation. SPH was first invented by Lucy (1977) [29] and Gingold and Monaghan (1977) [30] to model astrophysical problems. This technique has been applied in many areas due to its ability to model complicated cases which involve large-displacement deformations, such as the modeling of fast landslides in solid mechanics [4,31,32]. Good reviews can be found in the texts of Liu and Liu (2004) [33] and Li and Liu (2003) [34].
In this model, these 3D problems are transformed into 2D forms by applying a depthintegrated model. As a result, it provides an excellent combination of accuracy and computational time. This technique has been successfully applied to debris flows by Pastor et al. (2014) [10] and Cascini et al. (2014Cascini et al. ( , 2016 [35,36], and debris avalanches by Cuomo et al. (2014Cuomo et al. ( , 2016 [37,38].
To discretize the propagating mass in the SPH method, the first step is to present them as a set of nodes, as depicted in the figure, having individual material properties. In this paper, the mathematical equations are discretized while using the two-phase SPH technique to deal with two phases of flow involving solid and water particles. Therefore, two sets of nodes are introduced, one to represent the solid particles' movement and another to represent the movement of the fluid particles, as depicted in Figure 2. Then, an interpolation process calculates the relevant properties on each node over neighboring nodes through a kernel function W I J = W x I − x J , h without defining any element. The formulation can be expressed as: where the function f (x I ) is approximated at a position vector x in space.
Once the information is stored into nodes or particles, the integral interpolation of SPH kernel approximation is replaced with a discretized form of summation over all the particles within the region of compact support of kernel function. It is expressed as follows: where the infinitesimal volume dx J of the continuous integral representations is replaced by the volume of a neighboring particle ω J .
In the GeoFlow-SPH code, two alternative methods are applied to approximate the pore pressures inside a landslide. In the first alternative, similarly to the most depthintegrated models, the consolidation equation is solved using a quarter-cosine shape function that fulfills boundary conditions with a zero value at the surface and zero gradients at the basal surface to approximate the vertical distribution of pore-water pressure. By neglecting the porosity variations at the basal surface, the following first alternative consolidation equation can be obtained: where for simplification, we introduce β = c v π 2 /4h 2 .
The second  [26] extended this one-phase model to a two-phase model in order to fully approximate the pore pressures inside a landslide. In this technique, a finite-difference mesh, incorporated at each SPH node representing solid particles, is used to discretize pore-water pressures along the vertical axis, as depicted in Figure 3. The consolidation equation given in Equation (14) can be rewritten in a way that is more suitable for the formulation of FDM, as follows: where the consolidation equation is discretized by making two changes: (i) height variations and (ii) porosity variations. In this paper, we apply the first alternative model for the cases that are not needed to take into account basal surface permeability and porosity variations of debris flows.

Case Studies
This section presents two real case studies and discusses the results obtained from back-analyses performed by the two programs presented in the last sections. The first benchmark exercise was the prehistoric Acheron rock avalanche, which has been backanalyzed by Mergili et al. (2017) [11] using r.avaflow code. This real case was used to assess the validity and performance of the GeoFlow-SPH code. Another validation exercise was the Sham Tseng San Tsuen debris flow, which has been back-analyzed by Pastor et al. (2021) [26] using two-phase SPH modeling implemented in GeoFlow-SPH code. The latter case study was also used to assess the performance of the sets of code by comparing their numerical results.

Acheron Rock Vvalanche
The Acheron rock avalanche occurred in Canterbury, New Zealand, approximately 1100 years BP [39,40] based on radiocarbon dating analysis. An earthquake may have triggered the landslide, as it is located close to fault zones. In Figure 4, an overview of the Acheron rock avalanche is shown. Validation of sets of code relies on the availability of reliable information, including the impact area or deposition area of the event; and estimated deposition thickness, velocities, and final run-out.
The Acheron Rock Avalanche has been back-analyzed by Mergili et al. (2017) [11] using the r.avaflow code. The same back-calculated parameters were used here for both sets of code to simulate the rock avalanche. Flow parameters and coefficients required for both sets of code are given in Table 1. Table 1. Material parameters used in the analysis of the avalanche.

Parameter (Symbol) GeoFlow-SPH r.avaflow
Density of the of solid particle (ρ s ) 2700 kg/m 3 2700 kg/m 3 Density of the of water particle (ρ w ) 1000 kg/m 3 1000 kg/m 3 Mixture density (ρ) 2350 kg/m 3 -Internal friction angle (δ) The densities of solid and fluid particles were taken as ρ s = 2700 kg/m 3 and ρ w = 1000 kg/m 3 , respectively. In the GeoFlow-SPH code, the initial porosity is computed based on densities. On the other hand, the r.avaflow code uses the defined heights of fluid and solid phases to obtain porosities. For both sets of code, an initial porosity of 0.2, for which ρ = 2350 kg/m 3 , was considered.
In the previous numerical analysis [11], the best results were obtained with the dynamic basal friction angle of 17 • . In the r.avaflow code, the simulation was based on Mohr-Coulomb-type plastic flow and was performed with an internal angle of friction of 35 • and a bed friction angle 17 • . On the other hand, the GeoFlow-SPH code uses the voellmy rheological model to analysis case studies. Concerning its rheological parameters, the basal friction angle of 17 • and no turbulence coefficient were considered for this particular case. In both sets of code, the fluid phase is modeled based on the Manning formula. However, the fluid friction coefficient was neglected in this study.
The output of r.avaflow is discretized on the basis of geographic information system (GIS) coordinates and consists of raster maps. One of the advantages of the r.avaflow code is performing the simulation model in a GIS in order to manage geo-referenced data and facilitate the use of the results. In the GeoFlow-SPH code, post-processing tools are used to deal with the visualization of output. The main result is the evolution of solid and fluid flow heights, and the outputs are the values of velocity and depositional height computed at each time step. Optionally, total flow pressures and kinetic energies are generated in the r.avaflow code; relative pore-water pressure and total pore-water pressure are generated in the GeoFlow-SPH code. Figure 5 illustrates the comparison between the numerical results using the r.avaflow and GeoFlow-SPH sets of code, with the parameter values given in Table 1, at different time steps. As shown in the figure, with the same geotechnical parameter values, the two sets of code gave relatively similar deposit distributions, velocities, runout distances, and final deposition thickness values. Both sets of code show the traveling of small portions of the volume in the wrong direction (northern). It should be recalled that the event is an ancient rock avalanche, and perhaps there exist some unrecognizable deposition volumes in the field. The deposition volume was considered to be the field deposit estimate of 8.9 × 10 6 m 3 for both sets of code using the back-calculated erosion rate given in Table 1. As described in Sections 2.2 and 2.3, the erosion laws implemented in both sets of code are of empirical type.
During the propagation stage, a significant difference was observed in the numerical results of maximum height evolutions. However, the r.avaflow and GeoFlow-SPH sets of code had similar results regarding the shape of the final deposit and depositional height distribution. The maximum deposit thicknesses calculated by r.avaflow and GeoFlow-SPH sets of code were 31.8 and 27.3 m, respectively.
As shown in Figure 5, both sets of code reached a satisfactory approximation of the depositional area shape. However, the r.avaflow code gave a larger lateral mass spread than GeoFlow-SPH. The excessive lateral spreading might be solved by choosing appropriate combinations of internal friction angle δ and basal friction angel φ B . The GeoFlow-SPH code successfully modeled the final deposit, although an isotropic state of stress was assumed, and the earth pressure coefficient (k-value) was considered equal to 1. Lateral spreading also depends on the evolution of pore-water pressure, which the GeoFlow-SPH is capable of considering.
As depicted in Figure 6a, the r.avaflow code overestimated the height values with respect to the GeoFlow-SPH code until 95 s (time of reach). After that, the numerical result of maximum deposit thickness obtained by the r.avaflow decreased during the lateral spreading of the flowing mass until the numerical results of both sets of code became relatively close in terms of the shape of the final deposit and depositional height distribution.  Figure 6b shows the differences in maximum velocity values obtained from the GeoFlow-SPH and r.avaflow programs at the same time steps evaluated during the whole simulation. Both sets of code indicated that maximum velocity peaked at the time step of 20 s and achieved the maximum velocity of about 71 ms −1 . The results show that the r.avaflow and GeoFlow-SPH programs give approximately the same maximum velocities values until the time step of 95 s (time of reach). It is interesting to note that the most significant differences in computed velocities were during the spreading over the deposition area, where the GeoFlow-SPH code predicted that all the particles deposited at the time step of 95 s; on the other hand, at the same time the r.avaflow code indicated that some particles were still traveling at high speed-31 ms −1 .
The numerical analyses clearly show that the obtained calibrated values of rheological parameters can be transferred from one code to another, as both programs are based on the same mathematical formulation, although they employ different empirical laws and numerical solutions.

Sham Tseng San Tsuen Debris Flow
This debris flow took place in Hong Kong on 23 August 1999. This case was proposed as a benchmarking exercise by the Hong Kong Geotechnical Office with the objective of evaluating the precision of numerical and constitutive models [41]. They also provided a detailed digital map of terrain elevation, including the mobilized mass's original position and the reservoir's final position.
The event took place after heavy rain. During the 24 h before the event, the precipitation was 479 mm. The debris flow consisted of 600 m 3 of channeled material, according to the report. It originated from three shallow landslides, of which one was much larger than the others. In fact, the data provided included only this main source. The debris traveled downhill along an incised rocky stream course with negligible debris entrainment, and it destroyed several one-story buildings at the end of the canal. Figure 7 represents a map of the terrain showing the position of landslide A, the track, Sham Tseng San Tsuen village's location, and the house labeled as "No. 38", which was destroyed. More information can be found in a geotechnical report provided by Wilson et al. (2005) [42].
This simulation aimed to consider excess pore-water pressures in the two-phase model implemented in GeoFlow-SPH. In the previous simulation, it was assumed that the permeability was very large, and the excess pore pressure terms were neglected.
Densities of soil particles, pore fluid, and mixture were taken as 2400 Kg/m 3 , 1000 Kg/m 3 , and 2000 Kg/m 3 . The drag laws that have been implemented in both sets of code were based on the Anderson and Jackson law [27] model. In this case, the terminal velocity (V T ) and m were taken as 10 −4 ms −1 and 1, respectively.
According to the report [42], the triaxial tests show that the material did not present cohesion, and the friction angle varied from 37 • to 38 • , for which we found tan φ B = 0.75 − 0.78. Interestingly, the report states that "the distal end of the debris deposition had a travel angle of about 24 • ", which is much smaller than the friction angle found in laboratory tests. For the two-phase model developed by Pastor (2021) [26], the values of the friction angle and travel angle coincided well with those obtained in the reports by using the following equation: where φ B is the friction angle and φ app the travel angle. The rest of the parameters were similar to those in the previous case study.  In Figure 8, we can observe the larger velocities of the model simulated in the GeoFlow-SPH code, capable of considering excess pore-water pressures. The main reason for these differences was the pore-water pressure's evolution in the GeoFlow-SPH code. As described in Section 2.3, the consolidation equation was employed in the two-phase modeling of debris flow to compute the quantity of excess pore-water pressure, which has been implemented in the balance of momentum equation.
Besides, the frictional rheological equation implemented in the numerical code is capable of considering the effects of the pore pressure dissipation and generation. As can be seen in Equation (11), the generation of pore-water pressure is similar in effect to reducing the friction angle, and the dissipation of pore-water pressure is equivalent to increasing the friction angle. To obtain similar numerical results, it is necessary to use a lower value of basal friction angel in the r.avaflow code. Therefore, improvements to the numerical results obtained by the r.avaflow code can be achieved by modifying the values of rheological parameters.
In Figure 9, we provide the numerical results obtained from the GeoFlow-SPH code regarding the basal excess pore-water pressures relative to liquefaction. At the start, the relative pore-water pressure p rel w was taken as 1, indicating that the flowing material was completely liquefied. Then, the pore-water pressure dissipated during the propagation stage until the debris-flow reached the deposition area. In the authors' opinion, the back-analysis of both case studies indicates that the numerical results of r.avaflow are sensitive to variations in the basal friction angle φ B and the initial solid fraction n 0 ; and on the other hand, the numerical results of GeoFlow-SPH are more sensitive to the consolidation coefficient c v obtained by using parameters of terminal velocity V T and stiffness of the mixture k v ; consequently, the characteristic times of propagation and consolidation are the key aspects of the GeoFlow-SPH model.
There exist many different structural countermeasures to reduce the impacts of landslides. Bottom drainage screens [2] could be used as an effective control to work against channelized flows with narrow widths, such as the Sham Tseng San Tsuen debris flow case. Artificial barriers and walls [3] could be used as protection against rock and debris avalanches, such as the Acheron rock avalanche, with widespread deposition.

Conclusions
This paper described two sets of code, r.avaflow and GeoFlow-SPH, which are similar in mathematical formulation and differ by numerical implementation method. The two programs have been used to back-analyze two real cases that occurred in New Zealand in ancient history and Hong Kong in 1999. In previous works, the former has been backanalyzed with the r.avaflow code to validate the implemented depth-integrated, two-phase model, and the latter has been back-analyzed through the GeoFlow-SPH code to evaluate the two-phase model while considering the space-time evolution of pore-water pressure. The best fit values of geotechnical parameters, which were achieved through one of the programs in the previous works, were used in another program in this study.
The aims of the research were (1) to evaluate the prediction capabilities of programs, which are based on a depth-integrated two-phase model, by comparing the numerical results with available field data, and (2) to investigate whether the same geotechnical parameters calibrated in the one numerical code can be transferred to other numerical sets of code.
The numerical results were compared to the field data, and they evidenced that both sets of code were capable of reproducing the run-out distance, deposition thickness, and shape for the real cases, although some parameters are considered in one numerical code but not in another, such as internal friction angle, earth pressure coefficients, and excess pore pressure evolution. Thus, the results of both sets of code can be considered satisfactory. The first numerical analyses showed that the same values of rheological parameters can be used for different sets of code based on similar mathematical formulae. However, the mathematical formulations used for the second case study were different due to the implementation of a consolidation equation in the Geoflow-SPH code for considering the space-time evolution of pore-water pressure.
There were also considerable differences between the numerical results obtained by the two programs. In the first case study, the GeoFlow-SPH code underestimated the maximum propagation thickness compared to the r.avaflow code. In the second case study, the r.avaflow code underestimated the maximum propagation velocity compared to the GeoFlow-SPH code, which considers the pore-pressure evolution. These results confirm that using more than one simulation code is necessary to predict specific events to assess potential risks and design possible countermeasures.