Numerical Implementation of the Barcelona Basic Model Based on Return-Mapping Integration

Abstract

This paper implemented the Barcelona basic model (BBM) into the OpenGeoSys (OGS) platform for numerical modeling of the coupled hydro-mechanical (HM) behavior of unsaturated soil. Within the implicit integration approach in the OGS, the integration rule of the BBM was developed first and the integration form of the BBM under a return mapping algorithm was built. The closest point projection method was used for calculating the return mapping directions with the associative flow rule. The numerical simulation results show that the BBM is feasible in fitting the experimental results. The numerical integration algorithm can reflect the elastic–plastic mechanical behavior of materials, and improve the calculation accuracy. The material exhibits obvious elastic–plastic characteristics during numerical simulation and experiment, and the water absorption process can lessen the mixture’s compression stiffness while enhancing its recovery stiffness.

1. Introduction

The hydro-mechanical (HM) performance of unsaturated swelling clay soil is a serious concern for many engineering applications including infrastructure, transportation, and engineering buffers [1,2,3]. Take the high-level radioactive waste (HLW) repository as an instance, the typical swelling material, bentonite–sand mixture, is preferred as the backfill or buffer material to seal the repository with an attempt to close the excavating fracture and isolate nuclide migration by its swelling behavior and low permeability [4,5,6]. Evidence from the laboratory suggests a complex HM coupling process in unsaturated swelling soil [7,8,9], and great effort has been devoted to the study of the constitutive model for unsaturated swelling soil. Among these constitutive models, the Barcelona basic model (BBM) and its expanded formulation, named the Barcelona expansive model (BExM), obtained great success in describing the features of the HM coupling processes in unsaturated swelling soil [10,11,12,13].

Numerical simulations are commonly used to analyze the coupling process of multiple physical fields, as carrying out a long-term experimental investigation is often not possible [14]. Many studies have focused on numerical simulations of the mechanical properties of bentonite-based materials during water absorption and expansion, and have made many achievements. By taking the ideal elastic model, the Drucker–Prager plastic model, and the wetting expansive model, the stress–strain curve of bentonite was described in a simulation using ABAQUS [15]. According to the finite element method, the nonlinear swelling pressure development of bentonite under saturation conditions is predicted [16]. Ballarini et al. [17] simulated and analyzed the thermo-hydro-mechanical (THM) coupling process of buffer materials by OpenGeoSys, and the sensitive parameters were identified, which has guiding significance for the quantitative selection of materials in construction. Teams from four countries coded with ABAQUS, FRACON, THAMES, and ROCMAS to conduct numerical simulation studies on the coupling process of bentonite-based buffer materials in the international cooperation project DECOVALEX II [18,19,20,21].

Since the BBM can describe many typical features of the mechanical behavior of unsaturated soils [22], including the swelling or collapse strain caused by wetting, the BBM is widely used in the simulation of HLW buffer materials in countries such as Europe and Japan [23]. Xu et al. [24] used BExM and implemented in TOUGHREACT-FLAC3D, which is a multi-phase reactive transport and geomechanics simulator, to simulate the coupled thermal-hydraulic-mechanical-chemical (THMC) processes, with a focus on the influence of chemical changes on mechanical process over time. Lee et al. [25] employed the thermo-elasto-plastic version of the BBM to characterize the mechanical behavior for the numerical simulations of full-scale engineered barriers experiment (FEBEX) and obtained approximate results in simulations including temperature, relative humidity, total stress, and displacement.

Although research related to the mechanical constitutive model of the coupling process was widely performed, studies that evaluated the influence of integral algorithms on numerical simulation are limited. The integration method of constitutive equation is directly related to the accuracy and stability of calculation [26,27,28]. Different from the general unsaturated soil constitutive model, the integration algorithm of expansive unsaturated soil is more complicated due to the consideration of plastic characteristics and the effect of plastic deformation on water retention characteristics [29,30]. At present, implicit integration algorithms [31,32] and explicit integration algorithms [33] are commonly used in the field of numerical simulation. TOUGHREACT-FLAC3D uses an explicit integration algorithm, which does not require iteration and has high calculation efficiency. However, it does not strictly meet the consistency conditions in the integration process, which ultimately leads to the stress point possibly deviating from the yield surface, resulting in error accumulation and limiting the algorithm’s accuracy [32,34,35]. In contrast, the implicit integration algorithm corrects plastic stress using a nonlinear equation iterative method, which is widely used due to its high calculation accuracy and efficiency [33,36]. The application of the return mapping algorithm allows the deviation of the initial stress to the yield surface to be predicted elastically, and the implicit integration algorithm can correct the deviated stress plastically and pull the deviated stress back to the yield surface to complete the calculation [37,38,39,40].

In this paper, an implicit algorithm is used to study the numerical calculation of the elastic–plastic mechanical model of expansive unsaturated soil. The elastic–plastic mechanical response of bentonite/quartz sand mixture in the HM coupling process is analyzed by numerical simulation, and the applicability of the theoretical model is verified by experiments. The study of the elastic–plastic mechanical properties of bentonite/quartz sand mixtures in the HM coupling process of the HLW repository can provide theoretical and model support for the study of HM coupling and other multi-physical field coupling processes in the geological disposal of HLW and the implementation of engineering. At the same time, it has a certain guiding significance for the future application of this kind of materials in geological engineering and geotechnical slope engineering.

2. Implicit Integration of the BBM

As a traditional model to characterize the elastic–plastic mechanical behavior of unsaturated soil in geotechnical engineering, the BBM was developed based on the CAM clay model, proposed by Alonso et al. [41]. Based on plasticity theory, the critical state model, and the modified CAM clay model, the BBM is suitable for describing the mechanical behavior of unsaturated soils, expansive clay, low plastic sand, silty clay, and other materials [42,43]. Taking pore suction as a variable, the BBM was used to study the influence of soil suction on soil elastic–plastic mechanical behavior, and can be used to characterize the shear strength and consolidation brought on by hydroscopic expansion, collapsibility, and pore suction change in unsaturated rock and soil [44]. The BBM, like other elasto-plastic constitutive models, divides soil deformation into elastic and plastic deformation, employs three stress variables $(p,q,s)$, and employs yield surfaces $(LC,SI,CLS)$ to divide stress space into plastic and elastic space. Figure 1 presents the three-dimensional representation of the yield surface in the thermo-elasto-plastic BBM three-dimensional yield surface in $p-q-s$ space and $p^{\prime }-q-T$ space, where $p^{\prime }$ is the net mean stress (i.e., total stress minus gas-phase pressure), q is the deviatoric stress (or shear stress), s is suction, and T is temperature [45].

