Modelling large deformation problems in unsaturated soils

Applications of unsaturated soil mechanics often involve large deformations and displacements. This is the case of collapse behaviour of low density soils or the unrestrained swelling of expansive clays. Rain induced instability of unsaturated slopes is a further example of large displacements that cause important damages around the world every year. Since standard lagrangian Finite Element methods are not well suited to model large deformations, particle-based methods are under development. This is the case of the Material Point Method (MPM), which offers an interesting alternative. Recently, the MPM has been extended to model unsaturated soil problems, where the soil is understood as a unique medium integrated by three distinct phases (solid, liquid and gas). In this paper, the MPM computational cycle for unsaturated soils is described. In addition, a validation of the 3-phase MPM approach is presented by means of the modelling of a one-dimensional infiltration problem. Finally, the applicability of MPM to solve slope instabilities is presented. The simulation of the whole instability process of an embankment subjected to rain infiltration is analysed in detail.


Introduction
Many applications associated with unsaturated soils involve large deformations in history-dependent constitutive models.This is the case of landslides events and slope instabilities which are often triggered by heavy rainfall.They cause important damage around the world every year.Swelling problems in expansive clays or collapse behaviour of low density materials are other examples of problems involving large deformations.
Since standard lagrangian Finite Element methods are not well suited to model large deformations due to limitations on mesh tangling, particle-based methods are under development.This is the case of the Material Point Method (MPM) [1], which offers an interesting alternative.MPM discretizes the media into a set of lagrangian material points which move attached to the material and carry all the soil properties.The governing equations are solved incrementally at the nodes of a computational grid that remains fixed throughout the calculation.This dual description of the media prevents mesh distortion problems hence re-meshing techniques are not required.
Recently, Yerro et al. [2] extended MPM to model problems in unsaturated soils by means a coupled formulation that considers 3 phases (solid, liquid and gas) within each material point.In this work, the validation of such formulation in an infiltration problem is presented.Finally, the whole instability process (from failure initiation to final stabilisation) of an embankment slope subjected to rain infiltration will be analysed.

Basis of MPM
Material Point Methods (MPM) are rapidly evolving in the geotechnical field.This is due to their capability to analyse large deformation problems in a continuum framework.
Inspired by particle methods, MPM [1] discretizes the media by means of a set of material points.Each point moves attached to the portion of continuum which carries all the material information such as strain, stress, and mass.The main governing equations are solved at the nodes of a computational mesh that covers the whole domain and it is maintained fixed throughout the calculation.Domain integration is very similar to the operational rules used in FEM but integration points are now the material points instead of the familiar Gauss points in FEM.Linear shape functions are typically used in order to transfer information between both spatial discretizations (nodes and material points).
In this manner, most of the expertise learned from FEM can be extended to the MPM with the additional improvement that large deformations can be simulated without the requirement of re-meshing algorithms.
MPM was initially developed for fluid mechanics [3].However, the necessity of solving hydro-mechanical problems, such as undrained loading, consolidation or changes in water table, has prompted several authors to extend the method to solve multiphase problems.Two different approaches can be distinguished to solve the interaction between solid skeleton and fluids within a porous media or even with free water.The first one is a "single-layer" strategy which is based on representing the porous medium as a unique continuum by means of one set of material points that moves attached to the solid skeleton.It has been adopted by several authors to solve geotechnical problems under saturated conditions (twophases) [4][5][6].It has also been extended to unsaturated problems in [2] (three-phases).The information of pore fluids is carried as an internal variable at each material point.This approach is appropriate to simulate seepage problems but it is not suitable to model free liquid water.A second strategy was developed for modelling fully saturated soils [7,8].It is a "multi-layer" approach in which each constituent -grains and water in saturated soils-is described separately by means of different sets of material points.This MPM formulation has the capability of modelling both water within the pores and free water as a unique continuum which allows the simulation of fluid-structure interactions.However, the number of material points is much higher and the approach requires an additional computational cost.

