The importance of permeability in modelling soil-atmosphere interaction

. Soil-atmosphere interaction has been attracting increasing interest as the seasonal variation of pore water pressures (pwp) has been linked to a variety of geotechnical problems (e.g. slope stability and serviceability, foundation subsidence or swelling, desiccation cracking etc.) or has been identified as part of the solution of geotechnical problems (e.g. in sustainable urban drainage systems). Prediction of how the pwp will change within soils of low permeability under the combined effect of evapotranspiration and precipitation requires adequate knowledge of the soil permeability and how it varies spatially (e.g with depth) and temporally (e.g. with suction or degree of saturation, void ratio or due to the opening and closing of desiccation cracks). Nonetheless, in-situ measurements of permeability that satisfy both the spatial and temporal variation are difficult. In order to clarify the importance of variable permeability in predicting pwp variations under atmospheric loads, a series of one-and two-dimensional finite element analyses was performed, where the permeability model and the variation of permeability were parametrically studied. The results demonstrated that the variation of permeability, as well as the model employed in the analysis, e.g. allowing or not for desiccation cracking, influenced the values of suction calculated as well as the pwp profile with depth, highlighting the importance of estimating the spatial and temporal variation of permeability with some level of confidence.


Introduction
Mechanical soil properties, such as stiffness and strength, have traditionally monopolised the interest of practicing geotechnical engineers and researchers. The effect of pore water pressures (pwp) and of their variation on the stability and serviceability of geotechnical structures [1][2][3] and natural slopes [4,5], and on the coefficient of earth pressure at rest [6], has also long been recognised. Although, by association, the influence of hydraulic parameters, such as permeability and soil-water retention, is widely accepted in theory [7,8], in practice it is rare that complete and comprehensive laboratory and field data are dedicated to determining hydraulic parameters with accuracy and confidence, perhaps owing to the difficulty of the task.
Tsiampousi et al. [9] discussed the effects of hydraulic parameters on the serviceability of an infrastructure cut slope, primarily focusing on the soilwater retention (SWR) model and demonstrating its relevance on capturing numerically the expected behaviour of the slope. They also concluded that modelling permeability as a function of suction or of degree of saturation made little difference to the results. This paper studies the influence of the spatial and temporal variation of permeability on the computed pwp, via means of Finite Element (FE) analysis, with the ultimate goal of demonstrating its importance and emphasising the practical need to quantify it with reasonable accuracy. * Corresponding author: aikaterini.tsiampousi@imperial.ac.uk

Problem Definition
The analysis follows the case study of a slope cut in London clay (LC) in Newbury, South England, reported in Smethurst et al. [10], with the aim of utilising the field measurements of pore water pressures for assessing the FE analysis and its results. Smethurst et al. documented monthly rainfall data and daily potential evapotranspiration data for the calendar years 2003-2008, which match the pwp monitoring period. Pwp were presented for two vertical sections along the slope. For brevity, only the field data from a section close to the middle of the slope were considered here for comparison with the numerical data.
Two types of analyses were performed. Initially, a column of soil was considered under conditions that imposed 1D water flow, in order to study the effect of permeability and of the permeability model under simplified conditions. Subsequently, the slope itself was modelled allowing more realistic 2D conditions, which included the combined effect of soil-atmosphere interaction and pore water pressure dissipation following the excavation of the slope. Fully coupled hydro-mechanical analyses, accounting for the unsaturated nature of the weathered layers of LC [e.g. 11,12], were carried out in all cases with the FE program PLAXIS 2D [13]. While in some of the analyses only standard features of PLAXIS were employed, in others user-defined flow models were adopted which overwrite the permeability calculated by the program.

