Soil-atmosphere interaction in unsaturated cut slopes

Interaction between atmosphere and soil has only recently attracted significant interest. Soil-atmosphere interaction takes place under dynamic climatic conditions, which vary throughout the year and are expected to suffer considerable alterations due to climate change. However, Geotechnical Analysis has traditionally been limited to simplistic approaches, where winter and summer pore water pressure profiles are prescribed. Geotechnical Structures, such as cut slopes, are known to be prone to large irreversible displacements under the combined effect of water uptake by a complex vegetation root system and precipitation. If such processes take place in an unsaturated material the complexity of the problem renders the use of numerical analysis essential. In this paper soil-atmosphere interaction in cut slopes is studied using advanced, fully coupled partially saturated finite element analyses. The effect of rainfall and evapotranspiration is modelled through sophisticated boundary conditions, applying actual meteorological data on a monthly basis. Stages of low and high water demand vegetation are considered for a period of several years, before simulating the effect of vegetation removal. The analysis results are presented with regard to the serviceability and stability of the cut slope.


Introduction
It has long been shown that pore water pressure variation affects the stability and serviceability of infrastructure slopes (cut and embankment slopes). Pore water pressures may vary due to consolidation or due to seasonal variations in the water balance under the combined effect of evapotranspiration and rainfall infiltration. It is possible that the two mechanisms act simultaneously in a number of cases. The complexity of the problem is such that use of numerical tools is essential.
Additional complications arise when the soil is unsaturated, as the behaviour of unsaturated soils is fundamentally different to the behaviour of fully saturated soils. The implications of partial saturation may be significant and Geotechnical structures in unsaturated soils may exhibit the reverse behaviour of what is expected in saturated soils. For example, Tsiampousi et al. [1] have shown that the factor of safety against failure of slopes excavated in unsaturated soils does not generally decrease with time, as in saturated soils (see fig. 1). Due to their unpredictability, careful numerical analyses are of essence when unsaturated soils are involved.
The numerical analysis of a slope excavated in an unsaturated soil is presented in this paper. Pore water pressures change in the analysis due to the combined effect of consolidation and soil-atmosphere interaction. Different stages of vegetation type and density are considered and the associated effects on the stability and serviceability of the slope are studied. The finite element (FE) code ICFEP [2], which adopts a modified Newton-Raphson solution technique with an error controlled substepping stress point-algorithm, was used. First, boundary conditions relevant to soil-atmosphere interaction are briefly explained. The boundary value problem to be solved is then presented in detail and the results of the analysis are discussed. The analysis results suggest that dense, high water demand vegetation enhances the stability of the slope in expense of serviceability. Vegetation clearance, however, may induce slope instability.

Figure 1.
Variation of the factor of safety against failure of a slope cut in saturated and unsaturated soil (after [1]).

Vegetation
Extraction of water due to vegetation is modelled in ICFEP with the root water uptake model (RWUM), described by Nyambayo & Potts [3]. The potential evapotranspiration rate, T p , needs to be defined for each increment of the FE analysis. From this, a sink term, S max , is calculated and incorporated in the continuity equation of fluid flow [2]. S max is assumed to vary linearly with depth, r, as shown in Figure 2 (a), according to the following relationship: where r max is the root depth, which is an input parameter and below which root water uptake is zero. The above equation relates to how much water could potentially be extracted if the supply of moisture in the ground were constant. To calculate the actual evapotranspiration, which is a fraction of the potential one, the sink term is multiplied by a suction-dependent parameter α [4]. The variation of α with suction is illustrated in Figure 2 (a). Suctions S1 (anaerobiosis point), S2, S3 and S4 (wilting point) are input parameters.

Precipitation
The precipitation boundary condition is a dual boundary condition which enables the simulation of rainfall infiltration and runoff; it may operate as an infiltration condition, by specifying a constant inflow rate, q n , or as a constant prescribed pore water pressure condition, p fb [5]. Values of q n and p fb have to be prescribed for every increment of the analysis that the boundary condition is to be active. At the beginning of such increments the pore water pressure at boundary nodes is compared to p fb : if found to be more tensile, an infiltration boundary condition is imposed, with a flow rate equal to q n ; if found to be more compressive, a constant pore water pressure equal to p fb is applied. As conditions may change during an analysis increment and a switch from one boundary condition to the other may be required, an automatic incrementation algorithm (described in Smith, 2003) is available in ICFEP.

