Finite element modelling of unsaturated soils using a modified version of the Barcelona Basic Model

. In this paper, unsaturated soil behaviour is numerically simulated using a modified version of the Barcelona Basic Model. The model has been implemented as a User-Defined Soil Model (UDSM) in the finite element code PLAXIS. Unlike the original formulation based on net stress ([1]), this model is formulated based on the well-known Bishop’s effective stress, with suction as an additional stress variable. Furthermore, a non-associated flow rule and the dependency on Lode ’ s angle proposed by [2] are incorporated. The Soil Water Characteristic Curve (SWCC) follows the Van Genuchten model ([3]). First, the simulation of laboratory tests along different stress paths and levels of suction are compared against experimental data. Then, an unsaturated soil embankment subject to a time-dependent precipitation is simulated by means of a fully coupled flow-deformation analysis. The results show the capability of the model to reproduce the main features of unsaturated soil behaviour both at material point level and in boundary value problems.


Introduction
Classic soil mechanics, based on Terzaghi's effective stress concept, is formulated for fully saturated soils. However, a large part of the soils on the planet surface, as well as those used in earth constructions, is partially saturated. Therefore, in a wide variety of geotechnical engineering applications, such as (shallow) foundations, tunnelling, or earth dams and embankments, unsaturated soil mechanics should be considered to better capture the effect of partial saturation and suction on the mechanical behaviour of soils.
Many constitutive models have been proposed in the past decades for unsaturated soils. As one single effective stress often fails to represent the partially saturated soil behaviour, most of the models in the literature extend the stress space using two separate stress variables. In this way, there are two main categories of models. In the first category a 'net stress and suction' stress space is used, as in the case of the well-known Barcelona Basic Model (BBM) [1], while in the second category an 'effective stress and suction' stress space is employed (e.g., [4][5][6][7][8]). While the former presents an advantage for laboratory calibration, the latter allows a smooth transition from saturated to unsaturated modelling concept and is thereby more efficient for numerical implementation and modelling. Further developments of the models have also included the introduction of Lode's angle dependency [9,10] to improve the modeling of stress triaxiality.
In this context, a modified version of the BBM was implemented as a User-Defined Soil Model (UDSM) in * Corresponding author: anita.laera@bentley.com PLAXIS finite element code to simulate the behaviour of unsaturated soils. This PLAXIS BBM (PBBM) model incorporates the fundamental features of the original BBM, such as the capability to reproduce collapse upon wetting, together with the use of Bishop's effective stress combined with the suction stress variable, which facilitates numerical stability. Furthermore, compared to the original BBM, this modified version also includes Lode's angle dependency in its plastic formulation. The model was implemented using a fully implicit stress integration scheme with adaptive sub-stepping to ensure both numerical accuracy and stability.
In this paper, after describing the constitutive equations, the model is calibrated based on suctioncontrolled laboratory tests on a clayey silt [11] and used to simulate an unsaturated soil embankment subject to a time-dependent precipitation. The results show the capability of the model to reproduce the main features of unsaturated soil behaviour both at material point level and in boundary value problems.

Constitutive equations
In the following, soil mechanics convention is used, i.e. stresses are positive in compression. Vectorial and tensorial quantities are denoted in bold, while scalars are in normal character.

Principal concepts
Soil is modelled as an isotropic elasto-plastic material of which the total strain rate is decomposed into an elastic and a plastic part (Eq. 1):

̇=̇+̇
(1) As aforementioned, deformation is assumed to be governed by two stress variables. The first one is the Bishop's effective stress tensor ′ defined by [12]: and the second one is suction defined by: where σ is the total stress tensor, the second-order identity tensor, pw and pa the pore water and air pressures, respectively, and Seff the effective saturation, defined by: with and the residual and saturated saturation, respectively. The strain rate tensors (elastic and plastic) can be decomposed into a volumetric part ̇ and a deviatoric part ̇: In a similar way, the effective stress tensor can be decomposed into a hydrostatic part ′ and a deviatoric part , such as:

Elastic strains
The volumetric elastic strain rate can be calculated as: while the deviatoric elastic strain rate writes: where 0 is the initial void ratio, is the slope of the swelling line in (ln p', e) space, is the slope of the swelling line in (ln s, e) space, is the atmospheric pressure, and is the shear modulus calculated as: with ′ the bulk modulus and ′ the Poisson's ratio.

Plastic strains
The plastic strain rate is computed based on the standard plasticity theory with a plastic flow rule and a yield function satisfying Kuhn-Tucker's conditions: where is the plastic potential while represents the plastic multiplier. The yield function writes as: : is the Von-Mises shear stress, with the stress ratio at critical state for a fully saturated compressive triaxial state, related to the friction angle at critical state φ according to the following equation: In Eq. (12), the suction-induced tensile strength p't is computed by: while the unsaturated preconsolidation pressure p'c is characterized by: in which ′ 0 represents the preconsolidation at full saturation, considered as a hardening variable of the model. 0 is the slope of the isotropic normal compression line in the (ln p', e) plane at saturated state, while λ(s) refers to the same quantity but at a given value of suction s: In Eq. (15) and (16), pr, r and β are model parameters characterizing the shape of the Loading-Collapse (LC) curve, i.e., the yield surface in the (p', s) plane (see Fig  1). Eq (15) ensures that ′ ≥ ′ 0 , which realistically indicates a positive influence of suction on soil mechanical behaviour.
In Eq. (12) the ( ) function models the dependency on Lode's angle , proposed by [2]: where J2 and J3 are the second and third invariant of the deviatoric stress tensor , respectively. The introduction of ( ) represents an improvement with respect to both the Modified Cam-Clay model and the original BBM. In fact, it accounts for stress triaxiality and thereby avoids a strength overestimation in triaxial extension. It is important to note that in this case the yield function in the deviatoric plane shows a very close approximation of the well-known Matsuoka-Nakai's yield contour [10].
Similar to the original BBM, a non-associated flow rule is adopted: with α a model parameter. Finally, the isotropic hardening rule is classically defined as:

Model implementation
During each time step, elastic suction strain (see Eq. (7)) is subtracted from the total strain increment, and the resulting strain increment is then used to update effective stresses. To ensure both numerical accuracy and stability, a fully implicit stress integration with adaptive substepping is used.

Parameter calibration and simulation of laboratory tests
In this section, the PBBM model is calibrated based on laboratory test results available for Jossigny silt [11], consisting of triaxial and isotropic compression tests at different suction and confining pressures.
Only one triaxial test has been performed at full saturation, specifically at a confining pressure of 100 kPa under drained conditions. Having identified the yield point and the point at critical state (Fig. 2), the preconsolidation pressure p'0 and the stress ratio at critical state M are assumed as equal to 165 kPa and 1.26, respectively. The latter corresponds to a friction angle of 31.5°. Based on the assumed elastic region, we can estimate a Poisson's ratio equal to 0.3 ( [13]) and a slope of the swelling line κ equal to 0.012. As for the plastic behaviour, the evolution of the plastic volumetric strains during the test allows determining the slope λ0 as equal to 0.04, while the coefficient α for the non-associated flow rule has been determined through a trial-and-error procedure and found to be 0.35. The simulated triaxial compression test results show very good agreement with the experimental data in terms of the evolution of both deviatoric stress and volumetric strain (Fig. 3).
To simulate the tests on unsaturated samples, the Soil Water Characteristic Curve (SWCC) parameters are first to be calibrated. The residual degree of saturation is assumed equal to 0.57, as indicated in [13] for Jossigny silt. The measured degree of saturation is transformed into the effective degree of saturation for the tests performed at a constant suction of 200, 400, 800 and 1500 kPa (Table 1). In PLAXIS, the SWCC curve is calibrated based on Van Genuchten's model [3]: requiring, therefore, the calibration of two fitting parameters ga and gn, gc being (1 − )⁄ . In this case, ga and gn are fitted to be 0.2 and 1.5, respectively (Fig.  4, Table 1).

Fig. 3. Comparison of experimental (E) and simulated (S)
TXD results at a confining pressure s3 equal to 100 kPa.

Fig. 4. Numerical interpolation (S) of the measured (E) degrees of saturation (S and Seff).
Based on the SWCC parameters, the stress paths of other triaxial tests performed at a constant suction of 200 kPa with a total confining pressure s3 of 50, 100 and 200 kPa are represented in Fig. 5. The critical state points lie on a line that is parallel to the critical state line determined for the fully saturated case but shifting to the left, supporting the model's assumption on the shear and tensile strength parameters M and ks and allowing to calibrate ks=0.12 in this case.
In addition, the preconsolidation pressure p'c has been determined based on the isotropic compression tests performed at a constant suction of 200, 400, 800 and 1500 kPa. The parameters pr, r and β are calibrated such that the LC curve characterising the model, interpolates the preconsolidation pressure points at different suction values (Fig. 6). A summary of the model parameters is reported in Table 2.   For all the simulated laboratory tests, the initial void index e0 is set equal to the value measured in the laboratory for each test ( Table 3). Note that in the case of the test performed at a suction and confining pressure of 200 kPa, the initial state coincides with the one at 100 kPa and an isotropic compression up to p equal to 200 kPa is simulated prior to the shear phase. In Fig. 7, the results of the tests simulated at a confining pressure of 100 kPa and increasing suction (0, 200 and 400 kPa) are compared to the experimental results, while in Fig. 8   The results show a stiffer and stronger response with increasing suction (Fig. 7a), well simulated by the model, while, in terms of volumetric behaviour, no measurements are available for the tests on the two unsaturated samples at confining pressure of 100 kPa (Fig. 7b).
The increase in confining pressure has a similar effect in terms of stiffness and strength (Fig. 8a). The stress-strain behaviour at the confining pressure of 50 kPa is very well simulated with the PBBM model, as well as the case of s3 equal to 100 kPa. At the highest confining pressure, the model response is stiffer than what is observed experimentally, although similar maximum deviatoric stress seems to be reached in both cases. Larger volumetric strains with increasing confining pressures are also predicted by the model although the magnitude is generally overestimated when s3 is 50 kPa and underestimated for s3 equal to 200 kPa (Fig. 8b).
Finally, isotropic compression tests at constant suction (200, 400, 800 and 1500 kPa) have been simulated using the same set of parameters. The comparison with the experimental results (Fig. 9) shows a realistic volume change trend during the tests which are qualitatively comparable with experimental data.

Simulation of a road embankment construction
In this section, the construction of a road embankment and the effect of the water table change due to a precipitation event are simulated with the model PBBM using the finite element code PLAXIS 2D [14]. Due to symmetry, only half of the embankment geometry is modelled (Fig. 10) The numerical analysis consists of 5 phases. The stresses are initialised in a slightly overconsolidated (OCR = 1.2) homogeneous soil deposit having set K0 equal to 0.5. The deposit is 20 m deep and the water table is at a depth of 8 m from the ground surface. The soil deposit is modelled using the PBBM model, whose parameters have been calibrated in the previous section based on the laboratory tests on Jossigny silt. In addition, the slope κs has been set to 0.08 for the deposit, and zero for the embankment. The coefficients of permeability in both vertical and horizontal directions have been set equal to 0.1 m/day. The unsaturated and saturated unit weights of soil are equal to 18 and 20 kN/m 3 , respectively. After the stress initialisation, the first 2 m of the embankment are built in 5 days, followed by another phase building the remaining 2 m in additional 5 days. The embankment layers are assumed to have an initial degree of saturation of 70%. The two construction phases are simulated by performing a fully coupled flow-deformation analysis, after which the embankment is left to consolidate for further 30 days.
In the last phase, the effect of intense precipitation is simulated (Table 4).  . 11 shows the vertical displacement evolution in a node at the crest (x = 0 m and y = 4 m) and just underneath the embankment and at a distance of 3 m from the axis of symmetry. A large amount of settlement is observed during the construction of the embankment due to self-weight. During the consolidation phase, due to the dissipation of excess pore pressures, the embankment settles a bit further, reaching a maximum displacement in the selected nodes of 0.65 m at the top and 0.14 m at the bottom.
The precipitation induces a change in the distribution of suction in the embankment and the soils. It is observed (Fig.12) that suction is reduced in both two selected nodes due to wetting in the first 17 days of rain. As a consequence, two deformation mechanisms are expected: a contraction (collapse upon wetting) or an elastic swelling. Which effect is dominant depends on the stress path at each point. In this case, we observe only a small contraction followed by a large amount of swelling due to wetting. When the rain stops, suction also stops decreasing. Moreover, at the top node, water inflow from rain stops supplying to the node, but water outflow keeps moving downward due to permeability, leading to a recovery of suction.

Conclusions
This paper presents a modified version of the Barcelona Basic Model implemented in PLAXIS. The main modifications include the use of Bishop's effective stress instead of net stress, and the Lode's angle dependency. These modifications are motivated by both a more realistic prediction of material response, as well as numerical implementation and modelling convenience.
The calibration procedure of the model parameters was demonstrated using laboratory tests available for Jossigny silt. The simulation of the different triaxial and isotropic compression tests with a single set of parameters has shown the model's capability of successfully reproducing the behaviour of both saturated and unsaturated soils at different levels of suction.
The model was applied to simulate a road embankment construction process using fully coupled flow-deformation analysis. The effects of settlement during and after the construction, as well as rainfallinduced deformation, were analysed. The results are consistent with physical observation, demonstrating the capability of the model to simulate the complex behavior of unsaturated soil during soil-atmosphere interaction.