Modelling of Fluidised Geomaterials: The Case of the Aberfan and the Gypsum Tailings Impoundment Flowslides

The choice of a pure cohesive or a pure frictional viscoplastic model to represent the rheological behaviour of a flowslide is of paramount importance in order to obtain accurate results for real cases. The principal goal of the present work is to clarify the influence of the type of viscous model—pure cohesive versus pure frictional—with the numerical reproduction of two different real flowslides that occurred in 1966: the Aberfan flowslide and the Gypsum tailings impoundment flowslide. In the present work, a depth-integrated model based on the v-pw Biot–Zienkiewicz formulation, enhanced with a diffusion-like equation to account for the pore pressure evolution within the soil mass, is applied to both 1966 cases. For the Aberfan flowslide, a frictional viscous model based on Perzyna viscoplasticity is considered, while a pure cohesive viscous model (Bingham model) is considered for the case of the Gypsum flowslide. The numerical approach followed is the SPH method, which has been enriched by adding a 1D finite difference grid to each SPH node in order to improve the description of the pore water evolution in the propagating mixture. The results obtained by the performed simulations are in agreement with the documentation obtained through the UK National Archive (Aberfan flowslide) and the International Commission of large Dams (Gypsum flowslide).


Introduction
Flowslides are rapid flows, either saturated or unsaturated, where the material has a high compaction tendency, a low density and is characterised by a metastable structure. Since flow failures experience a sudden loss of strength followed by a very rapid development of deformations, their effects are often much more dramatic and devastating than other types of landslides. Thus, the prediction of flowslides' propagation distances, velocity and pore water pressure will be of great

Solid Skeleton Pore Fluid Coupled Mathematical Model
The starting point for the mathematical model considered in the present work, in order to describe the behaviour of the mixture that propagates within the landslide, is the well known v − p w Biot-Zienkiewicz model [28] which is composed of the balance momentum of the mixture (1) and the combined pore fluid balance of mass and linear momentum (2) ρ dv dt = ρb + div σ (1) div(−k w grad p w ) + div v = 0 Here v is the velocity of the solid skeleton, ρ is the mixture density, σ is the total Cauchy stress tensor while b is the gravity acceleration vector. The symbol d dt in (1) is the material time derivative with respect to the solid skeleton. It is defined by dv dt = ∂v ∂t + grad v · v where grad means gradient operator with respect to the spatial position. In both equations, the operator div is the divergence that in (2) is associated with grad. Moreover, in (2), p w is the pore fluid pressure and k w is the permeability coefficient. Incompressibility of pore fluid and solid skeleton particles has already been assumed in (2).
Most of the landslides are shallow, i.e., if L is a characteristic length of the landslide and H is a characteristic depth of the sliding mass, then ε = H/L 1. This geometric characteristic suggests, after a dimensionless analysis is performed, that the decomposition of the velocity field v [29] in a vertical consolidation component v c and a propagation component v p is possible and it is done as follows allowing the description of the landslide evolution as the combination of two different physical phenomena, namely, pore pressure evolution and propagation. Therefore, the mixture initially describe by (1) and (2), can be expressed by the following set of equations ρ dv p dt = ρb + div σ div v p = 0 (4) div v c = k w ∂ 2 p w ∂x 2 3 (5) where (4) governs the propagation of the mixture as an incompressible viscous fluid while (5) governs the evolution of the pore water pressure through a vertical consolidation process. The complete set of governing Equations (4) and (5) can be applied with confidence to the case of flowslides, where the mixture presents a medium permeability and consolidation while propagation phenomena are developed with a similar order of magnitude in time.
The principal component of volume changes within a spreading mixture governed by the set of Equations (4) and (5) is due to a vertical consolidation process. This behaviour, specific to flowslides, might not be applied to debris flows, where very high permeabilities may be observed and other sources of volume changes should be considered [30]. For the case of mudflows, where very low permeabilities are observed and the time required to develop consolidation is much larger than the time required for the propagation of the mixture, incompressible behaviour is usually observed. Therefore, if a mudflow is to be modelled, Equation (5) can be neglected keeping (4) as the main set of governing equations.
In order to be able to solve the Equations (4) and (5), appropriate initial and boundary conditions are needed. Regarding the boundary conditions (Figure 1), two different types of boundaries are considered: a no-slip condition at basal surface, defined by x 3 = Z(x 1 , x 2 , t), and a free surface, defined by of governing equations.
In order to be able to solve the Equations (4) and (5), appropriate initial and boundary conditions are needed. Regarding the boundary conditions (Figure 1), two different types of boundaries are considered: a no-slip condition at basal surface, defined by x h x ,x ,t Z x ,x ,t   . Figure 1. Reference system and notation used in the analysis [28].
As per the previous definition, the basal surface Z varies with time. Under this assumption, erosion can be considered in the mathematical description of the landslide by defining the erosion rate R e to be The system of partial differential Equations (4) and (5) are solved in the present work by numerical techniques. But first, in the following sections, system (4) is integrated into depth while Equation (5) is transformed into a diffusion-like equation.

