WCSPH with Limiting Viscosity for Modeling Landslide Hazard at the Slopes of Artificial Reservoir

Abstract

This work illustrated an application of the FOSS code SPHERA v.8.0 (RSE SpA, Milano, Italy) to the simulation of landslide hazard at the slope of a water basin. SPHERA is based on the weakly compressible SPH method (WCSPH) and holds a mixture model, consistent with the packing limit of the Kinetic Theory of Granular Flow (KTGF), which was previously tested for simulating two-phase free-surface rapid flows involving water-sediment interaction. In this study a limiting viscosity parameter was implemented in the previous formulation of the mixture model to limit the growth of the apparent viscosity, thus saving computational time while preserving the solution accuracy. This approach is consistent with the experimental behavior of high polymer solutions for which an almost constant value of viscosity may be approached at very low deformation rates near the transition zone of elastic–plastic regime. In this application, the limiting viscosity was used as a numerical parameter for optimization of the computation. Some preliminary tests were performed by simulating a 2D erosional dam break, proving that a proper selection of the limiting viscosity leads to a considerable drop of the computational time without altering significantly the numerical solution. SPHERA was then validated by simulating a 2D scale experiment reproducing the early phase of the Vajont landslide when a tsunami wave was generated that climbed the opposite mountain side with a maximum run-up of about 270 m. The obtained maximum run-up was very close to the experimental result. Influence of saturation of the landslide material below the still water level was also accounted, showing that the landslide dynamics can be better represented and the wave run-up can be properly estimated.

1. Introduction

Water represents a key factor in the mechanics of failure and in the dynamics of sloping rock–soil masses. In particular, intense rainfall events can be directly or indirectly involved in the trigger and run-out of landslides.

Rainfall events cause water infiltration at slopes, leading to an increase of both volumetric water content and pore pressure in the shallow soil layers that worsen the slope stability and may directly activate landslide movement [1].

Theoretical models are available to predict the influence of hydrological processes on the slope stability [2]. To evaluate landslide susceptibility, reliable prediction of groundwater pressure evolution during any meteorological event is fundamental: the simplest models assume steady or quasi-steady modification of the groundwater table with slope-parallel groundwater flow [3,4], while more advanced formulations include pressure change in the normal direction to the hill-slope caused by transient rainfall [5,6,7].

Beside the stability analysis for evaluating landslide potential, the dynamic analysis of the activated landslide is of great concern for defining related hazard.

Focusing on the important category of shallow landslides, triggered by heavy rainfalls and characterized by low displacement rates at a narrow shear zone just below the ground surface [8], their dynamics can be effectively described with finite element approach using few constitutive parameters that can be estimated through conventional geotechnical tests [9]. Preliminary analysis of shallow landslides can be performed adopting a simple sliding block model and relating the landslide movement directly to rain recordings, thus allowing future scenarios to be predicted from expected rainfall scenarios [10,11] and appropriate measures to be defined for slope stabilization [12].

Recent climate trends affecting rainfall [13] should be considered for proper definition of the rainfall characteristics as input for the above-mentioned models to obtain reliable estimate of landslide susceptibility.

In addition, rainfall events may indirectly trigger the landslide because intense rainfall can cause significant and rapid level rise of the water in artificial reservoirs with consequent seepage into the basin’s slopes and rise of groundwater surface. Subsequent rapid drawdown of the reservoir inverts the groundwater flow leading to high pore water pressure in the saturated natural slopes. The change in water content and pore pressure, in combination with the geotechnical properties of involved soil, affects the shear strength leading to possible landslides.

The interaction of the falling soil mass and the stored water induces formation of dangerous impulsive waves, usually referred to as a class of tsunami waves similar to open sea tsunami [14]. Empirical formulations relating the principal features of this class of waves with the landslide characteristics were defined in an experimental work [15] and applied to real cases, including the Vajont landslide that occurred in 1963, providing a reliable estimate of the order of magnitude of the expected consequences in terms of maximum wave height and wave run-up [16].

During the post-failure phase, sophisticated predictive models are required for reliable numerical investigation of the hazard associated to landslide run-out and its interaction with water in a complex scenario. These models should account for those physical aspects influencing fast dynamics, including the landslide water content affecting both sliding velocity and generated wave.

Moreover, complex landslides involve uncertainties related, for instance to angle of internal friction, pore water, and failure surface, which may significantly affect the model results [17].

Despite their intrinsic limitations, deterministic models represented by fast-running (in computational time) numerical models may be helpful in such a way they allow supporting risk analysis and assessment to address uncertainty inherent in landslide hazards [18]. By creating probabilistic maps of risk variables, a multi-disciplinary decision support system (DSS) for natural hazard risk reduction (NHRR) and management could be developed to explore many potential risk-reduction options with the help of numerical modeling [19,20].

The Smoothed Particle Hydrodynamics (SPH) method can be effectively adopted to set up a reliable numerical model for simulating the fast dynamics of landslide generated wave in artificial reservoir.

This Lagrangian meshless method proved to be successful for the engineering analysis of rapidly varied free-surface flows involving spillway hydraulics [21], hydraulic jump oscillations at an abrupt drop [22], dam break flow hazards in complex topography [23], wave propagation and impact pressure on rigid structures [24,25], and solid body transport in free surface flows [26].

Furthermore, SPH showed great potential in simulating multiphase problems, including non-cohesive sediment interacting with rapidly varied free-surface water flows with moving interfaces and large displacements: scouring behind seawall due to continuous tsunami overflow [27], impulsive dynamics of submerged non cohesive sediments involving buried explosion [28], scour over movable bed dam break flow [29] and downstream fluvial structures [30], shock wave propagation [31], surge waves generated by deformable landslides [32] and submarine deformable landslides [33]. SPH was successfully applied to the three-dimensional simulation of the Vajont landslide which was assumed to move as rigid body [34].

In this paper, an application of the FOSS software called SPHERA v.8.0 (RSE SpA) [35] that is based on the weakly compressible SPH method (WCSPH) is shown. SPHERA holds a mixture model which is consistent with the packing limit of the Kinetic Theory of Granular Flow. This model was previously tested for simulating erosional dam breaks [36].

In this work, SPHERA was applied to model the coupled fast dynamics of non-cohesive sediment and water. Even if this field is rather explored with SPH method [37,38,39], this work represented the first attempt to apply the Kinetic Theory of Granular Flow to simulate, at the slope of an artificial basin, the fast dynamics of a complex landslide undergoing large deformation and involving a sediment volume comparable to the stored water volume.

The model includes the effect of interstitial water in the saturated portion of the landslide below the still water level: as expected, pore water influences the rapid landslide movement and, hence, the generated wave and its maximum run-up.