MPM formulation for unsaturated soils
This work is based on the MPM approach presented by Yerro et al. [2].Here the soil, understood as a unique continuum, is a mixture of three distinct phases: solid (s), liquid (l) and gas (g) (Figure 1).In order to facilitate computations, only one component is considered within each phase (for instance, no air in liquid and no water in gas).The main governing equations are the momentum balances of the gas, the liquid and the mixture.These are integrated into the domain and discretised at the nodes of the computational mesh.At the beginning of each time step, information carried by the material points is mapped to the mesh in order to calculate nodal mass, nodal velocities, internal and external forces and dragging terms.Because it is a fully dynamic formulation, it leads to a system of equations in which liquid, gas and solid nodal accelerations are the principal unknowns of the problem.
An additional term can be included in the system in order to reduce dynamic effects and numerical instabilities.This is defined as a damping force proportional to the corresponding unbalanced force (: proportional factor) and opposite to the phase velocity [2].In dynamic problems, where the accelerations have an important role in the course of the calculation, this factor should be very small ( = 0-0.05).
Once the nodal accelerations are calculated, an Euler-Cromer scheme is used to update velocities, displacements and strains at the material points by means of shape functions.
Afterwards, in order to ensure mass conservation of each phase, the mass balance equations for the solid, liquid and gas are imposed in the material points.Considering the solid grans incompressible, the mass balance of solid becomes the material derivative of the porosity n: where t is time and v s is solid velocity of a material point.
Taking into account (1) and assuming that liquid and gas pressures (P l and P g ) are the state variables, liquid and gas mass balances can be written as follows, v l and v g being liquid and solid velocities at a material point, S l the degree of saturation, and ρ l and ρ g liquid and gas densities.Equations ( 2) and ( 3) provide the relationships to find liquid and gas pressures rates (dP l and dP g ).
This MPM approach is formulated in terms of two stress variables: net stress σ and isotropic suction (s).Writing such variables in the following convenient manner (4) and (5), in which σ is the total stress tensor and I is the identity matrix, the general form of a suitable stress-strain relationship can be written incrementally as follows in equation ( 6): where Δε is the strain increment vector.D and h are, respectively, the tangent matrix and a constitutive vector.
Then, the stress and also other soil properties are updated at the material points; for instance, the degree of saturation, intrinsic permeability and porosity.Finally, the material points carry all the updated information and time is updated.

Infiltration problem
The equation that represents the movement of water in unsaturated soils is the Richard's equation [9], which is derived from the water mass balance and the Darcy's equation (quasi-static liquid momentum balance).However, because the nonlinearities of soil hydraulic parameters (for instance, permeability depends on degree of saturation, and degree of saturation depends on fluid pressures), it is not easy to obtain an analytical solution to describe the unsaturated flow.

E-UNSAT 2016
With the aim of validating the hydro-mechanical MPM formulation outlined in the previous section, an analytical solution for the one-dimensional infiltration problem is required.To do that, and according to Alonso and Lloret [10], Richard's equation has been simplified following these assumptions: • vertical liquid flow • deformability of the solid skeleton is neglected • incompressible solid grains • constant permeability • constant liquid density • constant total stress field • linear water retention curve Linearizing the retention curve according to (7), a s being a constant, the mathematical expression which describes the vertical water flow within an unsaturated soil can be written as equation ( 8):   The previous expression is essentially the diffusion equation, where z is the infiltration direction and C i corresponds to k l and γ l being the intrinsic permeability and specific weight of liquid, respectively.
A dimensionless time T can be defined depending on C i as follows, where h is the infiltration length: At this stage, an MPM model has been developed in order to compare numerical and analytical solutions.A soil column 1 m high is considered.The material is linear elastic and all properties are summarised in Table 1.Whereas contours are permeable for gas, liquid can only flow vertically.Initially, the sample is unsaturated with a constant suction value (s 0 =0.5 MPa).Because the bottom is impervious for liquid, the infiltration length is 1 m.Gravity is neglected and the solid skeleton cannot deform.The wetting process starts when suddenly s=0 kPa (saturation condition) is imposed at the upper boundary.This condition is maintained throughout the calculation.

