FEM Analyses for T-H-M-M Coupling Processes in Dual-Porosity Rock Mass under Stress Corrosion and Pressure Solution

Abstract

The models of stress corrosion and pressure solution established by Yasuhara et al. were introduced into the 2D FEM code of thermo-hydro-mechanical-migratory coupling analysis for dual-porosity medium developed by the authors. Aiming at a hypothetical model for geological disposal of nuclear waste in an unsaturated rock mass from which there is a nuclide leak, two computation conditions were designed. Then the corresponding two-dimensional numerical simulation for the coupled thermo-hydro-mechanical-migratory processes were carried out, and the states of temperatures, rates and magnitudes of aperture closure, pore and fracture pressures, flow velocities, nuclide concentrations and stresses in the rock mass were investigated. The results show: the aperture closure rates caused by stress corrosion are almost six orders higher than those caused by pressure solution, and the two kinds of closure rates climb up and then decline, furthermore tend towards stability; when the effects of stress corrosion and pressure solution are considered, the negative fracture pressures in near field rise very highly; the fracture aperture and porosity are decreases in the case 1, so the relative permeability coefficients reduce, therefore the nuclide concentrations in pore and fracture in this case are higher than those in case 2.

1. Introduction

The rock mass below thousands of meters from the ground surface, which is dual-porosity medium with pore and fracture as the conduits of transporting, will be the site for the recovery of energy resources and minerals, and for the safe isolation and storage of high level radioactive wastes, CO₂, and so forth. So, the changes in the ambient stress and temperature conditions may affect the permeability characteristics of these conduits combined effects. Stress corrosion [1–4] and pressure solution [5–8], which are corresponding to the fractures, may result in the sealing and degradation of permeability through compaction driven by fracturing (or crushing) of the propping asperities and by dissolution at contacting asperities, respectively.

When local tensile stresses result from the compressive loading of contacting asperities, ‘‘subcritical’’ or ‘‘quasistatic’’ cracking may occur, leading to a time-dependent or progressive failure. Specifically, sub-critical crack growth in the presence of water is believed to be facilitated by chemical reaction, and the resulting process is termed stress corrosion. Pressure solution incorporates three serial processes: mineral dissolution at stressed contacts, diffusive transport of this material along the intervening thin film of water, and ultimate deposition of the mineral matter at the pore wall.

Dove [9] rigorously investigated the dissolution kinetics of quartz under the wide range of temperature and pH conditions and defined an empirical expression of mode I crack velocity resulted from chemical dissolution. Based on the experimental data, Yasuhara and Elsworth [10] investigated the evolution of fracture aperture within a sample of novaculite containing a natural fracture, and they also presented the models which separately account for stress corrosion and pressure solution to describe this response. Taron and Elsworth [11] introduced a kind of coupled thermo-hydro-mechanical-chemical model of dual-porosity medium, in which the influence of pressure solution, shrinkage and dilation of T-H-M, and precipitation and dissolution of mineral on the opening and closure of apertures was considered by simplification. Subsequently, on the basis of modifying permeability and porosity, the characteristics and change mechanisms of permeability within a rock mass containing natural fractures with TOUGHREACT and FLAC^3D were investigated. Taron and Elsworth [12, 13] developed a new model of pressure solution and applied it for numerical simulation of coupled mechanical and chemical processes in engineered geothermal reservoirs with dynamic permeability. Using simplified expression developed by Min et al. [14], the author improved the FEM code of T-H-M coupling of dual-porosity medium, modified the aperture timely and established the evolution of fracture permeability with pressure solution. Aiming at a hypothetical nuclear waste repository in an unsaturated dual-porosity rock mass as the calculation example, the relative numerical simulation [15] was carried out for three cases with different apertures. However, the concentration field of solute was not involved. It is well known that the leakage and diffusion of nuclide from nuclear waste repository are required to study [16]. Consequently, it is imperative to improve the existing model and FEM code, and to perform the analyses of T-H-M-M coupling based on the above work.

The author introduced primarily the models of stress corrosion and pressure solution by Yasuhara into the governing equations presented in [15], and the concentration of solute was involved as well. That is, in the dual-porosity rock mass, the stress field and the temperature field were single, but the water pressures both in pore and fracture are different as well as the concentrations. Therefore, the corresponding simulation for T-H-M-M coupling was constructed. And then, aiming at a hypothetical model for geological disposal of nuclear waste in an unsaturated dual-porosity rock mass, two computation cases were designed: (1) the fracture apertures were changed with the stress corrosion and pressure solution (the porosity of intact rock was also a function of stress); (2) the fracture aperture and the porosity of matrix rock were constants. The corresponding FEM analyses were performed under certain initial conditions of temperature, pore water pressure, in situ stress, and nuclide release intensity, and both the distributions and the changes of temperatures, pore pressures, flow rates, saturations, nuclide concentrations, and stresses in the near field of repository were investigated. Some conclusions were obtained.