Depth-Integrated Mathematical Model
In what follows, superindex p in (4) is dropped for brevity. Resolution of the above defined 3D boundary value problem is a formidable task if no reasonable simplifications are considered. In this context, depth-integrated models are a convenient simplification of 3D models, providing an acceptable compromise between computational cost and accuracy [30].
In order to obtain the depth integrated version of the system (4), the quasi-material derivative d dt is first introduced by the expression where j v is the j component of the average value of the velocity v over the flow depth h defined by Figure 1. Reference system and notation used in the analysis [28].
As per the previous definition, the basal surface Z varies with time. Under this assumption, erosion can be considered in the mathematical description of the landslide by defining the erosion rate e R to be The system of partial differential Equations (4) and (5) are solved in the present work by numerical techniques. But first, in the following sections, system (4) is integrated into depth while Equation (5) is transformed into a diffusion-like equation.

Depth-Integrated Mathematical Model
In what follows, superindex p in (4) is dropped for brevity. Resolution of the above defined 3D boundary value problem is a formidable task if no reasonable simplifications are considered. In this context, depth-integrated models are a convenient simplification of 3D models, providing an acceptable compromise between computational cost and accuracy [30].
In order to obtain the depth integrated version of the system (4), the quasi-material derivative d dt is first introduced by the expression where v j is the j component of the average value of the velocity v over the flow depth h defined by Integrating into depth, taking into account Leibniz's rule for differentiating integrals and considering the boundary conditions at the basal and free surface, the depth-integrated balance of mass reads as dh dt where e R is the erosion rate defined by (6), while div v = ∂v 1 ∂x 1 + ∂v 2 ∂x 2 . If a depth integration approach is applied to the balance of linear momentum and if an internal viscosity force within the incompressible mixture is neglected, as compared with the viscous resistance opposed by the basal surface to the displacement of the mass, the following expression is obtained [30] ρh Materials 2017, 10, 562

of 21
where g is the acceleration of gravity, ρ is the mixture density, h is the depth of the flow, v is the average velocity over h, Z is the basal surface, e R is the erosion rate and τ b is the basal shear stress that depends on the rheological law considered. Also, in (10),

Pore Pressure Evolution in Landslides
In what follows, superindex c in (5) is dropped for brevity. Bearing in mind the relation div v = dε v dt , where ε v is the volumetric deformation, it is assumed that the time rate of change of the volumetric deformation ε v can be related to the time rate of variation of effective confining pressure p by where K v is a suitable stiffness modulus, p = − 1 3 tr σ , and σ is the effective stress tensor. If the skeleton is elastic, K v is the elastic volumetric stiffness ratio. For example, for oedometric conditions, the volumetric stiffness ratio is considered as the oedometric modulus, K v = E m . Taking into account that p = p + p w , Equation (5) can be rewritten as In order to solve this equation, the landslide mass will be decomposed into differential elements of volume having a height h and a differential cross section dA at each time t, as shown in Figure 2.
The changes of the total mean confining pressure p are mainly caused by the height variation. So, for the differential volume in Figure 2, the total stress will vary as follows: Then, the total stress σ 3 , which depends on h, varies with it as Concerning the effective stress, it is observed that and by considering the relation dp = αdσ 3 , it is found that (16) which is the equation describing the evolution of the pore pressure along 3 x . Equation (16) has to be complemented with an initial and boundary condition at 3  x Z and 3   x Z h. For example, it can be considered zero at the surface and zero flow at the bottom.
Summarising, the mathematical model considered in the present work to reproduce the principal features of a landslide consists of: 2. A pore pressure evolution equation Before performing the numerical approach, a rheological model should be established in order to clarify the expression for the basal shear stress b  and its possible relation with the basal excess pore pressure. This will be clarified in the following section.

