Numerical modelling of slope–vegetation–atmosphere interaction: an overview

The equilibrium conditions within slopes, and hence their level of stability, depend on several factors (Terzaghi 1950), which include the boundary conditions at the ground surface. These generate exchanges of water between the topsoils and the atmosphere, through processes such as rainfall infiltration, water evaporation from the soil pores and transpiration through vegetation, which, as a whole, are referred to as ‘slope–atmosphere interaction’. Changes with time of the soil porewater pressures, and hence of the soil stress–strain conditions, are consequent to such interaction. In particular, the porewater pressure changes generate, on one hand, variations of the mobilized shear strengths and, on the other, variations of the available shear strengths, and, as such, may bring about either the onset or the progression of slope failure, with eventual final slope instability. Failure results from the slope–atmosphere interaction especially within slope portions whose stability is marginal and may be lost as a result of even small increases in porewater pressure. Otherwise, failure may occur within shallow slope portions whose stability is strictly dependent on the presence of significant suctions in the soil pores (i.e. above the water table), which are reduced upon rainfall infiltration. In any case, in general the slope stability varies with time as a result of climate.

Several observations are reported in the literature of damage to structures and infrastructure interacting with slope movements connected to the slope–atmosphere interaction, which hence represents a source of risk for society. Such risk may apply to either man-made slopes or natural slopes. For man-made soil embankments and fills, made of partially saturated compacted soils and location of either road or railway infrastructure, movements due to variations in the soil degree of saturation as a result of the soil–atmosphere interaction have been widely observed to jeopardize the serviceability and safety of the infrastructure (e.g. Fleming & Taylor 1980; Alexander 1986; Hungr et al. 1999; Guzzetti et al. 2003; Schmertmann 2006; Dijkstra & Dixon 2010; Klose et al. 2014). With regard to natural slopes, the damage to buildings and infrastructure interacting with unstable slopes is increasing, mostly as a result of the increase of structures built in areas of high landslide hazard connected to climate, and hence of structures put at risk by wrong planning schemes (Cascini et al. 2005; Cotecchia et al. 2010, 2011). The safety of railways and roads is often threatened by debris flows involving soil cover on the steep flanks of infrastructure. Such flows show the highest rates and are the most damaging landslides (Cruden & Varnes 1996; Picarelli et al. 2004; Jakob & Hungr 2005; Cascini et al. 2010b; Hungr et al. 2014), and are triggered by drops in soil suction at shallow depths due to rainfall infiltration. Furthermore, despite the fact that the highest soil–atmosphere exchange flows take place at shallow depths, recent research has demonstrated also that deep slope movements found to jeopardize the safety of roads and railways may be triggered by slope–atmosphere interaction. The deep movements may relate to the variations in seepage conditions at depth consequent to the slope–atmosphere interaction, as the whole seepage domain in the slope reacts to the phenomena taking place at the top boundary. Seasonal excursions of the piezometric heads of even 2 – 3 m have been measured from 30 m down to 50 m depth in clayey slopes (Cotecchia et al. 2014) and shown to be due to the seasonal climatic processes. Such seasonal cycles in porewater pressure have been shown to bring about yielding of deep weak clays and consequent seasonal deep movements (Cotecchia et al. 2015). Furthermore, in slopes that are the location of old landslide bodies of marginal stability, seasonal cycling of porewater pressure due to climate has been shown to trigger seasonal reactivation of sliding (Cotecchia et al. 2008). Also, seasonal reactivations of damage to buildings and infrastructure interacting with such slopes have been recorded.

The work presented in this paper aims at contributing to the dissemination of knowledge about the modelling strategies that can be employed to assess the effects of slope–atmosphere interaction on the stability of engineered slopes and has been conducted within the COST programme ‘Impact of climate changes on engineered slopes for infrastructure’. The assessment of the processes activated in the soils by climate is a necessary background for the evaluation of the instability risk for slope–structure systems and has to be based on interdisciplinary analyses. These entail knowledge in various scientific fields, mainly soil hydromechanics, hydrology, meteorology, agriculture and thermodynamics. The source of the disciplinary broadness of the subject and of the complexity of the modelling is the coupling of several phenomena taking place in the slope as a result of the soil–vegetation–atmosphere interaction. Over time, the capacity of models to simulate such phenomena has increased. The paper presents a review of the possible strategies to model these phenomena, from the variation of the piezometric conditions, to the evolution of the slope stability with time, and the corresponding soil displacements across the slope.

The paper reports first an outline of the phenomena to be treated in an advanced modelling, with a review of the possible mathematical formulations of the processes, according to different levels of advancement. In this way a methodological framework is portrayed as the basis for the classification of modelling strategies of different levels of accuracy. Several case histories are then considered as examples of the application of the different modelling strategies.

The stability of either natural or artificial slopes depends on thermo-hydro-mechanical processes taking place in the soil, which are connected to both the climatic and vegetation conditions at the ground surface (Fig. 1). The climatic factors representing the atmospheric conditions are, primarily, rainfall, relative humidity, temperature, net solar radiation and wind speed, which, together with the vegetation, determine the top boundary conditions for the seepage taking place through the soil pores. When the water table is below the ground surface, the soils above the water table may be in unsaturated conditions. Hence, seepage phenomena within the soil pores in the slope may involve both liquid and gas transport and may be also affected by thermodynamic processes taking place within the pores. The atmospheric conditions vary with time and, as such, determine a variable boundary condition that causes variations of the porewater pressure distribution across the whole slope and, in turn, variations of the available soil strengths and slope stability (Lu & Likos 2004; Gens 2010; Fredlund et al. 2012; Lu & Godt 2013).

The main variables defining the thermo-hydro-mechanical state of the slope soils, assuming the porewater to be fresh (with either zero or constant low solute concentration) or disregarding any chemical process, may be set as (Table 1) the pore liquid pressure (scalar P_l), the pore gas pressure (scalar P_g), the temperature (scalar T) and the solid phase displacement (vector u). The values attained by these variables are essentially controlled by the balance equations (Olivella et al. 1994; Olivella 1995; Gens 2010; see Table 1): (1) mass balance of liquid, variable P_l (kPa); (2) mass balance of gas, variable P_g (kPa); (3) internal energy balance for the medium, variable T (degrees); (4) momentum balance for the medium, variable u (m). Because the thermal, hydraulic and mechanical processes in the soil are coupled, the variables of the different balance equations are related to each other and the balance equations should be solved accounting for such coupling. Nevertheless, different levels of coupling are accounted for in the different modelling strategies, resulting in different levels of accuracy in the assessment, or prediction, of the slope conditions.

The physical variables involved in the processes are controlled by physical laws, such as the law expressing the liquid and the gas transport, the diffusion law, the change of phase law and the constitutive law of the soil skeleton. In the theoretical framework of the hydraulics of equivalent porous media, the laws that may be used to simulate some of the basic processes involving the pore fluids are Darcy's law, to describe either the liquid or the gas advective flux, and Fick's law, for the vapour and the air non-advective fluxes. Furthermore, Henry's law may be used to predict the air dissolution in the liquid phase, the psychrometric law to correlate suctions and relative humidities at a given temperature and predict the transition of water from the liquid to the gas phase, and Fourier's law to describe conductive heat fluxes through either the fluid or the solid skeleton (see Table 1). Moreover, soil constitutive laws are required in the integration, which are typically: the soil water retention curve (SWRC) expressed as the variation of the volumetric water content, θ, or the degree of saturation, S, with the soil porosity, n, and the soil suction, s; the hydraulic conductivity function, k(s, n, T); the mechanical constitutive law; the thermal conductivity law, which depends on T, S and n; the air conductivity law; and the diffusion and the dispersion coefficients, both dependent on T, S, n and P_g.

In the framework of continuum thermo-hydro-mechanics and equivalent porous medium, the numerical modelling will require a spatial discretization of the slope (Fig. 1). Therefore, the variables, or unknowns, of the analysis are obtained as simulation output at a discrete number of points. The following should be noted:

1. Usually the shallow layers, just below the ground surface, which are location of partial saturation conditions and host the vegetation roots, require finer discretization due to the more highly variable gradients with time triggered by the slope–atmosphere interaction (as shown in Fig. 1).

2. The variability of the slope soil properties, assessed by means of site investigations and laboratory studies, has to be included in the numerical model through appropriate discretization strategies.

3. Three-dimensional analyses may be required, especially when the hydraulic, morphological and mechanical features of the slope are spatially highly variable in the out-of-plane direction.

The boundary conditions for the balance equations in Table 1, at the lateral and base boundaries of the model (see Fig. 1) may be as follows: (1) the liquid pressures and/or fluxes; (2) the gas pressures and/or fluxes; (3) the temperature and/or thermal flux; (4) the solid displacements and/or forces.

