The Effects of Slope Initialization on the Numerical Model Predictions of the Slope-Vegetation-Atmosphere Interaction

Abstract

Deep slope movements and, eventually, slope failure, have been often interpreted to be due to slope-vegetation-atmosphere interaction on slopes formed of clayey materials in the Italian Southern-Eastern Apennines, as reported in the literature. Such slopes are generally formed of flysch, within which clay is the main lithotype. Such clays are characterized by a disturbed meso-fabric, as an effect of the intense tectonics. The paper presents the results of coupled hydromechanical numerical analyses of the slope-vegetation-atmosphere interaction for a clay slope representative for the geomechanical scenario where such climate-induced deep slope movements have been repeatedly recorded. In the analyses, different model initialization procedures and parameter values were adopted. The comparison of the numerical results with the site data is aimed at assessing the effects of the soil-vegetation-atmosphere interaction taking place in the top strata of the slope, on the stress-strain conditions across the whole slope, and on the slope stability. The comparison between the numerical results of the analyses carried out entailing different initialization stages are intended to evaluate the influence of such a stage on the model predictions. It is found that only when the slope model initialization accounts for the slope loading history, developed over geological time, the numerical predictions get close to the site observations. In such case, the numerical results confirm that deep movements consequent to progressive failure may take place in clay slopes due to the slope-vegetation-atmosphere interaction.

1. Introduction

The scientific literature provides evidence of a recent increase in research concerning the influence of climate on the inception of landslides [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17]. Such an influence develops through the soil-vegetation-atmosphere (SVA) interaction, which mostly involves the topsoils of slopes.

The SVA interaction is the combination of the mass (liquid or gas) and the energy exchanges between the atmosphere, the vegetation, the soil skeleton and the pore fluids. The phenomena bringing about such exchanges are of either thermodynamic, or hydraulic, or mechanical or chemical nature (e.g., [18,19,20,21,22,23,24,25,26]). With regards to water, such phenomena determine either its infiltration as a liquid, or its evaporation and transpiration as a vapor. The consequent transient seepage through the soils may generate variations over time of the pore water pressures, from the ground surface to depth in the slope, and consequent changes in the effective stresses and in the available shear strength, which impact the slope stability [21,27]. On the whole, the SVA interaction at the slope top boundary, together with the processes that it determines in the slope at any depth, represent the slope‒vegetation‒atmosphere (SLVA) interaction.

The literature reports studies of SLVA interaction, mostly for cases in which this brings about either serviceability problems, or the mobilization of shallow landslides such as, according to Cruden and Varnes [28], soil slipping, earthflows and flowsliding [11,29,30,31,32,33,34,35,36,37,38,39]. However, the SVA interaction may also determine measurable effects on the state of the soils at a large depth (the depth of deep landslides according to Cruden and Varnes [28], i.e., more than 30 m), where it can give rise to soil yielding [10,14,40,41], especially if the soils at depth are weak. This is found to be the case in the Italian Southern-Eastern Apennines, either for slopes location of deep pre-existing shear bands, due to past landsliding [42,43], or for slopes formed of fissured clays [44,45,46,47], as effect of tectonic shearing [3,21,27,48]. The present paper reports a numerical modelling of the SLVA interaction for clay slopes in such a geomechanical scenario. Two-dimensional (2D) fully coupled hydromechanical (HM) analyses of a uniform slope, made of an intensely fissured clay, are presented. The slope model presented here represents a prototype of several slopes that are location of landslide activity due to SLVA interaction in the Southern‒Eastern Apennines. In particular, several model features and parameter values have been derived from the Pisciolo slope case study [3,27].

In the HM numerical analyses, the SVA interaction is preprocessed according to the FAO Penman-Monteith method [18]. This provides an estimate of the transient flux at the slope top boundary, computed as net rainfall, equal to the difference between the total rainfall, R, the runoff, RO, and the evapotranspiration rate [19,22], ET, at the ground surface.

2. The Prototype Slope

The Pisciolo slope is representative of several clay slopes of the Southern-Eastern Apennines within which the current activity of deep landslides has been found to relate to climate [15,21,27,49,50,51,52]. As such, it has served as a reference in the design of the slope numerical model.

The Pisciolo slope is located at the orographic right of the Ofanto River Valley (Figure 1a). The slope is formed of sedimentary successions deposited in a pre-orogenic marine basin (Cretaceous‒Miocene), thereafter involved in the Apennine orogenesis. These successions are formed of intensively fissured clays (scaly clay), which may interbed isolated fractured rocks. The flysch in the Pisciolo slope is either Red Flysch (RF, Figure 1b,c), or the overlying Paola Doce Flysch (PD in Figure 1b,c). At the top of the PD succession, quartz sandstones of the Numidian Flysch (N) occur. After the deposition, the sedimentary successions were folded and faulted during the Apennine orogenesis. Consequently, a NW-SE anticline is located in the slope, with RF clays at its core and N rocks at its limbs, as shown in the geological map and section (i.e., section I‒II) in Figure 1c. The anticline is crossed by a W-E normal fault, along which the Pisciolo stream flows. Several slow landslides involve both PD and RF clays, as shown in Figure 1b.