Introduction
Deformation of a soil column [28].
Summarising, the mathematical model considered in the present work to reproduce the principal features of a landslide consists of:

2.
A pore pressure evolution equation Before performing the numerical approach, a rheological model should be established in order to clarify the expression for the basal shear stress τ b and its possible relation with the basal excess pore pressure. This will be clarified in the following section.

Introduction
In this section, the behaviour of fluidised geomaterials in fast landslides will be discussed. Once failure has been triggered, the behaviour of the soil mass is closer to that of fluids than to solids. This is why rheological models are used to describe the behaviour of such hazards. There are many types of materials involved in fast landslides, from assemblies of rock blocks to mixtures of clay and water.
Because of the computational cost of full 3D models, researchers have favoured the use of simpler depth-integrated models, as described previously, where the flow structure is lost and the basal friction is obtained from the depth average velocity.
The purpose of this section is then to describe three rheological models, two of them widely used and a new one based on Perzyna viscoplasticity for the study of landslides' propagation. The scope is to understand how these models provide an expression of basal friction using the hypothesis of simple shear Infinite Landslide Model, once the depth average velocity is known.

Pure Cohesive Viscoplastic Fluid: Bingham Model
The Bingham model includes two material parameters, the yield stress below which the material does not flow, and the viscosity. It was introduced by Bingham [31] in 1922. The expression for the Bingham model (where τ y is the yield stress) is written as: Depending on the fluid phase viscosity, mudflows, lahars and debris flow can be modelled as viscoplastic fluids with Bingham-like models. Considering a Bingham fluid initially at rest and increasing the shear stress, the fluid will start moving only when the shear stress reaches τ y . This behaviour creates what is generally called a "plug" or a zone where the velocity is constant and the rate of deformation is zero.
Concerning the bottom friction, it is assumed that it can be approximated under the hypothesis of simple shear flow conditions. As described in [15], the shear stress at the bottom τ b can be related to the depth-averaged velocity with the following expression:

Pure Frictional Viscoplastic Fluid
Frictional viscoplastic fluids are used to model fast landslides where friction is important. If the cohesion is assumed to be zero and using it is easy to obtain h is the total height, τ b = ρgh sin θ is the basal shear stress with density ρ and s(z) = −ρ d g(h − z) cos θ tan ϕ is the strength along z, being ρ d = ρ − ρ w . The velocity profile can then be obtained as where v h is the velocity at the surface. The basal shear stress becomes then where s b is the shear strength at the bottom.

Perzyna-Based Rheological Model for Frictional Materials
Viscoplastic models were found to provide a suitable and more economic approach than classical plasticity models when computing failure loads and mechanisms [32]. In the case of soils, viscoplastic models have been applied both to cohesive [33][34][35] and frictional materials. They have been found to reproduce well slow landslide movements [32,35,36].
and the Perzyna elasto-viscoplastic models. In the latter models, the relation between the effective stress and the rate of deformation tensor is given by Above, D e is the elastic constitutive tensor, d is the rate of deformation tensor and d vp is the viscoplastic component. The viscoplastic component of the rate of the deformation tensor is given by Perzyna [37,38] as where represents the Macaulay brackets and γ is the fluidity parameters. n g is a unit norm tensor characterising the direction of the plastic flow and Φ(F) is an arbitrary function. The function chosen here has the following form: where N is a model parameter and F a function describing a convex surface in the stress space. The value F 0 corresponds to the value of stress below which no viscoplastic flow occurs. If F is chosen to be equal to τ and F 0 is chosen to be the cohesive-frictional strength s, then the rate of viscoplastic strain can be rewritten as ∂v where γ = 1/µ 1 /m , N = 1 /m and the elastic contributions are neglected. Assuming a simple shear Infinite Landslide Model where τ = ρg(h − z) sin θ and s = ρg(h − z) cos θ tan ϕ, then Equation (31) becomes which results in a linear velocity profile.
A new simple rheological law based on Perzyna viscoplasticity for frictional materials is then easy to find, being the expression of the basal friction τ b where s b = σ 3b tan ϕ is the basal shear strength (z = 0), µ is the viscosity [s] and v is the depth-integrated velocity. More details on the model can be found in [28][29][30].

