3D effects of soil-atmosphere interaction on infrastructure slope stability

. It has long been established that pore water pressure (pwp) variations affect the stability and serviceability of slopes. Pwp variations may be due to consolidation/swelling processes within soils of low permeability or may be due to seasonal evapotranspiration and precipitation processes. The simultaneous study of such phenomena and of hydro-mechanical coupling in sloping ground requires use of advanced numerical methods. Often, infrastructure slopes are considered in plane strain (two-dimensional) conditions for simplicity and to date this is the case for most numerical analyses considering soil-atmosphere interaction. However, this approach makes it impossible to study the longitudinal extent of a possible slip surface forming following vegetation removal. Three-dimensional, fully-coupled numerical analyses of a cut slope were performed herein to study the effect of vegetation clearance on stability and explore numerically ways of implementing effective vegetation management.


Introduction
Engineered infrastructure slopes [1][2][3], i.e, slopes formed by excavating or compacting soil to accommodate highways, railways etc., and natural slopes [4,5] have long been known to be affected by interaction with the atmosphere due to the seasonal variation of pore water pressure [6], which affects the stresses, stiffness and the coefficient of earth pressure at rest [7]. For example, [2] observed that modestly vegetated slopes are prone to deep seated failures, primarily after wet periods, whereas densely vegetated slope suffer from severe serviceability problems, which are likely to occur towards the end of dry periods.
The effects of vegetation and of vegetation management on the stability and serviceability of cut slopes in unsaturated soils were studied by [8]. Conducting hydro-mechanically coupled, plane strain numerical analyses, they showed that following a relatively quick excavation, its stability, as indicated by the Factor of Safety, slightly reduced with time under the competing effect of swelling and soil-atmosphere interaction, but without ever reaching failure, when the slope was covered in shrubs. When trees were allowed to establish themselves on the slope, its stability steadily increased with time but on the expense of serviceability, as the onset of differential displacements along the bottom of the excavation, where the highway or railway asset would be located, started being observed. Finally, vegetation clearance led to a quick loss of stability and during the following wet season it became marginal, highlighting the fact that careful vegetation management is required to avoid both serviceability and stability problems.
* Corresponding author: aikaterini.tsiampousi@imperial.ac.uk As the analyses in [8] were plane strain, it was not possible to study either the out-of-plane extent of the failure mechanism formed after vegetation clearance, or vegetation management that included spacing different vegetation types and species in the out-of-plane direction. A preliminary study to this effect is presented here. A slope cut in London clay (LC) was considered in 3D fully coupled hydro-mechanical analyses, accounting for the unsaturated nature of the weathered layers [e.g. 9,10]. All analyses were carried out with the FE program PLAXIS 2D [11].