2. Modification of Fracture Permeability

2.1. Effect of Stress Corrosion on Aperture

Assume that the asperity contacts of brittle materials, schematically shown in Figure 1, within a fracture are in Hertzian contacts, and that a circumferential crack at or outside the contact may be induced by the tensile stress 𝜎_𝑡. And this crack is described as stress corrosion. Mode I crack velocity for quartz was defined by Dove, given as [9]𝑣_Si−O=𝐴_H2Oexp−Δ𝐻^H2O𝑅𝑇exp𝑏^∗H2O+𝐴_OH−𝐾₁𝜃^H2OSi−Oexp−Δ𝐻^OH−𝑅𝑇exp𝑏^∗OH−𝐾₁𝜃^OH−_Si−O,(2.1) where 𝑣_Si−O is mode I crack velocity caused by chemical dissolution; 𝐴_H2O and 𝐴_OH⁻ are the experimentally-determined factor related to temperature; Δ𝐻_H2O and Δ𝐻_OH⁻ are activation enthalpies; 𝑅 is the gas constant; 𝑇 is temperature; 𝑏^∗_H2O and 𝑏^∗_OH− are the experimentally determined constants derived from the geometry of crack tip; 𝐾₁ is the stress intensity factor; 𝜃^H2OSi_−O and 𝜃^OH−_Si−O always satisfying 𝜃^H2OSi_−O+𝜃^OH−_Si−O=1 are the fraction of Si–O reacting with molecular water or hydroxyl ions, and there will be 𝜃^H2OSi_−O=1 and 𝜃^OH−_Si−O=0 at low pH, and 𝜃^H2OSi_−O=0 and 𝜃^OH−_Si−O=0 at the high one. Consequently, the closure rate of fracture mechanical aperture due to stress corrosion, given by Yasuhara and Elsworth [10], is as follows: 𝑑𝐸^𝑠𝑑𝑡=−1−𝑅^𝑐⋅𝑣_Si−O𝐾₁,𝜎≈𝜎_𝑡√2𝜋𝑟,_𝑡=−(1−2𝜇)2𝜎^𝑎𝜎_𝑎,=𝑅𝑅_𝑐𝜎,(2.2) where 𝑟 is the distance parallel to the long axis direction of mode I crack caused by 𝜎_𝑡, and it is assumed to be infinitesimal as well as initial length of crack; 𝜇 is the Poisson’s ratio of material; 𝜎_𝑡 is the tensile stress induced by 𝜎_𝑎 which reaches the maximum value just at the edge of the contact; 𝜎_𝑎 is the real stress exerted over the contact area; 𝜎 is average macroscopic effective stress. 𝑅 is the nominal area of the fracture (taking unit value); 𝑅_𝑐 is the contact-area ratio, and 𝑅_𝑐≤𝑅.

Figure 1 

Schematic of fracture compaction induced by microcrack propagation [10].

𝑅_𝑐 can be calculated via the expression below: 𝐸^𝑠=𝐸_𝑟+𝐸₀−𝐸_𝑟exp−𝑅_𝑐−𝑅_𝑐0𝑎,(2.3) where 𝐸^𝑠 and 𝐸_𝑟 are the mean and residual apertures caused by stress corrosion, respectively; 𝐸₀ is the initial aperture; 𝑅_𝑐0 is the relative contact-area ratio at the reference stress; 𝑎 is empirical constant.

Therefore, the evolution of fracture mechanical aperture derived from stress corrosion is𝐸^{𝑠𝑡+Δ𝑡}=𝐸^𝑠𝑡+𝑑𝐸^𝑠𝑑𝑡Δ𝑡.(2.4)

2.2. Effect of Pressure Dissolution on Aperture

The dissolution defined by Yasuhara and Elsworth [10] is expressed as 𝑑𝑀_diss𝑑𝑡=3𝜋𝑉^2𝑚𝜎_𝑎−𝜎_𝑐𝑘₊𝜌_𝑔𝑑^2𝑐4𝑅𝑇,(2.5) where 𝑑𝑀_diss/𝑑𝑡 is the rate of addition of dissolved mass into solution at the interface; 𝑉_𝑚 is molar volume of the solid; 𝜎_𝑐 is the critical stress that defines stress state where the compaction will effectively halt and reach equilibrium while 𝜎_𝑎 is equal to 𝜎_𝑐; 𝑘₊ is the dissolution rate constant of the solid; 𝜌_𝑔 is the density of solid; 𝑑_𝑐 is the diameter of the asperity contact.