Depth-Integrated SPH Model Coupled with a Finite Difference Scheme for Pore Pressure
The Smoothed Particle Hydrodynamics (SPH) is a meshless method, which has been applied to a large variety of problems, introduced independently by Lucy [39] and Gingold and Monaghan [40]. The SPH is a numerical technique able to simulate the propagation of fast landslides, which are treated as fluidised masses. It is based on the approximation of given properties and their spatial derivatives by integral approximation defined in terms of smoothed functions or kernel functions. An interpolation process calculates the relevant properties of each "particle" over neighbouring "particles". Therefore, the SPH is based on introducing a set of nodes together with a set of nodal variables. For landslides' problems, these variables are the height of the landslide at node I, the depth-averaged 2D velocity, the surface force vector at the bottom and the pore pressure at the basal surface. Details of the formulation can be found in [41].
In the present work, the previous development of the GEOFLOW-SPH code [11] has been enriched adding a 1D finite difference grid to each SPH node, in order to improve the description of the pore water evolution in the propagating mixture. The Finite Difference (FD) scheme chosen is explicit and is centred in space and forward in time (FCTS). At every node and time step, critical step times of FD and SPH are compared, and, if necessary, the pore pressure equation is solved using internal sub-steps.

Introduction
In this section, a depth-integrated model based on the v − p w Biot-Zienkiewicz formulation, enhanced with a diffusion-like equation to account for the pore pressure evolution within the soil mass, is applied to the Aberfan flowslide and Gypsum tailings impoundment flowslide that both occurred in 1966. In the case of the Gypsum flowslide, a pure cohesive viscous model-the Bingham model-is considered, while for the Aberfan flowslide, a frictional viscous model based on Perzyna viscoplasticity is developed.

East Texas Gypsum Tailings Failure (1966)
Tailings impoundments involve very loose materials. Failure of the dam results in the propagation of the fluidised materials, which behave like cohesive-viscous fluids. A representative case for which there is available information is that of East Texas Gypsum tailings impoundment, which failed in 1966. It has been described by Jeyapalan et al. [42], and Pastor et al. [12], who modelled the problem using a depth-integrated finite element model assuming that the material behaved as a Bingham fluid.
The purpose of this section is to show how the problem can be modelled using a depth-integrated SPH model.
According to the description provided in [42], the impoundment was rectangular; the tailings having reached a depth of 11 m at the time the failure took place. The failure affected a length of the dyke of 140 m. The material propagated some 300 m beyond the dyke before stopping. The average velocity was in the range of 2.5-5 m/s, with a propagation time close to 60-120 s.
One key point is the rheological model and its parameters. Here, we have used a Bingham model, with a yield stress of 750 Pa and a viscosity of 35 Pa·s, obtained by back analysis. The tailings were non-plastic silts, according to Jeyapalan, with D50 of 0.07 mm, a density of the particles of 2450 kg/m 3 and a density of the mixture of 1400 kg/m 3 .
We provide the height of the soil results (in meters) of the analysis in Figures 3 and 4, where we have depicted the propagation of the tailings at a series of time stations (t = 0, 30, 60, 80 and 120 s). The results agree well with the observations, the runout onto the plane being approximately 300 m and the movement being close to zero at time 60 s. The vertical scale has been enlarged by a factor of 10.
In this SPH simulation, we have observed that the failure propagated inside the impoundment for a distance longer than the 110 m provided in the bibliography.
Regarding the computation, we have used normalisation of the tailings height close to the dykes, in order to avoid the particle deficiency problem found close to boundaries. We have used 2485 SPH nodes in the analysis.

10.
In this SPH simulation, we have observed that the failure propagated inside the impoundment for a distance longer than the 110 m provided in the bibliography.
Regarding the computation, we have used normalisation of the tailings height close to the dykes, in order to avoid the particle deficiency problem found close to boundaries. We have used 2485 SPH nodes in the analysis.

