*Article* **The Hydromechanical Interplay in the Simplified Three-Dimensional Limit Equilibrium Analyses of Unsaturated Slope Stability**

**Panagiotis Sitarenios 1,\* and Francesca Casini <sup>2</sup>**


**Abstract:** This paper presents a three-dimensional slope stability limit equilibrium solution for translational planar failure modes. The proposed solution uses Bishop's average skeleton stress combined with the Mohr–Coulomb failure criterion to describe soil strength evolution under unsaturated conditions while its formulation ensures a natural and smooth transition from the unsaturated to the saturated regime and vice versa. The proposed analytical solution is evaluated by comparing its predictions with the results of the Ruedlingen slope failure experiment. The comparison suggests that, despite its relative simplicity, the analytical solution can capture the experimentally observed behaviour well and highlights the importance of considering lateral resistance together with a realistic interplay between mechanical parameters (cohesion) and hydraulic (pore water pressure) conditions.

**Keywords:** unsaturated slope; Ruedlingen field experiment; lateral resistance; limit equilibrium solution

#### **1. Introduction**

It is well known that partial saturation can play a key role in stabilizing both natural and artificial slopes. The interparticle forces due to the formation of water menisci in the macrostructure and/or intermolecular forces due to absorptive and osmotic phenomena increase the shear strength of geomaterials, improving stability conditions. Loss of this beneficial contribution can significantly hinder stability conditions and trigger slope failures. A characteristic example is rainfall induced landslides which comprise one of the most common geotechnical failures [1–4].

Geotechnical failures occurring due to changes in hydraulic conditions highlight the significance that unsaturated soil mechanics can have in engineering problems. Unfortunately, despite its significant importance in many geotechnical problems, unsaturated soil mechanics still struggles today to find a place in everyday geotechnical practice. The reasons for this are usually related to the relative complexity of the material laws involved (i.e., water retention behaviour), the increased complexity of the required laboratory tests (i.e., suction controlled shear strength determination) and the lack of simple calculation tools accessible to everyday practitioners with fundamental knowledge of soil mechanics.

In this respect, the present paper presents and evaluates a simple slope stability limit equilibrium solution for translational failure modes. The proposed solution is initially obtained under two dimensions (2D) and it is later extended to three dimensions (3D) by additionally examining the potential development of shear strength along the sides of a sliding block with a finite lateral dimension. Similar solutions for unsaturated soil conditions exist in the international literature [1,5–8], where some refer to much simpler geometries (i.e., [6]) and others to more complex geometries requiring the use of the method of slices (i.e., [8]). Those closer to the mechanism adopted in this study usually involve assumptions which prevent them from remaining valid under saturated conditions,

**Citation:** Sitarenios, P.; Casini, F. The Hydromechanical Interplay in the Simplified Three-Dimensional Limit Equilibrium Analyses of Unsaturated Slope Stability. *Geosciences* **2021**, *11*, 73. https:// doi.org/10.3390/geosciences11020073

Academic Editors: Jesus Martinez-Frias, Roberto Vassallo, Piernicola Lollino, Luca Comegna and Roberto Valentino Received: 28 December 2020 Accepted: 25 January 2021 Published: 8 February 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

especially when positive pore water pressures appear. The most common reasons are either the use of stress parameters which cannot recall Terzaghi's effective stress upon saturation (i.e., [7]) or assumptions with respect to the lateral in situ stress and its evolution with partial saturation, which are incompatible with classical soil mechanics principles (i.e., [1,5]).

The proposed stability solution is a natural extension of classical soil mechanics which is key in facilitating the adoption of such solutions in engineering practices [4]. To evaluate the proposed solution, this is applied to the analyses of a well instrumented field experiment where a slope failure was triggered by means of artificial rainfall—the Ruedlingen experiment [1,9]. This proves that, despite its simplicity, the proposed mechanism can very satisfactorily represent the experimentally observed behaviour both qualitatively and quantitatively. The presented discussion highlights the importance of considering the lateral development of the failure mechanism (3D) and of the natural transition to the saturated domain in capturing the behaviour.

#### **2. A Slope Stability Limit Equilibrium Solution for Translational Modes in Unsaturated Soils**

This section discusses the limit equilibrium stability calculation of a simple translational slope stability mechanism. Firstly, the assumed geometry is introduced and simple actions such as the weight of the sliding block are calculated. Then, the shear resistance along the failure plane is calculated considering the effect that partial saturation has on soil shear strength. Finally, the Factor of Safety (FoS) against failure is formulated by additionally considering the lateral resistance developing at the sides of a 3D failing soil block.

#### *2.1. The Assumed Failure Mechanism*

Figure 1 shows the assumed failure mechanism. A soil layer of thickness *z* rests on top of another material—the latter was assumed to be a priory of infinite strength. Both the soil surface and the interface between the two materials is described by two parallel planes with an inclination *θ* with the horizontal axis. The interface between the two materials acts as a potential failure surface (plane) and a simple translational failure mechanism develops where the upper layer tends to slide along the interface. 

 **Figure 1.** The analysed translational failure mechanism along a planar failure surface with slope *θ*.

 To examine such a typical translational failure mechanism, we can isolate a block of soil with a length of *l* = 1 m along the *x* axis which is the one parallel to the failure plane as shown in Figure 1. To obtain an initial solution for a two-dimensional (2D) failure mechanism, an analyses per running meter was performed. The soil was assumed to be homogeneous and isotropic, while suction (*s*) and degree of saturation (*Sr*) were both assumed to be constant along the entire soil layer with no variation of depth in favour of simplicity. Moreover, for the sake of generality of the produced solution, a vertical surcharge equal to *q* was assumed on the surface of the soil.

<sup>௫</sup> <sup>௬</sup>

<sup>௫</sup> <sup>௬</sup>

The unit weight of the soil varies with degree of saturation (*Sr*) following Equation (1) as: = ⋅ ,

=

<sup>௦</sup> <sup>௪</sup>

ௌ +⋅ 1+ ⋅ ௪ ,

 = ⋅ ⋅ , <sup>௫</sup> = ⋅ = ⋅ ⋅ ⋅ , <sup>௬</sup> = ⋅ = ⋅ ⋅ <sup>ଶ</sup> ,

$$\gamma = \frac{G\_{\rm S} + e \cdot S\_{\rm r}}{1 + e} \cdot \gamma\_{\rm w.} \tag{1}$$

<sup>௪</sup> = 9.81/<sup>ଷ</sup>

where *G<sup>s</sup>* is the specific gravity and *e* the void ratio of the soil, while *γ<sup>w</sup>* is the unit weigh of the water usually assumed to be equal to *γ<sup>w</sup>* = 9.81 kN/m<sup>3</sup> . Having defined the unit weight of the soil, the weight of the soil block (*W*) can be calculated and analysed as two components, *W<sup>x</sup>* and *Wy*, along the *x* and *y* axes, respectively, as shown in Figure 2a: ௫

$$\mathcal{W} = \gamma \cdot z \cdot \cos \theta,\tag{2}$$

$$\mathcal{W}\_{\mathbf{x}} = \mathcal{W} \cdot \sin \theta = \gamma \cdot z \cdot \cos \theta \cdot \sin \theta \,, \tag{3}$$

$$\dot{W}\_y = \dot{W} \cdot \cos\theta = \dot{\gamma} \cdot z \cdot \cos^2\theta,\tag{4}$$

**Figure 2.** (**a**) Analyses of the weight of the block and of the applied surcharge; (**b**) equilibrium of the block along *x* axis (parallel to the failure plane) and *y* axis (perpendicular to the failure plane).

To further consider the effect of the uniform surcharge, we calculated the corresponding force (*Q*) and its two components, *Q<sup>x</sup>* and *Qy*, along the *x* and *y* axes, respectively, as shown in Figure 2a:

$$Q = q \cdot \cos \theta \,\,\,\,\,\tag{5}$$

$$Q\_{\mathbf{x}} = Q \cdot \sin \theta = q \cdot \cos \theta \cdot \sin \theta\_{\prime} \tag{6}$$

$$Q\_y = Q \cdot \cos \theta = q \cdot \cos^2 \theta. \tag{7}$$

Finally, as portrayed in Figure 2b, the sum of the forces due to: (a) the weight of the soil block and (b) the surcharge is considered; initially on the *x* axis where *F<sup>x</sup>* plays the role of the destabilizing action trying to slide the considered block along the failure surface,

$$F\_{\mathbf{x}} = Q\_{\mathbf{x}} + \mathcal{W}\_{\mathbf{x}} = q \cdot \cos \theta \cdot \sin \theta + \gamma \cdot z \cdot \cos \theta \cdot \sin \theta = (q + \gamma z) \cdot \cos \theta \cdot \sin \theta,\tag{8}$$

and along the *y*-axis where *F<sup>y</sup>* is equal to the reaction force from the bottom layer:

$$N = F\_y = Q\_y + \mathcal{W}\_y = q \cdot \cos^2 \theta + \gamma \cdot z \cdot \cos^2 \theta = (q + \gamma z) \cos^2 \theta. \tag{9}$$

#### *2.2. Shear Resistance under Unsaturated Conditions*

The shear resistance at the interface is assumed to be entirely controlled by the shear strength parameters of the soil block and is described by the Mohr–Coulomb failure

envelope. The effect of partial saturation on shear strength is considered by using Bishop's (average) skeleton stress [10]. Bishop's average skeleton stress (*σ* ′ ) is calculated according to Equation (10) as the sum of the total stress (*σ*) plus suction (*s* = −*uw*, assumed air pressure: *u<sup>a</sup>* = 0); the latter is scaled for parameter *χ*.

$$
\sigma' = \sigma + \mathbf{s} \cdot \chi \tag{10}
$$

Note that although the prime symbol is used for Bishop's average skeleton stress, this should not be considered an effective stress variable; in fact, a second stress variable must be defined to adequately describe the behaviour under unsaturated conditions—e.g., suction [11,12]. If *χ* is properly selected to correspond to *χ* = 1 upon full saturation (*S<sup>r</sup>* = 1.0) and to *χ* = 0 for dry states (*S<sup>r</sup>* = 0.0), then Bishop's skeleton stress recalls Terzaghi's effective stress upon saturation which is a desirable feature for the unified description of saturated and unsaturated material states.

The Bishop's skeleton stress combined with common failure criteria for soils proves adequate when the interest focuses exclusively on the evolution of shear strength with partial saturation. This is supported by an ensemble of different studies [13–16], all of them suggesting that *χ* must be a function of *S<sup>r</sup>* which ensures that *χ* = 1 for *S<sup>r</sup>* = 1.0 and *χ* = 0 for *S<sup>r</sup>* = 0.0, as previously discussed.

To this end, for the examined sliding block failure mechanism the shear resistance along the interface is described by the Mohr–Coulomb failure criterion combined with Bishop's skeleton stress:

$$
\pi = \mathfrak{c} + \{\sigma + \mathfrak{s} \cdot \chi\} \cdot \tan \phi \tag{11}
$$

where *τ* is the shear strength of the material described by the well known shear strength parameter cohesion (*c*) and angle of internal friction (*φ*). The failure criterion was then applied at the failure surface of the examined stability mechanism to calculate the shear resistance force (*T*) (see also Figure 2b), which can be calculated as (also see Appendix A):

$$T = c + (q + \gamma z) \cos^2 \theta \cdot \tan \phi + s \cdot \chi \cdot \tan \phi \tag{12}$$

Finally, the safety factor against sliding of the block can be evaluated as the ratio of the forces resisting potential sliding and those trying to destabilize the examined block in the direction of the *x* axis. After some algebraic calculations for the 2D failure mechanism, the Factor of Safety (FoS) can be calculated according to:

$$FoS = \frac{T}{F\_{\chi}} = \frac{\frac{c}{\cos^2 \theta} + \frac{s \cdot \chi \cdot \tan \phi}{\cos^2 \theta} + (q + \gamma z) \cdot \tan \phi}{(q + \gamma z) \tan \theta} \tag{13}$$

#### *2.3. Extension to 3D-Lateral Shear Resistance*

The simple 2D translational slope failure mechanism presented hereinbefore and quantified through Equation (13) can be used to evaluate the stability of planar slope instabilities where failure extends enough to reasonably assume that stability can be evaluated by performing a "per running meter analyses" where a slice of the falling mass with a thickness of 1m at the centre of the failure area is isolated and studied. However, in many cases, slope failures have rather finite dimensions where the aforementioned assumption becomes problematic.

In such cases, the main shortcoming of the 2D mechanism is that it neglects the contribution of the shear resistance developing at the sides of the failing mass. Neglecting their contribution proves conservative and usually this is the justification accompanying the assumption of simple 2D slope stability methods in the design of slopes and cuts. However, this is not true when it comes to the back analyses of slope failures where neglecting the actual dimensions of the failure mass and the contribution of the lateral resistance can significantly hinder the outcome of the analyses (i.e., overestimation of soil shear strength parameters and/or suction contribution).

To this end, this section extends the discussed translational failure mechanism in the 3D domain by additionally considering the contribution of lateral resistance. For this, the geometry presented in Figure 3 is considered. The failing block has a width (*b*) which limits its lateral development, while the depth of the failure mass is (*z*) and the length is *l* = 1 m, which as a horizontal distance translates to *a* = cos *θ*, similar to the 2D mechanism. ሺ) = 1, = cos