And 𝑘₊=𝑘⁰⁺exp−𝐸^𝑎𝜎𝑅𝑇_𝑐,=𝐸^𝑚1−𝑇/𝑇_𝑚4𝑉_𝑚,(2.6) where 𝑘⁰⁺ is constant factor; 𝐸_𝑎 is the activation energy; 𝐸_𝑚 and 𝑇_𝑚 are the heat and temperature of fusion, respectively.

The closure rate of fracture mechanical aperture caused by pressure solution is 𝑑𝐸^𝑝𝑑𝑡=−𝑑𝑀𝑑𝑡^diss⋅𝜌1_𝑔⋅1−𝑅^𝑐(𝜋/4)𝑑^2𝑐=−3𝑉^2𝑚𝑘⁰⁺1−𝑅_𝑐𝜎_𝑎−𝜎_𝑐𝑅𝑇exp−𝑅𝐸𝑇^𝑎.(2.7)

And the evolution of fracture mechanical aperture due to pressure solution can be expressed as 𝐸^{𝑝𝑡+Δ𝑡}=𝐸^𝑝𝑡+𝑑𝐸^𝑝𝑑𝑡Δ𝑡.(2.8)

2.3. Fracture Permeability

The fracture spacing in rock mass is assumed to be 𝑠, and then the total mechanical aperture for a single fracture at the time of 𝑡+Δ𝑡 is expressed as 𝐸_𝑡+Δ𝑡=𝐸₀+_Δ𝑡𝑑𝐸^𝑠+_Δ𝑡𝑑𝐸^𝑝.(2.9)

So, the hydraulic aperture for a single fracture is [17] 𝑒_𝑡+Δ𝑡=𝐸^2𝑡+Δ𝑡JRC^2.5,(2.10) where JRC is the roughness coefficient of fractures.

Consequently, the equivalent permeability coefficient of fracture in rock mass is [18]:𝐾_𝑡+Δ𝑡=𝑔𝑒^3𝑡+Δ𝑡12𝑣𝑠,(2.11) where 𝑔 is gravitational acceleration (9.81 m/s²) and 𝑣 is kinematics viscosity (the magnitude relative to purified water at 20°C is 1.0×10⁻⁶ m²/s).

2.4. Effect of Stress on Permeability of Rock Matrix

According to the empirical expression presented by J. P. Davies and D. K. Davies [19], the porosity and permeability of the rock matrix, when the stress in rock matrix changes, can be modified as𝑘𝜙==𝑘₀𝜙_𝑟+𝜙₀−𝜙_𝑟exp𝑓⋅𝜎^𝑚,exp𝑐⋅𝜙𝜙₀−1=𝐹_𝜙𝑘𝑘₀,(2.12) where 𝜙₀ and 𝑘₀ are the porosity and permeability of rock matrix at the stress state of zero, respectively; 𝜙_𝑟 is the residual porosity of rock matrix at a high stress state; 𝜎^𝑚 is average effective stress; 𝑓 and 𝑐 are the experimentally determined parameters, respectively; 𝐹_𝜙𝑘 is the modification factor of pore permeability.

3. Thermo-Hydro-Mechanical-Migratory Coupling Equations for Dual-Porosity Medium

For the dual-porosity medium shown in Figure 2, it can be thought that there exist pore water pressure and fracture water pressure, pore concentration and fracture concentration, but stress field and temperature field are single in the medium. So, one kind of three-dimensional model for coupled thermo-hydro-mechanical-migratory process is created. By omitting the complex deriving steps, the governing equations are given as follows.

Figure 2 

Porous-fractured media.

3.1. Equilibrium Stress Equation