The Aberfan Flowslide
Here, the event of the Aberfan flowslide will be analysed and the results of the simulations obtained using the mathematical and constitutive model described before will be shown. Moreover, a sensitivity analysis of the main parameters of the model is made in order to show how drastically the results change by changing the parameters model.
Aberfan is today a former coal mining village in South Wales (UK). In 1966, a flowslide of coal waste occurred, propagating onto the village itself and provoking 144 fatalities. Information about the failure mechanism and material properties have been provided by Bishop [43,44] and Hutchinson [14]. Other raw material is also available at the UK National Archive.
The Aberfan colliery waste was tipped on the side of a hill (Tip 7) facing the village. The triggering mechanisms of the flowslide lay in the hydrogeology of the site. Due to heavy rain, in fact, artesian pore pressure rose up in the sandstone beneath the less permeable glacial deposit at the toe of the slope, causing the liquefaction of the loose waste material dumped.
Tip 7 was about 67 m in height from the toe of the slope at the moment the slide occurred on 21 October and the underlying terrain had a slope of 12 degrees. The slide moved for 275 m before

The Aberfan Flowslide
Here, the event of the Aberfan flowslide will be analysed and the results of the simulations obtained using the mathematical and constitutive model described before will be shown. Moreover, a sensitivity analysis of the main parameters of the model is made in order to show how drastically the results change by changing the parameters model.
Aberfan is today a former coal mining village in South Wales (UK). In 1966, a flowslide of coal waste occurred, propagating onto the village itself and provoking 144 fatalities. Information about the failure mechanism and material properties have been provided by Bishop [43,44] and Hutchinson [14]. Other raw material is also available at the UK National Archive.
The Aberfan colliery waste was tipped on the side of a hill (Tip 7) facing the village. The triggering mechanisms of the flowslide lay in the hydrogeology of the site. Due to heavy rain, in fact, artesian pore pressure rose up in the sandstone beneath the less permeable glacial deposit at the toe of the slope, causing the liquefaction of the loose waste material dumped.
Tip 7 was about 67 m in height from the toe of the slope at the moment the slide occurred on 21 October and the underlying terrain had a slope of 12 degrees. The slide moved for 275 m before dividing itself into two lobes. The larger south lobe travelled for a distance of 500 m before impacting Aberfan buildings and stopped 100 m after, for a total propagation length of 600 m with estimated velocities in the range of 4.5 − 9 m s . Table 1 summarises the characteristics of the flowslide while in Figure 5 it is possible to see a photograph of the disaster.
So far, only simple 1D simulations of the Aberfan flowslide have been made as shown in the work of Pastor et al. [11,30]. In this paper, the authors want to present the 3D depth-integrated model of Aberfan flowslide, showing that the patterns observed are reproduced.
In order to do so, a proper 3D topographic mesh and an SPH mesh representing the initial mass is needed. The authors have built both of them by using topographic maps which are possible to consult in the UK National Archive. Figure 6 shows the topographic map that has been used in order to create the 3D topographical mesh (Figure 7).      The input for the topographic mesh is a Digital Terrain model with a Finite Element format with a total of 8733 nodes.
After identifying the edge of breakaway ( Figure 8) and its correspondent height of the sliding portion of Tip 7 that generated the flowslide, the SPH mesh has been properly created, as shown in Figure 9, with 1761 nodes with an average spacing of 3 m. It has been found that an average of 1700 nodes with a spacing no larger than 4 m, reproduces well the particular phenomena.
In Table 2, the parameters used to model the Aberfan flowslide which give the best agreement with field observations are presented. Erosion has been taken into account through the erosion coefficient of the Hungr erosion law [45]. In fact, with a careful reading of the report written immediately after the disaster and available from the UK National Archives, it is possible to see that erosion is widely mentioned by the author [43]. Moreover, rel w p represents the initial pore water pressure at the basal surface, varying between 0 and 1; 1 corresponding to liquefaction. Finally, the relative height of the basal saturated layer rel w h was assumed to be 0.4 times the total height of the flowslide at the beginning. The input for the topographic mesh is a Digital Terrain model with a Finite Element format with a total of 8733 nodes.
After identifying the edge of breakaway ( Figure 8) and its correspondent height of the sliding portion of Tip 7 that generated the flowslide, the SPH mesh has been properly created, as shown in Figure 9, with 1761 nodes with an average spacing of 3 m. It has been found that an average of 1700 nodes with a spacing no larger than 4 m, reproduces well the particular phenomena.
In Table 2, the parameters used to model the Aberfan flowslide which give the best agreement with field observations are presented. Erosion has been taken into account through the erosion coefficient of the Hungr erosion law [45]. In fact, with a careful reading of the report written immediately after the disaster and available from the UK National Archives, it is possible to see that erosion is widely mentioned by the author [43]. Moreover, p rel w represents the initial pore water pressure at the basal surface, varying between 0 and 1; 1 corresponding to liquefaction. Finally, the relative height of the basal saturated layer h rel w was assumed to be 0.4 times the total height of the flowslide at the beginning.

