Numerical study of the application of capillary barrier systems for prevention of rainfall-induced slope instabilities

. Slope instability is often caused by decreases in suction due to heavy and prolonged rainfall. In this study, the application of capillary barrier systems (CBSs) for suction control and slope stabilization purposes (i.e. reducing the risk of rainfall-induced slope instabilities) is analysed, due to their capacity to limit the percolation of water into the underlying soil. The behaviour of two slopes was studied numerically: a bare slope made of fine-grained soil and the same slope covered by a capillary barrier system. The time evolution of suction in the slopes subjected to realistic atmospheric conditions was studied by performing numerical finite element analyses with Code_Bright. In particular, multi-phase multi-physics thermo-hydraulic analyses were performed, modelling the soil-atmosphere interaction over periods of many years. Suction and degree of saturation distributions obtained from these analyses were then exported to the software LimitState GEO, which was used to perform limit analysis to assess the stability of the slopes. The CBS was able to limit the percolation of water into the slope and was shown to be effective in increasing the minimum values of suction attained in the underlying ground, resulting in improved stability of the slope.


Introduction
In unsaturated conditions, the presence of matric suction s, defined as the difference between pore-liquid pressure p l and pore-gas pressure p g (i.e. s=p l -p g ), imparts higher strength to the soil, compared to fully dry and fully saturated conditions. However, suction may vanish or greatly decrease after prolonged and heavy rainfall. The stability of slopes is often guaranteed by the effect of suction. In these cases, a heavy rainfall event may cause significant reductions in suction and shear strength in the soil and induce slope instability [1].
Capillary barrier systems (CBSs) are geotechnical structures made of an upper finer layer (F.L.) overlying a lower coarser layer (C.L.), placed over the ground with the aim of preventing the percolation of water into the underlying soil (U.S.) [2]. The coarser layer is typically at very low degree of saturation and, consequently, the corresponding unsaturated hydraulic conductivity may be several orders of magnitude lower than that of the finer layer. Thus, prior to significant water breakthrough to the coarser layer, rainwater is stored in the finer layer whereas the coarser layer acts almost as an impermeable barrier. This water can then be removed by evapotranspiration [3] and, if the barrier is sloped, by lateral drainage [4]. The barrier fails when the amount of water stored in the F.L. is so high that the suction at the interface between F.L. and C.L. decreases down to the "bulk water-continuity value" of the coarser layer, at which the hydraulic conductivity of the C.L. starts increasing significantly [5]. In this condition, water breakthrough occurs from the F.L. to the C.L., and eventually into the underlying soil.
CBSs have been primarily employed as landfill covers [6]. More recent research efforts have addressed the applicability of CBSs for suction control purposes and slope stability [7], given their ability to prevent water percolation into the underlying soil. However, more advanced numerical studies are required in order to clarify different aspects of the problem (e.g. the effect of complex soil-atmosphere interaction, the effect of the different parameters and long-term behaviour).
An advanced numerical modelling procedure was developed in this study, which links thermo-hydraulic finite element analyses (including modelling of soilatmosphere interaction) and mechanical limit analyses. This modelling procedure was applied to the study of the behaviour of a bare slope and a slope covered with a CBS, both subjected to realistic weather conditions.

Modelling procedure
The long-term behaviour of a bare slope (BS) and of the same slope covered by a CBS (CS), both subjected to real atmospheric conditions, was studied numerically. The modelling procedure was divided into three steps. 1.
Advanced finite element (FE) coupled thermo-hydraulic analyses were performed to establish temporal and spatial variations of suction s and degree of saturation S l . 2.
Values of the product of suction and degree of saturation s·S l were then interpolated to a new grid. 3.
The new grid values of s·S l were imported into computational limit analysis (LA) software to perform stability analyses considering the effect of unsaturated conditions on shear strength.