As a new feature, this work presented a modified rheological model where the limiting viscosity, early introduced to account for the experimental behavior of high polymer solutions approaching almost constant value of viscosity at low deformation rates [40], is implemented for the first time in the reference mixture model [36] and is used as a numerical parameter for results optimization. By assuming a suitable value for the reference permitted error (with respect to reference data), a significant reduction of the computational time was obtained for the reference simulation, thus providing a fast-running numerical model for possible coupling with a decision support system with the aim of exploring many potential landslide risk reduction options.

Some preliminary tests were performed by simulating 2D dam break flow over movable bed; these tests permitted identifying the reference simulation and showed that the limiting viscosity allows obtaining significant reduction of the computational time with negligible effects on the numerical solution.

It is worth noting that the number of intermediate simulations necessary to configure the reference simulation depends on many variables and is not involved in the estimate of the computational time. Further, once the optimum limiting viscosity is assessed for a fixed test case, all the variants of the same test case do not need any further optimization and can refer to the limiting viscosity value previously found.

Furthermore, once a sufficient number of test cases is addressed (on the same application field), it will be possible to suggest narrower interval of values for the ratio between the limiting viscosity and the maximum viscosity (these values, by definition, must be neither higher than unity nor lower than zero, in any case). From this point of view, the optimization procedure for the limiting viscosity is analogous to the optimization procedure for the Courant–Friedrichs–Lewy number (CFL): only the optimum value of the CFL is considered when assessing the computational time, neglecting the overall time resulting from all the simulations necessary to assess the optimum CFL.

Subsequently, SPHERA was validated with a test case simulating the coupled dynamics of a deformable landslide, occurring at the slope of an artificial basin, and the stored water. Validation of SPHERA was carried out by simulating the 2D scale experiment carried out in 1968 at the Hydraulic Laboratory of Padua University. A scale physical model was built reproducing the characteristic cross section of the Vajont basin close to the dam to investigate the early phase of the catastrophic event that occurred in 1963 when a landslide estimated in the order of 270 million cubic meters fell into the Vajont artificial basin located in the northeastern part of Italy. The stored volume at the date of the disaster was about 115 million cubic meters with water level at 700 m above mean sea level.

Padua experiment aimed at assessing the relation between the landslide run-out velocity and the maximum run-up of the generated wave climbing on the opposite slope of the basin.

To provide an idea of this catastrophic event, it should be known that the landslide was generated by the sudden collapse of an enormous portion of the northern side of the Toc Mountain defining the left-hand slope of the artificial lake near the dam.

The landslide slid into the hydroelectric basin with a run-out duration that was estimated from seismic recordings in about 20–25 s before the landslide front came to a standstill after impacting against the opposite mountain side.

The peculiar kinematics of the landslide, related to the complex three-dimensional nature of the involved geological structures [41], and the resulting high amount of kinetic energy of the detached soil mass were both responsible for a huge water wave to be generated that in the early phase climbed on the opposite mountain side with a maximum run-up of about 270 m above the dam crowning [42].

In the subsequent descending phase part of the rundown wave covered the landslide body where the residual lake of Massalezza was created in the central area; another portion of the wave propagated upstream the Vajont valley toward Erto village; and the remaining part of the wave overtopped the dam and flooded Longarone and some downstream villages. About 2000 casualties were caused by the disaster.

The catastrophic effect of the Vajont landslide, related to the height of the generated tsunami, was drastically underestimated before the tragic event for the reasons explained in the following.

A 3D experimental campaign was carried out before the disaster occured by realizing the scale physical model of the basin and landslide at the hydraulic research center of Nove (Italy).

During the experimental campaign, various catastrophic scenarios were considered, involving different run-out velocities and falling modalities with the aim of quantifying the actual risk depending on the landslide shape and velocity in relation with the level of the stored water.

According to the information collected in the early studies, it was stated that the falling time to be adopted in the experiments should have been ≥60 s; therefore, the maximum landslide velocity was strongly underestimated because the velocity associated to a run-out duration of 20 s results three times the maximum velocity adopted in the early experiments.

Consequently, the maximum landslide kinetic energy considered in the Nove experiments was about 1/9 of the actual one, leading to underestimate the maximum wave height which was found to be in the worst condition (i.e., run-out duration of 60 s) of about 27 m with stored water level at 700 m above mean sea level (a.m.s.l.).

Thus, the early experimental campaign of Nove led to the wrong conclusion that a stored water level of 700 m a.m.s.l. could be considered as a safe limit with respect to the worst scenario that could be taken into consideration a priori. Unfortunately, this estimate was disavowed by the facts.

The post-event two-dimensional experimental campaign carried out in 1968 at the Hydraulic Laboratory of Padua University aimed at demonstrating that the run-out duration played a key role in the wave generation and that the early experiments before the catastrophic event: (i) were properly conceived; (ii) were able to reproduce the relevant features of the Vajont landslide on a scale model using non-cohesive granular material with a controlled sliding movement; and (iii) would have provided more reliable prediction of the hazard if proper run-out time, closer to the actual one, had been considered.

The experimental model of Padua was reproduced with SPHERA: landslide run-out and maximum wave run-up were compared with experimental data showing good agreement. Influence of saturation on the landslide dynamics was also evaluated.

2. Numerical Aspects and Mathematical Details

The Smoothed Particle Hydrodynamics (SPH) technique was initially developed for astrophysical purposes [43] and was subsequently extended to the simulation of hydrodynamic problems [44]. SPH belongs to the category of numerical methods called meshfree particle methods [45] that are suitable to simulate large deformations of the material continuum without degradation of the numerical result due to mesh distortion.

Following a Lagrangian approach, any fluid continuum is discretized through a finite number of material points or particles with constant mass, whose dynamics respond to the conservation equations of the classical physics [46]. In contrast to traditional grid-based methods, these particles are not connected by any topological grid (or mesh). Consequently, the numerical approximation of the partial differential governing equations is performed through a weighting (or smoothing) function as described in the following.

A generic physical quantity f_i (e.g., velocity or density) can be estimated at a given continuum particle, denoted by index i, based on the values assumed at the N neighboring particles, each denoted by index j, according to the following discrete formulation (particle approximation):

$$\langle f_i\rangle ={\displaystyle \sum _{j=1}^N\frac{m_j}{{\rho }_j}{f}_j{W}_{ij,h}}$$

(1)

In Equation (1), m_j and ρ_j are the mass and density of the j-th particle. The neighboring particles are contained within the interaction domain that is centered on particle i and has a circular or spherical shape according to the problem dimension. The interaction domain represents the compact support of the kernel function W_ij,h used to carry out the weighted estimate. The kernel is a central weighting function whose radius is proportional to the so-called smoothing length h. The value assumed by the kernel depends on the relative distance between particles i and j.

The particle approximation of a function derivative can be obtained by replacing f_i with 99 $\nabla f_i$ in Equation (1) and performing some manipulations [47] leading to:

$$\langle \nabla f_i\rangle ={\displaystyle \sum _{j=1}^N\frac{m_j}{{\rho }_j}{f}_j\nabla W_{ij,h}}$$

(2)

