Hydro-mechanical modelling of MX-80 bentonite : one dimensional study

As a first step towards modelling the coupled Thermo-Hydro-Mechanical-Chemical (THMC) behaviour of bentonite, the Barcelona Basic Model (BBM) has been implemented into Numerrinfinite element code. This model has beenfully coupled with the single phase flow equation for unsaturated soils which models liquid water transport. Suction obtained from solving the flow equation is used as an input for the BBM model and the volumetric deformations from the mechanical analysis are used to update the pore water pressure field. As an alternative, BBM is used alongside the Kröhn’s model which assumes that bentonite re-saturation is mainly driven by water vapour diffusion. The paper simulates one dimensional infiltration test on MX-80bentonite with both implemented modelsfor water transport and compares the results withthose from a laboratory experiments based on X-ray tomography.The numerical results of the simulations are similar despite taking into account different physical phenomena.


Introduction
The paper initially describes implementation of the wellknown Barcelona Basic Model (BBM) [1]into Numerrin numerical solver [2].The implementation of BBM bases on standard explicit sub-stepping scheme with error control for stress integration [3][4].In addition in the subsequent finite element simulations two hydraulic models are tested:the single phase flow equation [5] and the vapour diffusion model (Kröhn's model) [6].The single phase flow equation is fully coupled withBBM.TheKröhn's model,on the other hand,does not support coupling directly and is used alongside BBM.The paper gives a brief description of the mechanical and hydraulic models before using them in the simulation of a 1D wettingexperimentof MX-80 bentonite under constant volume conditions.

Barcelona Basic Model
The Barcelona Basic Modelbases on the Modified Cam Clay model [7] and can be viewed as its extensionwhich accounts for unsaturated soil behaviour.BBM uses net stress =  tot -u a and suction s = u w -u a as independent variables, where  tot , u w and u a are the total stress, pore water pressure and pore air pressure respectively.In this paper u a is taken as zero.

BBM in elasticity
The rate of total strain ε  in BBM is decomposed into elastic e ε  and plastic part p ε  : During elasticity, the plastic part of the total strain ratevanishes 0  whereD e is the elastic stiffness matrix : with   in the above formulationK and G are the bulk and the shear modulus respectively, e stands for thesoil void ratio, for the soil swelling index, is the Poisson's ratio andp = tr()/3is the isotropic net pressure.The strain rate relatedto suction changeis: where s is the soil swelling index with respect to suction variation and p atm is the atmospheric pressure.

BBM in plasticity
Figure 1 shows a 3D representation of the yield surface ofBBM which equation is: The symbol q stands for the deviatoric stress: where    and  are the principal stresses.Thetrace of the yield surface in the plane q=0 is the Loading Collapse (LC) curve.
The plastic strain direction is determined using a plastic potential function: The factor  which allows forthe recovery oftheJaky'sapproximationof coefficient of at rest pressure K 0 in 1D compression isdefined as: where is the slope of normal consolidation line.BBM assumes that the soil shear strength increases linearly with product of suction and model constant k.The soil preconsolidation pressure p o is also considered to be suction dependent through the equation: with The preconsolidation pressure at full saturation is indicated by * o p .Equations( 10) and (11) introduce p c ,  and r as new BBM parameters.These parameters are used to define the effect of suction on the preconsolidation pressure and the post-yielding stiffness.
During plasticity the net stress rate is calculated as:

σ D ε D ε m ε
Te s where the plastic strain rate p ε  is determined by the flow rule: The plastic multiplier   can be derived based on plasticity theory [4]: , where:

  
Substituting equations ( 5) and ( 13) into equation( 12) yields the net stress rate: Finally, the preconsolidation pressure is updated as: where * oi p is the initial preconsolidation pressure and

BBM implementation into Numerrin code
BBM is implemented into Numerrinand used inthe Finite Element Method simulations.Stress is integrated with a standard explicit sub-stepping scheme with error control [3][4].In such a scheme if plasticity is detected then an automatic sub-stepping of the strain increment is initiated.The final number of sub-steps is dependent on the required accuracy by the user.
In order to validate the implementation, the data provided in literature [1] has been reproduced.For the sake of brevity only two numerical tests are shown:(a) wetting of unsaturated soil under isotropic condition and(b) shearing under constant confining pressure but at different suction levels.Both testsused model parametersgiven in Table 1.
Table 1.BBM parameters as used in the verification tests, after [1].  2. Numerical results versus reference data: isotropic loading followed by wetting.

Isotropic loading followed by wetting
A soil sample with an initial suction of 200 kPa is isotropically loaded up to a confining pressure of 350 kPa.At this confining pressure the soil is wetted till full saturation, which leads to soil collapse.After the full wetting the soil is isotropically loaded in fully saturated conditions up to 600 kPa. Figure 2shows impeccable match between the numerical results and the reference data.

Shear at different suction levels
The model implementation is furthertested by replicating the shear paths (Figure 3).The soil is sheared under constant confining pressures at three different suction levels (s=100, 200 and 300 kPa).The numerical results and the reference data match perfectly.

Single phase flow formulation
For the unsaturated water flow, the water mass conservation equation can be written as: wheren,  w and S r stand for soil porosity, water density and degree of saturation respectively.The average water velocityv = q / nwhere qis the specific dischargegiven by Darcy's law: wherek(h) is the suction dependent hydraulic conductivity and h is the total head being the sum of the pressure head  and the geodetic head z.Equation (18) shows that the net water flux expressed by the left hand side should be equal to the water mass stored or expelled by the soil expressed by the right hand side.On expandingequation (18)it becomes: The flow equation(20)includes mechanical couplingby acknowledging the effects of volumetric deformations v  on the water flow field.It also takes into account the effect of water compressibility through the parameter w .

Implementation into Numerrin code
The Finite Element discretization of equation ( 20)over domain Ωyields: whereN and u are the element shape functions and the nodal displacement respectively, andstands for gradient.
The soil specific moisture capacity Cis: By solving Equation (21) with suitable boundary conditions, one will be able to determine the suction variation overtime.These values of suction are passed to BBM for the deformation analysis.The resultingdisplacement fieldis passed back to the flow equation with the required volumetric deformations to update the flow field and to progress into the next time step in solving the coupled hydro-mechanical system.The basic input data needed to solve the unsaturated flow equation are the soil water retention curve (SWRC) and the soil permeability curve.These curves identify how the soil water content and permeability are changingwith suction variation.In most case some points of the SWRC are experimentally determined and then fitted to one of the common SWRC models.For example van Genuchten [8] or Brooks-Corey [9]modelscan be used for fitting the soil water retention data.The permeability curve is more difficult to establish and is usually based on the SWRC data.

Extended vapour diffusion model
The extended vapour diffusion model or Kröhn's model [6] is based on a phenomenological description of bentonite re-saturation.The model considers that the saturation process is driven by the vapour diffusion in bentonite pores and by the inter-lamellar water diffusion.The fundamental balance equation for this model is: The extended vapour diffusion model implementation has been used alongside BBM.The water content whas been converted to suction using the SWRC which affected the deformations.However, those mechanical deformationshaveno effect on the water transport as indicated inequation (23).

Simulation of 1D wetting test of MX-80 bentonite 5.1 Experimental procedure
One-dimensional wetting and swelling of a compacted cylindrical MX-80-bentonite sample confined in a constant volume were measured using a non-invasive method based on X-ray micro-tomography [10].The experiments were carried out using a table-top X-ray tomographic scanner (SkyScan 1172) and the sample holder shown in Figure 4.The method is based on comparison of X-ray tomographic images of the sample in the original unwetted state and in the wetted and deformed state.Full details of the experimental procedure and data analysis can be found in [11].The initial dry density and water content of the bentonite sample were  d = 1600 kg/m 3 and w = 12.1%, respectively.The duration of the measurement was 7.2 days.With the used sample holder geometry the measurement yields time evolution of the axial distribution of dry density, water content and swelling pressure during the wetting process.

Finite Element model
Figure 5a shows the finite element model used to simulate the wetting test.The model consists of 80 rectangular 4-noded elements with 4 integration points per element.To simulate a constant volume conditions the model is constrained in the directions normal to its boundaries as it clear in Figure 5b.The hydraulic boundaries along both sides of the sample are closed while keeping the top boundary open for discharge.The bottom boundary is kept under fixed zero pressure head to simulate continuous wetting, see Figure 5c.The simulationmonitors the variation of pressure head along the sample height in time as well as the development of the swelling pressure at a control point at the top of the sample (seeFigure 5a).The experiment is simulated twice:the first timethe BBM is coupled with the single phase flow equation(20), the other time with the Kröhn's model(23).

Mechanical parameters
The BBM mechanical parameters have been calibrated using the experimental data [10].The calibration processwas performed in such way that parameters are kept in the ranges of typical values used for MX-80 bentonite [12][13].Table 2 lists the parameters as used in the analyses.During the calibration process it was noticed that the values of the preconsolidation pressure * o p and theswelling index with respect to suction  s have the biggest impact on fitting the swelling pressure data.In fact, the preconsolidation pressure can be estimated depending on the pressing pressure during the preparation of the bentonite sample.However such information is not available for the current study, therefore the preconsolidation pressure was assumed and calibrated in light of data available in literature [14][15].

Hydraulic parameters
The soil water retention curve was estimated using the data from literature [16][17] as well as measurements [18] for the tested bentonite.The approximation of the datawith van Genuchtenmodel is shown in Figure 6.As the initial water content is 12.1% the SWRC gives a corresponding initial suction estimate of about 70 MPa.The saturated hydraulic conductivity was assumed ask sat = 4.0E-14 m/s.Similarly to [13] the soil relative permeability curve is takenas:

Using coupled BBM
Figure 7 shows that the coupled BBM results fit well the experimental measurements for the variation of water pressure head in time.However, the prediction for the evolution of swelling pressure at the control point is less accurate, even though the final swelling pressure value is matched (Figure 10).In the simulation,stresses follow a nonlinear elastic path upon wetting till the suction reaches a value of approximately 10 MPa.Below that suction value plastic deformation occurs and the stress path follows the LC yield curve.

Using decoupled BBM-Kröhn's model
Figure 8shows that the decoupled BBM-Kröhn's model can reproduce the variation of suction during bentonite re-saturation process.The final value of swelling pressure is well captured as shown in Figure 10.

Discussion
The paper presented finite element simulations of a simple wetting experiment.In the calculations BBM model coupled with two different formulations for water transport was used.Both obtained predictionsare very similar no matter which formulation for flow is chosen (Figure 9 and Figure 10).That is puzzling as the single phase flow equation and the extended vapour diffusion equation are based on quitedifferent physical processes.
In fact,bentonite re-saturation is most likely happening due to both water transport and vapour diffusion.However, based on this short study it seems that simply adding these two effects togethermaynot lead to improvement in the simulation results.Deeper links to the physicalprocesses behind wetting of bentonite must be included in the future, which may allow for some judgment on the weight of each transport mode during re-saturation.Furthermore, the mechanical model needs to be improvedin order to better capture the swelling pressure development in time.This might be achieved by incorporating the role of the bentonite micro-structural level intothe constitutive modelling framework.

.
The rate of elastic strain e ε  is decomposed further intothe elastic strain rate due to net stress σ e ε  and the elastic strain rate related to suction s e ,1,1,0,0,0}.Thus, the elastic net stress rate can be calculated as:

Figure 1 .
Figure 1.3D representation of the yield surface in BBM.
23) where v and  d arethe water vapour density and the bentonite dry density respectively.Diffusion is controlled by the vapour diffusion coefficient D v and the interlamellar water diffusion coefficient D w.

Figure 4 .
Figure 4.Schematic illustration (a) and design (b) of the sample holder used in X-ray tomographic experiments of 1D axial wettingof MX-80 bentonite [10].

Figure 5 .
Figure 5.Finite element model for the constant volume test on MX-80 bentonite.
sat and S res are the degree of saturation at full saturation (usually S sat = 1) and at the residual state respectively, calibration gave a value of S res = 0.1.ForKröhn's model the calibration yieldeda value ofD v = 10.0E-6m 2 /s for the vapour diffusion coefficient and D w = 4.0E-10 m 2 /s for the inter-lamellar water diffusion coefficient.

Figure 7 .
Figure 7.Pressure head in time: coupled BBM versus measurements.

Figure 9 .
Figure 9.Coupled BBM and decoupled BBM-Kröhn's model predictions for pressure head variation in time.

Figure 10 .
Figure 10.Coupled BBM and decoupled BBM-Kröhn's model predictions for swelling pressure in time.