#### Abstract

A set of thermo-hydro-mechanical coupling control equations was established, and the healing process of the joints between Gaomiaozi bentonite (GMZ01) buffer material blocks and the influence of the joint parameters were numerically simulated. The calculations consider the effect of joints on solute migration, the permeability and thermal conductivity of the buffer material, and the evolution of the healing effect. The effects of the joint design parameters, including the type, number, width, splicing form, and average dry density of the joint, are investigated. Studies show that, under an external water head, the joints will become hydraulic priority channels due to their higher permeability, which will shorten the saturation time of the blocks. As the bentonite gradually saturates, the swelling force compresses the joint material. This action improves the overall uniformity of the buffer material and reduces the priority channel effect. Meanwhile, the final average permeability and diffusion coefficient of the buffer material are found to mainly depend on the average dry density of the buffer material. The higher the average dry density of the buffer material is, the lower the final average permeability and diffusion coefficient are, whereas the distribution of joints and the block splicing are less affected by the average dry density of the buffer material. The findings of this study can provide a reference for the design of bentonite buffer material blocks in the repository.

#### 1. Introduction

In the geological disposal of high-level radioactive waste, highly compacted bentonite blocks are often used to fill the space between waste canister and the surrounding rock because bentonite has low permeability, good thermal conductivity, high sorption capacity, and durability in a natural environment [1, 2]. Joints inevitably occur during construction between the canister and blocks, between the surrounding rock and blocks, and between blocks [3, 4]. These joints accelerate the propagation of groundwater and decrease the effectiveness of high-level nuclear waste repository. Therefore, studying the self-sealing process and permeability characteristics of bentonite blocks and joints is the key to testing the reliability of bentonite as an engineering barrier.

Researchers have carried out experimental studies to elucidate the mechanism of joint action. Imbert and Villar [5] studied the hydraulic response of bentonite particle mixtures: during the saturation process, the microstructure of bentonite is rearranged, the swelling force reaches a steady state, and the mixture becomes homogeneous after full saturation. A study by the FEBEX project [6] in Europe shows that the presence of joints delays the action of the expansion force and produces a preferential flow hydraulic channel. Marcial et al. [7] varied the sample density and joint width in a gap sealing experiment to determine the water transmission law between metal canister and bentonite blocks. Some researchers have suggested backfilling construction joints with clay material [8, 9]: the swelling pressure of the blocks compacts the joint material, thus healing joints between bentonite blocks. Chijimatsu et al. [10] found that joints affect the mechanical properties, permeability properties, and water solute nuclide migration law of blocks. However, the experimental study on the multiphysical field evolution law of the blocks and joints and their influence on solute migration is still limited.

The long research cycle of buffer material necessitates the use of numerical calculations for repository design. A large-scale test device named China-Mock-Up was used to investigate transfer and swelling during the THM process of buffer material in [11, 12], and the study results were analysed. In addition, the application of new methods improves the accuracy of numerical calculation. Samaniego et al. [13] presented an alternative approach to treat the coupled problems purely based on the deep energy method (DEM). In order to further quantify the influence of input parameters on the model output, Vu et al. [14] proposed a sensitivity analysis method. Morteza [15, 16] revealed that the Damköhler number in the finite fracture affects the breakthrough of the solute much more considerably compared to that in the matrix and presented an analytical expression for the shear dispersion during solute transport in a coupled system comprised of a capillary tube and a porous medium. However, few studies have been performed on how joints between blocks affect the THM coupling process and their influence on solute migration.

In the present study, the use of GMZ01 bentonite buffer material with joints and bentonite powder as a joint filling material [17] was investigated. A framework for integrated analysis of the thermal, hydraulic, solute migration, and mechanical behaviors is presented in Figure 1. First, by establishing a set of thermo-hydro-mechanical coupling control equations, the healing process of the joints between Gaomiaozi bentonite (GMZ01) buffer material blocks is simulated. Second, the influence of the joint design parameters, such as the type, number, width, splicing form, and dry density, on the properties of the bentonite buffer material was analysed. Then, the quantitative relational expressions of the influence of the joint design parameters on the working performance of buffer materials are summarized. The present study aims to reveal the healing process of the joints and the influence of the joint parameters on the buffer material, which can provide a reference for the design of the repository.