The approximation of a function derivative in Equation (2) is valid when the interaction domain is entirely contained within the inner domain of the discretized continuum. Near its boundary, some corrections should be introduced as kernel truncation occurs [48].

The SPH particle approximations in Equations (1) and (2) represent the basic formulations that should be properly corrected to improve conservation and consistency properties [46].

The basic concepts on the representative volume for multiphase mixture models are described in [49]. From a microscopic point of view, a porous medium is characterized by the presence of a solid matrix and void spaces which are distributed throughout it; the void spaces may be occupied by some fluid phases separated from each other by an interface. The mathematical description of the dynamics in porous medium subjected to some imposed excitation is translated into a macroscopic level where the (material) point identifies a discretization volume that contains, in general, both a solid phase and void spaces filled by some fluid phases. This macroscopic description implies that each phase (fluid or solid) is represented by a material continuum and that different continua may overlap at the material points within the computational domain. The values of the relevant field variables at every point of these continua represent the average of the corresponding phase taken over the discretization volume centered at that point.

It must be pointed out that the mixture model was adopted here for simulating the fast dynamics of a landslide in the post failure phase. For this reason, the motion of the saturated solid matrix undergoing large deformation is treated following the fluidic approach [17]. Therefore, in subsequent development of the balance equations, both the granular (solid) material and the water that may fill the grain–grain interspaces are considered as two distinct and immiscible components.

2.1. Mixture Model for Dense Granular Flows

The mathematical model and the SPH numerical model presented in [36] are used in this study as the reference mixture model for dense granular flows (i.e., the “packing limit” of the Kinetic Theory of Granular Flow (KTGF)). The reference mixture model represents both fluid and solid granular material and is synthesized in the following.

The SPH model presented in the following was based on the assumption that the involved fluids are slightly compressible. Therefore, pressure represents an explicit thermodynamic variable that is computed solving an equation of state. Such a weakly compressible approach (WCSPH) offers some advantages concerning, for instance, the relative simplicity of the solution algorithm, but suffers from some limitations due to the introduction of artificial speed of sound leading to high-frequency pressure oscillations and time step restrictions. An alternative approach is based on the incompressible formulation (ISPH) that was successfully applied to the analysis of both Newtonian and non-Newtonian free-surface flows [50], as well as viscosity-dominated flows [51].

According to the weakly compressible SPH approach, some alternative expressions for the function derivative of Equation (2) may be derived [45]. Based on these expressions, the Lagrangian form of the governing equations of the fluid dynamics can be discretized. The following discretized form of the mass and momentum balance equations for a weakly compressible fluid (WCP) are used:

$$\begin{array}{l}\langle \frac{d{\rho }_i}{dt}\rangle =-{\rho }_i{\displaystyle \sum _{j=1}^N\frac{m_j}{{\rho }_j}\left(u_j-u_i\right)\cdot \nabla W_{ij,h}}+{{\rm B}}_{\rho }\\ \langle \frac{d{u}_i}{dt}\rangle =-\frac{1}{{\rho }_i^{}}{\displaystyle \sum _{j=1}^N\frac{m_j}{{\rho }_j}\left(p_i+p_j+\frac{{\Pi }_{ij}}{{\rho }_i^{}}\right)\nabla W_{ij}}+\frac{\nabla \cdot {\mathsf{\tau }}_i}{{\rho }_i}+g+{{\rm B}}_u\end{array}$$

(3)

In Equations (3), u is the mixture velocity vector, p is the mixture pressure of a particle, g is the gravity acceleration and ρ represents the mixture density. The term Π_ij denotes the well-known Monaghan artificial viscosity [52] that is used exclusively for the fluid phase.

The symbol τ denotes the mixture deviator viscous stress tensor that is considered in the momentum balance of the granular phase solely; its definition is provided in the following.

The terms B_ρ and B_u denote the boundary contributions to be included in the mass and in the momentum balance Equations (3) of those particles that are close to the boundary, respectively. The semi-analytic approach [24] is adopted to formulate these terms of boundary contributions. Therefore, these two terms represent integral contributions extended to the portion of the boundary that intersects the particle’s interaction domain; it is assumed that this portion of the boundary is continuously filled by fluid of the same phase as the considered particle and therefore treated as a material continuum with suitable distribution of velocity, density or pressure. The integral terms B_ρ and B_u can be computed analytically and their value depends only on the problem geometry, as explained in detail in [24].

The sound speed in both water and granular phase were selected to be at least 10 times the speed of bulk flow because this assures relative density variations of at most 1%. Therefore, the adopted state equation for barotropic flows was obtained by linear approximation to the well-known Tait’s equation about the reference density [24].

Time integration scheme followed a second-order Leapfrog scheme where the velocity of each particle is calculated at mid time-step with respect to both position and density. Details about this algorithm are provided in [26]. The numerical stability is based on the concepts discussed in [53] and the time step is selected in accordance to the following criterion:

$$dt={\min }_0\left\{CF{L}_{vis}\frac{2h^2}{\nu };CFL\frac{2h}{c+\left|u\right|}\right\}$$

(4)

The coefficient CFL_vis is set to 0.05, the CFL number is equal to 0.1 in the simulations of Section 3, v = µ/ρ is the kinematic viscosity of the mixture, and c represents the modulus of the highest sound speed among the values assumed by the particles all over the domain.

In the subsequent mathematical notation, the superscript α is introduced to distinguish those field variables that are representative of the liquid component when α = f, and those that belong to the solid component when α = s.

For every component of the fluid phase, the corresponding volume fraction ε^α is defined, along with its velocity u^α and the density ρ^α of pure component.

The mixture density ρ is defined as the sum of concentrations c^α = ε^αρ^α of both components:

$$\rho ={\epsilon }^f{\rho }^f+{\epsilon }^s{\rho }^s$$

(5)

while u is the mixture (or mass weighted) velocity given by the linear combination of the velocity u^α of each component with coefficients given by the corresponding mass fraction ε^αρ^α/ρ:

$$u=\frac{{\epsilon }^f{\rho }^f{u}^f+{\epsilon }^s{\rho }^s{u}^s}{\rho }$$

(6)

The Lagrangian formulation of mass balance for the mixture assuming a Weakly Compressible approach reads:

$$\frac{d\rho }{dt}=-\rho \nabla \cdot u$$

(7)

The Lagrangian formulation of the linear momentum balance for the mixture is given by:

$$\frac{d\left(\rho u\right)}{dt}=-\nabla \left(p^f+p^s\right)+\nabla \cdot \left({\mathsf{\tau }}^f+{\mathsf{\tau }}^s\right)+\rho g$$

(8)

According to [54], the pressure of the solid component is defined as the effective mean normal stress:

$$p^s=\frac{\mathrm{tr}{\mathsf{\sigma }}^{\prime }}{3}$$

(9)

where σ’ is the effective spherical stress tensor. The shear stress divergence for the solid component in Equation (8) can be expressed as a function of the strain rate tensor D through the apparent viscosity µ_fr:

$$\begin{array}{l}\nabla \cdot {\mathsf{\tau }}^s=\nabla \cdot \left(2{\mu }_{fr}D\right)\\ {\mu }_{fr}=\frac{{\sigma }^{\prime }\sin \varphi }{2\sqrt{I{I}_D}}\end{array}$$

(10)

In Equation (10) the symbol II_D is the second invariant of the strain rate tensor D, and ϕ denotes the angle of internal friction of the solid component.

For the Newtonian liquid component, the shear stress divergence in Equation (8) is given by:

$$\nabla \cdot {\mathsf{\tau }}^f={\epsilon }^f\left[\frac{{\mu }_f}{3}\nabla \left(\nabla \cdot u\right)+{\mu }_f{\nabla }^{2}u\right]$$

(11)

The mixture viscosity µ is defined as follows:

$$\mu ={\epsilon }^f{\mu }_f+H_{\left({\epsilon }^s-{\epsilon }^{s,p}\right)}{\mu }_{fr}$$

(12)

where H is the Heaviside step function and ε^s,p is the volume fraction in the packing limit.

Balance Equations (7) and (8) can be discretized based on Equations (3), as explained in detail in [36].

As can be seen from the second part of Equation (10), the apparent viscosity of the granular material µ_fr depends on the inverse of the second invariant II_D of the strain rate tensor; therefore, the apparent viscosity decreases as the shear rate grows, whereas µ_fr increases indefinitely as the rate of shear approaches zero, similar to a pseudoplastic fluid (blue dashed curve in Figure 1). This behavior affects the mixture viscosity µ in Equation (12) and, hence, the size of the time step dt that depends on the inverse of µ through Equation (4). In particular, at shear rates close to zero, the mixture viscosity assumes very high values and the time step needed to assure numerical stability becomes very small.

In this model, the unbounded growth of µ_fr is prevented by the threshold µ_thr for the mixture viscosity that has a physical meaning and acts when approaching zero shear rate. Mixture particles with a viscosity µ higher than µ_thr are considered in the elastic–plastic regime of soil deformation where the frictional regime of the packing limit in the KTGF does not apply; therefore, the threshold viscosity µ_thr is imposed to these particles that are excluded from the SPH computation and are held fixed with a pressure field derived from the calculated pressure of the neighboring mixture particles whose viscosity is below µ_thr. As a numerical consequence, the threshold viscosity µ_thr reduces the computational time.

The value of threshold viscosity should be properly selected so that it does not influence the model results. This goal is obtained through a convergence procedure: during this procedure subsequent simulations are performed by progressively increasing the threshold viscosity until the numerical solution converges because the only solid particles that should be excluded from the computation are those ones whose rate of deformation is so small that it does not affect significantly the global dynamics of the mixture.

The threshold viscosity holds a physical meaning as it mimics a sharp transition between (static) elastic–plastic regime and the (dynamic) frictional (or plastic–frictional) regime, similar to a Bingham fluid that exhibits a yield stress. In the present model, the apparent viscosity was not a constant in the frictional regime, similar to a pseudoplastic fluid [55].

2.2. Implementation of the Limiting Viscosity

The actual behavior of pseudoplastic fluids shows that an almost constant value of the apparent viscosity may be approached at very low deformation rates, as in the case of experimental tests with high polymer solutions [40]. This constant value is generally referred to as limiting viscosity and it is denoted by µ₀.

In [50], the non-Newtonian fluid has been simulated with Cross model that includes a similar rheological parameter µ₀ denoting the viscosity at very low shear rate. This constant parameter µ₀, which may be related to the Bingham yield stress under suitable condition, allows obtaining the gradual transition from the flow condition to the plug flow or no flow at all, and vice versa. Instead, adopting Bingham model, an abrupt transition occurs owing to the presence of the yield stress that induces the unbounded grow of the apparent viscosity when the shear rate becomes infinitesimal. Therefore, in the Cross model, the apparent viscosity varies continuously with the shear rate and approaches a constant value as the shear rate tends to be infinitesimal, thus preventing numerical divergence. In [50], the parameter µ₀ has been set equal to 10³ µ_∞, where µ_∞ is the viscosity at very high shear rates. In the present model, instead, the limiting viscosity parameter µ₀ was adopted as a parameter for numerical optimization.

Starting from the mixture model of Section 2.1, the limiting viscosity is presented hereafter as a new numerical feature of the reference numerical model and it is implemented in SPHERA [35]. It should be noticed that, starting from the physical meaning of the limiting viscosity, in this model, it was used as a numerical parameter to reduce the computational time, as illustrated in the following. In some SPH applications dealing with landslide run-out [17] and sediment scouring [56], a similar viscosity parameter was successfully adopted working as a numerical threshold.

The numerical role of the limiting viscosity is shown by the red curve in Figure 1: µ₀ acts at low deformation rates within the frictional regime, which is characterized by the frictional regime of the packing limit in the KTGF, by keeping constant the mixture viscosity and, hence, reducing the computational time.

The limiting viscosity µ₀ represents a numerical parameter that should be evaluated for the investigated problem to benefit from the reduction of the computational time while maintaining sufficiently accurate and reliable results. Therefore, some preliminary runs have been carried out, as illustrated in Section 3.1, by considering a simplified representative test. When performing the convergence procedure of the threshold viscosity µ_thr explained above, the limiting viscosity was deactivated by keeping µ₀ > µ_thr.

Once it was detected, the lower value µ_thr of the threshold viscosity above which the influence on the simulation results is negligible, a second set of simulations was carried out by reducing progressively the limiting viscosity below the value µ_thr: in this way, it is possible to detect the lower value µ₀^min of the limiting viscosity which matches the imposed numerical precision (i.e., allows obtaining the reference permitted error with respect to reference results). Setting µ₀ = µ₀^min assured a reliable compromise between model accuracy and computational time.

The above-mentioned procedure for the optimization of numerical results through the limiting viscosity was only carried out for the reference spatial resolution dx. In the case where it could be interesting to optimize the limiting viscosity for coarser resolutions, it is shown that µ₀ should increase with increasing spatial resolution to obtain the reference permitted error: this is because part of the error is already due to the increased spatial resolution.

3. Numerical Results

This section describes the simulations of some mixture multiphase free-surface flows carried out with SPHERA v 8.0. All the simulations performed concern the analysis of the fast dynamics of free surface flows involving the interaction between the water and non-cohesive sediment.

Section 3.1 describes the early tests carried out for studying the convergence of results by varying the threshold viscosity µ_thr and the subsequent calibration for the limiting viscosity µ₀ to find the minimum value µ₀^min assuring a suitable agreement between accuracy and computational cost. To reach such a goal, a simplified 2D erosive dam break test is adopted simulating the scouring induced by a shallow wave that is originated by the sudden collapse of a water column and that propagates over a thin bed of non-cohesive sediment.