**Figure 3.** The 3D failure mechanism in cut and planar view.

௦ ′ = ⋅ ′<sup>௩</sup> To calculate the shear resistance *T<sup>s</sup>* developing along each of the two lateral sides of the block, the Mohr–Coulomb failure criterion combined with Bishop's average skeleton stress is used again. The horizontal stress at the two lateral interfaces is assumed to be equal to *σ* ′ *<sup>h</sup>* = *K* · *σ* ′ *<sup>v</sup>* where *K* is the coefficient of lateral earth pressure. Applying the lateral earth pressure coefficient to Bishops' stress as opposed to total stress that other solutions assume (i.e., [1]) makes the derived solution applicable under positive pore water pressure regimes as well, while recent experimental evidence on soil anisotropy and partial saturation [17] also support such an assumption. The vertical Bishop's skeleton stress is calculated in the middle of the depth (*z*) as:

$$
\sigma'\_v = \sigma\_v + \mathbf{s} \cdot \boldsymbol{\chi} = q + \boldsymbol{\gamma} \cdot \frac{\boldsymbol{z}}{2} + \mathbf{s} \cdot \boldsymbol{\chi} \tag{14}
$$

and the corresponding horizontal stress as:

<sup>௦</sup> = <sup>௦</sup>

$$
\sigma\_h' = \mathbf{K} \cdot \sigma\_v' = \mathbf{K} \cdot \boldsymbol{q} + \mathbf{K} \cdot \boldsymbol{\gamma} \cdot \frac{\boldsymbol{z}}{2} + \mathbf{K} \cdot \mathbf{s} \cdot \boldsymbol{\chi} \tag{15}
$$

2

The stress as calculated in the middle of the layer is assumed uniform along the lateral sides of the excavation and the shear strength of the soil is then calculated as:

$$
\pi\_{\rm s} = \mathcal{c}\_{\rm s} + \sigma\_{\rm h}^{\prime} \cdot \tan \phi\_{\rm s} \tag{16}
$$

<sup>௦</sup> <sup>௦</sup> where *c<sup>s</sup>* and *φ<sup>s</sup>* are the cohesion and frictional resistance, respectively, along the lateral sides of the failure block. They are assumed to be different from their soil counterparts in favour of generality of the solution, by allowing different shear strength parameters to be assigned to the lateral resistance. Having defined the shear strength along the two lateral sides, we can calculate the corresponding shear resistance, where after some algebraic calculations (see Appendix B) the following expression is derived:

$$T\_s = c\_s \cdot z \cdot \cos\theta + \mathbf{K} \cdot \boldsymbol{\eta} \cdot z \cdot \cos\theta \cdot \tan\phi\_s + \mathbf{K} \cdot \boldsymbol{\gamma} \cdot \frac{z^2}{2} \cdot \cos\theta \cdot \tan\phi\_s + \mathbf{K} \cdot \mathbf{s} \cdot \boldsymbol{\chi} \cdot z \cdot \cos\theta \cdot \tan\phi\_s \tag{17}$$

ଶ

Finally, using the additional resistance along the sides of the 3D mechanism, the safety factor against sliding is reformulated:

$$FoS = \frac{T + \frac{2T\_s}{b}}{F\_x} \tag{18}$$

Note that the additional resistance at the two sides (2*Ts*) is normalized over the width (*b*) of the block to combine it mathematically with the already determined resistance *T* and action *F<sup>x</sup>* for the 2D mechanism. Finally, by using Equations (8), (12) and (17) to substitute for *Fx*, *T* and *T<sup>s</sup>* , respectively, in Equation (18), we can calculate that the Factor of Safety (FoS) against sliding for the 3D mechanism as:

$$\mathrm{FoS} = \frac{\frac{c}{\cos^2 \theta} \left(1 + \frac{2\rho\_c \cos \theta}{b}\right) + \frac{s\chi t \tan \theta}{\cos^2 \theta} \left(1 + \frac{2\rho\_b K z \cos \theta}{b}\right) + q \cdot \tan \theta \left(1 + \frac{2\rho\_b K z}{b \cdot \cos \theta}\right) + \gamma z \cdot \tan \theta \left(1 + \frac{\rho\_b K z}{b \cdot \cos \theta}\right)}{(q + \gamma z) \cdot \tan \theta} \tag{19}$$

In Equation (19), and in order to obtain a more elegant solution, the shear strength parameters used for the resistance along the inclined failure surface and along the sides of the block are correlated using the following factors:

$$
\rho\_s = \frac{c\_s}{c} \tag{20}
$$

$$
\rho\_{\varphi} = \frac{\tan \phi\_s}{\tan \phi} \tag{21}
$$

#### **3. Evaluation Using an Actual Case**

To evaluate the planar slope stability mechanism previously presented and to demonstrate the potential importance of lateral resistance in analysing slope instabilities, the results from the Ruedlingen full scale experiment [1,18] will be used and compared with the predictions of the proposed stability mechanism.

#### *3.1. The Ruedlingen Full Scale Experiment*

The Ruedlingen full scale field test [1,18] was performed to study the response of a steep forested slope, subjected to artificial intense rainfall within the context of the multidisciplinary research programme on the "Triggering of Rapid Mass Movements in steep terrain" (TRAMM). The experiment was carried out in northern Switzerland in a forested area near Ruedlingen village. The selected experimental site is located on the east-facing bank of the river Rhine, with an average slope angle of approximately 38◦ . An orthogonal area, with a length of 35 m and a width of 7.5 m, was instrumented with a wide range of devices including earth pressure cells, piezometers, tensiometers, time domain reflectometers (TDRs), and acoustic and temperature sensors, organized in clusters, as seen in Figure 4b. A landslide triggering experiment was conducted in March 2009 where the slope was subjected to artificial rainfall with an average intensity of 20 mm/h at the upper one-third (approximately) of the slope and with 7 mm/h in the remaining lower part leading to a shallow slope failure after 15 h.

The failure concentrated at the upper part of the instrumented area (see Figure 4b) and extended within a surficial layer of a medium to low plasticity (average PI ∼10%) silty sand (ML) [19]. The developed failure mechanism can be characterized as a surficial planar failure. The mobilized area had a length of ∼17 m, a width of ∼7 m while the failure surface can be practically assumed to be parallel to both the soil surface and the underlying soil–bedrock interface with its maximum depth exciding 1 m [18]. This soil layer rests on top of a fissured Molassic bedrock (sandstone and marlstone—see Figure 4a) [20,21]. Sitarenios et al. [9] numerically analysed the triggered slope instability and, by combining field measurements with numerical results, they concluded that failure happens under saturated conditions and is triggered by water exfiltration from the bedrock which results in considerable pore water pressure build-up at the soil–bedrock interface.

**Figure 4.** (**a**) Geomorphology of the greater experimental area (after [20]); (**b**) a schematic of the bedrock topography and of the instrumented area of the experiment (after [22]).

#### *3.2. 2D Solution—Neglecting Lateral Resistance*

In this section, the proposed limit equilibrium mechanism is used to evaluate the stability of the Ruedlingen slope by neglecting any contribution of lateral resistance or, in other words, by assuming an infinite failure mechanism perpendicular to the examined plane. When utilizing the mechanism, it is important to first establish the hydraulic and mechanical parameters of the Ruedlingen soil; to enable comparison of the results with the study of Sitarenios et al. [9], the same experimental data will be examined and similar assumptions will be put forward.

 <sup>௪</sup> ,௫ , = /ሺ1 − ) <sup>௪</sup> Starting with the hydraulic behaviour and according to what was previously discussed, calculation of the soil shear strength under unsaturated conditions necessitates knowledge of the degree of saturation as the latter is the most frequently made assumption for parameter (*χ*). To calculate the degree of saturation and thus parameter (*χ*), the Water Retention Model (WRM) of Equation (22) is used where *s* is the suction level, *b<sup>w</sup>* is a model parameter controlling the shape of the reproduced Water Retention Curve (WRC) and *Sr*,*max* and *Sr*,*min* are the maximum and residual degrees of saturation, respectively. It is the a Van Genuchten type WRM [23] further including dependence of the behaviour of the void ratio (*e* = *n*/(1 − *n*)) through parameter *P*, as described in Equation (23), where *P*<sup>0</sup> and *n*<sup>0</sup> are reference values and parameter *a<sup>w</sup>* controls the rate of change of the WRC with the void ratio.

$$\chi = \mathcal{S}\_r = \mathcal{S}\_{r, \text{res}} + \left(\mathcal{S}\_{r, \text{max}} - \mathcal{S}\_{r, \text{res}}\right) \left[1 + \left(\frac{s}{P}\right)^{\frac{1}{1 - b\_W}}\right]^{-b\_W} \tag{22}$$

$$P = P\_0 \cdot \exp[\alpha\_w (n\_0 - n)]\tag{23}$$

= ⋅ expሾ௪ሺ − )ሿ = 0.9 ⋅ It should be highlighted that the water retention behaviour is key to the accurate representation of the unsaturated shear strength. Consequently, in case of significant hysteresis in the water retention behaviour, attention must be given to the utilization of this branch of the curve (drying vs. wetting) which best represents the physically analysed problem. This study addresses a wetting problem (rainfall water infiltration) and, in this respect, the water retention behaviour is calibrated based on a set of suction controlled wetting laboratory tests [24]. The calibrated behaviour can be seen in Figure 5a and the parameters used are shown in Table 1. Figure 5b shows the WRC curve that is used in the analyses, which corresponds to an initial void ratio of *e* = 0.9 as this has been measured in the field and also used in previous studies [9]. It is accompanied by the evolution of the *s* · *S<sup>r</sup>* term which in fact reflects the increase in the shear strength of the material with suction.

⋅ **Figure 5.** (**a**) The calibrated Water Retention Model (WRM) (laboratory data from [24]); (**b**) the Water Retention Curve (WRC) used in the analyses with the corresponding evolution of the *s* · *S<sup>r</sup>* term.



=0 = 32 <sup>௦</sup> = 2.65 = 38 =0 With the WRM calibrated, Equation (11) is used to calculate the safety factor of the Ruedlingen slope under different combinations of suction and slope heights. With respect to the shear strength parameters involved in Equation (11), cohesion is set to zero (*c* = 0) and the friction angle equal to *φ* = 32 [9,19], while the unit weight (*γ*) of the slope follows the evolution of saturation through Equation (1), with specific gravity equal to *G<sup>s</sup>* = 2.65. Finally, the inclination of the failure surface is set equal to the average inclination of the slope (*θ* = 38) [1] and no surcharge is considered on the soil surface (*q* = 0).

Figure 6 presents the results of a parametric study. Figure 6a portrays the evolution of the FoS with suction for various slope heights ranging from 0.5 to 2.75 m. The reason why the height of the slope is handled parametrically is related to the in situ variability of the thickness of the soil layer which was found, in general, to range from 0.5 to 4.5 m [18,20], while in the part of the slope where failure occurred the variation was smaller, with depths up to 2.75 m. In Figure 6b, similar results are presented with depth in the horizontal axis and the various curves represent different suction levels from zero (0) up to 10 kPa, following the in situ initial suction values as those have been reported in [18].

⋅

**Figure 6.** (**a**) Evolution of Factor of Safety (FoS) with suction for different slope heights; (**b**) evolution of FoS with slope height for different suction levels; results from a 2D mechanism (no lateral resistance).

The results of Figure 6 clearly demonstrate the dominant effect that partial saturation can have on stabilizing a planar slope. It is characteristic that when no suction is present (*s* = 0), a cohesionless slope (*c* = 0) has a safety factor which follows the *tan*(*φ*) *tan*(*θ*) ratio; the latter, for the examined slope geometryand soil shear strength is equal to *tan*(32) *tan*(38) = 0.8. In both Figure 6a,b, it can be observed that even a small suction increase significantly alters stability conditions with various combinations of suction and slope heights exhibiting stable behaviours (FoS ≥ 1).

To better understand this interplay between the height of the slope and suction, Figure 7 presents stability charts correlating the two quantities, with suction represented by the horizontal axis and slope height by the vertical axis, respectively. Figure 7a suggests that there is an almost linear relation between the increase in depth and the required increase in suction for the slope to remain stable (FoS > 1), while in Figure 7b the same relation is presented using a logarithmic scale with suction. Figure 7b helps to understand that the increase in the stable slope height with suction in fact follows the trend of the increase in the shear strength with suction as the latter is controlled by the increase in the *s* · *S<sup>r</sup>* term in Figure 5b. ⋅

**Figure 7.** Stability chart; (**a**) stable vs. unstable suction and slope height combinations; (**b**) same as (**a**) in a semilogarithmic plot; results from a 2D mechanism (no lateral resistance).