The physical and mechanical properties of the slope materials have been thoroughly discussed by Cotecchia et al. [3,48] and Pedone [53]. The clays are dilatative and characterized by peak strength parameters that go from c_p′ = 5 kPa − ϕ_p′ = 14°, to c_p′ = 20 kPa − ϕ_p′ = 18° (Figure 2), if a few higher-strength soil samples are excluded [3]. The maximum rate of the landsliding occurs in late winter and it may reach tens of centimetres per year [3]. Landslides have severely damaged both an Apulian Aqueduct pipeline (lying few metres below ground level) and a road, both located at the base of the slope (Figure 1).

Landslide bodies C9, C and A (Figure 1) have a common toe and form a roto-translational multiple landslide, which was triggered in the past by the evolution of the hydrological network and concurrent morphological changes of the hillslope, determining an increase in water infiltration, as discussed by Cotecchia et al. [3]. The landslide toe is intercepted by inclinometer I12, whose monitoring has given evidence to a shear band at about 19 m depth (Figure 1c).

The displacement rates monitored along this shear band (inclinometric measurements) and those monitored at the ground surface (through the GPS sensor S2, Figure 1b,c), are reported in Figure 3, together with the piezometric levels recorded over time by means of the piezometers at 15 m (i.e., electric piezometer) and 36 m (i.e., Casagrande piezometer) depths down borehole P7 (Figure 1b,c). The concurring fluctuations of the piezometric levels and of the displacement rates, shown in Figure 3, match a seasonal (180-day period) trend. The maximum values of the piezometric head and of the displacement rates occur at the end of winter or early spring, whereas the minimum values correspond to the end of summer. The piezometric fluctuations concur with the 180-day cumulative rainfalls, also shown in the figure [3,27,53].

Therefore, the piezometric regime in the slope and the activity of the medium to deep landslides, C9 and C (Figure 1), appear to be controlled by the yearly climate. A similar influence of the yearly climate on both the piezometric regime and the deep landslide activity has been recognized in several other slopes in the Southern-Eastern Apennines [27,48,49]. The recurrence of such observations has motivated the HM modelling of the SLVA interaction, presented in the following.

3. Fully Coupled Hydromechanical Modelling

The predictions of the effects of the SLVA interaction on the stress‒strain field in the slope may be obtained if fully coupled hydromechanical (HM) analyses are carried out [21]. Such analyses involve the integration of the momentum‒balance, together with the mass balance equations (for both the gas mass and the liquid mass) in a single system. The prediction of the displacement field allows for the assessment of the serviceability of structures interacting with slopes subjected to the SLVA interaction [21,54,55,56].

The HM numerical analyses reported in this paper have been carried out using the finite element (FE) program PLAXIS 2D 2016 [57].

3.1. Governing Equations

Galavi [58] reports the FE formulation for the HM analyses implemented in Plaxis 2D 2016. Such a formulation was developed based on the Biot consolidation theory [59]. It involves two sets of governing equations: (i) the momentum balance equations; (ii) the water mass balance equations (which use the Darcy law for either the liquid or the vapour), whose unknowns are the displacement vectors and the pore water pressures.

In the case of partially saturated conditions, the integration requires the definition of three fundamental constitutive laws: the soil-water retention curve (SWRC), the relative hydraulic conductivity function and the mechanical constitutive law. For the first and second laws, the analyses presented have implemented the widely used van Genuchten-Mualem model [60].

The mechanical constitutive law has been modelled according to the single stress variable framework that, for partially saturated soil conditions, uses the generalized formulation of the normal effective stress [61]:

$${\sigma }^{\prime }=\left(\sigma -u_a\right)+\mathsf{\chi }\left(u_a-u_w\right),$$

(1)

where χ is the effective stress parameter, function of the degree of saturation, S; ${\sigma }^{\prime }$ and $\sigma$ are the effective and the total normal stress, respectively; $u_a$ and $u_w$ are the pore air pressure and the pore water pressure, respectively; $\left(p-u_a\right)$ is the net stress; and $\left(u_a-u_w\right)$ is the matric suction, s. The laws defining the effective stress state and the soil stress‒strain behaviour are the same across the whole system, where χ is 1 for the saturated part and lower than 1 for the partially saturated part. The modelling approach used neglects both the soil collapse upon wetting and that upon net loading. Such a choice is consistent with the experimental evidence that the clays in the reference slopes do not exhibit either of such collapses, which are instead typical for coarse soils.

Experimental studies have shown that the parameter χ depends on factors such as the wetting and drying history, the void ratio, e, and the soil structure [62]. Several formulations of χ are available in the literature [63,64,65,66]. In Plaxis 2016 the χ formulation is:

$$\mathsf{\chi }=S_{eff}=\frac{S-S_{res}}{S_{sat}-S_{res}},$$

(2)

