Wetting-induced instabilities: triggering mechanism and predisposing factors

. This paper addresses the topic of unsaturated soil stability focusing on wetting-induced instability mechanisms occurring before the attainment of the classical shear failure. The study was conducted by simulating a series of imbibition tests on elementary soil volumes by controlling water content and net stresses. The behaviour of an ideal coarse-grained soil was described by combining the WR2-Unsat model for the solid skeleton and the Van Genuchten relationship for the definition of the water retention curve. The model reproduces the soil response upon wetting, modelling the activation mechanism of the volumetric instability and identifying the factors that most contribute to its triggering.


Introduction
Natural soils close to the ground surface are often in an unsaturated state. The interaction with the atmosphere and, in particular, the rainfall infiltration and the following water redistribution can cause the triggering of instability phenomena.
The most common instabilities develop as shear failures occurring when the applied shear stress exceeds the available shear strength. These instabilities are due to the reduction of matrix suction and, thus, shear strength resulting from water infiltration [1][2]. Additional hydraulically induced instabilities can occur before attaining shear failures for positive values of the hardening modulus [3][4][5]. These instabilities manifest as sudden and significant collapses of the solid skeleton accompanied by uncontrolled growth of pore water pressures and soil saturation. These phenomena can play a crucial role either in triggering shallow landslides or in the development of displacements of foundations and earthworks [6][7][8].
This paper is devoted to the analysis of the second instability mechanism, focusing on: i) the causes activating the unstable process; ii) the identification of the main predisposing factors to the instability.
The numerical study was conducted using the finite element code Abaqus/Standard. The ideal wetting tests were performed on elementary soil volumes under the usual geotechnical testing conditions by controlling water content and total stresses. The soil behaviour was described using the WR2-Unsat model [9], consisting of an extension to unsaturated conditions of a recently proposed constitutive law for saturated structured soils [10]. The retention behaviour was described by employing the well-known Van Genuchten model, establishing a oneto-one relationship between suction and saturation * Corresponding author: giada.rotisciani@gmail.com degree. The numerical model reproduces the soil response during wetting, modelling the activation of the volumetric instability. The simulations cover a wide range of initial conditions and soil properties, allowing the identification of the factors that most contribute to triggering the instability.
The paper is organized as follows. In Section 2, the constitutive laws adopted for the solid and fluid phases are described together with the followed calibration procedure. Section 3 is devoted to the analysis of the activation mechanism of the volumetric instability. Section 4 presents the results of a parametric study aimed at investigating what extent soil properties and initial stress state can favour the occurrence of unstable processes.

Constitutive model
The hydro-mechanical soil behaviour was described by combining the WR2-Unsat model [9] and the Van Genuchten relationship for the definition of the Water Retention Curve (WRC).
The WR2-Unsat model is a new critical state-based hardening plasticity model obtained by extending to unsaturated conditions the model proposed by Boldini et al. [10]. The model reproduces the main features of the response of structured unsaturated soils under monotonic loading conditions, adopting a single-surface elasto-plastic formulation. Thanks to its hierarchical formulation, the various features of the model can be switched on or off, depending on the specific geomaterial to be modelled. In the following, a simplified version of the model is presented, specializing the constitutive laws to i) unstructured unsaturated soils and ii) axisymmetric testing conditions. Major details on the generalized model formulation can be found in [9,10], the last work only for the case of dry or fully saturated soils.

Model formulation
The WR2-Unsat model is formulated in terms of Bishop's effective stress , assuming χ equal to the saturation degree : , being the net stress, s the matrix suction and the Kronecker delta.
The reversible soil behaviour is described by the hyperelastic formulation proposed by Houlsby et al. [11] accounting for the non-linear dependency of the elastic stiffness on effective stress. The adopted free energy potential W has the following expression: where and are the volumetric and deviatoric elastic strain invariants, p a is a reference pressure corresponding to the atmospheric pressure, n e , k, and g are dimensionless parameters to be calibrated experimentally.
The yield locus YS delimiting the elastic domain follows the expression proposed by Bigoni and Piccolroaz [12]: where Φ is defined as: Φ=p′∕p c ; p′ and q are the mean effective and deviatoric stress invariants;  and m control the shape of the curve and N the pressure sensitivity; p c is the hardening variable. The plastic potential function G is characterised by an equation similar to that adopted for the yield surface: where is defined as: = p′∕̂ ; ̂ identifies the dimension of the G-locus passing through the current stress state; , , and define the shape of the plastic potential in the stress invariant plane.
It is worth noticing that selecting different values for the shape parameters of the YS-and G-loci implies the assumption of not associated flow rule.
The evolution of the hardening variable depends on plastic deformations and degree of saturation via the hardening rule. The latter differs from the one originally proposed by Boldini et al. [10] for the introduction of an explicit dependence of p c from S r : * (4) where is the volumetric plastic strain increment, * is the modified compression index evaluated in the bilogarithmic compressibility plane, b controls the mechanical effects induced by the variations of the degree of saturation. The latter parameter has a clear physical meaning since it defines the relative distance in the compressibility plane between the Normal Consolidation Line (NCL) related to the current degree of saturation and the saturated one.
The retention behavior was described using the well known Van Genuchten WRC: where is related to the air-entry value, n w and m w control the shape of the curve with 1 1 ⁄ ; S r,res is the residual saturation degree. Eq. (5) defines an uncoupled and not hysteretic retention behavior.

