Numerical Simulation of the Wormhole Propagation in Fractured Carbonate Rocks during Acidization Using a Thermal-Hydrologic-Mechanics-Chemical Coupled Model

Abstract

Acidizing is a widely adopted approach for stimulating carbonate reservoirs. The two-scale continuum (TSC) model is the most widely used model for simulating the reactive process in a carbonate reservoir during acidizing. In realistic cases, there are overburden pressure and pore pressure at present. When the injected acid reacts with the rock, the dissolution of the rock and the consumption of the acid in the pore will break the mechanical balance of the rock. Many experimental studies show that cores after acidizing have lower strength. However, it is still not clear how the deformation of rocks by the change of ground stress influences the acidizing dynamics. For fractured carbonate reservoirs, fractures play a leading role in the flow of injected acid, which preferentially flows into the fractures and dissolves the fracture walls. The effect of the combined action of rock mechanical balance broken and fracture wall dissolution on the formation of wormholes in fractured carbonate reservoirs remains to be studied. To address the above-mentioned issues, a thermal-hydrologic-mechanical-chemical coupled model is presented based on the TSC model for studying the wormhole propagation in fractured carbonate reservoirs under practical conditions. Linear and radial flow cases are simulated to investigate the influences of fracture distribution, reaction temperature, and effective stress on acidizing dynamics. The simulation results show that more wormhole branches are formed by acidizing if the fractures are perpendicular to the flow direction of acid. Temperature is a key parameter affecting the acidification dissolution patterns, so the influence of temperature cannot be ignored during the acidification design. As the effective stress of the formation increases, the diameter of the wormhole gradually decreases, and the branching decreases. More acid is needed for the same stimulation result under higher effective stress.

1. Introduction

Acidizing is one of the most common enhanced oil recovery (EOR) methods for oil and gas reservoirs. In matrix acidizing [1], the wormhole [2] penetrates the contaminated area near the wellbore, connecting the formation and the wellbore to increase production. There are typically five types of dissolution patterns for acidizing carbonate cores according to laboratory studies [3,4,5,6,7]. They are named face dissolution, conical dissolution, wormhole, ramified wormhole, and uniform dissolution. The patterns of dissolution depend on injection rates. Face dissolution takes place at a low injection rate, but uniform dissolution is formed when acid is injected at a high rate. Ramified wormholes and conical wormholes can be observed at other injection rates.

