J. Cent. South Univ. (2017) 24: 720-725

DOI: 10.1007/s11771-017-3473-1

Effects of porosity heterogeneity on chemical dissolution-front instability in fluid-saturated rocks

ZHAO Chong-bin(赵崇斌)¹, Peter SCHAUBS², Bruce HOBBS³

1. Computational Geosciences Research Centre, Central South University, Changsha 410083, China;

2. CSIRO Mineral Resource National Research Flagship, P. O. Box 1130, Bentley, WA 6102, Australia;

3. School of Earth and Environment, The University of Western Australia, Crawley, WA 6009, Australia

Central South University Press and Springer-Verlag Berlin Heidelberg 2017

Abstract: Homogeneity and heterogeneity are two totally different concepts in nature. At the particle length scale, rocks exhibit strong heterogeneity in their constituents and porosities. When the heterogeneity of porosity obeys the random uniform distribution, both the mean value and the variance of porosities in the heterogeneous porosity field can be used to reflect the overall heterogeneous characteristics of the porosity field. The main purpose of this work is to investigate the effects of porosity heterogeneity on chemical dissolution front instability in fluid-saturated rocks by the computational simulation method. The related computational simulation results have demonstrated that: 1) since the propagation speed of a chemical dissolution front is inversely proportional to the difference between the final porosity and the mean value of porosities in the initial porosity field, an increase in the extent of the porosity heterogeneity can cause an increase in the mean value of porosities in the initial porosity field and an increase in the propagation speed of the chemical dissolution front. 2) An increase in the variance of porosities in the initial porosity field can cause an increase in the instability probability of the chemical dissolution front in the fluid-saturated rock. 3) The greater the mean value of porosities in the initial porosity field, the quicker the irregular morphology of the chemical dissolution front changes in the supercritical chemical dissolution systems. This means that the irregular morphology of a chemical dissolution front grows quicker in a porosity field of heterogeneity than it does in that of homogeneity when the chemical dissolution system is at a supercritical stage.

Key words: porosity heterogeneity; chemical dissolution; front instability; computational simulation; porous rocks

1 Introduction

Homogeneity and heterogeneity are two totally different concepts in nature. The state of being homogeneous is referred to as homogeneity, while the state of being heterogeneous is referred to as heterogeneity. Thus, a homogeneous substance is comprised of same or similar constituents, while a heterogeneous substance is comprised of different or dissimilar constitutes of a conglomerate mass, which are connected, and can be distinguished from one another. At the particle length scale, rocks exhibit strong heterogeneity in their constituents and porosities. However, in terms of physical and mechanical properties (such as porosity and permeability), rocks can also be viewed as homogeneous substances at either the laboratory sample length scale or even larger length scales. The theoretical basis behind this viewpoint is that from the homogenization theory, a heterogeneous substance can be treated as an equivalent homogeneous one at an appropriate length scale. This is the reason why the concept of the representative elementary volume (REV) is widely used in the classical continuum mechanics [1]. From this point of view, the main advantage of using the continuum mechanics approach is that the physical and mechanical description of the rock can be significantly simplified, so that it is possible to obtain analytical solutions for some kinds of scientific and engineering problems. Nevertheless, the main disadvantage of using the continuum mechanics approach is that only can the macroscopic response of a system be obtained, but the microscopic response cannot be obtained at the particle length scale. Thus, in the case of studying the formation mechanisms of orebodies and mineral deposits, as well as hydrocarbon accumulation and reservoir characterization (at smaller length scales), it is necessary to consider the heterogeneity of the rocks.

Based on the continuum mechanics approach, the chemical dissolution-front instability problem in fluid-saturated rocks has been theoretically and computationally studied over a few decades [2-20]. Many factors, such as: 1) the reactive surface area of dissolvable rock particles; 2) the mineral dissolution ratio of the rock; 3) the mechanical dispersion of reactants in the rock; 4) the anisotropy of rocks; 5) the compressibility of both the pore fluid and rock; 6) the use of different permeability models; and 7) the temperature variation in the rock, were considered in both the theoretical analyses and the numerical simulations of chemical dissolution-front instability problems in fluid-saturated rocks [11, 20]. The related theoretical and numerical results of chemical dissolution- front instability problems in fluid-saturated rocks have cast new light on the dynamic mechanisms of the ore- forming processes and the karst stratum formation within the upper crust of the Earth.