Model calibration
The calibration of the constitutive functions is based on the literature data related to coarse-grained soils with a non-plastic fine.
The set of constitutive parameters adopted for the numerical simulations is listed in Table 1.
The dimensionless stiffness modulus g was determined according to the expressions of the maximum shear stiffness modulus available in literature, placing n e =0.5 and p a =1 kPa ( [13]).The dimensionless bulk modulus k is expressed as a function of g and Poisson's ratio, here assumed equal to 0.2.
The distorted shape of the yield locus shown in Fig.  1 was defined by back-analyzing the results of undrained compression triaxial tests carried out on nonplastic materials, including sands and sandy silts ( [14]). The plastic potential function was assumed to be coincident to the ellipse of the Modified Cam Clay model. The stress ratios at the maximum deviatoric stress of the YS-and G-loci are  p =0.9 and N g =1.331, respectively; the latter corresponds to the slope of the Critical State Line (CSL) in the (q, p′) plane.
Soil compressibility was determined coherently to the experimental data related to sands, silty sands, and sandy silts ( [15]). The hydro-mechanical coupling parameter b was estimated by interpreting a large set of experimental data within the conceptual framework adopted in this model ( [16]).
Finally, the WRC employed in the study is characterized by low air-entry value, small slope, and negligible residual degree of saturation. The values for the WRC's parameters stem from the analysis of numerous wetting curves related to sands and silty sands ( [17][18]).

Activation mechanism of volumetric instability
The response of elementary soil volumes upon wetting was investigated by simulating a series of ideal imbibition tests with the finite element code Abaqus/Standard. The tests were performed under the usual axial-symmetric loading conditions by controlling water content and total stresses. More specifically, a water inflow with constant intensity was prescribed at the boundaries, keeping unchanged the applied total stresses. It is worth noticing that the imposed boundary conditions mimic the ones expected in situ during rainfall events. Conversely, they are not representative of the ones imposed in laboratory tests, where suction and net stresses are usually controlled.
The initial conditions for the imbibition tests and some soil properties are summarized in Table 2.
Consider now the reference analysis corresponding to Test 1. This test was performed under a total vertical stress of 22.5 kPa, corresponding to the weight of a homogenous soil column 1.25 m high. The initial suction of 10 kPa is representative of suction distribution observed in wet seasons in Southern Italy, just below the ground surface. The material is lightly overconsolidated with an overconsolidation ratio R of 1.1 and a stress ratio  0 of 0.8. R is defined as the ratio between p c and the intersection of the yield surface passing through the current stress state with the horizontal axis. Finally, the initial void ratio was determined once known the initial stress state and the position of the NCL related to the current S r (procedure described in detail in [19]).
The increase in pore saturation entails a gradual loss of suction and, thus, a reduction of the mean effective stress under constant deviatoric stress (0F, Fig. 2a). The initial shrinkage of the yield locus due to the increase in S r leads quickly the stress state in the elastic-plastic regime. Despite the contractive volumetric behaviour, the elastic domain continuously decreases in volume following the stress state (segment 0P). Table 1. Parameters adopted for the numerical simulations.