Column Analysis
The column analysis consisted of a 20m deep 1m wide FE mesh comprising a 17m deep layer of unweathered LC, overlain by 3m of weathered LC. The weathered layer was further subdivided to include a top 1.5m rooted zone and simulate the effect of plant roots on permeability [14].
The mechanical properties of the three layers were identical and modelled with a standard Mohr-Coulomb (MC) model, employing a Poisson ratio, μ, of 0.3 and a Young Modulus, Ε, of 30 MPa at the ground surface and increasing by 5 MPa for every metre of depth. Cohesion, c', was equal to 7 kPa and the angle of shearing resistance, φ', was 23 o . The increase of strength with suction was indirectly accounted for through Bishop's effective stress, while any potential collapse on wetting was disregarded. Stiffness was assumed to be independent of suction and constant throughout the analysis.
The unweathered LC, which has an entry value of suction significantly higher than the suctions expected in the analysis [15], was modelled as fully saturated and a constant permeability of 0.319E-3 m/day was assigned. The weathered and rooted layers, however, were modelled as unsaturated. The same van Genuchten [16] SWR curve was assumed for both layers, with α = 0.0044, n = 1.02, m = 2-1/n and a residual degree of saturation of 7%. To model the effect of desaturation on the value of permeability, the Mualem [17] approach was followed, 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 . Two analyses were considered, adopting two different values of (which were the same for the weathered and the rooted layers) as summarised in Table 1, Desat.A and Desat.B, enabling the effect of desaturation to be studied (see also Figure 1 (a)). Two further analyses were performed, where only the permeability of the rooted zone was altered to account for desiccation cracks, according to the model by Nyambayo [18]: where is the saturated value of permeability sustained until the tensile minimum total principal stress reaches the input value of 1 , and is maximum possible value of permeability reached when = 2 . This is not a standard feature of PLAXIS and was implemented into it as a user-defined flow model. The model effectively relates tensile stresses to the opening/closing of desiccation cracks by appropriately increasing/decreasing the value of permeability. Note that this value is then subjected to desaturation, based on Mualem's model. Two scenarios regarding desiccation were considered. The relevant parameters are shown in Table 1 and the obtained variations in Figure 1 Stresses were initialised in the analysis based on a unit weight of 20 kN/m 3 for the three layers, an assumed initial ground water table (GWT) at 1m of depth and an initially hydrostatic pwp variation with depth, allowing suctions to develop above the GWT. Monthly rainfall and evapotranspiration rates from Newbury were then applied in a combined form (i.e. net inflow/outflow) on the ground surface on a monthly basis simulating soilatmosphere interaction during the calendar years from 2003 to 2008. The infiltration boundary condition available in PLAXIS was employed, allowing a maximum possible value of suction of 1500 kPa and a minimum of 0 kPa. 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. The horizontal displacements at the two vertical boundaries and the horizontal and vertical displacements at the bottom boundary were fixed. The two vertical boundaries were considered impermeable effectively simulating 1D flow through the column.

Slope Analysis
The plane strain slope analysis follows the geometry of the Newbury cut slope (16 o and 27 m wide at the crest). The presence of Lambeth group clay and sand and of Thanet sand underneath the LC, which was neglected in the column analysis, was included in the slope analysis. The FE mesh adopted is shown in Figure 2.
The mechanical and hydraulic parameters of the LC layers were the same as in the column analysis. The mechanical behaviour of the underlying layers was also simulated with a standard MC model with μ = 0.3 and E = 30 MPa (increasing by 1.75MPa/m depth) for the Lambeth group clay, and 220 MPa (constant) for the Lambeth group and Thanet sands. c' = 8 and 1 kPa, and φ' = 28 o and 30 o for the clay and sands, respectively. Lambeth group clay was modelled as consolidating material with a constant permeability of 0.2096E-4 m/day, whereas the sands were considered to be drained.
Stresses were initialised in the analysis based on a unit weight of 20 kN/m 3 for all layers, an assumed initial ground water table (GWT) at 1m of depth and an initially hydrostatic pwp variation with depth, allowing suctions to develop above the GWT. The undrained excavation of the cut slope was performed first. As the slope was cut in 1997, and field measurements commenced in 2003 and continued until 2008, 12 years of soil-atmosphere interaction in a hydro-mechanically coupled analysis were subsequently simulated. The same rainfall and evapotranspiration data used in the column analysis were applied monthly on the initial ground behind the crest on the excavation and on the newly formed face of the slope (cycled twice to extrapolate from 6 to 12 years). Seepage was allowed at the bottom of the excavation, where a highway was constructed (i.e. water could flow out of -but not intothe FE mesh as the excess pore water pressures generated during the undrained excavation dissipated, and could not pond on the FE mesh). At the beginning of the soil-atmosphere interaction phase of the analysis, the properties of the top 3m of the sloping ground were altered to simulate a rooted and a weathered zone of 1.5m depth each, assuming that some degree of weathering occurred as the face of the excavation was exposed and vegetation established itself. Smethurst et al. reported a 3m deep zone of weathered LC at the crest of the slope, but it is unlikely that the face of the slope remained unaffected by the presence of vegetation and its exposure to the elements. Although it is expected that weathering at the slope would be less pronounced than weathering of the initial ground surface, which was exposed for a substantially longer period, the simplification was made here that weathering had the same effect on the permeability irrespective of location. Four analyses were run making the same assumptions as for the column analyses regarding desaturation and desiccation (see Table 1).
The horizontal displacements at the two vertical boundaries and the horizontal and vertical displacements at the bottom boundary were fixed throughout the analysis. The left-hand-side vertical boundary, which is an axis of symmetry, was impermeable. The right-hand-side vertical boundary was also impermeable allowing pwp along it to change in response to the applied inflow/outflow rate. Seepage was allowed at the bottom of the mesh where the Thanet sand interfaces with the underlying chalk.

