Unsaturated structured soils: constitutive modelling and stability analyses

The paper presents a new single-surface elasto-plastic model for unsaturated cemented soils, formulated within the critical state soil mechanics framework, which should be considered as an extension to unsaturated conditions of a recently proposed constitutive law for saturated structured soils. The model has been developed with the main purpose of inspecting the mechanical instabilities induced in natural soils by bond degradation resulting from the accumulation of plastic strains and/or the changes in pore saturation. At this scope, the constitutive equations are used to simulate typical geotechnical testing conditions, whose results are then analysed in light of the controllability theory. The results of triaxial tests on an ideal fully saturated cemented soil and on the corresponding unsaturated uncemented one are first discussed, aiming at detecting the evidence of potentially unstable conditions throughout the numerical simulations. This is followed by similar analyses considering the combined effects of both the above features. For each analysed case, a simple analytical stability criterion is proposed and validated against the numerical results, generalizing the results, and highlighting the crucial role of state variables and model parameters on the possible occurrence of failure conditions.

Specific volume at p 0 = 1 kPa under isotropic conditions n Porosity n e Parameter controlling nonlinear elasticity n w Shape parameter of Van Genuchten retention curve (n/S r ) lim Threshold value of n/S r for partially saturated soils (n/S r ) lim,s Threshold value of n/S r for partially saturated structured soils p a Reference pressure for nonlinear elasticity p c Isotropic yield compression pressure p c,dummy Parameter controlling the dimension of the f g -locus passing through the current stress state p 0 Mean

Introduction
Natural soils can be structured, according to Burland [10] or Leroueil and Vaughan [37], and partially saturated. This is the case of structured soils located near the ground surface interacting with the atmosphere. The study of the hydro-mechanical behaviour of these deposits is relevant to many civil engineering applications, from the design of shallow and deep foundations to the stability analysis of slopes or road and railway embankments [4,22,64]. Despite the rather different origins, structure and partial saturation share some similar mechanical effects as both induce higher initial stiffness and strength which, at the microstructural scale, should be qualitatively related to the enhanced interparticle bonding due to cementation or capillary menisci, respectively. In many geotechnical boundary value problems, structure and partial saturation might act separately or in combination, leading to a wide range of different possible responses of the involved soil deposit.
In structured soil deposits, the enhanced interparticle bonding can be related to different processes as, for example, percolation of calcium carbonate or weak lithification. As mentioned above, its main engineering effect is the enhancement of the stiffness and strength of the soil.
Qualitatively, it can be assumed that, upon application of external loads, the stresses at certain interparticle contacts attain the bond strength, triggering the mechanical bond degradation (debonding) process. Bond degradation is an irreversible phenomenon which, based on experimental observations, results as mainly controlled by plastic strain accumulation (e.g. [2,3,34,48,57]). It is characterized by an initially more intense effect, which reduces with the accumulation of plastic strains, bringing the material towards its final fully destructured state.
In unsaturated soils, the capillary menisci existing at the inter-granular contacts have a stabilizing effect, providing an additional component of normal forces. Thus, the capillary forces induce a sort of hydraulic bonding, whose effects are similar to those activated by diagenetic bonds [39,63]. The experimental evidence shows that, starting from saturated conditions, an increase in suction causes a stiffening of the soil response and, in heavily over-consolidated soils, an increase in shear strength [66,69]. Furthermore, under compression tests, the unsaturated specimens are less compressible than the saturated ones, exhibiting an increased apparent preconsolidation pressure. Upon wetting, saturation degree increases and water menisci progressively vanish, triggering a hydraulic degradation process. Such a process is reversible and mainly controlled by the changes in the degree of saturation that can occur even for stress states within the elastic domain [16,18,20,72].
Despite the similar engineering effects induced by cementation and partial saturation, remarkable differences in the two phenomena exist. On one hand, cementation is progressively damaged by plastic straining of the soil, leading to a permanent and irreversible debonding process, which, at the macroscale, reflects into a shrink of the elastic domain in the stress space. Conversely, changes in the degree of saturation due to either mechanical loads, rainfall infiltration or evapotranspiration modify the capillary bonds between particles, resulting into either a bondinglike or debonding-like reversible process, which at the macroscale corresponds to an expansion or a shrink of the elastic domain, respectively.
In recent years, several constitutive models have been proposed to describe the mechanical behaviour of natural soils. Many of them are based on critical state soil mechanics and are formulated within the single-surface isotropic hardening-plasticity framework (e.g. [7,34,38]).
Their key ingredient is in the isotropic hardening rule, which typically accounts for strain-induced structure degradation by means of one or more negative exponential damage-type terms. Along the same track, a further improvement can be obtained accounting for early irreversibility, either by adopting a mixed isotropic-kinematic hardening multi-surfaces formulation and/or by implementing a bounding surface plasticity framework (e.g. [5,31,47,54,60]).
Despite the above-mentioned considerable advances in constitutive modelling, there are still few models able to capture the combined effects of structure-related bonding and partial saturation [9,24,36,42,46,62,70]. For this, in the following, a relatively simple model accounting for both structure degradation and partial saturation is proposed. The model is formulated extending to unsaturated conditions a recently proposed versatile constitutive formulation for saturated structured soils [7]. The model aims at reproducing the main features of the hydro-mechanical response of unsaturated cemented soils under monotonic loading conditions, adopting a single-surface isotropic hardening elasto-plastic formulation based on a single stress variable to account for partial saturation effects.
Similarly to the original model, the hardening laws describing the destructuration phenomena can easily be activated or deactivated depending on the observed experimental response of the soil under study.
In the following, the effects of structure degradation and partial saturation on the response of idealized natural soils are numerically investigated adopting the proposed model, separating first and then combining these two key features. In the proposed element test simulations, particular emphasis is placed on the failure conditions, with the aim of investigating the stability of the material by detecting the, so-called, loss of controllability domain. In fact, here we interpret the occurrence of actual or potential instabilities by means of the controllability theory, originally proposed by Nova [43] and then extended in Buscarnera and di Prisco [15], along the line of what discussed in Buscarnera and Nova [11,14], Buscarnera [11], Marinelli and Buscarnera [40], Buscarnera and Dattola [12] and in Rotisciani et al. [52]. Moreover, simple analytical criteria for capturing soil instabilities are proposed and validated against numerical results. Although these criteria only refer to the proposed constitutive formulation, they have the merit of clearly highlighting the fundamental role of the soil parameters and initial state of the material on the occurrence of instability conditions, as such deserving some discussion in the text.
The paper is organized as follows. The formulation of the constitutive law is presented in Sect. 2 while the basic ideas of the controllability theory are briefly summarized in Sect. 3. The numerical results of ideal triaxial test simulations are discussed in Sect. 4 to highlight the impact of structure degradation and partial saturation on the predicted hydro-mechanical response. For each simulated condition, an analytical criterion is provided emphasizing the factors that mainly affect soil stability. In the final section, some conclusions are addressed, together with possible future developments of the research.