General
The analysis loosely follows the geometry, stratigraphy and soil properties of a slope cut in London clay (LC) in Newbury, South England, reported in [12] and previously studied by [13].
A 50m deep layer of LC was considered. The top 3m were weathered. The 10m deep excavation was 50m wide at the crest and 30m at the bottom, giving rise to a 2:1 slope, which is typical of slopes cut in LC. The Finite Element (FE) mesh employed in the analyses is shown (after the excavation was performed) in Figure 1.
The mechanical properties of the weathered and unweathered LC were modelled as linear elastic perfectly plastic, employing a standard Mohr-Coulomb (MC) model. A Poisson ratio, μ, of 0.3 was assumed and a Young Modulus, Ε, of 30 MPa was set at the ground surface, increasing by 5 MPa per metre depth. Cohesion, c', was 7 kPa and the angle of shearing resistance, φ', was 23 o .
As the unweathered LC has an air entry value of suction which can exceed 1MPa [14], i.e. is significantly larger than the suctions expected in the analysis, it was modelled as fully saturated. A constant value of permeability of 0.3E-3 m/day was employed for the whole layer. This was a simplification in comparison to [13], where permeability was a function of mean effective stress, made in order to reduce the computational cost of the 3D analyses. The value of permeability in the current study was representative of the value of permeability at 10m depth (i.e., equal to the excavation depth) in the earlier study.
In contrast, the weathered LC, was modelled as unsaturated, with a van Genuchten [15] soil-water retention curve defined by parameters α = 0.0044, n = 1.02, m = 2-1/n and a residual degree of saturation of 7%. Its permeability was equal to 3.7E-3 m/day at full saturation and reduced due to desaturation according to the Mualem [16] approach, with the difference that the effective degree of saturation was not raised to the power of ½ but to a power equal to an input parameter which was set equal to 0. The unit weight of LC was 19.1 kN/m 3 . The initial ground water table (GWT) was assumed to be 1m deep and pwp varied hydrostatically with depth, allowing suctions to develop above the GWT. The coefficient of earth pressure at rest was 1.2.
The excavation was performed at the beginning of the analysis in five undrained phases, in each of which 2m of soil were removed. Next, five years of soilatmosphere interaction were considered: both the flat ground behind the excavation and the slope were covered in trees. Monthly rainfall and evapotranspiration rates from [13] were used to establish a typical year which was cycled five times to simulate the five years of soil-atmosphere interaction. They were then applied in a combined form (i.e. net inflow/ outflow) through the infiltration boundary condition available in PLAXIS. A maximum possible value of suction of 1500 kPa and a minimum of 0 kPa were allowed. The infiltration boundary condition is a dual boundary condition changing automatically from prescribed inflow/outflow to prescribed pore water pressure (head) when the above suction limits are exceeded. It was assumed that a drainage system capable of capturing and removing 50% of the rainfall was installed on the slope, similar to [13]. At the bottom of the excavation seepage was allowed, so that water could flow out of the FE mesh as the excess pore water pressures generated during the undrained excavation dissipated, while water could not pond on the boundary. The vertical boundaries were impermeable (axes of symmetry), with the exception of the out-of-plane righthand-side vertical boundary, where seepage was allowed allowing pwp along it to change in response to the applied inflow/outflow rate at the top boundary of the FE mesh. No change in pwp was prescribed at the bottom boundary where the interface with the permeable chalk is found. The horizontal displacements at the four vertical boundaries and the horizontal and vertical displacements at the bottom boundary were fixed. In the subsequent phases clearance of the vegetation was simulated, considering different scenarios as detailed subsequently.  Clearance along multiple 10m strips 1. 16 5 Clearance in a check pattern 1.11 6 Clearance along a horizontal strip 1.2

Vegetation clearance
Vegetation clearance was modelled by altering the infiltration rates in the analysis at selected areas on the slope to simulate different scenarios. Leaving the monthly rainfall rates unaltered, the monthly evapotranspiration rates were reduced to 10% of the rates corresponding to trees, following the approach by [13], to obtain rates more representative of grass. Six different scenarios of vegetation clearance were considered, as summarised in Table 1 and Figure 2. In the first one, vegetation was removed from whole of the slope, resembling a plane strain analysis (although it was still carried out in 3D). In the second scenario, a 10m wide vertical strip centred halfway along the slope in the out-of-plane direction (y coordinate of 50m) was cleared. Next, in scenario 3, the width of the strip was increased to 20m (centred at y of 55m). In scenario 4 multiple 10m strips were cleared, leaving 10m of untouched vegetation in between them, and in 5 clearance followed the check pattern shown in Figure 2. Finally, in scenario 6, a horizontal strip of 4m width was cleared.
The factor of safety, , was then calculated following a phase where ′-′ reduction was performed until failure was achieved. The corresponding failure mechanism was established by plotting vectors of incremental displacements at the end of the analysis (e.g., Figure 3). During the calculation pwp were no longer allowed to change: time was effectively frozen and the value of was calculated for the corresponding combination of stresses and pwp, as in [13]. Table 1 summarises the values of calculated for the six vegetation clearance scenarios studied.
was also calculated at the end of the 5 th year of soil-atmosphere interaction, before vegetation removal, as a benchmark.

