J. Cent. South Univ. (2012) 19: 2333-2339

DOI: 10.1007/s11771-012-1279-8

FEM analyses for influences of pressure solution on thermo-hydro-mechanical coupling in porous rock mass

ZHANG Yu-jun(张玉军)¹, YANG Chao-shuai(杨朝帅)^{1, 2}

1. State Key Laboratory of Geomechanics and Geotechnical Engineering

(Institute of Rock and Soil Mechanics, Chinese Academy of Sciences), Wuhan 430071, China;

2. China Railway Tunnel Group Co. Ltd., Louyang 471009, China

? Central South University Press and Springer-Verlag Berlin Heidelberg 2012

Abstract: The model of pressure solution for granular aggregate was introduced into the FEM code for analysis of thermo-hydro- mechanical (T-H-M) coupling in porous medium. Aiming at a hypothetical nuclear waste repository in an unsaturated quartz rock mass, two computation conditions were designed: 1) the porosity and the permeability of rock mass are functions of pressure solution; 2) the porosity and the permeability are constants. Then the corresponding numerical simulations for a disposal period of 4 a were carried out, and the states of temperatures, porosities and permeabilities, pore pressures, flow velocities and stresses in the rock mass were investigated. The results show that at the end of the calculation in Case 1, pressure solution makes the porosities and the permeabilities decrease to 10%-45% and 0.05%-1.4% of their initial values, respectively. Under the action of the release heat of nuclear waste, the negative pore pressures both in Case 1 and Case 2 are 1.2-1.4 and 1.01-1.06 times of the initial values, respectively. So, the former represents an obvious effect of pressure solution. The magnitudes and distributions of stresses within the rock mass in the two calculation cases are the same.

Key words: pressure solution; porous medium; thermo-hydro-mechanical coupling; FEM analysis

1 Introduction

If the sedimentary rockmass in deep stratum is in the environment strongly affected by thermo-hydro- mechanical-chemistry coupling (such as a high level radioactive nuclear waste geological repository), the mineral particles comprising rock mass will generate pressure solution (as shown in Fig. 1), which makes the porosity and permeability decrease. Pressure solution incorporates three serial processes [1]. The first is mineral dissolution at high stress contacts of grains; The second is diffusive transport of mineral matter from contacts to the pore space; Finally, this material is precipitated on the lowly stressed surfaces. A great number of studies on pressure solution have been performed. For example, sealing was reported in some tests at elevated temperature (>300 ℃) in sandstone made by TENTHOREY et al [2] and in granite performed by MOORE et al [3], at modest temperature (50-150 ℃) in tuff carried out by LIN et al [4] and    in novaculite executed by POLAK et al [5], and under low effective stress (0.2 MPa) with an acidic solution done by DURHAM et al [6]. Based on the test data,

Fig. 1 Schematic diagram of pressure solution for granular aggregate [1]

YASUHARA and ELSWORTH [7-8] presented a mathematics model which could show three processes of pressure solution and simulate the reductions of porosity (and fracture ratio) and permeability for rock-soil material. This is an academic achievement with much practical value, but it is not applied to numerical analysis of FEM for porous rock mass so far [9-11].

Therefore, the model pressure solution by YASUHARA was introduced into the FEM code of thermo-hydro-mechanical coupling for rock-soil media developed by ZHANG et al [12-13]. In this model, the porosity evolution with pressure solution was established, and the porosity and permeability had to be changed dynamically with the temperature and stress during the computation. Taking a hypothetical nuclear waste repository in an unsaturated rock mass as the calculation example, two calculating cases, in which the initial conditions of temperature, pore water pressure and stress in rock mass are the same, are designed: 1) the porosity and the permeability of rock mass are functions of pressure solution; 2) the porosity and the permeability are constants. Then, the FEM numerical simulation for thermo-hydro-mechanical coupling is carried out, and the states of temperatures, porosities and permeabilities, pore pressures, flow velocities and stresses in the rock mass are investigated. Some new understandings are obtained.