Material property Value
Porosity a s [kPa -1 ] 0.001 Applying such boundary conditions, the analytical solution that comes out from equation ( 8) is equivalent to the one-dimensional consolidation problem in saturated media, the well-known Terzaghi's solution [11].
The analytical solution for such problem is shown in Figure 2, in which the suction profile is presented along the depth z of the sample for different times during the infiltration process.It is clear that suction decreases as the sample gets wet.Different simulations have been carried out in order to analyse the effect of number of material points placed within each element.In Figure 3, the numerical results considering 1 and 4 material points are presented.In both cases, the numerical solution approaches the analytical solution.However, even if the case with 4 points fits perfectly well the analytical solution (Fig. 3a), some numerical oscillations are observed when only 1 material point is considered (Fig. 3b).Note that in this cases no artificial damping is considered (α = 0.0).
On the other hand, two more calculations are presented in Figure 4, considering 4 points in each element and different damping factors: α = 0.05 (Fig. 4a) and α = 0.75 (Fig. 4b).See also Figure 5, in which the suction evolution of a material point located at z = 0.49 m is presented for different damping factors (α = 0, α = 0.05, α = 0.75).
From Figures 4 and 5, it is clear that the inclusion of artificial damping (α) has a direct effect on the infiltration rate.Large values of α lead to slowing the infiltration process.For small values, α = 0 and α = 0.05, numerical results fits the analytical expression (see Fig. 3b and Fig. 4a).
Note that over-damping (α = 0.75) leads to a gross error in the calculated spatial and temporal evolution of suction.A small damping (α = 0.05) is able to stabilize the MPM solution and it reproduces accurately the analytical solution.Finally, two additional calculations have been carried out in order to emphasise the importance of taking adequate assumptions in the mass balance calculations when unsaturated fluxes are modelled.Consider, in particular, the hypothesis of neglecting the spatial variations of water and air masses.Assuming that gradients of liquid and gas masses can be neglected equations ( 2) and ( 3) can be rewritten as follows: The results presented in Figure 6 have been calculated assuming that equations ( 11) and ( 12) are correct.In particular Figure 6a shows the results maintaining the same properties considered in Figure 3b.It is clear that for this analysis, numerical solution moves away from the analytical one, especially for T=0.2 and T=0.7, which means that spatial variations of fluid masses are relevant and cannot be neglected.
One of the constitutive equations that play a significant role in the spatial distribution of the fluids masses is the water retention curve.In this example, it is a linear function that depends on a s .According to equation ( 7), a small value of a s implies that small changes in suction correspond to small change of S l .All parameters used in calculating Figure 6a coincide with values listed in Table 1 (a s = 0.001 kPa −1 in particular), whereas in Figure 6b a s is a smaller value (a s =0.0001 kPa −1 ).Comparing these two plots, it can be concluded that the larger the a s , the larger the committed error.Therefore, advective fluxes due to spatial variations of fluids masses cannot be neglected in materials such as sand, although such terms will be less relevant in clay.
Another parameter that influences the spatial distribution of fluid mass is the porosity.If materials with very different porosities are in contact, the spatial gradients of fluid mass along the contact cannot be neglected.

Rainfall effects on an embankment slope
The slope stability problem presented here is inspired by a real case in which several road embankments of medium height (6-8 m) became unstable and shallow failures developed after a heavy rainfall, immediately after the end of construction.The estimated run-out of the slides was 2-4 m.The embankments were built in summertime and the soil, a low to medium plasticity sandy clay, was compacted dry of optimum.
The numerical model presented here is a 7 m high slope with an angle of 32.5º.
The elastoplastic constitutive model presented in [2]  First terms (c' and φ') are the values for saturated conditions, which in this analysis are considered c'=1kPa and φ'=20º.Second terms in equations ( 13) and ( 14) introduce suction effects and provide an additional strength.It has been accepted that cohesion increases from c' to a maximum value c'+Δc max .P atm is the atmospheric pressure and B is a constant parameter that controls the rate of apparent cohesion.Although changes in friction angle are typically less relevant, it is considered that it has linear dependence with suction depending on A. In the model Δc max =15 kPa, B=0.7 and A=0.1.
The water retention curve considered in modelling is based on field measurements [12].Other material properties are summarized in Table 2.