Constitutive model
The mechanical behaviour of the solid skeleton is here described by employing a modified version of the WR2 model proposed by Boldini et al. [7], named WR2-Unsat. The original WR2 is a critical state-based hardeningplasticity model, aimed at reproducing the main aspects of the mechanical response of structured soils under monotonic loading paths. In this study, the constitutive law was modified following the approach proposed by Jommi [29] to account for the effects of capillary forces exerted on the solid skeleton in the partial saturation regime.
In the following, the basic constitutive equations are briefly outlined, referring to the triaxial plane where all the stress and strain paths investigated in this study are confined. Details on the original multiaxial formulation for fully saturated conditions and its predictive capabilities in reproducing experimental results can be found in Boldini et al. [7].
The model is formulated in terms of Bishop's effective stress r 0 ij , assuming v equal to the saturation degree, S r : where r ij is the net stress, s is the matrix suction defined as the difference between pore air u a and water u w pressures, and d ij is the Kronecker delta. Equation (1) coincides with the classical definition of Terzaghi's effective stress under fully saturated conditions. The yield surface YS delimiting the elastic domain is given by Bigoni and Piccolroaz [6]: where U is defined as: U ¼ p 0 =p c ; p 0 and q are the mean effective and deviatoric stress invariants; a and m define the shape of the curve in the meridian section; N controls the pressure-sensitivity and it is assumed to be constant; p c is the hardening variable governing the current size of the yield locus. Within the elastic domain, the soil behaves as a hyperelastic material with stiffness moduli increasing nonlinearly with both stress invariants. The elastic constitutive equations are deduced by differentiating the free energy potential W defined by Houlsby et al. [26,27]: where e e v and e e s 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.
The plastic potential function f g obeys to the following expression, similar to the one adopted for the yield locus: where U g is defined as: p 0 =p c;dummy ; p c,dummy has no physical meaning and only identifies the fictitious dimension of the f g -locus such that it passes through the current stress state; a g , m g , and N g define the shape of the plastic potential in the (q, p 0 ) plane. In this model, thus, the flow rule can be associated or not, depending on the values assigned to the shape parameters. The formulation of both yield surface and plastic potential function is very flexible and can be adapted to reproduce well-known criteria, such as non-associated Mohr-Coulomb or Modified Cam Clay (MCC) models (see [7]).
The initial dimension of the yield locus results from the combination of the effective stress and bond histories. Subsequent plastic straining and/or variation of the degree of saturation trigger a homothetic evolution of the yield domain, here accounted for by an isotropic hardening law. This latter differs from the one proposed by Boldini et al. [7] for the introduction to an explicit dependency on the saturation degree: where _ e p v and _ e p s are the increments of the plastic strain invariants, e d v is the controlling variable for the volumetric strain-induced de-structuring and is defined as: dt, e p s is the deviatoric plastic strain, k Ã is the modified compression index evaluated in the bi-logarithmic compressibility plane, g v ; n v ; g s , and n s are the volumetric and deviatoric destructuration parameters, _ S r is the increment in the saturation degree and b defines the relative position on the normal compression line (NCL) related to the current S r with respect to the saturated one.
In Eq. (5), the first term coincides with the hardening law of the MCC model for saturated soils expressed as a function of the modified compression index. The second and third terms describe the strain-induced structure degradation using two separate negative exponential damage-type forms. The last term represents the simplest form of the hydraulic hardening that allows capturing most of the observed phenomena for unsaturated soils as, for example, the increase of the preconsolidation pressure with suction, and the wetting-induced collapse [17,35,49,51]. It worth noting that this expression is valid only when S r is not far from the unity and water menisci are homogeneously distributed over the soil volume [13], as assumed in the following.
The modelling strategy employed for capturing hydraulic bonding effects is similar to that used in cemented soils. This choice is justified by the similarities observed experimentally between the mechanical behaviour of unsaturated and bonded soils. In this framework, the mechanical bond degradation is an irreversible process controlled by plastic strain accumulation, whose intensity exponentially decreases with the development of the plastic strains. Instead, the hydraulic bonding is a reversible process controlled by the changes in degree of saturation, which is considered as a measure of the overall distribution of the capillary forces responsible for the hydraulic bonding.
In this model, the effects of saturation and cementation are uncoupled. This is of course a simplification of the real soil behaviour and may represent a limitation in terms of predictive capabilities of the model when applied to those soils in which the above coupling plays a non-negligible role, as, for example, the weakly bonded unsaturated sapriolite discussed in [36].
The expression suggested by Van Genuchten [65] is adopted for describing the soil water retention curve (WRC): where a w is linked to the air-entry value, n w and m w control the shape of the curve, and S r,res is the residual saturation degree. For the sake of simplicity, suction is assumed to be a unique function of the degree of saturation, neglecting both hydraulic hysteresis and the dependence on the current void ratio.
The retention model given by Eq. (6) can be reformulated in a more convenient way for the analysis of the possible occurrence of instabilities as follows: where _ s Ã is the incremental hydraulic stress variable conjugate to _ S r in the definition of the second-order work d 2 W, _ s Ã ¼ n _ s [15], n is the porosity and _ s is the increment of matrix suction.