Factor of safety calculation
To assess the stability of the cut slope, the factor of safety (FoS) against failure was calculated with the strength reduction technique explained in Potts & Zdravkovic [6]. Furthermore, the variation of the FoS with time was obtained following a similar approach to the one used by Tsiampousi et al. [1]: secondary analyses were initiated at appropriate increments of the main fully coupled analysis but consolidation was switched off and pore water pressures were not allowed to change (drained analyses). In this way, the FoS was obtained for a particular time instance. To obtain the variation of FoS with time the same process was repeated at various, appropriately selected time instances.

Geometry
A 10 m deep excavation in an unsaturated silty clay was considered ( fig. 3a). The geometry of the cut slope was discretised using isoparametric, quadrilateral 8-noded solid elements, resulting in the FE mesh shown in Figure  3 (b). The problem was considered in plane strain in a fully coupled consolidation analysis. Two displacement degrees of freedom were assigned at each node of each element of the FE mesh. Pore water pressure degrees of freedom were only assigned to corner nodes.

Soil properties and model parameters
The Barcelona Basic model [7], with modifications by Georgiadis et al. [8] and Tsiampousi et al. (2013-HV), was used to simulate the mechanical behaviour of the soil (model parameters in Table 1).The same properties were employed in a similar study by Tsiampousi et al. [9] and correspond to a silty soil tested by Estabragh & Javadi [10]. For the hydraulic behaviour, a van-Genuchten (1980) type soil-water retention curve (SWRC) was used, in combination with a permeability model in which the logarithm of permeability varies linearly with suction, from an initial saturated value, k sat , corresponding to suction s A , to a final value k min , corresponding to suction s B (Table 2). Finally, the actual transpiration properties adopted are summarised in Table 3.

Initial stresses
At the commencement of the analysis, soil stresses were initialised employing a unit weight of 19.1 kN/m 3 above and below the ground water table (GWT). The coefficient of earth pressure at rest K 0 was 2.1 at the surface, reducing to 0.6 at 15m below ground level. The GWT was 1 m deep (see also fig. 3a) and the pore water pressure profile was hydrostatic, with suctions developing above the phreatic surface.

Climate data
The necessary precipitation and potential evapotranspiration data were obtained from a meteorological station in Greenwich, London ( Figure 4). They refer to average long-term monthly data obtained from 1971 to 2000.
As discussed in the following section, the analysis is divided into five phases, each considering a different stage of vegetation, varying from low to high water demand, and finally simulating clearance of the vegetation on the slope. For this reason, the vegetation boundary condition was altered appropriately in each phase, with regard to both the root depth and the potential evapotranspiration rates. Similarly, the precipitation boundary condition varied between phases (most notably during the excavation phase). In total, the precipitation and vegetation boundary conditions (BC) were employed in four different manners (precipitation and vegetation BC1 to BC4 in Table 4). For simplicity, the climate data in Figure 4 formed the basis for all the rates applied in the five phases, as summarised in Table 4, e.g. the potential evapotranspiration rates T p in vegetation BC2 were reduced to 50% of those in Figure 4 in order to simulate low water demand plants.

Initialisation
The first phase of the analysis refers to the time prior to the excavation, assuming that the ground surface was covered in dense, high water demand vegetation. Its purpose is to establish a pore water pressure and stress regime which would have been affected by the precipitation and evapotranspiration, prior to the excavation. Figure 5 (a) illustrates the boundary conditions applied on all four boundaries of the FE mesh during this phase. The area to be excavated subsequently is highlighted in grey. The climate boundary conditions were applied on the top horizontal boundary of the mesh. Information about the vegetation and precipitation boundary conditions can be found in Table 4. It should be highlighted that the precipitation and evapotranspiration rates correspond to those in Figure 4. They were applied monthly and the typical year, starting in April (as in Figure 4), was repeated 9 times to reproduce 9 years. This time was found to be adequate for an annual pore water pressure variation, repeating itself year after year, to be established.

Excavation
The excavation was performed at the beginning of the 10 th year of the analysis, in a period of time short enough for essentially undrained conditions to be maintained. The hydraulic boundary conditions were altered to those illustrated in Figure 5 (b). The changes in the precipitation boundary condition are of particular interest, notably the fact that q n was switched to zero (see also  Table 4). Additionally, the vegetation boundary condition was inactive during this phase. The applied boundary conditions ensure that any further changes of the pore water pressures during the excavation were solely due to unloading and not to soil-atmosphere interaction. T p As in Fig.  4 50% of rates in Fig. 4 As in Fig.  4 10% of rates in