Finite element (FE) model
The long-term behaviour of a bare slope (BS) (see Fig.  1(a)) and of the same slope covered by a CBS (CS) (see Fig. 1(b)) was studied numerically using the FE code Code_Bright [8]. Thermo-hydraulic, multi-phase and multi-physics two-dimensional simulations were performed. For the mass balance, advective liquid flow and diffusion of water vapour within the gas phase were included. Gas pressure was considered uniform and constant, equal to the atmospheric pressure, p g =100 kPa. For the energy balance, conductive heat flow and convective heat flow (i.e. heat flux associated to mass fluxes of liquid and vapour) were modelled. No mechanical behaviour was modelled (all materials were assumed rigid). The slope was 6 m high with a slope angle of 35°. The CBS, made of a 40 cm-thick finer layer overlying a 20 cm-thick coarser layer, covered the slope in all its height. A drain was modelled at the toe of the sloping CBS in order to simulate the eventual collection of water diverted by the CBS.
The underlying soil (U.S.), the coarser layer (C.L.) and the finer layer (F.L.) were representative of a silt, a gravelly sand and a fine sand, respectively. The principal laws and parameters used to model the materials are shown in Table 1. Default laws and parameters were used for the physical properties (e.g. liquid viscosity µ l , liquid density ρ l , gas density ρ g , water mass fraction in the gas phase ω w g ) and more details can be found in the Code_Bright User's manual. Unsaturated hydraulic properties of the different materials (soil water retention curve SWRC and soil hydraulic conductivity curve SHCC) are shown in Fig. 2. SWRC and SHCC of the U.S. were modelled using the Van Genuchten-Mualem model [9] whereas an improved model incorporating the influence of liquid film flow at low degree of saturation was used for the C.L. and F.L. [5]. Water retention hysteresis was included in all the models using a hysteretic bounding surface water retention model [5]. Lateral and bottom boundaries were modelled as impermeable to liquid and heat flux and they were placed sufficiently far from the slope in order not to affect the results. At the top boundary, "atmospheric" boundary conditions were applied. For the mass transfer, the atmospheric boundary conditions included: rain P, runoff R (occurring when p l at the ground surface is equal to the atmospheric gas pressure), and evaporation E. For the energy transfer, the atmospheric boundary conditions included: radiation R n , sensible heat flux (advection) H s and latent heat flux H c (convection).
The evaporation E was modelled as [10]: where k is Von Karman's constant (k=0.4), z a is the screen height (z a =1.5 m in this case), v a is the wind speed at the screen height, ψ is the stability factor (ψ=1), z 0 is the roughness length (z 0 =0.001 m for a surface covered by short grass), ρ va is the absolute humidity of the atmosphere at the screen height and ρ v is the absolute humidity at the soil surface (i.e. boundary nodes). ρ va is a is a function of atmospheric air temperature T a , relative humidity RH a and atmospheric gas pressure p ga, whereas ρ v is a function of soil surface temperature T, pore-liquid pressure p l and pore-gas pressure p g . These relationships are governed by the psychrometric law. The sensible heat flux H s was modelled as [10]: where ρ ga is the atmospheric gas density, C a is the specific heat of the gas, T a is the atmospheric temperature at the screen height and T is the soil surface temperature.
The atmospheric boundary condition applied in this study was representative of the climate of Cagliari (Italy), located in a relatively dry area of Europe but subjected to sporadic intense rainfall events. Historical weather data for years 1981-2010 were available for the weather station of Cagliari Elmas [11]. The corresponding monthly averages are represented by the histograms in Figs. 3(a-e) whereas the solid lines represent yearly sinusoidal distributions fitted to the monthly averages. Moreover, Fig. 3(f) shows daily rainfall values for a wet 10-year period (1984)(1985)(1986)(1987)(1988)(1989)(1990)(1991)(1992)(1993).
The numerical simulations involved three different intervals. Atmospheric temperature T a , atmospheric relative humidity RH a , radiation R n and wind speed v a were modelled using the yearly sinusoidal distributions shown in Figs. 3(a-d) throughout all the intervals. The applied rain was different for the different intervals. In the first interval, lasting 20 years, the CBS was not built yet, rainfall was applied to the underling soil surface in both BS and CS models using the yearly sinusoidal distribution shown in Fig. 3(e). The purpose of this first interval was to set up realistic initial distributions of s, S l and T in the underlying soil. In the second interval, lasting 10 years, the CBS was constructed for the CS model and the same yearly sinusoidal rainfall distribution was applied to the CBS surface with the Table 1. Constitutive laws and parameters used to model the materials in the numerical FEM analyses.  purpose of setting up initial realistic distributions of s, S l and T in the CBS. For the BS model, the same rainfall distribution was again applied to the underlying soil. In the third interval, lasting 10 years, the daily rainfall distribution shown in Fig. 3(f) was applied to the top boundaries (underlying soil for the BS model, CBS for the CS model) with the purpose of assessing the response of the two models to a particularly wet 10-year period, modelled with daily accuracy. Fig. 3. Atmospheric data from Cagliari Elmas weather station: (a) air temperature T a , (b) wind speed v a , (c) relative humidity RH a , (d) radiation R n , (e) average monthly precipitation P, (f) daily precipitation P for years 1984-1993.