2 Pressure solution and evolution of porosity and permeability

2.1 Pressure solution

The granular aggregate is presented in Fig. 2, and the dissolution rate of the granular aggregate due to pressure is described by YASUHARA et al [1,14]:

               (1)

where dM_diss/dt is the flux of dissolved mass into solution at the grain contact; V_m is molar volume of the solid; σ_c means the critical stress that defines stress state where the compaction will effectively halt and reach equilibrium while σ_a is equal to σ_c; k₊ represents the dissolution rate constant of the solid; ρ_g denotes the density; d_c is the diameter of the grain contact; R is the gas constant; T is absolute temperature.

Fig. 2 Proposed geometric model of grain-to-grain contact [1]

And

                        (2)

                           (3)

                           (4)

where  is constant factor; E_k+ is the activation energy; E_m and T_m are the heat and temperature of fusion, respectively; A₀ and A_c are the maximum areas of grain section and contact, respectively.

Integrating a circular contact of radius r in the range h≤r≤d_c/2, it can be expressed by Fick’s first law:

               (5)

                        (6)

where dM_diff/dt is the diffusive mass flux; ω is the thickness of the water film trapped at the interface; D_b represents the diffusion coefficient, and the term D₀ is a constant; E_D denotes the activation energy; h is infinitesimal length which is substituted to avoid a singularity in integrating the Fick’s first law (h was taken as 1/1 000 of the diameter of grain contact by YASUHARA);  and  are the mineral concentrations of fluid in the interface and pore space, respectively.

It can be defined by the difference between the mineral concentration in the pore space and the equilibrium concentration [14-15]:

              (7)

                        (8)

where  is the precipitation mass flux at the grain face; A_pore is the nominal area of the void (It was taken as  by YASUHARA et al [14]); k_- represents the precipitation rate constant of the mineral; C_eq is the equilibrium solubility of the dissolved mineral; m is the reaction order, which is experimentally constrained;  is the constant factor;  is activation enthalpy.

During the calculation, it is taken as [7,    14]

(9)

where V_p is the volume of the pore space, and

                             (10)

                           (11)

where Q is the flow rate.

2.2 Evolution of porosity and permeability

According to the test by ELIAS et al [16], the strain rate can be calculated by the following expression:

                           (12)

where  and  are porosities at the time t and 0, respectively.

It is convenient to relate dissolution mass fluxes to strain rates by YASUHARA et al [1]. So, the strain rate can be expressed as

                (13)

We made an improvement by taking contribution of the precipitation mass flux at the grain face into the strain rate, and approximately took it as d=d₀. Thus, the strain rate can be modified as

      (14)

So, the total porosity at time of t can be presented by the following expression:

          (15)

Accordingly, the permeability at time t can be calculated [9]:

                        (16)

where k₀ is the original permeability.

The above model of pressure solution for granular aggregate was introduced into the FEM code of thermo-hydro-mechanical coupling for porous media by ZHANG and ZHANG [11].

3 Computation example

As shown in Fig. 3, it is assumed that a canister filled with the vitrified radioactive nuclear waste is disposed at the depth of 2 000 m beneath the ground surface, and the surrounding rock mass is unsaturated quartzite. 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 margin of the vitrified waste to the 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= 53.4 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 and temperature are constant with values of -4.59 MPa and 20 ℃, respectively; the original porosity of the rock mass is 0.3. The environment of coupled T-H-M is to act as the role of pressure solution on the grains of rock mass, which gives rise to changes of porosities and permeabilities. The relative calculating parameters are tabulated in Table 1 and Table 2 (most of the parameters are consulted in Refs. [1, 15]). At the initial state, temperature of rock mass is the uniform value of 20 ℃. The waste continuously releases heat with a constant power of    1 000 W during a period of 4 a.

Fig. 3 Computation model

Table 1 Main computation parameters