Examining the stability chart of Figure 6 in relation to the observed behaviour during the Ruedlingen field experiment, it seems that the suggested evolution of stability with suction and slope height does not match the observed behaviour. If we consider an average slope height of 1.5 to 2 m, Figure 7 suggests that a minimum suction value of 10 to 15 kPa is required for stability. This may partially explain why the slope remains stable under the initial hydraulic regime which suggests a value of around 10 kPa, but at the same time it is obvious that even the smallest decrease in suction should cause failure. However, both in situ measurements [1] and numerical analyses [9] clearly indicate that even complete loss of suction is not enough to trigger the Ruedlingen failure. In this respect, it seems reasonable to assume that additional factors contribute to stability with respect to those examined with the present 2D failure mechanism.

#### *3.3. 3D Solution—Including Lateral Resistance*

Lateral resistance can be one of the factors contributing to improved stability conditions with respect to those of Figure 7. To consider the effect of lateral resistance in the examined problem, the analysis of the previous section is repeated by additionally considering the lateral resistance at the two sides of a failing block with a width equal to b = 7.5 m. For the side resistance, the same friction angle as for the soil–bedrock interface is assumed (*ρφ*=1.0), while for the lateral earth pressure coefficient an arbitrary value of K = 0.5 is chosen, close to the value suggested by Jaky's formula for *φ* = 32 ⇒ *K* = 1 − sin 32 = 0.47.

Figure 8 presents results similar to those of Figure 6; however, Figure 8 includes the additional contribution of lateral resistance. The results clearly demonstrate the beneficial effect of lateral resistance in increasing FoS under given combinations of slope height and suction. It is interesting to observe that Figure 8a seems to suggest that the effect of suction in improving stability conditions becomes less profound as the slope height increases. Moreover, Figure 8b indicates that under relatively small slope height values (i.e., less than 0.6 m), the influence of suction is quite significant, while as the slope height increases its effect becomes less profound, reaching a minimum for a given depth and then slightly increasing again.

= 32 ⇒ = 1 − sin 32 = 0.47

థ

**Figure 8.** (**a**) Evolution of FoS with suction for different slope heights; (**b**) evolution of FoS with slope height for different suction levels; results from a 3D mechanism (with lateral resistance).

 This effect is more profound in the results shown in Figure 9, presenting a stability chart similar to Figure 7 with the addition of lateral resistance. It is observed that as the height of the slope increases, positive pore water pressure values are needed to cause instability. Qualitatively, this trend in behaviour is much better aligned with the behaviour of the Ruedlingen slope, where failure happened following considerable water pressure build-up which reduced the available shear strength at the soil–bedrock interface. Moreover, both postfailure in situ observations and numerical results suggest that the depth of the failure surface ranged from 1.0 to 2.0 m. This is in good agreement with Figure 9, where within the same range of *z* values, the stability chart exhibits a peak coinciding with the most unfavourable stability conditions for the given assumptions.

**Figure 9.** Stability chart; (**a**) stable vs. unstable suction and slope height combinations; (**b**) as (**a**) in a semilogarithmic plot; results from a 3D mechanism (with lateral resistance).

5

=1

= 1.0

=5

= 0, 1, 2, 3, 4

However, from a quantitative perspective and within the same range of slope heights, the required suction for stability suggests that the slope should have failed before reaching full saturation contrary to what was observed in situ. This discrepancy can be attributed to several factors which can result in better stability conditions than those described by the analysis used for this study. Amongst others, three of the most likely potential candidates for the different behaviours can be (a) the presence of some cohesion either due to fines or due to the contribution of vegetation and roots [25]; (b) a less extended mechanism in the third direction (smaller b) and (c) a different WRC. However, it should be recognized that the WRC can only affect the behaviour in the unsaturated regime and still cannot explain why a higher positive pore water pressure may be needed to trigger the instability.

#### *3.4. 3D Solution—The Effect of Cohesion and 3D Mechanism Development*

The previous section demonstrated how the inclusion of the lateral resistance in the analyses can justify an improvement of the stability conditions for the analysed Ruedlingen field experiment. However, it was not feasible to quantitatively approach the behaviour satisfactorily. As mentioned before, either the presence of cohesion or a narrower than assumed mechanism can serve as the most probable explanations and this hypothesis is put to test in this section.

To examine the effect of different assumptions regarding the aforementioned factors, a parametric analysis is conducted to check the sensitivity of the stability charts on the presence of cohesion (c) and on the width (b) of the 3D mechanism. Figure 10a presents how the stability curve differentiating between stable and unstable conditions changes with an increase in cohesion. Five different cohesion values, namely *c* = 0, 1, 2, 3, 4 and 5 kPa (from the blue to the red curve), are examined and the same cohesion values also control the lateral resistance by assuming that *ρ<sup>c</sup>* = 1.0 (Equation (20)). The plotted charts indicate a significant contribution of cohesion to the stability conditions. Even the smallest nonzero value of *c* = 1 kPa is enough to justify a stable slope under the presence of the smallest suction and irrespective of its height. This is in accordance with the experimentally observed behaviour, while the presence of roots can justify even higher values of cohesion [25] which, if adopted (i.e., *c* = 5 kPa), can explain why significant pore pressures up to almost 10 kPa were needed in the soil–bedrock interface for failure to occur in the Ruedlingen experiment.

**Figure 10.** (**a**) Effect of cohesion on the stability chart; (**b**) effect of the 3D mechanism width on the stability chart; results from a 3D mechanism (with lateral resistance).

 = 7.5 Figure 10b offers a similar discussion by explaining the effect of the width of the assumed mechanism. The main assumption put forward in the analyses presented hitherto is that the width of the mechanism is *b* = 7.5 m, coinciding with the width of the experimental area (dashed curve on Figure 9b). However, the in situ postfailure observations suggest that the failing mass had a smaller average width (see also Figure 4b), and the

=7

=

=

=

=

=

=1

results of Figure 10b indicate that such a reduction in the lateral development of the failure mechanism is quite beneficial for stability. In more detail, Figure 10b presents the evolution of the stability curves with the width of the mechanism, namely from *b* = 7 m down to *b* = 1 m (from the blue to the red curve). It is interesting that the effect of lateral development was qualitatively different compared to that of cohesion, with the slope height playing a key role. It is observed that as the slope becomes higher and the mechanism narrower, the effect of the lateral resistance on improving stability becomes higher. This is a clear reflection of the fact that the ratio of the area of the lateral side over the area of the soil–bedrock interface becomes higher, increasing the contribution of the side resistance and limiting the contribution of the failure plane resistance, thus progressively making the lateral resistance more and more dominant.

Regarding the in situ behaviour, it is clear that a narrower mechanism cannot explain the in situ behaviour alone; as for the slope height relevant to the Ruedlingen soil (*z* = 1 to 2 m) and the range of the lateral development of the mechanism (*b* = 5 to 7 m), the slope is still unstable, even though there is a noticeable improvement in stability conditions with respect to the refence scenario (*b* = 7.5 m). In this respect, the best candidate to justify the experimental behaviour is cohesion. From the results of Figure 10, it seems that reasonable assumptions regarding a small presence of cohesion in the order of *c* = 3 to 4 kPa, combined with a relatively smaller average mechanism *b* = 5 to 7 m, can explain very satisfactory the experimentally observed behaviour during the Ruedlingen experiment, both qualitatively and quantitively.

To test this assumption, four combinations of cohesion and slope width within the aforementioned range of values are finally analysed Figure 11b. In the same graph, the grey shaded area indicates the range of values of slope height and "suction" (negative suction values correspond to positive pore water pressure) most representative of the in situ failure conditions. In more detail, Figure 11a presents the failure mechanism (after [1]) where the average depth of the failure surface was more than 1 m. It also contains a plot of the evolution of the in situ and of the numerically determined pore water pressure at the failure area (after [9]) indicating that an average pore water pressure in the order of 5 kPa was present in the slope mass at failure. The very good match of the stability curves with the shaded area is a clear indication that the relatively simple 3D translational failure mechanism used in this study can provide a very good insight into the Ruedlingen failure.

**Figure 11.** (**a**) The Ruedlingen slope failure (after [1]) and the evolution of pore water pressure in the slope mass (after [9]); (**b**) stability curves for different combinations of cohesion and slope widths which can accurately capture the observed instability (grey shaded area).

l = 1 m

=⋅+ ሼ ⋅+⋅⋅ሽ ⋅ l = 1m

= + ⋅ + ⋅ ⋅

#### **4. Conclusions**

In this paper, a simplified limit equilibrium slope stability solution for infinite translational failure mechanisms was introduced. The solution uses Bishop's average skeleton stress to account for the increase in the shear strength under unsaturated conditions. It was extended in 3D dimensions by considering a potential limited lateral development of the failure mechanism by accounting for the lateral shear resistance that develops across the sides of the failing mass. Contrary to other similar attempts, the proposed solution assumes that the lateral earth pressure coefficient (K) describes the ratio of the horizontal over the vertical Bishop's stress and not only the net stress ratio. By doing so, the solution remains valid under fully saturated conditions where Bishop's average skeleton stress recalls Terzaghi's effective stress.

The mechanism was evaluated using a set of actual data from a full scale field experiment—the Ruedlingen artificial rainfall induced landslide. The proposed stability solution, despite its simplicity, can capture very satisfactory the experimentally observed behaviour. Additionally, its simplicity and very limited calculation costs enable an easy evaluation of different scenarios with respect to slope geometries and soil parameters which can prove very useful in back analysis of slope failures. Such simple validated tools for the analyses of engineering problems under unsaturated conditions, formulated as a natural extension of similar classical soil mechanics solutions, are very important in facilitating a wider adoption of unsaturated mechanics in everyday geotechnical practice.

**Author Contributions:** Conceptualization, F.C.; Investigation, P.S.; Methodology, P.S.; Writing original draft, P.S.; Writing—review & editing, F.C. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A**