#### 2. THM Coupled Governing Equation

The same convection control equation is used to model the blocks and joints, which are both composed of Gaomiaozi bentonite (GMZ01). However, different dry densities are used to model the blocks and joints, resulting in different calculation parameters. Two assumptions are made in this study: heat transfer induced by air flow is neglected, and the temperature of the soil particles and water are considered to reach equilibrium, such that there is only one temperature variable *T*. The heat transfer equation for unsaturated soil can be written as [3, 11]where ∇ is the Hamiltonian operator, *H* is the heat capacity per unit volume, *λ*_{k} is the heat transfer coefficient, *ρ*_{l} is the density of liquid water, *C*_{l} is the specific heat capacity of water, *u* is the water flow rate vector, and *q* is the heat source intensity.

The thermal conductivity of bentonite, *λ*_{k} (W/(m∙°C)), is approximately linear in the dry density *ρ*_{d} (g/cm^{3}) and the saturation [18, 19]:

Bentonite is characterized by a high density and slow seepage and can therefore be modelled by the Richards equation for seepage [11, 20]:where *θ*_{l} is the water volume fraction, is the gravity acceleration vector, *E*_{l} is an evaporation term, is a source term, *p* is the water pressure, and *K*_{l} is the permeability coefficient.

At a constant air pressure, the suction *s* depends on the saturation, temperature, and dry density of the soil, i.e., [21–23]

#### 3. Stress Path Analyses

The Barcelona basic model is used to analyse the stress paths for the blocks and joints during soil wetting. The maximum swelling force at constant volume produced by bentonite after complete saturation depends on the dry density, initial moisture content, montmorillonite content, and test methods. The result from a test carried out by Ye et al. [8] for the relationship between the swelling pressure and dry density of GMZ01 bentonite with an initial water content of 11% and an initial temperature of 20°C is given below:

The dry densities of the block and joint are 1.7 g/cm^{3} and 1.1 g/cm^{3}, respectively. Using Equation (5), *P*_{sMj} of the joint material and *P*_{sMb} of the block are 0.068 MPa and 5.83 MPa, respectively.

The deformation of the block caused by the change in the suction under low confining pressures is always elastic and can be written as [18, 24]where *s* is the suction, *κ*_{s} is the elastic deformation index for the suction, and *P*_{v0} is a fitting parameter.

Equation (7) is the fitted relationship between the free swelling ratio and the suction of a block with a dry density of 1.7 g/cm^{3}, as measured by Niu et al. [24], from which *κ*_{s} = 0.10 and *P*_{v0} = 5.0 MPa are obtained. Increasing the external pressure reduces the amount of swelling. Denoting the average external pressure on the soil by *P* and replacing *P*_{s} with *P*_{s}−*P* in equation (7) yields the stress-strain relationship for the block during the swelling process under load as [11]where and *s*_{b0} are the specific volume and the suction, respectively, of the block in the initial state.

The compression index of the joint material is taken as *λ*(0) in the numerical calculation, yielding the stress-strain relationship of the joint material during compression aswhere is the reference stress, *N* (0) is the reference specific volume, which can be measured in compression tests [25], and is the specific volume corresponding to the preconsolidation pressure.

Assuming that the soil is linearly elastic for > ,where *p*_{cj} is the preconsolidation pressure of the joint material.

Setting *p'* = 0 in Equation (9) results in = . The following result is thus obtained:

Denoting the elastic compression index of the block as *κ* yields [26, 27],

The values of the parameters used in the diffusion and mechanical models are shown in Table 1 and the specific derivation process can be seen in Appendix A. The symbols used in the calculation are listed in Table 2.

#### 4. Solute Migration Governing Equation

Assuming that solute migration is affected by the THM field but that solute migration does not affect the THM field corresponds to a unidirectional coupling process. The governing equation for solute migration in unsaturated soil is as follows:where is the volumetric water content of the porous medium, is a diffusion term, is a convection term, *u* is the water flow rate vector, represents an adsorption term, is the dry density, *C* is the solute concentration (mol/m^{3}), *C*_{s} is the quantity of solute adsorbed on the solid framework (mol/kg), *qC*_{R} is a source term, *q* is the source volume flow rate per unit volume, *C*_{R} is the concentration of the source solution (mol/kg), denotes the rate of reaction degradation or decay of the considered substance (mol/m^{2}s), and is the rate at which other substances form the considered substance through reaction or decay (mol/m^{2}s):