Controllability theory
Material stability is here discussed with reference to the controllability theory proposed by Nova [43]. This theory combines the Hill's stability criterion with the concept of test controllability and demonstrates that the crucial factor determining the occurrence of instability is the way in which the perturbation is actually imposed to the specimen [28,43].
According to the Hill's criterion, positive values of the second-order work d 2 W under any imposed strain increment guarantee the stability of the material response. Null or negative values reveal, instead, a latent instability, not necessarily related to the occurrence of an uncontrolled deformation process. More notably, the evidence of d 2 W B 0 points out the existence of at least one loading program under which the control of the incremental response is lost: instability takes place when the imposed loading program exactly coincides with the one along which the control is lost.
According to Nova [43], the loss of both uniqueness and existence of the incremental response occurs when: where X is the control matrix linking the control variables to the response ones. The condition for the loss of controllability (Eq. (8)) corresponds to the vanishing of d 2 W [43]. Beyond inspecting the stability under the experimentally adopted set of control variables, this theory can be employed to track the susceptibility of the material to specific modes of failure and to identify the conditions for their activation [41]. This strategy, known as latent instability analysis, requires the definition of control matrixes corresponding to sets of control variables different from those actually imposed in a test, and the monitoring of their determinant. This approach allows to inspect the consequences on soil stability of possible changes in statickinematic constrains.
In this study, the above strategy is employed for investigating failure in both hardening and softening regimes for fully and partially saturated soils. For this purpose, the control matrixes X s and X w were defined and monitored in the numerical simulations, the first one being checked in drained/suction-controlled tests, whereas the second in undrained/constant water content tests. X s is the hydro-mechanical stiffness matrix that links the control variables _ r Ã and _ s Ã to the response variables _ e and _ S r : where _ e is the principal strain rate vector, _ r Ã is the principal stress measure conjugate to _ e in the definition of d 2 W, defined as: is the mechanical elastoplastic stiffness matrix defined as: where D e is the elastic stiffness matrix, of =or 0 and og=or 0 are the partial derivatives of the yield surface and plastic potential function with respect to r 0 , H is the hardening modulus: and H c is the critical hardening modulus: D wr is a null vector since the WRC is assumed to be independent of void ratio; D ww is a scalar term resulting from the Van Genuchten expression and D rw describes the mechanical effects induced by changes in S r : where of =op c is the partial derivative of the yield surface with respect to p c . Under drained conditions, X s coincides with the mechanical elasto-plastic stiffness matrix, and the vanishing of its determinant points out the attainment of shear failure.
X w can be derived from Eq. (9) assuming as control variables the total stress increments _ r and the variations in water content _ e w = 1 þ e ð Þ, and as response variables _ e and _ s: By monitoring the determinant of X w , the occurrence of potential instabilities due to, for example, volumetric collapse, static liquefaction or undrained shear failure can be captured. Under fully saturated conditions, X w is given by:

Stability analyses under axisymmetric loading paths
In this section, the material response is analysed under the usual geotechnical testing conditions along axisymmetric loading paths. A series of compression triaxial tests was numerically simulated using the commercial finite element code Abaqus/Standard, where the proposed model was implemented adopting an explicit integration scheme with automatic sub-stepping and error control, modifying the Fortran 90 routines developed by Boldini et al. [7] for its fully saturated version. The ideal tests were carried out on a single eight-node axisymmetric element (denoted as CAX8P in the Abaqus code) keeping the suction constant (corresponding to drained conditions for fully saturated soils) or preventing water outflow (corresponding to undrained conditions for fully saturated soils). The simulations were run controlling the axial strain and radial stress. Such mixed stress-strain control guarantees the stability of the soil response along the entire loading path, preventing the development of any instability phenomenon which, conversely, can occur under stress-controlled test conditions. In the mixed control numerical tests illustrated in the following, potential instability conditions can be captured by monitoring the determinant of the appropriate control matrixes, i.e. X s and X w , depending on the imposed hydraulic boundary conditions.
All the simulations were carried out by the WR2-Unsat model illustrated in Sect. 2. This model relies on a flexible and hierarchical formulation suitable for describing soil instabilities caused by non-associated flow rule, structure degradation induced by the accumulation of plastic strains and the changes in pore saturation. This study focuses on the effects on soil stability of the last two factors, assuming for the sake of simplicity an associated flow rule throughout the whole set of numerical experiments. Details on the role of non-associated conditions on the stability of saturated granular media can be found in Nova [44]. In this section, the effects of structure degradation and partial saturation are first inspected separately and then combined, to mimic more realistically the state of many natural soils.
The set of hydro-mechanical constitutive parameters used in the numerical simulations is summarized in Table 1. The considered material is an ideal fine-grained soil characterized by low plasticity (see [7,50]), whose hydro-mechanical properties were selected to emphasize the instability phenomena triggered by the loss of structure and suction.
The study was conducted by employing a relatively simple version of the WR2-Unsat model, characterized by the classical MCC yield surface (a = 1, m = 2, N = M, where M is the slope of the critical state line in the p 0 and q plane), an associated flow rule (a g = a, m g = m, N g-= N) and possible structure degradation. The adopted WRC follows the typical Van Genuchten expression with a low air-entry suction and a non-zero value for the residual saturation degree (Fig. 1).
A series of 16 compression triaxial tests was numerically performed on fully and partially saturated soil specimens (see Table 2). The set of initial conditions was defined with the main aim of analysing the influence of the over-consolidation ratio (considering lightly and heavily over-consolidated soil samples), degree of saturation, and void ratio on soil stability.
In all the cases, the initial stress state is isotropic and the isotropic yield compression pressure is initially equal to 200 kPa. The over-consolidation ratio R p , defined as the ratio between p c and p 0 , is fixed equal to 1.3 to 8. Only for the case of fully saturated structured soils, an additional value for R p of 3.3 is considered. The initial void ratio is determined consistently with the position of the NCL related to the current S r , for a fixed p c (Fig. 2, procedure described in detail in [53]), while it is assumed constant in the elastic domain, given the low compressibility of the material.
In the first set of numerical analyses (Tests 1 to 4), the stress-strain response of fully saturated soils was inspected in absence of destructuration processes. In the second set (Tests 5 to 10), the parameters n v , g v , n s and g s controlling the debonding processes were activated. In the third set (Tests 11 to 14), partially saturated soil specimens, with no cementation bonds, were brought to failure preventing any water outflow. Finally, in the last set (Tests 15 to 16), both sources of instability were combined by considering the effects of mechanical and hydraulic bonding.
For each test, the stress-strain response was analysed to identify the loss of controllability domain in the (q, p 0 ) plane. In addition, for each set of tests a specific analytical criterion is also proposed, to express the stability conditions in the elasto-plastic regime highlighting the role of the soil properties and initial conditions on the occurrence of possible failure mechanisms.