Low water demand vegetation
The subsequent phase refers to years 10 to 19 of the analysis, i.e. it extends for 10 years following the excavation. The vegetation and precipitation boundary conditions were reactivated, as illustrated in Figure 5 (c). Dense vegetation was preserved at a distance of 10m from the slope crest (vegetation BC1 in Table 4). The remaining area and up to a vertical distance of 1.5m from the bottom of the excavation was assumed to be covered in low water demand vegetation (e.g. shrubs and bushes). This corresponds to vegetation BC2 in Table 4.
The precipitation boundary condition on the horizontal ground beyond the crest of the slope was restored to its form prior to the excavation (i.e. precipitation BC1). On the slope only 50% of the precipitation rates in Figure 4 were applied (precipitation BC4), assuming that a drainage system was in place, capable of capturing and removing 50% of rainfall. Finally, at the bottom of the excavation the precipitation rate was set to zero. No further changes were made to the precipitation boundary condition for the remainder of the analysis.

High water demand vegetation
At the beginning of the 20 th year of the analysis, growth of high water demand trees on the slope is simulated. The potential evapotranspiration rates applied on the slope are equal to those in Figure 4, with the root depth increasing gradually from 0.5 to 2m over a period of 5 years (vegetation BC3 in Table 4). This and the remaining boundary conditions are schematically explained in Figure 5 (d). Overall, this phase of high water demand vegetation on the slope lasted for a total of 10 years (i.e. years 20 to 29 of the analysis).

Vegetation clearance
Finally, at the end of March of Year 29 the vegetation on the slope was cleared and replaced by low and sparse vegetation (vegetation BC4 in Table 4).The total duration of the final stage of the analysis was once more 10 years (i.e. years 30 to 39 of the analysis).

Stability
The FoS was calculated bi-annually, at the end of March and August, for the 30 years following the excavation. Its variation with time, measured in months from the end of excavation, is shown in Figure 6, where the symbols correspond to results obtained from FoS analyses.
With reference to Figure 6 and for the low water demand vegetation stage, the values of FoS for August are systematically higher than the values for March. However, an overall decrease in the FoS is observed, which can be attributed to the combined effect of slopeatmosphere interaction and the dissipation of tensile excess pwp's developed during undrained excavation. By the end of this phase of the analysis, the values of FoS have more or less stabilised and oscillate about the value of 1.4. This is possibly due to the fact that average monthly rainfall and potential evapotranspiration rates were applied year after year. Although evapotranspiration data do not vary greatly from year to year, this is not so for rainfall data, which exhibit a large annual variation. It is reasonable to expect that a particularly dry year would have enhanced the stability of the slope, whereas a particularly wet year could have brought it to failure. This cannot be captured when average rainfall data are employed.
The FoS visibly increased during the next phase, where the slope was covered in high water demand vegetation. The enhanced stability can be justified by the larger suctions that have developed in comparison to the previous phase, as expected, and which contribute to a strength increase. Additionally, the values of FoS for August are systematically higher than the values for March, however, the overall increasing trend observed signifies that the beneficial effect of suction build-up during the summer months is not entirely lost during the subsequent winter months.
When the vegetation was cleared from the slope, a steady and fast decrease in the FoS was obtained. Rainfall rates are now high enough to quickly eliminate the suctions that have developed in the previous phase. It is interesting to note that the FoS no longer increases temporarily during summer and within five years the stability of the slope becomes marginal.
The formation of a slip surface about 4.5 m deep is indicated by the horizontal displacement profiles along a vertical line at mid-slope, shown in Figure 7. These displacements are sub-accumulated from the end of the previous phase and therefore represent the displacements developed due to vegetation clearance.    There is little difference between the August and March displacement profiles for both Year 19 and 29. Furthermore, although the displacements are not uniform across the bottom of the excavation the maximum difference is only about 5 mm. Serviceability problems have been reported in highly vegetated slopes excavated in fully saturated soils. The discrepancy may be due to the increased stiffness unsaturated soils generally exhibit.

Conclusions
This paper presents the results of a numerical study on the effect of soil-atmosphere interaction on the behaviour of a slope cut in an unsaturated soil. Under the combined effect of vegetation, precipitation and excess pore water pressure dissipation, the factor of safety exhibited an overall decrease when the slope was covered with low water demand vegetation. A seasonal fluctuation in the factor of safety was nonetheless visible. A significant increase in the overall factor of safety was obtained when high water demand vegetation was established on the slope. This type of vegetation is known to cause serviceability problems in fully saturated soils. However, differential settlements in the present analysis remained small, possibly due to the increased stiffness unsaturated soils generally exhibit. Finally, vegetation clearance was shown to cause stability issues, with the factor of safety reducing fast and a deep seated failure mechanism developing within five years.