By means of the linear stability theory, the rock involved in a chemical dissolution system is first treated as a homogeneous medium, and then a small perturbation is commonly applied to the system for the purpose of investigating the stability of the system. Since the source of such a small perturbation may come from the heterogeneity of rocks, it is desirable to investigate how the rock heterogeneity can affect the chemical dissolution-front propagation in fluid-saturated rocks. This issue has not been addressed up to date, so that it becomes the main purpose of carrying out this study. The outcome of this study will answer the following realistic question: how can the rock heterogeneity affect the chemical dissolution-front propagation in fluid-saturated rocks when the chemical dissolution system is in a supercritical state?

2 Brief mathematical description of problem

The chemical dissolution-front instability problem in a fluid-saturated (rigid) rock can be treated as a coupled nonlinear problem involving porosity evolution, pore-fluid flow and reactive single-solute transport in the fluid-saturated (rigid) rock. Based on the related scientific principles and the mass conservation of the pore fluid and the chemical species (i.e. the dissolvable mineral), the governing partial differential equations, which are used to describe the chemical dissolution-front instability phenomenon, have been derived in the previous studies [18, 19]. For the sake of completeness of this paper, they are briefly given below. However, if one is interested in the detailed deriving process of these governing partial differential equations, please refer to the related open publications [18, 19].

Under isothermal conditions, the dimensionless governing partial differential equations that are used to describe the chemical dissolution-front instability phenomenon can be expressed as follows [18, 19]:

                      (1)

      (2)

                     (3)

where f is the porosity of the rock; and are the dimensionless pore-fluid pressure and the dimensionless concentration of the solute (i.e. the dissolvable mineral in the rock), respectively; ε is the mineral dissolution ratio of the chemical dissolution system; τ is the dimensionless time to describe the slowness of the chemical dissolution that takes place in the chemical dissolution system; f_f is the final (i.e. maximum) porosity of the rock after the completion of dissolvable mineral dissolution; I is a unit second-order tensor (i.e. a 2 by 2 unit matrix for a two- dimensional problem); is the dimensionless second-order tensor due to pure dispersion;is a dimensionless quantity that has the same definition as used in the previous studies [18, 19].

From the previous study [18, 19], the dimensionless second-order tensor due to pure dispersion can be expressed in the following form:

      (4)

where and are the dimensionless transversal and longitudinal dispersivities of the chemical species (i.e. the dissolved mineral in the pore space of the rock); is a dimensionless quantity, which can be expressed as follows:

(under isothermal conditions)(5)

where k(f) is the permeability of the rock; k(f_f) is the maximum permeability of the rock in correspondence to the final porosity, f_f.

For the chemical dissolution-front propagation problem of a rectangular computational domain, as shown in Fig. 1, the corresponding boundary conditions can be expressed in the following dimensionless forms:

(at                                (6)

(at        (7)

(at and  (8)

whereis the dimensionless pore-fluid pressure gradient at the entrance of the chemical dissolution system; is the dimensionless pore-fluid pressure at the exit of the chemical dissolution system.

The initial conditions for this problem are: expect at where and expect at where Note that f₀ is the initial porosity field. The physical meaning of these initial and boundary conditions is that the chemical dissolution system in the computational domain is initially in an equilibrium state, and the fresh pore-fluid (i.e. water) is injected at the entrance of the chemical dissolution system, so that the chemical equilibrium state is broken and the chemical dissolution reaction starts from the entrance of the chemical dissolution system. As a result, the chemical dissolution- front gradually propagates from the entrance to the exit of the chemical dissolution system. From the previous studies [18, 19], such a propagation process of the chemical dissolution-front is fully controlled by the governing partial differential equations that are expressed by Eqs. (1) to (3). Specifically, the location of a chemical dissolution-front can be expressed by that of the moving porosity front in the chemical dissolution system.

3 Effects of porosity heterogeneity on chemical dissolution-front propagation in fluid-saturated rocks

In the process of dealing with the computational simulation of chemical dissolution-front propagation in fluid-saturated rocks, an amalgamation procedure consisting of using the finite element and finite difference methods has been developed and well verifiedover the past decade [18, 19]. In the developed amalgamation procedure, the finite element method [21] is used to discretize the geometry of a computational domain, while the finite difference method is used to discretize the computational time of a computational model. Due to the versatility and robustness of the developed amalgamation procedure, it can be also used to investigate the effects of porosity heterogeneity on chemical dissolution-front propagation in fluid-saturated rocks.

For the computational simulation of chemical dissolution-front propagation in a fluid-saturated rock, it is necessary to know the specific form of the dimensionless quantity represented by (in Eq. (2)) during the computational simulation of the corresponding numerical solutions [18, 19]. In order to consider the porosity dependence of solute diffusion in a fluid-saturated rock, it is commonly assumed that where D₀ is the solute diffusion coefficient in pure water [1, 2]. In addition, the dimensionless longitudinal dispersivity (i.e. of the chemical species is assumed to be 0.2, while the dimensionless transversal dispersivity (i.e. of the chemical species is assumed to be 0.02 in the forthcoming computational simulation.

In order to simulate the chemical dissolution-front propagation problem shown in Figs. 1 and 2, the computational domain of the problem is discretized into 29601 four-node square elements with 30000 nodal points in total. The Kozeny–Carman formula is used to express the porosity-permeability relationship in the computational simulation. The final porosity (i.e. which is associated with the complete dissolution of the dissolvable mineral in the chemical dissolution system, is equal to 0.2. For the purpose of considering supercritical chemical dissolution systems, the Zhao number of the chemical dissolution system is equal to 5 (i.e. Zh=5.0 or The mineral dissolution ratio is equal to 0.001, while the dimensionless time-step length is also equal to 0.001 for all the forthcoming computational simulations. For the purpose of investigating the effects of porosity heterogeneity on chemical dissolution-front propagation in fluid-saturated rocks, four random initial porosity fields (shown in Fig. 2) of uniform distribution are considered in the four different computational models. Note that the assumed four random initial porosity fields are used to approximately represent some potential heterogeneous distributions of porosities in the real rocks. The porosity variation ranges of these four random initial porosity fields arein case 1, in case 2, in case 3 and in case 4, respectively. This means that the mean values of the porosities are respectively 0.1005, 0.105, 0,12 and 0.14, while the variances of the porosities are respectively andfor the four random initial porosity fields. The variance of porosities in a heterogeneous porosity field can be used to express the extent of porosity heterogeneity. It is obvious that with the increase of the extent of the porosity heterogeneity (i.e. from case 1, case 2, case 3 to case 4), there is an increase in the propagation speed of the chemical dissolution front. The reason for this phenomenon is due to the fact that from the previous theoretical result [19], the propagation speed of a chemical dissolution front is inversely proportional to the difference between the final porosity (which is associated with the complete dissolution of the dissolvable mineral in the chemical dissolution system) and the mean value of porosities in the initial porosity field. Since the mean value of porosities in the initial porosity field increases from 0.1005, 0.105, 0,12 to 0.14 when the case number increases from 1, 2, 3 to 4, the corresponding difference between the final porosity (i.e. and the mean value of porosities in the initial porosity field is equal to 0.0995, 0.095, 0.08 and 0.06 respectively. Therefore, among the four cases under consideration, the chemical dissolution front propagates fastest in case 4, while it propagates lowest in case 1. Another interesting phenomenon that can be observed from the numerical results (shown in Fig. 3) is that an increase in the variance of porosities in the initial porosity field can cause an increase in the instability probability of the chemical dissolution front in the fluid-saturated rock. This means that the chemical dissolution front in case 4 is more unstable that that in case 1 when it propagates in the chemical dissolution system. To explain this phenomenon theoretically, it is necessary to approximately calculate the critical Zhao number of the chemical dissolution system (in the case of using the following formula [19]:

 (9)

where is the critical Zhao number of the chemical dissolution system; is the permeability of the rock in correspondence to the mean value of porosities in the initial porosity field, ; is the maximum permeability of the rock in correspondence to the final porosity,

Fig. 1 Geometry and boundary conditions of chemical dissolution problem in a fluid-saturated rock

Fig. 2 Four different random distributions of porosities in initial heterogeneous porosity field:

Fig. 3 Effects of porosity heterogeneity on chemical dissolution-front propagation in fluid-saturated rock (porosity):

By substituting both the related mean value of porosities in the initial porosity field and the final porosity value into Eq. (9), the corresponding critical Zhao number can be evaluated as 1.7733, 1.8193, 2.0239 and 2.4926 for cases 1, 2, 3 and 4, respectively. Since the Zhao numbers of these four cases are identical, the small value of the critical Zhao number of a system means that the chemical dissolution front of the system is easy to become unstable. However, from the linear stability theory [19], the porosity heterogeneity can be viewed as a kind of the initial perturbation (that is applied to the chemical dissolution system). When a chemical dissolution system is in a supercritical state, the morphology of the chemical dissolution front depends not only on the magnitude of the applied initial porosity perturbation, but also on the dimensionless growth rate of the applied initial porosity perturbation in the chemical dissolution system. Generally, the applied initial porosity perturbation can be expressed as where is the perturbation ratio, which is defined as the ratio of the applied initial porosity perturbation to the initial porosity of the chemical dissolution system. Thus, the modified initial porosity before running a computational model is: , so that the mean value of porosities in the initial porosity field (i.e. ) can be expressed as . This implies that the mean value of porosities in the initial porosity field is directly proportional to the perturbation ratio. By means of the above-mentioned relationship between the perturbation ratio and the mean value of porosities in the initial porosity field, the corresponding perturbation ratios are equal to 0.01, 0.1, 0.4 and 0.8 for cases 1, 2, 3 and 4, respectively.

From the previous theoretical study [19], the dimensionless growth rate of an applied small perturbation in the chemical dissolution system can be expressed in the case of the fundamental mode (i.e. in the case of the dimensionless wavenumber being unity) as follows:

 (10)

where is the dimensionless growth rate of an applied small perturbation in the chemical dissolution system; Zh is the Zhao number of the chemical dissolution system. Other quantities have the same meanings as those mentioned previously.

Figure 4 shows the variations of the dimensionless growth rate with mean values (denoted by Phi_m) of porosities in the heterogeneous porosity fields due to three different Zhao numbers, namely Zh=1, 5 and 10 respectively. It is obvious that with the increase of the Zhao number, the dimensionless growth rate of the applied initial porosity perturbation increases accordingly. In the case of Zh=5, the dimensionless growth rate of the applied initial porosity perturbation is almost unchanged when the mean value of porosities in the initial porosity field varies from 0.1005 (in case 1) to 0.14 (in case 4). For example, the dimensionless growth rates of the applied initial perturbation are equal to 23.48, 23.55, 23.14 and 20.44 for cases 1, 2, 3 and 4 respectively. Since the perturbation ratios are equal to 0.1 and 0.4 for cases 2 and 3, the irregular morphology of the chemical dissolution front changes quicker in case 3 than in case 2, because the growth of the irregular morphology depends largely on the product of the perturbation ratio and the value of This indicates that the greater the mean value of porosities in the initial porosity field, the quicker the irregular morphology of the chemical dissolution front changes in the supercritical chemical dissolution systems. Thus, the irregular morphology ofthe chemical dissolution front changes quicker in a heterogeneous porosity field than it does in a homogeneous porosity field when the chemical dissolution system is at a supercritical stage.

Fig. 4 Variations of dimensionless growth rate with mean values of porosities in heterogeneous porosity fields due to different Zhao numbers

4 Conclusions

The effects of porosity heterogeneity on chemical dissolution-front instability in fluid-saturated rocks have been investigated by computational simulation methods. When the heterogeneity of porosity obeys the random uniform distribution, both the mean value and the variance of porosities in a heterogeneous porosity field can be used to reflect the overall characteristics of the heterogeneous porosity field. Such a treatment can provide a theoretical foundation for understanding how the porosity heterogeneity affects chemical dissolution-front instability in fluid-saturated rocks.

The related theoretical results demonstrated that: 1) since the propagation speed of a chemical dissolution front is inversely proportional to the difference between the final porosity and the mean value of porosities in the initial porosity field, an increase in the extent of the porosity heterogeneity can cause an increase in the mean value of porosities in the initial porosity field and an increase in the propagation speed of the chemical dissolution front. 2) An increase in the variance of porosities in the initial porosity field can cause an increase in the instability probability of the chemical dissolution front in the fluid-saturated rock. 3) The greater the mean value of porosities in the initial porosity field is, the quicker the irregular morphology of the chemical dissolution front changes in the supercritical chemical dissolution systems. This means that the irregular morphology of a chemical dissolution front grows quicker in a heterogeneous porosity field than it does in a homogeneous porosity field when the chemical dissolution system is at a supercritical stage.

References

[1] BEAR J. Dynamics of fluids in porous media [M]. Amsterdam: Elsevier, 1972: 636.

[2] CHADAM J, HOFF D, MERINO E, ORTOLEVA P, SEN A. Reactive infiltration instabilities [J]. IMA Journal of Applied Mathematics, 1986, 36: 207-221.

[3] CHADAM J, ORTOLEVA P, SEN A. A weekly nonlinear stability analysis of the reactive infiltration interface [J]. IMA Journal of Applied Mathematics, 1988, 48: 1362-1378.

[4] CHEN J S, LIU C W. Numerical simulation of the evolution of aquifer porosity and species concentrations during reactive transport [J]. Computers and Geosciences, 2002, 28: 485-499.

[5] CHEN J S, LIU C W, LAI G X, NI C F. Effects of mechanical dispersion on the morphological evolution of a chemical dissolution front in a fluid-saturated porous medium [J]. Journal of Hydrology, 2009, 373: 96-102.

[6] COHEN C E, DING D, QUINTARD M, BAZIN B. From pore scale to wellbore scale: impact of geometry on wormhole growth in carbonate acidization [J]. Chemical Engineering Science, 2008, 63: 3088-3099.

[7] GOLFIER F, ZARCONE C, BAZIN B, LENORMAND R, LASSEUX D, QUINTARD M. On the ability of a Darcy-scale model to capture wormhole formation during the dissolution of a porous medium [J]. Journal of Fluid Mechanics, 2002, 457: 213-254.

[8] HINCH E J, BHATT B S. Stability of an acid front moving through porous rock [J]. Journal of Fluid Mechanics, 1990, 212: 279-288.

[9] KALIA N, BALAKOTAIAH V. Modeling and analysis of wormhole formation in reactive dissolution in carbonate rocks [J]. Chemical Engineering Science, 2007, 62: 919-928.

[10] KALIA N, BALAKOTAIAH V. Effect of medium heterogeneities on reactive dissolution of carbonates [J]. Chemical Engineering Science, 2009, 64: 376-390.

[11] LAI K H, CHEN J S, LIU C W, YANG S H. Effect of permeability-porosity functions on simulated morphological evolution of a chemical dissolution front in a fluid-saturated porous medium [J]. Hydrological Processes, 2014, 28: 16-24.

[12] ORMOND A, ORTOLEVA P. Numerical modeling of reaction- induced cavities in a porous rock [J]. Journal of Geophysical Research, 2000, 105: 16737-16747.

[13] ORTOLEVA P, CHADAM J, MERINO E, SEN A. Geochemical self-organization II: The reactive-infiltration instability [J]. American Journal of Science, 1987, 287: 1008-1040.

[14] RENARD F, GRATIER J P, ORTOLEVA P, BROSSE E, BAZIN B. Self-organization during reactive fluid flow in a porous medium [J]. Geophysical Research Letters, 1998, 25: 385-388.

[15] SHERWOOD J D. Stability of a plane reaction front in a porous medium [J]. Chemical Engineering Science, 1987, 42: 1823-1829.

[16] ZHAO C, POULET T, REGENAUER-LIEB K. Replacement of annular domain with trapezoidal domain in computational modeling of nonaqaeous-phase-liquid dissolution front propagation problems [J]. Journal of Contral South University, 2015, 22: 1841-1846.

[17] VERMA A, PRUESS K. Thermohydrological conditions and silica redistribution near high-level nuclear wastes emplaced in saturated geological formations [J]. Journal of Geophysical Research, 1988, 93: 1159-1173.

[18] ZHAO C, HOBBS B E, HORNBY P, ORD A, PENG S, LIU L. Theoretical and numerical analyses of chemical-dissolution front instability in fluid-saturated porous rocks [J]. International Journal for Numerical and Analytical Methods in Geomechanics, 2008, 32: 1107-1130.

[19] ZHAO C. Physical and chemical dissolution front instability in porous media: Theoretical analyses and computational simulations [M]. Berlin: Springer, 2014: 354.

[20] ZHAO C. Advances in numerical algorithms and methods in computational geosciences with modeling characteristics of multiple physical and chemical processes [J]. Science China: Technological Sciences, 2015, 58: 783-795.

[21] ZIENKIEWICZ O C. The finite element method [M]. London: McGraw-Hill, 1977: 536.

(Edited by YANG Hua)

Cite this article as: ZHAO Chong-bin, Peter SCHAUBS, Bruce HOBBS. Effects of porosity heterogeneity on chemical dissolution-front instability in fluid-saturated rocks [J]. Journal of Central South University, 2017, 24(3): 720-725. DOI: 10.1007/s11771-017-3473-1.

Foundation item: Project(11272359) supported by the National Natural Science Foundation of China

Received date: 2015-10-25; Accepted date: 2016-01-20

Corresponding author: ZHAO Chong-bin, Professor; Tel: +86-731-88830028; E-mail: chongbin.zhao@iinet.net.au