Type Value
Elastic parameters ne = 0.5 pa = 1 kPa k = 6666.7 g = 5000 Upon exceeding the stress ratio  p , the trend switches, and the yield locus starts to increase in size (segment PF). This response results from the development of positive plastic volumetric strains that counterbalance the softening induced by S r changes.
The soil response is shown in Fig. 2b in terms of the evolution of pore water pressure u w and volumetric strains  v with e w , being e w the water ratio defined as the ratio between the water and solid volumes. The water ratio is related to the water content and is, thus, a control variable of the test. The abrupt change in the slope of u w marks the transition from elastic to elasto-plastic response. The relationship between e w and u w keeps linear up to u w -values around -3 kPa. Upon exceeding these values, the curve changes continuously in slope, and u w starts to increase quickly. The rapid and uncontrolled growth in the pore water pressures points out the loss of stability of the material response, causing the immediate stop of the simulation.
The material exhibits a contractive volumetric behaviour. Appreciable volumetric strains start to develop when the soil is in the elasto-plastic regime.  v increases gradually with e w up to the e w -value of 0.67, corresponding to u w -levels ranging from -3 and -2 kPa. Upon exceeding these values, the strain rate starts to increase rapidly, making evident the loss of control of the test.
The similarities observed in the curves (u w , e w ) and ( v , e w ) stem from the mutual interaction between the two response variables. The S r increase due to the imposed water inflow and the consequent suction removal cause skeletal deformations that, in turn, affect the current values of suction and saturation degree ( [5]). The wetting process becomes unstable when the soil cannot sustain anymore the applied loading. This circumstance is pointed out by continuous changes in soil particle configuration and ever-increasing pore water pressures.  The occurrence of volumetric instability is strictly connected to the constraints imposed in the test. Under the control of net stress and suction, the same material never exhibits an unstable response. Pore saturation grows following the increase of suction, never going out of control ( [4][5]). The instability remains, thus, in the latent form. This evidence justifies the choice to perform the laboratory wetting tests under the control of suction and net stresses.
It is important to note that the instability activated by the suction removal can trigger static liquefaction once attained fully saturated conditions. Static liquefaction is a well-known instability phenomenon occurring in undrained conditions under stress states lying within the potential liquefaction domain. The latter gathers all stress states with  included between  p and N g . This phenomenon manifests with a sudden decrease in the soil strength accompanied by the development of large pore water pressures that reduce up to vanish the effective stresses in the soil.

Predisposing factors to instability
The reference analysis shows the activation of a wetting instability before the complete removal of suction. This form of instability depends on the interaction between hydraulic and mechanical variables that can lead to a loss of control of pore saturation and soil deformations.
The set of numerical analyses summarized in Table  2 aims to identify what factors most contribute to the occurrence of this instability. To this goal, the ideal imbibition tests were repeated by modifying, one at a time, the parameters controlling the hydro-mechanical coupling and the initial state of the material. More specifically, the selected b-and  * -values cover the ranges typically observed in silty sands: (1, 5) for b, (0.01, 0.05) for  * . The total vertical stress varies between 22.5 and 36 kPa, corresponding to the weight of a soil column 1.25-2 m high. The initial stress ratio ranges from 0.6 to 1; such interval includes the stress ratio at the peak of the yield locus. The suction, degree of saturation, and overconsolidation ratio were kept unaltered in all the tests. To the shortage of space, the test results related to the influence of hydraulic variables and their relationship (WRC) sit outside of this publication. Figure 3 shows the impact of the b-value on the evolution of u w and  v with e w /v 0 . The latter variable defines the variations of water ratio normalized to the initial specific volume, corresponding to the applied perturbation.

Soil properties
As revealed by Table 2, Tests 1, 2, and 3 differ only in the initial specific volume that increases with the bvalue, coherently to the increasing distance between the saturated NCL and the one related to S r0 .
Regardless of the b-value, the instability takes place in all the tests with an uncontrolled growth of pore water pressures and volumetric strains before the complete removal of suction. Nevertheless, the hydro-mechanical parameter has a marked influence on i) the moment of activation of the instability and ii) the magnitude of the external perturbation triggering the loss of control.
The increase in the b-value introduces an everstronger interplay between mechanical and hydraulic variables, favouring instability. In highly collapsible soils, the injection of a water volume involves significant volumetric collapses that magnify the effects of the external perturbation in terms of S r changes, leading to early instability. As confirmed by the results depicted in Fig. 3, the increase in the b-value entails: i) the occurrence of the volumetric instability under lower u w -levels, that is, more distant to full saturation; ii) the reduction of the injected water volume triggering the instability. The curves reported in Fig. 4 refer to Tests 1, 4, and 5, differing in the value of  * . The latter controls the slope of the NCL in the bi-logarithmic compressibility plane, thus the compressibility of the solid skeleton in the elasto-plastic regime. Keeping unchanged the initial stress state, the reduction of this parameter entails the increase of the initial specific volume.
The soil response is unstable for  * -values greater than 0.03, whereas it remains stable when the  * -value is equal to 0.01. In the latter case, indeed, the simulation ends when the soil is fully saturated. In the other cases (Tests 1 and 4), the simulation ends when the changes in pore water pressure and volumetric strains are no longer controllable. This circumstance manifests i) under u wlevels that decrease as the  * -value increases and ii) after the injection of a water volume lessening as soil compressibility rises.
These results confirm the crucial role played by the *-value in the occurrence of volumetric instability. The increase in soil compressibility favours the wetting collapse, strengthening the coupling between the hydraulic and mechanical variables. This feature, combined with the reduction of the initial specific volume, increases the risks associated with the occurrence of an unstable process, reducing the magnitude of the triggering perturbation and lowering the threshold values beyond which the control is lost (suction, pore saturation, specific volume). Finally, it is worth noticing that whenever the control is lost, the current stress state lies in the potential liquefaction domain. This means that in the analysed cases, the instability induced by the suction removal could give rise to static liquefaction once attained full saturation.