The water retention curve of porous medium conforms to the van Genuchten model, that is:

             (17)

where α=3.86×10^{-6 m-1} and β=1.41; γ=1-1/β; ψ is the water potential head; s_ws and s_wr are the maximum saturation with a magnitude of 1.0 and the minimum one with a magnitude of 0.19, respectively.

Table 2 Parameters for pressure solution

The relationship between relative permeability and saturation is

                                 (18)

The thermal water diffusivities of the rock matrix is taken as

D_t=2.5×10^{-12 m2}/(s·℃)                        (19)

The reaction order m in Eq. (7) is approximately taken as 1 during computing.

For the two cases mentioned above, the changes and distributions of temperatures, porosities and permeabilities, pore pressures, flow velocities and stresses in the rock mass were calculated. The main results and analyses are as follows.

The temperature fields in the calculation domain for the two cases are almost the same. Taking Case 1 for instance, the temperatures versus time at the nodes of 432, 433, 434 and 435 are plotted in Fig. 4. It can be seen that in the early 0.2 a, the temperature of buffer rises fast, and then increases slowly. At the end of computing, the temperatures at the nodes of 432, 433, 434 and 435 are 77.8, 61.9, 52.6 and 45.7 ℃, respectively. Temperature contours at 4 a in calculation domain are presented in  Fig. 5.

Fig. 4 Temperatures versus time at some nodes

Fig. 5 Temperature contours in calculation domain at 4 a (?C)

Contours of porosity in rock mass within a range of 2 m×2 m around the canister at 4 a for Case 1 are shown in Fig. 6. At this time, the porosity values at the four points, which have the distances of 0.3, 0.7, 1.1 and   1.5 m from the canister centre in horizontal symmetry axis, are 0.13, 0.08, 0.06 and 0.04, respectively, and are about 43%, 27%, 20% and 13% of their initial values, respectively. It is indicated that the further it is away from canister centre, the more the porosity reduces. The reason of response is that the temperature near the vitrified waste is higher, which gives rise to the thermo effect to diminish the difference between effective stress σ_a and critical stress σ_c slightly, and that the dissolution mass flux, dM_diss/dt, is inversely proportional to the term temperature T. Therefore, porosity in this part decreases slowly, and vice versa. Porosities at the four points mentioned above versus time are shown in Fig. 7. It can be found that the porosities in rock mass linearly descend with time under the effect of pressure solution, but non-synchronously rise with the temperature increasing. Contours of permeability in rock mass at 4 a are shown in Fig. 8, and permeabilities versus time at the four points above in rock mass at 4 a for Case 1 are plotted in Fig. 9, from which it can be seen that the evolutions of permeability are similar to those of porosity, but obviously represent the nonlinearity with time (determined by Eq. (8)). The permeabilities at these four points are 0.60×10^-14, 0.16×10^-14, 0.05×10^-14 and 0.01×10^{-14 m/s, and are about 5%, 1.3%, 0.4% and 0.1% of their initial values, respectively. It is to be noted that the reduction of porosity in rock mass during computing is mainly generated by dM}_diss/dt in the first process of pressure solution while it is quite slightly influenced by the precipitation of dM_prec/dt in the third process (The values of the latter are almost four orders lower than those of the former).

Fig. 6 Contours of porosity in rock mass at 4 a for Case 1

Fig. 7 Porosities versus time at four nodes (x=0 at center of vitrified waste)

Fig. 8 Contours of permeability in rock mass at 4 a for Case 1

Fig. 9 Permeabilities versus time at four nodes (x=0 at center of vitrified waste)

Pore pressures versus time at the four nodes of 432, 433, 434 and 435 for two cases are illustrated in Fig. 10. From Fig. 10, it can be found that, for the effect of pressure solution in Case 1, pore pressures rise substantially with time due to the reduction of porosities and permeabilities. Pore pressures at these four nodes are -6.45, -6.43, -6.31 and -6.17 MPa in Case 1, which are about 1.4 times of their initial values (-4.59 MPa), and are -4.90, -4.89, -4.87 and -4.84 MPa in Case 2, which are about 1.06 times of their initial values, respectively. Contours of pore pressures in rock mass at 4 a for two cases are shown in Fig. 11.