Limit analysis (LA) model
Limit analyses (LA) were performed to assess the stability of the BS model and the CS model at different times using LimitState GEO [12]. This software adopts a discontinuity layout optimization method [13] to work out potential failure mechanisms and the corresponding factors of safety [12]. The geometry of the models was the same presented for the FE models. Lateral and bottom boundaries were modelled as fixed (i.e. no displacements allowed). Materials were modelled as rigid-perfectly plastic with Mohr-Coulomb yield criterion and associative plastic flow. Material properties (i.e. unit weight γ, effective cohesion c' and friction angle ϕ') used for the LAs are shown in Table 2. It should be noticed that the friction angles of F.L. and C.L., representative of a fine sand and a gravelly sand, are much higher than that of the U.S., representative of a silt. Small values of c' were introduced to avoid numerical instabilities. Unsaturated conditions were included in the yield criterion using the Bishop stress [14] with χ=S l , which although it is not able to represent properly all aspects of the mechanical behaviour of unsaturated soils, is adequate to model the shear strength [15]. The resulting shear strength τ can be expressed as: where σ is the total normal stress.
Values of suction s and degree of saturation S l were exported from the FE mesh nodes and converted to values of the product s·S l , which were then interpolated in a new regular grid which was imported into LimitState GEO.

Finite element (FE) analyses results
In this section, results concerning the response of the BS and the CS FE models to the applied weather conditions of the third interval of the analysis (i.e. daily rainfall distribution, years 1984-1993) are presented. Fig. 4 shows the degree of saturation contours of the two models at the end of the most critical rainfall event, which occurred 1 year and 57 days after the start of the analyses (i.e. 7 th of March 1985). In the absence of the CBS (see Fig. 4(a)), the soil surface was fully saturated down to a depth of approximately 80 cm, where a sharp wetting front separated the overlying saturated zone from the underlying zone at lower degree of saturation (S l ≈0.27). In the presence of the CBS (see Fig. 4(b)), the soil underlying the CBS was maintained at lower degree of saturation even during this most critical rainfall event, unlike the lateral areas not covered by the CBS. Within the CBS, the C.L. remained at very low degree of saturation (i.e. no breakthrough), whereas the lower part of the F.L. was at high degree of saturation (S l > 0.99). Fig. 5 shows the suction profiles in the U.S. at the end of the most critical rainfall event, for the BS and CS models, in three sections: top, middle and toe (indicated in Fig. 4). In agreement with what was observed for the degree of saturation, in all sections the presence of the CBS had a beneficial effect on the suction profiles in the underlying soil. Very low values of suction were attained close to the soil surface for the BS model (i.e. between 0 and 8 kPa) whereas much higher suction values were maintained in the CS model (i.e. higher than 1.7 MPa).
The beneficial effect of the CBS is also confirmed by the results presented in Fig. 6, where time histories of s and S l are plotted at the underlying soil surface of the middle section (see points O indicated in Fig. 4). In the BS model, although suction attained very high values during summer due to the effect of the evaporation (see Fig. 6(a)), zero suction values are attained several times E3S Web of Conferences 195, 01027 (2020) E-UNSAT 2020 https://doi.org/10.1051/e3sconf/202019501027 after intense rainfall events. By contrast, minimum suction values attained with CS model during winter are permanently higher than 1 MPa (see Fig. 6(b)). Similar concepts apply to the degree of saturation which, in the BS model, experienced high yearly fluctuations (i.e. very low during summer and very high during winter) reaching full saturation several times. By contrast, in the CS model, S l had a much smoother variation with values always under 0.35.  In order to understand how the CBS limited the percolation of water into the underlying soil, absolute liquid velocity v l and degree of saturation S l profiles in the CBS at the end of the most critical rainfall event are shown in Fig. 7 for top, middle and toe sections. It can be seen that, proceeding from the top to the toe of the slope, a higher amount of water is transported laterally down the slope by the F.L. of the CBS, in particular close to the interface with the C.L. as indicated by the increase of S l and v l . In contrast, no water is transported by the C.L. which acts only as an impermeable barrier. The water diverted laterally was ultimately collected by the drain placed at the toe of the sloping CBS. This mechanism qualitatively agrees with previous studies on the lateral water diversion capacity of sloping CBSs [4].