Due to the low permeability of bentonite materials, the seepage–stress coupling time is long during the encapsulation process [46]. The numerical calculation time step is an important factor affecting the accuracy of the calculation results and the time cost. Therefore, the implicit integration algorithm based on the return mapping is used to establish the numerical integration expressions of the BBM, which is embedded in the OpenGeoSys numerical simulation software to reduce the calculation steps and improve the calculation accuracy. In this section, the equations of the BBM were introduced firstly, afterward, the integration form of the BBM was developed according to the approach of Borja and Lee [47]. The consistent tangent modulus is a very important part of the constitutive model integration, and the determination of Tangent Modulus of Stress–Strain Curve is given in Appendix A.

2.1. Equations of the BBM

The stress state in the BBM is described by the net mean stress in the form of

(1)${\sigma }^{\prime }=\sigma -u_a\mathit{I}$

in which the volumetric stress p and deviatoric stress q are defined as

(2)$p=\frac{1}{3}tr\left({\sigma }^{\prime }\right),q=\sqrt{\frac{3}{2}}\parallel \xi \parallel ,$

here, $\xi ={\sigma }^{\prime }-1/3tr\left({\sigma }^{\prime }\right)\mathit{I}$.

The elastic strain of volumetric and shear components are expressed as

(3)$d{ϵ}_v^e=\frac{\kappa }{v}\frac{dp}{p}+\frac{{\kappa }_s}{v}\frac{ds}{\left(s+p_{a}t\right)}$

(4)$d{ϵ}_s^e=\frac{dq}{3G}$

It is notable that the BBM has two yield surfaces in the form of

(5)$f_1\left(p,q,s,p_0^{*}\right)\equiv q^2-M^2\left(p+p_s\right)\left(p_0-p\right)=0$

and

(6)$f_2\left(s,s_0\right)\equiv s-s_0=0$

where

(7)$p_s=k_{s}s$

(8)$\frac{p_0}{p^c}={\left(\frac{p_0^{*}}{p^c}\right)}^{\left[\lambda \left(0\right)-\kappa \right]/\left[\lambda \left(s\right)-\kappa \right]}$

(9)$\lambda \left(s\right)=\lambda \left(0\right)\left[\left(1-r\right)exp\left(-\beta s\right)+r\right]$

The yield surface of the BBM is determined by the hardening parameters $p_0^{*}$ and $s_0$, which depend on the total plastic volumetric strain increment $d{ϵ}_v^p$ in the form of

(10)$\frac{d{p}_0^{*}}{p_0^{*}}=\frac{v}{\lambda \left(0\right)-\kappa }d{ϵ}_v^p$

(11)$\frac{d{s}_0}{s_0+p_{at}}=\frac{v}{{\lambda }_s-{\kappa }_s}d{ϵ}_v^p$

2.2. Equations of Associative Flow Rule

The evolution of plastic strain can be expressed with the associative flow rule in the form of

(12)${\dot{\epsilon }}_{vm}^p=\dot{\varphi }\frac{\partial f_1}{\partial {\sigma }^{\prime }}$

where $\dot{\varphi }$ is a factor. Submitting Equation (2) into Equation (12), we get

(13)$\Delta {ϵ}_v^p=\Delta \varphi tr\left(\frac{\partial f_1}{\partial {\sigma }^{\prime }}\right)=\Delta \varphi \frac{\partial f_1}{\partial p}$

The term of $\partial f_1/\partial {\sigma }^{\prime }$ can be obtained by the chain rule as Borja and Kavazanjian [48]

(14)$\frac{\partial f_1}{\partial {\sigma }^{\prime }}=\frac{\partial f_1}{\partial p}\frac{\partial p}{\partial {\sigma }^{\prime }}+\frac{\partial f_1}{\partial q}\frac{\partial q}{\partial {\sigma }^{\prime }}=\frac{1}{3}\left(\frac{\partial f_1}{\partial p}\right)\mathit{I}+\sqrt{\frac{3}{2}}\left(\frac{\partial f_1}{\partial q}\right)\widehat{\mathit{n}}$

where $\widehat{\mathit{n}}$=$\xi /\parallel \xi \parallel$. The derivation of $f_1$ respecting to p, q, $p_0$, and $p_s$ can be deduced from Equation (5) as

(15)$\frac{\partial f_1}{\partial p}=M^2\left(2p+p_s-p_0\right)$

(16)$\frac{\partial f_1}{\partial q}=2q$

(17)$\frac{\partial f_1}{\partial p_s}=M^2\left(p-p_0\right)$

(18)$\frac{\partial f_1}{\partial p_0}=M^2\left(p+p_s\right)$

2.3. Closest Point Projection of the Constitutive Relation

The returning mapping tensor of stress can be defined as

(19)${\sigma }_{n+1}^k={\sigma }_{n+1}^{tr}-{\mathbb{C}}^e:\delta {\epsilon }^p$

where ${\sigma }_{n+1}^{tr}$ is the trial stress. For the volumetric and deviatoric parts of ${\sigma }_{n+1}^k$, which has a formulation of

(20)$p\equiv p_{n+1}^k=\frac{1}{3}tr\left({\sigma }_{n+1}^k\right)=p_{n+1}^{tr}-K\Delta {ϵ}_v^p$

(21)$q\equiv q_{n+1}^k=\sqrt{\frac{3}{2}}\parallel {\xi }_{n+1}^k\parallel =\sqrt{\frac{3}{2}}\parallel {\xi }_{n+1}^{tr}-2\mu \Delta {\gamma }^p\parallel$

where K and $\mu$ are the elastic moduli defined as

(22)$K=\frac{1+e}{\kappa }p$

(23)$\mu =\frac{3K\left(1-2\upsilon \right)}{2\left(1+\upsilon \right)}$

According to Borja and Lee [47], p and q can be determined by

(24)$p=p_{n+1}^{tr}-K\Delta \varphi \frac{\partial f_1}{\partial p}=p_{n+1}^{tr}-K\Delta \varphi M^2\left(2p+p_s-p_0\right)$

(25)$q=q_{n+1}^{tr}-3\mu \Delta \varphi \frac{\partial f_1}{\partial q}=q_{n+1}^{tr}-6\mu \Delta \varphi q$

The updation of hardening parameters $p_0^{*}$ and $s_0$ is complemented by the hardening equations of Equations (10) and (11). The integration of Equations (10) and (11) over a finite time increment yields

(26)$\begin{array}{cc}\hfill {\left(p_0^{*}\right)}_{n+1}& ={\left(p_0^{*}\right)}_{n}exp({\vartheta }^{*}\Delta {ϵ}_v^p)\hfill \\ & ={\left(p_0^{*}\right)}_{n}exp({\vartheta }^{*}\Delta \varphi \frac{\partial f_1}{\partial p})\hfill \\ & ={\left(p_0^{*}\right)}_{n}exp\{{\vartheta }^{*}\Delta \varphi M^2\left(2p+p_s-p_0\right)\}\hfill \end{array}$