where $S_{eff}$ is the effective saturation, $S_{sat}$ is the degree of saturation for $s =0kPa$; $S_{res}$ is the degree of saturation at residual state [67]. In the present modelling, the χ function provided by Equation (2) has been calibrated to represent the fissured clay response through the selection of the $S_{res}$ value, as discussed later in the paper.

3.2. Geometry, Discretisation and Hydromechanical Parameters of the Slope Model

Figure 4 reports the FE mesh of the slope model. Between the base and the top horizontal ground levels, the numerical model has been set to have a constant slope, selected to be about the average of slopes found to be location of climate-induced landslide activity in the Southern-Eastern Apennines (12°; e.g., section III‒IV in Figure 1c). The mesh includes 136,101 nodes in 16,918 15-node triangular elements, each having nine stress points. The total model extension, L, is 1850 m, whereas its height is 425 m upslope, HT, and 200 m downslope, HB (Figure 4), after the slope initialization.

Both the left and right vertical boundaries are set to be impervious, since they represent a watershed and the centre of the river valley, respectively, according to the in situ conditions. Consistently, at both these boundaries, the horizontal displacements are set to null. The bottom boundary has been set as impervious and of null total displacement. The distance of the bottom boundary from the top of the ground surface in the slope model has been chosen after a sensitivity analysis [54], showing that both the seepage regime and the stress‒strain states in the model do not change any more, for depths of the model bottom boundary higher than 200 m below the ground level at the slope toe.

The slope model is made of homogeneous clay since, in the geomechanical reference scenario, slope failure has been found to start and progressively develop mostly through the clays part of the flysch [3,48,68]. In this respect, the numerical modelling is intended to analyse the extent to which the SLVA interaction may cause first failure in slopes made of rather weak clays of the type present in the Southern‒Eastern Apennines, disregarding the presence of pre-existing shear bands, despite these are recurrent in the slopes of the region. Hence, it is worth highlighting that the model results will provide evidence of the possible displacements caused by climate-induced first failure in the slopes of the region, whereas displacements monitored in situ are often the effect of climate-induced re-activation of failure, e.g., see the displacement rates in Figure 2 for Pisciolo. Therefore, the numerical results are meant to be compared with the real scale observations, mainly for the morphology and propagation style of the slope failure, rather than the size of the displacement vectors.

The constitutive model used to simulate the mechanical behaviour of the clay is the linear elastic-perfectly plastic Mohr-Coulomb model. All the model parameters implemented in the analyses are reported in

Table 1. Hydromechanical parameters characterizing the model clay.

Property	Symbol	Value
Saturated unit weight	γsat	19 [kN/m3]
Coefficient of saturated permeability	ksat	3 × 10−9 [m/s]
Effective Young modulus	E’	20,000 [kPa]
Effective Poisson’s ratio	ν′	0.3
Effective cohesion intercept	c’	20 [kPa]
Effective friction angle	ϕ’	23 [°]
Dilation angle	ψ	0 [°]
Unsaturated unit weight	γunsat	17 [kN/m3]
Saturated degree of saturation	Ssat	1
Residual degree of saturation	Sres	0.45
Van Genuchten parameter	gn	1.7
Van Genuchten parameter	ga	0.0095 [1/m]
Van Genuchten parameter	gl	0.5

. The strength parameter values, c’ = 20 kPa and ϕ’ = 23°, have been selected to be slightly higher than the average values recorded for the clays in the reference slopes (e.g., those recorded for the clay samples taken in the Pisciolo slope, Figure 2), in order to account for the local interbedding of coarser strata and rock inclusions within the clays in the real slopes.

Table 1 also reports the parameters for the partially saturated conditions. The unsaturated weight, ${\gamma }_{unsat}$, applies to the material above the water table, irrespective of the variability of the degree of saturation in this slope portion; as such, it has been selected to be slightly higher than the measured ${\gamma }_{dry}$ of the clay. The parameters of the van Genuchten model for the SWRC have been calibrated based upon both the results of drying tests on Pisciolo clay samples [3,27] and the need to accomplish, through the integration of the governing equations, the simulation of the clay volumetric straining measured in drying‒wetting paths. The specific way in which, in the present work, the χ function of the Bishop effective stress has been calibrated to predict the straining of the clay when in partially saturated condition is explained in the following. At the same time, the interplay between this calibration and that of the Van Genuchten parameters is clarified.

The calibration of the χ function was conceived to reproduce the volumetric straining of the clay in drying‒wetting paths above the water table, using the Mohr-Coulomb constitutive law. Such a strategy was selected in light of the knowledge of the major impact that the cited volumetric straining has on the prediction of the water mass balance over time, hence on the prediction of the water infiltration and of the SLVA interaction. Therefore, the χ function (Equation (2)) was calibrated to provide values of the invariant $p^{\prime }$, i.e., the mean normal effective stress:

$$p^{\prime }=\left(p-u_a\right)+\frac{S-S_{res}}{S_{sat}-S_{res}}\left(u_a-u_w\right),$$

(3)