Column Analysis
The column analysis results are shown in Figures 3 to 6 in terms of pwp plotted against depth, in comparison to the field data from the Newbury cut. Only the analysis and field data for the end of dry periods are considered here (September or October, depending on when the field measurements were obtained) whereas the data for the end of wet periods (usually May) have been disregarded for visual clarity. During the wet periods, suctions almost disappeared, hence the associated results are of less interest. As monthly rainfall and evapotranspiration rates were input in the analysis, the numerical results that are presented in the figures correspond to the end of month that fits the field data better (e.g. end of October for measurements taken on the 21 st of October, and end of September for measurements taken on the 28 th of September). Although the dates of field measurements and the analyses dates do not coincide, given the simplifications made in the analyses with respect to the rainfall and evapotranspiration rates, the numerical results are indicative of the field measurements.
All analyses gave generally satisfactory results in that they captured the range of suctions measured in the field and the annual variation in suctions observed. What was less successfully captured was the variation of suction in the top 1.5m in relation to the next 1.5m (from a depth of 1.5 to a depth of 3m). The field data from 21-Oct-03 show a very clear change of slope at 1.5 of depth. A change of slope can also be observed for the field data from 23-Oct-05 and from 28-Sept-06, indicating that the depth of the rooted layer of LC in the analysis was appropriately chosen. The suctions for 03-Oct-03 and 22-Sep-07 are smaller than the previously mentioned field data and there is no marked change in slope.
Comparing between the numerical results for the two desaturation cases, it can be argued that Desat.B, for which a faster reduction of permeability with desaturation was simulated in comparison to Desat.A (Figure 1), produced better results, certainly for a depth range between 1.5 and 3m: whereas Desat.A overpredicted the suctions at these depths for all but the 22-Sep-07 case, Desat.B overpredicted only the 03-Oct-06 data. At shallower depths, predictions were not as accurate. In particular for the shallowest depth at which pwp were measured, Desat.B overpredicted the field measurement taken on 21-Oct-03 but underpredicted all other measurements.
Similar observations can be made for the desiccation analyses, with Desic.B, where desiccation had a smaller effect on permeability for the values of suction obtained (Figure 1), giving overall better results at shallower depths.

Slope Analysis
The slope analysis results are compared to the same field data in Figures 7 to 10. Interestingly, the picture is reversed from what was observed for the column analyses, in that it is now Desat.A and Desic.A, rather than Desat.B and Desic.B, which resulted in a better simulation of the measured suctions at shallower depths. Overall, it can be argued that it is Desic.A that gives the best match with the field observations. Whereas flow was 1-D in the column analysis, it was 2-D in the slope To demonstrate the potential effect of anisotropic permeability on the analysis results, Desic.B was repeated assuming an anisotropy ratio of 100 between the horizontal and vertical permeability of London clay [19]. For the rooted zone, where the desiccation model was adopted, it was assumed that desiccation cracks were predominantly vertical. This meant that the horizontal permeability remained constant throughout the analysis, whereas the vertical permeability varied with the tensile total principal stress in the same manner as for Desic.B. This required a modification of the Nyambayo [18] model used in Desic.A and Desic.B (where both the vertical and horizontal permeabilities changed by the same amount during desiccation and behaviour remained isotropic), so that the horizontal and vertical permeabilities were calculated independently. The prediction of suctions ( Figure 11) at shallow depths improved considerably, demonstrating the importance of anisotropic permeability in 2D flow conditions.

Conclusions
The column analyses results indicated that faster reduction of permeability with desaturation improved the numerical simulation of suctions within the lower part of the weathered layer. They also indicated that the numerical reproduction of suctions within the rooted zone was better for the case where desiccation had a smaller effect on permeability. In any case, simulating the effect of desiccation improved the suction prediction within the rooted zone. The slope analyses results highlighted the potential importance of anisotropic permeability. Although there are several features of this simplified analyses that could be improved, such as the time resolution from months to days and employing a hysteretic SWR model (e.g. [20]), it is undisputed that permeability, its value and its variation during the analysis plays a major role on the computed results. This numerical exercise underlines the importance of determining with some level of accuracy not only the bulk permeability of a soil layer subjected to soilatmosphere interaction, but also the spatial and temporal variation of permeability due to desaturation and the effect of vegetation and weathering.