(27)$\begin{array}{cc}\hfill {\left(s_0\right)}_{n+1}+p_{at}& =\left({\left(s_0\right)}_n+p_{at}\right)exp({\vartheta }^s\Delta {ϵ}_v^p)\hfill \\ \hfill & =\left({\left(s_0\right)}_n+p_{at}\right)exp({\vartheta }^s\Delta \varphi \frac{\partial f_1}{\partial p})\hfill \\ \hfill & =\left({\left(s_0\right)}_n+p_{at}\right)exp\{{\vartheta }^s\Delta \varphi M^2\left(2p+p_s-p_0\right)\}\hfill \end{array}$

in which

(28)${\vartheta }^{*}=\frac{v}{\lambda \left(0\right)-\kappa }$

(29)${\vartheta }^s=\frac{v}{{\lambda }_s-{\kappa }_s}$

(30)$\int (\nabla q)wd\Omega =n\int (\nabla qw)d\Omega -n\int \nabla q\nabla wd\Omega$

3. Numerical Implementation

To validate the proposed numerical algorithm, this section builds a numerical model based on compression experiments of bentonite/quartz sand mixture. The elastic–plastic mechanical model of expansive unsaturated soil embedded in the OGS was used in the simulation. The THMC processes in porous media are numerically simulated using the finite element approach, which is based on a flexible numerical framework [49,50]. Using the BBM, the elastic–plastic mechanical behavior of a bentonite/quartz sand mixture during compression was simulated, and the findings of the numerical simulation were compared with the experimental results.

3.1. Experiment Benchmark

According to the compression experiment of MX-80 bentonite/quartz sand mixture completed by Wang et al. [51], the ratio of MX80 bentonite/quartz sand mixture used in the experiment was 70:30, and the sample suction was adjusted by dialysis method. Three standard samples with diameters of 38 mm and heights of 10 mm were prepared in the experiment and named SO-02, SO-03, and SO-04. In order to consider the influence of the gap between the sample and the container on the compression process, sample SO-01 with a diameter of 35.13 mm and a height of 10 mm was also prepared.

The specimens were placed in the dialysis device, as shown in Figure 2, and initially loaded 0.1 MPa pressure to ensure full contact between the specimens and the stress sensor, followed by a steam cycle to allow the specimens to reach a preset suction force. The initial dry density of samples SO-02~04 was 1670 kg/m${}^3$, and steam circulation was used to achieve suction forces of 4.2, 12.6, and 38 MPa. Before compression, the samples were allowed to expand fully until stability under the premise of maintaining the axial load of 0.1 MPa. The initial density of sample SO-01 was 1970 kg/m${}^3$, whose diameter was smaller than that in the dialysis container. After absorbing water and filling voids, the dry density of sample SO-01 was close to the initial dry density of the other three samples of 1670 kg/m${}^3$. Sample size and initial suction parameters in the compression test are shown in Table 1.

After the specimen reached the preset suction, the axial compressive stress was gradually loaded up to 50 MPa by stress control at a loading rate of 0.5 MPa/h, and then we unloaded the stress at the same rate. The change of axial displacement with time was measured automatically by the experimental system during the experiment. Through calculation, the change data of volume strain with time in the process of experimental compression recovery is obtained, and then the change data of void ratio with time can be calculated. The relation between void ratio and volumetric strain is

(31)$K=\frac{1+e}{\kappa }p$

3.2. Software Utilization

OGS is a scientific open-source project developed in the 1980s, which is dedicated to solving the coupling problems of temperature–seepage–stress and other physical fields in HLW packaging, geothermal development, energy storage, and other fields. It has been applied in many fields involving the coupling problems of temperature–seepage–stress chemistry and other physical fields. The integration algorithm of the unsaturated soil elastic–plastic model was used in the numerical simulation, and the software implementation of the BBM with implicit integration has been shown in Figure 3.

According to the experimental process described above, the numerical calculation model adopts an axisymmetric model. As shown in Figure 4 for specimen SO-01, the vertical displacement of the bottom side of the specimen was fixed, and set the bottom and side as the water inflow boundary condition until the model reached the preset suction. Set the water inflow boundary condition at the bottom of the model and fix the lateral displacement of the lateral boundary when the sample expands to the pore closes. The suction of the sample would reach to 0.01 MPa when saturated. When the water absorption of sample SO-02~04 reaches the preset suction, the vertical displacement of the bottom boundary and the lateral displacement of the side boundary of the model are fixed, and the bottom boundary of the model is set as the water inflow boundary, as shown in Figure 5. In the process of compression and recovery, the boundary conditions of samples SO-01~04 are the same. The axial pressure is gradually loaded to 50 MPa on the top boundary of the sample for 100 h and then unloaded the pressure within the same time. The time step for both compression and recovery is 1 hour. In the process of water absorption, the final axial strains of samples SO-01~04 are approximately 18%, 6.8%, 5.4%, and 1.2%, respectively. Therefore, in the compression experiment, the initial porosity of samples SO-01~04 is 0.97, 0.75, 0.73, and 0.69, respectively.

In the compression experiment under suction control, parameters related to the numerical calculation of water absorption and expansion process are shown in Table 2, and parameters related to the numerical calculation of the compression and recovery process are shown in Table 3. In the process of water absorption and expansion, the initial stress value of the model was calculated to reach the fixed suction force. The compression coefficient and recovery coefficient in the process of compression and recovery are obtained by fitting the experimental data, which correspond to the slope of elastic V-ln (P) and the slope of saturated V-ln (P) in Table 3. In order to calculate the void ratio during the compression experiment, the displacement and stress at the node of the symmetric axis on the top boundary of the model during the numerical calculation are output. The model volume change was calculated according to the output displacement, and then the void ratio axial stress curves of the samples SO-01~04 were calculated.

3.3. Result Analysis

Experimental results and numerical simulation results are shown in Figure 6. It can be seen from Figure 6 that the four groups of samples show obvious plastic processes during compression and recovery. The yield stress of samples SO-01~04 increased gradually, indicating that the water absorption process reduced the yield stress of the material. According to the fitted compression parameters during the compression process, the compression parameters of bentonite decreased gradually during the water absorption process. It can be seen from the elastic parameters fitted during the recovery process that the absorption of bentonite reduces the linear elastic modulus of the material.

Figure 7 and Figure 8 show the displacement diagram of samples SO-01 and SO-02 respectively. As can be seen from Figure 7 and Figure 8, the axial displacement of the model in the numerical calculation is small within the compression pressure of 0.5 MPa. When the compression load reaches 10 MPa, the model has significant displacement change and reaches the maximum displacement at the compression load of 50 MPa, and the model is mainly manifested as the deformation of the upper part of the model under compression. Then the displacement of the upper part of the sample decreases gradually with axial pressure decreasing. By comparing Figure 7 and Figure 8, it can be found that the sample with smaller suction has a larger axial displacement, indicating that the water absorption process reduces the stiffness of the bentonite sample.