Section 3.2 shows the validation of the reference numerical model with the new feature of limiting viscosity by simulating the 2D experiment that was set up at the Hydraulic Laboratory of the Padua University in 1968 to perform the post-event analysis of the huge landslide occurred on the left-hand slope of the Vajont artificial basin in 1963.

3.1. Two-Dimensional Erosive Dam Break

The two-dimensional erosive dam break of a water column collapsing over an erodible bed has been performed with SPHERA to evaluate the threshold viscosity and to assess the optimum value for the limiting viscosity.

The problem is illustrated in Figure 2 showing on the left-hand panel the water and solid bed distribution at the initial time; this demonstrative test represents a modification of an erosional dam break tutorial available in SPHERA v8.0 package [35] and shares some relevant features with similar laboratory tests that have been adopted for numerical model validation, including weakly compressible SPH [57,58] and incompressible SPH [29].

In the performed test, the initial depth of the water column is much greater than the thickness of the bed sediment, resulting in a higher and faster excavation at the wave front, as shown by the representative frames on the right-hand side of Figure 2. Moreover, the horizontal dimension of the computational domain is rather limited and therefore the generated wave suddenly impacts against the downstream wall originating a vertical jet.

The initial length and height of the water column are L_w = 1.50 m and H_w = 1.30 m, respectively. The density and bulk modulus of the liquid component are ρ^f = 1000 kg/m³ and K^f = 1.30 × 10⁶ Pa. The non-cohesive granular bed at initial time has uniform thickness of H_s = 0.20 m and the longitudinal length is L_s = 3.00 m. The geotechnical parameters of the solid component are: median diameter d₅₀ = 1.0 mm, density ρ^s = 1500 kg/m³, bulk modulus K^s = 1.95 × 10⁶ Pa, angle of internal friction ϕ = 20° and porosity ε^f = 0.5. The saturated mixture density is obtained through Equation (5) ρ = 1250 kg/m³.

The adopted spatial resolution is dx = 0.005 m while the smoothing length is assumed h = 1.3 dx. The total number of particles in all dam break simulation is 11,300. The coefficient of artificial viscosity is set at α_M = 0.1, while the speed of sound for both components of the mixture is c_s = 36.1 m/s.

The water column, initially at rest, is confined on the right-hand side by the vertical wall of the tank, while the vertical left-hand boundary of the water is free to collapse. A shallow water wave propagates from the toe of the column toward the left-hand vertical wall of the tank; the mean effective stress represents the main cause of mobilization of the solid grains (as observed also in [36]) and formation of a bed-load transport layer that moves along with the water wave. The shape of the bed-load transport layer is such that a sloping surface develops below the front of the water wave that is therefore pushed upward. A concentric plunging wave, made of an outer water coating and an inner sediment layer, develops and subsequently impacts the vertical downstream wall at about t = 0.29 s. After the impact, a vertical water jet climbs on the left-hand wall of the tank while a sediment wedge accumulates at the left corner below the water jet.

The first set of simulations has been carried out to evaluate the convergence of the numerical results with respect to the threshold viscosity.

Table 1. Summary of the runs performed for numerical convergence of results with respect to threshold viscosity.

Run N	µthr (kPa s)	µ0 (kPa s)	Aver. Err. % (-)	Total Elapsed Time (s)
1	10	10.7	5.97	12,478.41
2	20	20.7	4.97	24,223.50
3	40	40.7	4.92	44,777.07
4	80	80.7	-	77,299.08

summarizes this set of runs showing the value of threshold viscosity and limiting viscosity adopted in each run. The threshold viscosity is progressively increased from 10 kPa s to 80 kPa s. The limiting viscosity is excluded by assuming in each run µ₀ > µ_thr. The last column of Table 1 shows the total elapsed time in every run.

The time evolution of the water free surface and of the interface between water and the bed load transport layer have been monitored for quantitative comparison among the different runs to obtain a measure of the influence of µ_thr on the mixture dynamics and to evaluate numerical convergence of results.

During the early phase of the water collapse the differences between the results of the simulations in Table 1 are almost negligible; this is reasonably because the mixture particles are interested by an elevated rate of deformation, far from the elastic–plastic regime. On the contrary, the results become different after the impact on the downstream wall because stagnation occurs; some solid grains at the lower left-hand corner of the tank enter the elastic–plastic regime owing to their low rate of deformation thus becoming fixed.

For the above reasons, the comparison of results from the runs is made at a reference time of t = 0.5 s. Figure 3 compares the free surface (fs) and the mixture–water interface (bls) close to the impact side for all the runs in Table 1. It can be seen that, as the threshold viscosity approaches the maximum value of 80 kPa s, the maximum rise of the bed-load sediment layer reduces and stabilizes at 0.3 m. The slope of the bottom sediment layer is less affected by changes in the threshold viscosity, as there are minimal differences between the four mixture–water interfaces at x > 0.05 m.

In Figure 3, it can be seen that the water free surface is poorly affected by changes in µ_thr, probably because the changes in the sediment layer profile are small, in percentage, with respect to the average water depth. The evolution of the free surface and of the bottom interface is quite regular in every run, as well as the pressure field that does not show significant instability phenomena at the transition interface between the water and the bed-load layer. The maximum pressure on the impacting area of the downstream vertical wall is around 10 kPa.

Based on the discussed results, it can be concluded that Run 2 allows obtaining an average profile and maximum run-up of the sediment layer which is rather close to the ones obtained in Runs 3 and 4 with higher values of the threshold viscosity; therefore, solution convergence can be obtained by assuming µ_thr^min = 20 kPa s and, in addition, this choice allows saving approximately 69% of computational time with respect to Run 4 and about 46% with respect to Run 3 (see Table 1).

To obtain a quantitative representation of the numerical optimization, the percentage average errors of the mixture–water interfaces are also reported in Table 1. These values have been obtained by calculating, near the left-hand corner of the tank, the average of the percentage differences between the mixture–water interface computed in the corresponding run and the mixture–water interface of Run 4. The average percentage error reduces about 1/6 from Run 1 to Run 2; indeed, it remains almost the same from Run 2 to Run 3. These results confirm the above conclusion.

A second set of runs was performed to investigate the effects of the limiting viscosity on saving the computational time and preserving numerical convergence of results.

Table 2. Summary of the runs performed for numerical optimization with respect to limiting viscosity (Runs 5 and 6). Run 7 tests an alternative rheological behavior with SPH mixture particles held fixed in case their viscosity values are higher than µ₁ (green line in Figure 1).

Run N	µthr (kPa s)	µ0 (kPa s)	µ1 (kPa s)	Aver. Err. % (-)	Total Elapsed Time (s)
2 (ref.)	20	20.7	-	-	24,223.50
5	20	15	-	5.4	19,163.10
6	20	5	-	5.7	7,089.25
7	-	-	5	-	7,098.34