At the ground surface, the rainfall (R) represents a positive flux, whereas the evapotranspiration (ET) and the runoff (RO), resulting from the combination of the climatic variables, vegetation and soil state and behaviour, are negative fluxes. The difference R − ET − RO represents the infiltration water flux at ground surface (Fig. 1). As will be discussed below, whereas the rainfall will be always an input flux of the model at the ground surface, the evapotranspiration flux may be a result of the model simulation, if thermohydraulic coupling is accounted for, while it has to be introduced as an input boundary condition in the analyses disregarding such coupling. As specified in the section on thermohydraulic modelling, thermohydraulic analyses require the characterization of the slope material properties also in terms of thermal conductivity and volumetric heat capacity and account for the direct interference in the thermohydraulic conditions of the soil from all the climatic variables, such as air temperature, net solar radiation, wind speed, relative humidity and rainfall. However, to date not all the climate effects can be yet modelled by simulating the specific physical phenomena taking place at the top boundary of the slope and some empirical laws, or some simple assumptions, have to be used (e.g. for the effects of wind). Common to most models, either solely hydraulic or coupled, runoff processes are modelled in a simplified way (Blight 1997). They are often assumed to occur when the rainfall rate exceeds the saturated permeability of the soil (Krahn 2003). At this stage, the very top porewater pressure is set to zero, hence the top boundary condition turns into a pressure value from being a flux value. Another established approach assumes that runoff occurs when the very top porewater pressure becomes zero (Smith et al. 2008). Also in this case, at the zeroing of the top porewater pressures the top boundary condition starts by being provided as a pressure value, until the very top porewater pressure becomes negative. In both cases, the onset of ponding can be allowed by specifying positive values of porewater pressure corresponding to the desired height of the pond (Smith et al. 2008). However, a more advanced hydrological assessment of the runoff flow may be implemented in the modelling at the ground surface. Probabilistic analyses may offer a rational way for engineers to implement the variability of the meteorological conditions at the top boundary of the slope model within hazard analyses (e.g. Li & Lumb 1987; Griffiths & Fenton 1993; Rezaur et al. 2002).

In the following, available numerical strategies to model the effects of the slope–atmosphere interaction on the thermal, hydraulic and mechanical state of the slope are discussed, which are of different levels of complexity. The discussion will proceed from the simplest to the most advanced ones. The simplest integrate only part of the algorithms simulating the processes taking place in the slope and may even disregard their coupling. This is the case for the hydraulic models, which simulate the transient seepage in the slope evolving with the atmospheric conditions by accounting solely for the fluid mass-balance equations and disregarding the effects of the variations in temperature within the soil and the deformation of the soil skeleton. The models that take account also of the thermal processes within the soil and of the thermohydraulic coupling are more accurate than the hydraulic ones and will be presented afterwards. Finally, the models that account also for the coupling of the fluid mass balances with the soil skeleton deformations will be discussed.

The assessment of the slope stability conditions with time requires the prediction of the variations of the porewater pressures across the slope in relation to the atmospheric conditions. If the porewater pressure distribution is derived based upon the liquid and gas mass-balance equations (Table 1), but is not coupled with both the momentum and the energy balance (and hence does not account for the skeleton deformations and the effects of the variations in temperature within both the soil skeleton and the pore fluids), the modelling will be defined as purely hydraulic. In this case, all the terms referring to the mechanical and thermal behaviour of the soil skeleton and the thermal behaviour of either the liquid or the gas phase are not accounted for.