Figure 9 and Figure 10 reflect the strain cloud diagrams of samples SO-02 and SO-03 under six different compression pressures, respectively. It can be seen from Figure 10 that the strain generated by the model at 0.5 MPa is very small, and when the compression stress reaches 2.5 MPa, obvious strain occurs in the area near the symmetry axis at the upper end of the model. When the compression load reaches 10 MPa, the strain in the upper part of the model near the symmetry axis is large, while the strain in the bottom part of the model is close to 0. In the process of gradually unloading the compression load, the position near the symmetry axis at the bottom of the sample has obvious deformation. When the compressive stress is released to 0.5 MPa, sample SO-03 still maintains a large strain in the center area at the bottom, indicating that the model has a large stiffness under the condition of a large suction force.

4. Results

Based on the implicit integration algorithm of return mapping, the numerical integration expression of the BBM is established, and implement the BBM into the OGS platform for numerical modeling of the coupled HM behavior of unsaturated soil. The numerical basis for studying the elastic–plastic mechanical behavior of bentonite/quartz sand mixture in the coupled process of seepage and stress is provided. The corresponding experimental research is conducted and compared the results with the numerical simulation, the results show that:

• (1). Numerical calculation of void ratio along with the change of axial compression stress results can fit the experimental results well, showing that the numerical calculation described the elastic–plastic mechanics response of bentonite/quartz sand mixture parameters with rationality, and the unsaturated soil elastic–plastic model of the numerical integration algorithm can well reflect the bentonite/elastic–plastic mechanical behavior of the quartz sand mixture.

• (2). In the initial stage of the compression process, the strain distribution of the numerical model is uniform, showing an obvious elastic stage, and then the local strain increases, showing an obvious plastic process.

• (3). Both numerical calculation and experimental results show that the water absorption process can reduce the compression stiffness of bentonite/quartz sand mixture and increase the recovery stiffness.

Author Contributions

Conceptualization, H.Y. and Z.J.; methodology, G.L.; software, G.L.; formal analysis, H.Y.; investigation, G.L. and Z.J.; data curation, G.L.; writing—original draft preparation, H.Y. and Z.J.; writing—review and editing, H.Y. and B.L.; visualization, H.Y.; supervision, H.Y. and Z.W.; funding acquisition, H.Y. and Z.W. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The data supporting this study are included in the paper.

Conflicts of Interest

On behalf of all authors, the corresponding author states that there is no conflict of interest.

Abbreviations

${\sigma }^{\prime }$	effective stress tensor
$u_a$	pore pressure
$\sigma$	total stress tensor
$\xi$	deviatoric stress tensor
p	volumetric stress
$p_{at}$	atmospheric air pressure
q	deviatoric stress
s	suction
$s_0$	maximum suction in experiment
$p_0^{*}$	saturated preconsolidation
$p_0$	unsaturated preconsolidation
$p_c$	reference stress
$\kappa$	elastic stiffness for changes in stress
${\kappa }_s$	elastic stiffness for changes in suction
$\lambda \left(s\right)$	stiffness parameters in stress
${\lambda }_s$	stiffness parameter for changes in suction
${\delta }_{ij}$	Kronecker delta
$\mathit{I}$	indetity
r	parameter defining the maximum stiffness
$\beta$	parameter controlling the increment of stiffness with suction
$\epsilon$	strain tensor
f	plastic potential
v	specific volumn
$\upsilon$	Poisson’s ratio

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Tables

Figure 1. Three-dimensional representation of the yield surface in the thermo-elasto-plastic BBM.

Figure 2. Experimental equipment of suction controlled compression.

Figure 3. Numerical simulation implementation of the model.

Figure 4. Numerical modeling approach of SO-01 in suction controlled compression.

Figure 5. Numerical modeling approach of SO-02–04 in suction controlled compression.

Figure 6. Experimental and numerical void ratio vs. vertical stress of suction controlled compression. (a–d) corresponding to samples SO-01 to 04.

Figure 7. Displacement contour of SO-01 at (a) 0.5 MPa, (b) 2.5 MPa, (c) 10 MPa, and (d) 50 MPa in compression, and (e) 25 MPa and (f) 0.5 MPa in recovery.

Figure 8. Displacement contour of SO-02 at (a) 0.5 MPa, (b) 2.5 MPa, (c) 10 MPa, and (d) 50 MPa in compression, and (e) 25 MPa and (f) 0.5 MPa in recovery.

Figure 9. Vertical strain contour of SO-02 at (a) 0.5 MPa, (b) 2.5 MPa, (c) 10 MPa, and (d) 50 MPa in compression, and (e) 25 MPa and (f) 0.5 MPa in recovery.

Figure 10. Vertical strain contour of SO-03 at (a) 0.5 MPa, (b) 2.5 MPa, (c) 10 MPa, and (d) 50 MPa in compression, and (e) 25 MPa and (f) 0.5 MPa in recovery.

Table 1

Physical parameters and dimensions of sample for suction controlled compression.

Tests	${\mathit{\rho }}_{\mathit{d}\mathbf{I}}$ (Mg/${\mathbf{m}}^3)$	D (mm)	s  (MPa)
SO-01	1.97	35.13	0.0
SO-02	1.67	38.00	4.2
SO-03	1.67	38.00	12.6
SO-04	1.67	38.00	38.0

Table 2

Numerical parameters for hydromechanical modeling of bentonite/sand mixture.

Mechanical Parameters
M (-)	1.07
$P_{os}$ (MPa)	0.66
$P_{ref}$ (MPa)	0.01
$\lambda$ (-)	0.1272
${\lambda }_s$ (-)	0.03
$kappa$	0.081
$kapp{a}_s$	0.008
r (-)	0.75
$\beta$ (MPa${}^{-1}$)	0.25
$e_o$ (-)	0.67
k (-)	0.1
$s_o$ (MPa)	20.2
Seepage Parameters
$n_e$ (-)	0.389
a (-)	1.2
b (-)	6.2
$K_{int}$ (m${}^2$)	$7\times {10}^{-21}$
$\alpha$ (-)	22
A (-)	2
B (-)	4

Table 3

Numerical parameters in suction controlled compression.

Paramter	Unit	SP01	SP02	SP03	SP04
M	-	1.07	1.07	1.07	1.07
$P_{0s}$ (-)	MPa	0.67	1.052	1.052	1.052
$P_c$ (-)	MPa	0.01	0.01	0.01	0.01
${\lambda }_0$	-	0.1352	0.104	0.106	0.12
$Kappa$	-	0.0687	0.034	0.0308	0.024
$Kapp{a}_s$	-	0.03	-	-	-
r	-	0.75	0.75	0.75	0.75
$\beta$	-	$0.25\times {10}^{-6}$	$0.25\times {10}^{-6}$	$0.25\times {10}^{-6}$	$0.25\times {10}^{-6}$
$e_o$	-	0.97	0.75	0.73	0.69
$k_s$	-	0.1	0.1	0.1	0.1
${\sigma }_x$	MPa	−1.52	−1.68	−1.443	−1.22
${\sigma }_y$	MPa	−1.52	−1.68	−1.443	−1.22