consistent with the prediction (through the Mohr‒Coulomb model) of the volumetric straining measured in the laboratory during free drying tests. In particular, the calibration has concerned the parameter $S_{res}$ and made reference to laboratory drying paths starting from fully saturated conditions. The calibrated value of $S_{res}$ (Table 1) was then imposed in the Van Genuchten model. The free drying test data used for the $S_{res}$ calibration are those measured in tests on the Pisciolo clay samples [3].

The $S_{res}$ calibration did not consist of a mere fitting of the measured drying test data through the numerical modelling; rather, it was carried out according to the framework of the response of clays to drying (Figure 5) reported by Cafaro and Cotecchia [69]. Such an approach was selected because, within the vast literature concerning the drying behaviour of partially saturated soils, most contributions concern compacted soils, whereas the framework proposed by Cafaro and Cotecchia [70] is one of the few referring to the response to drying of clays deposited and consolidated in saturated conditions (at the base of a water column). Moreover, it is worth highlighting that the reference clays in the work by Cafaro and Cotecchia [69] exhibit retention properties close to those of the reference clays in the present slope modelling.

In particular, through the analysis of several test data, Cafaro and Cotecchia [69] demonstrated that the clay response to drying depends on its current void ratio, volumetric stiffness, and overconsolidation ratio. In particular, in this framework the state of an initially saturated overconsolidated clay (like the clay assumed to form the slope model, e.g., the Pisciolo clay) moves along a recompression line (in first approximation coincident with the swelling line, of gradient κ in e-ln p’ plane; Figure 5a) when it is either isotropically compressed by external loading and keeps being saturated (dashed path A-Y in Figure 5a), or when it is compressed due to s increase upon free drying, but still keeping $\mathrm{S}=1$ (Figure 5, continuous thin path A-GAE, (a) in the legend). In both cases, $\mathsf{\chi }=1$ (Equation (2)); however, in the first case $u_w$ is positive while p increases, whereas p equals zero and $u_w$ decreases to become negative in the second case. Cafaro and Cotecchia [69,70] discuss how this is also the case when $\mathrm{S}$ reduces from 1 to 0.9 with drying, i.e., the soil state becomes quasi-saturated (Figure 5, dash and dot path A-GAE, (b) in the legend). The maximum value of s by which $\mathrm{S}\ge 0.9$ is defined as $s_{des}$ (Figure 5b,c). When s increases beyond $s_{des}$, the clay starts experiencing major desaturation. Cafaro and Cotecchia defined the onset of major desaturation at $s =s_{des}$ as Gross Air Entry (GAE; Figure 5b,c), and found that an immediate limited drop in void ratio generally occurs at GAE (Figure 5b). Beyond GAE, instead, the major $\mathrm{S}$ reduction occurs with minor reduction in the clay void ratio (path GAE-B in Figure 5b,c), since the clay moves towards its shrinkage limit ($w_{sr}$) at about constant void ratio (Figure 5b). As shown by Cafaro and Cotecchia [70,71], most often for overconsolidated clays GAE corresponds to the gross yield state in isotropic compression, $p{\prime }_{yis}$, so that $p{\prime }_{GAE}$ is about $p{\prime }_{yis}$ (Figure 5a).

Summarizing, according to the framework, for s <$s_{des}$ S can drop to a limited extent, not below 0.9, and the $\mathsf{\chi }$ function must provide p’ values consistent with the prediction of a clay straining along the path A-GAE in Figure 5a, controlled by the pre-yield compressibility. For s about $s_{des}$ there is an immediate limited volumetric collapse and for $s$ >$s_{des}$, S < 0.9 and the $\mathsf{\chi }$ function must provide $p\prime$ values allowing for the prediction of zero volumetric strain increments. Therefore, in the present modelling the $S_{res}$ calibration in free drying has complied with the following conditions: (i) for s ≤ $s_{des}$, and corresponding p’ ≤ $p{\prime }_{GAE}$, the prediction of void ratio variations controlled by the parameters E’ and ν’ (Table 1), i.e., by a linear elastic compression law; (ii) for s > $s_{des}$, the prediction of zero increase in $p\prime$, keeping p’ = $p{\prime }_{GAE}$; (iii) the assumption of $p{\prime }_{GAE}$ = $p{\prime }_{yis}$ = 980 kPa, as measured for the Pisciolo clays. The value of $S_{res }$ achieved through the calibration was $S_{res}= 0.45$.

Figure 6a,b reports the free drying test data for a sample of Pisciolo clay; its behaviour has been seen to be coherent, and hence representative of several other free drying tests carried out on the Pisciolo clay samples [69]. Figure 6c shows the e-p’ data predicted by the present modelling for the free drying test and the corresponding measured data. The figure shows that HM modelling is able to provide a prediction of the volumetric straining of the clay upon drying that is very close to the measured one.