Under ideal conditions, that is, a constant temperature, a saturated porous medium, a steady-state porosity, and no reaction or decay of the solute, ; in addition, can be defined in the adsorption model, and a retardation coefficient can be introduced as . Equation (13) can be used to obtain the following equation for isothermal solute migration in saturated porous media [30]. The specific derivation process can be seen in Appendix A and the symbols are listed in Table 2:

#### 5. THM Coupling Calculations and Results

##### 5.1. Physical Model and Boundary Conditions

Referring to the China-Mock-Up experiment [1, 8], the overall repository model is a cylinder with a radius of 450 mm, and the central section comprises a heat source with a radius of 150 mm that is surrounded by rock. The buffer material is formed by splicing blocks and joint material. The initial dry densities of the block and the joint material are 1.70 g/cm^{3} and 1.10 g/cm^{3}, respectively. As the blocks are placed at equal intervals, a 1/6 model is used for the calculation. Symmetrical boundary conditions are used at the left and right sides and the upper and lower sides of the model. The model is 300 mm high, and the temperature of the central section is 90°C; the initial temperature of the buffer layer and the ambient temperature are both 20°C. The initial water content of the blocks and joints are both 0.11 and the inner side is set to “no flow.” The deformation of the heat source and the surrounding rock are not considered, and the inner and outer sides of the soil are set to “roller support.” The program COMSOL is used for numerical calculation and the model mesh contains 3456 domain elements, 3480 boundary elements, and 1208 edge elements.

Figure 2 illustrates the specific size of the physical model and the points used to record data during the calculation, where B1, B2, B3, and B4 are points in the block, J1 and J2 are points in a transverse joint, and J3 and J4 are points in a radial joint. In the numerical calculation, the block material is modelled by a nonlinear elastic stress-strain relation, in which the Poisson ratio is fixed and the swelling force is applied as an external stress; the joint material is modelled using an elastoplastic model, the Mohr–Coulomb yield criterion, and an associated flow rule. Table 1 shows the values of the material parameters obtained from the literature. Each physical field is discretized by the quadratic method. The “MUMPS” solver is applied to solve the transient problem. The backward difference formula is used for time stepping with an initial step size of 10^{−8} d.

During the calculation, the outer boundary of the material is set as the normal fixed boundary, and the risk factors such as shear failure and hydraulic fracturing are not considered. Meanwhile, the time effects such as creep and stress relaxation are not taken into consideration.

##### 5.2. Healing Process

The strain in bentonite is directly reflected by the dry density, which determines the thermal conductivity, seepage, and strength characteristics of soil. The calculated results show that, under an external water head, both the joint and block materials gradually saturate from the outside to the inside, as evidenced by the increased rate of saturation at B2 and B4 compared to that at B1 and B3. The final saturation of the block at B1 is 0.72 without joints, which is lower than the corresponding value with joints (0.99). Water absorption by the block creates a swelling pressure that squeezes the joint. After complete saturation of the buffer material, the dry density of both the inner and outer layers of the block reaches 1.62 g/cm^{3}, the dry density of the joint material after compression reaches approximately 1.47 g/cm^{3}, and the final dry densities of blocks at B1 and B2 without joints are 1.71 g/cm^{3} and 1.69 g/cm^{3}, respectively, which are higher than the corresponding values with joints (1.62 g/cm^{3} and 1.61 g/cm^{3}, respectively). The dry density of the joint material increases significantly compared to that of the initial state (1.1 g/cm^{3}), and the overall buffer material becomes more uniform.

The analysis presented above shows that the presence of joints shortens the saturation time of blocks for the following reasons: the joint has a considerably higher permeability than the block, resulting in the formation of a hydraulic channel, and swelling increases the porosity permeability coefficient of the block. The calculation results also show that the permeability of the joint is reduced from 8.0 × 10^{−18} m^{2} in the initial state to approximately 5.0 × 10^{−20} m^{2} after saturation, and the outer joints heal before the inner joints. The permeability of the block is increased from 2.0 × 10^{−21} m^{2} in the initial state to approximately 6.0 × 10^{−21} m^{2} after saturation. The function of the buffer material is to delay the penetration of pollutants [31]. Therefore, there should be less than an order of magnitude difference in the permeability between the joint and the block for effective healing. The data analysis shows that the joint has a significant healing effect and the overall average permeability of the block is considerably reduced after the joints are healed.