Supposing there are 𝑛 sets of fractures in a fractured porous rock mass, the equilibrium stress equation can be written in the global coordinate system as below:𝑑𝜎𝑑𝑡=𝐷𝑑𝑑𝜀𝑡−𝑚𝐶¹−13𝐾_𝑠𝑠^𝑤1+𝐷_𝑠1𝑝_𝑤1𝑑𝑝^𝑤1−𝑚𝐶₂𝑑𝑡−13𝐾_𝑠𝑠^𝑤2+𝐷_𝑠2𝑝_𝑤2𝑑𝑝^𝑤2𝑑𝑡−𝑚𝛽3^𝑆𝑑𝑑𝑇𝑡,(3.1) where 𝜎 and 𝜀 are the total stress and total stain, respectively; 𝐷=(𝐶₁+𝐶₂)⁻¹ is the elastic matrix; 𝑚^𝑇=[110] is the unit normal column matrix; 𝐾_𝑠, 𝛽_𝑆, and 𝑇 are the bulk modulus, synthesized thermal expansion coefficient, and temperature of the fractured porous rock mass, respectively;𝑠_𝑤1, 𝑝_𝑤1, 𝐷_𝑠1, 𝐶₁ and 𝑠_𝑤2, 𝑝_𝑤2, 𝐷_𝑠2, 𝐶₂ are the saturation degree, water pressure, specific moisture content, and flexibility matrix of rock matrix and fractured network, respectively; 𝑡 is the time.

3.2. Continuity Equation for Groundwater

On the basis of the principle of mass balance, the water volume flowing into an object during a time increment of 𝑑𝑡 is equal to the rate of water accumulation within the object. Assuming that the seepage of water can be expressed by Darcy law, the continuity equation for rock matrix is expressed by−∇^𝑇𝑘¹𝑘_𝑟𝑤1𝜇_𝑤∇𝑝_𝑤1+𝛾_𝑤𝑧+𝛼𝑘₁𝑘_𝑟𝑤1𝜇_𝑤𝑝_𝑤1+𝐴₁−𝑝_𝑤2𝜕𝜀𝜕𝑡+𝐵¹𝜕𝑝_𝑤1𝜕𝑡+𝐸¹𝜕𝑝_𝑤2𝜕𝑡+𝐹¹𝜕𝑇𝜕𝑡−∇^𝑇𝐷_𝑡1∇𝑇=0,(3.2) where 𝑘₁ and 𝑘_𝑟𝑤1 are the intrinsic permeability tensor and relative permeability of rock matrix, respectively; 𝜌_𝑤, 𝜇_𝑤, and 𝛾_𝑤 are the density, dynamic viscosity and unit weight of water, respectively; 𝑧 is the head above some arbitrary datum; 𝛼 is a parameter determined by the aperture and geometry of fracture; 𝐷_𝑡1 is the thermal water diffusivity of rock matrix; and 𝐴₁, 𝐵₁, 𝐸₁, and 𝐹₁ are the constant matrixes.

For the fractured medium, the continuity equation of groundwater is:−∇^𝑇𝑘²𝑘_𝑟𝑤2𝜇_𝑤∇𝑝_𝑤2+𝛾_𝑤𝑧−𝛼𝑘₁𝑘_𝑟𝑤1𝜇_𝑤𝑝_𝑤1+𝐴₂−𝑝_𝑤2𝜕𝜀𝜕𝑡+𝐵²𝜕𝑝_𝑤2𝜕𝑡+𝐸²𝜕𝑝_𝑤1𝜕𝑡+𝐹²𝜕𝑇𝜕𝑡−∇^𝑇𝐷_𝑡2∇𝑇=0,(3.3) where 𝑘₂ and 𝑘_𝑟𝑤2 are the intrinsic permeability tensor and relative permeability of fractured medium, respectively; 𝐴₂, 𝐵₂, 𝐸₂ and 𝐹₂ can be obtained by replacing subscripts 1 and 2 in expressions of 𝐴₁, 𝐵₁, 𝐸₁ and 𝐹₁ with subscripts 2 and 1; 𝐷_𝑡2 is the thermal water diffusivity of fractured medium.

3.3. Energy Conservation Equation

In accordance with the principle of energy conservation, the rate of heat flowing into an object equals the increase of the internal energy within the object. The temperature field is single, and the energy conservation equation takes the form below:−∇^𝑇𝜆∇𝑇++s1_w1−𝜙𝜑₁V^𝑎1₁+s_w2𝜙₂V^𝑎2𝜌_𝑤𝐶_𝑤∇^𝑇𝑇−1−𝜑₁𝐶_𝑠𝑇𝐾𝜌^𝑠_𝑠𝜑₁𝐶_𝑤𝑇𝐾𝜌^𝑤_𝑤𝑠^𝑤1+𝐷_𝑠1𝑝_𝑤1𝜕𝑝𝜕𝑡^𝑤1+1𝐶_𝑠𝑇𝜌_𝑠𝛽_𝑠+𝜑₁+𝜑₂𝐶_𝑤𝑇𝜌_𝑤𝛽_𝑤−1−𝜑₁𝜌_𝑠𝐶_𝑠+𝜑₁+𝜑₂𝜌_𝑤𝐶_𝑤𝜕𝜕𝑇𝑡21−𝜙¹𝛽_𝑠𝑇𝜕𝜕𝑡𝑢^𝑖,𝑗+𝑢_𝑗,𝑖𝛿_𝑖𝑗=0,(3.4) where 𝐶_𝑤 is the specific heat of water; 𝐶_𝑠, 𝜌_𝑠, and 𝜆 are the specific heat, density and thermal conductivity matrix of fractured porous rock mass, respectively; 𝑉^𝑎1 and 𝑉^𝑎2 are the apparent flow velocities of pore water and fracture water, respectively; 𝑢_𝑖 and 𝑢_𝑗 are the displacement components; 𝛿_𝑖𝑗 is the Kronecker’s delta.