Figure 7 shows the SWRC and the hydraulic conductivity function of the HM modelling. $S_{res}= 0.45$ is slightly higher than the $S_{res}$ reached in most drying tests on the Pisciolo clay samples. However, throughout the year, the in situ saturation degree of the clay has never been measured to be lower than 0.6–0.7 [3], so that the erroneous portion of the modelled SWRC, i.e., at very low saturation degrees, is not expected to impact the mass balance predictions. Conversely, through the strategy presented above, the HM modelling predicts realistic volumetric straining at low to medium suction; moreover, at high suction, it does not predict unrealistic volumetric elastic strains, which would otherwise be predicted by adopting lower values of $S_{res}$.

At the same time, though, as is well known, the Mohr‒Coulomb model underestimates the volumetric straining of the clay when saturated and isotropically compressed to $p\prime >p{\prime }_{yis}$ [71], since the model solely implements the yielding in shearing. Also, the model does not simulate the clay strain softening in shearing. Such model features are expected to underestimate the progression of yielding and failure in the slope.

4. Initialization of the FE Slope Model

The initialization stage of the modelling is aimed at providing a distribution of the stress–strain states across the slope as close as possible to the site conditions. As such, it represents a challenging task, requiring the numerical simulation, in the framework of elastoplasticity, of the most important loading processes that have involved the slope soils in the geological history. However, with natural slopes these processes are often known only to a limited extent. In the case of engineered slopes, (e.g., of embankment slopes or slopes cut in horizontal strata), instead, the soil loading history is controlled; therefore, the numerical initialization is quite straightforward. Moreover, for the natural slopes location of deep failures, such as the slopes of interest here, the slope dimensions (i.e., height and width) are usually much larger than for the engineered slopes. Hence, the degree of uncertainty in the definition of the loading history for the different portions of the slope [72,73] is correspondingly higher. Furthermore, the dependency of the numerical predictions on the slope initialization is greater the more advanced the adopted soil mechanical constitutive law [21].

Despite an awareness of the impact of the slope model initialization procedure on the numerical predictions, the literature on such initialization procedures is limited, especially for natural slopes. Potts et al. [74] extensively discussed the effects that the model initialization procedure may have on the prediction of progressive failure, for engineered slopes resulting from cuts in London clay. They imposed oedometric (referred to as k₀ later) conditions across a horizontal clay stratum, before running an excavation stage to generate the clay slope. The authors show that the k_{0_initial} value (i.e., the ratio of the principal horizontal effective stress to the principal vertical effective stress in oedometric conditions) impacts the numerical prediction of progressive failure in the slope after excavation, in particular when using a strain-softening constitutive law. The k_{0_initial} value not only affects the morphology of the shear bands forming after excavation, but also the timing of progressive failure.

According to this premise, the modelling being presented was anticipated by a work aimed at reconstructing the loading history of several clay slopes location of landsliding in the region of interest, based upon geological surveys and field studies [3,48,49,50,68]. Thereafter, the numerical modelling was carried out for different initialization stages, either accounting or not for such reconstructions. In this way, the work is intended to provide a novel contribution to the assessment of the effects that the model initialization procedure may have on the prediction of the slope response to SLVA interaction, simulated after the completion of the initialization stage.

Table 2. Initialization stages of the numerical analyses discussed in the paper.

Analysis	Initialization Procedure	Ratio of the Horizontal to the Vertical Effective Stress (k0)	Poisson’s Ratio (ν′)
A	Gravity loading	k0 = kela = 0.428	ν′ = 0.3
B	Gravity loading	k0 = kela = 0.65	ν′ = 0.39
C	k0 procedure	k0 = 0.428	ν′ = 0.3
D	k0 procedure	k0 = 0.65	ν′ = 0.39
E	k0 procedure	k0 = 0.65	ν′ = 0.3
F	k0 procedure	k0 = 1	ν′ = 0.3
G	k0 procedure	k0 = 1.5	ν′ = 0.3
H	k0 procedure	k0 = 2	ν′ = 0.3

summarizes all the numerical analyses that have been conducted, with different initialization procedures.

The loading history of the natural reference slopes in the present work has mainly been determined by the tectonic processes relating to orogenesis and involving overconsolidated clays, and by successive valley erosion [3,75]. Therefore, the modelled slope has been assumed to be the result of lateral compression taking place during the orogenetic uplifting, followed by river erosion [3]. With the aim of simulating such a slope loading history, for one set of numerical analyses the “k₀ procedure” was adopted in the initialization stage (Table 2; analyses C to H), in which the initial ground level is horizontal and the stress‒strain conditions are oedometric, with a value of k_{0_initial} set by the user. Hence, the numerical integration has been run, in drained conditions, to reach equilibrium in the whole model, resulting in a final k₀ value depending on the selected soil constitutive law and parameter values. After the initialization stage, the profile of the slope (Figure 4) was obtained by running the drained excavation of nine soil layers, each of 25 m depth, in order to simulate the river erosion. In the excavation phases, the steady state seepage through the slope model was run, accounting for the following hydraulic boundary conditions: (i) zero pore water pressure at the horizontal ground surface upslope, coherent with the presence of springs and ponding; (ii) zero pore water pressure at the ground surface downslope, simulating the presence of the Ofanto River; (iii) a suction of 40 kPa along the sloping ground level, according to the average suction monitored in situ at the ground surface [3]; (iv) impervious condition along both the lateral vertical boundaries and the bottom horizontal one.