Flow rate vectors of pore water in calculation domain at 4 a for Cases 1 and 2 are displayed in Fig. 12. The ratio of the former to the latter is 1-10. Compared with Case 2, the pore rates decrease by one order for pressure solution in Case 1, and their vector distributions are also quite different. Particularly, the magnitudes of them near the canister are much larger than those in other parts. Taking the nodes 432, 433, 434 and 435 for instance, the pore flow rates are 6.96×10^-10, 3.90×10^-10, 2.02×10^-10 and 1.09×10^{-10 m/s for Case 1, and are} 16.82×10^-10, 24.02×10^-10, 30.34×10^-10 and 36.03×10^{-10 m/s for Case 2, respectively.}

Fig. 10 Pore pressures versus time at 4 nodes for two cases:  (a) Case 1; (b) Case 2

Fig. 11 Contours of pore pressures in rock mass at 4 a for two cases (MPa): (a) Case 1; (b) Case 2

Fig. 12 Flow vectors of pore water in calculation domain at 4 a for two cases: (a) Case 1; (b) Case 2

Generally, the magnitudes and distributions of stresses in calculation domain in two cases are the same due to ignoring the impact of negative pore pressure in rock mass [17]. Taking Case 1 for instance, normal stress contours in calculation domain at 4 a are presented in  Fig. 13. It can be found that the stress fields, influenced by the existence of the vitrified waste and the effect of radiating heat, are different from those only caused by the gravity of overburdened rock mass (the contours of the latter are the horizons). At the time of 4 a, the horizontal stresses at the nodes of 432, 433, 434 and 435 are -12.7, -15.2, -16.8 and -18.1MPa, and the vertical stresses at the nodes of 432, 433, 434 and 435 are -37.5, -39.4, -40.8 and -41.8 MPa, respectively. So reduction range of initial stress, caused by effect of thermal expansion resulting from vitrified waste, tends to be little with the distance from the canister to the beyond.

Fig. 13 Normal stress contours in calculation domain at 4 a for two cases (MPa): (a) Case 1; (b) Case 2

4 Conclusions

1) Temperatures for the two cases are almost the same, and at the time of 4 a, they can reach 30.0-80.0 ℃ in the near field.

2) Pressure solution makes the porosities and the permeabilities decrease to 10%-45% and 0.05%-1.4% of their initial values, which is mainly generated in the first process of pressure solution but slightly influenced in the third process.

3) Under the action of the release heat of nuclear waste, pore pressures in rock mass rise with time, and the values of them at 4 a for Case 1 and Case 2 are 1.2-1.4 and 1.01-1.06 of their initial values, respectively. So the former represents an obvious effect of pressure solution. The pore flow vectors are quite different in the two cases.

4) The magnitudes and distributions of stresses within the rock mass in the two calculation cases are the same for ignoring the impact of negative pore pressure in rock mass.

References

[1] YASUHARA H, ELSWORTH D, POLAK A. A mechanistic model for compaction of granular aggregates moderated by pressure solution [J]. J Geophys Res, 2003, 108(B11): 2530.

[2] TENTHOREY E, COX S, TODD H. Evolution of strength recovery and permeability during fluid-rock reaction in experimental fault zones [J]. Earth and Plane Sci Lett, 2003, 206: 161-172.

[3] MOORE E, LOCKNER A, BYERLEE D. Reduction of permeability in granite at elevated temperatures [J]. Science, 1994, 265: 1558- 1561.

