Calibration of a double structure constitutive model for unsaturated compacted soils

This paper describes a calibration procedure for the double structure constitutive model ICDSM (Imperial College Double Structure Model), developed for highly expansive clays, when the model is applied to MX-80 bentonite. Firstly, the model calibration process is discussed and organised in a number of hierarchical steps. These steps involve the estimation of the macrostructural parameters that can be derived from oedometer, isotropic and triaxial laboratory data. Estimation of the microstructural parameters is more challenging due to the limited knowledge of an expansive clay’s fabric and of the physico-chemical phenomena that control its evolution upon wetting. Nevertheless, this paper discusses the available sources of data and identifies the appropriate information that is needed to characterise the micro-structural behaviour of the bentonite. Finally, through the simulation of a swelling pressure test on a bentonite plug, the hydration of the material is studied as a hydro-mechanical coupled process. Particular attention is devoted to the evolution of the stress state of the sample, which is compared to the experimental measurements in order to demonstrate that the constitutive model accurately reproduces the expansive behaviour of MX-80 bentonite.


Introduction
The constitutive model ICDSM has been developed at Imperial College for the modelling of unsaturated, highly expansive clays ( [1]). The complexity of the model requires a hierarchical calibration process. The interpretation of experimental evidence through a double porosity structure framework is a delicate matter requiring knowledge of the micro-structure of compacted clays which is still largely unknown. However, as discussed in [2] for the Barcelona Expansive Model (BExM), the double structure should be accounted for in order to capture the expansive behaviour of compacted clays. In this paper a calibration process, specific to the IC DSM, is proposed for MX-80 bentonite. This process does not follow the same patterns suggested over the years for the Barcelona Basic Model (BBM, [3]), for example in [4]. In particular, the calibration procedure is divided into four different steps: isotropic, oedometric, triaxial and microstructural characterisation. Subsequently, the model is employed to study a radial swelling pressure test carried out in [5], which consisted of constant volume saturation followed by a swelling phase. The numerical results obtained with the ICDSM are compared with the experimental data and the results given by the Imperial College Single Structure Model (ICSSM, [6]).