Appendix A. Determination of Tangent Modulus of Stress–Strain Curve

In order to calculate the stress tensor through normal stress and deviatoric stress, the stress tensor consists of the superposition of normal stress and deviatoric stress, and the stress tensor can be expressed as (A1)$${\sigma }_{n+1}^k=\frac{1}{3}tr\left({\sigma }_{n+1}^k\right)\mathbf{I}+∥{\xi }_{n+1}^k∥\widehat{\mathbf{n}}=p\mathbf{I}+\sqrt{\frac{2}{3}}q\widehat{\mathbf{n}}$$

In the stress–strain curve, the mechanical response of strain increment can be expressed as (A2)$$\Delta \epsilon ={\epsilon }_{n+1}^k-{\epsilon }_n$$

Then, in the new iterative cycle step, the normal tangent modulus can be obtained from the partial derivative of the new stress tensor with respect to the new strain tensor (A3)$${\mathbf{c}}_{n+1}^k=\frac{\partial {\sigma }_{n+1}^k}{\partial {\epsilon }_{n+1}^k}$$

Substitute Equation (A1) into Equation (A3), we get (A4)$${\mathbf{c}}_{n+1}^k=\frac{\partial {\sigma }_{n+1}^k}{\partial {\epsilon }_{n+1}^k}=\mathbf{I}\otimes \frac{\partial p}{\partial {\epsilon }_{n+1}^k}+\sqrt{\frac{2}{3}}q\frac{\partial \widehat{\mathbf{n}}}{\partial {\epsilon }_{n+1}^k}+\sqrt{\frac{2}{3}}\widehat{\mathbf{n}}\otimes \frac{\partial q}{\partial {\epsilon }_{n+1}^k}$$

Corresponding to the partial variables in the above equation, where (A5)$$\frac{\partial p}{\partial {\epsilon }_{n+1}^k}=\frac{\partial p_{n+1}^{\pi }}{\partial {\epsilon }_{n+1}^k}-K{M}^2\left[\left(2\frac{\partial p}{\partial {\epsilon }_{n+1}^k}-\frac{\partial p_0}{\partial {\epsilon }_{n+1}^k}\right)\Delta \varphi +\left(2p+p_s-p_0\right)\frac{\partial \Delta \varphi }{\partial {\epsilon }_{n+1}^k}\right]$$ (A6)$$\frac{\partial p_0}{\partial {\epsilon }_{n+1}^k}=p_0\vartheta \left(s\right)M^2\left[\left(2\frac{\partial p}{\partial {\epsilon }_{n+1}^k}-\frac{\partial p_0}{\partial {\epsilon }_{n+1}^k}\right)\Delta \varphi +\left(2P+P_s-P_0\right)\frac{\partial \Delta \varphi }{\partial {\epsilon }_{n+1}^k}\right]$$ (A7)$$\frac{\partial q}{\partial {\epsilon }_{n+1}^k}=\frac{\partial q_{n+1}^{\pi }}{\partial {\epsilon }_{n+1}^k}-6\mu \Delta \varphi \frac{\partial q}{\partial {\epsilon }_{n+1}^k}-6\mu q\frac{\partial \Delta \varphi }{\partial {\epsilon }_{n+1}^k}$$

In Equations (A5) and (A6), let (A8)$$B=\left(2\frac{\partial p}{\partial {\epsilon }_{n+1}^k}-\frac{\partial p_0}{\partial {\epsilon }_{n+1}^k}\right)\Delta \varphi +\left(2P+P_s-P_0\right)\frac{\partial \Delta \varphi }{\partial {\epsilon }_{n+1}^k}$$

Combined with (A6) and (A8), we can see (A9)$$B=\frac{1}{P_0\vartheta \left(s\right)M^2}\frac{\partial p_0}{\partial {\epsilon }_{n+1}^k}$$

By substituting Equation (A9) into Equation (A5), we can see (A10)$$\frac{\partial p_0}{\partial {\epsilon }_{n+1}^k}=\frac{p_0\vartheta \left(s\right)}{K}\left(\frac{\partial p_{n+1}^n}{\partial {\epsilon }_{n+1}^k}-\frac{\partial p}{\partial {\epsilon }_{n+1}^k}\right)$$

By substituting Equation (A10) into Equation (A5), we can see that (A11)$$\frac{\partial p}{\partial {\epsilon }_{n+1}^k}=\frac{\partial p_{n+1}^{\pi }}{\partial {\epsilon }_{n+1}^k}-K{M}^2\Delta \varphi \left[2\frac{\partial p}{\partial {\epsilon }_{n+1}^k}-\frac{p_0\vartheta \left(s\right)}{K}\left(\frac{\partial p_{n+1}^n}{\partial {\epsilon }_{n+1}^k}-\frac{\partial p}{\partial {\epsilon }_{n+1}^k}\right)\right]$$ (A12)$$-K{M}^2\left(2p+p_s-p_0\right)\frac{\partial \Delta \varphi }{\partial {\epsilon }_{n+1}^k}$$

The above formula can be obtained (A13)$$\frac{\partial p}{\partial {\epsilon }_{n+1}^k}=a_{1}K\mathbf{I}+a_{2}K\frac{\partial \Delta \varphi }{\partial {\epsilon }_{n+1}^k}$$