Fully saturated soil with no interparticle bonding
Drained and undrained triaxial compression tests were simulated on fully saturated soil samples by increasing the axial strain while keeping fixed the total radial stress. The results of the tests are shown in Fig. 3 in terms of loading path (a), stress-strain response (b), and evolution of the determinant of the control matrix X s (c). This latter is not directly associated to the mixed control test conditions adopted, as being the one related to the corresponding loadcontrolled tests. However, monitoring the evolution of the detX s allows to detect the latent instabilities related to a potential change in boundary conditions, i.e. revelling what Hydraulic parameter would have been the stability conditions of the soil element if it were subjected to load-controlled test.
As expected, the lightly over-consolidated soil element exhibits a ductile behaviour along the imposed drained stress path (see Fig. 3a, b). The deviatoric stress invariant increases monotonically up to critical state line (CSL), where shear failure occurs. In Test 2, the soil response is brittle and is characterized by a peak in the q-e a curve and a subsequent drop in the available shear strength. For stress paths attaining the yield surface on the dry side, plastic dilation occurs accompanied by a softening-induced shrink of the elastic domain. Related to this, the peak of the deviatoric stress coincides with the yielding point.
In both cases, the simulation ends upon attaining stationary conditions while the soil response remains stable along the entire loading path. In Test 1, the positive definiteness of the matrix X s guarantees the stability of the soil response irrespectively of the selected control variables (see Fig. 3c). Instead, in Test 2 the negative values of detX s reveal potentially unstable conditions that would give rise to a loss of control upon switching to an axial stress-controlled test condition.
These results suggest that, under the above drained conditions, the control can only be lost in the elasto-plastic regime under stress-controlled test conditions and the related loss of controllability domain simply coincides with the dry side of the yield surface (area coloured in pale grey in Fig. 3a).
Starting from the same initial conditions, two undrained compression triaxial tests were numerically simulated by preventing any water flow. As before, both tests were performed under mixed stress-strain control by monitoring at each calculation step the determinant of X w . Once again, this latter matrix refers to stress-controlled loading conditions, thus highlighting the susceptibility of the material to undrained shear failure under such different constrains. Figure 4a shows the stress paths followed by both soil specimens. In Test 3, upon yielding the imposed kinematic constrain causes the rapid build-up of excess pore water pressure and the related substantial reduction of p 0 . Similarly to what observed under drained conditions, the soil response is ductile and the simulation ends upon attaining stationary conditions (Fig. 4b). In Test 4, the combination of two competitive factors, namely the excess pore water pressure-induced increase in the effective confinement pressure and the negative volumetric plastic strain-induced  softening, results in a mildly brittle response. However, the small reduction in shear strength is, in this case, far less abrupt as compared to that of Test 2 and takes place for stress ratios close to that of the critical state line M. Under undrained conditions, the response of lightly over-consolidated samples would always be stable even under stress-controlled conditions, as confirmed by the positive values of detX w , and failure conditions coincide with the critical state ones (Fig. 4c). In Test 4, the attainment of the deviator peak corresponds to the vanishing of detX w and marks the onset of a latent instability (empty circle in Fig. 4c). In this case, besides identifying the potentially unstable conditions, the controllability theory allows to exactly detect the available maximum undrained shear strength of the material, whose value depends on both soil properties and state variables. The resulting failure domain (coloured area in Fig. 4a) covers only a portion of the stress zone related to softening.

Stability criteria in analytical form
The results shown in the previous section suggest that under drained conditions the loss of controllability domain coincides with the softening region. Therefore, the control can be lost along any loading path attaining the yield surface on the dry side.
Under undrained conditions, only highly over-consolidated soils exhibit potentially unstable response, occurring under stress ratios close to that of critical state. In this case, the identification of the loss controllability domain is not readily apparent since the deviator peak does not coincide with the yielding point.
For the above case, shear failure occurs under axisymmetric loading paths when: where D pq , D qq are the terms of the elasto-plastic stiffness matrix linking the increment in the deviatoric stress invariant _ q to _ e v and, _ e s . Assuming an associated flow rule, Eq. (16) is satisfied for stress ratios g ranging from M to g MAX (see Appendix, case 1): where D e pp is the equivalent elastic bulk modulus, of =op 0 is the partial derivative of the yield surface with respect to p 0 , which depends on the shape of the elastic domain. Equation (17) was deduced neglecting the off-diagonal terms of the elastic stiffness matrix. This assumption leads to simple analytical solutions, whose accuracy was preliminary verified by comparing them to the numerical results.
The resulting failure domain, which never involves the wet side, is characterized by a dimension that strongly depends on the assumed soil properties, most notably, on k * , k, m, and a. The influence of each parameter was analysed using the same set of constitutive parameters employed in the numerical study (Table 1). The failure domain increases in size for decreasing values of k* (Fig. 5a). The influence of the soil compressibility is magnified by the increase of p c . This outcome is not surprising, as the reduction of k*, for increasing values of p c , enhances the strain-softening behaviour triggered by postyielding plastic volumetric deformations. A similar effect on the instability domain is observed for decreasing values of the dimensionless parameter k, controlling the elastic bulk modulus (see Fig. 5b). Under undrained conditions, this latter influences the stress changes deriving from the imposed kinematic constrain. The smaller the value of k, the smaller the increase in mean effective stress and, thus, the more likely is the onset of shear failure.
The impact of the shape of the yield surface on the failure domain was investigated by changing the values of the parameters m and a. The increase of m leads to a distorted-shape of the elastic domain, characterized by different stress ratios for the maximum deviatoric stress (see Fig. 6b). Upon increasing the values of m, the region gathering the unstable stress states rotates in a counterclockwise direction (M increases) and its degree of opening decreases. The latter is a direct consequence of the elongated shape of the yield surface along the critical state line. The influence of this parameter is amplified for increasing values of p c , as already observed in the previous cases (see Fig. 6a).
Finally, a lower value of a leads to an asymmetric shape of the yield surface with the apex moving rightward, resulting in a less curved dry side of the yield surface (Fig. 7b). The corresponding failure domain spans a smaller portion of the stress plane, just above the critical state line, showing the already observed increase in size for increasing values of p c (Fig. 7a).