The k_{0_initial} values set in the analyses carried out adopting the “k₀ procedure” have been: 0.428 (C), 0.65 (E), 1 (F), 1.5 (G) and 2 (H), with ν’ = 0.3; k_{0_initial} = 0.65, with ν’ = 0.39 (D). In particular, the analyses implementing high k_{0_initial} values, F (k_{0_initial} = 1), G (k_{0_initial} = 1.5), and H (k_{0_initial} = 2), account for the high horizontal stresses caused by both the overconsolidation of the clay and the lateral compression due to tectonics. The alternative procedure used in the initialization of the slope model was “gravity loading.” This procedure sets the initial stress states by applying the soil self-weight to the slope model, which is set to have the final geometry from the start. The initial ratio of the normal effective stress on the vertical plane to the normal effective stress on the horizontal plane (which are not principal planes) is automatically set to the k₀ value in oedometric condition for an elastic material:

$$k_0=k_{0\_el}=\frac{\mathsf{\nu }\prime }{\left(1-\mathsf{\nu }\prime \right)}.$$

(4)

After the “gravity loading,” a plastic nil step was carried out to bring the system to equilibrium. The analyses adopting the gravity loading initialization procedure were run (analyses A and B, Table 2) for two different values of $\mathsf{\nu }\prime$: $\mathsf{\nu }\prime$= 0.3 (Figure 8a) and $\mathsf{\nu }\prime$= 0.39 (Figure 8b), where $\mathsf{\nu }\prime$= 0.3 is about the average for the reference clays, whereas $\mathsf{\nu }\prime$= 0.39 is not as representative. However, where the $k_{0\_el}=0.428$ (Equation (4)) corresponding to $\mathsf{\nu }\prime$ = 0.3 is quite low with respect to the site conditions in the context of the reference, for $\mathsf{\nu }\prime$ = 0.39 $k_{0\_el}$= 0.65, which is closer to reality. Values of k₀ higher than 0.65 cannot be achieved using a reasonable value of $\mathsf{\nu }\prime$, i.e., $\mathsf{\nu }\prime$ ≤ 0.4. Therefore, the gravity loading initialization could not be run for $k_{0\_el}$ higher than 0.65. The seepage conditions in analyses A and B were set to be equivalent to those set in the other analyses.

To assess the effects on the numerical predictions of the differences in the sole initialization procedure of the modelling, i.e., keeping the same values of the model parameters (Table 1), the results of the analyses with gravity loading initialization, A and B, need to be compared with analyses C and D (Table 2), respectively.

Discussion of the Initialization Stage Results

In Figure 8 the plastic points (i.e., nodes in which the stress state is at yield) at the end of the initialization stage are shown (red dots) for the different analyses (for the gravity loading, Figure 8a,b; for the k₀ procedure, Figure 8c–h). The corresponding shear strain fields are shown in Figure 9.

Neither of the analyses initialized through the gravity loading procedure (i.e., analyses A and B, Figure 8a,b and a,b) reached numerical convergence by the initialization stage. Hence, the results shown correspond to their last integration step before the loss of convergence. In both these analyses, the distribution of the plastic points and of the corresponding shear strains do not appear to configure shear bands compatible with the landsliding of a slope portion. Rather, the shear stresses mobilized by the loss of convergence are about to reach the maximum shear strength in large portions of the slope model. Therefore, in these two cases the numerical initialization provides a model with loading conditions incompatible with equilibrium, almost all the way through the model. Once initializing the model through the k₀ procedure and using the same parameter values of either analysis A or B, the numerical modelling achieved convergence. This was the case even if the predicted fields of plastic points and of shear strain in the analyses C and D were rather similar to those resulting from the corresponding analyses A and B, respectively. However, in the analyses using the k₀ initialization procedure, by the end of the initialization stage the mobilized shear stresses were lower than those mobilized in analyses A and B in large portions of the slope model. By the end of the initialization stage, the slope achieved a safety factor of F = 1.135 in analysis C and F = 1.140 in analysis D (

Table 3. Factor of safety obtained through the “c, ϕ reduction” technique.

Analysis	Initialization Procedure	Ratio of the Horizontal to the Vertical Effective Stress (k0)	Safety Factor
A	Gravity loading	k0 = kela = 0.428	No convergence
B	Gravity loading	k0 = kela = 0.65	No convergence
C	k0 procedure	k0 = 0.428	1.135
D	k0 procedure	k0 = 0.65	1.140
E	k0 procedure	k0 = 0.65	1.142
F	k0 procedure	k0 = 1	1.136
G	k0 procedure	k0 = 1.5	1.138
H	k0 procedure	k0 = 2	1.136

). All the reported safety factor values (Table 3) have been computed through the shear strength reduction technique [76,77]. It can be concluded that the strain fields predicted through the analysis using the Mohr‒Coulomb model are not so sensitive to the differences in initialization procedure, as already reported by Griffiths et al. [78]. Nonetheless, the size of the mobilized shear stress depends on the initialization procedure. Evidently, the c’ and ϕ’ values representative for the reference slopes (Table 1) are too low to guarantee equilibrium in the slope, once the initialization procedure disregards the loading history of the slope, as is the case with gravity loading.