Introduction
In this section, the results obtained with the setup described in the previous sections and the model parameters of Table 2 will be shown. Moreover, in order to prove the validity of the model itself and to justify the used parameters, a parametric study is conducted. This parametric study has the objective of showing the non-negligible differences that arise when it comes to choosing the right parameters which allow the description of the phenomena. In particular, the results shown and considered of most importance by the authors are: • the height of propagation • the propagation profile In order to do so, the parameters that will vary in the simulations will be the angle of friction ϕ and viscosity factor µ, specific to the Perzyna-based model, and the erosion rate of growth.

Simulation Results
In Figure 10, the pore pressure contours evolution is presented at 0, 2, 5 and 10 s. Please note that in order to improve readability, the saturated layer has been expanded and now it occupies the whole mass. This is possible because we are considering the depth of the basal saturated layer proportional to the one of the landslide.

Introduction
In this section, the results obtained with the setup described in the previous sections and the model parameters of Table 2 will be shown. Moreover, in order to prove the validity of the model itself and to justify the used parameters, a parametric study is conducted. This parametric study has the objective of showing the non-negligible differences that arise when it comes to choosing the right parameters which allow the description of the phenomena. In particular, the results shown and considered of most importance by the authors are:  the height of propagation  the propagation profile In order to do so, the parameters that will vary in the simulations will be the angle of friction  and viscosity factor  , specific to the Perzyna-based model, and the erosion rate of growth.

Simulation Results
In Figure 10, the pore pressure contours evolution is presented at 0, 2, 5 and 10 s. Please note that in order to improve readability, the saturated layer has been expanded and now it occupies the whole mass. This is possible because we are considering the depth of the basal saturated layer proportional to the one of the landslide. The results of the propagation and height of the soil obtained with the parameters described in Table 2 are shown in Figure 11 at time 5, 10, 25, 35 and 50 s. Results satisfactorily reproduce the flowslide. The legend in the picture refers to the height of the soil in meters at 50 s. It is possible to note that the final height of the soil of the left lobe at 50 s is almost 10 m which matches with the real height reported in [43]. Furthermore, the SPH program reproduces well the division of the flowslide The results of the propagation and height of the soil obtained with the parameters described in Table 2 are shown in Figure 11 at time 5, 10, 25, 35 and 50 s. Results satisfactorily reproduce the flowslide. The legend in the picture refers to the height of the soil in meters at 50 s. It is possible to note that the final height of the soil of the left lobe at 50 s is almost 10 m which matches with the real height reported in [43]. Furthermore, the SPH program reproduces well the division of the flowslide into two lobes. Results of the soil height also match well with the one-dimensional results obtained by Pastor et al. [11,30]. into two lobes. Results of the soil height also match well with the one-dimensional results obtained by Pastor et al. [11,30].  Table 2. Colors refer to the height of the soil.

Parametric Study
After showing the results of the simulation, a parametric study is now shown in order to understand the difference that other sets of parameters might have on the simulation results. Table 3 shows the parameters chosen for every simulation for a total of six simulations to compare, including the original one displayed above. At first, the difference in the soil height is shown. Profiles are displayed at 10 s, 20 s, 30 s and 50 s in Figures 12 and 13. In Figure 6, line A-A' shows the profiles' perspective. In order to improve the results' readability, Table 4 summarises the values of the maximum height reached at every stage.   Table 2. Colors refer to the height of the soil.