The calculated results presented in Figure 3 show that although the joint has a slightly smaller thermal conductivity than the block, the joints do not significantly affect the overall temperature distribution. The heat transfer coefficients of the block and joint both increase with the saturation, and the considerable increase in the thermal conductivity of the joints meets the International Atomic Energy Agencies recommended requirement of exceeding 1 W/m°C.

##### 5.3. Effect of Joints on the Seepage Process

The analysis presented above shows that joints significantly affect the overall permeability of the buffer material, which is the key factor affecting the solute migration rate. The design parameters of the buffer material include the number and width of joints, the splicing form of the blocks, and the dry density of the materials used. An analysis of the influence of the joint design parameters on heat water migration and solute migration is useful for the design of the geological disposal schemes for high-level radioactive waste. The same initial and boundary conditions are used in this section, which were also presented in Section 5.1.

Figure 2 shows three types of joints (circumferential, radial, and transverse) in the buffer material. The effect of the joint type on seepage is assessed by building models with only one joint type and the same overall size of the buffer material: the calculation results are shown in Figure 4. In these models, the width of the circumferential (a) joint is 10 mm, the width of the circumferential (b) joint has a width of 5 mm, the radial joint angle is 2°, and the transverse joint width is 10 mm. A comparative analysis of the calculation results for different joint widths shows that the samples with wider joints require a relatively short time for the centre of the block to reach 95% saturation and have a lower dry density. The average permeability of the block at = 95％ is similar for the radial and transverse joints but less than that of the circumferential (a) joint. It can be concluded that the main function of the circumferential joint is to provide space for the block to swell, which reduces the dry density of the block and thus increases the permeability. The radial and transverse joints lie perpendicular to the surface of the tank. Consequently, these joints act as hydraulic channels, in addition to providing space for the block to swell. The normal joints have a larger influence on seepage than the circumferential joints.

The influence of the number of joints on the permeability was investigated: considering the widespread use of circumferential joints, only the number of radial joints is changed in the calculations, and the influence of transverse joints is neglected. The calculations are performed for 3, 6, 8, and 12 radial joints. A plan view of the buffer layer with 6 radial joints is shown in Figure 5. *ρ*_{dj} and *ρ*_{db} are the dry densities of the joints and blocks at *r* = 0.2 m after saturation, respectively. *ρ*_{da} is the average dry density of the buffer material, which is obtained as a volume-weighted average of the dry densities of the joints and blocks. *K*_{I} and *K*_{F} represent the average permeability of the section at *r* = 0.2 m in the initial state and after full saturation, respectively, where *K*_{F} reflects the overall permeability of the block. Figure 5 shows the calculated results for a radial joint angle of 2°: as the number of joints increases, the joint dry density decreases and the average permeability increases. Therefore, the number of radial joints should be minimized to ensure the barrier performance of the buffer material. Figure 6 shows the calculation model used to investigate the influence of the radial joint width on the permeability. The influence of the radial joint width is investigated by changing the radial joint angle for 6 radial joints and a constant circumferential joint width. Figure 6 shows that, as the radial joint angle increases, the final dry density of the joints decreases and the average permeability increases. Therefore, the precision during construction should be carefully controlled to ensure that the joints are not overly wide.

Figure 7 shows the results obtained using a model with a constant number and width of joints while varying the dry density of the blocks and joints to simulate the effect of a buffer material with a variable dry density. As the joint has a significantly higher permeability than the block, the joint dry density has a more important effect on the initial average permeability *K*_{I} of the buffer layer than the block dry density. The higher the joint dry density is, the lower the initial average permeability is. During self-healing, the block dry density determines the average permeability *K*_{F} after full saturation, that is, the average permeability *K*_{F} decreases as the dry density increases. The expansion force increases with the block dry density, thus squeezing the joint material into a denser state and significantly reducing the overall permeability. Therefore, increasing the dry density of the block can significantly enhance the antiseepage characteristic of the buffer layer.