Numerical convergence was also achieved by the end of the initialization stage for the other analyses adopting the k₀ procedure (i.e., analyses E, F, G and H), whose factors of safety are reported in Table 3. For all the analyses initialized through the k₀ procedure (analyses C to H), the values of the factor of safety are slightly higher than 1.1.

In analyses C to H, the results show that when k_{0_initial} <1, the propagation of failure is characterized by an advancing mode. It starts at the top of the slope and propagates downslope, but does not reach the toe by the end of the initialization (Figure 8c–e and Figure 9c–e). On the contrary, when k_{0_initial}≥1, failure starts at the toe and retrogresses.

However, only for k_{0_initial} = 1 (i.e., analysis F), the shear band appears to acquire a morphology compatible with the progressive development of a deep roto-translational sliding mechanism. For k_{0_initial} >1, instead, the shear band tends to deepen towards the bedrock of the slope model (Figure 9g,h). Moreover, the bigger the k_{0_initial}, the higher the deepening of the shear band. The deepening of the shear bands has been previously observed by Potts et al. [74], although when using a strain softening Mohr‒Coulomb model.

Both Figure 9d,e show that, with lower k_{0_initial}values, the highest shear strains are achieved in the upper portion of the slope, where the yielding may be largely a tension yield (i.e., controlled by tension cutoff).

In Figure 10 the stress paths during the excavation process for all the analyses implementing $\mathsf{\nu }\prime$ = 0.3 are reported for two points in the slope, M and N in Figure 4. In particular, Figure 10a reports all the stress paths for both stress points, M upslope and N downslope. The stress paths for the stress point M are also zoomed in Figure 10b. For k_{0_initial} >0.65 (i.e., analyses F, G and H), the stress paths at N follow similar trends. They are all characterized by an initial increase in deviatoric stress associated with a decrease in mean normal effective stress, which corresponds to a first stage of the valley excavation. After the stress paths approach the yield envelope, they move along it, during the unloading, while remaining at yield. It is worth clarifying that the stress paths in Figure 10 seem not to join the yield envelope, but rather follow lines at a small distance from it, merely due to a numerical tolerance adopted in the calculation, since point N in Figure 4 is predicted always to reach yield. Furthermore, the analysis predicts that the onset of yielding at N delays with reducing k_{0_initial}.

For k_{0_initial} equal to 0.428 and to 0.65 (analyses C and E), the first part of the stress path complies with an initial reduction in deviatoric stress upon the vertical unloading due to excavation. Thereafter, the approach of the isotropic stress state, the stress path starts following a trend similar to that of the stress paths for k_{0_initial} >0.65 described before.

The stress paths computed at point M (Figure 4), are chaotic (Figure 10b). In any case, point M reaches yield in the analyses C, D, F, G and H in Table 2, but not in the analysis E, although also in this case point M is surrounded by plastic points (Figure 8e). It may be concluded that yield is reached upslope, although here the shear strains are much less than at the toe.

On the whole, the analyses of the model results, reported in Figure 8 and Figure 9, show that when the unloading caused by the excavation is major, as in the case of k₀ = 1.5 and k₀ = 2, the shear strains localized in the shear band are much larger than for 0.65 ≤ k₀ ≤ 1, and the yielding is more diffuse. However, the morphology of the shear band predicted by the analyses in these cases is not such as to predispose the mobilization of a roto-translational landslide body. Conversely, the analysis with k_{0_initial} = 1 predicts the generation of a shear band predisposing the onset of roto-translational landsliding.

5. The Predicted SLVA Interaction

After the initialization stage, the hydraulic condition along the sloping portion of the model top boundary was changed from a hydraulic head condition to an infiltration condition, in order to run the SLVA interaction (Figure 11). In particular, the 365 days of net rainfall computed for the climatic year 2006–2007 were input 10 times. The net rainfall was determined as difference between the total daily rainfall and the daily evapotranspiration (Figure 11), from 1 September 2006 to 31 August 2007. The evapotranspiration rate was computed by means of the FAO Penman‒Monteith method [18] using the same procedure, input data and parameter values, extensively reported by Cotecchia et al. [27]. Only a fraction of the net rainfall effectively infiltrates through the outcropping soil layer, since runoff may occur; the latter is computed automatically by the FE code switching the infiltration boundary condition at the ground surface into an hydraulic head condition when a user-input hydraulic head is reached.

The need for modelling the partially saturated behaviour of the clay arises from the strong impact the state of the outcropping material has on the overall hydraulic balance at the ground surface. In particular, the soil water retention properties, together with the reduction in hydraulic conductivity due to the suction level, strongly modify the amount of water that can infiltrate in the soil in this calculation stage.