To understand the mechanism of wormhole propagation in carbonate reservoirs, extensive research has been conducted [3,4,8,9,10,11,12,13,14,15,16,17,18,19,20,21]. The two-scale continuum (TSC) model developed by Panga et al. (2005) [22] is the most commonly used one due to its computational efficiency. In the past decade, the TSC model has been continuously extended to better simulate wormhole propagation. Up to now, the TSC model can examine the influence of the reservoir conditions (temperature [23,24,25], pressure [25], domain geometry [26], rock properties (permeability [18], heterogeneity [18,27,28,29], anisotropy [27], minerals proportions [30], fracture geometry [31] and acid system characteristics (reaction kinetics and mass transfer [22], reaction rate and viscosity [32] on the wormhole formation. In addition, carbonate reservoirs usually have natural fractures. The injected acid penetrates the formation guided by the fractures, increasing the permeability of the formation. The productivity enhancement effect of acidization on fractured carbonate reservoirs is usually greater than that on nonfractured reservoirs [33,34,35]. Natural fractures have considerable effects on the acid distribution and wormhole propagation. Hence, it is important to study the acidizing process of naturally fractured carbonate reservoirs.

The models used to simulate the reactive flow process in fractured porous media can be divided into four categories. (1) the single fracture model [36,37,38,39,40,41,42,43]; (2) the fracture network model [44,45]; (3) the pseudo-fracture model [26,46,47] and (4) the discrete fracture model (DFM) [33,34,48,49,50,51,52]. The single fracture model only studies the reaction and flow in a single fracture. The fracture network model assumes that the acid flows only along the fractures while ignoring the fractures that are not directly connected to the wellbore. The DFM method treats the fractures explicitly and can characterize the effect of each fracture on fluid flow and solute transport. Due to the complexity of fracture geometry, an unstructured mesh is needed. However, when the distance or the angle between the adjacent fractures is very small, tiny mesh elements are required to discretize the area between fractures, leading to computational difficulties. The pseudo-fracture model treats fractures as high-porosity areas, which is simple in concept and can effectively model fractures in the dissolution process. Therefore, the pseudo-fracture model is adopted in the current study.

In realistic formation, the rock bears the internal fluid pressure in the pore and the overburdened rock pressure. When the injected acid reacts with the rock, the rock is constantly dissolved and the consumption of the acid in the pore will break the mechanical balance of the rock, resulting in the change of effective stress. When the effective stress changes, the pore structure [53] of the reservoir is deformed, which leads to the change of geophysical properties, hence the reservoir is stress sensitivity [54,55]. The stress sensitivity of fractured carbonate rocks is greater than that of non-fractured carbonate rocks and results mainly from that of fractures in the reservoir. Therefore, stress plays an important role in the acidification process of fractured carbonate rocks. It is necessary to study the formation mechanism of wormholes with formation pressure.

This paper establishes a Thermal-Hydrologic-mechanics-Chemical (THmC) coupled model to study the dissolution dynamics of acidizing fractured carbonate rocks. We use the capital letters T, H, and C to represent the full coupling of the thermal-hydrologic-chemical process, and lowercase (m) to indicate simplified mechanics. The model is focused on the Thermal-Hydrologic-Chemical processes and their coupling, with additional mechanical processes considered within a simplified geomechanically model. The main aim of this paper is to establish a comprehensive and practical simulation method for simulating the acidizing process in fractured carbonate reservoirs. The model takes into account the effect of reaction kinetics, rock heterogeneity, stress, and temperature effects on wormhole formation during acidization. Using realistic reservoir conditions, the dissolution process can be investigated by our model accurately.

This paper is organized as follows. In Section 2, the mathematical model is presented to describe the reaction involving the chemical system, fluid flow, mass and energy transport, rock dissolution, etc. In Section 3, the dimensionless mathematical model is derived. In Section 4 and Section 5, the numerical method is presented. In Section 6, a sensitivity analysis of fracture parameters, reaction exotherm, and effective stress is conducted. In the last section, the conclusions of the paper are summarized.

2. Mathematical Model

The assumptions regarding the thermal-hydrologic-mechanics-chemical coupled model in carbonate rocks are as follows. First, both fluid and solid materials are assumed to be incompressible, while the rock itself can still deform under stress due to a modification of the pore space. Second, the thermal equilibrium between fluid and rock can be achieved instantaneously. Third, ${\mathrm{CO}}_2$ generated from the acid reaction is dissolved in the liquid phase.

2.1. Chemical System

The main mineral component of the carbonate rock is ${\mathrm{CaCO}}_3$ and ${\mathrm{CaMg}(\mathrm{CO}}_3{)}_2$. Therefore, hydrochloric acid (HCl) [56] is appropriate for stimulating the formation. Hydrochloric acid can be completely ionized into hydrogen ions ${\mathrm{H}}^{+}$ and chloride ions (${\mathrm{Cl}}^{-}$) in aqueous solution, and the reaction between the hydrochloric acid and calcite, occurring at the fluid-solid surface can be expressed as:

$$2{\mathrm{H}}^{+}{+\mathrm{CO}}_3^{2-}\to {\mathrm{H}}_2\mathrm{O}+\mathrm{C}{\mathrm{O}}_2$$

(1)

2.2. Darcy Scale Model

The governing equations of Darcy flow in porous media are as follows:

$$\mathrm{v}=-\frac{\mathrm{K}}{\mu }·\nabla P$$

(2)

$$\frac{\partial \varphi }{\partial t}+\nabla ·\mathrm{v}=0$$

(3)

where $\mathrm{v}$ denotes the velocity, t is the time, $\mathrm{K}$ is the permeability tensor of the formation, P is the pressure of formation, $\varphi$ is the formation porosity and $\mu$ is the fluid viscosity. The material balance equation for the solute transport in the liquid phase is

$$\frac{\partial \left(\varphi C_f\right)}{\partial t}+\nabla ·\left(\mathrm{v}{C}_f\right)=\nabla ·(\varphi {\mathrm{D}}_e·\nabla C_f)-k_c{a}_{\mathrm{v}}(C_f-C_s)$$

(4)

where $C_f$ denotes the concentration of acid in the liquid phase, $C_s$ is the concentration of acid at the fluid–solid surface, ${\mathrm{D}}_e$ is the dispersion tensor, $k_c$ is the mass transfer coefficient and $a_v$ is the specific surface area. In general, a chemical reaction occurs on the liquid-solid surface and the amount of acid liquid per unit area transmitted to the liquid-solid surface is the same as that consumed by the surface reaction, expressed as:

$$k_c(C_f-C_s)=R\left(C_s\right)$$

(5)

where $R\left(C_s\right)$ denotes the reaction power, $R\left(C_s\right)$ is equal to $k_s·C_s$ when the chemical reaction rate is proportional to the concentration of the reactants. $k_s$ is the surface reaction rate constant that is dependent on temperature, calculated by the Arrhenius equation as:

$$k_s=A·exp(-\frac{E_a}{RT})$$

(6)

where A denotes the pre-exponential factor, $E_a$ is the activation energy and the value of both of them can be determined experimentally [57]; R is the universal gas constant; T is the temperature (in Kelvins). Let $k_s^0$ indicate the rate at temperature $T_0$. The rate at other temperatures can be calculated by the following formula

$$k_s=k_s^0·exp\left(-\frac{E_a}{R}\left(\frac{1}{T}-\frac{1}{T^0}\right)\right)$$

(7)

Porosity is changed by chemical reactions over time as:

$$\frac{\partial \left(\varphi \right)}{\partial t}=\frac{R\left(C_s\right)a_v\alpha }{\rho }$$

(8)

where ${\rho }_s$ denotes the density of solid and $\alpha$ is the power of the acid dissolving.

2.3. Pore Scale Model

Coefficient $k_c$ indicates the mass transfer rate from the liquid phase to the liquid-solid interface. Its value is related to the velocity and physical properties of the fluid along with the pore structure. The dimensionless form of the mass transfer coefficient is named the Sherwood number $Sh$ and is written as:

$$Sh=\frac{2k_c{r}_p}{D_m}=S{h}_{\infty }+0.7{Re}_p^{1/2}S{c}^{1/3}$$

(9)

where $r_p$ is the pore radius; $S{h}_{\infty }$ is the asymptotic Sherwood number; The value of $S{h}_{\infty }$ is taken as 3 in this article. ${Re}_p$ is the pore Reynolds number expressed as ${Re}_p=2r_{p}u/\nu$, where $\nu$ denotes the kinematic viscosity; $S_c$ is the Schmidt number defined as $S_c=\nu /D_m$ and $D_m$ is the molecular diffusivity.

The fractal dimension for pore spaces $D_f$ and tortuosity $D_T$ can be expressed respectively as:

$$D_f=d_E-\frac{ln\varphi }{ln\frac{{\lambda }_{min}}{{\lambda }_{max}}}$$

(10)

$$D_T=\frac{ln\left(1-\frac{\varphi }{2}+\frac{\sqrt{1-\varphi }}{4}+\frac{\left(\varphi +1+\sqrt{1-\varphi }\right).\sqrt{9-5\varphi -8\sqrt{1-\varphi }}}{8\varphi }\right)}{ln\left(\frac{D_f-1}{\sqrt{D{}_f}}\frac{{\lambda }_{max}}{{\lambda }_{min}}\sqrt{\frac{1-\varphi }{\varphi }\frac{\pi }{8-4D_f}}\right)}$$

(11)

where $d_E$ is the Euclidean dimension. ${\lambda }_{min}$ and ${\lambda }_{max}$ are the smallest and largest pore diameter, respectively. Sensitivity analysis shows that the value of ${\lambda }_{min}$ and ${\lambda }_{max}$ has little influence on the fractal dimension [28]. In this paper, the ratio of ${\lambda }_{max}/{\lambda }_{min}$ is set to be a typical value 0.01 following in Liu et al., 2017 [28]. Permeability is expressed by:

$$K=\frac{{\left(\pi D_f\right)}^{\frac{\left(1-D_T\right)}{2}}{\left(8-4D_f\right)}^{\frac{\left(1+D_T\right)}{2}}{d}^2}{128\left(3+D_T-D_f\right)}{\left(\frac{\varphi }{1-\varphi }\right)}^{\frac{3+D_T}{2}}$$

(12)

During the chemical reaction, the variation of permeability is correlated with the change in pore radius. The permeability, pore radius, and rock-acid interfacial surface area are expressed as follows:

$$\frac{K}{K_0}=\frac{{\left(\pi D_f\right)}^{\frac{\left(1-D_T\right)}{2}}{\left(8-4D_f\right)}^{\frac{\left(1+D_T\right)}{2}}\left(3+D_{T0}-D_{f0}\right)}{{\left(\pi D_{f0}\right)}^{\frac{\left(1-D_{T0}\right)}{2}}{\left(8-4D_{f0}\right)}^{\frac{\left(1+D_{T0}\right)}{2}}\left(3+D_T-D_f\right)}\frac{{\left(\frac{\varphi }{1-\varphi }\right)}^{\frac{3+D_T}{2}}}{{\left(\frac{{\varphi }_0}{1-{\varphi }_0}\right)}^{\frac{3+D_{T0}}{2}}}$$

(13)

$$\frac{r_p}{r_0}=\sqrt{\frac{K}{K_0}\frac{{\varphi }_0}{\varphi }}$$

(14)

$$\frac{a_v}{a_0}=\frac{\varphi }{{\varphi }_0}\frac{r_0}{r_p}$$

(15)

where $K_0$ indicates the initial permeability of the matrix, the subscript 0 is the initial value of the corresponding parameter, $\beta$ is the constant that relates with the structure of the medium, $r_0$ and $a_0$ are the initial values of average pore radius and interfacial area, respectively.

The diffusion coefficient following the exponential Arrhenius relation is given by:

$$D_m=D_m^{0}exp\left(-\frac{E_D}{R}\left(\frac{1}{T}-\frac{1}{T^0}\right)\right)$$

(16)

where $E_D$ indicates the activation energy of diffusion, and its value can be determined experimentally [58], varying from 12.6 to 28.1 kJ/mol. R is the universal gas constant.

The dispersion tensors depend on the pore structure, the pore-scale flow, and the properties of the fluid in the fluid flow. The relative ratio of convection to diffusion is expressed by the Peclet number, defined as:

$$P{e}_p=\frac{\left|\mathrm{v}\right|d_h}{\varphi D_m}$$

(17)

where $\left|\mathrm{v}\right|$ denotes the magnitude of the Darcy velocity, $d_h$ is the pore diameter defined as $d_h=2r_p$. When pore connectivity is good, Panga [22] gave the calculation formula of diffusion coefficient as:

$$D_{eX}={\alpha }_{os}{D}_m+D_m{\lambda }_{X}P{e}_p$$

(18)

$$D_{eT}={\alpha }_{os}{D}_m+D_m{\lambda }_{T}P{e}_p$$

(19)

where the subscripts X and T indicate acid injection direction and vertical direction respectively. The X denotes the x direction, and T is the y and z direction in the 3D orthogonal coordinate system. However, in 3D radial coordinate system, the X is the r direction and T is the $\theta$ and z direction; ${\alpha }_{os}$, ${\lambda }_X$ and ${\lambda }_T$ are the constants related to the pore structure. In this paper, we set typical values for them, which are 0.5, 0.5, and 0.1 respectively.

2.4. Energy Transport

First, the chemical reaction between HCl and carbonate rocks releases heat, resulting in a change in formation temperature. Second, the chemical reaction and acid transport rates are influenced by temperature. Third, there is energy transfer between the injected acid and the formation. The energy transport equation is written as:

$$\begin{array}{cc}\hfill & \frac{\partial }{\partial t}\left(\varphi {\rho }_f{C}_f^{P}T+(1-\varphi ){\rho }_s{C}_s^{P}T\right)+{\rho }_f{C}_f^P\nabla ·\left(\mathrm{v}T\right)\hfill \\ \hfill & =\left(\varphi {\mathrm{\Gamma }}_f+(1-\varphi ){\mathrm{\Gamma }}_s\right)\nabla ·(\nabla T)+2(-\Delta H_r\left(T\right))a_v{k}_e{C}_f\hfill \end{array}$$

(20)

where ${\rho }_f$ indicates the density of the liquid phase; ${\rho }_s$ is the density of the solid phase, $C_f^p$ is the specific heat capacity of the liquid phase; $C_s^p$ is the specific heat capacity of the solid phase. T is the temperature; ${\mathrm{\Gamma }}_f$ is the coefficient of heat conduction in the liquid phase, ${\mathrm{\Gamma }}_s$ is the coefficient of heat conduction in the solid phase; $\Delta H_r\left(T\right)$ is the heat released by reaction at temperature T, defined as the heat generated per mole of ${\mathrm{H}}^{+}$ consumed and expressed as:

$$\Delta H_r\left(T\right)=\sum |\Delta H_{resultants}\left(T\right)|-\sum |\Delta H_{reactants}\left(T\right)|$$

(21)

Equation (21) indicates the enthalpy difference between resultants and reactants. Kirchhoff’s law determines the value of enthalpy at different temperatures. At temperature T, the exothermic heat of reaction Equation (1) derived as follows:

$$\Delta H_r\left(T\right)=-6846+8.038T-3.22\times {10}^{-3}{T}^2-870.3T^{-1}$$

(22)

3. Dimensionless Model

The governing equations are rewritten into their dimensionless forms using the following parameters:

$$\begin{array}{c}{x}_D=\frac{x}{L},y_D=\frac{y}{L},z_D=\frac{z}{L},r_D=\frac{r}{L},r_{wD}=\frac{r_w}{L},\mathrm{U}=\frac{\mathrm{v}}{v_0},t_D=\frac{t}{L/v_0},r_{pD}=\frac{r_p}{r_0},a_{vD}=\frac{a_v}{a_0},\\ T_D=\frac{T}{T^o},\kappa =\frac{\mathrm{K}}{k_0},D_{mD}=\frac{D_m}{D_m^0},C_D=\frac{C_{}}{C_0},C_{fD}^P=\frac{C_f^P}{C_s^P},P_D=\frac{P-P_e}{\left(\mu v_{0}L\right)/K_0},h_T^2=\frac{2K_s\left(T^0\right)r_0}{D_m^0},\\ D_a=\frac{K_s\left(T^0\right)a_{v0}L}{u_0},N_{ac}=\frac{\alpha C_0}{{\rho }_s},{\rho }_{fD}=\frac{{\rho }_f}{{\rho }_s},A_E=\frac{E_a}{R{T}^0},D_{mD}=\frac{D_m}{D_m^0},a_{vD}=\frac{a_v}{a_{v0}},\\ P{e}_L=\frac{u_{0}L}{D_m^0},\eta =\frac{2r_0}{L},H_T^2=\frac{K_s\left(T^0\right)a_{v0}{L}^2}{D_m^0},P{e}_H=\frac{L{u}_0{\rho }_s{C}_s^p}{f{\mathrm{\Gamma }}_{ef}+(1-f){\mathrm{\Gamma }}_{es}},H_D=\frac{-\Delta H_r{c}_0}{T_0{\rho }_s{C}_s^p}.\end{array}$$

where subscript D indicates that the variable is dimensionless; L is the characteristic length set to be the core length for linear flow and the core radius for radial flow; $\mathrm{U}$ is the dimensionless velocity; $\kappa$ is the dimensionless permeability. $\eta$ is the ratio of the pore-to-domain scale. $A_E$ is the dimensionless activation energy; $H_D$ is the dimensionless exothermic reaction; $H_T^2$ and $h_T^2$ are Thiele modulus, which is the reaction velocity divided by the diffusion velocity at the initial pore volume. $Da$ is the Damkohler number defined as the reaction rate divided by the convection velocity; $P_{eL}$ and $P_{eH}$ are the Peclet number defined on a core scale and heat transport scale, respectively. However, $P_{eL}$ is the convective transport rate divided by the diffusive transport rate. $P_{eH}$ is the ratio of heat transfer by convection to that by thermal conduction. The acid capacity constant $N_{ac}$ is defined as the volume of dissolved solids per unit volume of consumed acid. Therefore, the governing equations can be rewritten as:

$$\mathrm{U}=-\kappa ·\nabla P_D$$

(23)

$$\frac{\partial \varphi }{\partial t_D}+\nabla ·\left(-\kappa ·\nabla P_D\right)=0$$

(24)

$$\frac{\partial \left(\varphi C_D\right)}{\partial t_D}+\nabla ·\left(\mathrm{U}{C}_D\right)=\nabla ·\left({\mathrm{D}}_{eD}\nabla C_D\right)-\frac{Sh·D_{mD}·Da·a_{vD}·C_D·exp\left(A_E(1-1/T_D)\right)}{Sh·D_{mD}+h_T^2{r}_{pD}·exp\left(A_E(1-1/T_D)\right)}$$

(25)

$$\frac{\partial \varphi }{\partial t_D}=\frac{Sh·Da·a_{vD}·C_D·D_{mD}·N_{ac}·exp\left(A_E(1-1/T_D)\right)}{Sh·D_{mD}+h_T^2{r}_{pD}·exp\left(A_E(1-1/T_D)\right)}$$

(26)

$$\begin{array}{cc}& \frac{\partial T_D}{\partial t_D}\left(\varphi {\rho }_{fD}{C}_{fD}^P+(1-\varphi )\right)+{\rho }_{fD}{C}_{fD}^P\nabla ·\left(\mathrm{U}{T}_D\right)\hfill \\ & =\frac{1}{P{e}_H}\nabla ·(\nabla T_D)+\frac{H_D·Sh·D_{mD}·Da·a_{vD}·C_D·exp\left(A_E\left(1-\frac{1}{T_D}\right)\right)}{Sh·D_{mD}+h_T^2{r}_{pD}·exp\left(A_E\left(1-\frac{1}{T_D}\right)\right)}\hfill \end{array}$$

(27)

$$D_{eD}^X=\frac{\varphi {\alpha }_{os}Da{D}_{mD}}{H_T^2}+{\lambda }_X\left|\mathrm{U}\right|r_{pD}·\eta$$

(28)

$$D_{eD}^Y=\frac{\varphi {\alpha }_{os}Da{D}_{mD}}{H_T^2}+{\lambda }_Y\left|\mathrm{U}\right|r_{pD}·\eta$$

(29)

The dimensionless form of the initial and boundary conditions can be written as:

$$-\kappa ·\nabla P_D·\mathrm{n}=1,\mathrm{on}\partial {\Omega }_{in}$$

(30)

$$C_D-\left(\frac{\varphi {\alpha }_{os}Da{D}_{mD}}{H_T^2}+{\lambda }_X\left|U\right|r_{pD}·\eta \right)\nabla C_D·\mathrm{n}=1,\mathrm{on}\partial {\Omega }_{in}$$

(31)

$$T_D=T_{fD},\mathrm{on}\partial {\Omega }_{in}$$

(32)

$$P_D=0,\mathrm{on}\partial {\Omega }_{out}$$

(33)

$$\nabla c_D·\mathrm{n}=0,\mathrm{on}\partial {\Omega }_{out},\mathrm{and}\partial {\Omega }_{tran}$$

(34)

$$\nabla T_D·\mathrm{n}=0,\mathrm{on}\partial {\Omega }_{out},\mathrm{and}\partial {\Omega }_{tran}$$

(35)

$$\nabla P_D·\mathrm{n}=0,\mathrm{on}\partial {\Omega }_{tran}$$

(36)

$$c_D=0,T_D=T_{sD},\mathrm{at},t=0$$

(37)

4. Numerical Method

The governing equations are solved sequentially and described as follows. First, the permeability, mass transfer coefficient, and diffusion coefficient are calculated based on initial conditions. Then, they are brought into Equation (24) to solve the pressure distribution at $t^{n+1}$. Next, the pressure distribution is brought in Equation (23) to obtain the Darcy velocity. The reactive transport procedure is divided into two parts, i.e., a physical transport process and a process involving chemical reactions. In the physical transport process, the distribution of temperature and acid concentration is calculated by physical transport, followed by their values modified by reactions, as shown below:

$$\left\{\begin{array}{cc}\hfill & \frac{\partial T_D}{\partial t_D}\left(\varphi {\rho }_{fD}{C}_{fD}^P+(1-\varphi )\right)+{\rho }_{fD}{C}_{fD}^P\nabla ·\left(\mathrm{U}{T}_D\right)=\frac{1}{P{e}_H}\nabla ·(\nabla T_D)\hfill \\ \hfill & \frac{\partial \left(\varphi C_D\right)}{\partial t_D}+\nabla ·\left(\mathrm{U}{C}_D\right)=\nabla ·\left({\mathrm{D}}_{eD}·\nabla C_D\right)\hfill \\ \hfill & \frac{\partial \varphi }{\partial t_D}=0\hfill \end{array}\right.$$

(38)

$$\left\{\begin{array}{cc}\hfill & \frac{\partial }{\partial t_D}\left[\left(\varphi {\rho }_{fD}{C}_{fD}^P+(1-\varphi )\right)T_D-\frac{H_D\varphi }{N_{ac}}\right]=0\hfill \\ \hfill & \frac{\partial }{\partial t_D}\left(\varphi C_D+\frac{\varphi }{N_{ac}}\right)=0\hfill \\ \hfill & \frac{\partial \varphi }{\partial t_D}=N_{ac}·\frac{Sh·Da·a_{vD}·C_D·D_{mD}·exp\left(A_E(1-1/T_D)\right)}{Sh·D_{mD}+h_T^2{r}_{pD}·exp\left(A_E(1-1/T_D)\right)}\hfill \end{array}\right.$$

(39)

The finite volume method is employed to discretize these governing equations. Orthogonal and radial grids are used for the linear and radial cases, respectively. The diffusion terms are discretized by the central difference scheme and the TVD scheme is used to discretize the convection terms. The backward Euler scheme is used for time integration to help the convergence of solutions at larger time steps. For the linear case, the simulation is performed on a $5\times 2$ cm rectangle core. This rectangular domain is discretized with a $500\times 200$ orthogonal grid. For the radial case, the simulation is performed on an annular domain, which is discretized with a $160\times 640$ radial grid system (160 cells in r direction and 640 cells in $\theta$ direction). The fractures are introduced by the following procedure:

(1)

Generate randomly distributed discrete fractures.

(2)

Calculate the starting point and end point coordinates of each fracture.

(3)

Loop all faces of the grid cell to check whether the fracture segment intersects the face segment.

(4)

If the fracture intersects the face segment of the grid cell, set the two cells sharing the face as the fracture cell, and change its porosity to 0.99.

The porosity value 0.99 used to characterize fractures is obtained by comparing the Darcy formula with the laminar flow formula, through which the cubic law can be derived as:

$$K=\frac{b^2}{12}$$

(40)

In the current work, we set the fracture aperture b = 1 mm, which means that the permeability of fracture is $8\times {10}^7$ md. For a porous medium with an average porosity of 0.35 and an average permeability of 100 md, it can be calculated with Equation (13) that the porosity of fracture is 0.99.

5. Effect of Stress

Prior to the development of oil and gas reservoirs, the pressure of overburden formation is balanced by the pressure of pore fluid and the support of the rock skeleton. With the dissolution of rocks, the porosity increases consequentially. If the influence of stress is ignored, the variation of porosity can be obtained through the procedure described in Section 4. However, stress exists in the actual formation. The influence of stress should be considered on the basis of the original porosity increment induced by rock dissolution. The influence of stress on porosity increment can be expressed as follows:

$$\Delta \varphi ={\Delta \varphi }_0{e}^{-a_0\sigma }$$

(41)

where $\Delta \varphi$ indicates the increment of porosity under the current effective formation stress; $\Delta {\varphi }_0$ is the porosity increment induced by rock dissolution; $a_0$ is the pore compression coefficient of rock at the original pressure and its value is 0.08; $\sigma$ is the effective stress of rock and its value can be calculated by

$$\sigma =P_0-P$$

(42)

where $P_0$ is the pressure of the primordial formation; P is the current pore pressure of the rock.

6. Results and Discussions

This section presents simulation results for the sensitivity analysis of fracture parameters, temperature, and stress. The values of the parameters are as in

Table 1. List of values of parameters used in the simulation.

Parameters	Values
L	5 cm
H	2 cm
$r_e$	10 cm
$r_w$	0.5 cm
${\varphi }_0$	0.35
$K_0$	100 md
$r_0$	$5\times {10}^{-6}$
$k_s$	0.2
$\upsilon$	0.0286
$S{h}_{\infty }$	3
$Sc$	1000
$N_{ac}$	0.1
$\eta$	$1\times {10}^{-6}$
$v_0$	0.004 cms
${\alpha }_{os}$	0.5
${\lambda }_T$	0.1
${\lambda }_X$	0.5
m	1
$\beta$	1
$h_T^2$	0.07

unless otherwise stated.

6.1. Effect of Fracture Parameters

6.1.1. Effect of Fracture Orientation

The effects of fracture orientation on the dissolution structure are studied from both radial and linear aspects (Figure 1), while keeping the other parameters fixed. It is evident in Figure 1 that the wormhole is formed along the fractures which become part of the wormhole. However, the direction of the wormhole is not along fractures if the fractures are perpendicular to the injection direction (Figure 1 row(c), column (A)). The calculated $P{V}_{BT}$ is plotted against fracture inclinations (Figure 2). It can be seen from Figure 2 that as the fracture inclination increases, the breakthrough pore volume increases, but then decreases. When the dip angle between the fracture and the flow direction is small, the fracture becomes part of the wormhole, and the breakthrough can be achieved by dissolving a small amount of rock. When fractures are perpendicular to the flow direction, wormholes spread along both the fracture direction and fluid flow direction. Furthermore, the presence of the fractures makes the wormholes produce more branches, which means more acid is needed to dissolve the rock. So, the breakthrough pore volume reaches its maximum when the fracture inclination is ${90}^{\circ }$. However, fracture orientation has little effect on wormhole structures in radial acidizing simulation as shown in Figure 1 column (B).

6.1.2. Effect of Fracture Length

The distribution of natural fractures is an important factor for acidization Figure 3 and Figure 4 show the influence of fracture length on dissolution structure. $L_L$ and $L_R$ are the average fracture lengths of linear and radial acidizing simulations, respectively. The values are placed under the picture. Related parameters are shown in Table 1. In Figure 3A and Figure 4A, the wormholes are more irregular, wider, and more branched. Obviously, as the fracture length is increased, wormholes become thinner with fewer branches (Figure 3 and Figure 4). When the fractures are longer, they can be interconnected to allow acid advances faster through fractures and prevent the formation of ramified wormholes. The relationship is shown between $P{V}_{BT}$ against fracture length in Figure 5. The $P{V}_{BT}$ for the longest fractures is the smallest, demonstrating that longer fractures can significantly improve the efficiency of wormhole propagation.

6.1.3. Effect of Fracture Density

The density of natural fractures has a significant influence on wormhole propagation. The reactive flow simulation of cores with a different number of fractures was carried out for both the linear and radial cases (Figure 6 and Figure 7). It is obvious that when the fracture density is small, wormholes are created with more branches. As the increase of the number of fractures, competition appears between fractures. When the number of fractures is large, the injected acid mainly flows along the fractures. This is attributable to the fact that the permeability of the fracture is much larger than that of the matrix. The injected acid preferentially flows along the path with the least resistance during wormhole propagation, resulting in a smaller breakthrough volume of injected acid and much higher efficiency in wormhole formation.

6.2. Effect of Reaction Temperature

In this section, the acid is injected into the rock radially or linearly at the same formation temperature (375 K). The effect of temperature on the dissolution pattern is investigated. The acid reaction with carbonate minerals releases heat and increases the temperature of the formation, which in its turn accelerates the reaction. To investigate the effect of this heating on the dissolution pattern, in this section, the formation temperature is set to a constant value of 375 K. The dissolution structures are shown in Figure 8 and Figure 9. It is obvious that different wormhole patterns are generated at the breakthrough stage under the influence of temperature. This is because acid reaction becomes faster at a higher temperature. The diffusivity and reactivity also increase with higher temperatures. Hence, different dissolution modes are produced at the same injection rate.

6.3. Effect of Stress Sensitivity

Stress sensitivity is a common characteristic in oil and gas reservoirs. Many theoretical, numerical, and experimental studies have been conducted to study the effect of stress sensitivity [59,60,61,62]. The distribution, variation, and stress sensitivity of fractures along with the hydraulic properties of the matrix are regarded as important parameters for thermal–hydraulic–mechanical modeling [63,64,65].

Although many studies have been carried out on the stress sensitivity of fractured reservoirs, they mainly focus on the relationship between reservoir permeability and stress sensitivity but neglect the effect on the dissolution patterns. The stress sensibility of the reservoir results mainly from the presence of fractures. The influence of the effective stress on wormhole patterns is analyzed in fractured carbonate rocks from the radial and linear aspects. This is evident in Figure 10 and Figure 11 that the range of the wormhole propagation is reduced at low injection rates. While at high injection rates, the range of the wormhole propagation is extended (Figure 10d). It is worth noting that at a low injection rate, the change of effective stress has little effect on wormhole structures in Figure 11d. The rapid injection and incomplete reaction allow some reactants to propagate further downstream, resulting in a more uniform dissolution. The effective stress has little effect on the dissolution patterns of wormhole in linear simulation (Figure 10). However, the dissolution pattern is changed by the conical dissolution of the wormhole in radial simulation (Figure 11 row (b)). The reason is that the value of the formation porosity decreases as the effective stress increases. The value of $P{V}_{BT}$ is plotted against $1/Da$ in Figure 12. It could be observed in Figure 12 that $P{V}_{BT}$ increases with higher effective stress regardless of the injection rate.

7. Conclusions

In this study, we developed a Thermal-Hydrologic-mechanics-Chemical coupled model to simulate the dissolution process during acidizing in carbonate reservoirs. The following conclusions are made after analyzing the influence of fracture parameters, temperature, and stress on the dissolution process.

1. Wormhole propagation is affected by the orientation of fractures. The wormhole grows along the fracture direction in both linear and radial cases. When acid is injected into the formation radially, no matter which direction the fracture is, the injected acid can always flow along the fracture direction. Therefore, the direction of wormhole propagation is determined by the fracture distribution of the formation.

2. With the increase in fracture length, thinner wormholes with fewer branches are obtained. If there is a long fracture in the formation, the wormhole will expand along this fracture. Therefore, to make all zones near the wellbore be stimulated sufficiently, long fractures need to be plugged before the acid is injected into the formation.

3. The presence of natural fractures enhanced the efficiency of acidizing. As increasing natural-fracture density, wormholes are formed with fewer branches. The injected acid preferentially selectively flows into fractures with the least resistance, resulting in much higher efficiency in wormhole propagation.

4. The influence of temperature on the dissolution pattern cannot be ignored, and hence should be considered in the acidizing design. With the increase in temperature, not only the acid reaction becomes faster, but also the diffusivity and reactivity also increase. Hence, different dissolution patterns are produced at the same injection rate by Comparing the non-isothermal and isothermal dissolution structures.

5. As the effective stress increases, the increment of porosity decreases, and the permeability enhancement is inhibited, resulting in the increase of $P{V}_{BT}$. Therefore, the mechanical compaction effects should be considered when designing the acidizing stimulation treatment.

It should be noted that in this work we are mainly concerned with coupling in Thermal-Hydrologic-Chemical(THC) processes and only local effects are considered for the mechanics, so only limited simulation result is presented. In practice, the mechanics should be coupled to THC processes fully. A fully coupled Thermal-Hydrologic-Mechanics-Chemical model for simulating the wormhole formation during acidizing will be developed in our future work.