Material property Value
Porosity The computational mesh is formed by tetrahedrons.The initial distribution of material points is presented in Figure 7a.A damping factor α = 0.05 is adopted.
The lower boundary is fixed and horizontal displacements along vertical contours are prevented.Lateral and bottom contours are impervious for the liquid phase.A constant zero gas pressure in excess of atmospheric pressure is prescribed in all the boundaries (P g =0 kPa).The initial stresses and pore pressures of the slope are in equilibrium with the gravity force and the prescribed suction (s 0 =800 kPa) distributed along the slope surface, which is in contact with the atmosphere (see initial suction distribution in Fig. 7a) The heavy rainfall is modelled by applying a reduction of suction on the ground surface from 800 to 0 kPa during 10 seconds.Afterwards, the saturated boundary condition is maintained constant during the entire simulation.
The embankment response is presented in Figures 7, 8 and 9 in terms of suction, deviatoric strain and total displacements respectively.Figure 7 shows the evolution of suction as a result of the imposed wetting at 3 different times.In Figure 8 the contours of deviatoric plastic strain at times t = 15s and t = 200s are plotted.High shear strains begin to develop at the slope toe soon after the beginning of wetting due to a strength softening.A shear band defining a potential shallow failure surface at an average depth of 1.5 m is already defined at t = 15s (Fig. 8a), although the slope remains stable.A few seconds later, the slope becomes unstable, large deformations are involved and a failure mechanism is well defined.In Figure 9 the total displacements of 4 material points (indicated in Fig. 7a) are presented.Material points located in the central part of the slope (S2, D2) slide down and pass over the material points located in the slope toe (S1) which experiences small displacements.The lowest point, D2, remains motionless because it is located below the shear band.
The final run-out, defined as the distance between the initial toe of the slope and the toe of the final geometry, is 2.5 m.

Conclusions
A three-phase formulation (solid, liquid and gas) of the MPM is described in the paper.Material points are assumed to carry all the necessary information for the three phases The MPM formulation for unsaturated soils presented in [2] has been validated by modelling an infiltration problem.A complete analysis has been carried out and the effect of some numerical aspects (number of material points, artificial damping term) and modelling assumptions was discussed.
In the second part of the paper an embankment slope instability induced by heavy rains has been simulated.An elastoplastic suction dependent Mohr-Coulomb model formulated in terms of two stress fields (net stress and suction) has been used in order to model strength variations due to suction changes.Suction decrease in the slope results in a marked strength softening.The result is a complex motion, which simulates observations in rainfall-induced instabilities.
Other large deformation problems, such as wetting induced collapse or swelling may be analysed by the same approach but they will require the consideration of a different constitutive model.However, the general 3phase MPM algorithm will remain unchanged.

Figure 1 .
Figure 1.Scheme of the MPM discretisation for unsaturated soils.

Figure 2 .
Figure 2. Analytical solution.Suction evolution along depth z.Different simulations have been carried out in order to analyse the effect of number of material points placed within each element.In Figure3, the numerical results considering 1 and 4 material points are presented.In both cases, the numerical solution approaches the analytical solution.However, even if the case with 4 points fits perfectly well the analytical solution (Fig.3a), some numerical oscillations are observed when only 1 material point is considered (Fig.3b).Note that in this cases no artificial damping is considered (α = 0.0).On the other hand, two more calculations are presented in Figure4, considering 4 points in each element and different damping factors: α = 0.05 (Fig.4a) and α = 0.75 (Fig.4b).See also Figure5, in which the suction evolution of a material point located at z = 0.49 m is presented for different damping factors (α = 0, α = 0.05, α = 0.75).From Figures4 and 5, it is clear that the inclusion of artificial damping (α) has a direct effect on the infiltration rate.Large values of α lead to slowing the infiltration process.For small values, α = 0 and α = 0.05, numerical results fits the analytical expression (see Fig.3band Fig.4a).

Figure 5 .
Figure 5. Suction evolution in a material point located at z = 0.49 m.Comparison between analytical solution and numerical results obtained different damping factors and 4 material points per element.

Figure 6 .
Figure 6.Numerical results when advective fluxes due to spatial variations of fluids masses are neglected.Comparison between different a s values: (a) a s =0.001 kPa −1 ; (b) a s =0.0001 kPa −1 .
is considered in this work to model the effect of suction in the stability of the embankment.The shear strength is written according to the Mohr-Coulomb criterion and strength parameters (c and φ) are be written as follows,