3.4. Percolation-Migration Equation

The percolation-migration equation in [20] was improved by us with the new meaning of adding the solute exchange between rock matrix and fractured network due to concentration difference. The new percolation-migration equation is derived from the old one as follows:𝑅_𝑖𝜃_𝑖𝜌_𝑤𝜕𝑐_𝑖𝜕𝑡=∇^𝑇𝜃_𝑖𝜌_𝑤𝐷_𝑖∇𝑐_𝑖−𝜃_𝑖𝜌_𝑤𝑉_𝑖∇𝑐_𝑖−𝑅_𝑖𝜃_𝑖𝜌_𝑤𝜒𝑐_𝑖+(−1)^𝑖+1𝜔𝜃_𝑖𝜌_𝑤𝐷₁𝐶₁−𝐶₂−𝑄_𝑐𝑖,(3.5) where 𝑖=1,2 correspond to rock matrix and fractured network, respectively; 𝑅_𝑖 is the retardation coefficient and is defined as 𝑅_𝑖=𝑉_𝑖/𝑉^∗𝑖=(1+(𝜌_di/𝜃_𝑖)𝐾_di); 𝑉_𝑖 is the apparent velocity of groundwater; 𝑉^∗𝑖 is the transport velocity of radioactive nuclide; 𝜌_di is the dry density of rock matrix or fractured network; 𝐾_di is the distribution coefficient for saturated media; 𝑉_𝑖 is the apparent velocity of groundwater; 𝜃_𝑖 the volumetric water content; 𝐷_𝑖 is the diffusion tensor; 𝑐_𝑖 is the concentration of solute; 𝑉_𝑖 is the apparent velocity vector of groundwater; 𝜒 is the radioactive decay constant; 𝜔 is the coefficient which depends on the fracture aperture and geometry; 𝑄_𝑐𝑖 is the source term.

The diffusion tensor can be given by 𝐷_𝑖𝛼𝛽=𝛼_𝑖𝑇||𝑉_𝑖||𝛿_𝛼𝛽+𝛼_𝑖𝐿−𝛼_𝑖𝑇𝑉^𝑖𝛼𝑉_𝑖𝛽||𝑉_𝑖||+𝛼^𝑖𝑚𝜏_𝑖𝛿_𝛼𝛽,(3.6) where 𝛼_𝑖𝑇 is the transversal dispersivity; 𝛼_𝑖𝐿 is the longitudinal dispersivity; |𝑉_𝑖| is the absolute value of the apparent flow velocity; 𝛼_𝑖𝑚 is the molecular diffusion coefficient; 𝜏_𝑖 is the tortuosity coefficient; 𝛿_𝛼𝛽 is the Kronecker’s delta.

The discretizations both in space and time domains are carried out for the equilibrium equation, the continuity equation, the energy conservation equation, and the percolation-migration equation by Galerkin method, and then the FEM pattern can be obtained.

The models of stress corrosion and pressure dissolution developed by Yasuhare et al. were introduced into the governing equations above for T-H-M coupling in dual-porosity rock mass by the author, and the corresponding algorithm was consulted in [21, 22].

4. Computation Example

The computation model in laboratory scale is shown in Figure 3. A canister filled with the vitrified radioactive nuclear waste is disposed at a certain depth beneath the ground surface, and the surrounding rock mass is quartzite which is also an unsaturated dual-porosity medium. As an approximate simplification, it is treated to be a plane strain problem. A computation region with a horizontal length of 4 m and a vertical length of 8 m is taken. There are 800 elements and 861 nodes in the mesh. From the midpoint at the right margin of the vitrified waste to right, the node numbers are 432, 433, 434, 435, and 436, respectively. The boundary conditions are as follows.

