A theoretical and numerical framework for modeling gypsum cavity dissolution

Abstract

In this paper, a large-scale diffuse interface model is used to describe the evolution of a gypsum cavity formation induced by dissolution. The method is based upon the assumption of a pseudo-component dissolving with a thermodynamic equilibrium boundary condition. A methodology is proposed based on numerical computations with fixed boundaries in order to choose suitable parameters for the diffuse interface model, and hence predict the correct dissolution fluxes and surface recession velocity. Additional simulations were performed to check which type of momentum balance equation should be used. The numerical results did not show a strong impact of this choice for the typical initial boundary value problems under consideration. Calculations with a variable density and Boussinesq approximation were also performed to evaluate the potential for natural convection. The results showed that the impact of density driven flows was negligible in the cases under investigation. The potential of the methodology is illustrated on two large-scale configurations: one corresponding to a gypsum lens located strictly within a porous rock formation and the other to an isolated pillar in a flooded gypsum room and pillar quarry. boundary condition for the small-scale dissolution problem) for limestone, calcite, gypsum, or salt follow a similar form expressed aswhere k s is the surface reaction rate coefficient, c s is the mass concentration of the dissolved species at the surface, and c eq the corresponding equilibrium concentration or solubility. It is worth noticing that the reaction order, n, can lead to a highly nonlinear reaction rate. For instance, it ranges from 1 to 4.5 approximately in the case of gypsum dissolution, depending on c s , as obtained by Jeschke et al. [21]. Also, if the surface Damköhler number is very large, for instance a very large k s , this boundary condition tends to the classical equilibrium condition c s D c eq at the solid surface. This latter condition is often used for salt dissolution, for instance. Dissolution models will be based on the coupling of the transport equations and the interface recession. From a numerical viewpoint, some models have been developed to simulate the early stage of conduit evolution, that is, the dissolution enlargement of preexisting fractures in rock matrix, mainly for hydrology analyses. Early studies used 1D pipe models to analyze the single conduit development in limestone [24][25][26]. A main conclusion of these studies was that the positive feedback of mutual enhancement of conduit enlargement and flow rate governed the early stage of karst genesis. 2D pipe network models representing interconnected conduits were developed in [27][28][29] for more complex limestone structures, and flow types were extended from laminar to turbulent. It was reported in these works that both flow patterns and hydraulic boundary conditions have significant impact on cave patterns, with laminar flow favoring the development of single passages, while turbulent flow triggering the developme...

Modeling of salt and gypsum dissolution: applications, evaluation of geomechanical hazards

“…Natural geological systems are often highly heterogeneous, involving multiple scales ranging from pore-scale of millimeter size, to large-scale of tens of meters (Guo, Quintard, & Laouafa, 2015;Guo, Laouafa, & Quintard, 2016), as schematically illustrated in Figure 2. Consequently, the numerical simulation of subsurface dissolution of porous media or cavities can be conducted at different scales.…”

“…Several important effective parameters may be introduced in the macro-scale models to link the averaged behaviour with the pore-scale features, and they are obtained by solving the so-called closure problems developed in the upscaling process. Among various effective parameters, the effective dispersion tensor (Guo et effective mass exchange coefficient (Coutelieris, Kainourgiakis, Stubos, Kikkinides, & Yortsos, 2006;Guo et al, 2015Guo et al, , 2016 are of great significance. In addition, non-classical parameters may also have some influences on the accuracy of the macro-scale models, as discussed in Guo et al (2015).…”

“…In addition, non-classical parameters may also have some influences on the accuracy of the macro-scale models, as discussed in Guo et al (2015). Darcy-scale non-equilibrium models have been successfully used in Guo et al (2016), Soulaine et al (2017), Zhao, Hobbs, Ord, and Peng (2010) to describe the moving front between fluid and porous domains, with the front represented by a transient zone characterised by a porosity gradient.…”

“…When the effective mass exchange coefficient is large enough, the front tends to be very thin and local equilibrium condition is recovered. The ability of DIM to capture the recession of solid formations, even coupled with hydrodynamic instabilities such as natural convection, has been demonstrated in several works (Luo et al, 2012(Luo et al, , 2014Guo et al, 2016).…”

“…In this article, we are interested in the potential to use DIM in conditions corresponding to real industrial situations. Since the framework to develop the volume averaged equations and to use them for the modeling of given soluble rocks has been provided in several previous works (Guo et al, 2015(Guo et al, , 2016Luo et al, 2014), this question will not be detailed in the article which is organised as follows. We will first go briefly through the theoretical part for reader's comprehension.…”