The pore water pressures predicted in the 10 years of analysis for the model points V and W in Figure 11, corresponding to the piezometers installed down the borehole P7 in Figure 2, are shown for the analysis F (k_{0_initial} = 1) in Figure 12. Point V is at 15 m depth and point W at 36 m depth. The results of analysis F are very close to those achieved in all the other analyses initialized with the k₀ procedure (C to H).

The numerical predictions in Figure 12 show that the modelling is successful in predicting the transient seepage connected to the SLVA interaction for the reference slopes. Indeed, the numerical modelling predicts seasonal piezometric fluctuations consistent with those recorded in situ (Figure 3) down to a large depth, despite the slope is formed of uniform clay (constant parameter values across the whole model). It is expected that a better simulation could be achieved by inputting values of both the soil stiffness and the coefficient of saturated permeability varying with depth, as generally is the case in situ.

For the analyses E to H, Figure 13 shows the shear strains cumulated during the 10 years of net rainfall, which therefore represent the strain increments due to the sole SLVA interaction. The numerical results demonstrate that such strain increments depend on the loading history rehearsed in the initialization stage of the elastoplastic modelling. The shear strain increment fields suggest that, only in the case of k_{0_initial} = 1 (analysis F), the shear strain increments localize in a shear band that represents the further development of the shear band formed by the slope toe at the end of the river erosion stage. In this case, the maximum shear strain increment is 16.3%; in the other cases, instead, the shear strain increments are much lower. Evidently, for k_{0_initial} values equal to 0.65, 1.5 and 2, the stress strain history of the soils across the slope (modelled through the initialization stage) has caused major straining and yielding in portions of the slope that do not appear to have significant further loading as an effect of the SLVA interaction. Hence, this interaction is not capable of generating important new shear straining.

Figure 14 and Figure 15 show respectively the fields of cumulated shear strain and total displacement at the end of the first, second, third, fourth, fifth and 10th years of SLVA interaction for analysis F. Figure 14f shows that the SLVA interaction causes a total maximum displacement of about 0.49 m by the end of the 10th year, with a maximum shear strain increment at the toe of the slope. However, most of both the cumulated displacement and the cumulated shear strain increment develop by the second year of SLVA interaction (Figure 14a,b and Figure 15a,b). Thereafter, the progressive failure caused by the SLVA interaction seems to slow down. This is due, on one side, to the redistribution of the stresses across the whole system with cyclic loading, but also to the elastic recovery of significant portions of the strain increments in each unloading stage. Such recovery is largely unrealistic and relates to the overestimation of the size of the elastic domain of the clay within the Mohr‒Coulomb elastoplastic model.

Following the above, the numerical results discussed so far show that it is the loading history of the slope that makes it more or less prone to becoming a location of progressive failure at large depth as an effect of the SLVA interaction. The initialization simulating such history that is found to be most suited to model the response to the SLVA interaction of slopes like those of reference in the work appears to be the k₀ procedure, assuming k_{0_initial} = 1. This is consistent with a geological history in which the soils have been overconsolidated, then subjected to significant tectonic-induced lateral loading and, later, to lateral unloading due to river erosion.

6. Conclusions and Future Research Perspectives

The study has provided evidence of the impact that different procedures of initialization of the stress state in the slope model can have on the model predictions of the SLVA interaction. As anticipated by Potts et al. [74] and Griffiths et al. [78], this impact is expected to be more important when adopting more advanced constitutive models (e.g., implementing hardening), in which case the initialization of the slope requires more care. However, the present work has shown that also when using a non-hardening elastoplastic constitutive model, the slope initialization influences the results.

In particular, the success of the numerical predictions has been shown to depend on the modelling of the main stages of the geological history of the natural slope. For natural slopes in the geological reference context (i.e., Southern‒Eastern Apennines), the gravity loading procedure should be abandoned in favour of the k₀ procedure, using k_{0_initial} = 1, in order to account for the important lateral loading that the slope soils have been subjected to in their history, and for the following unloading due to the excavation of the valley caused by the river erosion.

For all the slope models implementing the k₀ initialization procedure, the SLVA interaction has been shown to cause seasonal piezometric fluctuations down to a large depth. For k_{0_intial} = 1, these have been shown to determine a slope progressive failure capable of generating, in the long term, the base portion of the shear band of a possible deep roto-translational landslide.

The size of the cumulated displacements predicted by the modelling by the 10th year of SLVA interaction (tens of centimetres) is already such as to cause serviceability problems to infrastructures interacting at the toe. The model, though, does not predict the onset of landslides by the end of the analysis and predicts the reach of a stage of zero displacement increments, even if the SLVA interaction persists. However, such long-term predictions may underestimate the progression of failure in the slope, due to the limitations of the adopted soil constitutive model. This is because the model does not implement strain softening and overestimate the size of the elastic domain. Therefore, the research should be further developed in order to implement more advanced constitutive models for the slope soils, as well as an accurate simulation of the slope stratigraphy, which is disregarded in this work.