The data for the average dry density *ρ*_{da} and average permeability *K*_{F} presented in Section 5.3 are regressed to yield the following correlation:

The calculated results presented in Figure 8 show that the final average permeability is approximately exponential in the average dry density in the presence of multiple types of joints. The average permeability decreases as the average dry density increases, which is similar to the trend observed for the existing test results.

#### 6. Effect of Joints on the Solute Migration

An important function of the buffer material is to block the migration of pollutant solutes. This study is based on the assumptions that solute migration does not affect the THM field (via a one-way coupled process) and occurs on a considerably longer time scale than heat conduction and seepage [32, 33]. Therefore, the calculated results for solute migration are analysed in a separate section. The governing equations used in the calculation are shown in Section 3, and the calculation parameters are shown in Table 3.

The splicing mode of the block and dry density of the material are changed, the law of solute migration after self-healing is analysed, and the relationship between the design parameters and the final average diffusion coefficient is determined.

##### 6.1. Evolution of Solute Migration

Figure 9 shows a calculation model based on the following assumptions: when the waste container is damaged, the bentonite is completely saturated under the long-term action of groundwater and the external water head, the temperature of the waste container and the solute concentration at the inner boundary remain constant [28], *C*_{0} = 1 mol/m^{3}, and the initial soil saturation is 1.

Figure 10 shows the solute concentration distribution along the radial joint at different times, where the time unit is years (denoted as a). The solute concentration diffuses from the inside to the outside, and the solute concentration at the same distance gradually increases with increasing diffusion time. These results are in line with expectations.

##### 6.2. Influence of the Block Splicing Form on Solute Migration

Figure 11 is a plan view of the block splicing, in which the inner and outer radial joints have an average width of 10 mm. The staggered and aligned arrangement corresponds to the scheme used for the FEBEX project [8]. A transverse joint has a similar effect as a radial joint and is therefore not considered in this section. The calculated results are shown in Table 4.

*D*_{j} and *D*_{b} are the diffusion coefficients of the joints and blocks at *r* = 0.2 m after deformation, where *D*_{f} is the final average diffusion coefficient for the section at *r* = 0.2 m, and the diffusion coefficient *D* is given by Equation (A.25) in Appendix A. *t*_{c10} is the time for the average solute concentration in the outer seam of *R* = 0.45 m to reach 0.1*c*_{0}, which reflects the overall barrier performance of the buffer material. An analysis of the data in Table 4 shows that the alignment of the outer and inner joints has little effect on solute migration for a fixed joint volume proportion. This result shows that, as the soil is fully saturated at the beginning of the calculation, the joints on the inner and outer sides have healed, the diffusion coefficient of the material is relatively uniform, and the concentration distribution during solute migration is relatively uniform.

##### 6.3. Effect of Dry Density of Joints and Blocks on Solute Migration

The results presented in Figure 12 show that increasing the dry density of the block or joint decreases the diffusion of the solute concentration to the outer boundary, indicating that increasing the dry density of the material can enhance the barrier performance. The data presented in Table 5 are regressed to obtain the correlation for the average dry density *ρ*_{da} and the final average diffusion coefficient *D*_{f} shown in Figure 13. The logarithm of the final average diffusion coefficient decreases approximately linearly with increasing average dry density. The fitting relationship is as follows:

#### 7. Summary and Conclusions

Through the establishment of THM coupling model, the healing process of the joints, and the influence of different joint design parameters, such as joint type, number, width, splicing form, dry density of block, and joint, on the heat transfer, seepage and solute migration of buffer material are studied. The main conclusions drawn are as follows:(1)The presence of joints in the buffer material significantly increases the overall material permeability and reduces the swelling stress to some extent but has little influence on the thermal conductivity.(2)The high swell ability of bentonite facilitates healing by joints during the saturation process and improves the uniformity of the buffer material.(3)Radial and transverse joints considerably affect the permeability of the buffer material; thus, these joints should be excluded from designs to the greatest possible extent.(4)The dry densities of the block and joint determine the final average permeability and diffusion coefficient of the buffer material. At a fixed average dry density, the splicing pattern of the blocks has little effect on the overall permeability and solute migration.