Figure 3 

Computation model.

The free displacement is allowed for the top of computation domain over which the vertical distributed load of 𝜎_𝑣=26.7 MPa is exerted; both the left and right sides are fixed horizontally; the bottom face is fixed vertically; on all the boundary faces the pore pressure, fracture pressure, and temperature are constant with values of −4.59 MPa, −0.46 MPa, and 20°C, respectively. There exist one set of horizontal fractures and one group of vertical ones in the rock matrix, separately. The state of coupled T-H-M is to act as the role of stress corrosion and pressure solution on the fracture aperture. The relative calculating parameters are tabulated in Tables 1, 2 (The parameter values in these two tables are assumed by the authors, but they have reasonable orders.), and Table 3 (All the parameters in this table are taken from [10] except 𝑅_𝑐0, 𝑅, and JRC.). The saturations of rock matrix and fracture system are 0.44 and 0.01, respectively, and the temperature of rock mass is the uniform value of 20°C at the initial state. The waste continuously releases heat with a constant power of 1000 W during a period of 4 years [23].

Table 1 

Main computation parameters.

Table 2 

Parameters for fracture sets used in calculation.

Table 3 

Parameters for stress corrosion and pressure solution.

The water retention curves of both porous and fracture media conform to the Van Genuchten model, that is, 𝑠_𝑤=𝑠_𝑤𝑠−𝑠_𝑤𝑟1+||𝛼𝜓||^𝑛^−𝑚+𝑠_𝑤𝑟,(4.1) where 𝛼=3.86×10⁻⁶ m⁻¹, 𝑛=1.41 for the rock matrix; 𝛼=5.26×10⁻⁴ m⁻¹, 𝑛=2.55 for the fracture system; 𝑚=1−1/𝑛; 𝜓 is the water potential head; 𝑠_𝑤𝑠 is the maximum saturation with a value of 1.0 while 𝑠_𝑤𝑟 is the minimum saturation of which the values are 0.19 for the rock mass and 0.01 for the fracture system, respectively.

The relationship between relative permeability and saturation degree is 𝑘_𝑟𝑤=𝑠^𝑤2.0.(4.2)

Both the thermal water diffusivities of the rock matrix and fracture system are taken as 𝐷_𝑡=2.5×10⁻¹⁰m²/s^∘C.(4.3)

The vitrified waste is the source term with a diffusive mass flux of radioactive nuclides 𝑄_𝑐1=1.44×10⁻¹⁰ mol·kg/m³·s⁻¹. The constants used in the computation concerned with the percolation-migration of nuclide are supposed as follows: the tortuosity coefficients 𝜏₁, 𝜏₂ are 0.4 and 0.8, respectively; the dispersivities in the longitudinal direction 𝛼_1𝐿 and 𝛼_2𝐿 are 1.0 m and 2.0 m, respectively; the dispersivities in the transversal direction are 𝛼_𝑖𝑇=𝛼_𝑖𝐿/10; the molecular diffusion coefficients 𝛼_1𝑚 and 𝛼_2𝑚 are 1.0×10⁻⁹ m²/s and 2.0×10⁻⁹ m²/s, respectively; the distribution coefficients 𝐾_𝑑1 and 𝐾_𝑑2 are 8.0 mL/g and 5.3 mL/g, respectively; the dry densities 𝜌_𝑑1 and 𝜌_𝑑2 are 23.0 kg/m³ and 21.0 kg/m³, respectively; the parameter 𝜔 is 100.0 m⁻²; the radioactive decay constant 𝜒=ln2/𝑇_half, where 𝑇_half is the half life of radioactive nuclide and is taken as 1000 years in the computation. The waste radiates continuously heat with a power of 1000 W during a period of 4 years, and the time step is taken as 100000 s.

For the two cases with different evolutions of fracture aperture above, the change and distribution of the temperatures, displacements, pore pressures, nuclide concentrations and stresses in the rock mass are studied. The analyses of the main computation results are as follows.

The changes of temperatures in calculation region for case 1 and 2 are basically the same. Taking case 1 for instance, the temperatures versus time at nodes 432, 433, 434, and 435 are shown in Figure 4. In the early 0.1 a, the temperature of buffer increases fast, then it grows slowly. At the termination of computation, the temperatures of nodes 432, 433, 434, and 435 are 77.8°C, 62.0°C, 52.5°C, and 45.7°C, respectively.

Figure 4 

Temperatures versus time at some nodes for case 1.