[4] LIN W, ROBERTS J, GLASSLEY W, RUDDLE D. Fracture and matrix permeability at elevated temperatures [C]// Workshop on Significant Issues and Available Data, Near-field/Altered-zone Coupled Effects Expert Elicitation Project. San Francisco, 1997.

[5] POLAK A, ELSWORTH D, YASUHARA H, GRADER A, HALLECK P. Permeability reduction of a natural fracture under net dissolution by hydrothermal fluids [J]. Geophys Res Lett, 2003, 30(20), 2020, doi: 10.1029/2003GL017575.

[6] DURHAM W, BOURCIER W, BURTON E. Direct observation of reactive flow in a single fracture [J]. Water Resour Res, 2001, 37: 1-12.

[7] YASUHARA H, ELSWORTH D. Evolution of permeability in a natural fracture: Significant role of pressure solution [J]. J Geophys Res, 2004, 109, B03204, doi: 10.1029/2003JB002663.

[8] YASUHARA H, ELSWORTH D. Compaction of a rock fracture moderated by competing roles of stress corrosion and pressure solution [J]. Pure Appl Geophys, 2008, 165: 1289-1306.

[9] TARON J, ELSWORTH D. Coupled mechanical and chemical processes in engineered geothermal reservoirs with dynamic permeability [J]. International Journal of Rock Mechanics and Mining Sciences, 2010, 47: 1339-1348.

[10] TARON J,　ELSWORTH D. Constraints on compaction rate and equilibrium in the pressure solution creep of quartz aggregates and fractures: Controls of aqueous concentration [J]. Journal of Geophysical Research, 2010, 115: B07211. doi:10.1029/ 2009JB007118.

[11] ZHANG Yu-jun, ZHANG Wei-qing. 3D thermo-hydro-mechanical model and finite element analyses of dual-porosity fractured medium for ubiquitous-joint rock mass [J]. Sci China Tech Sci, 2010, 53(8): 2172-2182.

[12] ZHANG Yu-jun, YANG Chao-shuai. Coupled hermo-hydro- mechanical–migratory model for dual-porosity medium and numerical analysis [J]. Journal of Central South University of Technology, 2011, 18(4): 1256-1262.

[13] ZHANG Yu-jun, ZHANG Wei-qing, YANG Chao-shuai. FEM analyses for influences of stress corrosion and pressure solution on THM coupling in dual-porosity rock mass [J]. Sci China Tech Sci, 2011, 54(7): 1748-1756.

[14] YASUHARA H, ELSWORTH D, POLAK A, JISHAN L, GRADER A, HALLECK P. Spontaneous switching between permeability enhancement and degradation in fractures in carbonate: Lumped parameter representation of mechanically-and chemically-mediated dissolution [J]. Transport in Porous Media, 2006, 65: 385-409.

[15] LEE D, ELSWORTH D, YASUHARA H, WEAVER J, RICKMAN R. Experiment and modeling to evaluate the effects of proppant-pack diagenesis on fracture treatments [J]. Journal of Petroleum Science and Engineering, 2006, 74: 67-76.

[16] ELIAS P, HAJASH A. Change in quartz solubility and porosity change due to effective stress: An experimental investigation of pressure solution [J]. Geology, 1992, 20: 451-454.

[17] CHIJIMATSU M, KURIKAMI H, ITO A, SUGITA Y. Implication of THM coupling on the near-field of a nuclear waste repository in a homogeneous rock mass [R]. Tokyo: DECOVALES III-Task3-Bench Mark Test 1(BMT1)-Subtask BMT1-B, 2002: 1-43.

(Edited by HE Yun-bin)

Foundation item: Project(2010CB732101) supported by the National Basic Research Program of China; Project(51079145) supported by the National Natural Science Foundation of China; Project(2009BAK53B03) supported by the National Key Technology R&D Program of China

Received date: 2011-06-23; Accepted date: 2011-10-17

Corresponding author: ZHANG Yu-jun, PhD; Tel: +86-27-87198482; E-mail: yjzhang@whrsm.ac.cn