Initial conditions
In this sub-section, the influence of the initial stress state on the soil response is investigated by modifying the ratio between the two stress invariants, namely  0 (Tests 1, 6, and 7), and the magnitude of single stress components, by acting on the total vertical stress  v0 (Tests 1, 8, and 9).
Keeping unchanged  v0 , the reduction of  0 entails the increase of the mean effective stress and the reduction of the deviator stress. The changes in the initial effective stress state involve opposite variations in the initial preconsolidation pressure p c (R is kept equal to 1.1) and specific volume. Figure 5a shows the wetting-induced stress path and the final yield locus related to each test. In all the selected cases, the stress path moves horizontally and ends before attaining the classical shear failure. The final stress ratio  f differs in the three tests: the higher  0 , the higher  f . In any case,  f exceeds  p .
The responses of the soil samples are shown in Fig.  5b in terms of variations of pore water pressures and volumetric strains with the applied perturbation. In all the tests, the injection of water volume provokes an increase in pore saturation, the collapse of the solid skeleton, and the growth of pore water pressure. Nevertheless, there are significant differences in the shape of the curves. Starting from high stress ratios (Tests 1 and 6), the curves exhibit ever-increasing gradients until  v , and u w are no longer controllable. This evidence reflects the occurrence of unexpected instability modes taking place at u w -levels decreasing with increasing  0 . Conversely, in Test 7, u w and  v increase more slowly, and their control is never lost.
The simulation ends just after the complete removal of suction.
These results can be explained keeping in mind that the soil response results from the combination of two factors: the reduction of p′ due to the suction removal and the softening induced by the increase in S r . While the current stress ratio keeps being lower than  p , these factors have opposite effects on volume changes. When  exceeds  p , both factors force the development of positive plastic volumetric strains: the higher the stress ratio, the stronger the tendency to collapse upon wetting, and the higher risks of the volumetric instability. Tests 1, 8, and 9 differ in  v0 that increases from 22.5 to 36 kPa. Keeping unchanged  0 , the rise of the vertical component involves a proportional increase of the horizontal one. Passing from Test 1 to Test 9, thus, both stress invariants increase as well as the preconsolidation pressure, whereas the initial specific volume lightly decreases.
Similarly to the previous cases, the soil responses are reported in the planes (u w ,e w /v 0 ) and ( v ,e w /v 0 ) in Fig.   6. The impact of  v0 is appreciable only in the final parts of the tests when the curves are no longer overlapped. The increase of the vertical stress improves soil stability, limiting the wetting-induced volumetric collapse and the resulting changes in pore water pressures. Similarly to the reference test, Test 8 ends before completing the saturation process. Nevertheless, the delay in the occurrence of instability and the increase of the injected water volume make readily apparent the beneficial effect of the  v0 increase. In Test 9, this effect is even more evident since the soil response remains stable up to fully saturated conditions. These results can be interpreted taking into account that the increase in stress ratio induced by suction removal lessens as total and, thus, effective stresses grow. Therefore, increasing  v0 , the soil exhibits an ever-weaker tendency to collapse that makes unlikely the activation of instability in the unsaturated regime.

Concluding remarks
This paper presents some results of a numerical study devised to reproduce the activation mechanism of wetting instabilities occurring before the attainment of the classical shear failure. The study was conducted by employing the WR2-Unsat model for reproducing the behaviour of the solid skeleton and Van Genuchten model for describing the WRC.
The instability manifests as an uncontrolled growth of pore water pressures and skeletal deformations that causes an abrupt acceleration of the saturation process. Depending on soil properties, state of the material and drainage conditions, this instability can trigger a second unstable mechanism under fully saturated conditions, known as static liquefaction.
The results of this study point out that: -The instability is driven by the strong interplay between hydraulic and mechanical variables; -The activation of this phenomenon is markedly dependent on the soil properties controlling the hydro-mechanical coupling. The higher the b-or  * -values, the higher the risks triggering volumetric instability; -The activation of this mechanism is significantly dependent on the initial stress state. High stress ratios favour instability, exacerbating the tendency to the collapse upon wetting. Conversely, the increase in the magnitude of the stress state has a beneficial effect on soil stability.
These results suggest paying attention when water infiltration involves highly collapsible soils located just below the ground surface under high deviatoric stresses. These conditions are often satisfied in unsaturated slopes consisting of shallow silty sands layers with slope angles close to the critical state friction angle.