shows the details of the second set of runs carried out. Run 2 from previous analysis has been assumed as a reference for comparison. Runs 5 and 6 adopt the same value of the threshold viscosity used in the reference run µ_thr = 20 kPa s while their limiting viscosity has been reduced to 15 kPa s and 5 kPa s, respectively, for evaluating the lower value µ₀^min of the limiting viscosity that reduces the computational time while preserving the results accuracy with respect to the reference run.

Run 7 was carried out to test an alternative rheological behavior (green line in Figure 1) with SPH mixture particles held fixed in case their viscosity values are higher than µ₁ = 5 kPa s.

As previously done, the obtained results have been analyzed at the reference time t = 0.5 s at which the maximum discrepancies were detected in first set of runs.

The obtained free surface (fs) and mixture–water interface (bls) profiles are compared in Figure 4. In this case, the water profile is also poorly affected by the changes in the limiting viscosity: the free surface for Runs 5–7 is practically coincident with the shape obtained in the reference Run 2.

When considering Runs 5 and 6, the interface between the bed-load sediment layer and the water is weakly sensitive to the decrease of the limiting viscosity. In these runs, the maximum height at which the suspended sediment can rise is approximately the same as the maximum height obtained in the reference Run 2. This means that numerical optimization can be achieved by assuming the limiting viscosity equal to 5 kPa s thus obtaining almost the same average percentage error as in Run 5 with a significant decrease of the computational time that for Run 6 is about 71% lower than the value obtained in the reference Run 2 (see Table 2).

The same benefit of Run 6 in term of saving computational time can be obtained by excluding the effect of the limiting viscosity and setting the viscosity µ₁ equal to 5 kPa s, as it can be seen from Run 7 in Table 2. The mixture dynamics deviates significantly from the reference run because the maximum run-up of the bed-load sediment layer is considerably underestimated. This can be seen in Figure 4, where the profile of the bed-load sediment layer remains below the reference profile in the interval x < 0.1 m.

To obtain clear representation of this behavior, Figure 5 shows a magnification of the mixture flow around the lower left-hand corner of the tank for both Run 6 (left-hand panel) and Run 7 (right-hand panel).

It can be seen that in Run 7 the run-up of the bed-load plume is smaller than Run 6: such behavior is due to the presence of a cluster of solid particles forced to remain fixed (colored in light brown in Figure 5b) whose apparent viscosity is higher than µ₁ = 5 kPa s. For this reason, these particles are excluded from the computation and kept fixed, thus neglecting part of the kinetic energy of the bed-load layer that prevents the bed-load plume from reaching the same height of Run 6 where the threshold viscosity is instead much greater to avoid the formation of a cluster of solid particles in the elastic–plastic regime (Figure 5a).

This result confirms that introduction of limiting viscosity allows obtaining a suitable physical representation of the phenomenon saving computational time. Further reduction of threshold viscosity to obtain similar benefit in term of computational time leads to a less consistent representation of the mixture flow.

3.2. Validation

This section describes the test case that has been performed for validating the mixture model implemented in SPHERA on a fast landslide. The test reproduces a two-dimensional experiment that was set up in 1968 at the Hydraulic Laboratory of Padua University for the analysis of the huge landslide occurred in 1963 on the left-hand slope of the Vajont artificial basin.

The 2D experimental setup and its functioning scheme are widely explained in [17] to which the reader is referenced for the details.

The left-hand panel in Figure 6 replicates the central part of the two-dimensional 1:500 scale model schematically reproducing the geometry of a representative cross section of the Vajont basin close to the dam. As the landslide dynamics was driven mainly by gravity, Froude similitude was adopted to convert both height and time to the full scale. Several experiments were carried out in the post-event campaign by testing: different run-out durations, ranging between 0.7 s and 22.5 s, corresponding to 15.7 s and 503.1 s, respectively, at full scale, different gravel sizes for the landslide, and two strokes of the rigid plate pushing the landslide (i.e., 0.5 and 0.8 m).

In the experimental campaign it was observed that the size of the granular material exerted scarce influence on the maximum run-up as, for a given stroke and run-out duration, the experimental points showed negligible dispersion. The larger the stroke, the higher the maximum run-up obtained for a given run-out duration. The maximum height reached for a given plate stroke increases as the run-out duration decreases.

Owing to the partial lack of information on the Padua experiment, only the toe of the landslide has been reproduced in the 2D numerical model. At the initial time, the still water level corresponds at full scale to a height of 700 m a.m.s.l.

The landslide is pushed into the reservoir by a rigid vertical wall that in the numerical model was assumed to be at the abscissa x = 0.0 m (not shown in the Figures) for simulating the effect of the rigid plate in the experimental facility. The movement of the wall is controlled by imposing its velocity components that follows a linear ramp to prevent incorrect numerical effects due to instantaneous change in the velocity at the boundary of the solid particles. During the simulations it was observed that few solid particles penetrate the wall and accumulate at its back. Even if the percentage of lost particles is so small that the effect on the landslide dynamics can be reasonably neglected, in future studies, an alternative approach available in SPHERA that may avoid particles loss will be tested: it consists of simulating rigid plate through a solid body with prescribed law of motion.

The water wave is generated by the landslide entering the basin, and the maximum run-up on the opposite slope is monitored and compared with the experimental result Figure 6, right).

The experimental test No. 48 is here simulated with run-out duration of 0.8 s (corresponding to 17.6 s at full scale) and observed maximum run-up of 0.712 m (corresponding to 863 m a.m.s.l.), which are close to the characteristics of the Vajont catastrophic event.

It must be pointed out that reliable prediction of the experimental run-up is obtained with the same run-out time of 0.8 s and reduced stroke of 0.36 m with respect to the test No. 48 where a stroke of 0.5 m was instead used. This may be partly related to the fact that only the toe of the landslide was here modeled owing to the lack of available information; therefore, the plate should accelerate an amount of granular material which is considerably smaller than in the Padua experiment.

Table 3. Summary of the relevant model parameters for sensitivity analysis (symbol - denotes parameter invariance with respect to reference Run V1).

Run N	ρs (kg/m3)	Ks (Pa)	αM (-)	ϕs (°)	Saturated (-)	µthr (kPa s)	µ0 (kPa s)	εf (-)	d50 (mm)
V1 (ref.)	2650.0	4.24 × 106	0.075	35.0	FALSE	320.0	5.0	0.35	4.0
V2	-	-	-	25.0	-	-	-	-	-
V3	-	-	-	-	-	-	10.0	-	-
V4	-	-	-	-	TRUE	-	-	-	-

summarizes the runs performed to evaluate model sensitivity to some relevant parameters. In all runs the same spatial resolution dx = 0.005 m was adopted and the total particles number was 12,309. The first run denoted as V1 is assumed as reference: the values assigned to geotechnical parameters of the landslide are physically consistent with the characteristics of the poorly graded sand with 3–4 mm diameter adopted in the experimental test number 48; in particular, the angle of internal friction is set equal to 35° following the guidelines on the non-cohesive soils [59]. In Runs V1–V3, the landslide material was assumed dry, thus neglecting the influence of pore water in the landslide portion below the still water level. This simplification was also assumed in a previous SPH model reproducing the test No. 48 of Padua experiment [17].