Limit analysis (LA) results
Limit analyses were performed to assess the slope stability considering three models: the bare slope model (BS), the slope covered by the CBS in which the stability of the CBS is assessed (CS -CBS) and the slope covered by the CBS in which the stability of the underlying soil is assessed (CS -U.S). In the latter, only the underlying soil was modelled with s and S l distributions obtained from the CS model. The stability is indicated by the factor of safety FoS, here defined as: where ϕ lim and c lim are respectively friction angle and effective cohesion which would cause collapse. Fig. 8 shows the failure mechanisms and the corresponding FoSs at the end of the most critical rainfall event for three models analysed. The BS model  5). The corresponding FoS is lower than 1 suggesting that the slope was unstable under the most critical rainfall event. The use of a CBS significantly improved the stability of the underlying soil, as shown in the CS -U.S. model (see Fig. 8(c)), due to the higher suction values attained. The failure line was much deeper and the corresponding FoS was enormously higher than 1 (FoS=78.8). For the CS model, the stability of the CBS became more critical than that of the U.S. Indeed, the failure mechanism of the model CS -CBS involved the F.L. of the CBS but the corresponding FoS was now higher than 1 (FoS=1.37) which suggests a stable condition. Fig. 9 shows the FoSs of the different models in 10 critical rainfall events identified in the 10 years analysed. It can be seen that in two events (i.e. 1.18 years and 3.10 years) the FoS of the BS was lower than 1. Introducing the CBS, the underlying soil was now permanently stable with corresponding FoSs always very high. The stability of the CBS is almost unaffected by the weather conditions, indicated by an almost constant trend of the FoS, always higher than 1.

Conclusions
Advanced FE coupled thermo-hydraulic multi-physics multi-phase numerical simulations considering the longterm soil-atmosphere interaction were coupled with mechanical limit analyses to assess the effect of the use of a capillary barrier system (CBS) as a slope cover. Due to its ability to prevent water ingress into the underlying soil and to divert water laterally, the CBS was able to maintain high values of suction and low values of degree of saturation in the underlying soil throughout the analysis, when subjected to realistic weather conditions (Cagliari, Italy). Consequently, the CBS was highly efficient in preventing rainfall-induced slope instability. The problem becomes thus controlled by the stability of the CBS, typically made of materials with high friction angle and affected by fewer uncertainties (i.e. predictable shear strength parameters).