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

Authors: Yu-Jun Zhang; Chao-Shuai Yang; Gang Xu
Journal: Journal of Applied Mathematics
 (2012)

---

## 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.


---

## Body

## 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, CO2, 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 FLAC3D 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 Figure1, within a fracture are in Hertzian contacts, and that a circumferential crack at or outside the contact may be induced by the tensile stress σt. And this crack is described as stress corrosion. Mode I crack velocity for quartz was defined by Dove, given as [9](2.1)vSi–O=AH2Oexp⁡⁡(-ΔHH2ORT)exp⁡⁡(bH2O*K1)(θSi–OH2O)+AOH-exp⁡⁡(-ΔHOH-RT)exp⁡⁡(bOH-*K1)(θSi–OOH-),
where vSi–O is mode I crack velocity caused by chemical dissolution; AH2O and AOH- are the experimentally-determined factor related to temperature; ΔHH2O and ΔHOH- are activation enthalpies; R is the gas constant; T is temperature; bH2O* and bOH-* are the experimentally determined constants derived from the geometry of crack tip; K1 is the stress intensity factor; θSi–OH2O and θSi–OOH- always satisfying θSi–OH2O+θSi–OOH-=1 are the fraction of Si–O reacting with molecular water or hydroxyl ions, and there will be θSi–OH2O=1 and θSi–OOH-=0 at low pH, and θSi–OH2O=0 and θSi–OOH-=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:
(2.2)dEsdt=-(1-Rc)⋅vSi–O,K1≈σt2πr,σt=-(1-2μ)2σa,σa=R̅Rcσ̅,
where r is the distance parallel to the long axis direction of mode I crack caused by σt, and it is assumed to be infinitesimal as well as initial length of crack; μ is the Poisson’s ratio of material; σt is the tensile stress induced by σa which reaches the maximum value just at the edge of the contact; σa is the real stress exerted over the contact area; σ̅ is average macroscopic effective stress. R̅ is the nominal area of the fracture (taking unit value); Rc is the contact-area ratio, and Rc≤R̅.Figure 1
Schematic of fracture compaction induced by microcrack propagation [10].Rc can be calculated via the expression below:
(2.3)Es=Er+(E0-Er)exp⁡⁡[-(Rc-Rc0)a],
where Es and Er are the mean and residual apertures caused by stress corrosion, respectively; E0 is the initial aperture; Rc0 is the relative contact-area ratio at the reference stress; a is empirical constant.Therefore, the evolution of fracture mechanical aperture derived from stress corrosion is(2.4)Et+Δts=Ets+dEsdtΔt.


### 2.2. Effect of Pressure Dissolution on Aperture

The dissolution defined by Yasuhara and Elsworth [10] is expressed as
(2.5)dMdissdt=3πVm2(σa-σc)k+ρgdc24RT,
where dMdiss/dt is the rate of addition of dissolved mass into solution at the interface; Vm is molar volume of the solid; σc is the critical stress that defines stress state where the compaction will effectively halt and reach equilibrium while σa is equal to σc; k+ is the dissolution rate constant of the solid; ρg is the density of solid; dc is the diameter of the asperity contact.And(2.6)k+=k+0exp⁡(-EaRT),σc=Em(1-T/Tm)4Vm,
where k+0 is constant factor; Ea is the activation energy; Em and Tm are the heat and temperature of fusion, respectively.The closure rate of fracture mechanical aperture caused by pressure solution is(2.7)dEpdt=-dMdissdt⋅1ρg⋅1-Rc(π/4)dc2=-3Vm2k+0(1-Rc)(σa-σc)RTexp⁡(-EaRT).And the evolution of fracture mechanical aperture due to pressure solution can be expressed as(2.8)Et+Δtp=Etp+dEpdtΔt.


### 2.3. Fracture Permeability