Induced by the stress corrosion and pressure dissolution, the aperture closure rates of the horizontal fracture and the vertical fracture at the midpoint on the right edge of the canister versus time are plotted in Figures 5 and 6, respectively. It can be seen that both the rates due to these two factors climb up, then decline after the peak, and furthermore tend towards stability slowly. The aperture closure rates caused by stress corrosion are almost six orders higher than those caused by pressure solution. This response is similar with the conclusions presented in [10]. Meanwhile, the aperture closure rates of horizontal fractures are larger than those of vertical fractures, and the reason is that the vertical stresses are higher than the horizontal ones in rock mass. The apertures and the asperity contact-area ratios of the horizontal fracture and the vertical fracture at the midpoint mentioned above versus time are presented in Figures 7 and 8, respectively. For the former, the apertures decrease from the original value and then tend towards the residual value. The contact-area ratios of asperities increase also from initial value then towards the nominal value (unit value), and the changes of the values corresponding to the horizontal fractures are more significant. It can be seen in Figure 9 that stress intensity factor ratio on vertical crack is much larger than that on horizontal crack at this midpoint, and both of them reduce over time. It is shown in Figure 10 that at this midpoint, the critical stresses of horizontal fracture and vertical fracture are equal. They decline rapidly at the beginning, and then tend towards constant. This phenomenon is just due to the combined effects of temperature, stress, and chemistry.

Figure 5 

|𝑑𝐸𝑠/𝑑𝑡| caused by stress corrosion versus time at middle point of right margin of vitrified waste.

Figure 6 

|𝑑𝐸𝑝/𝑑𝑡| caused by pressure solution versus time at middle point of right margin of vitrified waste.

Figure 7 

Fracture apertures versus time at middle point of right margin of vitrified waste.

Figure 8 

Contact-area ratios versus time at middle point of right margin of vitrified waste.

Figure 9 

Stress intensity factor ratios versus time at middle point of right margin of vitrified waste.

Figure 10 

Critical stresses versus time at middle point of right margin of vitrified waste.

Pore and fracture pressures at nodes 432, 433, 434, and 435 versus time for case 1 and case 2 are presented in Figures 11 and 12, respectively. It can be seen that negative pore and fracture pressures rise higher for case 1 than those for case 2. Particularly, at node 432 where the effects of stress corrosion and pressure solution are the most intense, the negative fracture pressure reaches a quite large value. The reason of this response is that the reduction of stress corrosion and pressure dissolution on the fracture apertures and the change of pore permeability with time are considered for case 1, while the fracture apertures and the pore permeability remain constants for case 2. The negative pore and fracture pressures at node 432 at 4 a are −12.25 MPa, −7.95 MPa for case 1 and −6.03 MPa, −0.66 MPa for case 2, respectively. Contours of pore and fracture pressures within a range of 2 m × 2 m around the canister at 4 years for case 1 and case 2 are described in Figures 13 and 14, respectively. It is found that the fracture pressures affected by the stress corrosion and pressure dissolution for case 1 have a significant growth around the canister as compared with case 2.

(a) Pore pressure

(b) Fracture pressure

Figure 11 

Pore and fracture water pressures versus time at some nodes for case 1.

(a) Pore pressure

(b) Fracture pressure

Figure 12 

Pore and fracture water pressures versus time at some nodes for case 2.

(a) Pore pressure

(b) Fracture pressure

Figure 13 

Contours of pore and fracture pressures in a 2 m × 2 m area at 4 years for case 1 (MPa).

(a) Pore pressure

(b) Fracture pressure

Figure 14 

Contours of pore and fracture pressures in a 2 m × 2 m area at 4 years for case 2 (MPa).

The flow vector distributions of pore and fracture water in calculation domain at 4 a for the two cases are presented in Figure 15. The fracture flow vectors for case 1, on which the effects of stress corrosion and pressure dissolution are considered, are quite distinguished from those for case 2, especially in the vicinity of canister. Taking the node 432 for instance, flow velocities of pore and fracture are 3.40×10⁻⁸ m/s, 1.52×10⁻⁸ m/s for case 1 and 2.32×10⁻⁸ m/s, 2.77×10⁻⁸ m/s for case 2, respectively.

Figure 15 

Flow vectors of pore and fracture water in calculation domain at 4 years.