Scenarios 1, 2 and 3
The factor of safety, , at the end of the 5 th year of soilatmosphere interaction was equal to 3.17. Although this is not directly comparable to the computed by [13], as different stages were considered in the earlier analysis in relation to the present one (e.g., soil-atmosphere interaction prior to the slope excavation and 10 years of the slope being vegetated with shrubs before trees established themselves), the value falls within the range of expected values. In [13], after 5 years of treedominated soil-atmosphere interaction, increased by about 1.2. Considering that it was 2.1 after the excavation (which, however, started from pwp already depressed because of the previous soil-atmosphere interaction), the value of 3.17 in the present analysis is not unreasonable.
This value reduced to 1.06 in scenario 1, where vegetation was cleared from everywhere on the slope (i.e., in an equivalent plane strain manner), in good agreement with [13]. When a 10m wide strip was cleared of its vegetation in scenario 2, reduced to 1.56 and when the strip was 20m wide, was 1.1, indicating that the latter scenario approached the plane strain equivalent scenario whereas the former one did not. Figure 4 presents the normalised incremental horizontal displacements in the x-direction (see Figure  1 for directions of axes) along a vertical line through the centre of the cleared area, i.e., halfway up the slope and at half the width of the area stripped of vegetation, for scenarios 1 to 3, at the last converged increment of the analysis. Zero incremental displacements signify stable areas, whereas non-zero incremental displacements signify areas of movement and therefore the plot is indicative of the depth of the failure mechanism at this particular vertical line. Although their relative magnitudes are indicative of the failing mass, the actual magnitudes of incremental displacements are irrelevant as failure occurs and the FE analysis becomes unstable. Therefore, they have been presented in Figure  4 normalised by the maximum value obtained in each analysis. It can be observed that the failure mechanism at this particular location for the plane strain equivalent analysis (scenario 1) was about 5m deep, which is similar to [13]. For scenarios 2 and 3 (10 and 20m wide strips cleared of vegetation, respectively) the depth of the failure mechanism was about 9 and 6m, i.e., once again scenario 3 approached the plane strain condition whereas scenario 2 did not.
The horizontal extent of the failure mechanism can be investigated by examining the incremental horizontal (in the x-direction) displacements along the toe of the slope. They are presented for scenarios 2 and 3 in Figure  5, normalised by the maximum value obtained in each analysis. Naturally, in the plane strain equivalent analysis of scenario 1 the extend of the failure mechanism was equal to the extent of the FE mesh (100m). In scenario 2, the centre of the clearance strip was at a y-coordinate of 50m, which coincides with the location of the maximum incremental horizontal displacement along the toe. The centre of the 20m wide strop was at a y-coordinate of 55m (see also Figure 2), which also coincided with the location of the maximum incremental horizontal displacement for this analysis (hence the two analyses have their maxima at different locations). Interestingly, although the failure mechanism was just short of 20m wide in scenario 3, approximately coinciding with the width of the cleared strip (20m), in scenario 2 the width of the failure mechanism exceeded the width of the strip (10m), signifying that areas where the vegetation had been left untouched participated in the failure mechanism. Figure 6 explores this further. The normalised incremental displacements and the suctions along the toe of the slope for scenario 2 (10m wide strip) have been plotted in the same graph for comparison. It is very clear that the suctions were zero from y = 45m to y = 55m, i.e., for the extent of the cleared strip. It is also clear that the failure mechanism extended beyond the cleared area and mobilised the soil outside the cleared strip, where suctions were as high as 18 kPa. This is further demonstrated in Figure 7, where the contour of zero suction and the vectors of incremental displacements at failure are superimposed on the 3D geometry of the numerical model.
This was not the case for scenario 3 (20m wide strip), as illustrated in Figure 8: the extend of the horizontal failure mechanism is shown to be just short of the width of the cleared strip where zero suctions were obtained.
The results indicate that clearing a 20m wide strip of land across the slope reduced the factor of safety to levels comparable to clearing the whole slope of its vegetation. The failure mechanism that occurred was similar in depth to the plane strain equivalent analysis and its horizontal extent did not exceed the width of the cleared area, i.e., this was wide enough to fit the whole extent of the failure mechanism.
When a 10m strip of land was cleared across the slope, the factor of safety was not reduced to critical levels. Nonetheless, was failure to occur, the analysis results indicated that the likely failure mechanism would not be limited by the presence of trees on either side of the cleared strip. Instead, it would potentially extend horizontally to areas beyond the cleared zone, where suctions would resist its formation, leading to mobilisation of soil strength deeper into the ground in comparison to the 20m wide strip. This explains the deeper failure mechanism observed in Figure 4.