#### Appendix

#### A. The Specific Derivation Process of the Governing Equation

##### 1.1. THM Coupled Governing Equation

The same convection control equation is used to model the blocks and joints, which are both composed of GMZ01 bentonite. However, different dry densities are used to model the blocks and joints, resulting in different calculation parameters. Two assumptions are made in this study: heat transfer induced by air flow is neglected and the temperature of the soil particles and water are considered to reach equilibrium, such that there is only one temperature variable *T*. The heat transfer equation for unsaturated soil can be written as [3, 11]where ∇ is the Hamiltonian operator, *H* is the heat capacity per unit volume, *λ*_{k} is the heat transfer coefficient, *ρ*_{l} is the density of liquid water, *C*_{l} is the specific heat capacity of water, *u* is the water flow rate vector, and *q* is the heat source intensity.

The thermal conductivity of bentonite, *λ*_{k} (W/(m∙°C)), is approximately linear in the dry density *ρ*_{d} (g/cm^{3}) and the saturation [18, 19]:

The heat capacity per unit volume *H* (J/(m^{3}∙°C)) is given as follows:where *n* is the porosity and *C*_{ps} and *C*_{pl} represent the specific heats of the soil particles and water, respectively.

Bentonite is characterized by a high density and slow seepage and can therefore be modelled by the Richards equation for seepage [11, 20]:where *θ*_{l} is the water volume fraction, is the gravity acceleration vector, *E*_{l} is an evaporation term, is a source term, *p* is the water pressure, and *K*_{l} is the permeability coefficient [11, 36]:where *k*_{r} = *S*_{ω}^{4} is the relative permeability.

The permeability *K*_{int} is approximately exponential in the bentonite dry density [11]:where *T*_{rm} is the room temperature in the laboratory and *μ*_{l} is the dynamic viscosity of water, which depends on the temperature (°C) as given below [21, 37].

at a constant air pressure, and the suction *s* depends on the saturation, temperature, and dry density of the soil, i.e., [21–23]where *σ(T)* reflects the effect of the surface tension coefficient of water on the suction and can be expressed as

Here, the relation between *s*_{r} and *S*_{ω} can be described by the Van Genuchten model:

The model parameters *a* and *b* are obtained by fitting the test results of a soil-water characteristic curve for a compacted bentonite block at 20°C obtained by Bai et al. [38]. The change in the air pressure is neglected. Figure 14 shows the bentonite SWCC curve for different temperatures calculated using Equation (A.9).

##### 1.2. Stress Path Analyses

The Barcelona basic model is used to analyse the stress paths for the blocks and joints during soil wetting. The maximum swelling force at constant volume produced by bentonite after complete saturation depends on the dry density, initial moisture content, montmorillonite content, and test methods. The result from a test by Ye et al. [8] for the relationship between the swelling pressure and dry density of GMZ01 bentonite with an initial water content of 11% and an initial temperature of 20°C is given below:

The dry densities of the block and joint are 1.7 g/cm^{3} and 1.1 g/cm^{3}, respectively. Using Equation (A.10), *P*_{sMj} of the joint material and *P*_{sMb} of the block are 0.068 MPa and 5.83 MPa, respectively.

The swelling pressure *P*_{s} of bentonite increases with decreasing suction and reaches a maximum *P*_{sM} at a suction force of 0 [39–41]:

The deformation of the block caused by the change in the suction under low confining pressures is always elastic and can be written as [18, 24]where *s* is the suction, *κ*_{s} is the elastic deformation index for the suction, and *P*_{v0} is a fitting parameter.

Equation (A.13) is the fitted relationship between the free swelling ratio and the suction of a block with a dry density of 1.7 g/cm^{3}, as measured by Niu et al. [24], from which *κ*_{s} = 0.10 and *P*_{v0} = 5.0 MPa are obtained. Increasing the external pressure reduces the amount of swelling. Denoting the average external pressure on the soil by *P* and replacing *P*_{s} with *P*_{s}−*P* in Equation (A.13) yields the stress-strain relationship for the block during the swelling process under load as [11]where and *s*_{b0} are the specific volume and the suction, respectively, of the block in the initial state.