The fracture spacing in rock mass is assumed to bes, and then the total mechanical aperture for a single fracture at the time of t+Δt is expressed as
(2.9)Et+Δt=(E0+∑ΔtdEs+∑ΔtdEp).So, the hydraulic aperture for a single fracture is [17]
(2.10)et+Δt=Et+Δt2JRC2.5,
where JRC is the roughness coefficient of fractures.Consequently, the equivalent permeability coefficient of fracture in rock mass is [18]:(2.11)Kt+Δt=get+Δt312vs,
where g is gravitational acceleration (9.81 m/s2) and v is kinematics viscosity (the magnitude relative to purified water at 20°C is 1.0×10-6 m2/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(2.12)ϕ=ϕr+(ϕ0-ϕr)exp⁡(f⋅σm′),k=k0exp⁡[c⋅(ϕϕ0-1)]=Fϕkk0,
where ϕ0 and k0 are the porosity and permeability of rock matrix at the stress state of zero, respectively; ϕr is the residual porosity of rock matrix at a high stress state; σm′ is average effective stress; f and c are the experimentally determined parameters, respectively; Fϕk is the modification factor of pore permeability.


## 2.1. Effect of Stress Corrosion on Aperture

Assume that the asperity contacts of brittle materials, schematically shown in Figure1, within a fracture are in Hertzian contacts, and that a circumferential crack at or outside the contact may be induced by the tensile stress σt. And this crack is described as stress corrosion. Mode I crack velocity for quartz was defined by Dove, given as [9](2.1)vSi–O=AH2Oexp⁡⁡(-ΔHH2ORT)exp⁡⁡(bH2O*K1)(θSi–OH2O)+AOH-exp⁡⁡(-ΔHOH-RT)exp⁡⁡(bOH-*K1)(θSi–OOH-),
where vSi–O is mode I crack velocity caused by chemical dissolution; AH2O and AOH- are the experimentally-determined factor related to temperature; ΔHH2O and ΔHOH- are activation enthalpies; R is the gas constant; T is temperature; bH2O* and bOH-* are the experimentally determined constants derived from the geometry of crack tip; K1 is the stress intensity factor; θSi–OH2O and θSi–OOH- always satisfying θSi–OH2O+θSi–OOH-=1 are the fraction of Si–O reacting with molecular water or hydroxyl ions, and there will be θSi–OH2O=1 and θSi–OOH-=0 at low pH, and θSi–OH2O=0 and θSi–OOH-=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:
(2.2)dEsdt=-(1-Rc)⋅vSi–O,K1≈σt2πr,σt=-(1-2μ)2σa,σa=R̅Rcσ̅,
where r is the distance parallel to the long axis direction of mode I crack caused by σt, and it is assumed to be infinitesimal as well as initial length of crack; μ is the Poisson’s ratio of material; σt is the tensile stress induced by σa which reaches the maximum value just at the edge of the contact; σa is the real stress exerted over the contact area; σ̅ is average macroscopic effective stress. R̅ is the nominal area of the fracture (taking unit value); Rc is the contact-area ratio, and Rc≤R̅.Figure 1
Schematic of fracture compaction induced by microcrack propagation [10].Rc can be calculated via the expression below:
(2.3)Es=Er+(E0-Er)exp⁡⁡[-(Rc-Rc0)a],
where Es and Er are the mean and residual apertures caused by stress corrosion, respectively; E0 is the initial aperture; Rc0 is the relative contact-area ratio at the reference stress; a is empirical constant.Therefore, the evolution of fracture mechanical aperture derived from stress corrosion is(2.4)Et+Δts=Ets+dEsdtΔt.


## 2.2. Effect of Pressure Dissolution on Aperture

The dissolution defined by Yasuhara and Elsworth [10] is expressed as
(2.5)dMdissdt=3πVm2(σa-σc)k+ρgdc24RT,
where dMdiss/dt is the rate of addition of dissolved mass into solution at the interface; Vm is molar volume of the solid; σc is the critical stress that defines stress state where the compaction will effectively halt and reach equilibrium while σa is equal to σc; k+ is the dissolution rate constant of the solid; ρg is the density of solid; dc is the diameter of the asperity contact.And(2.6)k+=k+0exp⁡(-EaRT),σc=Em(1-T/Tm)4Vm,
where k+0 is constant factor; Ea is the activation energy; Em and Tm are the heat and temperature of fusion, respectively.The closure rate of fracture mechanical aperture caused by pressure solution is(2.7)dEpdt=-dMdissdt⋅1ρg⋅1-Rc(π/4)dc2=-3Vm2k+0(1-Rc)(σa-σc)RTexp⁡(-EaRT).And the evolution of fracture mechanical aperture due to pressure solution can be expressed as(2.8)Et+Δtp=Etp+dEpdtΔt.


## 2.3. Fracture Permeability

The fracture spacing in rock mass is assumed to bes, and then the total mechanical aperture for a single fracture at the time of t+Δt is expressed as
(2.9)Et+Δt=(E0+∑ΔtdEs+∑ΔtdEp).So, the hydraulic aperture for a single fracture is [17]
(2.10)et+Δt=Et+Δt2JRC2.5,
where JRC is the roughness coefficient of fractures.Consequently, the equivalent permeability coefficient of fracture in rock mass is [18]:(2.11)Kt+Δt=get+Δt312vs,
where g is gravitational acceleration (9.81 m/s2) and v is kinematics viscosity (the magnitude relative to purified water at 20°C is 1.0×10-6 m2/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(2.12)ϕ=ϕr+(ϕ0-ϕr)exp⁡(f⋅σm′),k=k0exp⁡[c⋅(ϕϕ0-1)]=Fϕkk0,
where ϕ0 and k0 are the porosity and permeability of rock matrix at the stress state of zero, respectively; ϕr is the residual porosity of rock matrix at a high stress state; σm′ is average effective stress; f and c are the experimentally determined parameters, respectively; Fϕk 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 Figure2, 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 aren sets of fractures in a fractured porous rock mass, the equilibrium stress equation can be written in the global coordinate system as below:(3.1)dσdt=D[dεdt-m(C1-13Ks)(sw1+Ds1pw1)dpw1dt-m(C2-13Ks)(sw2+Ds2pw2)dpw2dt-mβS3dTdt],
where σ and ε are the total stress and total stain, respectively; D=(C1+C2)-1 is the elastic matrix; mT=[110] is the unit normal column matrix; Ks, βS, and T are the bulk modulus, synthesized thermal expansion coefficient, and temperature of the fractured porous rock mass, respectively;sw1, pw1, Ds1, C1 and sw2, pw2, Ds2, C2 are the saturation degree, water pressure, specific moisture content, and flexibility matrix of rock matrix and fractured network, respectively; t 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 ofdt 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(3.2)-∇T{k1krw1μw∇(pw1+γwz)}+α̅k1krw1μw(pw1-pw2)+A1∂ε∂t+B1∂pw1∂t+E1∂pw2∂t+F1∂T∂t-∇TDt1∇T=0,
where k1 and krw1 are the intrinsic permeability tensor and relative permeability of rock matrix, respectively; ρw, μw, and γw are the density, dynamic viscosity and unit weight of water, respectively; z is the head above some arbitrary datum; α̅ is a parameter determined by the aperture and geometry of fracture; Dt1 is the thermal water diffusivity of rock matrix; and A1, B1, E1, and F1 are the constant matrixes.For the fractured medium, the continuity equation of groundwater is:(3.3)-∇T{k2krw2μw∇(pw2+γwz)}-α̅k1krw1μw(pw1-pw2)+A2∂ε∂t+B2∂pw2∂t+E2∂pw1∂t+F2∂T∂t-∇TDt2∇T=0,
where k2 and krw2 are the intrinsic permeability tensor and relative permeability of fractured medium, respectively; A2, B2, E2 and F2 can be obtained by replacing subscripts 1 and 2 in expressions of A1, B1, E1 and F1 with subscripts 2 and 1; Dt2 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:(3.4)-∇Tλ∇T+(sw1ϕ1V1a+sw2ϕ2V2a)ρwCw(∇TT)+[(1-φ1)CsTρsKsφ1CwTρwKw](sw1+Ds1pw1)∂pw1∂t-{(1-φ1)CsTρsβs+(φ1+φ2)CwTρwβw-[(1-φ1)ρsCs+(φ1+φ2)ρwCw]}∂T∂t+12(1-ϕ1)βsT∂∂t(ui,j+uj,i)δij=0,
where Cw is the specific heat of water; Cs, ρs, and λ are the specific heat, density and thermal conductivity matrix of fractured porous rock mass, respectively; V1a and V2a are the apparent flow velocities of pore water and fracture water, respectively; ui and uj are the displacement components; δij 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:(3.5)Riθiρw∂ci∂t=∇TθiρwDi∇ci-θiρwVi∇ci-Riθiρwχci+(-1)i+1ω̅θiρwD1(C1-C2)-Qci,
where i=1,2 correspond to rock matrix and fractured network, respectively; Ri is the retardation coefficient and is defined as Ri=Vi/Vi*=(1+(ρdi/θi)Kdi); Vi is the apparent velocity of groundwater; Vi* is the transport velocity of radioactive nuclide; ρdi is the dry density of rock matrix or fractured network; Kdi is the distribution coefficient for saturated media; Vi is the apparent velocity of groundwater; θi the volumetric water content; Di is the diffusion tensor; ci is the concentration of solute; Vi is the apparent velocity vector of groundwater; χ is the radioactive decay constant; ω̅ is the coefficient which depends on the fracture aperture and geometry; Qci is the source term.The diffusion tensor can be given by(3.6)Diαβ=αiT|Vi|δαβ+(αiL-αiT)ViαViβ|Vi|+αimτiδαβ,
where αiT is the transversal dispersivity; αiL is the longitudinal dispersivity; |Vi| is the absolute value of the apparent flow velocity; αim is the molecular diffusion coefficient; τi 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].


## 3.1. Equilibrium Stress Equation

Supposing there aren sets of fractures in a fractured porous rock mass, the equilibrium stress equation can be written in the global coordinate system as below:(3.1)dσdt=D[dεdt-m(C1-13Ks)(sw1+Ds1pw1)dpw1dt-m(C2-13Ks)(sw2+Ds2pw2)dpw2dt-mβS3dTdt],
where σ and ε are the total stress and total stain, respectively; D=(C1+C2)-1 is the elastic matrix; mT=[110] is the unit normal column matrix; Ks, βS, and T are the bulk modulus, synthesized thermal expansion coefficient, and temperature of the fractured porous rock mass, respectively;sw1, pw1, Ds1, C1 and sw2, pw2, Ds2, C2 are the saturation degree, water pressure, specific moisture content, and flexibility matrix of rock matrix and fractured network, respectively; t 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 ofdt 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(3.2)-∇T{k1krw1μw∇(pw1+γwz)}+α̅k1krw1μw(pw1-pw2)+A1∂ε∂t+B1∂pw1∂t+E1∂pw2∂t+F1∂T∂t-∇TDt1∇T=0,
where k1 and krw1 are the intrinsic permeability tensor and relative permeability of rock matrix, respectively; ρw, μw, and γw are the density, dynamic viscosity and unit weight of water, respectively; z is the head above some arbitrary datum; α̅ is a parameter determined by the aperture and geometry of fracture; Dt1 is the thermal water diffusivity of rock matrix; and A1, B1, E1, and F1 are the constant matrixes.For the fractured medium, the continuity equation of groundwater is:(3.3)-∇T{k2krw2μw∇(pw2+γwz)}-α̅k1krw1μw(pw1-pw2)+A2∂ε∂t+B2∂pw2∂t+E2∂pw1∂t+F2∂T∂t-∇TDt2∇T=0,
where k2 and krw2 are the intrinsic permeability tensor and relative permeability of fractured medium, respectively; A2, B2, E2 and F2 can be obtained by replacing subscripts 1 and 2 in expressions of A1, B1, E1 and F1 with subscripts 2 and 1; Dt2 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:(3.4)-∇Tλ∇T+(sw1ϕ1V1a+sw2ϕ2V2a)ρwCw(∇TT)+[(1-φ1)CsTρsKsφ1CwTρwKw](sw1+Ds1pw1)∂pw1∂t-{(1-φ1)CsTρsβs+(φ1+φ2)CwTρwβw-[(1-φ1)ρsCs+(φ1+φ2)ρwCw]}∂T∂t+12(1-ϕ1)βsT∂∂t(ui,j+uj,i)δij=0,
where Cw is the specific heat of water; Cs, ρs, and λ are the specific heat, density and thermal conductivity matrix of fractured porous rock mass, respectively; V1a and V2a are the apparent flow velocities of pore water and fracture water, respectively; ui and uj are the displacement components; δij 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:(3.5)Riθiρw∂ci∂t=∇TθiρwDi∇ci-θiρwVi∇ci-Riθiρwχci+(-1)i+1ω̅θiρwD1(C1-C2)-Qci,
where i=1,2 correspond to rock matrix and fractured network, respectively; Ri is the retardation coefficient and is defined as Ri=Vi/Vi*=(1+(ρdi/θi)Kdi); Vi is the apparent velocity of groundwater; Vi* is the transport velocity of radioactive nuclide; ρdi is the dry density of rock matrix or fractured network; Kdi is the distribution coefficient for saturated media; Vi is the apparent velocity of groundwater; θi the volumetric water content; Di is the diffusion tensor; ci is the concentration of solute; Vi is the apparent velocity vector of groundwater; χ is the radioactive decay constant; ω̅ is the coefficient which depends on the fracture aperture and geometry; Qci is the source term.The diffusion tensor can be given by(3.6)Diαβ=αiT|Vi|δαβ+(αiL-αiT)ViαViβ|Vi|+αimτiδαβ,
where αiT is the transversal dispersivity; αiL is the longitudinal dispersivity; |Vi| is the absolute value of the apparent flow velocity; αim is the molecular diffusion coefficient; τi 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 Figure3. 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σv=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 Rc0, R̅, 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.
PropertyRock massVitrified wasteDensity,ρ (kN·m−3)26.725.0Porosity,ϕ10.110.0Permeability,k1/μw (m2·Pa−1·s−1)1.24×10-131.0×10-27Young’s modulus,E (GPa)37.053.0Poisson’s ratio,μ0.30.25Specific heat,C (kJ·kg−1·°C−1)1.00.7Thermal expan. coeff.,β (°C−1)8.8×10-61.0×10-5Thermal conductivity,λ (W·m−1·°C−1)2.85.3Table 2
Parameters for fracture sets used in calculation.
ParameterHorizontal fractureVertical fractureSpacing,S (m)0.30.3Continuity ratio,l11Dip angle,θ(°)090Normal stiffness,kn  (MPa/m)1000.02000.0Shearing stiffness,ks  (MPa/m)500.01000.0Porosity,ϕ20.010.01Permeability,k2/μw (m2/Pa·s)9.7×10-99.7×10-9Table 3
Parameters for stress corrosion and pressure solution.
ParameterUnitValuesEmpirical constant,a—5.0Origin asperity contact-area ratio,Rc0—0.1Nominal asperity contact-area ratio,R̅—1.0Roughness coefficient of fracture, JRC—2.5Factor,AH2Oms−11.12×10-4 TFactor,AOH-ms−12.51×103 TOrigin aperture,E0m0.0125Residual aperture,Erm0.0025Constant,bH2O*N−1 m3/22.69×10-5 (quartz)Constant,bOH-*N−1 m3/21.78×10-5 (quartz)Activation energy,EaJ·mol−17.0×104 (quartz)Heat of fusion,EmJ·mol−18.57×103 (quartz)Activation enthalpy forH2O, ΔHH2OJ·mol−16.6×104Activation enthalpy forOH-, ΔHOH-J·mol−18.27×104Reference dissolution rate constant,k+0mol·m−2 s−11.59Infinitesimal distance from crack tip,rm1.0×10-6Gas constant,RJ·mol−1 K−18.31Temperature of fusion,TmK1883 (quartz)Molar volume,Vmm3 mol−12.27×10-5 (quartz)Fraction of Si–O reacting withH2O, θSi–OH2O—0.99921 (pH = 7)Fraction of Si–O reacting withOH-, θSi–OOH-—0.00079 (pH = 7)The water retention curves of both porous and fracture media conform to the Van Genuchten model, that is,(4.1)sw=(sws-swr)(1+|αψ|n)-m+swr,
where α=3.86×10-6 m−1, n=1.41 for the rock matrix; α=5.26×10-4 m−1, n=2.55 for the fracture system; m=1-1/n; ψ is the water potential head; sws is the maximum saturation with a value of 1.0 while swr 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(4.2)krw=sw2.0.Both the thermal water diffusivities of the rock matrix and fracture system are taken as(4.3)Dt=2.5×10-10m2/s∘C.The vitrified waste is the source term with a diffusive mass flux of radioactive nuclidesQc1=1.44×10-10 mol·kg/m3·s−1. The constants used in the computation concerned with the percolation-migration of nuclide are supposed as follows: the tortuosity coefficients τ1, τ2 are 0.4 and 0.8, respectively; the dispersivities in the longitudinal direction α1L and α2L are 1.0 m and 2.0 m, respectively; the dispersivities in the transversal direction are αiT=αiL/10; the molecular diffusion coefficients α1m and α2m are 1.0×10-9 m2/s and 2.0×10-9 m2/s, respectively; the distribution coefficients Kd1 and Kd2 are 8.0 mL/g and 5.3 mL/g, respectively; the dry densities ρd1 and ρd2 are 23.0 kg/m3 and 21.0 kg/m3, respectively; the parameter ω̅ is 100.0 m−2; the radioactive decay constant χ=ln2/Thalf, where Thalf 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 Figure4. 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 Figures5 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
|dEs/dt| caused by stress corrosion versus time at middle point of right margin of vitrified waste.Figure 6
|dEp/dt| 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 Figures11 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.Pore and fracture water pressures versus time at some nodes for case 1.
(a)
Pore pressure(b)
Fracture pressurePore and fracture water pressures versus time at some nodes for case 2.
(a)
Pore pressure(b)
Fracture pressureContours of pore and fracture pressures in a 2 m × 2 m area at 4 years for case 1 (MPa).
(a)
Pore pressure(b)
Fracture pressureContours of pore and fracture pressures in a 2 m × 2 m area at 4 years for case 2 (MPa).
(a)
Pore pressure(b)
Fracture pressureThe flow vector distributions of pore and fracture water in calculation domain at 4 a for the two cases are presented in Figure15. 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-8 m/s, 1.52×10-8 m/s for case 1 and 2.32×10-8 m/s, 2.77×10-8 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 Figures16 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−3 mol/m3). 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.Nuclide concentrations versus time at some nodes for case 1.
(a)
Pore concentration(b)
Fracture concentrationNuclide concentrations versus time at some nodes for case 2.
(a)
Pore concentration(b)
Fracture concentrationContours 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 concentrationContours of nuclide concentration in a 2 m × 2 m area at 4 years for case 2 (10−3 mol/m3).
(a)
Pore concentration(b)
Fracture concentrationThe 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.Normal stress contours in calculation domain at 4 years for case 1 (MPa).
(a)
Horizontal stress(b)
Vertical stress


## 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.



---

*Source: 983718-2012-05-29.xml*