Parametric Study
After showing the results of the simulation, a parametric study is now shown in order to understand the difference that other sets of parameters might have on the simulation results. Table 3 shows the parameters chosen for every simulation for a total of six simulations to compare, including the original one displayed above. At first, the difference in the soil height is shown. Profiles are displayed at 10 s, 20 s, 30 s and 50 s in Figures 12 and 13. In Figure 6, line A-A' shows the profiles' perspective. In order to improve the results' readability, Table 4 summarises the values of the maximum height reached at every stage.   Table 3.  Table 3.  Table 3.  Table 3.
The height of the flowslide at the end of its propagation was observed to be 10 m when hitting the school building of the village [43]. It can be observed that by changing the parameters of the model (Table 3), the final height varies substantially. In the case of a higher viscosity value (N.2), zero erosion (N.4) or higher friction angle (N.6), the maximum height at the end of the simulation is far lower than what it is supposed to be. On the other hand, by lowering the viscosity value (N.3) or lowering the friction angle (N.5), the height reached exceeds the reference value. In the case of N.3, the 10 m height is reached at 45 s, while for the N.5 case, the 10 m height is reached at 42 s. Another clear observation that one can find in [43] is the fact that the erosion takes place during the hazard in the form of a channel excavated during the flowslide. By not taking into account the erosion rate, the model does not successfully reproduce the same pattern, as visible in simulation N.4.
Moreover, it is also interesting to see how the propagation profiles change with respect to the different parameters used as in Figures 14 and 15.

Conclusions
A depth-integrated model, based on thew v p Biot-Zienkiewicz formulation enhanced with a diffusion-like equation to account for the pore pressure evolution, is presented. This paper clarified the influence of the selection of the viscous model and proposed a rheological model based on Perzyna viscoplasticity. In a landslide, the pore pressure-shear stress interaction cannot take place when a pure cohesive viscous model is used. The difference between pure cohesive and frictional rheological behaviour is discussed. This approach, which couples solid skeleton and pore fluid, has been applied in order to simulate two case studies: the Aberfan flowslide and Gypsum tailings impoundment flowslide that

Conclusions
A depth-integrated model, based on thew v p Biot-Zienkiewicz formulation enhanced with a diffusion-like equation to account for the pore pressure evolution, is presented. This paper clarified the influence of the selection of the viscous model and proposed a rheological model based on Perzyna viscoplasticity. In a landslide, the pore pressure-shear stress interaction cannot take place when a pure cohesive viscous model is used. The difference between pure cohesive and frictional rheological behaviour is discussed. This approach, which couples solid skeleton and pore fluid, has been applied in order to simulate two case studies: the Aberfan flowslide and Gypsum tailings impoundment flowslide that For simulation N.2, N.4 and N.6, it is clearly visible that the propagation profile does not reach its extension either in longitude (N.4) or width (N.2). For simulation N.3 and N.5, the profiles are much more similar to the original, especially the N.3, but they tend to exceed the propagation extension (N.5). Moreover, as comparable in Figure 5, the propagation profile is somehow homogeneous, while for both N.3 and N.5, several small branches are visible which do not replicate the real hazard or the division into two lobes as well as the parameters used for the original simulation.

Conclusions
A depth-integrated model, based on the v − p w Biot-Zienkiewicz formulation enhanced with a diffusion-like equation to account for the pore pressure evolution, is presented. This paper clarified the influence of the selection of the viscous model and proposed a rheological model based on Perzyna viscoplasticity. In a landslide, the pore pressure-shear stress interaction cannot take place when a pure cohesive viscous model is used. The difference between pure cohesive and frictional rheological behaviour is discussed.
This approach, which couples solid skeleton and pore fluid, has been applied in order to simulate two case studies: the Aberfan flowslide and Gypsum tailings impoundment flowslide that both occurred in 1966. By combining the SPH technique with a set of Finite Differences, it has been possible to gain insight into the pore water pressure developed during the hazard due to changes in height, vertical consolidation and changes in total stresses.
The application of the methodology proposed and the results obtained show its suitability to be applied in studies of the propagation phase of fast landslides. The results of the Gypsum tailings impoundment flowslide agree well with the observations reported in terms of runout. Aberfan flowslide turns out to be an excellent example where the trajectory of the flowslide matches well with that which occurred, especially when the final bifurcation takes place. Moreover, the final height of the simulation also matches with the one reported in [43] and satisfactory profiles of the pore water pressure evolution varying through time are reached.