In Run V1 (reference run), the minimum value for the threshold viscosity has been selected by following the procedure described in previous section, resulting in µ_thr^min = 320 kPa s. The limiting viscosity has been set equal to the value of µ₀ = 5 kPa s as in the 2D erosive dam break test.

Run V2 aims at evaluating the effects on the landslide–water coupled dynamics of lowering the angle of internal friction to 25°.

Run V3 shows the influence of increasing the limiting viscosity to µ₀ = 10 kPa on the numerical results.

The effect of pore water has been accounted in Run V4, where it has been assumed that landslide portion below the stored water level was fully saturated.

Table 4. Maximum run-up obtained in the numerical simulations and comparison with the experimental result.

Run N	Total Elapsed Time (s)	Max Run-Up (m)	Time (s)	Zr−u (m a.m.s.l.)	Zexp (m a.m.s.l.)	Δη %
V1 (ref.)	21,395.15	0.718	1.35	870.3	863.0	5
V2	22,705.37	0.733	1.35	888.5		16
V3	43,814.04	0.708	1.35	858.2		3
V4	21,737.65	0.701	1.35	849.7		8

summarizes some relevant results from the performed runs: total elapsed time (for 1.6 s of simulated time on 16 cores @ 2.3 GHz 128 GB RAM), calculated maximum run-up and time at which it is reached. The maximum run-up obtained in each run N is translated into the corresponding maximum height $Z_{r-u}^{\left(N\right)}$ reached by the generated wave at full scale (with respect to mean sea level) and then it is compared with the maximum height Z_exp from the experimental test No. 48. The absolute percentage error in the last column is computed for each run N as:

$$\Delta {\eta }_{\%}^{\left(N\right)}=\frac{\left|Z_{r-u}^{\left(N\right)}-Z_{\exp }\right|}{Z_{\exp }-Z_{still}}100$$

(13)

where Z_still = 700 m a.m.s.l. is the stored water level.

As it can be seen in Table 4, the angle of internal friction ϕ does not affect significantly the computational time: in Run V2 the total elapsed time is about 6% greater than the reference Run V1. On the contrary, the maximum run-up is notably influenced by a reduction of ϕ from 35° to 25° leading to noticeable increment of the absolute percentage Δη % that becomes almost three times the one obtained for Run V1.

Such behavior may be related to the dependence of the apparent viscosity µ_fr from the sine of ϕ through the second of Equations (10): as the angle of internal friction decreases, also µ_fr reduces in a non-linear manner and lowers the intensity of the shear stresses in the granular material through the first of Equations (10). Therefore, a reduced part of kinetic energy is dissipated during the run-out, thus allowing the landslide to undergo higher deformation and pushing progressively its front more into the basin as ϕ reduces.

This is confirmed by the analysis of Figure 7 that compares the landslide and water profiles for both Run V1 and Run V2 at the instant of maximum run-up that is t = 1.35 s. As the angle of internal friction decreases the landslide front gets closer and touches the opposite slope (dash-dot brown line) providing an increased thrust on the stored water that reaches a height of maximum run-up (dash-dot blue line) significantly greater than Run V1 (continuous blue line).

From the comparison of both water profiles in Figure 7, it can be seen that the surface in Run V2 (dash-dot blue line) is, on average, above the surface obtained in Run V1 (continuous blue line). This means that the reduction of the angle of internal friction influences not only the shape of the landslide front, but it affects also the soil volume entering the basin: in particular, this volume increases because the displaced water volume is greater (it is worth noting that, at the time of maximum run-up, the velocity modulus is almost zero everywhere in the climbing water tongue, as shown in the following).

Dashed lines refer to experimental data: black line indicates the height of maximum run-up; blue line shows the water surface profile that is well reproduced by the Run V1 at abscissa between 0.7 m and 1.1 m, while above it is slightly overestimated probably because particle resolution cannot reproduce the actual thickness of the water blade.

Figure 8 shows the results from Run V3 (dash-dot lines) that are obtained for an increment of the limiting viscosity to 10 kPa s.

The comparison with the reference Run V1 (continuous lines) shows that this change in the limiting viscosity does not affect significantly the profile of the landslide. It can be also noticed that the water profile of Run V3 (dash-dot blue line) is almost overlapped to the corresponding profile obtained in Run V1 (continuous blue line). Therefore, the water surface profile from Run V3 also matches the experimental one (blue dashed curve) at abscissa between 0.7 m and 1.1 m. Even if the maximum wave run-up from Run V3 is a little bit closer to the experimental one, it seems scarcely affected by the considered change of limiting viscosity because the absolute percentage error Δη % remains of the same order of Run V1.

On the contrary, doubling the limiting viscosity induces a significant increase of the computational time of about 105%.

From these considerations follow that assuming µ₀ = 5 kPa s allows optimizing the simulation by halving the required computational time while preserving the accuracy of obtained results regarding landslide dynamics and wave run-up.

Figure 9 shows the comparison between the profiles from Run V1 (continuous lines) and Run V4 (dash-dot lines) at the same instant of maximum run-up (i.e., t = 1.35 s). Run V4 differs from Run V1 because in the former it has been assumed saturated condition for that portion of landslide below the still water level, as shown in Table 3.

The amount of soil penetrated into the reservoir is almost the same as in Run V1: this is confirmed by the fact that the water profiles are, on average, very close each other for both runs and hence the displaced water volume is about the same (the two extreme water volumes comprised between the water profiles of Run V1 and Run V4 balance each other). In the case of Run V4, the water surface profile is rather overlapped to the experimental one (blue dashed curve) below abscissa 1.1 m, including the left-hand front of the water surface that is closer to the experimental profile above the back of the landslide.

The two landslide fronts have different spatial distributions. In particular, the saturated layer below the still water level undergoes a higher deformation, as expected. Therefore, the landslide front in Run V4 (dash-dot brown line) appears more sharpened and touches the opposite side of the valley, while the maximum run-up reached by the generated wave is slightly lower than Run V1.

This behavior suggests that the way the landslide deforms and the resulting shape of its front influence the amount of kinetic energy transferred to the water and the maximum wave run-up: if the front is sharpened and sloping (as in Run V4) the maximum wave run-up decreases even if the protrusion of the landslide toward the opposite slope increases.

Conversely, if the landslide front is close to the vertical direction, as in the case of Run V2 in Figure 7, more kinetic energy can be transferred to the stored water resulting in an increased maximum run-up with respect to the reference run.

The results from Run V4 show the influence of pore water in the portion of the landslide below the still water level. As expected, the interstitial water contributes to lessen the effective stresses and to increase the soil susceptibility to deformation; therefore, the lower part of the landslide front is more prominent with respect to the case of Run V1, touching the opposite slope and showing a behavior that is similar to the model described in [17].