Starting from the Mohr–Coulomb criterion, Equation (11) is multiplied with the area of the failure surface, which is equal to *l* = 1 m, by 1 m for the distance perpendicular to the examined 2D problem (analyses per running meter. Thus, the shear resistance as a force is derived

$$T = c \cdot l + \{\sigma\_n \cdot l + s \cdot \chi \cdot l\} \cdot \tan \phi \text{ (for } l = 1 \text{ m) or }\tag{A1}$$

$$T = \mathcal{c} + \mathcal{N} \cdot \tan \phi + \mathbf{s} \cdot \boldsymbol{\chi} \cdot \tan \phi \tag{A2}$$

In Equation (A2), *N* represents the reaction force at the failure interface as this has been calculated with Equation (9) and so we can finally write:

$$T = c + (q + \gamma z) \cos^2 \theta \cdot \tan \phi + s \cdot \chi \cdot \tan \phi \tag{A3}$$

#### **Appendix B**

In Equation (16), which represents the shear strength along the side of the examined mechanism, Bishop's average skeleton stress is described through Equation (15) to obtain:

$$\tau\_{\rm s} = c\_{\rm s} + \sigma\_{\rm h}^{\prime} \cdot \tan \phi\_{\rm s} = c\_{\rm s} + \left(\mathbf{K} \cdot \boldsymbol{\eta} + \mathbf{K} \cdot \boldsymbol{\gamma} \cdot \frac{\boldsymbol{z}}{2} + \mathbf{K} \cdot \mathbf{s} \cdot \boldsymbol{\chi}\right) \tan \phi\_{\rm s} \tag{A4}$$

Focusing on one side of the 3D mechanism, the shear strength above develops along the lateral side of the mechanism which is portrayed in Figure 1 by the light grey shaded

area. It is a parallelogram with an area equal to *z* · *a* = *z* · *cosθ*. Multiplying Equation (A4) with *z* · *cosθ*, the shear resistance along the side is:

$$T\_s = c\_s \cdot z \cdot \cos\theta + \left(\mathbf{K} \cdot \boldsymbol{\eta} + \mathbf{K} \cdot \boldsymbol{\gamma} \cdot \frac{z}{2} + \mathbf{K} \cdot \mathbf{s} \cdot \boldsymbol{\chi}\right) z \cdot \cos\theta \cdot \tan\phi\_\mathbf{s} \text{ or }\tag{A5}$$

$$T\_s = c\_s \cdot z \cdot \cos\theta + K \cdot q \cdot z \cdot \cos\theta \cdot \tan\phi\_s + K \cdot \gamma \cdot \frac{z^2}{2} \cdot \cos\theta \cdot \tan\phi\_s + K \cdot s \cdot \chi \cdot z \cdot \cos\theta \cdot \tan\phi\_s \tag{A6}$$

#### **References**


## *Article* **Experimental and Numerical Investigations of a River Embankment Model under Transient Seepage Conditions**

**Roberta Ventini 1,\*, Elena Dodaro <sup>2</sup> , Carmine Gerardo Gragnano <sup>2</sup> , Daniela Giretti <sup>3</sup> and Marianna Pirone <sup>1</sup>**


**Abstract:** The evaluation of riverbank stability often represents an underrated problem in engineering practice, but is also a topical geotechnical research issue. In fact, it is certainly true that soil water content and pore water pressure distributions in the riverbank materials vary with time, due to the changeable effects of hydrometric and climatic boundary conditions, strongly influencing the bank stability conditions. Nonetheless, the assessment of hydraulic and mechanical behavior of embankments are currently performed under the simplified hypothesis of steady-state seepage, generally neglecting the unsaturated soil related issues. In this paper, a comprehensive procedure for properly defining the key aspects of the problem is presented and, in particular, the soil characterization in partially saturated conditions of a suitably compacted mixture of sand and finer material, typical of flood embankments of the main river Po tributaries (Italy), is reported. The laboratory results have then been considered for modelling the embankment performance under transient seepage and following a set of possible hydrometric peaks. The outcome of the present contribution may provide meaningful geotechnical insights, for practitioners and researchers, in the flood risk assessment of river embankments.

**Keywords:** riverbank; unsaturated soils; water retention curve; unsaturated permeability curve; transient seepage; slope stability

#### **1. Introduction**

Riverbanks are passive defense works, consisting of earthen structures built along the edges of a stream or river channel to prevent flooding of the adjacent land. In recent years, the risk assessment of river embankments has received great attention as consequence of the considerable socio-economic damages caused by flooding. Since 2000, floods in Europe have caused at least 700 deaths, the evacuation of about half a million people, and at least 25 billion euro in insured economic losses [1]. The economic impact of climate-related events varies considerably across countries. In absolute terms, the highest economic losses in the period 1980–2019 were registered in Germany, followed by Italy, then France [2].

Nevertheless, investments in earth retaining structures are relatively limited compared to other hydraulic engineering works, such as water regulation and canalization. In many countries, including Italy, minor embankments are not even subjected to any formal legislation or published technical guidelines and this has a severe impact on the performance of levees and safety of territories [3].

The design of riverbanks and their individual components depends on many factors, including functionality under different hydraulic loadings, feasibility, but also environmental constraint and availability of the materials. Numerous cross-sectional variations exist, each with the primary objective of reducing flood risk within the surrounding area. The simplest type of section is the homogeneous earthfill, composed of granular and cohesive

**Citation:** Ventini, R.; Dodaro, E.; Gragnano, C.G.; Giretti, D.; Pirone, M. Experimental and Numerical Investigations of a River Embankment Model under Transient Seepage Conditions. *Geosciences* **2021**, *11*, 192. https://doi.org/10.3390/ geosciences11050192

Academic Editors: Jesus Martinez-Frias, Roberto Vassallo, Luca Comegna and Roberto Valentino

Received: 31 March 2021 Accepted: 25 April 2021 Published: 29 April 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

soil, obtained from site excavation or borrow sources, though riverbank may also present an impermeable core (zoned embankment) or be composite with superstructures, waterside or landside berm and internal structures [4].

Riverbank failure, defined as the inability to prevent the inundation of the area surrounding a watercourse, may occur due to various hydraulic and structural causes, such as overtopping, seepage flow through foundation layer, internal and external erosion, slopes instability, scour and liquefaction [5]. Some possible failure mechanisms are reported in Figure 1. One of the most critical is represented by the overall macro-instability of the slopes triggered by the unfavorable seepage of water through the riverbanks, due to rivers persistently running at unusually high levels. An increasing number of studies have recognized changes in pore water pressures within a riverbank as a fundamental factor in determining conditions of instability and in the triggering of bank failures [6–9]. Moreover, burrowing animals are acknowledged by agencies responsible for earthen dams and riverbanks to have a significant role in substantial and costly damages [10].

**Figure 1.** Overview of main riverbank failure mechanisms in earthen structures, modified from Kok et al. [11]: (**a**) Overtopping; (**b**) Uplift and piping; (**c**) Erosion of waterside slope; (**d**) Macro-instability in landside slope; (**e**) Macro-instability in waterside slope.

The development of high pore water pressures, due to a persistent groundwater flow, actually leads to a decrease in the soil shear resistance within the riverbank that may result in catastrophic slope failure on the landside [12,13]. This type of collapse strongly depends on the retention characteristics of the filling material such as suction variations with time, also due to seasonal variations of the meteorological conditions. In fact, the recurring soil water content variations in the embankment induced by hydrometric fluctuations and rainfalls occurrences may dramatically reduce the suction contribution to shear strength and eventually modify the microstructure of the soil, even creating a cracking state in the riverbank that could form preferential flow channels with general increase of the soil permeability and decrease of the shear strength [14]. As an example, the study of the breach that occurred along the river Secchia, Northern Italy, on 19 January 2014, provides clear evidence of a landside slope failure, caused by a series of simultaneous river stage raising combined to rainfall on the embankment surface and the wildlife activity along fluvial systems [15]. For the failure of the landslide, it is well known that toe drainage by gravelly materials is able to effectively lower the phreatic surface and to prevent water from seeping out of the 'landslide slope'.

The failure of the riverside is relatively different and could be triggered by external erosion of the toe and foundation of the embankment, primarily due to the action of waves, currents or turbulence; filling material is so removed from the waterside surface, leading to loss of stability and consequently to deep or shallow slope failure. A recurrent phenomenon causing riverside slope failure is also the rapid drawdown, typically occurring when the hydrometric level undergoes a rapid decrease, after a high stage, as the bank material

is still near a saturated condition and the confining pressure of the river decreases to zero [8,16,17]. A frequently adopted solution for averting the aforesaid failure generally consists in protecting the riverside slope of the embankment using armourstone, gabions and concrete slabs. In general, different types of interventions, including the realization of components, such as impervious layers, cut-offs barriers, drains, and reinforcement layers within the embankment sections, are indeed considered by management and streambank control agencies to improve flood protection [4].

The identification of the main factors adversely affecting riverbank stability and triggering failure conditions is still a fundamental step for the optimal design and verification of these linear infrastructures, together with the individuation of those critical sectors, where structural improvements are mostly needed. The hydraulic and mechanical behavior of a riverbank during high water events can be generally investigated by a suitable combination of transient seepage and stability analyses, taking into account appropriate unsaturated soil geotechnical characteristics. In this framework, the assessment of hydromechanical parameters of the soil constituting the body of the embankment and foundation in partial and total saturated conditions is necessary in order to carry out proper predictive analyses of the collapse mechanism.

This paper describes a simple procedure for experimentally estimating optimal soil properties for riverbanks and their implementation in numerical analyses. The study case and materials are selected as representative for the riverbank systems of the alpine and Apennines tributaries of river Po (generally 5–10 m high above the ground level), which recently experienced various sudden overall collapses (Figure 2). Results hereafter discussed, besides constitutes a preliminary benchmark for the design of physical models that, in further studies, will be tested with the final aim to develop diagnostic and operational tools for the assessment of river embankments resistance to flood waves and support local governments to draw their resilience action plans.

**Figure 2.** Recent riverbank breaches occurred in: (**a**) San Matteo (MO), Secchia river, 2014; (**b**) Lentigione (RE), Enza river, 2017; (**c**) Villafranca (FC), Montone river, 2019; (**d**) Nonantola (MO), Panaro river, 2020.

#### **2. Materials and Methods**

#### *2.1. Case Study*

The geometry considered in present study has been defined as suitable for replicating a typical failure for the embankments of tributaries of the Po river and concurrently reproducible in scaled physical models. The case study is thus referred to a simplified geometry (sketched in Figure 3) where the riverside and landside are equally sloped (45◦ );

the height of the embankment, which has a simple trapezoidal shape, is 7.5 m above the surrounding level and the crown is 3.0 m wide.

**Figure 3.** Geometry of the section of the riverbank analyzed.

The embankments of tributaries of river Po are prevalently constituted by a heterogeneous mixture of sands and silt, sometimes clayey soils in the above ground part, while the subsoil frequently consists of clayey and silty deposits of floodplain environment. Similarly, the soil here selected for the earthfill is constituted by a compacted mixture of 70% Ticino sand (TS) and 30% Pontida clay (PON), while a foundation characterized by a homogeneous consolidated layer of Pontida clay is taken as reference.

The whole model is supposed to be contained in a rigid box and this leads to peculiar boundary conditions for the foundation layer, which will be detailed in Section 4.1.1.

Particular attention will be devoted to riverside slope instability due to rapid drawdown after high water stages and to the landside instability, due to a consistent saturation of the embankment body.

#### *2.2. Soil Physical Properties*

Ticino sand (TS) is a well-known Italian sand which has been subjected to intensive conventional and advanced laboratory tests [18–22]. It is a uniform coarse to medium sand made of angular to sub rounded particles and composed of 30% quartz, 65% feldspar and 5% mica [18–23]. Pontida Clay (PON) is a low plasticity kaolinitic clayey silt obtained from a quarry of fine material located in Pontida (Italy), deposited in a post-glacial lake environment. X-ray diffraction tests indicate the following mineralogical composition: 25–30% illite, 20–25% kaolinite, 20–25% quartz, 12% calcite and dolomite, less than 10% feldspar, and less than 5% chlorite.

The grain size distributions for each single component (TS and PON) and the considered mixture (TS70%-PON30%) are given in Figure 4. Following a preliminary stage of laboratory investigations, soil index and physical main properties have been determined and are listed in Table 1.

**Figure 4.** Grain size distributions of the tested materials.

γ γ **Table 1.** Average on index properties of the tested materials (γmin maximum value of soil unit weight, γmax minimum value of soil unit weight, emin minimum value of void ratio, emax maximum value of void ratio, G<sup>s</sup> specific soil density, D<sup>50</sup> mean particle size, Uc uniformity coefficient, LL limit liquid, PL plastic limit, PI plasticity index).


<sup>1</sup> Average of two measurements. <sup>2</sup> Average of three measurements.

#### *2.3. Specimen Preparation and Initial Condition*

Soils commonly used as embankment fill are compacted to a dense state to obtain satisfactory engineering properties such as shear strength and/or permeability. The degree of compaction and the molding water content needed to achieve the required compaction degree can be determined through standard laboratory tests. Above all, the optimum water content and maximum dry unit weight values of a soil can be determined under standard Proctor energy (SPE) and modified Proctor energy (MPE), which represent a quite reproducible compaction procedure and typical for embankments and retaining structures. The SPE allows to obtain a degree of compaction substantially lower than the modified maximum dry unit weight. Since more energy is applied for compaction using MPE, the soil particles are more closely packed and the optimum moisture content is lower. Higher energy results in greater shear strength, lower air voids and decreased permeability. Here, the modified Proctor compaction test has been performed on TS70%-PON30% following ASTM standard [24] using an automatic compactor. The specimen properties so obtained are reported in Table 2.


**Table 2.** Index properties and initial suction of the TS70%-PON30% specimens compacted at the optimum Modified Proctor (γ<sup>d</sup> dry unit weight, w natural water content, n soil porosity, e void ratio, Sr degree of saturation).

Compacted soils used in earth structures are often in unsaturated conditions. Since the Modified Proctor compaction test reproduces the staged construction of earth embankments that is expected to be in unsaturated conditions, thus, immediately after this test the matric suction was measured in order to check the initial conditions (Table 2).

After instrument saturation, the porous stone of a small-tip tensiometer (Figure 5a) was coated with a slurry of Pontida clay, to provide adhesion between the ceramic cup and the soil in which it was inserted. Then, with the help of a screw of the same diameter of the porous cup, a hole was created into the soil contained into the Proctor mold. The bottom of the hole laid at a distance of 5 cm from the surface of the material, then, when the porous cup was finally inserted into the soil, the superficial and disturbed layer was removed and the measurement returned a suction value that can be assumed representative of the entire compacted material contained into the mold. When the hole was excavated and the porous cup inserted, the sample and the small-tip tensiometer were hermetically covered, to avoid any input of air that would compromise the accuracy of the measurements. Considering the grain size distributions involved and the type of the instrument used, the equilibrium was reached in few hours within 100 and 200 min. The tests have been carried out with two small-tip tensiometers to verify the reliability of the measurements. At the equilibrium, the two instruments gave two overlapping trends of the suction values for all the specimens, as shown in Figure 5b. The measured suction values are in general low, ranging between 4.5 kPa and 8 kPa, mainly due to high percentage of TS in the mixture.

**Figure 5.** Suction measurements in samples after compaction: (**a**) Soil moisture Probe Parts; (**b**) Matric suction vs. Time.

#### **3. Soil Hydro-Mechanical Characterization**

*3.1. TS70%-PON30%*

#### 3.1.1. Saturated Permeability from Permeameter Test

The saturated permeability was determined through a constant head test in a permeameter by measuring the volumes of water exchanged by the saturated specimen. This test relies on the Darcy's law [25], reported in Equation (1):

$$
\Delta \mathbf{V} / (\mathbf{i} \cdot \mathbf{A}) = \mathbf{k}\_{\text{sat}} \cdot \Delta \mathbf{t} \tag{1}
$$

where ksat is the saturated permeability, A is the cross-section area of the specimen, ∆V is the volume of water that flowed through the specimen in the time interval ∆t, and i is the hydraulic gradient. The volume ∆V is equal to the average of the volume of water flowing inside and outside of the specimen for each time step.

The test was performed on saturated soil specimens using the setup presented in [26]. The equipment is composed of a permeameter in which the sample is placed, two graduated burettes, two water reservoirs, two air pressure regulators and a pressure transducer. Two porous plates (previously saturated) are placed on the top and bottom of the specimen, connecting the soil water to the drainage grooves machined in the sealing acrylic caps. Two O-rings are used to seal the caps to prevent any leakage of water. For the same reason, the acrylic caps are tied together with a well tightened screw system. Finally, filter papers for sandy soils are placed on top and bottom of the specimen to avoid the blockage of the porous plates.

The soil specimen was taken from the Proctor mold immediately after suction measurements were carried out. The specimen is contained in a cylindrical steel ring 7.21 cm in diameter and 6.08 cm in height (total volume of about 248.7 cm<sup>3</sup> ). The saturation is obtained by flushing distilled and deaerated water into the specimen. Then, a small and constant pressure head gradient (5 kPa between bottom and top) is applied to the specimens to induce a water flow in the upward direction. The saturation process is controlled by continuously monitoring the volumes of water coming in and out of the specimen; it is assumed that the saturation process is complete when a stationary condition is reached, and the incoming and outgoing flow rates become equal [26]. In order to quantify ksat, the variable E' has been plotted against the time (Figure 6), with E' defined as the volume of water discharge per unit bulk cross sectional area and hydraulic gradient (Equation (2) [25]):

$$\mathbf{E}' = \Delta \mathbf{V} / (\mathbf{i} \cdot \mathbf{A}) = \mathbf{k}\_{\text{sat}} \cdot \mathbf{t} \tag{2}$$

where ksat is equal to the slope of the linear regression in the t–E' space.

Since this type of test could be affected by random errors, due to the operator's readiness during the readings from two burettes at a time, four measurements have been made for each of the four compacted specimens, as reported in Figure 6.

The values determined in permeameter evidence a strong variability in results, despite of similar porosities. Some dispersion can be surely related to heterogeneity that can occur during the compaction among all the sampled specimens. The saturated permeability values determined for each specimen of TS70%-PON30% mixture are reported in Table 3, with the mean logarithmic (of base 10) value equal to 1.29 × 10−<sup>7</sup> m/s.

**Figure 6.** Saturated permeability measurements: (**a**) MP2; (**b**) MP3; (**c**) MP4; (**d**) MP5 specimens.


**Table 3.** Saturated permeability (ksat) of the TS70%-PON30% specimens and hydraulic parameters of the Mualem–van Genuchten (MVG) retention model.

3.1.2. Soil Water Retention Curve (SWRC) and Hydraulic Conductivity Function (HCF)

In this study, the device described by [26] was used to identify the SWRC of the mixture along the main drying branch (Figure 7).

**Figure 7.** Apparatus for the evaporation tests, modified from Nicotera et al. [26]: (**a**) sampler changer; (**b**) specimen ring and sampler holder.

A commercial apparatus (ku-pF by Umwelt-Geräte-Technik GmbH) developed for performing evaporation tests according to the simplified procedure of [27] was adopted. After completely saturating the specimen in the permeameter, it was mounted in the ku-pF apparatus (Figure 7a) for carrying out the evaporation test. The specimen has been enclosed in a metallic ring which includes two holes to house the tensiometers installed 3 cm apart (respectively 1.5 and 4.5 cm from the top). Each pair of pressure sensors is connected to a conditioning unit arranged upon each sample holder (Figure 7b) that is fitted with an electronic assembly that digitalizes the tensiometer readings at selected intervals and sends them to the data logger via a two-wire line. The connection to the data logger takes place through contacts in the support rods on the rotating carrier.

The tensiometers were saturated before each test and the calibration checked. Once calibrated, the tensiometers were installed in the holes of the saturated specimen, previously made with a special guide. The material resulting from the drilling was kept aside, weighted, and then put into the stove for the evaluation of water content. Then, the specimen enclosed in the ring was placed on a perforated base covered with filter paper while the top was covered with clingfilm and a metal cup. At this point, the specimen was placed on the closed base of the basket of the ku-pF apparatus. The basket was installed in a rotating changer arm and was weighed every 10 min on a precision scale with a resolution of 0.01 g to register variations in water content.

At first, the test was conducted with top clingfilm and top cap still mounted on the specimen to reach the hydrostatic condition inside the specimen. The equilibrium was reached when the difference between the readings of top and bottom tensiometers was 0.3 kPa (the tensiometers are placed 3 cm apart). Then, when the equilibrium is reached, top clingfilm and top cap are removed and weighted. The actual evaporation process began and lasted until the tensiometers reached the cavitation value (about 80 kPa). During the evaporation test, the water content decrease of the sample was measured at set intervals by a digital balance logged into an acquisition system. After the evaporation test, the soil specimen was moved in the Pressure Plate apparatus [28,29] to determine additional experimental points for suction values between 90 kPa and 1 MPa. The results of tests on MP3 and MP5 specimens in the ku-pF apparatus are presented in Figure 8a on the soil water retention plane. In particular, the mean measured matric suction has been coupled to the average water content estimated by the measured variation in the soil weight during the evaporation test. On the same figure the points determined in the Pressure Plate at higher suction are overlapped.

**Figure 8.** (**a**) **Soil Water Retention Curves** (SWRCs) and (**b**) **Hydraulic Conductivity Functions** (HCFs) determined by inverse analyses in Hydrus 1D (black continuous line) overlapped on measurements obtained in the ku-pF apparatus, in the Pressure Plate and in the permeameter (red symbols).

σ The Mualem–van Genuchten (MVG) model [30,31] is then considered in the interpretation of laboratory data through the software HYDRUS-1D [32], which numerically solves the Richards' equation using standard Galerkin-type linear finite element schemes. HYDRUS-1D is able to simulate the water movement in one-dimensional variably saturated media.

The relationship between the volumetric water content (θ) and the matric suction (s) in the MVG model is here described by Equation (3) [31], the relationship between the effective degree of saturation (Se) (Equation (4) [31]), and the hydraulic conductivity (k) is described by Equation (5) [30]:

$$\boldsymbol{\Theta} = \boldsymbol{\Theta}\_{\rm r} + (\boldsymbol{\Theta}\_{\rm sat} - \boldsymbol{\Theta}\_{\rm r}) \cdot \left\{ \frac{1}{\left[ \overline{\boldsymbol{\Upsilon}} + (\boldsymbol{\alpha} \cdot \mathbf{s})^{\text{n}\_{\rm V}} \right]} \right\}^{m} \tag{3}$$

$$\mathbf{S}\_{\mathbf{e}} = \frac{\boldsymbol{\Theta} - \boldsymbol{\Theta}\_{\mathbf{r}}}{\boldsymbol{\Theta}\_{\text{sat}} - \boldsymbol{\Theta}\_{\mathbf{r}}} \tag{4}$$

$$\mathbf{k} = \mathbf{k}\_{\rm sat} \cdot \mathbf{S}\_{\rm e}^{\rm l} \cdot \left[ 1 - \left( 1 - \mathbf{S}\_{\rm s}^{\rm 1/m} \right)^{\rm m} \right]^2 \tag{5}$$

where θ<sup>r</sup> is the residual volumetric water content and θsat represents the volumetric water content at saturation; α, n<sup>v</sup> and l are fitting parameters of the models while m = 1 − 1/nv; ksat is the saturated permeability and S<sup>e</sup> is the effective degree of saturation.

The model parameters were estimated via inverse analysis of the experimental measurements by means HYDRUS-1D [32]. According to Nicotera et al. [26], application of the inverse method consists of a numerical analysis simulating a lab test, and of determination of the values of the model parameters for which differences between observed and simulated flow variables are minimized; the estimated values of the parameters are those that optimize the model response. Therefore, the best-fitting parameters of the MVG model were obtained by means of inverse analysis of the evaporation tests carried out in ku-pF apparatus. The parameter vector associated with the main drying (θd, θsat, α, n, l) was determined using the software HYDRUS-1D [32], in which the objective function used to fit the data was minimized using the Levenberg–Marquardt nonlinear minimization method [33,34]. The initial condition was fixed as the hydrostatic pressure distribution estimated from initial suction measurements. The water flow occurring within the specimen was vertical, with a null flux at the bottom and an upward flux at the upper boundary equal to the weight change in the specimen in the given time range. Data sets and respective

weights in the objective function were composed of (i) matric suction measurements at the tensiometers during the monitoring time, with a weight of 1; (ii) the pair (s, θ) obtained from the pressure plate, with a weight of 1; and (iii) the pair (s, θ) corresponding to the AEV, with a weight of 1. The AEV was roughly identified as the point of maximum curvature on the SWRC curve, obtained by coupling the mean measured matric suction to the average water content estimated by the measured variation in the soil weight during the evaporation test.

Suction and water content numerically determined are overlapped to the experimental results in Figure 8a for tests on MP3 and on MP5 specimens. In Figure 8b, the permeability curve is shown where the value of saturated permeability measured is also overlapped. The fitting parameters of the MVG model determined by inverse analyses for each sample are reported in Table 3.

Lastly, it is important to keep in mind that, in this study, the hydraulic hysteresis is neglected and the hydraulic soil behavior has been modelled using a single branch (drying) of the SWRC. According to the literature [34,35], the main wetting branch and the scanning paths are located below the main drying curve on the SWR plane, these are characterized by hydraulic conductivities generally lower than those detected along drying curves. Therefore, neglecting hysteresis could cause errors in the mass flux calculations and soil shear strength determination.

#### 3.1.3. Mechanical Soil Properties at Full Saturation

Monotonic triaxial compression tests were performed to study the stress–strainstrength behavior of the TS70%-PON30% mixture. Six monotonic consolidated drained (CD) triaxial tests were carried out using different confining pressures (σ'c) on TS70%- PON30% specimens with 38 mm diameter and 76 mm height (Table 4). Specimens were compacted at the optimum modified Proctor (MPE). The weighed materials were moistened uniformly and then two modified Proctor compaction tests were done to prepare two sets of three specimens from each Proctor mold.


**Table 4.** Characteristics of TS70%-PON30% specimens used in triaxial testing.

All monotonic tests were performed in drained conditions (CD) using a straincontrolled compression loading system. The prepared specimens have been saturated by adopting base-top drainage technique in the triaxial apparatus until Skempton's B-value exceeded 0.98. After achieving full saturation, the specimens were subjected to the predefined confining pressures and then the axial load was applied up to an axial strain of about 20%. The imposed axial displacement rate was 0.05 mm per minute. During the tests, axial stresses, axial strains and volumetric strains were recorded.

Figure 9a,c show the variations of deviatoric stress and volumetric strain versus axial strain for tested material. As expected, the peak deviator stress increases with increasing applied confining pressure and the higher the value of confining pressure, the greater the axial strain corresponding at peak deviator stress. The rapid decrease in deviatoric stress after peak indicates the brittle response of the dense mixture. All specimens showed a dilating trend in their volume change, as evident also in Figure 9d. Only four specimens reached critical state condition (Figure 9a–c). The functions of the critical state line (CSL) in the q-p' and e-p' plane are also reported in Figure 9b,d while the peak strength envelope in the q-p' plane is shown in Figure 9b.

ε ε ε **Figure 9.** Stress-strain behavior of TS70%-PON30% mixture: (**a**) q-εa; (**b**) q-p'; (**c**) εv-εa; (**d**) e-p' planes.

By processing the experimental data from CD triaxial tests, the angle of shearing resistance and cohesion values obtained are 45.5◦ and 5 kPa in peak condition and 34◦ at critical state condition, respectively.

#### *3.2. Pontida Clay*

Pontida Clay has been extensively used in the past to investigate the effect of temperature on the stress-strain behaviour of the clay skeleton and clay-water system [36,37].

To reproduce the state of a foundation layer deposited in a floodplain environment, a slurry of Pontida Clay was prepared at a water content w = 43%, equal to 1.8 time its limit liquid, (distilled water and dried Pontida clay were mixed under vacuum) and one-dimensionally compressed under a maximum vertical effective stress of 200 kPa.

− Saturated permeability of Pontida Clay was derived from an oedometer test at different vertical effective stress as a function of the primary consolidation coefficient c<sup>v</sup> and the compressibility coefficient mv, according to Terzaghi's formulation of one-dimensional consolidation. Figure 10a,b show the computed permeability and the void ratio as a function of the applied vertical stress. Given the vertical overburden stress under the levee, ranging from about 100 to 150 kPa, a reference value ksat = 6.67 × 10−<sup>10</sup> m/s, relating to last loading test of the oedometer test, has been assumed for simulating the state of the levee foundation layer.

**Figure 10. Pontida clay** (PON): (**a**) saturated permeability and (**b**) oedometric curve, from oedometer test.

ϕ The mechanical properties of PON have been obtained through a series of drained and undrained triaxial tests carried out on both isotropically and anisotropically consolidated samples. The samples were obtained from large specimens, 1-D compressed from a slurry in a consolidometer. At the end of consolidation in the triaxial cell, some samples were normally consolidated, some others had an over consolidation ratio in the range 1.5–16. All the samples reached the critical state (or constant volume conditions). The state of the samples at critical state, evidenced with black squares in the q-p' (Figure 11a) and e-p' (Figure 11b) planes, have been interpreted to derive the critical state line (CSL), which has a slope M = 1.33 (i.e., a shearing resistance angle at critical state φ'cs = 33◦ ). The critical state conditions are also plotted in Figure 11b in the e-p' plane and have been interpreted via the logarithmic equation reported in the figure. ϕ

**Figure 11. Pontida clay** (PON) mechanical properties: critical state conditions in the (**a**) q-p' and (**b**) e-p' planes from triaxial tests.

#### **4. Numerical Analysis on the Hydraulic Response and Stability Assessment**

#### *4.1. Seepage Analyses*

The finite element software SEEP/W [38], capable of solving complex saturated and unsaturated groundwater flow problems, has been used to perform both steadystate and transient seepage analyses with reference to a bidimensional river embankment model, based on the simplified geometry sketched in Figure 3. The software calculation is formulated on the base of Darcy's law [25], originally derived for saturated soil. In the case of unsaturated porous media, the hydraulic conductivity is assumed to be a function of volumetric water content and matric suction [39].

The partial differential equation governing two-dimensional transient unsaturated water flow through riverbank is the Richard's mass conservation law, which states that the difference between the flux entering and leaving an elemental volume at any instant of time is equal to the change in storage of the system [39,40]. As implemented in the code, it is given by (Equation (6) [40]):

$$\frac{\partial}{\partial \mathbf{x}} \left( \mathbf{k}\_{\mathbf{x}} \frac{\partial \mathbf{h}}{\partial \mathbf{x}} \right) + \frac{\partial}{\partial \mathbf{y}} \left( \mathbf{k}\_{\mathbf{y}} \frac{\partial \mathbf{h}}{\partial \mathbf{y}} \right) + \mathbf{Q} = \frac{\partial \mathbf{\theta}}{\partial \mathbf{t}} \tag{6}$$

where h is the total hydraulic head, k<sup>x</sup> and k<sup>y</sup> represent the hydraulic conductivities in the x and y direction respectively, Q is the sink term, θ is the volumetric water content and t is the time. For the present application, soil permeability has been assumed to be as isotropic. This model assumes that the soil porous medium is rigid under partially saturated conditions (no volume deformations, Equation (6)). However, SEEP/W is able to simulate bidimensional un-coupled consolidation in the fully saturated soil by specifying the coefficient of volume compressibility. When the hydraulic load is applied at soil surface, un-coupled approach could not result rigorous and a fully coupled stress/pore-water pressure analysis might be required. Especially in rapid drawdown phenomena, the un-coupled approach could overestimate pore water pressures respect to fully coupled analysis. This reflects in the prediction of lower safety factor, in favor of safety. Nevertheless, in this study, the coefficient of volume compressibility of mixture is set equal to zero; this hypothesis seems here reasonable due to the predominantly granular texture of the compacted mixture, which is characterized by significant stiffness. Therefore, the differences in terms of positive pore water pressures, and thus factor of safety, computed by the uncoupled approach, adopted here, and by the fully coupled one are expected to be small.

Furthermore, at the present stage of the analysis, soil-atmosphere interaction processes have been neglected in seepage modelling.

Under steady-state conditions, which may be reasonably referred to extended period of stationary river level, the water storage in the model is constant in time, so the right side of the differential Equation (6) reduces to zero.

An unstructured mesh, characterized by an irregular pattern, consisting of 5878 triangular and quadrangular elements, with an approximate element global size of 0.20 m, was adopted for the entire domain as shown in Figure 12. The mesh density was iteratively optimized by reducing the element size until no significant change in pore-water pressure was observed after refinement.

**Figure 12.** Initial pore water pressure distribution for transient seepage analysis (isoline increment = 10 kPa), boundary conditions and modelling mesh. Phreatic surface in blue bold line.

Soil permeability of the foundation layer was considered constant (only saturated conditions) and equal to its saturated value, while the soil forming the embankment was assumed as partially saturated considering the average parameters listed in Table 3 and Equations (3) and (5) as retention and hydraulic functions.

#### 4.1.1. Boundary and Initial Conditions

To simulate a stationary flow through the embankment for a persistent retained high water level, a steady-state seepage analysis was performed assuming a maximum hydrometric height of 6.75 m (90% of the maximum embankment height) on the riverside, with the landside water table located at the ground level. As mentioned, since the model is supposed to be contained in a rigid box, the bottom horizontal, right, and left vertical outlines of the foundation layer were assumed to be impermeable. A no flow condition with a potential seepage face review has been assigned, instead, to the crest of the embankment and the landside slope; this implies that the water flow will remain zero until the pore water pressures computed by the code are not positive, otherwise the boundary condition will be reviewed to zero pressure head.

Concerning transient seepage analysis, the definition of initial condition in terms of suction and pore water pressure represents a crucial, preliminary, aspect. In the present study, a steady-state regime associated to a 0.75 m retained water level upon the riverside ground level and a hydrostatic suction pattern above the phreatic line, is assumed as starting condition. Suction distributions measured on site are generally different from the hydrostatic pattern, e.g., [9], especially during the wet period, nevertheless this hypothesis has been considered admissible for a preliminary model stability evaluation.

With regard to boundary conditions in transient seepage analysis, a variable hydrometric condition was imposed to model river level fluctuations on the surfaces potentially affected by water action. In order to simulate a realistic series of river flood events, data from the sudden collapse occurred on 19 January 2014 on river Secchia (Figure 2a) were taken as reference. In this particular case, the hydrograph recorded from 25 December 2013 to 19 January 2014 was characterized by a rapid series of three hydrometric peaks, which occurred in about one month [41]. Similarly, in the present analysis, a synthetic hydrograph is so characterized by a series of three hydrometric peaks. Starting from a river stage of 0.75 m in the first two days, each high water event reaches a maximum level of 6.75 m, raising in one day, then returning to the minimum water elevation in three days after a day of hydrometric peak persistence. The time between two following events is six days, while the final high-water level lasts for a longer period (two days). The water level was modelled as time-dependent boundary condition, expressed with linear variations of total hydraulic head (Figure 13) with a total elapsed time of 30 days. Moreover, rainfall infiltration fluxes have not been considered to specifically focus on the transient effects of hydrometric variations.

**Figure 13.** Numerical total hydraulic head boundary condition assumed in the transient seepage analysis.

#### *4.2. Stability Analyses*

The pore water pressure distributions obtained from seepage analysis, both in transient and steady state conditions, were used to evaluate the stability of the embankment on both riverside and landside slopes, expressed in terms of Factor of Safety. Stability analyses were performed with the software SLOPE/W [42] and adopting the Morgenstern and price limit equilibrium method [43].

Circular slip surfaces were generated by setting geometric constraints, consisting in the ranges for the entry and the exit zones of the sliding mechanisms, selected at the top and the toes of the embankment (Figure 14).

**Figure 14.** Slip surfaces investigated for the riverside (**a**) and the landside slopes (**b**). Entry and exit ranges of sliding mechanism are drawn in red.

The failure criterion of Vanapalli et al. [44], implemented in the software, was used for considering unsaturated shear strength contribution and it can be expressed in the form (Equation (7)):

$$\boldsymbol{\pi} = \mathbf{c}\prime + (\boldsymbol{\sigma}\_{\rm n} - \mathbf{u}\_{\rm a})\tan\varphi\prime + (\mathbf{u}\_{\rm a} - \mathbf{u}\_{\rm W})\mathbf{S}\_{\rm e}\tan\varphi\prime\tag{7}$$

where τ is the shear strength, σ<sup>n</sup> is the total normal stress, u<sup>a</sup> is the pore air pressure, u<sup>w</sup> is pore water pressure, c' is the effective cohesion, ϕ′ is the friction angle in terms of effective stress, and S<sup>e</sup> is the effective degree of saturation. In saturated conditions, when the pore-air pressure, ua, is equal to the pore water pressure, uw, Equation (7) turns into the classical Mohr–Coulomb failure criterion and the shear strength is a function solely of the normal effective stress. The second part of the equation provides the shear strength contribution due to matric suction (u<sup>a</sup> − uw), that is correlated to the soil water retention curve by the term Se. The soil parameters adopted in the evaluation of the factors of safety, related to the sliding mechanisms hereafter investigated, are reported in Table 5.


τ = c′ + (σ<sup>୬</sup> − uୟ)tanφ′ + (u<sup>ୟ</sup> − u୵)Sୣ

(u<sup>ୟ</sup> − u୵)

߬ σ<sup>୬</sup> u<sup>ୟ</sup> u<sup>୵</sup>

φ′

tanφ′

**Table 5.** Soil parameters adopted in the river embankment stability analyses. Sୣ

#### *4.3. Results and Discussion*

Sୣ

The results of the steady-state analysis, performed taking the maximum hydrometric height, in terms of pore water pressure and total head distributions are showed in Figure 15. The phreatic surface reaches a considerable height on the landside slope (half of the overall embankment elevation) and the body of the embankment is almost fully saturated.

**Figure 15.** Results of the steady state seepage analysis: (**a**) pore water pressure distribution (isoline increment = 10 kPa) and (**b**) water total head distribution (isoline interval = 0.5 m). The blue bold line indicates the estimated location of the phreatic surface.

The factor of safety (FoS) associated to the steady-state regime investigated is just above one (1.05), which theoretically corresponds to an incipient condition of collapse. The fact that the FoS turned out to be slightly greater than one, despite the significant increase in pore water pressure and the substantial raise of the phreatic surface, is due to the high shear resistance of the TS70%-PON30% mixture, which constitutes the embankment body. It should be noted that the analyses in steady state seepage conditions, though not being representative for the specific case, may provide useful information with regard to the existing margin of safety of the river embankment, as the related FoS might be identified with the reference threshold for a limit condition. Nevertheless, only transient seepage analysis, taking into account the effective fluctuations of the river water level, can predict the actual safety conditions of an earthen retaining structure. Figure 16 summarizes the results of the transient seepage analysis in terms of pore water pressure contours,

with reference to four significant time steps: (i) at the beginning of the first hydrometric peak (t = 3 days), (ii) after the first drawdown (t = 7 days), (iii) at the end of the second hydrometric peak (t = 15 days) and (iv) in correspondence of the end of the third peak (t = 27 days), while the hydrometric level is again at its maximum value.

27 days), while the hydrometric level is again at its maximum value.

**Figure 16.** Pore water pressure distribution at different stages during the transient seepage analysis: (**a**) at the first peak of the hydrograph (t = 3 days), (**b**) after the first drawdown (t = 7 days), (**c**) during the second hydrometric peak (t = 15 days), (**d**) in correspondence of the end of the third peak (t = 27 days). The black arrows represent velocity vectors of the groundwater flow.

− − In the light of the transient seepage results plotted above, some meaningful observations can be made. The groundwater flow predominantly develops in the body embankment, due to the remarkable difference in hydraulic conductivities between this unit and the foundation layer (1.29 × 10−<sup>7</sup> m/s and 6.67 × 10−<sup>10</sup> m/s respectively). During the first high water event, about a third of the embankment turns out to be already in saturated conditions. This outcome can be fairly attributable to the retention properties, the significant hydraulic permeability which characterizes the TS70%-PON30% mixture, as well as to the magnitude of the flooding event on the embankment.

After the first drawdown, the phreatic surface within the riverbank does not come back to its initial configuration, on the contrary it remains at a considerable height in correspondence of the toe of the embankment on the riverside. It is worth noting that, while the drawdown occurs, the groundwater flow reverses its motion, heading from the bank towards the river. The progression of hydrometric peaks results in a gradual advancement of the phreatic surface towards the landside slope, which leads to a significant increasing of the pore water pressure in the entire riverbank section.

The correlation between river stages, the associated pore water pressures and their effects in terms of bank stability is evident observing the evolutions of the factor of safety over time, both for the riverside and landside slopes, presented in the following diagram together with the hydrograph considered (Figure 17).

**Figure 17.** Variation over time of Factor of Safety computed at each time step for the critical slip surface of riverside and landside slopes, compared to the numerical water level assumed.

Following the first river stage rise, a reduction in the shear strength of the riverbank material is induced principally by the increasing of positive pore water pressures and the decrement in the resistance contribution related to soil suction. Hence, the FoS associated to the landside slope gradually decreases during the succession of hydrometric peaks, reaching a minimum value at last step of analysis.

With reference to the riverside slope, the stabilizing effect provided by the river water level prevails over the destabilizing forces generated by the build-up of pore water pressures when the water stage rises and, consequently, the FoS tends to increase. A significant decrease of the safety factor is however observed during the persistence of the hydrometric peaks, as a consequence of increasingly higher pore water pressures near the riverside slope. Then, when a rapid drawdown occurs, following a high water level stage, the embankment material is still close to a saturated condition, but the confining pressure is no longer present to counterbalance the shear strength reduction and this results in a fast decay of the safety factor.

Slip surfaces obtained by Limit Equilibrium stability analyses for both river embankment's sides are reported in Figure 18, through a graphical representation called 'safety map' [45]. Considering all the valid trial slip surfaces determined by the numerical code, these can be grouped into homogeneous ranges of FoS identified by different colored bands in order to easily identify critical zones. The red color is here associated to the lowest factors of safety, while the green represents the sliding mechanisms related to the highest FoS values.

**Figure 18.** Riverbank safety maps for riverside (**a**) and landside (**b**) slopes, isolines of pore water pressure distribution for the most critical time steps and minimum FoS attained. The most critical slip surface is drawn in white.

Since the most critical slip surfaces for both riverside and landside slopes involve the upper portion of the foundation layer, starting from the crest of the riverbank, the predominant failure mechanism for the present embankment section might be the macroinstability of the slopes.

However, as the minimum value of the safety factor is approximately equal to 1.32 for the landside slope and to 1.41 for the riverside slope, the embankment results to be in a stable condition for the whole elapsed period of time. Although both values are far from unity, the riverside zone appears quite more secure than the landside one, as suggested by the color maps. Moreover, it is noted that the minimum value of safety factors is not attained at the same time in both slopes. It corresponds to the third hydrometric peak for the landside slope and to the third short drawdown for the riverside slope, as expected.

Comparing the results in terms of factor of safety between the steady-state and the transient seepage analyses carried out, it is evident that the steady-state condition is associated to a safety margin even lower, frequently resulting in a highly conservative river embankment stability assessment and, generally, as an excessively prudential assumption.

#### **5. Conclusions**

In the present contribution, a methodological approach for the investigation of the hydraulic response and the stability assessment of fluvial embankments, based on laboratory testing and numerical modelling has been thoroughly described. The entire procedure takes its beginnings from the synthetic representation of the main characteristics of riverbanks subjected to frequent fluctuations of retained water levels, where consequent variation in the positive pore water pressure and suction distribution are typical. A simplified geometry, eventually reproducible in scaled physical models and characterized by a trapezoidal embankment lying on a finer grained soil foundation, has been so taken as reference.

The laboratory setup starts with the identification of the optimal physical properties and the retention behaviour of the unsaturated compacted earthfill material and includes the hydraulic characterization of both layers involved in the study, as well as a series of triaxial compression tests in saturated conditions. Seepage and stability numerical simulations have been carried out through a combination of FE (for pore water pressure computation) and LE (for stability evaluation) analyses, firstly considering both the stationary and transient flow conditions and then focusing on two main global failure mechanisms of the riverside and landside slopes. The importance of carrying out transient seepage analyses, in the case of succeeding river flood events, has been highlighted, since steady-state seepage analyses performed at the maximum hydrometric height recorded, frequently results in excessively conservative stability assessments. Factor of safety variations with time have been determined, with respect to all these aspects, gathering key information on the effects of river level fluctuations on the embankment stability conditions, which can turn out very useful for further detailed experimental and numerical investigations on the river embankment response.

**Author Contributions:** Conceptualization, R.V. and M.P.; Data curation, R.V., E.D., C.G.G. and D.G.; Formal analysis, R.V. and E.D.; Investigation, R.V.; Methodology, R.V., E.D., C.G.G. and M.P.; Project administration, R.V.; Software, E.D., C.G.G. and M.P.; Supervision, R.V. and M.P.; Validation, R.V. and M.P.; Writing—original draft, R.V., E.D. and M.P.; Writing—review & editing, R.V., E.D., C.G.G., D.G. and M.P. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded under the scheme for "Research Projects of National Relevance" (in Italian: Progetti di Ricerca di Rilevante Interesse Nazionale—PRIN), Bando 2017, grant number 2017YPMBWJt, promoted by the Italian Ministry of Education, University and Research (in Italian: Ministero dell'Istruzione, dell'Università e della Ricerca—MIUR).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The authors confirm that all data used during the study appear in the submitted article.

**Acknowledgments:** This work has been supported by the project "REDREEF-Risk Assessment of Earth Dams and River Embankment to Earthquakes and Floods" funded by the MIUR Progetti di Ricerca di Rilevante Interesse Nazionale (PRIN) Bando 2017-grant 2017YPMBWJt.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


## *Article* **Capillary Barriers during Rainfall Events in Pyroclastic Deposits of the Vesuvian Area**

**Ciro Sepe 1,2, Domenico Calcaterra <sup>1</sup> , Manuela Cecconi 3,\* , Diego Di Martire <sup>1</sup> , Lucia Pappalardo <sup>4</sup> , Riccardo Scarfone 5,6, Enza Vitale <sup>1</sup> and Giacomo Russo <sup>1</sup>**


**Citation:** Sepe, C.; Calcaterra, D.; Cecconi, M.; Di Martire, D.; Pappalardo, L.; Scarfone, R.; Vitale, E.; Russo, G. Capillary Barriers during Rainfall Events in Pyroclastic Deposits of the Vesuvian Area. *Geosciences* **2021**, *11*, 274. https://doi.org/10.3390/ geosciences11070274

Academic Editors: Roberto Vassallo, Luca Comegna, Roberto Valentino and Jesus Martinez-Frias

Received: 12 May 2021 Accepted: 23 June 2021 Published: 29 June 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

**Abstract:** In the present paper, the capillary barrier formation at the interface between soil layers, which is characterized by textural discontinuities, has been analyzed. This mechanism has been investigated by means of a finite element model of a two-layer soil stratification. The two considered formations, belonging to the pyroclastic succession of the "Pomici di Base" Plinian eruption (22 ka, Santacroce et al., 2008) of the Somma–Vesuvius volcano, are affected by shallow instability phenomena likely caused by progressive saturation during the rainfall events. This mechanism could be compatible with the formation of capillary barriers at the interface between layers of different grain size distributions during infiltration. One-dimensional infiltration into the stratified soil was parametrically simulated considering rainfall events of increasing intensity and duration. The variations in the suction and degree of saturation over time allowed for the evaluation of stability variations in the layers, which were assumed as part of stratified unsaturated infinite slopes.

**Keywords:** pyroclastic soils; infiltration; capillary barriers; stability analysis

#### **1. Introduction**

Infiltration processes through layered soil deposits can be deeply affected by the presence of layers with contrasting hydraulic properties. In saturated conditions, the hydraulic conductivity of an upper, finer-grained layer overlying a coarser-grained layer can be higher than that of the coarser layer; however, in unsaturated conditions, the situation can be the opposite. The coarser layer can hinder the infiltration process through reduced unsaturated hydraulic conductivity, thus promoting the formation of a capillary barrier. This principle is applied in artificial layered soil covers to develop capillary barriers that prevent infiltration and promote lateral drainage (e.g., for landfill covers [1] and slope stabilization [2]).

Layered pyroclastic soil profiles are often complex, owing to the different eruptive phases and the deposition process of volcanic materials that follows. This is the case for the pyroclastic materials spread out over a large part of the mountainous areas surrounding the volcanic centers of Campania (southern Italy). Soil covers are indeed characterized by alternating layers of ashes, pumices, scoriae, and other volcanic products, which presents significant differences in terms of the geochemical nature, texture, density, hydraulic and mechanical properties [3,4].newpage Pyroclastic soil covers of Campania are frequently involved in rainfall-induced shallow landslides. Multiple factors can trigger landslides among these materials, including the site conditions (geomorphological characteristics,

rainfall patterns, etc.) and properties (e.g., porosity, particle size distributions) of the pyroclastic soils [5–8].

The spatial variations in the textural characteristics and the related hydro-mechanical properties of the pyroclastic layers are the key factors controlling the rates of water infiltration, with the potential to form capillary barriers at the interface between layers of different soil textures (e.g., [1,9]). Capillary barriers have been considered by some authors to be one of the factors influencing the possible instability of shallow landslides [10,11]. Several studies [10–14] have shown that the differences in texture between the layers of the pyroclastic covers of Campania may affect infiltration processes and, eventually, slope instability.

The in situ observation of capillary barrier formation is a difficult task, due to the transient nature of the phenomenon and the fact that it is highly dependent on the initial soil-moisture profiles and boundary conditions, as shown by several laboratory experiments in soil columns, lysimeters, and instrumented flumes [15–19]. Numerical simulations represent an alternative to experimental measurements when examining the performance of the capillary barrier system. Oldenburg and Pruess [20] comprehensively studied the behavior of capillary barriers through the use of numerical methods. Stormont and Morris [21] and Khire et al. [1] performed parametric numerical analyses on horizontal capillary systems to investigate the effects of different parameters on the water storage capacity and hydraulic response of capillary barrier systems when subjected to different environmental loads. Numerical analyses were conducted by Rahardjo et al. [22] to determine the thickness and length of the fine- and coarse-grained materials of the capillary barrier system. Scarfone et al. [23,24] used advanced numerical modeling to extensively study the effectiveness of capillary barrier systems in introducing water retention hysteresis into the hydraulic behavior of unsaturated soils, considering a new hydraulic conductivity model at a low degree of saturation, and the complex phenomenon of soil–atmosphere interaction in the long term.

In the present paper, the capillary barrier formation at the interface between soil layers, which is characterized by textural discontinuities, has been investigated by means of a finite element model. Two layers of soil stratification belonging to the pyroclastic succession of the "Pomici di Base" Plinian eruption (22 ka, Santacroce et al. [25]) of the Somma–Vesuvius volcano have been considered in the numerical model. The stratification is affected by shallow instability phenomena [26] that are likely caused by the progressive saturation during the rainfall events; this is a mechanism that is compatible with the formation of capillary barriers at the interface between layers of contrasting hydraulic properties during infiltration. The water retention properties and hydraulic conductivity function were derived for the two layers, and one-dimensional infiltration was parametrically simulated considering rainfall events of increasing intensity and duration. Variations in the suction and degree of saturation over time allowed for the evaluation of the stability of the stratification assumed as a stratified unsaturated infinite slope.

#### **2. Pyroclastic Soils**

The territory around the town of Palma Campania, about 20 km east of the city of Naples, Southern Italy, is located in the "Conca Napoletana" (Campanian plain) at the border of the Nola and Sarno plains; it is delimited to the East by the Apennine Mesozoic carbonate reliefs and to the West by the Somma–Vesuvius Quaternary volcanic complex (Figure 1).

**Figure 1.** Geological map of the southwestern portion of the Campanian plain: (1) alluvial deposits, (2) travertine, (3) loose ash-fall deposits, (4) mainly coherent ash-flow deposits, (5) lavas, (6) debris and slope talus deposits, (7) Miocene flysch, (8) middle Jurassic–Upper Cretaceous limestones, (9) lower Triassic–Middle Jurassic dolomites and calcareous limestones, (10) fault, (11) total isopach lines (m) of the most important Mt. Somma–Vesuvius explosive eruptions. This figure is modified from [5].

On the carbonate slopes surrounding Palma Campania, relatively recent small landslides of the translational sliding type [27] have been observed among pyroclastic materials. A detailed field survey allowed for the identification of a series of recent landslides, mainly located in the "Vallone Lupici", involving the Somma–Vesuvius pyroclastic deposits. The failure surface has been identified at the interface between the grey pumice and dark scoria layers (Figure 2). This is possibly due, in large part, to the substantial variation in the grain size distributions and hydraulic properties between these two levels, such that water infiltration is affected.

**Figure 2.** (**A**) A 3D satellite image of the studied area (from Google Earth); (**B**) the detachment zone of a landslide in Vallone Lupici; and (**C**) the cross-section of a typical landslide.

The studied volcanic succession (Figure 3) includes a few-meters-thick fallout deposit consisting—from the base to the top—of a basal, high-vesiculated, white pumice layer; an intermediate level of high to moderate vesiculated grey pumices; and an upper level of denser, black scoriae. This stratigraphic sequence corresponds to fallout deposits emplaced during the so-called "Pomici di Base" Plinian eruption (22,000 yrs BP [25]) of the Somma– Vesuvius volcano, which is well described in the literature [28–31].

**Figure 3.** Vallone Lupici pyroclastic deposits—the identification of the Pomici di Base Plinian Eruption layers: BS, black scoriae; GPs, grey pumices; WPs, white pumices.

The observed variations in lithological facies in the stratigraphy of this fallout deposit well-match the changes in eruption dynamics during the progression of the volcanic activity. In fact, the reduction in the clast vesicularity and size at the top of the sequence reflects the existence of a compositional gradient in the magma chamber. Whereas the first explosive phases tap the volatile-rich top of the felsic magma body (forming the first deposited materials at the base of the sequence), the following less-explosive phases sample partially degassed mafic magmas at the bottom of the reservoir (deposited later in the upper stratigraphic levels).

Bulk samples were collected at different stratigraphic heights, consistent with variations in lithological facies (i.e., grain size, color, and texture). Moreover, relatively undisturbed whole samples were retrieved using thin-walled tube samplers, from the zone close to the failure surface, specifically at the transition between the grey pumices and the black scoriae layer.

Grain size distributions of samples vary along the stratigraphic sequence. As reported in Figure 4—according to ASTM D422—black scoriae consist of 80% sand and 20% gravel, grey pumices consist of 80% gravel and 20% sand, and white pumices consist of 70% gravel and 30% sand.

**Figure 4.** Grain size distributions of the investigated units in stratigraphic succession.

The specific weight of solid grains (γs) is 27.8 kN/m<sup>3</sup> for black scoriae, 27.2 kN/m<sup>3</sup> for grey pumices, and 24 kN/m<sup>3</sup> for white pumices. Variations in the density of the solid phase are consistent with the geochemical changes described in the literature [29], from the dense trachybasaltic-latitic scoriae of the top layer to the light trachytic pumices at the bottom of the volcanic succession.

Bulk weight measurements, carried out on at least 100 single clasts of each representative class (namely 2 mm, 5 mm, and 10 mm), showed no significant variations within each stratigraphic unit [26]. The dry specific weight (γd) increased from the base to the top of the stratigraphic succession, with lower average values pertaining to basal white pumices (γ<sup>d</sup> = 6 kN/m<sup>3</sup> ), intermediate values pertaining to grey pumices (γ<sup>d</sup> = 12.4 kN/m<sup>3</sup> ), and larger values pertaining to black scoriae (γ<sup>d</sup> = 15 kN/m<sup>3</sup> ). Intergranular porosities, determined using a microstructural analysis performed by means of X-ray tomography on the reconstituted sample [26], were estimated to be equal to 33.8% for black scoriae and 36.5% for grey pumices.

#### **3. Water Retention Properties**

The water retention properties of the pyroclastic soil under study were determined using the model proposed by Arya and Paris [32], based on the particle size distribution of each formation in the stratigraphic sequence. In the model proposed by Arya and Paris, it is assumed that the solid grains are spherical, and the pore volume is approximated to that of

cylindrical capillary tubes. For each ith particle-size class, the pore radius (r<sup>i</sup> )—associated with the mean grain radius (R<sup>i</sup> )—is related to the suction value s<sup>i</sup> according to:

$$r\_i = R\_i \left[ 2en\_i^{\left(1-a\right)} / 3 \right]^{1/2} \text{ and } s\_i = \frac{2T\_w}{r\_i} \tag{1}$$

where *n*<sup>i</sup> , *e*, *T*w, and α are, respectively, the number of spherical grains, the void ratio, the surface tension of water (*T*<sup>w</sup> = 7.27 × 10−<sup>2</sup> N/m at 20 ◦C), and a constant parameter larger than unity. Thus, for a given grain size distribution, Equation (1) allows for the calculation of the value of the suction required to desaturate a given fraction of pores. In the following, α is assumed to be constant, as consistent with the one evaluated for other typical pyroclastic weak rocks and the coarse-grained soils of Central Italy [33,34].

The Arya and Paris model was applied to two layers of the considered stratigraphic succession, namely to the black scoriae (BS) and grey pumices (GP), which form the interface of interest for the possible formation of a capillary barrier during an infiltration process. From the application of the method described above, volumetric water content values were obtained as a function of suction, which represents the behavior of the materials during wetting. These values were expressed in terms of the effective degree of saturation (*Se*), assuming a residual degree of saturation for both soils equal to zero, and were interpolated using the function proposed by Van Genuchten [35]:

$$\mathbf{S}\_{\ell} = \left[\frac{1}{1 + \left(as\right)^{n}}\right]^{m} \tag{2}$$

where *a*, *n*, and *m* are fitting parameters, reported in Table 1 for both GP and BS, along with fitting parameters proposed in the literature for similar coarse layers. Soil water retention curves for both the BS and GP layers are reported in Figure 5.


**Table 1.** Fitting parameters of the Van Genuchten model.

**Figure 5.** Water retention curves of GP and BS soils, estimated using the Arya and Paris (AP) model and the Van Genuchten (VG) model.

#### **4. Infiltration Process**

The implemented numerical model, aimed at simulating the formation of a capillary barrier at the interface between two layers, is characterised by textural discontinuities. A parametric study of infiltration was carried out by varying the intensity and duration of the considered rainfall events.

The hydraulic constitutive Van Genuchten model for unsaturated soils was implemented in the CODE\_BRIGHT finite-element code [36,37]. The numerical model was a vertical column of soil made of two layers: an upper layer, 2.5 m thick, representing the finer layer (BS) of a capillary barrier system, and a lower layer, 0.5 m thick, representing the coarser layer (GP) (Figure 6). The materials forming the two layers were each modeled by defining the soil water retention curve and the porosity. The intergranular porosities of reconstituted samples of BS and GP were assumed to be representative for each layer. The potential hydraulic effects of double porosity (i.e., intergranular and intragranular porosity) were not considered. Each of the two layers was treated as a uniform material. The hydraulic conductivity function was derived from the Mualem model [38].

**Figure 6.** The geometry of the two-layer numerical model, Van Genuchten water retention curves, and Mualem hydraulic conductivity.

During numerical simulations of one-dimensional infiltration, only isothermal liquid transport was considered, with the solid phase treated as nondeformable and the gas phase as nonmobile. Thus, constant and uniform values of temperature (T = 20 ◦C), displacements of the solid phase (u = 0 m), and gas pressure (u<sup>a</sup> = 0 kPa) were imposed.

The initial condition for the numerical analyses was a hydrostatic suction profile with *s* = 20 kPa at the bottom boundary, *s* = 50 kPa at the top, and a linear variation between these values, with null water pressure at depth equal to 2.0 m below the bottom boundary (*z*<sup>w</sup> = 5 m from the ground table). In this initial condition, the coarser layer was at a very low degree of saturation. For the bottom boundary condition, a free drainage condition was imposed. This consists of a zero water flow value if the suction is higher than 0.3 kPa, whereas a suction value of 0.3 kPa is imposed if the suction at the bottom boundary attains this value. The value of 0.3 kPa is the suction value at which the lower coarser layer GP starts draining water (see hydraulic properties in Figure 6). For the top boundary condition, a constant value of vertical water flux (the infiltration rate) was imposed. To assess the influence of the infiltration rate, different rainfall scenarios—in terms of rainfall intensity (5 mm/h and 10 mm/h) and duration (12 h and 24 h)—were considered. The selected rainfall events have been considered to be representative events at the site, as highlighted by the data shown in Figure 7, which was acquired across approximately two years of monitoring.

**Figure 7.** Rainfall events monitored at the site during the time interval from 05/2017–02/2019.

#### **5. Results and Discussion**

#### *5.1. Infiltration Process*

The infiltration analysis is carried out with reference to the suction and degree of saturation profiles, and by checking the water storage in the upper BS layer. These parameters are controlled over time, beginning at the start of the considered rainfall event.

For a rainfall event with an intensity of 5 mm/h and a duration of 12 h, it is observed that (Figure 8) the infiltration front proceeding downward causes a reduction in suction and an increase in the degree of saturation along the vertical column for periods of time that are longer than the duration of the rainfall event. When the interface is reached (depth of 2.5 m), the degree of saturation increases and reaches its maximum value in the upper layer, while it remains almost constant in the lower layer. Similarly, suction decreases at the interface in the BS layer and reaches the minimum value, while it remains at the initial value in the grey pumices. The conditions described correspond with the formation of a capillary barrier at the interface between black scoriae and grey pumices. The capillary barrier also persists for a longer time, as evidenced by both the suction and saturation degree profiles at observation time t = 480 h. Further confirmation of this permanence comes from the constant value over time of the water storage in the BS layer, as shown in Figure 9.

**Figure 8.** Infiltration process for a rainfall event with an intensity of 5 mm/h and a duration of 12 h: (**a**) the suction profile, and (**b**) the degree of saturation profile.

**Figure 9.** Water stored in the BS layer during the simulated infiltration processes.

For a rainfall event with an intensity of 5 mm/h and a duration of 24 h, the infiltration front close to the interface causes an increase in the degree of saturation (Figure 10a); however, this is coupled with suction reduction at the interface. The suction attained at the interface is so low that break-through has occurred in the capillary barrier, and water flow takes place toward the lower formation, with a progressive reduction in water storage in the upper black scoriae layer, as highlighted in Figure 9. As reported in Figure 11, break-through occurs when suction at the interface drops down from 25 kPa to 0.6 kPa, which is close to the theoretical value (s = 0.46 kPa) proposed by Lu and Likos [39].

**Figure 10.** Infiltration process for a rainfall event with an intensity of 5 mm/h and a duration of 24 h: (**a**) the suction profile, and (**b**) the degree of saturation profile.

**Figure 11.** Suction evolution at the interface between BS and GP during the rainfall event at an intensity of 5 mm/h and a duration of 24 h.

For a rainfall event with an intensity of 10 mm/h and a duration of 24 h (Figure 12), the infiltration front crosses the interface in less time than the duration of the rainfall event. The observed increase in the degree of saturation in the BS layer approaching the interface is due to the contrast in the hydraulic conductivities of the two layers. The reduction in water storage over time, beginning at the end of the rainfall event, is a further confirmation that the capillary barrier effect at the interface is rapidly lost for the simulated rainfall event.

**Figure 12.** Infiltration process for a rainfall event with an intensity of 10 mm/h and a duration of 24 h: (**a**) the suction profile, and (**b**) the degree of saturation profile.

#### *5.2. Stability Analyses*

This section refers to the results obtained from a numerical example of application, which was aimed at highlighting the hydraulic effects on slope stability conditions following a rainfall event lasting 24 h and at an intensity equal to 5 mm/h and 10 mm/h, respectively. The case study is a slope (slope angle α = 36◦ ) of black scoriae layer (ϕ = 35◦ , c' = 0, γ = 15 kN/m<sup>3</sup> ) 2.5 m thick superposed to a bottom layer of grey pumices (ϕ = 37◦ , c' = 0, γ = 12.4 kN/m<sup>3</sup> ) 0.5 m thick. The shear strength parameters of the two layers have been deduced by [5] and [40]. The water table is located at a depth of z<sup>w</sup> = 5 m from the ground table. The profiles of the saturation degree and soil suction during the infiltration events are shown in Figure 7, Figure 9, and Figure 10 at different times of observation (t = 12 h, 24 h, 240 h, and 480 h). These profiles, in turn, made it possible to estimate the soil shear strength at increasing depths, by adopting the failure criterion for unsaturated cohesionless soils proposed by Vanapalli et al. [41]. This leads to Equation (3), below:

$$\text{SF}(\mathbf{z}) = \frac{\mathbf{r}\_f}{\tau} = \left( 1 + \frac{\mathbf{S}\_r(z) \cdot \mathbf{s}(z)}{\gamma z \cos^2 \mathfrak{a}} \right) \cdot \frac{\tan \mathfrak{q}}{\tan \mathfrak{a}} \tag{3}$$

The decrease in the safety factor with depth, as the infiltration front proceeds, has been shown in Figure 13 as a consequence of the reduction in the positive contribution provided by the degree of saturation and suction in the expression (3). For a slope inclination close to the value of the internal friction angle of the shallow layer, it can be observed that for a rainfall event with an intensity (i) of 5 mm/h and a duration of 24 h, the factor of safety is progressively reduced as the infiltration front proceeds, but no instability of the slope occurs. At longer observation times, when the front reaches the interface between black scoriae and grey pumices, the formation of the capillary barrier (as described in Section 5.1 and shown in Figure 10) favors the further reduction in the safety factor—down to the instability condition—which is evident for the observation times of 240 h and 480 h. The failure surface is localized within the black scoriae layer, while the underlying grey pumices layer remains stable (SFmin = 1.07).

**Figure 13.** Safety factor (SF) vs. depth (z) for a rainfall event with an intensity (i) of 5 mm/h and a duration of 24 h, at increasing observation times.

For a rainfall event at a higher intensity (i = 10 mm/h) and a duration of 24 h, the infiltration front reaches the interface during the rainfall event, and the formation of the capillary barrier favors the onset of the unstable condition during the 24 h of a rainfall event.

A stability analysis of the infinite slope model (for stability conditions close to failure) showed how the formation and permanence of the capillary barrier can induce the instability of the shallow layer, confirming the in situ observation of the localization of instabilities in the shallow finer layers. The failure condition is triggered by rainfall events at intensity and duration levels favoring water storage in the upper layer and minimizing the beneficial contribution of suction to the stability of the slope.

#### **6. Conclusions**

A detailed field survey allowed for the identification of a series of recent landslides, mainly located in the "Conca Napoletana" (Campanian plain), involving the Somma– Vesuvius pyroclastic deposits. The failure surface has been identified at the interface between grey pumice and black scoriae layers, likely due to the contrasting hydraulic properties between these two levels, which promote the formation of capillary barriers. This hypothesis has been explored in the present paper by means of a numerical study of water infiltration, parametrically simulated considering rainfall events of increasing intensity and duration. The numerical analyses showed the effects of progressive water storage in the upper, finer-grained layer of black scoriae, which can give place to the capillary barrier at the interface with the coarser-grained layer of grey pumices. The time history of the infiltration process allowed the present analyses to highlight the capillary barrier water break-through as a consequence of progressive suction equalization between the two layers, as well as the gravity-driven downward flux of water over time. The impact

of infiltration, and the influence of the capillary barrier on the layers' stability, has been evaluated with reference to the unsaturated infinite slope model. The attainment of failure conditions of the black scoriae layer for a medium intensity, long duration rainfall event is consistent with the field observations.

**Author Contributions:** Conceptualization and methodology, M.C., G.R., D.C., L.P. and D.D.M.; investigation, C.S.; data curation, C.S., R.S., L.P. and E.V.; writing—original draft preparation, G.R. and M.C.; writing—review and editing, D.C., D.D.M., L.P., R.S. and E.V.; supervision, M.C., G.R. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Geosciences* Editorial Office E-mail: geosciences@mdpi.com www.mdpi.com/journal/geosciences

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34 Fax: +41 61 302 89 18

www.mdpi.com ISBN 978-3-0365-1875-6