In the above equation (A14)$$\left\{\begin{array}{c}{a}_1=\frac{1+p_0\vartheta \left(s\right)M^2\Delta \varphi }{a}\hfill \\ a_2=\frac{-M^2\left(2P+P_s-P_0\right)}{a}\hfill \\ a=1+2K{M}^2\Delta \varphi +p_0\vartheta \left(s\right)M^2\Delta \varphi \hfill \end{array}\right.$$

By combining Equations (A5) and (A8), we can see (A15)$$\frac{\partial p}{\partial {\epsilon }_{n+1}^k}=\frac{\partial p_{n+1}^j}{\partial {\epsilon }_{n+1}^k}-\frac{K\partial p_0}{p_0\vartheta \left(s\right)\partial {\epsilon }_{n+1}^k}$$

Substitute Equation (A14) into Equation (A6) (A16)$$\frac{\partial p_0}{\partial {\epsilon }_{n+1}^k}=p_0\vartheta \left(s\right)M^2\left[\left(2\frac{\partial p_{n+1}^n}{\partial {\epsilon }_{n+1}^k}-\frac{2K\partial p_0}{p_0\vartheta \left(s\right)\partial {\epsilon }_{n+1}^k}-\frac{\partial p_0}{\partial {\epsilon }_{n+1}^k}\right)\Delta \varphi +\left(2P+P_s-P_0\right)\frac{\partial \Delta \varphi }{\partial {\epsilon }_{n+1}^k}\right]$$

After finishing, we can obtain (A17)$$\frac{\partial p_0}{\partial {\epsilon }_{n+1}^k}=a_{3}KI+a_{4}K\frac{\partial \Delta \varphi }{\partial {\epsilon }_{n+1}^k}$$

In the above equation (A18)$$\left\{\begin{array}{c}{a}_3=\frac{2p_0\vartheta \left(s\right)M^2\Delta \varphi }{a}\hfill \\ a_4=\frac{p_0\vartheta \left(s\right)M^2\left(2P+P_s-P_0\right)}{a}\hfill \end{array}\right.$$

By sorting out Equation (A7), we can get (A19)$$\frac{\partial q}{\partial {\epsilon }_{n+1}^k}=2a_5\mu \widehat{\mathbf{n}}+2a_6\mu \frac{\partial \Delta \varphi }{\partial {\epsilon }_{n+1}^k}$$

In the above equation (A20)$$\left\{\begin{array}{c}{a}_5=\sqrt{\frac{3}{2}}{(1+6\mu \Delta \varphi )}^{-1}\hfill \\ a_6=\frac{-3q}{1+6\mu \Delta \varphi }\hfill \end{array}\right.$$

According to the associated flow rule and the chain rule (A21)$$F_1(\Delta \varphi )=F_1\left(\Delta \varphi \left({\mathit{z}}_{n+1}^k\right)\right)=q^2-M^2\left(p+p_s\right)\left(p_0-\mathrm{p}\right)\equiv 0$$

Thus we could get (A22)$$\frac{\partial F_1}{\partial {\epsilon }_{n+1}^k}=\frac{\partial F_1}{\partial p}\frac{\partial p}{\partial {\epsilon }_{n+1}^k}+\frac{\partial F_1}{\partial q}\frac{\partial q}{\partial {\epsilon }_{n+1}^k}+\frac{\partial F_1}{\partial p_0}\frac{\partial p_0}{\partial {\epsilon }_{n+1}^k}\equiv 0$$

The partial derivatives of the yield function $F_1$ with respect to $p,q,p_0$, and $p_s$ are (A23)$$\left\{\begin{array}{c}\frac{\partial F_1}{\partial P}=M^2\left(2p+p_s-p_0\right)\hfill \\ \frac{\partial F_1}{\partial q}=2q\hfill \\ \frac{\partial F_1}{\partial P_s}=M^2\left(p-p_0\right)\hfill \\ \frac{\partial F_1}{\partial P_0}=-M^2\left(p+p_s\right)\hfill \end{array}\right.$$

Substitute Equation (A23), (A12), (A16), and (A18) into Equation (A21) and calculate (A24)$$\frac{\partial \Delta \varphi }{\partial {\epsilon }_{n+1}^k}=b_1\mathbf{I}+b_2\widehat{\mathbf{n}}$$ where (A25)$$\left\{\begin{array}{c}{b}_1=\frac{-M^{2}K\left[a_1\left(2p+p_s-p_0\right)-a_3\left(p+p_s\right)\right]}{b}\hfill \\ b_2=\frac{-4a_5\mu q}{b}\hfill \\ b=4a_6\mu q+M^{2}K\left[p\left(2a_2-a_4\right)+p_s\left(a_2-a_4\right)-a_2{p}_0\right]\hfill \end{array}\right.$$

Return (A22) to Equations (A12), (A16), and (A18) to get (A26)$$\left\{\begin{array}{c}\frac{\partial p}{\partial {\epsilon }_{n+1}^k}=K\left(a_1+a_2{b}_1\right)\mathbf{I}+K\left(a_2{b}_2\right)\widehat{\mathbf{n}}\hfill \\ \frac{\partial p_0}{\partial {\epsilon }_{n+1}^k}=K\left(a_3+a_4{b}_1\right)\mathbf{I}+K\left(a_4{b}_2\right)\widehat{\mathbf{n}}\hfill \\ \frac{\partial q}{\partial {\epsilon }_{n+1}^k}=2\mu \left(a_6{b}_1\right)\mathbf{I}+2\mu \left(a_5+a_6{b}_2\right)\widehat{\mathbf{n}}\hfill \end{array}\right.$$

Finally, Equation (A25) is substituted into Equation (A4) to get (A27)$$\begin{array}{cc}\hfill {\mathbf{c}}_{n+1}^k=2\mu \xi \mathbf{I}& +\left[K\left(a_1+a_2{b}_1\right)-\frac{1}{3}2\mu \xi \right]\mathbf{I}\otimes \mathbf{I}\hfill \\ \hfill +& K{a}_2{b}_2\mathbf{I}\otimes \widehat{\mathbf{n}}+2\mu \sqrt{\frac{2}{3}}\left(a_6{b}_1\right)\widehat{\mathbf{n}}\otimes \mathbf{I}\hfill \\ \hfill & +2\mu \left[\sqrt{\frac{2}{3}}\left(a_5+a_6{b}_2\right)-\xi \right]\widehat{\mathbf{n}}\otimes \widehat{\mathbf{n}}\hfill \end{array}$$ (A28)$$\begin{array}{cc}\hfill \xi & =∥{\xi }_{n+1}^k∥/∥{\xi }_{n+1}^{*}∥=q/q^{tr}\hfill \end{array}$$

References

1. Siemens, G.; Blatz, J.A. Evaluation of the influence of boundary confinement on the behaviour of unsaturated swelling clay soils. Can. Geotech. J.; 2009; 46, pp. 339-356. [DOI: https://dx.doi.org/10.1139/T08-125]

2. Fattah, M.; Salim, N.; Irshayyid, E. Swelling behavior of unsaturated expansive soil. Transp. Infrastruct. Geotechnol.; 2021; 8, pp. 37-58. [DOI: https://dx.doi.org/10.1007/s40515-020-00112-z]

3. Zhang, J.; Sun, D.; Yu, H.; Jiang, J.; Xu, Y. Swelling of unsaturated GMZ07 bentonite at different temperatures. Bull. Eng. Geol. Environ.; 2020; 79, pp. 959-969. [DOI: https://dx.doi.org/10.1007/s10064-019-01595-y]

4. Komine, H. Predicting hydraulic conductivity of sand–bentonite mixture backfill before and after swelling deformation for underground disposal of radioactive wastes. Eng. Geol.; 2010; 114, pp. 123-134. [DOI: https://dx.doi.org/10.1016/j.enggeo.2010.04.009]

5. Cui, Y.J. On the hydro-mechanical behaviour of MX80 bentonite-based materials. J. Rock Mech. Geotech. Eng.; 2017; 9, pp. 565-574. [DOI: https://dx.doi.org/10.1016/j.jrmge.2016.09.003]

6. Nazir, M.; Kawamoto, K.; Sakaki, T. Properties of granulated bentonite mixtures for radioactive waste disposal: A review. GEOMATE J.; 2021; 20, pp. 132-145. [DOI: https://dx.doi.org/10.21660/2021.81.GX254]

7. Su, W.; Wang, Q.; Ye, W.; Deng, Y.; Chen, Y. Swelling pressure of compacted MX80 bentonite/sand mixture prepared by different methods. Soils Found.; 2021; 61, pp. 1142-1150. [DOI: https://dx.doi.org/10.1016/j.sandf.2021.06.005]

8. Agus, S.S.; Schanz, T.; Fredlund, D.G. Measurements of suction versus water content for bentonite–sand mixtures. Can. Geotech. J.; 2010; 47, pp. 583-594. [DOI: https://dx.doi.org/10.1139/T09-120]

9. Chen, Z.G.; Tang, C.S.; Ye, W.M.; Wang, D.Y.; Shi, B. Volume change characteristics of bentonite-sand mixture under hydro-mechanical coupling condition. Yantu Lixue/Rock Soil Mech.; 2017; 38, 1041–1051, 1059

10. Wang, Q.; Tang, A.M.; Cui, Y.J.; Barnichon, J.D.; Ye, W.M. Investigation of the hydro-mechanical behaviour of compacted bentonite/sand mixture based on the BExM model. Comput. Geotech.; 2013; 54, pp. 46-52. [DOI: https://dx.doi.org/10.1016/j.compgeo.2013.05.011]

11. Yamamoto, S.; Komine, H.; Kato, S. Development and validation of mechanical model for saturated/unsaturated bentonite buffer. Comput. Geotech.; 2010; 54, pp. 46-52.

12. Rawat, A.; Lang, L.; Baille, W.; Dieudonne, A.C.; Collin, F. Coupled hydro-mechanical analysis of expansive soils: Parametric identification and calibration. J. Rock Mech. Geotech. Eng.; 2020; 12, pp. 620-629. [DOI: https://dx.doi.org/10.1016/j.jrmge.2019.12.013]

13. Lee, J.; Yoon, S.; Kim, G. Mechanical Constitutive Models of Unsaturated Expansive Clays: A Review of BExM. Proceedings of the Korean Radioactive Waste Society Conference; Korea Institute of Science and Technology Information: Daejeon, Republic of Korea, 2018; pp. 231-232.

14. Wang, W.; Kosakowski, G.; Kolditz, O. A parallel finite element scheme for thermo-hydro-mechanical (THM) coupled problems in porous media. Comput. Geosci.; 2009; 35, pp. 1631-1641. [DOI: https://dx.doi.org/10.1016/j.cageo.2008.07.007]

15. Cai, G.Q.; Tian, H.; Li, J. Numerical simulation on water absorption by MX-80 bentonite for nuclear waste disposal. Appl. Mech. Mater. Trans. Tech. Publ.; 2014; 490, pp. 78-82. [DOI: https://dx.doi.org/10.4028/www.scientific.net/AMM.490-491.78]

16. Bosch, J.A.; Ferrari, A.; Laloui, L. Coupled hydro-mechanical analysis of compacted bentonite behaviour during hydration. Comput. Geotech.; 2021; 140, 104447. [DOI: https://dx.doi.org/10.1016/j.compgeo.2021.104447]

17. Ballarini, E.; Graupner, B.; Bauer, S. Thermal–hydraulic–mechanical behavior of bentonite and sand-bentonite materials as seal for a nuclear waste repository: Numerical simulation of column experiments. Appl. Clay Sci.; 2017; 135, pp. 289-299. [DOI: https://dx.doi.org/10.1016/j.clay.2016.10.007]

18. Borgesson, L. DECOVALEX II-Test Case 2C-1, THM Experiment of a Simplified Deposition Hole in Kamaishi Mine–thermo-hydro-mechanical Modeling of Yje Buffer Material Progress Report, Clay Technology AB and FEM-Tech AB. 1998.

19. Nguyen, T.S. DECOVALEX II T-H-M in-situ experiment at the Kamaishi Mine, Japan. Report on Task 2C: Prediction of T-H-M Response of Bentonite and Rock Mass, Technical Report Atomic Energy Control Board, Canada. 1999.

20. Sugita, Y.; Chijimatsu, M.; Fujita, T.; Duc, P.T. Coupled Thermo-hydro-mechanical Experiment at Kamaishi Mine Technical note 14-99-01. Verification of the Buffer Material Emplacement Technique; Japan Nuclear Cycle Development Instition: Tokai, Japan, 1999.

21. Rutqvist, J.; Börgesson, L.; Chijimatsu, M.; Nguyen, T.S.; Tsang, C.F. Coupled thermo-hydro-mechanical analysis of a heater test in fractured rock and bentonite at Kamaishi Mine—Comparison of field results to predictions of four finite element codes. Int. J. Rock Mech. Min. Sci.; 2001; 38, pp. 129-142. [DOI: https://dx.doi.org/10.1016/S1365-1609(00)00069-1]

22. Gens, A.; Alonso, E. A framework for the behaviour of unsaturated expansive clays. Can. Geotech. J.; 1992; 29, pp. 1013-1032. [DOI: https://dx.doi.org/10.1139/t92-120]

23. Rutqvist, J.; Zheng, L.; Chen, F.; Liu, H.H.; Birkholzer, J. Modeling of coupled thenno hydro mechanical processes with links to geochemistry associated with bentonite backfilled repository tunnels in clay fonnations. Rock Mechan. Rock Eng.; 2014; 47, pp. 167-186. [DOI: https://dx.doi.org/10.1007/s00603-013-0375-x]

24. Xu, H.; Zheng, L.; Rutqvist, J.; Birkholzer, J. Numerical study of the chemo-mechanical behavior of FEBEX bentonite in nuclear waste disposal based on the Barcelona expansive model. Comput. Geotech.; 2021; 132, 103968. [DOI: https://dx.doi.org/10.1016/j.compgeo.2020.103968]

25. Lee, C.; Lee, J.; Kim, G.Y. Numerical analysis of coupled hydro-mechanical and thermo-hydro-mechanical behaviour in buffer materials at a geological repository for nuclear waste: Simulation of EB experiment at Mont Terri URL and FEBEX at Grimsel test site using Barcelona basic model. Int. J. Rock Mech. Min. Sci.; 2021; 139, 104663. [DOI: https://dx.doi.org/10.1016/j.ijrmms.2021.104663]

26. Sheng, D.; Sloan, S.W.; Gens, A. Finite element formulation and algorithms for unsaturated soils. Part I: Theory. Int. J. Numer. Anal. Methods Geomech.; 2003; 27, pp. 745-765. [DOI: https://dx.doi.org/10.1002/nag.295]

27. Solowski, W.T.; Gallipoli, D. Explicit stress integration with error control for the Barcelona basic model. Part I: Algorithms formulation. Comput. Geotech.; 2010; 37, pp. 59-67. [DOI: https://dx.doi.org/10.1016/j.compgeo.2009.07.004]

28. Ma, T.T.; Wei, C.F.; Chen, P. Implicit scheme for integrating constitutive model of unsaturated soils with coupling hydraulic and mechanical behavior. Appl. Math. Mech.; 2014; 35, pp. 1129-1154. [DOI: https://dx.doi.org/10.1007/s10483-014-1859-6]

29. Simo, J.C.; Kennedy, J.G.; Govindjee, S. Non-smooth multisurface plasticity and viscoplasticity. Loading/unloading conditions and numerical algorithms. Int. J. Numer. Methods Eng.; 1988; 26, pp. 2161-2185. [DOI: https://dx.doi.org/10.1002/nme.1620261003]

30. Vaunat, J.; Cante, J.C.; Ledesma, A. A stress point algorithm for an elastoplastic model in unsaturated soils. Int. J. Plast.; 2000; 16, pp. 121-141. [DOI: https://dx.doi.org/10.1016/S0749-6419(99)00033-9]

31. Sloan, S.W. Substepping schemes for the numerical integration of elastoplastic stress-strain relations. Int. J. Numer. Methods Eng.; 2010; 24, pp. 893-911. [DOI: https://dx.doi.org/10.1002/nme.1620240505]

32. Sun, Z.; Wang, C.; Liu, H.; Yin, Y.; Xiao, Y. Bounding surface plasticity model for granular soil and its integration algorithm. Yantu Lixue/Rock Soil Mech.; 2020; 41, pp. 3957-3967.

33. Simo, J.C.; Taylor, R.L. Consistent tangent operators for rate-independent elastoplasticity. Comput. Methods Appl. Mech. Eng.; 1985; 48, pp. 101-118. [DOI: https://dx.doi.org/10.1016/0045-7825(85)90070-2]

34. Sołowski, W.T.; Hofmann, M.; Hofstetter, G.; Sheng, D.; Sloan, S.W. A comparative study of stress integration methods for the Barcelona Basic Model. Comput. Geotech.; 2012; 44, pp. 22-33. [DOI: https://dx.doi.org/10.1016/j.compgeo.2012.03.007]

35. Sheng, D.; Gens, A.; Fredlund, D.G.; Sloan, S.W. Unsaturated soils: From constitutive modelling to numerical algorithms. Comput. Geotech.; 2008; 35, pp. 810-824. [DOI: https://dx.doi.org/10.1016/j.compgeo.2008.08.011]

36. Chen, Z.; Huang, M.; Shi, Z. Application of a state-dependent sand model in simulating the cone penetration tests. Comput. Geotech.; 2020; 127, 103780. [DOI: https://dx.doi.org/10.1016/j.compgeo.2020.103780]

37. Sun, X.; Guo, X.X.; Shao, L.T. A return-mapping algorithm and implementation of thermodynamics-based critical state model. Yantu Lixue/Rock Soil Mech.; 2015; 36, pp. 85-93.

38. Jian, L.I.; Cheng-gang, Z.; Yan, L. Numerical implementation of a bounding surface model for unsaturated expansive clays. Chin. J. Rock Mech. Eng.; 2017; 36, pp. 2551-2562.

39. Sheng, D.; Smith, D.W.; Sloan, S.W.; Gens, A. Finite element formulation and algorithms for unsaturated soils. Part II: Verification and application. Int. J. Numer. Anal. Methods Geomech.; 2003; 27, pp. 767-790. [DOI: https://dx.doi.org/10.1002/nag.296]

40. Gens, A.; Sánchez, M.; Sheng, D. On constitutive modelling of unsaturated soils. Acta Geotech.; 2006; 1, pp. 137-147. [DOI: https://dx.doi.org/10.1007/s11440-006-0013-9]

41. Alonso, E.E.; Gens, A.; Josa, A. A constitutive model for partially saturated soils. Géotechnique; 1990; 40, pp. 405-430. [DOI: https://dx.doi.org/10.1680/geot.1990.40.3.405]

42. Wang, L.J.; Liu, S.H.; Fu, Z.Z.; Li, Z. Coupled hydro-mechanical analysis of slope under rainfall using modified elasto-plastic model for unsaturated soils. J. Cent. South Univ.; 2015; 22, pp. 1892-1900. [DOI: https://dx.doi.org/10.1007/s11771-015-2708-2]

43. De Gennaro, V.; Pereira, J.M. A viscoplastic constitutive model for unsaturated geomaterials. J. Cent. South Univ.; 2013; 54, pp. 143-151. [DOI: https://dx.doi.org/10.1016/j.compgeo.2013.06.005]

44. Wheeler, S.J.; Gallipoli, D.; Karstunen, M. Comments on use of the Barcelona Basic Model for unsaturated soils. Int. J. Numer. Anal. Methods Geomech.; 2002; 26, pp. 1561-1571. [DOI: https://dx.doi.org/10.1002/nag.259]

45. Gens, A. Constitutive Laws. Modern Issues in Non-Saturated Soils; Gens, A.; Jouanna, P.; Schrefler, B.A. Springer: Vienna, Austria, 1995; pp. 129-158. [DOI: https://dx.doi.org/10.1007/978-3-7091-2692-9_2]

46. Villar, A.L.V. Advances on the knowledge of the thermo-hydro-mechanical behaviour of heavily compacted “FEBEX” bentonite. Phys. Chem. Earth Parts A/B/C; 2007; 32, pp. 701-715.

47. Borja, R.; Lee, S. Cam-Clay plasticity, Part 1: Implicit integration of elasto-plastic constitutive relations. Comput. Methods Appl. Mech. Eng.; 1990; 78, pp. 49-72. [DOI: https://dx.doi.org/10.1016/0045-7825(90)90152-C]

48. Borja, R.; Kavazanjian, E. A constitutive model for the stress–strain–time behaviour of ’wet’ clays. Geotechnique; 1985; 35, pp. 283-298. [DOI: https://dx.doi.org/10.1680/geot.1985.35.3.283]

49. Kolditz, O.; Bauer, S.; Bilke, L.; Bttcher, N.; Zehner, B. OpenGeoSys: An open-source initiative for numerical simulation of thermo-hydro-mechanical/chemical (THM/C) processes in porous media. Environ. Earth Sci.; 2012; 67, pp. 589-599. [DOI: https://dx.doi.org/10.1007/s12665-012-1546-x]

50. Thien, B.; Kosakowski, G.; Kulik, D. COTHERM: Modelling fluid-rock interactions in Icelandic geothermal systems. Proceedings of the Egu General Assembly Conference; Vienna, Austria, 27 April–2 May 2014.

51. Wang, Q.; Tang, A.M.; Cui, Y.J.; Delage, P.; Barnichon, J.D.; Ye, W.M. The effects of technological voids on the hydro-mechanical behaviour of compacted bentonite–sand mixture. Soils Found.; 2013; 53, pp. 232-245. [DOI: https://dx.doi.org/10.1016/j.sandf.2013.02.004]