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

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 dualporosity medium developed by the authors. Aiming at a hypothetical model for geological disposal of nuclear waste in an unsaturated rock mass fromwhich 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.


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 2 , 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,

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 σ t .And this crack is described as stress corrosion.Mode I crack velocity for quartz was defined by Dove, given as 9 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: 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 ; R c is the contact-area ratio, and R c ≤ R. R c can be calculated via the expression below: where E s and E r are the mean and residual apertures caused by stress corrosion, respectively; E 0 is the initial aperture; R c0 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

Effect of Pressure Dissolution on Aperture
The dissolution defined by Yasuhara and Elsworth 10 is expressed as where dM diss /dt is the rate of addition of dissolved mass into solution at the interface; V m 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; d c is the diameter of the asperity contact.And where k 0 is constant factor; E a is the activation energy; E m and T m are the heat and temperature of fusion, respectively.The closure rate of fracture mechanical aperture caused by pressure solution is And the evolution of fracture mechanical aperture due to pressure solution can be expressed as 2.8

Fracture Permeability
The fracture spacing in rock mass is assumed to be s, and then the total mechanical aperture for a single fracture at the time of t Δt is expressed as So, the hydraulic aperture for a single fracture is 17 where JRC is the roughness coefficient of fractures.Consequently, the equivalent permeability coefficient of fracture in rock mass is 18 : where g is gravitational acceleration 9.81 m/s 2 and v is kinematics viscosity the magnitude relative to purified water at 20 • C is 1.0 × 10 −6 m 2 /s .

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 where φ 0 and k 0 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.

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.

Equilibrium Stress Equation
Supposing there are n sets of fractures in a fractured porous rock mass, the equilibrium stress equation can be written in the global coordinate system as below: where σ and ε are the total stress and total stain, respectively; D C 1 C 2 −1 is the elastic matrix; m T 1 1 0 is the unit normal column matrix; K s , β S , and T are the bulk modulus, synthesized thermal expansion coefficient, and temperature of the fractured porous rock mass, respectively; s w1 , p w1 , D s1 , C 1 and s w2 , p w2 , D s2 , C 2 are the saturation degree, water pressure, specific moisture content, and flexibility matrix of rock matrix and fractured network, respectively; t is the time.

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 dt 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 where k 1 and k rw1 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; D t1 is the thermal water diffusivity of rock matrix; and A 1 , B 1 , E 1 , and F 1 are the constant matrixes.
For the fractured medium, the continuity equation of groundwater is: where k 2 and k rw2 are the intrinsic permeability tensor and relative permeability of fractured medium, respectively; A 2 , B 2 , E 2 and F 2 can be obtained by replacing subscripts 1 and 2 in expressions of A 1 , B 1 , E 1 and F 1 with subscripts 2 and 1; D t2 is the thermal water diffusivity of fractured medium.

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: where C w is the specific heat of water; C s , ρ s , and λ are the specific heat, density and thermal conductivity matrix of fractured porous rock mass, respectively; V a 1 and V a 2 are the apparent flow velocities of pore water and fracture water, respectively; u i and u j are the displacement components; δ ij is the Kronecker's delta.

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: where i 1, 2 correspond to rock matrix and fractured network, respectively; R i is the retardation coefficient and is defined as V i is the apparent velocity of groundwater; V * i is the transport velocity of radioactive nuclide; ρ di is the dry density of rock matrix or fractured network; K di is the distribution coefficient for saturated media; V i is the apparent velocity of groundwater; θ i the volumetric water content; D i is the diffusion tensor; c i is the concentration of solute; V i is the apparent velocity vector of groundwater; χ is the radioactive decay constant; ω is the coefficient which depends on the fracture aperture and geometry; Q ci is the source term.
The diffusion tensor can be given by where α iT is the transversal dispersivity; α iL is the longitudinal dispersivity; |V i | 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 percolationmigration 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 .

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.
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 , n 2.55 for the fracture system; m 1 − 1/n; ψ is the water potential head; s ws is the maximum saturation with a value of 1.0 while s wr 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 Both the thermal water diffusivities of the rock matrix and fracture system are taken as    The vitrified waste is the source term with a diffusive mass flux of radioactive nuclides Q c1 1.44 × 10 −10 mol•kg/m 3 •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 m 2 /s and 2.0 × 10 −9 m 2 /s, respectively; the distribution coefficients K d1 and K d2 are 8.0 mL/g and 5.3 mL/g, respectively; the dry densities ρ d1 and ρ d2 are 23.0 kg/m 3 and 21.0 kg/m 3 , respectively; the parameter ω is 100.0 m −2 ; the radioactive decay constant χ ln 2/T half , where T 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.
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.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.
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 −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.
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.55and 9.36/1.76for 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/m 3 .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.
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.124MPa and −26.75 MPa, respectively.The compressive effect is not to be analyzed in this paper.

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.

Journal of Applied Mathematics
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.

2 . 1 wherev
Si-O is mode I crack velocity caused by chemical dissolution; A H 2 O and A OH − are the experimentally-determined factor related to temperature; ΔH H 2 O and ΔH OH − are activation enthalpies; R is the gas constant; T is temperature; b * H 2 O and b * OH − are the experimentally determined constants derived from the geometry of crack tip; K 1 is the stress intensity factor; θ H 2 O Si-O and θ OH − Si-O always satisfying θ H 2 O Si-O θ OH − Si-O 1 are the fraction of Si-O reacting with molecular water or hydroxyl ions, and there will be θ H 2

− 1 Figure 5 :
Figure 5: |dE s /dt| caused by stress corrosion versus time at middle point of right margin of vitrified waste.

Figure 6 :
Figure 6: |dE p /dt| caused by pressure solution versus time at middle point of right margin of vitrified waste.

Figure 7 :Figure 8 :
Figure 7: Fracture apertures versus time at middle point of right margin of vitrified waste.

5 )Figure 9 :Figure 10 :
Figure 9: Stress intensity factor ratios versus time at middle point of right margin of vitrified waste.

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

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

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

Figure 14 :
Figure 14: Contours of pore and fracture pressures in a 2 m × 2 m area at 4 years for case 2 MPa .

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

Figure 16 :
Figure 16: Nuclide concentrations versus time at some nodes for case 1.

Figure 17 :
Figure 17: Nuclide concentrations versus time at some nodes for case 2.

Figure 20 :
Figure 20: Normal stress contours in calculation domain at 4 years for case 1 MPa .

Table 2 :
Parameters for fracture sets used in calculation.

Table 3 :
Parameters for stress corrosion and pressure solution.