Fully saturated soil with destructuration
The numerical simulations reported in this section differ from the previous ones, in terms of input data, only for the activation of the volumetric and deviatoric destructuration parameters, as all the remaining material properties and initial state variables were kept unchanged (see Table 2). Two additional numerical analyses were performed on an over-consolidated soil specimen with R p = 3.3 under drained and undrained conditions. In the following, the tests results are illustrated in terms of stress-strain response and corresponding latent instability analyses, which were conducted monitoring the evolution of the control matrix X s or X w , depending on the specific drainage condition. Under drained conditions, the lightly over-consolidated sample exhibits a strain-hardening behaviour similar to the one observed in absence of destructuration (Test 5, Fig. 8a, b): the deviator stress increases monotonically up to the final stationary condition corresponding to critical state. The initial reduction of the elasto-plastic stiffness reflects the progressive breakage of interparticle cementation bonds. The soil response is stable irrespectively of the test control conditions, as indicated by the positive definiteness of matrix X s (Fig. 8c). In Test 6, upon yielding, the elastic domain shrinks, inducing a temporary stress path reversal, followed by a final smooth increase. In fact, despite the compaction regime, the structure degradation causes an initial loss of strength, followed by a positive hardening stage starting when the material is already partially destructured. The potential instability observed in the first loading stages is underlined by the negative values of detX s . The change in sign of this determinant corresponds to the turning point from strain-softening to strain-hardening behaviour (empty circle in Fig. 8c).
The highly over-consolidated sample exhibits a brittle response, as expected along a drained stress path intersecting the yield surface on the dry side. The breakage of the cementation bonds increases the rate of strength loss. At yielding the negative value of the detX s highlights the occurrence of shear failure, thus indicating that by that point it would be impossible to keep running a numerical simulation under stress-controlled conditions. In summary, the above drained tests show that adding the de-structuring process results in a larger loss of controllability domain as compared to the unbonded case. Moreover, the failure domain evolves with plastic strains: the one plotted in Fig. 8a refers to the initial conditions. It worth noticing that the stress states lying on the dry side of the yield surface are always potentially unstable (area coloured by dark grey), whereas those located on the wet side in the region coloured by pale grey can be either stable or potentially unstable, depending on the accumulated plastic strains. Figure 9 shows the results of the undrained compression triaxial tests performed on fully saturated soil samples starting from the same initial conditions as those assumed in the drained tests (Table 2). In Test 8, the undrained stress path denotes a potential instability in the compaction regime at stress ratios lower than M (Fig. 9a). Upon exceeding the deviator peak, the stress path tends asymptotically towards the critical state line, reaching the final stationary conditions once the bond degradation is completed. The latent instability is confirmed by the negative values of detX w in the post-peak branch (Fig. 9c). Both over-consolidated soil samples exhibit brittle responses. In Test 9, the yield point coincides with the maximum deviatoric stress that the soil can sustain. In Test 10, the loss in shear strength takes place after yielding, when the increase in the confining pressure is no longer sufficient to compensate for the softening due to plastic strains (empty circle in Fig. 9a). In both cases, the vanishing of the detX w reflects the initiation of failure conditions. In all the undrained cases, the progressive breakage of cementation bonds leads to significant strength reductions, as taking place while the stress path approaches the critical state line. Similarly to the drained case, even in the undrained one, structure degradation processes enhance soil fragility and result in a larger loss of controllability domain as compared to the unbonded case. For the case under study, characterized by the set of parameters listed in Table 1, the undrained failure domain is confined between g = 1.3 and 1.72 (coloured area in Fig. 9a). However, these threshold values depend on the degradation parameters n v and n s adopted in the numerical analyses.