Pore and fracture concentrations at nodes 432, 433, 434, and 435 versus time for the two cases are presented in Figures 16 and 17, respectively. Compared with case 2 in which all of the aperture, porosity, and pore permeability are constants, the nuclides both in fracture and pore are gathered largely in case 1 for the reason that both the reduction of aperture due to stress corrosion and pressure dissolution and the compression of porosity due to mean effective stress lead to decreasing the permeabilities of pore and fracture. The nuclide concentrations at nodes 432, 433, 434, and 435 at 4 a for the two cases are 20.18/6.82, 15.42/3.60, 12.06/2.55 and 9.36/1.76 for rock matrix, and 10.86/8.44, 8.39/6.82, 6.28/5.54 and 5.00/4.63 for fracture system, respectively, (the values in the left and right of “/” are for case 1 and 2, resp., and their units are 10⁻³ mol/m³). Contours of pore and fracture concentrations within a range of 2 m × 2 m around the canister at 4 years for case 1 and case 2 are described in Figures 18 and 19, respectively.

(a) Pore concentration

(b) Fracture concentration

Figure 16 

Nuclide concentrations versus time at some nodes for case 1.

(a) Pore concentration

(b) Fracture concentration

Figure 17 

Nuclide concentrations versus time at some nodes for case 2.

(a) Pore concentration

(b) Fracture concentration

Figure 18 

Contours of nuclide concentration in a 2 m × 2 m area at 4 years for case 1 (10−3 mol/m3).

(a) Pore concentration

(b) Fracture concentration

Figure 19 

Contours of nuclide concentration in a 2 m × 2 m area at 4 years for case 2 (10−3 mol/m3).

The differences between the magnitudes and distributions of stresses within the rock mass in the two cases are quite small for the reason that the impacts of negative pore pressure and negative fracture pressure on the mechanical balance are not considered [24]. For instance, normal stress contours in calculation domain at 4 a for case 1 are given in Figure 20. It can be known that the stress fields, influenced by the existence of the vitrified waste and the effect of radiating heat, are distinguished from those caused only by the gravity of rock mass (the contours of the latter are the horizons). At 4 a, the horizontal stress and vertical stress at the midpoint on the right edge of the canister are −0.124 MPa and −26.75 MPa, respectively. The compressive effect is not to be analyzed in this paper.

(a) Horizontal stress

(b) Vertical stress

Figure 20 

Normal stress contours in calculation domain at 4 years for case 1 (MPa).

5. Concluding Remarks

Based on the introduction of stress corrosion and pressure dissolution of fracture aperture as well as the concentration field of solute, the existing governing equations for T-H-M coupling in dual-porosity rock mass were developed to a model for T-H-M-M coupling. Taking a hypothetical model for geological disposal of nuclear waste with a nuclide leakage in an unsaturated dual-porosity rock mass as a calculation example, on the basis of the two cases whether the changes of fracture apertures with stress corrosion and pressure dissolution are considered or not (meanwhile whether the porosity of rock matrix is the stress function or not), the change and distribution of temperatures, rates and magnitudes of aperture closure, pore pressures, flow velocities, nuclide concentrations, and stresses in rock mass were investigated by the two-dimensional FEM simulation for the coupled T-H-M-M processes. It is shown from the computing results that the temperature differences between case 1 and case 2 are not large, and the temperature in near field can reach 30.0~80.0°C at the end of calculation (4 a); the aperture closure rates caused by stress corrosion are almost six orders higher than those produced by pressure solution, and the two kinds of closure rates rise and then reduce, and furthermore tend towards stability; the fracture apertures decrease from the original value and tend towards the residual value while the contact-area ratios of asperities increase from the original value and tend towards the nominal value; the tensile stress and critical stress exerted over cracks decline over time and then tend towards constants; the negative fracture pressures for the case in which the effects of stress corrosion and pressure solution are considered in near field rise more highly than those for the case in which the corresponding effects are not considered, and the differences of flow vectors between the two cases are quite large; the permeabilities of fracture and pore decline resulted from stress corrosion, pressure dissolution, and mean effective stress in case 1, while they are constants in case 2, so the concentrations both in fracture and pore for the former are larger than those for the latter. but the differences between the magnitudes and distributions of stresses within the rock mass in two cases are very small.

However, the models of stress corrosion and pressure solution by Yasuhara et al. are based on the laboratory test for small scale rock, and the application of them in large-scale rock mass engineering remain to be examined. They are applied on FEM analysis with THMM coupling for a hypothetical model of geological disposal of nuclear waste in an unsaturated rock mass by the authors, and the reliability of it is limited in a certain extent. The further research remains to be carried out in the future.

Acknowledgments

This paper is supported by the National Key Basic Research and Development Program of China (973 Project) (Grant no. 2010CB732101), the National Natural Science Foundation of China (Grant no. 51079145), and the National Key Technology R&D Program (Grant no. 2009BAK53B03).