Similar behavior can be obtained by lowering the angle of internal friction (Run V2) but the front shape is less sharpened resulting in an overestimation of the maximum run-up. The maximum water height obtained in Run V4 is instead closer to the experimental result Z_exp; the absolute percentage error Δη % is half the one obtained in Run V2 and is rather close to the absolute percentage error of Run V1 (Table 4).

This suggests that in the present study the angle of internal friction has a physical meaning and does not represent a calibration parameter.

Figure 10 shows the field of velocity magnitude for the reference Run V1 at two representative instants: t = 0.80 s when the pushing plate comes to a stop (upper frame) and at t = 1.35 s when the maximum run-up is reached (lower frame).

In the upper panel, a wedge within the landslide which is characterized by zero velocity magnitude can be noticed: the wedge, contained between the sliding surface and the rigid plate, has a sloping side toward the basin forming an angle of about 45° with respect to the horizontal direction; the velocity of its solid particles is set equal to zero because their apparent viscosity exceeds the threshold viscosity µ_thr. The landslide body beyond this fixed wedge is characterized by a regular distribution of velocity magnitude ranging from zero, in close proximity of the wedge frontier, to the highest value of about 1 m/s at the landslide upper surface. The landslide front pushing the stored water is characterized by lower values of the velocity modulus with a maximum of 0.65 m/s close to the free surface. During this phase, the stored water is pushed by the landslide front entering the basin and is accelerated over the right-hand slope where the tip of the water tongue shows a maximum value of the velocity modulus around 0.9 m/s.

The lower panel of Figure 10 shows the instant t = 1.35 s of maximum run-up of the generated wave; this wave is characterized by almost zero velocity magnitude and, after climbing the opposite slope, it reaches the maximum height of 0.718 m in the laboratory frame of reference, corresponding to a height of 870.3 m at the real scale.

At the left-hand side of the free surface, a small wedge of water particles with a velocity magnitude of about 0.23 m/s can be noticed. During the descending phase of the generated wave toward the reservoir, these water particles further increase their velocity magnitude and cover the back of the landslide in accordance with the actual behavior of the Vajont landslide that was described above.

Previous results have been obtained for the reference spatial resolution dx = 0.005 m and evaluating the corresponding optimized value of the limiting viscosity µ₀ that allows obtaining the reference permitted percentage error Δη % of the estimate of the wave maximum run-up (see Table 4).

In the case where one is interested in optimizing the parameter µ₀ for coarser resolutions (i.e., greater values of dx), two additional simulations have been carried out to show how limiting viscosity changes by varying the resolution to obtain the allowed absolute percentage error. These simulations are summarized in

Table 5. Runs with reduced particles resolution: maximum run-up and percentage error with the experimental result.

Run N	dx (m)	µ0 (kPa s)	Total Elapsed Time (s)	Max Run-Up (m)	Time (s)	Zr−u (m a.m.s.l.)	Zexp (m a.m.s.l.)	Δη %
V5	6.25 × 10−3	7.5	14,476.08	0.702	1.35	850.9	863.0	7
V6	7.50 × 10−3	10.0	9334.12	0.699	1.35	847.3		9

where the adopted values for the spatial resolution dx and for the limiting viscosity µ₀ in each run are indicated. All the other parameters of Runs V5 and V6 are the same as Run V4 (see Table 3).

Figure 11 shows the profiles of the landslide (in brown) and of the water surface (in blue) for Run V5 (continuous lines) and for Run V6 (dotted lines with plus marker). These profiles are almost overlapped and the computed maximum run-up is very close to the one obtained for the higher resolution Run V4 (Table 4). In addition, the experimental water profile is suitably reproduced with major discrepancy above abscissa 1.1 m for the reason previously discussed.

The calculated values of the absolute percentage error Δη % are also reported in Table 5: for the investigated range of coarser resolutions the results show that, to obtain analogous values of Δη % that are of the same order of the permitted error in the estimate of maximum run-up, the parameter µ₀ should be varied almost linearly with respect to the particle resolution dx.

This result seems consistent with theoretical considerations because, when increasing dx, the limiting viscosity should be higher (i.e., should be less influential) to obtain the reference permitted error because part of the error is already due to the lower spatial resolution.

4. Conclusions

This paper shows an application of SPHERA v.8.0 (RSE SpA) FOSS code based on SPH method. SPHERA holds a reference mixture model [36], consistent with the packing limit of KTGF, for simulating two-phase free-surface fast flows involving water-sediment interaction.

The reference mixture model adopts a viscosity threshold µ_thr that has a physical meaning and its value is obtained through a convergence procedure. The viscosity threshold controls the abrupt transition of the sediment particles between the frictional regime and the elastic–plastic regime at very low shear rates where the viscosity is set to µ_thr: consequently, it avoids the unbounded growth of the apparent viscosity at very low deformation rates that would lead to significant reduction of the computational time step to assure stability.

The reference mixture model has been here improved by introducing a parameter µ₀ that, starting from the physical concept of the limiting viscosity that is exhibited by high polymer solutions, plays the role of a numerical parameter for code optimization acting in the frictional regime at low deformation rates near the transition zone to the elastic–plastic regime.

An important numerical consequence of introducing the limiting viscosity is the considerable drop of computational time that can be achieved through numerical optimization, as illustrated in a simplified 2D erodible dam break test. It has been shown that proper selection of the value to be assigned to the limiting viscosity leads to a considerable drop of the computational time of about 70% while assuring the reference permitted error of numerical results with respect to reference data.

Finally, SPHERA has been validated by simulating the 2D experimental test that was performed at Padua University in 1968 for investigating the relation between the landslide run-out velocity and the maximum run-up of the wave generated in the catastrophic event occurred in 1963 at the Vajont artificial basin. This may represent the first validation of such a model on the run-out of a complex and fast landslide interacting with water (both pore water and stored water). The results have shown that the mixture model allows obtaining an accurate representation of the generated wave whose profile and maximum run-up are very close to the experimental results. The geotechnical parameters of the landslide do not require any calibration and their values are consistent with the geo-mechanical properties of the considered granular material.

The proposed numerical model allows accounting for the effects of soil saturation on the prediction of the landslide dynamics and the relevant features of the generated wave. In particular, the experimental maximum run-up can be suitably reproduced by simulating the motion of the saturated layer that is more susceptible to higher deformation due to pore pressure that lower the effective stresses in the solid matrix.

For these reasons, the reference mixture model [36], empowered by the present developments on the limiting viscosity, seems suitable for general applications to the engineering investigation of landslide hazard at the slopes of artificial reservoirs and may help in developing probabilistic maps of risk variables.

Furthermore, by saving computational time, the proposed model may be suitable for the integration within a decision support system for natural hazard risk reduction to explore a number of potential risk-reduction options.