Stability criteria in analytical form
Under drained conditions, the loss of controllability domain coincides to the softening one. Differently from the case with no bond degradation, the failure domain is not univocally defined since it evolves with the internal variables. Therefore, for structured soils the controllability theory can be a very useful tool to monitor step-by-step the stability, thus identifying the occurrence of potential instabilities.
For active de-structuring processes, failure can occur even in the compaction regime when: Focusing on the initial conditions (e d v ; e p s ¼ 0), the above expression simplifies into: The stress ratio delimiting the failure domain on the wet side ( _ e d v ¼ _ e p v ) can be deduced from the following expression: where d is the dilatancy defined as: d ¼ À _ e p v = _ e p s . Adopting as plastic potential function the MCC ellipse, Eq. (20) can be rewritten as: which identifies the loss of controllability domain upon yielding without providing any information on its evolution for increasing strains. Assuming for n v and n s the values reported in Table 1, the control can be lost along any drained stress paths intersecting the yield surface at g [ 1.06 (as confirmed by Fig. 8a). For values of n v larger than 15.8, the soil response becomes potentially unstable upon yielding, for any drained stress paths. As already discussed for the previous cases, such an instability can only be activated by the change in control variables. Under undrained conditions, the region gathering the stress states at failure can be found solving the following expression (see Appendix, case 2): Equation (22) was obtained neglecting the off-diagonal terms of the elastic stiffness matrix and the effects of plastic strains before attaining the deviator peak. Thanks to these assumptions, the initiation of undrained shear failure can be captured using a practical and sufficiently accurate criterion.
Starting from lightly over-consolidated samples, Eq. (22) is satisfied by stress ratios ranging from g MIN to M, where g MIN is given by: For highly over-consolidated samples, instead, failure occurs under stress ratios in the range M B g B g MAX , where: Based on Eqs. (23) and (24), it can be concluded that the loss of controllability domain covers only the portion of the softening domain included between g MIN and g MAX . The limit values depend on soil properties, most notably on soil compressibility, volumetric stiffness modulus, yield surface shape variables and destructuration parameters. Here, the attention will be focused only on the influence of the last factor since the role played by the others coincides with what outlined in Sect. 4.1.1. Figure 10 shows the domains where material instabilities can take place under undrained conditions for different combinations of n v and n s . The curves refer to an MCC-like model characterized by the constitutive parameters listed in Table 1. The activation of the destructuration processes causes an increase in size of the loss controllability domain. Because of the kinematic constrain on volumetric strains, the curves are scarcely affected by the n v value. By assuming n s = 0 and n v [ 0, the instability domain remains confined within the dry side of the yield surface, and its opening is only slightly larger than the one obtained in absence of inter-granular debonding. By activating both components of destructuration, i.e. n s [ 0 and n v [ 0, the control can be lost even for stress ratios lower than M. Finally, increasing values of p c enhance the debonding-induced strain-softening, leading to larger undrained shear failure domains.

Partially saturated soil with no interparticle bonding
The response of partially saturated specimens to mechanical perturbations was numerically investigated by simulating a series of compression triaxial tests (Tests from 11 to 14 of Table 2). The tests were performed at a constant rate of axial strain keeping unchanged the water content.
Herein, the results of triaxial tests at constant suction are not reported as being similar to those obtained on saturated samples under drained conditions. The tests differ from those reported in the previous sections in terms of initial values of S r , s and e (see Table 2 and Fig. 2). In Test 11, because of the low R p value, the imposed stress path reaches the yield locus on the wet side (Fig. 11a). Despite the contractive volumetric behaviour, the soil element exhibits a brittle response. More in detail, the initial hardening regime is followed first by a sharp reduction of the deviatoric stress, controlled by the increase in S r , and then by a recovery and increase of strength when saturation conditions are finally attained (Fig. 11b). The transition from unsaturated to saturated conditions is accompanied by the rapid growth of the pore water pressure and the corresponding reduction of the effective confining pressure (Fig. 11a, c), as typically observed in undrained triaxial tests. By continuously increasing axial strains, the soil attains the final critical state conditions. The latent instability is analysed by a step-by-step monitoring of the detX w,NORM , defined as:  the maximum absolute value of detX w reached in each test (Fig. 11d). The evolution of detX w,NORM confirms the latent instability in the post-peak branch resulting from the hydraulic debonding induced by the S r increase. At this stage, a change in the control variables (i.e. by imposing increasing axial total stresses) would result in uncontrolled soil deformations and pore saturation. The soil switches back to an unconditionally stable response once fully saturated conditions are approached (detX w,NORM [ 0). The attainment of both the maximum and minimum deviatoric stresses coincides with a vanishing detX w,NORM (empty circles in Fig. 11d). Test 12, carried out on an highly over-consolidated specimen, exhibits a strain-softening response similar to the one observed in drained tests performed on fully saturated highly over-consolidated soil elements (Fig. 11a). The peak of the stress-strain response corresponds to the yield point, while the loss of shear strength only ends upon attaining critical state conditions (Fig. 11b). Because of the dilatant soil behaviour, the deviatoric loading stage induces a continuous decrease of u w , that keeps the soil in a partially saturated state (Fig. 11c). The soil response is potentially unstable under the entire plastic loading stage, preventing the run of any axial load-controlled numerical analysis (detX w,NORM \ 0, Fig. 11d).
Tests 13 and 14 differ from the previous ones only for the initial conditions: keeping constant p c , S r is increased from 0.59 to 0.8 and, consequently, e decreases from 0.8 to 0.63 (the ratio n/S r decreases from 0.75 to 0.48).
The stress-strain response of the lightly over-consolidated specimen is characterized by a local peak upon yielding, a wetting-induced abrupt loss of strength and a final increase of q when the soil is fully saturated (Test 13, Fig. 12a, b). In this case, the combination of the low void ratio and high saturation level makes the hydro-mechanical coupling responsible for the potential instability. The development of the latent instability is confirmed by the negative values for detX w,NORM and is interrupted by the transition from unsaturated to fully saturated conditions (detX w,NORM = 0, empty circle in Fig. 12c, d). Once saturated, indeed, the specimen exhibits a strain-hardening response regardless of the way in which the test is controlled. Test 14 shows a response similar to the corresponding Test 12 (see Fig. 12a): the major difference is in the small increase in the deviatoric stress observed in Test 14 just after yielding (Fig. 12b), which results from the reduction of S r causing a homothetic expansion of the yield surface. Thus, in this case the hydro-mechanical coupling has a beneficial effect on the stress-strain response as it delays the onset of the shear failure. The latter is pointed out by the vanishing of the detX w,NORM (empty circle in Fig. 12d), marking the turning point from an stable to a potentially unstable response.
The numerical results demonstrate that the soil stability is markedly influenced by both the initial values of state variables and their evolution along the imposed stress path. This dependency does not allow to establish a direct correspondence between the stress state and soil stability and, as such, to identify a loss controllability domain in the stress space. In this perspective, some factors, as the initial n/S r value, the volume changes induced by the loading path and the resulting variations in S r , seem to be crucial for the prediction of soil stability. Most notably, in the elastoplastic regime, the response of soil samples with high initial n/S r values is potentially unstable at g [ M, whereas it is stable at g \ M provided that the wetting-induced volume changes keep limited. Starting from low n/S r values, latent instability is predicted at g \ M, whereas, on the dry side, the soil response remains stable until its behaviour is governed by the hydro-mechanical coupling.

Stability criteria in analytical form
The numerical results highlight a remarkable increase in complexity of the stability analyses when switching from saturated to unsaturated conditions. In tests at constant suction, the loss of controllability domain is confined to the dry side of the yield surface, where the soil exhibits a strain-softening response. Under undrained conditions (i.e. keeping fixed the water content), the failure domain is expected to vary during the deformation process because of the dependency of the soil stability on the incremental loading path and the evolving state of the material.
The instability of unsaturated soils under axisymmetric loading paths occurs when: where D qp , D qw are the terms of the elasto-plastic hydromechanical stiffness matrix linking _ q to _ e v , and _ S r , respectively. Equation (25) is satisfied when (see Appendix, case 3): By expressing: Equation (26) can be simplified into the following form: In Eq. (27), os oS r is the slope of the WRC, and os oS r S r þ s controls the variations of _ p 0 induced by the changes in pore saturation. This solution is obtained by neglecting the offdiagonal terms of elastic stiffness matrix, following an approach similar to the one adopted in the previous cases. The numerator in Eq. (28) evolves with both the current stress state and state variables, keeping, in all the investigated cases, a value close to zero. It follows that, under the compaction regime ( _ e v [0, _ S r [0), the soil response remains stable for n/S r larger than the limit value. By contrast, the dilatant behaviour ( _ e v \0, _ S r \0) leads to latent instability when n/S r exceeds the limit value. Therefore, the potential triggering of instability can be captured by comparing n/S r , here used as a compact descriptor of the state of the material, with its limit value that delimits the state under which the soil response is controlled by the hydro-mechanical coupling. It worth noting that the ratio n/S r used above is defined for degrees of saturation ranging over the open interval (0,1), where the endpoints representing dry and fully saturated conditions, respectively, are neglected.
These conclusions are confirmed by the numerical results shown in the previous section (Fig. 13). In Test 11, despite the hydraulic debonding induced by the increase in S r , the soil exhibits an initial strain-hardening response. The latent instability emerges when n/S r becomes lower than (n/S r ) lim (see the empty circle, Fig. 13a). This circumstance triggers the volumetric collapse that causes the temporary reduction in the deviatoric stress. In Test 13, the initial state is below the limit n/S r value. Under such a condition, the hydro-mechanical coupling governs the soil behaviour making the material response potentially unstable just after yielding.
An opposite trend was observed in Tests 12 and 14 performed on highly over-consolidated states (Fig. 13b).
The tests differ only for the initial value for n/S r . The response of the looser specimen is scarcely affected by the hydro-mechanical coupling and, thus, latent instability appears upon yielding (exhibits a strain-hardening response). By contrast, the denser specimen initially benefits from the S r reduction due to the dilatant soil behaviour. The positive effect of the hydraulic bonding is lost when n/S r exceeds the limit value, resulting in a potentially unstable soil response.
Based on these results, it can be concluded that, when n/ S r [ (n/S r ) lim , particular attention should be paid to the stress paths intersecting the yield surface on the dry side. These paths always lead to instability, which can actually take place or remain in latent form depending on the set of control variables. On the wet side, the instability is predicted only under loading paths characterized by significant volume changes. The loss of control is, hence, less likely.
When n/S r \ (n/S r ) lim , the soil response remains stable only along those loading paths intersecting the yield locus on the dry side and not inducing significant increment of volume. In all the other cases, the soil is susceptible to failure, as such it is important to identify the conditions for the activation of the instability phenomena.
In all the cases at hand, (n/S r ) lim is about equal to the product between the soil compressibility and the parameter governing the hydro-mechanical coupling. Small changes occur only when fully saturated conditions are approached in Test 11 and 13 (Fig. 13a).
Under general conditions, (n/S r ) lim depends on material properties (k * , b) and state variables (p 0 , p c , s, S r ), see Eq. (27). The influence of these factors is analysed in Figs. 14 and 15 by assuming an MCC-like yield locus, an associated flow rule, and a Van Genuchten WRC characterized by the hydraulic parameters listed in Table 1. Increasing values of b or k * enhance the hydro-mechanical coupling, leading to higher values of (n/S r ) lim (see Fig. 14a, b). For a fixed value of b or k * , (n/S r ) lim depends on both p c and the ratio p 0 /p c that defines the position of the current stress state on the yield surface. Increasing values of p c imply a decrease of the limit value on the dry side (p 0 / p c \ 0.5), and a slight increase on the wet side (p 0 /p c-[ 0.5). For fixed p c , (n/S r ) lim decreases for increasing p 0 / p c , assuming a constant value on the wet side.
Finally, Fig. 15 shows the influence of pore saturation on (n/S r ) lim . Figure 15a refers to stress states lying on the wet side at p 0 /p c = 0.6. Regardless of the p c value, no significant changes in (n/S r ) lim are observed under pore saturation ranging from 0.65 to 0.85. In this range, the moderate slope of the WRC augments the hydro-mechanical coupling, keeping (n/S r ) lim to the maximum value. By contrast, when S r approaches the residual value or full saturation, the slope of WRC appreciably increases causing significant reductions in the (n/S r ) lim value. On the dry side, the curve follows the opposite trend (see Fig. 15b related to p 0 /p c of 0.4). Indeed, (n/S r ) lim keeps at its minimum value when S r ranges from 0.65 to 0.85, whereas appreciably increases when S r is outside this range.
In conclusion, on the wet side, the instability is more likely in highly collapsing soils, characterized by large values of b and k * and small slopes of the WRC. By contrast, on the dry side, the risks associated with a potential loss of control are higher in hard soils having low b values and WRCs with small slopes.

Partially saturated soil with destructuration
The combined effect of partial saturation and bond degradation on the hydro-mechanical behaviour is here investigated by simulating two compression triaxial tests under constant water content (Tests 15, 16 of Table 2). These simulations differ from Tests 11 and 12 only for the activation of the destructuration processes. Herein, the results of triaxial tests at constant suction are not reported since they are similar to those obtained on bonded fully saturated samples under drained conditions. In Test 15, the deviator peak coincides with the yield point (see Fig. 16a). The brittleness in the soil response results from the combined effect of the strain-induced structure degradation and the homothetic contraction of the yield locus due to the increase in S r . After this initial softening, the deviatoric stress lightly increases again, as the impact of the debonding process reduces for increasing cumulated plastic strains. This trend is inverted once the hydro-mechanical coupling becomes the main factor controlling the soil behaviour (see Fig. 16b): at this stage, the increase in S r leads to a significant reduction in q down to extension triaxial stress conditions. Once the material is fully saturated, the stress path reverses again its direction and q starts increasing, initially along an almost p 0 path, to then deviate to the left towards the critical state line. This trend is due to the sudden build-up of positive excess pore water pressures related to the imposed kinematic constraint (see Fig. 16c). Before reaching the final stationary conditions, a further light decrease of the deviatoric stress is observed: this feature should be ascribed to the residual destructuring process still ongoing at that final stage.
Once again, the mixed stress-strain control guarantees the stability of the soil response. However, under unsaturated conditions, latent instabilities are observed just after yielding and before full saturation. The vanishing of detX w coincides with the attainment of the local minimum and maximum of the deviatoric stress. When fully saturated, the soil response is again potentially unstable once the critical state line is approached, as confirmed by the negative values for detX w .
In Test 16, the material exhibits a brittle response similarly to what observed in the absence of cementation bonds (see Test 12). The activation of the structure degradation makes the loss of shear strength more rapid as compared to that observed for the unstructured case. The dilatant soil behaviour keeps the material in an unsaturated regime and leads to increasing suctions (see Fig. 16c). The hydraulic hardening induced by the reduction of S r is not sufficient to balance the mechanical softening. The soil response is, hence, potentially unstable along the entire stress path (detX w \ 0, Fig. 16d).
These results indicate that the activation of the destructuration process causes an increase in size of the loss of controllability domain as compared to the case with no cementation bonds. Upon yielding, the domain gathering the potentially unstable stress states covers the entire yield surface. However, this outcome should not be generalized, as it depends on the specific test conditions and soil properties assumed in the proposed tests. Moreover, the domain evolves with the material state, as already found for unstructured unsaturated soils.

Stability criteria in analytical form
The results shown in the previous section confirm the ability of the controllability theory to capture potential/ actual instabilities in undrained triaxial tests performed on structured unsaturated soils. Under undrained conditions, just after yielding (e p v , e p s ¼ 0), the instability phenomena take place when (see Appendix, case 4): where in this case (n/S r ) lim,s is defined as follows: and l is a function of soil compressibility and destructuration parameters, and is given by: As in the previous examples, Eq. (29) was deduced by neglecting the off-diagonal terms of the elastic stiffness matrix and can only be employed to check the soil stability upon first yielding (i.e. for zero initial plastic strains). This criterion differs from the one proposed for the unsaturated soils with no interparticle bonds (Eq. (28)) only for the introduction of the factor l. Therefore, most of the comments on soil stability reported in the previous section are also valid in this case and summarized as: (i) control can be lost in stress-controlled tests when n/S r equals (n/S r ) lim,s , causing an uncontrolled growth of soil deformations and pore saturation; (ii) under the compaction regime, the soil response upon yielding is stable when n/S r is larger than the limit value, whereas it is potentially unstable when n/S r-\ (n/S r ) lim,s ; (iii) under the dilatant regime, the soil response at first yielding is stable when n/S r is lower than  Fig. 16 Results of water content-controlled triaxial tests performed on partially saturated bonded specimens the limit value, whereas it is potentially unstable when n/S r is greater than (n/S r ) lim,s .
The limit value for n/S r defined in Eq. (30) is a function of material properties (k * , b, n v , n s , shape parameters of YS and WRC) and state variables (p 0 , p c , s, S r ). In the following, only the influence of the destructuration parameters is investigated, as the impact of the other factors has already been analysed in absence of cementation bonds (see Figs. 14 and 15). Figure 17 shows the loss of controllability domains under undrained conditions for different combination of n v and n s , assuming the same set of constitutive parameters employed in the numerical analyses (Table 1).
In absence of inter-granular bond degradation (n v , n s-= 0), the limit value is roughly equal to the product between b and k* (Fig. 17a). Upon activating the structure degradation process, the limit value increases on the wet side and decreases on the dry one. In both cases, the region involving potentially unstable conditions expands, making more likely the occurrence of instabilities. Assuming a non-null value only for n v , the (n/S r ) lim,s value abruptly changes at p 0 /p c equal to 0.5, corresponding to the critical state conditions, passing from 0.38 on the dry side to 0.73 on the wet side (Fig. 17b). If destructuration process is assumed to depend on both volumetric and deviatoric strain components (n v , n s = 5), the limit value varies remarkably with p 0 /p c (Fig. 17c). On the wet side, (n/S r ) lim,s increases as p 0 /p c decreases and tends to infinite when p 0 /p c is equal to 0.65, corresponding to g = 1.06. Under p 0 /p c ranging from 0.5 to 0.65, the response is potentially unstable regardless of the n/S r value, because of the strain-softening due to the mechanical degradation process. On the dry side, (n/S r ) lim,s increases as p 0 /p c decreases, assuming null value at critical state. According to the numerical results, starting from n/S r of 0.75, the response of both lightly and highly over-consolidated specimens will be potentially unstable once the material yields.

Concluding remarks
In this paper, a single-surface elasto-plastic model for unsaturated cemented soils is presented. This model is formulated extending to unsaturated conditions the WR2 model proposed by Boldini et al. [7] for saturated bonded soils following the approach suggested by Jommi [29]. The latter exploits the modelling techniques typically adopted for bonded soils, considering the capillary forces due to menisci as a further cause of bonding effect on solid grains. The constitutive model was developed with the main purpose of inspecting the mechanical instabilities in natural soils induced by bond degradation due to the accumulation of plastic strains and changes in pore saturation. For this aim, the constitutive equations were combined with the mathematical concept of controllability [15,43], and the response of a single element was analysed under the usual geotechnical testing conditions.
The numerical simulations discussed in the previous sections show that: shape of the yield locus, as highlighted by the proposed stability criterion; • Regardless of drainage conditions, the loss of controllability domain for fully saturated structured soils is larger than the one observed in the unbonded case. Depending on the rate of the degradation process, the failure domain might result as no longer confined to the dry side of the yield locus. Its initial opening under drained and undrained conditions can be determined using the analytical criteria presented in this paper, once the soil properties and material state are known; • Under constant suction conditions, the loss of controllability domain of partially saturated soils with no interparticle bonding is confined within the dry side of the yield locus. Under constant water content conditions, the failure domain cannot be univocally defined since the soil response is strongly dependent on the initial values of the state variables and their evolution along the imposed loading path. Nevertheless, the proposed stability criterion can be fruitfully employed to identify the scenarios where the onset of instability phenomena is more likely. Such criterion highlights the crucial role played by the hydro-mechanical soil properties (k * , b) and the evolving state of the material (e, p 0 , p c , s, S r ) on the occurrence of uncontrolled deformation processes; • Similarly to what observed for fully saturated soils, the progressive breakage of the interparticle bonding in structured unsaturated soils leads to a larger loss of controllability domain compared to the case of unbonded unsaturated soils. As confirmed by the proposed analytical criterion, the structure degradation process makes more likely the occurrence of instability phenomena under generic loading paths.
In conclusion, the results of this study confirm the ability of the proposed constitutive model to capture the basic hydro-mechanical phenomena affecting the behaviour of cemented unsaturated soils. In particular, the model provides a good description of phenomena, i.e. the strain-induced structure degradation and wetting-induced collapse, that, under certain loading conditions, could give rise to instability mechanisms. It is worth recalling that the numerical results reported in this paper refer to single elements subjected to mechanical perturbations. In perspective, the ability of the model in reproducing the response of natural soils is left to be explored under drying and wetting paths and in boundary value problems, where the combined effect of structure and partial saturation might play a role in the response of a whole deposit of soil.