Scenarios 4, 5 and 6
Scenarios 4 and 5 were inspired by scenario 2 and explored different combinations of 10m wide areas of vegetation clearance. In scenario 4, repeated 10m wide strips, 10m apart, were considered. As the failure mechanism extended into the vegetated area in scenario 2 discussed above, it is not surprising that the distance of 10m between strips was not adequate in maintaining the factor of safety, which dropped to 1.16 from 1.56 in scenario 2. Indeed, Figure 9 demonstrates that the whole slope was involved in the failure mechanism and the contribution of the vegetated strips of land is limited to the slightly higher factor of safety obtained in comparison to the plane strain equivalent analysis (1.06).
This insight could be only have been provided by a method of analysis that allows the failure mechanism in 3D to be reliably obtained. For example, had a 10m wide strip of the slope been cleared in a field experiment where failure did not occur (e.g., because of a winter drier than normal), the extent of the failure mechanism would not have been revealed and it would perhaps be tempting to assume that repeated strips of similar width at regular intervals of the same width could be adopted in vegetation management, invertedly reducing the factor of safety beyond what was intended.
Similar conclusions can be derived from the analysis of scenario 5, where a check pattern of vegetation clearance was examined, each rectangle of the pattern being 10m in width. In the same manner as the horizontal spacing of the strips did not help maintain the factor of safety, the vertical spacing failed to stop the failure mechanism from progressing, perhaps owing to its significant depth. The corresponding was only 1.11. Scenario 6, where a horizontal strip of land close to the slope toe was cleared, gave the most satisfactory value of (1.2) from the three clearance scenarios. Figure 10 plots the normalised horizontal displacements in the x-direction along a vertical line at the centre of the slope (i.e. half way across and half way along the slope). It can be observed that a very deep failure mechanism was obtained, much deeper than the mechanism in scenario 1, which is also plotted for comparison. Figure  11 shows the suctions on the slope at failure. The pwp pressure was zero at the bottom of the excavation and at the toe and the zero pwp contour has climbed up the slope covering the whole vertical extent of the cleared area. Nonetheless, higher up the slope suctions remained high (up to 100kPa), pushing the slip surface of the failure mechanism deep into the soil. Scenario 6 can serve as the basis for future further studies on vegetation management, where not only the stability but also the serviceability of the slope should be taken into account.

Conclusions
The numerical study has provided preliminary but useful insight into 3D effects of soil-atmosphere interaction on slope stability which would have been difficult, time-consuming and costly to obtain through field and large-scale laboratory testing.
It was shown that clearing a 10m wide strip on the slope of trees may lead to a deep and wide failure mechanism, exceeding 10m in width and mobilising the soil in vegetated as well as in cleared areas, i.e. the width was not adequate for a failure mechanism to develop fully within it for the particular combination of geometry and soil properties considered. On the contrary, a strip of 20m was adequately wide to accommodate the full failure mechanism, with a depth and factor of safety approaching those for a plane strain analysis. It can be concluded that vegetation management should avoid 20m wide clearance strips and that perhaps 10m wide strips are more appropriate.
To avoid out-of-plane serviceability issues along the route of the railway or highway, a more generalised approach to vegetation management than clearing single strips should be adopted. The analysis results, however, demonstrated that neither spacing 10m wide strips at 10m horizontal distance between them nor following a check pattern of 10m wide areas were adequate in maintaining the factor of safety. Instead, clearing the toe of the slope of its vegetation may be a better approach and it could form the basis for further numerical investigation.
Future numerical studies should consider different geometries and soil properties, the possibility that a weathered layer develops with time on the newly formed face of the slope and its unsaturated nature, including hydraulic hysteresis (e.g. [17]) and unsaturated permeability. Numerical analysis can provide useful and cost-effective preliminary information about slope stability and guide field and large-scale laboratory investigations, which are necessary and irreplaceable, but require a lot more resources and therefore need to be well planned out in order to maximise their benefit.