Overview of the IC DSM
Two stress variables are used in the present framework. One is the equivalent stress σ = σnet+ sair, where σnet is the net stress, defined as σnet = σtot -patm, σtot being the total stress and patm being the atmospheric pressure. sair is the air entry value. The second stress variable is the equivalent suction, seq = s -sair , s being the matric suction.
The model accounts for the double-porosity structure, typical of compacted expansive clays, by introducing two interactive scales of structure: the macrostructure and the microstructure. The fabric is characterised through a model parameter labelled "the void factor", defined as the ratio of microstructural void ratio, em, over the total void ratio, e: For consistency e= em + eM, where eM is the macrostructural void ratio. Thus, the void factor expresses whether the fabric of the material is predominantly influenced by the microstructure or by the macrostructure. The macrostructure is governed by the ICSSM developed in [6]. Conversely, in accordance with experimental evidence, the microstructure is assumed to be saturated and its behaviour is purely volumetric. Consequently, the effective stress principle holds and the material behaviour at this level is controlled by the mean microstructural effective stress, defined as: where p is the mean equivalent stress. Similar to the assumption in the Barcelona Expansive Model (BExM, [7]), the change in microstructural effective stress triggers both reversible and irreversible strains. At the microstructural level the elastic strains are quantified as follows: where the microstructural bulk modulus is defined as: where κm is the elastic compressibility parameter of the microstructure. The material's elastic behaviour upon changes of equivalent stress and/or equivalent suction is defined by contributions from both the macrostructure and the microstructure. The void factor, VF, is used to weight the contribution of each term, hence obtaining the following total elastic strains due to changes in equivalent suction: K M s is the macrostructural bulk modulus with respect to changes in equivalent suction and it is defined as: κs being the elastic compressibility parameter of the macrostructure with respect to changes in equivalent suction. Moreover, the global bulk modulus with respect to changes in equivalent stress is: where K M σ is the macrostructural bulk modulus with respect to changes in equivalent stress, which is defined as: κ being the elastic compressibility parameter of the macrostructure with respect to changes in equivalent stress. Finally, Young's modulus is evaluated as: where μ is Poisson's ratio. With these components the elastic constitutive behaviour is defined. Microstructural plastic strains are assumed to be equal to the elastic microstrains multiplied by a scalar, fβ, which presents the interaction function: The expression for fβ given in [7] is adopted. In the case of micro-compression: or fβ = cc1 if pr/p0<0 (11) and in the case of micro-swelling: or fβ = cs1+ cs2 if pr/p0<0 (12) where cc1, cc2, cc3 and cs1, cs2, cs3 are shape coefficients and the ratio pr/p0 is an expression of the degree of openness of the structure ([1]), as it takes into account the distance of the current stress state from the yield surface. p0 is the isotropic yield stress and pr is related to the current stress state. At every loading step during a finite element analysis, one of the fβ functions from Equations (11) and (12) is active, depending on whether the microstructure is swelling or contracting, which in [7] is established according to the sign of the microstructural effective stress change. However, in a finite element analysis, this selection is problematic because the integration of the constitutive model equations to obtain the change in mean microstructural effective stress, Δp ' , from the changes in total strains, Δε, requires prior knowledge of its value (i.e. Δp ' is an unknown in the iterative solution of the governing equations). Consequently, an alternative criterion from that reported in the literature is adopted here, which relies on the total strain changes when selecting the active interaction function: where Δεvol is the change in volumetric total strain. This adjustment in the formulation allows the implementation of the model in a general finite element software, required to solve boundary value problems. Thus, the model allows the contribution of the microstructure to the irreversible behaviour to be included in a relatively simple conceptual manner. However, this complicates the formulation because the microstructural plastic mechanism cannot be associated with a proper yield surface and, hence, it cannot be defined within the classical plasticity theory. Details on how this issue has been overcome, as well as the model's implementation in a finite element code, are discussed in [1].
Irreversible changes due to loading imply the evolution of the material at both scales. The state of the microstructure is monitored by the void factor, VF, which, in the model formulation, acts as a hardening parameter. Its evolution law is defined as follows: On the other hand, the macrostructure follows the familiar hardening law of the critical state framework: where p * 0 is the equivalent fully saturated yield stress, v is the specific volume, λ(0) is the fully saturated compressibility coefficient and Δε p vol is the change in plastic volumetric strain. The latter is the sum of two contributions, the macrostructural one, indicated with the subscript LC, and the microstructural one, indicated with the subscript β: As a result, the macrostructural hardening ensures the coupling between the two structures. However, while the micro-scale influences the evolution of the macrostructure, the opposite does not apply. Finally, the model takes into account that the fabric of the material undergoes permanent changes upon hydration. This means that the double porosity structure disappears upon full saturation and the material assumes a single porosity structure ( [8]). Such change is irreversible. Consequently, once the suction has reduced below the air-entry value, the void factor is set to zero and the microstructure ceases to exist.

Calibration
Observed soil behaviour combines macrostructural and microstructural contributions, preventing the use of iterative calibration methods, such as the one proposed by [4] for the BBM ( [3]), when calibrating the IC DSM. Limited knowledge of the micro-structure behaviour represents a complication that, at present, cannot be fully resolved. It is important, therefore, to prioritise the calibration of parameters that can be derived directly from experimental data, i.e. those that characterise only one level of structure. Calibration is carried out in four steps: isotropic, oedometric, triaxial and microstructural characterisation, as discussed below for the case of compacted MX-80 bentonite.

Isotropic characterisation
As indicated by Equation (7), the elastic bulk modulus combines macro-and micro-structure contributions. Consequently, κ is no longer directly measurable through the slope of the Unloading Reloading Line (URL) and cannot be dissociated from κm. Without preliminary knowledge of κm, only a suitable range of values can be established for κ. As post-yield compressibility can be estimated with more certainty, it is prioritised.
Three isotropic tests, T1, T2 and T3, by [9], comprising a wetting phase and a loading phase on compacted MX-80 bentonite (Figure 1), are selected from the literature. Each loading phase is approximated by the URL and the VCL (Virgin Compression Line), yielding the slopes reported in Table 1. The ICL slope for T1 represents a direct measurement of the model parameter λ(0), and ICL slopes for T2 and T3 are used to fit the expression λ(seq)=λ(0)[(1-r) e -βSeq +r] ( [6]), as illustrated in Figure 2. The parameters obtained are: λ(0)=0.25, r=0.61, β=0.00007 1⁄kPa and κ=0.08. κ has been chosen within the range of measured URL slopes.
It is noted that, in T2 and T3, the post-yield slope is still changing at the end of the experiment implying that the axial loading stopped prematurely. Consequently, if the capacity of the laboratory equipment allows it, extending compression paths to higher loads would be beneficial for calibration purposes.

Oedometric characterisation
Oedometer data provide information on the compressibility coefficient with respect to changes in suction, κs. Nevertheless, as expressed in Equation (5), the elastic strains due to changes in suction also depend on κm, therefore the value of κs is chosen within a range of realistic values as outlined below. Two oedometer tests by [10] on compacted MX-80 bentonite, V1 and V2 in Figure 3, are selected. Focusing on the initial wetting phase, a range of values from approximately 0.05 to 0.09 can be defined for the elastic compressibility coefficient κs.

Triaxial characterisation
Since the microstructure is assumed to have a purely volumetric behaviour, any information about the triaxial characterisation of the material refers solely to the macrostructure, as characterised in [6]. Hence, in principle, the strength parameters, Mf and Mg, and the shape parameters of the yield and plastic potential surfaces in the (p,J) plane, i.e. αf, μf, αg and μg, can be directly estimated.
Two triaxial tests, one undrained and one drained, conducted on fully saturated samples of normally consolidated MX-80 bentonite by [11] are considered. The stress paths are shown in Figure 4(a) and Figure 4(b), respectively. In both figures the critical state line (CSL) is shown (dashed line) and it has a slope of Mf =0.5. It can be argued that the material seems to soften after reaching a peak deviatoric stress. Such behaviour would be typical of an over-consolidated sample rather than a normally consolidated one. Additional information about sample preparation could perhaps render the interpretation of these tests slightly less ambiguous. As only a limited amount of triaxial testing on unsaturated MX-80 bentonite can be found in the literature it is difficult to characterise the shape of the yield and plastic potential surfaces in the (p,J) plane. Consequently, despite the flexibility of the IC DSM which allows for virtually any shape ( [6]), it was decided to opt for an associated modified Cam-Clay ellipse. Therefore, the shape parameters are set to αf = αg = 0.4 and μf = μg = 0.9 and Mf = Mg is imposed.

Microstructural characterisation
The IC DSM defines three microstructural parameters: VF, κm and fβ.
The void factor, VF, weights the contributions of the micro-and macro-structure to the elastic behaviour, as shown in Equations (5) and (7). In order to determine its initial value, it is believed that a Mercury Intrusion Porosimetry (MIP) investigation would be suitable. MIP shows the distribution function of the diameter of the pores. Consequently by defining an arbitrary separation between micro-pore and macro-pore diameters, the ratio of micro-pores over total pores, i.e. the void factor, can be obtained. Nevertheless, MIP-based investigations are quite scarce and, more importantly, further research and collaboration between experimentalists and numerical analysts is required to define a reliable relationship between the void factor and the pore size distribution function.
The microstructural elastic compressibility parameter, κm, contributes to the elastic response of the material to changes in effective stress and equivalent suction. At present, to the authors' knowledge, there is no laboratory test that is designed to provide useful information for its estimation. In fact, the physico-chemical phenomena that take place within the micro-structure that are responsible for the swelling of single grain particles ( [8]), are still partly unknown. A suitable method for the measurement of κm can therefore be envisaged as more of a chemical issue than a mechanical one. At present, however, the adopted approach assumes simply that κm must be comparable to the macrostructural compressibility coefficients. Fig. 3. Experimental stress paths, [10]. Fig. 4. Estimation of the CSL slope from an undrained triaxial test (a, above) and a drained triaxial test (b, below).
The interaction functions, fβ, control the microstructural contribution to plasticity, according to Equation (10). Given the arbitrary nature of this assumption, it is impossible to associate fβ with a precise physical meaning and, therefore, with a measurement method. It can be envisaged that, upon improvement of the knowledge of the micro-structure, fβ will be conceptually reinterpreted. For the moment, coefficients cc1, cc2, cc3 and cs1, cs2, cs3 are set to values taken from the literature ( [7]): cc1= cs1=-0.1, cc2= cs2=1.1 and cc3= cs3=2.

Additional parameters
Additional parameters include: Poisson ratio, μ=0.4 ([12]), the air-entry value of suction, sair=1 MPa ( [13]) and the characteristic pressure, pc=1 MPa. These are based on engineering judgement as no direct measurements are available. The complete list of IC DSM parameters for MX-80 bentonite is reported in Table 2. From a hydraulic standpoint, the bentonite is characterised with a Van Genuchten ( [14]) retention curve fitted to data provided by [13]. The logarithm of permeability is assumed to vary linearly from the saturated, ksat, to a minimum, kmin, value of permeability, when suction changes from s1 to s2, according to the model by [15]. All parameters are listed in Table 3.

Numerical analysis and discussion
The radial swelling test by [5] is simulated in the bespoke finite element code ICFEP ( [15]) comprises the constant volume saturation of a cylindrical sample of compacted MX-80 bentonite, 46.8mm in diameter and 40mm high, followed by free expansion in the radial direction until a radial strain of 6.5% is reached. Hydration is imposed by providing the lateral boundary with free access to water. Once the swelling pressures stabilise, the sample is trimmed circumferentially, allowing for the outward radial strain of 6.5% upon hydration.
A coupled consolidation, axi-symmetric analysis of the test was performed, employing 54 8-noded elements with displacement degrees of freedom in each node and pore pressure degrees of freedom in the corner nodes. Initially, the soil has axial and radial stresses of 1.896 and 0.621 MPa, respectively and a suction of 48 MPa, corresponding to a water content of 13% under a dry density of 1.6 Mg/cm 3 .
During confined hydration horizontal and vertical displacements are set to zero along the vertical and horizontal boundaries of the mesh, respectively. The nodal suction at the right-hand side vertical boundary is progressively decreased until seq = 0, while the remaining boundaries are impermeable. Trimming of the sample is then simulated by excavating the appropriate elements, while the lateral displacement constrains on the righthand side vertical boundary are lifted. A suction decrease is re-established on the same boundary until the radial swelling reaches the desired value of 6.5%. At this point, the lateral confinement is also re-established under further reduction in suction.
The evolutions of total axial stress and total radial stress are shown in Figures 5 and 6, respectively. Results from analyses with the ICDSM and the ICSSM are shown. The results obtained using the new double-structure model are very encouraging, in spite of the uncertainties highlighted in the calibration process. The ICSSM, on the contrary, under-predicts the stress level. This is because the micro-structure in compacted clays is responsible for a considerable amount of the swelling.
The material seems to exhibit anisotropic behaviour, since the peak measured values of σr (10 MPa) and σa (12 MPa) do not coincide. The anisotropy could be caused by the one-dimensional compaction process operated by [5] or, alternatively, could be a material feature worth investigating experimentally. The model, on the other hand, predicts isotropic behaviour, as expected. Figure 7 shows two predictions for the radial swelling test, each of which is based on different input parameters and, in particular, on the choice of microstructural parameters. The one labelled "Cal#1" is the same as that shown in Figure 6 and was obtained by employing the model parameters in Table 2, while "Cal#2" was obtained by only changing the values of cc= cs=5, m=0.1 (as in [7]) and s=0.091 (which corresponds to the highest measured value). Both predictions show an excellent agreement with the test data and, at present, there is insufficient experimental evidence to choose one calibration over the other. Hence there are difficulties in characterising unequivocally the material.

Conclusions
The IC DSM has proven effective at reproducing the expansive behaviour of MX-80 bentonite during a radial swelling test. The double porosity structure considered in the formulation strongly improves the numerical predictions compared to predictions obtained from the single porosity IC SSM model. Nevertheless, the proposed calibration of the IC DSM has highlighted the scarcity of experimental data for a suitable estimation of the micro-structure's compressibility and plastic contribution. The importance of investigating the material fabric at the micro-scale, using for example the MIP technique, has been emphasised, as this would pave the way for an improvement of the IC DSM's conceptual representation of the micro-structure and its calibration. The work presented in this paper is funded by Wood PLC and Radioactive Waste Management Ltd., UK.