Hence, in the hydraulic modelling the porewater pressure computation results from the analysis of the transport processes (Darcy's law, Table 1), the diffusion processes (Fick's law, Table 1) and the dissolution of air in water (Henry's law, Table 1). The liquid and the gas densities are assumed to vary exclusively with the liquid and the gas pressure; that is, isothermal conditions are considered (Freeze & Cherry 1979; Fredlund & Rahardjo 1993; Gens 2010).

Further simplification can be introduced in the hydraulic modelling by disregarding the gas balance (e.g. by assuming the gas pressure as constant) and neglecting both the water vapour diffusion within the gas phase and the dissolution of air in water. In this case the model simulates only the liquid mass balance. If also the liquid density is assumed to be constant, the liquid mass-balance equation turns into a volume-balance equation of liquid water, which is commonly known as Richards’ equation (Richards 1931).

Given the uncoupled approach described above, in the hydraulic modelling the evapotranspiration of water from the soil must be an input boundary outflow at the ground surface. Therefore, evapotranspiration has to be estimated using phenomenological interpretations and semi-empirical laws, which must relate the climatic variables to such outflow and have to be calibrated based upon the outcropping soil conditions and the vegetation typology and state. Several approaches have been proposed to estimate evapotranspiration (Penman 1948; Thornthwaite 1948; Turc 1954; Monteith 1965). The FAO Penman–Monteith method (Allen et al. 1998) is at present one of the most commonly adopted to provide physically based estimates of the evapotranspiration rate. It uses empirical algorithms to represent the different physical processes that combine in giving rise to the outflow, accounting for the outcropping soil, the vegetation typology and the plant conditions during the year.

The hydraulic properties of the soils are the main internal parameters of the hydraulic model. For the top partially saturated soils, both the water retention curve and the hydraulic conductivity function must be implemented in the model, once derived from laboratory testing. A wide review of water retention curves and hydraulic conductivity functions was provided by Leong & Rahardjo (1997a,90b), but several have been the recent developments of the algorithms suited to represent such soil hydraulic properties (Wheeler 1996; Wheeler et al. 2003; Li 2005; Sun et al. 2007; Khalili et al. 2008; Nuth & Laloui 2008a; Pedroso & Williams 2010; Cafaro & Cotecchia 2015), which are not strictly empirical, but rather also based on a theoretical interpretation of the phenomena. It is worth highlighting that, when the Richards’ equation is used, the water retention curve is expressed in terms of the volumetric water content, θ(s) = Sn. In this case, the soil deformability inherent in the drying-wetting process and causing the variation in porosity, n, that contributes to the variation of θ upon drying–wetting is taken into account; hence the integration accounts for the soil skeleton deformation caused by the variations in suction. However, it disregards the strain compatibility across the soil system and the possible consequent redistribution of stresses and variations in stress–strain conditions (Tsaparas et al. 2002; Calvello et al. 2007; Cascini et al. 2010a; Tommasi et al. 2013; Cotecchia et al. 2014), hence it does not implement the full hydromechanical coupling. Below, the Pisciolo case study provides an example of hydraulic modelling of the evolution of the slope conditions and stability with time, implementing the ingredients noted above. When the SWRC is implemented by inputting the function S(s) and the porosity is assumed constant, the seepage problem is integrated totally disregarding the soil skeleton deformability (even the deformations upon drying–wetting).

Once the hydraulic regime within the slope is realistically simulated, the corresponding pore pressure distribution with time can be employed as input in slope stability analyses, performed, for instance, using limit equilibrium (Abramson et al. 1996), to derive the variation with time of the slope stability factor.

The variation in slope stability factor for different climatic conditions (e.g. in winter and summer) may be derived even using the results of steady-state seepage analyses. These represent even more basic approaches to estimate the differences in the slope seepage for different climatic conditions (e.g. in summer and winter). They simulate a steady-state slope seepage representative of the conditions at the time of year of interest, without exploring the transient seepage that takes place throughout the year. If saturated conditions are assumed to apply above the water table, the model will have to implement hydraulic boundary conditions (either as porewater pressures or as fluxes) able to provide a piezometric regime consistent with the regime at the time of interest in the analysis. In this case, the model will be realistic if piezometric data from the field monitoring are available and are used to calibrate the model boundary conditions to generate piezometric predictions across the slope close to the monitored data. Hence, the seepage modelling becomes a back-calculation of the surveyed slope seepage. Such an approach may be useful, for example, when comparing the stability of the slope at the end of the rainy season and of the dry season, once seepage is back-calculated for both stages (see the Volturino case history presented below). The approach is even more accurate if steady-state seepage models that implement partially saturated soil conditions above the water table are considered.

Strength parameters are input for the limit equilibrium analyses. When dealing with the stability of slope portions that are mostly partially saturated, the strength of the unsaturated soils can be represented using the Mohr–Coulomb failure criterion modified to account for partial saturation (Fredlund et al. 1978), as was done by Ng & Shi (1998), for example. They used the steady-state hydraulic approach to simulate the seepage in the slope and employed the revised version of the failure criterion and Bishop's simplified limit equilibrium method (Bishop 1955) to assess the variation in stability of a typical unsaturated cut slope in Hong Kong during different representative rainfall events.

The modelling is defined as thermohydraulic when the equations describing the mass balance of both the liquid and the gaseous phase are solved together with the equation describing the energy balance (e.g. Krahn 2003). Generally speaking, such a modelling approach allows the prediction of the liquid pressure, P_l, the gas pressure, P_g, and the temperature, T. However, either because of shortage of input data or to reduce the complexity of the numerical formulation, P_g is often assumed to be constant (i.e. equal to the atmospheric pressure) and the mass-balance equation of the gaseous phase is disregarded.

With T being calculated at each node, the thermohydraulic modelling has the advantage of predicting the water transition from the liquid to the vapour phase, which takes place when the latent heat of vaporization is supplied to the soil as a result of heat exchanges with the atmosphere. This means that the evapotranspiration fluxes can be derived from the numerical simulations, provided that all the climatic variables are defined (i.e. air temperature, net solar radiation, wind speed, relative humidity and rainfall), together with the main characteristics of the vegetation on the slope, and employed in the calculation of the liquid, gas and heat fluxes at the top boundary of the model. Boundary conditions need to be defined also on the lateral and the bottom boundaries of the mesh, either as liquid, gas and heat fluxes, or directly as imposed pressures (P_l and P_g) and temperatures (T).

Similarly to the transient seepage computed through the hydraulic modelling, for the transient thermohydraulic processes all the nodal variables need to be initialized. However, in this case, the initial conditions of the slope have to be defined not only in terms of liquid pressures, P_l, but also in terms of temperatures, T (and, eventually, in terms of gas pressures, P_g). To date, the required comprehensive characterization of the initial slope conditions (i.e. also in terms of temperature and gas pressure) still poses serious difficulties in the use of this modelling approach in practice.

The solution of the energy-balance equation introduces additional complexity, as the thermal conductivities and the volumetric heat capacities of the different slope materials also become essential model ingredients, whose definition requires additional tests (e.g. Woodside & Messmer 1961; Midttømme & Roaldset 1999). Moreover, because of the relevant role played by the unsaturated soil behaviour within the analyses, the variation of both the thermal conductivities and the volumetric heat capacities with the soil volumetric water content should also be defined.

Performing non-isothermal analyses also gives the possibility of introducing the dependence of all the material properties on T. The hydraulic conductivity, for instance, can be defined as a function of both void ratio and temperature (e.g. Lambe & Whitman 1969; Delage et al. 2000).

Because of the large amount of input data required and the model complexity, 2D thermohydraulic modelling of the slope–atmosphere interaction has rarely been performed so far (e.g. Cotecchia et al., in prep.), whereas 1D modelling is more frequently discussed in the literature (e.g. Rajeev et al. 2012).

The assessment of the slope displacements resulting from the slope–atmosphere interaction is crucial for the evaluation of the serviceability of structures interacting with slopes and requires the performance of mechanically coupled numerical analyses (Zienkiewicz & Taylor 1989; Potts & Zdravkovic 1999, 2001; Zienkiewicz et al. 1999). Such analyses entail the solution of the momentum-balance equation complying with the variation in time of the pore liquid and gas pressures. If the mass-balance equations and the momentum balance are solved at the same time as a single system, the analysis complies with a full hydromechanical coupling. If the solution of the mass-balance equations is carried out separately from that of the momentum balance, by a staggered approach, full coupling is not pursued and the accuracy of the results will depend on the size of the calculation step.

Irrespective of the adopted solution approach, any hydromechanical simulation requires the definition of the initial state of the soil in the slope, resulting from its geo-hydro-mechanical history. Hence, the initialization of the stress–strain conditions across the model represents a very challenging task, which can be rarely accomplished by simply applying an initial gravitational load. In particular, the model prediction will depend on the initial conditions especially when more advanced soil constitutive laws are adopted. Elasto-plastic laws are the most commonly used to reproduce the soil behaviour and have also been extended to partially saturated conditions.

In the Volturino case history presented below (Lollino et al. 2010, 2016), hydromechanical finite-element analyses have been performed adopting an elasto-plastic Mohr–Coulomb model for a slope initialized with the K₀ condition followed by excavation. These analyses are carried out for two permanent hydraulic steady-state conditions, respectively representative of the average condition for summer and winter. When the model is set to investigate the effects of the slope–atmosphere interaction with time, instead, the numerical simulations have to reproduce the stress–strain evolution associated with the pore pressure variations. In this case, the analysis of the deformation processes has to be coupled with the analysis of the transient seepage processes induced by the climatic conditions. The simplest strategy to represent these effects implements at the ground level, as a boundary condition of the model, suction variations close to those really occurring in situ in the topsoils as an effect of the climatic perturbation. Kovacevic et al. (2001) and Nyambayo et al. (2004), among others, adopted this modelling approach to investigate the stability with time of typical vegetated railway embankments in the UK, developing fully coupled hydromechanical finite-element simulations, but assuming the soil to be saturated above the water table. Those researchers employed an elastic strain-softening Mohr–Coulomb model, capable of simulating the key features of the mechanical behaviour of the materials forming the investigated slopes.

More recent studies have reported the results of fully coupled hydromechanical simulations in which the unsaturated hydromechanical soil behaviour is also reproduced (Fredlund & Rahardjo 1993; Wong et al. 1998; Smith 2003). The formulation of the basic equations originally developed for saturated conditions is modified to model the elasto-plastic deformations under partial saturation conditions. Moreover, the unsaturated formulation should also accommodate the definition of both the hydraulic conductivity and the water retention curve. The latter is here necessarily defined in terms of degree of saturation, as the soil deformations evolve with time and control the soil porosity. This approach has been adopted by Rouainia et al. (2009) to develop fully coupled modelling of the effects of future climatic scenarios on the stability of cuttings in London Clay, as discussed below. In this case, the pore pressure variations at ground level, in the unsaturated soils, have been calculated through an advanced modelling of infiltration, accounting for rainfall, runoff and evapotranspiration, all input at the model top boundary.

One of the first and still most commonly used constitutive laws for unsaturated soils was presented by Alonso et al. (1990). This law, known as the Barcelona Basic Model (BBM), is formulated using the first approach for the stress variables. It has been employed in the hydromechanical modelling used by Askarinejad (2013), reported in the case histories section below, to evaluate the stability of the banks of the River Rhine (Switzerland), with reference to both 2D and 3D conditions.

Several other constitutive laws, either stemming from the BBM or sharing the same choice of stress variables, have been presented in the literature (Wheeler & Sivakumar 1995; Chiu & Ng 2003; Sheng et al. 2008). However, recent constitutive models have been formulated in accordance with the second approach (Pereira et al. 2005; Russell & Khalili 2006; Khalili et al. 2008; Khalili & Zargarbashi 2010).

It is worth observing that advanced constitutive laws should be employed together with advanced formulations of the water retention curves (e.g. Gallipoli et al. 2003; Tsiampousi et al. 2013a; Cafaro & Cotecchia 2015), to take into account the intrinsic coupling between suction, void ratio, degree of saturation and soil mechanical properties. Such an approach has been proposed by Tsiampousi et al. (2013b) and employed by Pedone et al. (2016) to investigate the stability of a slope interacting with both hydraulic and transport infrastructure.

Desiccation cracking is the result of shrinking in clays and is commonly considered to initiate when, as a result of the volumetric contraction of the top strata where the main effects of evapotranspiration are observed, the tensile stresses exceed the soil tensile strength, which is generally very low. The morphology of observed soil cracking is generally found to depend on the vegetative cover (i.e. root reinforcement, water uptake potential, etc.), the soil plasticity index, the soil mesofabric and the drying conditions (potential evapotranspiration magnitude and duration). With changes in soil moisture, the soil can crack in a brittle mode under dry conditions, after a linear elastic stage, or in a ductile way, when the moisture content is high and the soil is soft (Hallett & Newson 2005).

The ability to capture the generation of discontinuities effect of desiccation is crucial for an accurate simulation of the slope response, as cracking reduces the soil shear strength and increases the water infiltration. The development of coupled hydromechanical models that include cracking in clayey soils represents a challenge for numerical modelling, given the difficulties inherent in the mathematical simulation of crack generation and flow through the cracks. The presence of both discrete deep cracks and extensive shallow crack patterns is generally the source of a more rapid response of the seepage domain to climatic changes, which may ultimately lead to slope failure.

Research into the modelling of the effects of cracking has increased in recent decades. Approaches that do not explicitly predict the cracking process but implement the variation in soil properties due to cracking include the estimation of the crack evolution with both depth and time (Fredlund & Rahardjo 1993), bimodal SWRCs, hydraulic conductivity and water storage functions (Fredlund et al. 2010) for the soil and the cracks, and have been applied also to full-scale slopes by Booth (2015).

Other research has focused on capturing the processes that initiate cracking, but does not attempt to incorporate discrete cracking into full-scale slope geometries. Researchers have followed two alternative numerical strategies, modelling cracking in either a continuum or a discontinuous medium. Both finite-element and finite-difference methods (FEM and FDM) have been employed and represent the simplest way within traditional slope stability modelling to reproduce the main effects of cracks. Meshing in FEM and FDM typically limits the ability to capture crack localization and propagation. Therefore, researchers have developed novel schemes of mesh fragmentation and discretization (Sanchez et al. 2014; Stirling 2014; Pouya et al. 2015). Alternatively, the extended finite-element method (XFEM), where additional degrees of freedom are included to represent cracks within finite elements, can be used for fluid-driven cracking, such as desiccation (Mohammadnejad & Khoei 2013), although at present the application of this strategy to soil desiccation is limited.

The application of linear elastic fracture mechanics has also been explored (Lachenbruch 1961) and included in a modelling strategy that uses a stepped approach to predict the moisture profile (with 1D FEM), then crack depth, through linear elastic fracture mechanics, and spacing, through linear elastic FEM (Konrad & Ayad 1997).

Alternatively, the discrete element method (DEM) has been used to simulate the medium subjected to cracking as an assemblage of single, interacting grains subject to contact laws (Peron et al. 2009; Amarasiri et al. 2011; Muslelak & Sliwa 2012). This approach allows the formation of discontinuities as a result of grain bond breakage and can even capture fluid flow through the cracks. However, DEM is limited in its ability to simulate the multi-phase processes within the matrix before cracking; that is, drying is often inferred through prescription of porewater pressure gradients, rather than by applying moisture variations, or heat or mass exchanges (Bui et al. 2015).

Vegetation has to be modelled when dealing with slope–atmosphere interaction, because it affects the water balance of the slope. As mentioned above, transpiration processes occurring in the topsoil layers result in a negative flux. Furthermore, leaves and steams shield the slope, intercept water and limit the infiltration of rainfall.

However, transpiration processes, although increasing the slope stability (Glendinning et al. 2009), may induce cyclic movements, as they occur on a seasonal basis. This cyclic loading can be detrimental for serviceability if it occurs in the vicinity of infrastructure. Moreover, when the cyclic pore pressure variations are pronounced and prolonged, they can also promote progressive failure phenomena (Take & Bolton 2011).

All these aspects can be predicted by advanced coupled numerical modelling, if a boundary condition simulating the presence of vegetation is appropriately implemented. This is the case in the numerical simulations presented by O'Brien et al. (2004), who modelled the hydromechanical behaviour of railway embankments covered by either grass or trees. Transpiration processes induced by grass were simulated by applying suctions at the top boundary of the model, whereas the presence of trees was simulated by applying an internal boundary condition acting in the area of the slope where the roots of the trees are concentrated. The numerical modelling was validated by means of a comparison with the available track monitoring data and proved that significant deformations can occur as a result of the seasonal pore pressure variations associated with transpiration. Similarly, Lees et al. (2013) have shown how the mechanical response of an old clay embankment was strongly influenced by the transpiration processes of some trees located in the lower part of the slope. The presence of the trees was modelled by applying the suctions induced by their roots and the results were verified by comparison with inclinometric and piezometric measurements.

The direct application of suctions as an internal boundary condition represents a simplified way of simulating the effects of transpiration and may lead to unrealistic porewater pressure predictions. Otherwise, root water uptake models should be employed (Nyambayo & Potts 2010), in which the transpiration fluxes represent the input ingredient, so that pore pressure variations can be estimated and not imposed.

It is worth observing that vegetation can also induce permeability variations of the topsoil layers, because high suction values can trigger the formation of cracks (as discussed above), which increase the soil permeability (Li et al. 2016). Similarly, if the vegetation has been removed from the slope, the outcropping soil layers should be considered more permeable, because the holes created by the roots will act as preferential pathways for the infiltration of rainwater (Smethurst et al. 2015).

Advanced numerical modelling should also account for the mechanical characteristics of the vegetation roots. This is because the tensile strength of roots within the soil mass can improve the capacity of the soil to resist the mobilized shear stresses. The maximum tensile strength, or pull-out resistance, of the roots together with their size and distribution with depth are the ingredients of the modelling of the reinforcement effects of vegetation for slope stability analyses. In particular, the experimental data obtained from direct shear tests performed on blocks of soil containing roots have shown that the presence of vegetation produces an increase in soil cohesion, leaving the friction angle unchanged (Wu et al. 1988; Faisal & Normaniza 2008). Such a mechanical effect is usually introduced in the Mohr–Coulomb failure criterion through an ‘apparent cohesion’ term, c_R, which adds to the soil effective cohesion, c′ (Waldron 1977; Wu et al. 1979).

Wu et al. (1979) studied the stability of natural slopes on Prince of Wales Island (Alaska), before and after the removal of a forest cover. They incorporated the apparent cohesion owing to the roots into simple limit equilibrium analyses, assuming infinite slope and steady-state seepage conditions. The results indicated that the additional strength provided by the tree roots is important for the stability of steep slopes and the loss of root strength following clear-cutting can seriously affect slope stability. Chok et al. (2004) and Gentile et al. (2010) analysed the mechanical effect of vegetation on the stability of ideal slopes, characterized by dimensions and material properties typical of highway embankments, using the FE method. The results show that vegetation plays an important role in stabilizing shallow-seated failure of slopes, indicating that the increment of the slope factor of safety (FoS) is more significant if the slope toe elements are treated as vegetated soil. The effect increases with the deepening of the root system, reaching the zones where the failure mechanism starts. Nevertheless, the factor of safety of the slope cannot increase indefinitely, as it reaches an asymptotic value with increasing root additional cohesion. Similar results were obtained by Ji et al. (2012), who performed 2D FE simulations to investigate the effect of root additional cohesion on the stability of two natural slopes in NW China.

Another example of numerical assessment of the root mechanical reinforcement is given by the work by Koda & Osinski (2011), who studied the stability of the Radiowo landfill located in the northwestern part of Warsaw (Poland). The stabilizing effects of the vegetation have been considered in the factor of safety calculation, making use of the General Greenwood Method (Greenwood 2006). The results of the stability analysis indicate that a 20% increase of the factor of safety can be reached along the vegetated slopes of the landfill, thus showing how slope reinforcement does not always require heavy engineering methods, but can be achieved with simple, environment-friendly and cost-effective techniques, such as the implementation and maintenance of a vegetation cover. Nevertheless, the mechanical improvement provided by the root systems is strongly controlled by the depth of the root zone. Therefore, the vegetation mechanical effects are less significant in slopes where deep-seated failure mechanisms are likely to occur.

This section presents a collection of real and ideal cases where the stability of natural and artificial slopes subjected to climatic perturbations has been investigated by some of the participants in the COST Action, using different modelling strategies, classified as outlined so far (Table 1).

The influence of the slope–atmosphere interaction on slope stability has been investigated implementing the results of hydraulic slope modelling in limit equilibrium analyses for the Pisciolo hill-slope (Melfi, Italy) (Cotecchia et al. 2014), which is the location of deep slow landslide movements typically observed in fissured clay slopes in the south of Italy. The slope has been thoroughly studied because important infrastructure has been severely damaged owing to its interaction with the slope movements. In particular, a pipeline, a road and a railway are located at the toe of the slope, where landslide activity is highest. The geological setting of the area has been carefully reconstructed based upon in situ surveys and underground investigations. The slope is mainly formed of clayey turbidites that are part of the Paola Doce formation and consist essentially of high-plasticity fissured clays with interbedded coarse and fractured rock inclusions. Paola Doce clays do not outcrop all over the slope, due to the presence of an anticline with an axis that crosses the hill-slope transversally. The central portion of the slope has been uplifted and deeper clays, belonging to the Red Flysch formation and characterized by even higher plasticity and fissuring degree than the Paola Doce clays, crop out about the fold axis. Also, more recent arenaceous and sandy materials, which are part of the Numidian Flysch, crop out on the fold flanks. This geostructural reconstruction has supported the generation of the geotechnical model of the slope, which implements much higher strength properties and higher permeabilities where the Numidian Flysch and the rock or sand inclusions are located, and a great reduction in strength parameters and permeability within the slope locations of the Paola Doce clays (c′_av= 0 kPa, ϕ′_av = 20.7°) and the Red Flysch clays.

Ten rotational–translational clay-slides have been identified to be active within the Pisciolo hill-slope, based upon stereoscopic analyses of past and recent aerial photographs, geomorphological surveys and the interpretation of inclinometer data. The landslides involve mainly the Paola Doce soils (clays and interbedded sands and rocks) and, locally, the Red Flysch. The slope has been subjected to extensive piezometric monitoring, conducted by means of Casagrande and electrical piezometers. Most piezometric levels reach a few metres below the ground level, also for very deep piezometers (down to 60 m depth). The landslides are a consequence of the relatively high slope mass permeability, arising from the fissuring of the clays and the occurrence of the coarse strata and the fractured rock inclusions within the clays. As well as the presence of the higher permeability inclusions, the average permeability of the clay portions (site measurements: k_sat > 10⁻⁹ m s⁻¹), which is higher than for unfissured clays by a few orders of magnitude, eases rainfall infiltration, which generates high piezometric heads that are detrimental to the available shear strengths. Weather data have been collected at a weather station near the slope. The 180 day cumulative rainfall evolves according to a seasonal trend, with maximum values at the end of winter and minimum values at the end of summer. The piezometric heads and the displacement rates of the most active landslides, located in the southern portion of the hill-slope, follow a similar fluctuation trend. In particular, at the toe of the most active rotational–translational multiple landslide, the displacement rates have been measured using inclinometers and global positioning system (GPS) sensors and have been found to fluctuate seasonally. The correspondence between the pore pressure excursions, the displacement rate variations and the 180 day cumulative rainfall variations confirms that the climatic regime represents the main triggering factor of the current activity of the Pisciolo landslides, and the relatively high mass permeability represents the internal factor that mainly predisposes the slope to landsliding, along with the relatively low available strengths of both the Paola Doce and the Red Flysch clays.

Two-dimensional transient seepage analyses have been performed with the finite-element code SEEP/W (Geo-Slope International 2004), which integrates the Richards’ equation, according to the method presented above. The analysed section, which crosses longitudinally the most active landslides, is reported in Figure 2, together with plots of the model results compared with some monitoring data. The section includes the different materials found in the slope: the Paola Doce clays, the Red Flysch clays and the more permeable fractured rock inclusions. The van Genuchten (1980) model has been used to fit the laboratory retention data for the fissured clays and to simulate the permeability variation with suction. Pore pressures constantly equal to zero have been imposed at the ground surface at the upper and lower portions of the slope ground surface, according to the corresponding presence of a spring and of the Ofanto river. Net rainfalls referring to the period September 2006–August 2007 have been cyclically applied along the rest of the slope ground surface. The net rainfalls have been computed as the difference between the total rainfalls and the evapotranspiration fluxes, the latter estimated by means of the FAO Penman–Monteith method (considering the presence of winter wheat and using the single crop coefficient approach). The numerical results have been compared with the piezometric heads measured along three verticals. As shown in Figure 2a for the vertical P7, the piezometric variations predicted by the numerical model are in good agreement with those measured in situ, especially once the more permeable inclusions are implemented in the model. When the pore pressure distributions resulting from the seepage analysis for different stages of the year are implemented in the limit equilibrium analysis, variations of the factor of safety from 5 to 20% are calculated for the landslide bodies in Figure 2b. Such safety factor variations reflect the changes in actions that cause the seasonal accelerations of landsliding on the Pisciolo slope, which results in the recurrence of damage to both the pipeline and the road present at the toe of the slope.

The failure of a peat dyke at Wilnis (Netherlands) during summer 2003 provides evidence of the effects of severe droughts on the stability of this type of engineered slope. The failure did not cause casualties, but its economic impact was important, as the resulting economic loss was about €20 M. Since then, interest in the evaluation of the safety of the 7000 km of peat dykes existing in the Netherlands has increased, as demonstrated by the ‘Knowledge for Climate’ programme (2008 – 2014) of the Dutch Ministry of Infrastructure and Environment, which has supported research into the behaviour of engineered slopes under changing climatic conditions. Within this research programme, the stability of an existing dyke near Boskoop (Netherlands), the Middelburgse kade peat dyke, has been studied by Van Esch (2012). This example is presented here with the aim of outlining the evolution of both seepage and stability applying to peat dykes and to assess the impact of climate changes on these engineering structures, for more efficient and effective advice for their maintenance and remediation. The dyke is characterized by a height of about 3 m and a longitudinal length of 40 m, with height-to-width ratio of 1:14. Figure 3a shows a schematic cross-section of the dyke and the peat layers that cover a sandy subsurface.

Previous research by Van Esch et al. (2007) has proposed a procedure to assess the consequences of droughts for peat dykes, using the finite-element software Plaxflow (Brinkgreve et al. 2003). This procedure complemented an agro-meteorological model to derive the evapotranspiration flow at the ground surface, based on the Penman–Monteith expression, and a groundwater flow model to predict the piezometric distribution based on the Richards’ equation (hydraulic modelling). Flow through cracks in unsaturated subsoil was simulated via a double-porosity approach and Scott's model (Scott et al. 1983) was adopted to simulate the drying and wetting behaviour of peat, as laboratory experiments could not be reproduced well by a single van Genuchten expression. It was found that the approach adopted to assess the slope piezometric and stability conditions was computationally very demanding; this approach is worth using if detailed information on material behaviour is available. Therefore, using Richards’ model it was found that is very important to account for anisotropy in the permeability tensor, because of the horizontal layering of peat. The horizontal component of the permeability tensor in the mostly saturated zone is about 10 times its vertical component. Hence, a Dupuit approximated approach (De Marsily 1986) has been used to simulate the flow through the saturated zone of peat dykes. The Dupuit approximation states that groundwater flow in aquifers of high permeability, such as sandy deposits, is oriented horizontally and the flow through the aquitards interbedding low-permeability clay layers takes place in the vertical direction. To verify the extent to which the simplified Dupuit approach is able to provide realistic predictions in practice, Van Esch (2012) employed the Penman–Monteith expression to generate boundary conditions for a Dupuit model, and used a finite-volume scheme, implemented using in-house software (Van Esch et al. 2013), to simulate transient groundwater flow. In the case under investigation, saturated flow through the peat dyke is oriented horizontally and unsaturated flow takes place in the vertical direction. The proposed model incorporates effective porosity and effective permeability in the unsaturated zone. Effective parameters have been obtained from a calibration procedure and have been computed more easily than model parameters that support the soil water retention curves. The simplified model is computationally highly efficient. Figure 3a shows the computed piezo-lines in the aquifer and the dyke, from which the water pressure distribution can be derived. The figure also shows the location of four observation wells, where groundwater heads have been monitored over a period of 2 years. The numerical predictions in terms of hydraulic head obtained at different observation points with the finite-volume Dupuit model have been compared with those obtained through the more advanced finite-element simulation using Richards’ model. Both groundwater flow models produce results that compare well, as shown in Figure 3b, thus demonstrating the good ability of the Dupuit approximation to simulate groundwater pressure fields. Figure 3c compares the calculated heads at different observation points with observed heads in the field, showing a good agreement between numerical predictions and measured data. The pore pressure distribution calculated with the Dupuit model is then used as input for a limit equilibrium analysis, which employs Spencer's method (Spencer 1967) to evaluate the evolution of the dyke stability with time. The critical slip surface is shown in Figure 3a, and Figure 3d presents the stability response of the dyke over the interval under observation. The figure indicates a stable condition throughout the simulation period as the stability factor is well above 1.0. In general, the stability of the dyke is found to decrease when its weight reduces. For peat dykes the weight can strongly decrease depending on the amount of water that leaves the system by evaporation. A tipping point analysis has, therefore, revealed that an increase of evapotranspiration by a factor of two would have led to the failure of the peat dyke (Fig. 3d).

The Po river (northern Italy) is the Italian main water course in terms of length and capacity. Past flooding events have often caused damage and casualities. An ad hoc public body, the Po River Management Authority ‘AIPO’, is in charge of the river safety issues and early warning system, manages emergencies, and finances and coordinates research activities addressing risk mitigation. In the past, flooding prevention required remediation of existing defence embankments, through height increase (even up to 10 m), cross-section enlargement, or improvement of the mechanical properties of the construction materials. Throughout the river network (about 2000 km), embankments are made of compacted soils taken from pits along the river banks. To investigate the hydraulic behaviour of both the fine-grained embankments located in the mid-course of the river and the coarse-grained ones along its upper course, research studies have been carried out on full-scale physical models at three sites: in Viadana (MN), for the study of the mid-course embankment's behaviour, and in Bormida (AL) and Motta dei Conti (VC) for the upper-course embankments. The main objective of the research study has been to define a design method that takes account of the seepage in the banks and in the foundation soils, under either partially or fully saturated conditions (Calabresi et al. 2013; Belardi et al. 2014).

In the case of Viadana, according to recommendations formulated by the Po River Management Authority, an embankment prototype was built on the floodplain, beside the existing embankment, forming a pond. The pond was filled to reproduce historical floods. The engineering properties of the foundation soils have been investigated through in situ tests and complemented by suction-controlled laboratory tests (Calabresi et al. 2013). The soil profile at Viadana includes different layers, consisting of medium to fine sands, clayey silts, medium to fine silty sands and sandy gravels. Porewater pressures have been measured in the embankment and in its foundation before, during and after the experimental simulation of two floods, performed to reproduce those that occurred in 1976 and in 2000. Atmospheric variables have been monitored at the site during the 6 months of experimental activity. The transient seepage processes occurring in unsaturated conditions within the prototype, given the ponding conditions and the slope–atmosphere interaction, have been investigated by Calabresi et al. (2013) through hydraulic analyses conducted with the finite-element code SEEP/W (Geo-Slope International 2004). The mesh adopted in the analyses is reported in Figure 4a. The water flows at the ground surface due to slope–atmosphere interaction (i.e. rainfall and evapotranspiration) have been computed separately and input into the model, according to an isothermal hydraulic modelling approach (see Table 1), integrating the Richards’ equation and neglecting both hydromechanical and hydrothermal coupling. The transient seepage stages following changes in ponding conditions have been analysed. With this aim the hydraulic boundary conditions have been changed over time, according to steps corresponding to the history of the impounding levels of the Po river. A ‘seepage surface’ condition has been applied along the unsaturated portion of the embankment boundary not subjected to impounding, where pore pressure proved to become either negative or zero. Along this boundary, soil–atmosphere interaction has been simulated. Setting the flow rate as the net balance between the outflow owing to evaporation and the inflow owing to rainfall, the prediction of the soil suction has been less satisfactory than that of the subsoil porewater pressures (see Fig. 4b). The transient stages have been found to be persistent and a significant delay in reaching steady-state conditions occurs, due to the low unsaturated soil permeability, consistent with the soil suctions varying according to the retention curve of the material. The water retention properties of the embankment materials are in reality inhomogeneous, given the randomly distributed variations in grain size and/or in compaction degree of the soils. As a consequence, suctions may have been more scattered than those predicted by the analyses performed assuming homogeneous embankment properties. With regard to the differences in the predictions when accounting for the slope–atmosphere interaction effects and when not, the comparison in Figure 4c shows that these differences are major for tensiometers close to the embankment boundary (see locations T1A, T2A and T6A). Within the embankment (see locations T3A and T5A), the effects of the slope–atmosphere interaction decrease. Hence, it emerges that modelling soil–atmosphere interaction is crucial to predict a reliable evolution and distribution of the suctions across embankments.

The distribution of the porewater pressures provided by the modelling has been input into embankment stability analyses, using Bishop's traditional limit equilibrium method. The safety factors computed in transient conditions are increasingly higher than those associated with steady-state conditions reached when impounding is maintained for a time sufficient to attain a new hydrostatic regime. There is little penetration, in the short term, of the saturation zone during severe flooding, as the vanishing of suction is limited to a thin zone near the wet boundary and, in some cases, along the foundation boundary. Therefore, given the presence of negative pore pressures in large portions of the embankment, realistic stability analyses of the embankment during flooding should be carried out by taking into account the strength contribution provided by suction.

The stress–strain conditions in a slope location of a deep slow-moving landslide in fissured marly clays, the Fontana Monte slope at Volturino (Southern Italy), have been investigated by Lollino et al. (2010, 2016) through hydromechanical finite-element analyses performed with the code PLAXIS 2D (2012) according to the method presented above. The landslide is observed to undergo seasonal reactivations, as the rate of movements increases in late winter–early spring and reduces in summer, as observed for several slopes formed of clay flysch in the Italian Southern Apennines. The landslide body has been found to involve mainly the Toppo Capuana Flysch, which is mostly made of marly clays with an average clay fraction of about 50% and plasticity index increasing with depth, from about 10% at shallow depth, to 40% at 45 m depth. Geomorphological analyses (Fig. 5a) had revealed that the Fontana Monte landslide was already active in late nineteenth century, so that the current slope activity can be considered to be largely the effect of further straining, relating to climate, along pre-existing shear bands. In particular, the seasonal accelerations of the landslide suggest that periodic variations of the shear strength occur in the slope, probably as a result of seasonal porewater pressure variations. Seasonal fluctuations of the piezometric levels, resulting from transient seepage processes generated by the slope–atmosphere interaction, have been monitored in the slope to depth, as reported by Lollino et al. (2016). The landslide activity in the upper portion of the slope is documented by cracks and fissures that affect buildings and roads located along and at the back of the rear scarp of the landslide, at the top of the slope, as well as along the lateral borders of the landslide body (see Fig. 5a), as consequence of differential settlements of foundations connected to the landslide movements. The analyses have explored the slope straining processes evolving from summer to winter in the slope, assuming the slope to be formed of undisturbed clays that include a deep pre-existing shear band, whose morphology has been derived from investigations of bored corings and from inclinometric monitoring. The variations in piezometric head from summer to winter have been input into the analyses through the simulation of both the summer and the winter steady-state seepage. The validity of the seepage simulations has been checked against piezometric site measurements at several locations along the slope. A standard linear elastic–perfectly plastic Mohr–Coulomb model with non-associated flow rule has been adopted for the materials, in accordance with the moderately brittle soil shear behaviour observed in the laboratory. A cohesion of 13 kPa, an effective friction angle equal to 20° and a coefficient of permeability equal to 5 × 10⁻⁹ m s⁻¹ have been adopted in the analyses. The analyses have been carried out assuming the shear band to be of either 37 m maximum depth or of 50 m maximum depth, to evaluate the influence of the depth of the pre-existing shear band on the current slope stability. Because the shear band is formed of pre-failed clays, strength parameters intermediate between peak and post-peak have been used for the shear band (i.e. c′ =8 kPa and ϕ′ = 18.7°). The summer and winter analyses have been carried out under plane-strain drained conditions, with pore pressure distributions resulting from steady-state seepage calculations fitting the site measurements, as explained above in the section on hydraulic modelling. The maximum localization of shear straining has been found to occur within the 50 m deep shear band when the piezometric heads increase from the summer values to the winter ones. Figure 5b and 5c report the contours of incremental shear strain in the slope, which provide indication of the shear strain concentrations in the deepest shear band. The slope is found to be stable in summer and to fail in winter (i.e. loss of numerical convergence). The corresponding contours of vertical incremental displacement are reported in Figure 5d: negative displacements (i.e. settlements) are recorded in the upper part of the slope and in the centre, whereas maximum positive displacements (i.e. heave), of about 20 mm just before loss of convergence, are observed at the base, where the toe of the landslide body is located. Maximum incremental settlements, of about 10 mm, are calculated in the scarp area, where buildings are periodically damaged. As such, the accumulation of settlement along the slope crest over few years provides a reason for the observed pattern of structural damage.

This numerical study proves the tendency of the deep Fontana Monte landslide to reactivate as an effect of seasonal piezometric fluctuations. Furthermore, it proves that this type of slope is most unstable when very deep pre-existing shear bands are present, with respect to shallower ancient shear bands. However, it is worth noting that a significant improvement in the displacement predictions would result from the adoption of a strain-softening model, especially for the undisturbed clays, outside the ancient shear band. Also, the use of 3D numerical modelling, implementing hydromechanical coupling within the simulation of the transient seepage, from summer to winter, would provide a great enhancement in the prediction of the slope response to seasonal rainfall infiltration and evapotranspiration during the year (e.g. Pedone et al. 2016).

The stability of a clay slope in the highly urbanized area of Veles (Macedonia), the Ramina slope along the left bank of the river Vardar, has been studied by Josifovski et al. (2013). The slope is the location of an active landslide since the nineteenth century, which has been reactivated several times, most recently in 1963, 1999 and 2002. The sliding mechanism is composed of two sliding masses, the upper one of length 350 m, width 110 m, height 60 m and depth up to 24 m, and a lower one of length 200 m, width 90 m, height 35 m and depth up to 18 m (Fig. 6a). The separation of the two masses is due to the presence of a rock buttress in the middle of the landslide (see Fig. 6a). The total area covered by the landslide is c. 37 600 m² and involves a mass of around 475 200 m³, making it the largest landslide in the Balkan countries and, possibly, in SE Europe. The geotechnical investigations indicate that the sliding zone is along the contact between low-plasticity clay-like materials and the bedrock, which is represented by Palaeozoic amphibolitic schists that are highly foliated, faulted and locally weathered. In the upper part of the slope, additional fill is present on the top of the clay layer and the water table occurs at 8 m depth, whereas in the lower part the water table is 2 m below the ground surface. The hydrogeological conditions indicate that the zones of higher water content correspond to the sliding zone. Also, the piezometric levels at depth are found to be higher than those at shallower depths in almost all boreholes (Jovanovski et al. 2003). The material properties have been determined through extensive laboratory and field investigations, performed from 1963 to 2002. According to site observations, the sliding process appears to be caused by heavy rainfall, as the rainfall water infiltration accesses the weak schistose rocks. The climate is semi-arid, characterized by gusty rainfall events from October to March, with an average monthly rainfall of around 100 mm. After the landslide occurrence in 2002, hydromechanical simulations have been conducted (according to the method presented above in the section on hydromechanical modelling) to evaluate the effects of slope–atmosphere interaction on the stability of the slope and estimate the potential further sliding, using the finite-element code PLAXIS 2D (2012). Based on the available data about the mechanical behaviour of the slope materials, the elasto-plastic Mohr–Coulomb model has been adopted. The van Genuchten hydraulic model has been employed in the analysis and calibrated against measured data, as presented in Figure 6b. The rainfall infiltration has been set as 10 mm h⁻¹ at ground surface, plus an additional 1 m³ h⁻¹ on the right side as inflow boundary. The runoff has been assumed to be not greater than 0.5 mm h⁻¹ owing to the surface characteristics at the site (which is fairly urbanized). In the model vegetation and temperature have not been taken into account and evapotranspiration has not been modelled. Mesh updating has been performed, due to the occurrence of large straining. A realistic simulation of the sliding mass behaviour has been obtained, as presented in Figure 6c: the sliding of the mass occurs from around 12 h after the start of the rainfall, with maximum suction of 169.7 kN m⁻² (Fig. 6d) and maximum displacement of 1.17 m (at the rear scarp). The results show a sliding mechanism whereby the upper part of the slope starts to move as a result of the rainfall event, which causes a significant rise of the water table and the porewater pressures. Thus, different ground water levels develop in the upper and lower part of the slope, representing a trigger for the sliding of the upper mass and the loading at the top of the lower mass. The slope deformations, recorded by means of inclinometers located along sections of the slope different from the analysed one (about 0.5 m displacement in the middle of the landslide), have been qualitatively compared with those simulated through the hydromechanical analysis (reported in Fig. 6c). The comparison appears to indicate that the hydromechanical simulation of the failure mechanism so far presented is reasonable.

An instrumented cutting on the A34 Newbury bypass in southern England has been studied by Rouainia et al. (2009) to investigate the effects of present and future climate on infrastructure slope stability. The slope is formed of London Clay overlying Lambeth Group deposits and the Upper Chalk (see Smethurst et al. (2006) for further details) and takes the form of an 8 m high, 28 m wide cut (Fig. 7a). In the modelling, the cutting surface is vegetated with grass, but in reality, over time larger shrubs have started to grow on the site (Smethurst et al. 2012).

The modelling has been undertaken by inputting currently monitored meteorological variables into a SHETRAN (Ewen et al. 2000) hydrological model of the site. SHETRAN uses weather data to estimate the water flow in the soil and the evapotranspiration from the slope surface. The hydrological model assumes that surface soil deformation is negligible and does not significantly affect the porewater pressures. In addition, future climate data have been created using a weather generator based on the UK Climate Impacts Programme (UKCIP02) projections (Kilsby et al. 2007) and applied to the same SHETRAN model to estimate porewater pressure response to future climate. Once the daily surface porewater pressures have been calculated with SHETRAN, they have been transferred to the surface of a finite-difference geotechnical model developed in FLAC (Itasca 2011) to simulate the coupled hydromechanical slope response according to the method presented in the modelling section above. In this way, the effects of present and future climate on the mechanical response of the cutting have been modelled. The van Genuchten model has been used to derive the soil water retention curve and the hydraulic conductivity function (as summarized in Fig. 7a). Through this modelling, changes in the slope porewater pressure have been calculated as response to the weather conditions with time. The modelled porewater pressure response has been calibrated to comply with published monitored porewater pressure data from Smethurst et al. (2006) through a parametric study on the effect of the assumed coefficient of permeability. It has been observed that the saturated coefficient of permeability required by the hydromechanical model to reproduce correctly the measured pore pressures is three orders of magnitude greater than the coefficient of permeability recorded within triaxial cell testing and two orders of magnitude greater than that estimated on site from bail-out tests. Field permeability has been therefore selected, as reported in Figure 7a, to allow for the macroscopic effects of fissures, sand lenses and other heterogeneous features, as well as to compensate for scale factors. The hydromechanical coupling has been accomplished using a mean effective stress-dependent stiffness, strain-softening Mohr–Coulomb model. The predicted mid-slope horizontal displacements and the porewater pressures at the ground surface for both the present and future climate conditions are plotted in Figure 7b. A number of readily apparent differences between the two responses are evident. The first concerns the magnitude of the suctions calculated for future climate, compared with the minor present-day suctions, where it can be seen that the peak summer suctions generated by the future climate are always larger than those for the present: the future minimum suction of 575 kPa exceeds the maximum suction for the present climate. The maximum future climate suctions in turn reach c. 1175 kPa in the summer of year 9 and approach or exceed 1000 kPa in eight of the modelled years. This increased variation in suction appears to be due to the combined effect of the more intense drying and the larger evapotranspiration that will be connected to the higher temperatures expected in future. In both cases, it can be seen that the suctions totally dissipate during winter and the corresponding larger effective stress fluctuations are reflected in larger shrink–swell cycles and permanent seasonal downslope movements. The larger deformations for the future climate lead to more rapid softening of the material at the toe of the slope and, ultimately, to the failure of the slope in only 22 years, whereas for the present-day climate failure occurs after 61 years. It can be concluded that the model has the potential to assist in the assessment of the future behaviour of clay cuttings and embankments subjected to modified climatic conditions.

Within the context of the multidisciplinary research programme ‘Triggering of Rapid Mass Movements in Steep Terrain’ (TRAMM, http://www.cces.ethz.ch/projects/hazri/tramm), a full-scale field test has been set up on a steep forested slope located on the east-facing banks of the river Rhine near the village of Ruedlingen, in the Schaffhausen canton, northern Switzerland. Here, several shallow landslides occur due to extreme rainfall. To analyse the effects of hydrological events on the stability of unsaturated silty sand slopes, a 38° steep forested slope has been monitored during the application of artificial intense rainfall events, after carrying out an extensive geo-hydro-mechanical characterization of the site. The selected experimental site, of 35 m length and 7.5 m width, is a small part of a slope on the east-facing bank of the river. Fifteen hours of artificial rainfall, with an average intensity of 11 mm h⁻¹, has been applied to the slope, and this resulted in the triggering of a landslide of 130 m³ volume. The lithostratigraphic profile of the slope consists mainly of Seawater Molasse from the Tethys Ocean and Freshwater Molasse. Several boreholes, as well as natural scarps, revealed the presence of horizontal layering of sedimentary rocks (fine-grained sandstones and marlstones; Tacher & Locher 2008). Fissures, with openings of more than some centimetres in size, have been mapped in the lower Freshwater Molasse and found to be parallel to the Rhine (Brönnimann 2011). The bedrock depth is between 0.5 and 5 m, according to dynamic probing (DPL) and electrical resistivity tomography (Lehmann et al. 2013). The bedrock on the right side of the field (looking upslope) is found to be shallower than on the left side. The soil can be classified as medium- to low-plasticity silty sand (ML), according to the Unified Soil Classification System (USCS). The clay fraction increases with depth, from 4% at shallow depths to 10% at about 2 m depth. The silt fraction also increases with depth, from 25 to 32%, whereas the sand fraction decreases from 67 to 56% (Casini et al. 2010). The results of conventional drained and undrained triaxial tests on saturated samples revealed a critical-state friction angle of ϕ′_CS = 32.5° (Casini et al. 2010, 2013). Moreover, a series of unsaturated constant shear drained triaxial tests (as proposed by Brand 1981) have been performed by Askarinejad (2013) to replicate the stress path that a soil element experiences during the rise of porewater pressure in the slope owing to rainfall infiltration. The results suggest that the stress ratio attained at failure (ϕ′ = 34.2°) is generally higher than the critical-state stress ratio determined by conventional drained or undrained triaxial tests, consistent with what was observed by Anderson & Sitar (1995), Zhu & Anderson (1998) and Casini et al. (2010). Several direct shear tests have also been performed on statically compacted unsaturated specimens and the results are indicative of a peak internal friction angle of 30° (Colombo 2009; Askarinejad et al. 2014). Wetting and drying branches of the SWRC of natural samples have been determined using the axis translation technique (Hilf 1956). The air-entry value of the drying branch is 1.7 kPa (Askarinejad & Springman 2015). In all the numerical analyses presented in this section, the unimodal van Genuchten formula has been used to fit the data for both wetting and drying branches.

Two-dimensional uncoupled hydraulic finite-element simulations have been carried out by Askarinejad et al. (2012) and Askarinejad (2013) to simulate the conditions of the pilot site slope before the artificial rainfall event. The mechanical features of the unsaturated soils and the reinforcing effects of vegetation have been implemented. The simulations have been performed using the SEEP/W and SLOPE/W modules of GeoStudio (Geo-Slope International 2004), according to the method presented above. The hydraulic boundary conditions imposed during the test at the ground surface have been applied to the finite-element (SEEP/W) model and the pore pressure distribution with time has been calculated. The factor of safety of the slope has been determined using SLOPE/W, based on the pore pressure distribution exported from SEEP/W for each time step. The pore pressure distribution has been calculated according to Richards’ equation for flow in unsaturated soil. The shear strength of the unsaturated soils has been simulated in SLOPE/W using the extended form of the Mohr–Coulomb criterion after Fredlund et al. (1978). Different bedrock topographies (derived from the DPL results) and inclinations of fissures have been examined. The bedrock has been considered to be elastic and impermeable, compared with the soil layer. The location and sizes of the fissures in the bedrock have been determined by monitoring the sequences of the changes in the degree of saturation in the bedrock. The fissures have been implemented in the geometry of the model and have been filled with the overlying soil, with similar hydromechanical properties. Figure 8a shows the geometry of the model for the right-hand sides of the test site (looking upslope) in the case of the shallower bedrock. Fifteen hours of rainfall, with constant but different intensities for seven zones, have been applied all over the slope, in accordance with the real rain event in the Landslide Triggering Experiment (LTE, March 2009). The seepage analyses show that the water applied in rain zone 7 is conveyed through a fissure in the upper part of the slope to the horizontal permeable layer and then it exfiltrates from the bedrock back into the soil layer, applying an uplifting hydraulic pressure (Fig. 8b). However, the majority of the infiltrated water drains through the fissures into the bedrock in the lower part of the slope. The evolution of the factor of safety as a function of time during the LTE indicates that the FoS decreases with the rain and attains a minimum value 2.6 h after the rain stops. The results of the limit equilibrium analyses suggest that the reduction in factor of safety occurs more rapidly in the case of the shallower bedrock as a result of the exfiltration process described in Figure 8b.

In addition, a series of 2D and 3D coupled hydromechanical finite-element analyses have been performed by Askarinejad (2013) to take the coupled response of the unsaturated soil into account, according to the method presented above. Only the results of the 2D coupled analyses are commented on here. The simulations have been performed using the finite-element software CODE_BRIGHT (Olivella et al. 1996). The constitutive modelling of the coupled hydromechanical response of the unsaturated soils has been achieved using the Barcelona Basic Model (BBM) (Alonso et al. 1990). The time dependence of the stresses and strains of the material, owing to viscosity, has been also taken into account for the plastic region, based on the formulation suggested by Perzyna (1966). Accordingly, the model used for these series of simulations is an elastic–viscoplastic model. The constitutive model parameters have been defined by comparing the yield function of the BBM with the modified yield function proposed by Desai & Zhang (1987) for unsaturated soils, as specified by Askarinejad (2013). It should be noted that the time dependence of the downslope movements is more related to the cyclic changes in suction and volumetric water content induced by the infiltrated rainfall than the viscosity of the material, which has been introduced in these simulations for numerical convenience. The equation proposed by Averjanov (1950) has been used to describe the unsaturated hydraulic conductivity function in CODE_BRIGHT, whereas the saturated hydraulic conductivity has been determined to be 5 × 10⁻⁵ m s⁻¹ (Askarinejad et al. 2010). Two analyses have been carried out to investigate the effect of exfiltration from the bedrock on the response of the overlying soil layer during the LTE. No exfiltration has been simulated in the first model (2D_D), whereas water exfiltrating from the bedrock at two locations has been considered in the second model (2D_DE). The geometry of the 2D_DE model, the location of the fissures in the bedrock and the hydraulic boundary conditions are shown in Figure 9a. The maximum dimension of the triangular elements used in the FE mesh is 0.3 m. The bedrock is assumed to be impermeable, with draining fissures that are located beneath clusters 1 and 2. The lower boundaries of the model have been set as ‘seepage face’ lines, where the water can drain out from the model and the pressure head is kept to zero. Exfiltration has been modelled at a simulation time of 5 h, by applying 9 kPa constant-head boundary conditions at the desired locations in the upper slope (cluster 3). The contours of total displacement before the failure for the 2D_DE model are shown in Figure 9b. The approximate shape of the failure surface of the Ruedlingen slope from the landslide triggering experiment is also depicted in this figure. Changes in porewater pressure at depth of 150 cm in cluster 3 are shown in Figure 10. The differences between the models with and without exfiltration are significant. The results of the 2D_DE model have been found to be in fairly good agreement with the field measurements in terms of the evolution of the porewater pressures (Fig. 10) and location of failure (Fig. 9b). However, the time of the failure did not match the in situ observations (i.e. failure occurred c. 13 h after the start of the rainfall). Nonetheless, the numerical analysis confirmed the importance of the hydraulic interaction of the bedrock with the overlying soil layers, in terms of drainage and exfiltration, and showed that these interactions can play a major role in the pattern of pore pressure distributions and in stabilizing or destabilizing the slope.

A theoretical slope, characterized by the typical geometry of a cutting or embankment for linear infrastructure, has been modelled by Arnold & Hicks (2011) to analyse the influence of soil heterogeneity on the slope stability during a rainfall event. The so-called random finite-element method (RFEM) has been employed (Griffiths & Fenton 1993), where a series of different statistically possible set-ups (realizations) of material properties have been simulated. Using a Monte Carlo framework, the data have been compiled into a statistical output of the factor of safety, leading to a description of the reliability of the slope. The model uses Bishop's effective stress and involves a saturated or unsaturated seepage analysis performed in a stepwise staggered manner (as described above); that is, a hydraulic analysis is undertaken at a single time step, with the results imported into a stability analysis. The strength reduction method (Griffiths & Lane 1999) has been used to determine the FoS at each time step. Partial saturation is included with both the SWRC and the unsaturated hydraulic conductivity following the van Genuchten formulation. A linear elastic–perfectly plastic soil model has been adopted in the simulations. Hydraulic conductivity and shear strength parameters can be correlated or independent and it is possible to incorporate anisotropy. The model has been designed to investigate the effect of soil heterogeneity on stability, and therefore embodies a number of more simplified features for other aspects; for example, no hardening or softening is incorporated in the constitutive model, fluid flow is single phase with air pressures assumed to be atmospheric, and vegetation and other surface effects are neglected, as are any thermal influences. The model has been implemented in an in-house finite-element code, coupled to a Random Field (RF) generator (Hicks & Samy 2002). The RF generator produces the different set-ups of material properties by use of the Local Average Subdivision (LAS) method (Fenton & Vanmarcke 1990), which utilizes as input the mean, standard deviation and scale of fluctuation (a measure of spatial correlation between values at different locations). This method generates fields with the same scale of fluctuation in both vertical and horizontal directions, which can be post-processed to produce an anisotropic field by squashing or stretching the isotropic field. Rainfall events in this model are simulated by fixing a zero pore pressure value on the nodes at the embankment surface. This approach is valid where the rainfall event is sufficient to fully saturate the upper layers of the slope. Any additional rain is assumed to run off, and is not included in the simulation. Care must be taken not to simulate excessive inflow, and it is possible to include mixed boundary conditions to ensure this automatically.

In the study by Arnold & Hicks (2011), a 45° slope of 5 m height resting on a 5 m basement layer has been analysed in two dimensions, assuming plane strain conditions and using 0.25 m wide eight-noded quadrilateral finite elements. Five soil properties have been selected as spatially varying in the slope: the effective cohesion and friction angle, the porosity, the saturated hydraulic conductivity and the inverse of the air-entry pressure. Their mean values have been set equal to 10 kPa, 25°, 0.4, 0.0036 m h⁻¹ and −1.0 m⁻¹, respectively, with coefficients of variation equal to 0.2, 0.3, 0.15, 1.75 and 0.9. The five properties have been assumed to be log-normally distributed. Based on test results found in the literature, a typical correlation matrix has been set up to define the covariance structure for this boundary value problem. In particular, positive correlations have been assumed between porosity, saturated hydraulic conductivity and the inverse of the air-entry pressure, whereas negative correlations have been assumed between friction angle, saturated hydraulic conductivity and the inverse of the air-entry pressure. In the case presented by Arnold & Hicks, the effective cohesion and friction angle have been considered constant, with only the SWRC parameters varying in space. A homogeneous case, for which the mean values of the above-mentioned soil properties have been adopted, and a series of different heterogeneous slopes have been investigated. In the heterogeneous case, a vertical scale of fluctuation of 2 m and a horizontal scale of fluctuation of 6 m have been assumed. The stability of the models has been investigated applying a continuous flux of 18.0 mm h⁻¹ for 48 h at the ground surface, representative of a heavy rainfall event, followed by 1.0 mm h⁻¹ for about 550 h.

For a given factor of safety, the reliability of the slope is equal to R = 1 − P_f = 1 − N_f/N, where P_f is the probability of failure, N is the total number of realizations and N_f is the number of realizations in which the slope fails. A sufficient number of realizations must be considered to produce a converged result, which increases with the higher reliability levels. In this case, 300 realizations have been found to provide a satisfactory response. The FoS calculated prior to rainfall, without taking into account suction, is 1.39, which increases to 1.75 taking into account suction stresses. The FoS is observed to fall during and after rainfall as the water front propagates through the soil. In the homogeneous case the FoS falls from 1.75 to 1.58 after c. 500 h. Figure 11 presents a comparison of the Bishop's stress contours for the homogeneous and a single realization heterogeneous slope, with the corresponding potential failure mechanisms shown by the deformed mesh. The results are presented after 48 h, immediately after the rainfall event. A reduction in Bishop's stress can be observed close to the surface in both cases owing to increased levels of saturation (and reduction in suction). In the heterogeneous case, the FoS is reduced from 1.90 to 1.42 over the same time period. The failure mode is also seen to be different, with a shallower failure surface noted for the heterogeneous case. It is also observed that the mean FoS of the heterogeneous slope is lower than that of the homogeneous case; that is, a reliability of 50% is c. 90% of the homogeneous FoS and a 99.75% reliability gives an FoS that is c. 70% of the homogeneous case, thus indicating how failure can be driven by the spatial variability of the effective shear strength parameters within the slope.

The stability analyses performed taking into account the interaction phenomena between the soil and the environment through a multidisciplinary approach have shown that changes in climate can affect the behaviour and stability of artificial and natural slopes, and engineers are required to assess their impact to minimize risk and financial costs. A range of approaches, from simple limit equilibrium solutions to advanced finite-element and finite-difference analyses, have been proposed in the literature to assess slope–atmosphere interaction and the related slope stability. In particular, the research related to thermo-hydro-mechanical coupling processes in soils has, nowadays, provided advanced modelling tools to take account of climate effects in slope stability and serviceability analyses. However, these tools are seldom used in practice, given their mathematical complexity and the in-depth knowledge of soil constitutive behaviour required for their application. The paper contributes to the dissemination of these modelling approaches by providing a schematic description of all the alternative procedures, accounting step by step for all the ingredients useful to represent the different slope conditions and processes. A review of the available methods to model the slope behaviour is presented, which allows, through the comparison of the methods and of their results, identification of their benefits and limitations and the different input requirements. In this perspective, the paper outlines the work by some of the participants in the COST Action TU1202 concerning the effects of climate and vegetation on slope deformation and stability performance. In all cases the review emphasizes the importance, for any reliable modelling of the slope behaviour, of simulating the fundamental soil hydromechanical properties, as derived from in situ monitoring and/or laboratory tests.

The scenarios of slope conditions that have been analysed are large in number and the classical soil mechanics formulation has been progressively generalized to account for a broader range of phenomena and soil behaviour, including thermohydraulic, hydromechanical and thermo-hydro-mechanical coupling. A common key aspect of all the numerical approaches described in the paper is the implementation of the climate effects, from the simplest way of simulating saturated conditions with suctions calibrated to represent real or monitored porewater pressures, to the simulation of the hydraulic behaviour of partially saturated soils (through SWRC and k(s)) and the implementation of rainfall and evapotranspiration. From these hydraulic analyses it is usually found that the soil permeability is the single most influential parameter in the infiltration and evapotranspiration processes in the slope. This finding highlights the importance of suitable methods to measure this property in situ. In conjunction with the hydraulic modelling, thermal coupling or mechanical coupling can be included in more sophisticated simulations of slope–atmosphere interaction problems. When considering the hydromechanical coupling, slope progressive failure can be predicted. For strain-softening soils, appropriate constitutive models implementing this aspect of the shear behaviour have to be used. Moreover, the soil permeability can significantly affect the rate of progressive failure, as higher permeabilities allow for a greater seasonal porewater pressure change at depth, thus leading to larger seasonal displacements and increased strain-softening. In addition, the effects of desiccation cracking and corresponding increase of permeability are not fully considered in most of the numerical approaches adopted in practice. Therefore, the development and practical application of advanced numerical tools for the prediction of swelling–shrinkage cracking at the slope surface represents a future research challenge.

Finally, the effects of temperature are usually not fully incorporated in the modelling. The temperature is not only a fundamental ingredient to analyse climatic processes at the interface between the slope and the atmosphere, but it can also affect the mechanical behaviour of soils. Therefore, thermo-hydro-mechanical coupling represents the most advanced numerical approach to solve the full set of balance equations. Coupled thermo-hydro-mechanical modelling has been successfully used for the analysis of deep geological storage and disposal of high-level nuclear waste (Gens 2010), but its application to slope stability problems is much less common (Nishimura et al. 2009) and does not appear to be well established yet. Nevertheless, its use should be promoted further and it is hoped that, with the increase of computational power and resources, coupled thermo-hydro-mechanical simulations of slope–atmosphere interaction will be more widely adopted in practice.

This paper is an output of Working Group 1 of EU COST Action TU1202 – Impacts of climate change on engineered slopes for infrastructure. TU1202 comprises four working groups: WG1, Slope numerical modelling; WG2, Field experimentation and monitoring; WG3, Soil/vegetation/climate interactions; WG4, Slope risk assessment. Outputs from each working group have been submitted to the Quarterly Journal of Engineering Geology and Hydrogeology and are intended to be read as a thematic set.