The compression index of the joint material is taken as *λ* (0) in the numerical calculation, yielding the stress-strain relationship of the joint material during compression aswhere is the reference stress, *N* (0) is the reference specific volume, which can be measured in compression tests [25], and is the specific volume corresponding to the preconsolidation pressure.

Assuming that the soil is linearly elastic for > yieldswhere *p*_{cj} is the preconsolidation pressure of the joint material.

Setting *p'* = 0 in Equation (A.15) results in = . The following result is thus obtained.

denoting the elastic compression index of the block as *κ* yields [26, 27]:

The modulus of elasticity of bentonite varies with the temperature. Zhang et al. [30] correlated the elastic modulus with the temperature. Zhang et al. [42] carried out a saturated expansion test of bentonite at different temperatures and found that the expansion rate decreases as the temperature increases: this behavior can be modelled by modifying the elastic modulus in the calculation model. Introducing a temperature correction factor *E*_{T} into the model yields [9, 11]

The Mohr–Coulomb yield criterion can be used to determine when the material yields based on the principal stress, cohesion *c*, and internal friction angle *φ*. Wood [29] reported that the cohesion *c* is proportional to the temperature and the dry density, whereas the internal friction angle *φ* mainly depends on the saturation. Thus, the following relationships can be written:where *C*_{0}, *C*_{d}, *C*_{T}, *φ*_{0}, and *φ*_{T} are obtained by fitting the test results of pure GMZ01 bentonite samples with high water contents [29]. The values of the parameters used in the diffusion and mechanical models are shown in Table 1.

##### 1.3. Solute Migration Governing Equation

Assume that solute migration is affected by the THM field, but that solute migration does not affect the THM field corresponds to a unidirectional coupling process. The governing equation for solute migration in unsaturated soil is as follows:where is the volumetric water content of the porous medium, is a diffusion term, is a convection term, *u* is the water flow rate vector, represents an adsorption term, is the dry density, *C* is the solute concentration (mol/m^{3}), *C*_{s} is the quantity of solute adsorbed on the solid framework (mol/kg), *qC*_{R} is a source term, *q* is the source volume flow rate per unit volume, *C*_{R} is the concentration of the source solution (mol/kg), denotes the rate of reaction degradation or decay of the considered substance (mol/m^{2}s), and is the rate at which other substances form the considered substance through reaction or decay (mol/m^{2}s).

In the absence of the injection of additional solution into the solution domain, *q* in the source term can be considered to be the change in the pore water volume caused by seepage, that is, , for a sufficiently small unit, and *C*_{R}=*C*. The quantity of solute adsorbed on the solid framework [29] can be expressed as , considering the effect of temperature. Then, Equation (A.20) can be written as follows:

Under ideal conditions, that is, a constant temperature, a saturated porous medium, a steady-state porosity, and no reaction or decay of the solute, ; in addition, can be defined in the adsorption model, and a retardation coefficient can be introduced as . Equation (A.21) can be used to obtain the following equation for isothermal solute migration in saturated porous media [30]:

The effective diffusion coefficient *D*_{eff} is composed of a mechanical diffusion coefficient *D*_{m} and a molecular diffusion coefficient *D*_{e}. As the velocity of seepage in bentonite is low, the influence of mechanical diffusion can be neglected, that is, *D*_{m} = 0. Molecular diffusion occurs via the random thermal motion of molecules and is positively correlated with the thermodynamic temperature [20]; thus, the saturation of the porous medium and the permeability *K*_{int} [34] can be written as follows:where *D*_{e0} is the reference value of the molecular diffusion coefficient and *K*_{int0} is the reference permeability.

Sodium-based GMZ bentonite adsorbs radionuclides through ion exchange reactions and surface coordination reactions [43, 44]; the adsorption of Eu(III) in GMZ can be described by the Freundlich equation:where *k*_{d} is an adsorption coefficient (m^{3}/kg) that can be obtained by fitting experimental data [35].

The diffusion coefficient *D* is often normalised by *R*_{d} to produce a direct input parameter for numerical calculation [11]. The expression for *D* is

#### Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.

#### Conflicts of Interest

The authors declare that there are no conflicts of interest regard this paper.

#### Acknowledgments

The authors would like to thank the National Natural Science Foundation of China (52078031) for financially supporting this study.