Numerical Study of Reactive Flow in Fractured Carbonate Rock

Acidizing technology is an effective reformation method of oil and gas reservoirs. It can also remove the reservoir pollution near wellbore zones and enhance the fluid transmissibility. The optimal injection rate of acid is one of the key factors to reduce cost and improve the effect of acidizing. Therefore, the key issue is to find the optimal injection rate during acid corrosion in fractured carbonate rock. In this work, a novel reactive flow mathematical model based on two-scale model and discrete fracture model is established for fractured carbonate reservoirs. The matrix and fracture are described by a two-scale model and a discrete fracture model, respectively. Firstly, the two-scale model for matrix is combined with the discrete fracture model. Then, an efficient numerical scheme based on the finite element method is implemented to solve the corresponding dimensionless equations. Finally, several important aspects, such as the influence of the injection rate of acid on the dissolution patterns, the influence of fracture aperture and fracture orientations on the dissolution structure, the breakthrough volume of injected acid, and the dynamic change of fracture aperture during acidizing, are analyzed. The numerical simulation results show that there is an optimal injection rate in fractured carbonate rock. However, the fractures do not have an impact on the optimal acid injection rate, they only have an impact on the dissolution structure.


INTRODUCTION
Since the 1980s, many scholars have conducted systematic research on matrix acidizing. Hoefner (Hoefner and Fogler, 1989;Fredd et al., 1997) and (Fredd et al., 1997) first undertook studies on the dissolution patterns of porous media. They injected the inorganic acid into limestone, then, injected the low melting point alloy into limestone after acid etching. They studied the dissolution structure of porous media by observing the shape of the alloy. (Hoefner and Fogler, 1988) used the network model to study the wormhole propagation and formation in the porous media and performed a series of experiments to analyze the mechanism of wormhole formation and many numerical simulations to study the characteristics of the mass-transfer limited regime and reaction-limited regime. These results suggested that the branch length of the main channel will never exceed the distance between the branch and the main channel. Next, based on this mathematical model, (Budek and Szymczak, 2012) improved the extended pore network model to qualitatively characterize the dissolution patterns at different Damköhler numbers and analyze optimal injection velocity. (Kim and Santamarina, 2015) explored how CO 2 dissolves into water and flows into the reservoir with a 2-D pore network model. (Wang et al., 1993) adopted the core displacement method to study the influence of temperature, acid concentration, and the velocity of the injected acid on the wormhole. Finally, (Frick et al., 1994) also found the optimal injection rate at different acid concentrations and temperatures in the radial core. These analysis result consistent with their mathematical model, which can predict the wormhole diameter and the permeability influences the breakthrough time in the cross-section of the core. (Bazin et al., 1995) proposed a pressure drop function based on morphology after acidizing. (Daccord et al., 1989) used water flowing on gypsum to simulate wormhole formation, three dissolution patterns were obtained: compaction, wormhole, and uniform. (Bazin, 2001) studied the optimal injection rate of injected acid using CT scanning technology and (Ziauddin and Bize, 2007) discusses the influence of the heterogeneity of porous media on the dissolution patterns based on CT scanning technology, NMR, and SEM technologies. (Zhang et al., 2017) studied the dissolution structure of acid-etched core by CT scanning technology, undertaking (He, 2009) an acid erosion experiment on carbonate rocks with undeveloped fractures, but the permeability is so low that the injected acid struggles to form a wormhole and break through the core end. (Yang and Pan, 2000) studied the influence of the injection rate and type of acid on dolomite based on the mechanism of acid filtration. (Daccord et al., 1989) described the dissolution structure based on fractal theory. (Golfier et al., 2001) studied the reactive flow combined with Brinkman. They found that the optimal condition was related to the acid concentration and the length of the core. This was the first confirmation that the acid capacity number affects the breakthrough time of injected acid.
A.D.Hill (Hill et al., 2009) built a model of wormhole formation considering acid filtration theory. (Liu et al., 1997) developed a simulator to describe the dissolution patterns. (Chen and Ying, 2006) deduced the equations of acid diffusion and acid surface reaction rate of reactive flow. (Catherine, 2004) injected the CO2-enriched water into the limestone to study the relationship between porosity and permeability, which are distinct at different dissolution stages in the core. (Andersen and Evje, 2016) considered flow through a fracture coupled with diffusion to the surrounding matrix, where the reaction occurring presented a reactive flow model in the fractured medium. (Nierode and Williams, 1971) established description equation for reactive flow based on the reaction kinetics, indicating that the mass transfer coefficient of solute in acid affects the result of acidizing treatment, and the effect of matrix acidizing was related to the injection rate of acid in the ground.
The above physical experiments and numerical simulation experiments did not consider the influence of natural complex fractures on the dissolution patterns and the dynamic change of fracture aperture during acid-rock reaction time. In response, this paper established a novel mathematical model of reactive flow in a fractured medium based on the two-scale model and the discrete fracture model. The numerical simulation of reactive flow in the fractured medium was realized.
In The Mathematical Model for Matrix of this paper, the twoscale mathematical model in matrix is presented. The Mathematical Model for Fracture outlines the mathematical model of reactive flow in a fractured medium based on the discrete fracture model. Then, in Dimensionless the governing equations and boundary conditions are written in a dimensionless formula. Validation of the Model verifies the theoretical model by examining the conclusions of previous studies. In Numerical Results and Discussion, the numerical simulation of reactive flow is studied in different conditions. Finally, the paper is summarized by conclusions in Conclusion.

Darcy Scale Model
The flow equation of injected acid in a porous medium is described by Darcy's law, the formula is written as: Where, v is the vector of Darcy's velocity, m/s; K is the tensor of permeability in study region, m 2 ; μ is fluid viscosity, mPa·s; P is fluid pressure, Pa. The continuity equation is derived from the law of mass conservation, the formula is as follows: Where, ϕ is the porosity of carbonate rock; t is reaction time, s. Q mf is the normal flow rate through the interface between the matrix system and the fracture system. Generally, when the acid is injected into the carbonate rock, the equation of acid transport in the rock is described by the following mass balance equation.
Where, C f is acid concentration in a liquid phase, mol/m 3 ; D e is the tensor of the diffusion coefficient, m/s 2 ; Φ mf is solute mass transmitted into the matrix system through the interface between matrix medium and fracture medium. R is reaction term, its physical significance is that the amount of acid transferred to the liquid-solid surface is equal to the consumption of surface reaction, the formula of reaction term is expressed as: Where, a v is specific surface area, m −1 ; k c is local transport coefficient of acid, m/s; C s is the acid concentration of liquidsolid surface in core pore, mol/m 3 , it is calculated by: The porosity will change as acid corrodes carbonate rock. The calculation formula of dynamic change of porosity is as follows: zϕ zt Where, α is the dissolving power of acid, kg/mol; ρ s is the density of carbonate rock, kg/m 3 .

Pore Scale Model
The porosity, radius of the pore, and specific surface area of rock are decided by pore structure. It is a dynamic process and the injected acid changes the pore structure of rock when it continuously erodes porous media. This means that the different kinds of physical parameters of carbonate rock will change dynamically. There are two methods to compute the physical parameters of carbonate rock in reactive flow (Panga et al., 2005;Kalia and Balakotaiah, 2009;Maheshwari et al., 2013;Ghommem et al., 2015;Liu et al., 2020). In the first method, these various physical parameters are computed based on the dynamic change of pore structure obtained through the laboratory. Another method is to calculate the various physical parameters of rock according to empirical or semi-empirical formulas. This paper adopts the modified Carman Kozeny equation to describe the dynamic change of porosity and permeability, the formulas are as follows: Where, r p is pore radius, m; a 0 is initial specific surface area, m-1; k 0 is initial permeability, m 2 ; ϕ 0 is initial porosity; β is a constant that depends on the structure of the medium. When solute of liquid phase flows in the pores, the velocity of solute in the acid through the pore flows from liquid phase to pore surface and contact with the pore surface is expressed by the rate of the mass transfer. It can be seen from Eq. 1.5 that the effect of the rate of mass transfer on reactive flow cannot be ignored when the acid system is certain, because its magnitude makes a difference to the chemical reaction. (Panga et al., 2010) defined a dimensionless number to the analysis of the associated factors, it is referred to the Sherwood number and represents the dimensionless mass transfer coefficient, it is given by: Where, r p is the average pore radius of porous media, m; D m is the molecular diffusion coefficient, m/s 2 ; Sh ∞ is asymptotic Sherwood number of the pore, b is a constant number related to the structure of porous media, b 0.7/m 1/2 , where, the m is the ratio of pore length to pore radius, Re p is pore Reynolds number, Re p 2ur p /], where, ] is hydrodynamic viscosity; Sc is Schmidt number, Sc ] k /D m , ] k is hydrodynamic viscosity.

Diffusion Coefficient
When the acid is transported in the homogeneous and anisotropic matrix, in 2D, the diffusion tensor is characterized by axial diffusion coefficient D eX and transverse diffusion coefficient D eT . There is only solute molecular diffusion in the liquid phase and the transverse diffusion coefficient equals the axial diffusion coefficient when the fluid does not flow, the formula is as follows: Where, D m is the molecular diffusion coefficient, m/s 2 ; α os is a constant related to the structure of porous media (such as tortuosity or connectivity between pores). The twodimensional diffusion coefficient of acid relates to the geometry of porous media, the flow pattern of pore scale, and the properties of acid. Researchers usually define a dimensionless Peclet number to describe the contrast between convection and diffusion. The expression formula of the Peclet number is found to be Where, |v| describes the magnitude of Darcy's velocity, m/s; d h is pore diameter of porous media, m. The formula of diffusion coefficient derived by Panga (Panga et al., 2010) is used to describe the diffusion of acid in porous media in this paper. the diffusion coefficient can be computed as: Where, the subscript of X and T denotes the injection direction and vertical transverse direction of acid injection; α os , λ T, and λ x are constants that have typical values of 0.5, 0.5,0.1 for a packedbed of spheres (Kalia, 2008;Kalia and Balakotaiah, 2009;Liu et al., 2017), respectively.

Initial Conditions and Boundary Conditions
The pressure at the outlet is constant and the outlet boundary is a constant pressure boundary. The upper and down boundaries are closed boundaries. The specific expression of boundary conditions is given by: Where, u 0 is initial injection velocity, m/s; C 0 is a concentration of injected acid at the inlet, mol/m 3 ; P e is a constant, it is boundary pressure at the outlet, Pa.
A random function f is defined to simulate the heterogeneity of rock. Its value varies from -Δϕ 0 to +Δϕ 0 and satisfies the random function of uniform distribution, the generated initial porosity field is shown in Figure 1. The heterogeneity of fracture aperture is defined in this way.

THE MATHEMATICAL MODEL FOR FRACTURE
The governing equation of the fracture system is written as: Where, v f is the velocity vector of injected acid in the fracture system, the subscript f denotes fracture, d is fracture aperture, ϕ f is fracture porosity; Q mf is the flow exchange between the matrix system and fracture system. According to the cubic law, the flow velocity of injected acid in fracture system is expressed by the Poiseuille equation, the expression equation of flow velocity is given as: The expression equation of permeability in fracture system is derived as: Similarly, the convection-diffusion equation in fracture system is obtained as: Where, C f is the concentration of injected acid in the fracture system, D e,f is the effective diffusion coefficient of solute in acid in the fracture system. Φ mf is solute mass transmitted into the fracture system through the interface between matrix medium and fracture medium. Fracture aperture changes when the acid flow in the fracture system alters. The dynamic change process of fracture aperture can be obtained as: The effective diffusion coefficient (Steefel and Lichtner, 1998) of solute in acid in fracture system is as follows: Where α f is constant, which is equal to 0.5 in this paper.

DIMENSIONLESS
Many factors affect the process of reactive flow. The above governing equations and boundary conditions are written in a dimensionless formula to make the numerical solution of reactive flow easier and numerical simulation results more clear. The expression equations are defined: Where the superscript p denotes dimensionless; L is a characteristic length of the research area of the mathematical model; x,y are parameters of the coordinate system; U p indicates that velocity vector is dimensionless; K p indicates that the permeability vector is dimensionless; k 0 is initial permeability of porous media; r 0 is average pore radius of porous media; a 0 is initial specific surface area, B indicates that fracture aperture is dimensionless; the subscript f denotes fracture; D f p indicates diffusion coefficient is dimensionless; K p indicates that permeability vector of fracture is dimensionless. Da represents the rate of reaction velocity to convection velocity at the core scale.
Based on the expression of the above-mentioned variables, dimensionless variables are substituted into governing equation. The dimensionless equations are obtained: The dimensionless equation of porosity changing is Frontiers in Earth Science | www.frontiersin.org November 2021 | Volume 9 | Article 765139 The dimensionless equation of diffusion coefficient is rewritten as: The above equations are mathematical models of matrix systems for reactive flow. The dimensionless mathematical equations in the fracture system for reactive flow are The dimensionless boundary conditions and initial conditions are arranged as:

VALIDATION OF THE MODEL
In this section, the correctness of the mathematical model of the matrix is verified. The numerical simulation results are compared with (Fredd and Fogler, 1999) physical experiment results. In the physical experiment, the diameter of the core is 3.8 cm and the length of the core is 10.2 cm. The comparison results of numerical simulation and physical experiment are shown in Figure 2.
Through the comparison results, there are five different dissolution patterns. The results of numerical simulation in this paper are consistent with physical experiment in Fredd's paper. A numerical example is used to test the accuracy of the mathematical model in Panga's paper (Panga et al., 2010). The length and width of the model are 5 and 2 cm, separately. The initial porosity is 0.2, the fluctuation range is 0.05, V 2 10 5, N ac 0.1, h T 2 0.07. The comparison cures are shown in Figure 3. It shows that the two cures are not coincident. The reasons are the difference between accurate value and solving algorithm.

NUMERICAL RESULTS AND DISCUSSION
In this section, the numerical simulation of reactive flow in the fractured medium is implemented based on the above dimensionless mathematical model. The parameters and values in the numerical simulation are shown in Table 1. Unless otherwise specified, all parameters remain unchanged. The corresponding dissolution patterns in complex fractured media are studied, then, the effect of fracture aperture and the fracture orientation on dissolution structure and breakthrough volume of injected acid are studied in a single fracture system. The effect of heterogeneity of initial fracture aperture and dynamic variation of fracture aperture with dissolution on dissolution patterns, dissolution structure, and breakthrough volume of injected acid are studied. In this case, Da 250.
The definition equation of the breakthrough volume of injected acid is as follows: Where, V acid is the consumption volume of injected acid; V ip is the initial pore volume of the core; Q acid is the flow rate of injected acid, T bt is waste time that injected acid breaks through core when the pressure rate of inlet end to outlet end is 0.01 at each time point.

Dissolution Patterns at Different Damköhler Number
The results of the numerical simulation are shown in Figure 4. There are five different dissolution patterns with Damköhler number changing. They are face dissolution, conical dissolution, wormhole, branch dissolution, and uniform dissolution. Figure 4A is the concentration field, Figure 4B is the porosity field. Column (A) is the concentration field and porosity field in the initial phase. (B) is the concentration field and porosity field when injected acid breaks through the core. It can be seen that there is a big difference among dissolution patterns. The reason is that the dissolution patterns are determined by the combined effect of convection velocity and reaction velocity of the solute in acid. When the convection velocity is smaller than the diffusion velocity and reaction velocity, the reaction between acid and rock is mainly controlled by the diffusion and reaction velocity of solute in acid. Therefore, the hydrogen ions transferred from the center of the pore to the surface of the pore will be completely consumed by the chemical reaction, and the rate of injected acid cannot promote the injected acid to penetrate the core deeper. In this case, the dissolution of acid in the core is named surface dissolution. As the rate of injected acid gets higher, so does the convection velocity of solute in acid, and the convection velocity is dominant for diffusion velocity and reaction velocity.
Frontiers in Earth Science | www.frontiersin.org November 2021 | Volume 9 | Article 765139 5 The solute in the acid flows away before it reacts on the pore surface of carbonate rock. This results in the core pores being filled with acid and leads to a uniform variation in pore porosity. There is no dominant channel in the core like wormhole, and the dissolution pattern is called uniform dissolution.
As can be seen from Eq. 1.4, the reaction rate is determined by the properties of acid and rock, and have nothing to do with the existence of fractures in the matrix. In this case, we can draw an important conclusion: in the case of the same acid-rock reaction system, whether there are fractures in the matrix has no effect on the dissolution patterns of the fractured medium.
In addition, the comparison results of numerical simulation show that the fracture is equivalent to the main advantage channel in the matrix. An injection rate that is too high or too low will not substantially change the dissolution patterns. This is because the injected acid has already reacted with the rock completely when the acid is injected into the core at a low injection rate and there is no excess acid to flow into the fracture. In the case of high injection velocity, the dominant role of fracture in the flow is negligible and may even play a negative role. In the above situations, the existence of fractures does not affect the dissolution patterns, however, the dissolution structure will change when the injection rate of acid is between the two. This is because the flow of injected acid is controlled by fracture conductivity, and the ultimate dissolution pattern is closely related to the distribution of fracture in the fractured medium.

The Effect of Fracture Orientations
In this section, the effect of fracture orientations on the wormhole structure is studied when Damkohler number is 250. Firstly, a fracture is located in the matrix and then the dissolution structure is obtained by numerical simulation. The results of numerical simulation are compared with the results obtained in a matrix
(A) is concentration field and porosity field of numerical simulation in the initial phase and (B) is the concentration field and porosity field of numerical simulation when injected acid breaks through the core. In Figure 5, it can be obtained from column (A), the existence and orientations of fracture in the matrix do not affect the formation of the initial wormhole, but the direction of wormhole formation is consistent with the injected port of fractures. (B) shows that the fracture orientations influence the orientation of the wormhole and the fracture gradually becomes the part of the wormhole. The generation trend of the wormhole is random and branching in the matrix without fracture. When the orientation of fracture and injection direction of acid is at an angle of 90°, as shown in row (d) of Figure 5, the formation direction of the wormhole is not controlled by the orientation of the fracture. Instead, it is consistent with the direction of injected acid. Different from  the wormhole in the matrix without fractures, the wormhole radius has increased and the direction of the wormhole has also changed in the fractured medium. This shows that the influence mechanism of fracture orientation on the formation direction of the wormhole is like to a high permeability zone. The breakthrough volume of injected acid in a matrix with different fracture orientations is drawn in Figure 6. It shows that when the orientation of fracture is consistent with the injection direction of acid or the orientation of fracture is similar with horizontal direction, the breakthrough volume of injected acid is smaller than the breakthrough volume in the matrix without fracture, the branching randomness of the wormhole is also getting weaker. When the direction of fracture is bigger than a certain angle, the breakthrough volume of injected acid becomes larger, and so does the consumption of acid. At this time, the effect of fracture inhibits the flow of acid in the fractured medium. The waiting time in the fractured medium becomes larger and then consumes more acid. When the direction of fracture is between 67.5°and 112.5°, the breakthrough volume of injected acid increases first and then decreases. It is because the conductivity of fracture has a smaller impact on the development and formation of the wormhole. When the direction of fracture is bigger than 90°, the breakthrough volume of acid is similar to volume when the angle is less than 90°, the slight difference is due to the heterogeneity of the core.
In this example, the breakthrough volume of injected acid is largest in the core with a single fracture and the orientation of fracture is 67.5. The consumption of injected acid is even larger than the matrix without fracture.   In this section, the effect of fracture aperture on wormhole structure is studied when Damkohler number is 250. This model contains a fracture. The numerical simulation results are compared with the matrix without fracture, and the result of the numerical simulation is shown in Figure 7.
(A) is concentration field and porosity field of numerical simulation in the initial phase and (B) is the concentration field and porosity field of numerical simulation when injected acid breaks through the core in Figure 7. In column (A) of the concentration field and porosity field in Figure 7, whether the fracture exists or not and no matter how the fracture aperture changes, it does not affect the initial wormhole structure.
From column (B), it can be concluded that fracture aperture don't affect the formation trend of the wormhole and the fracture gradually becomes a part of the wormhole. The fracture aperture has no significant effect on the dissolution structure. The relationship between the breakthrough volume of injected acid and the fracture aperture is plotted as a columnar graph shown in Figure 8. Figure 8 shows that the breakthrough volume of injected acid is the smallest when the fracture aperture is largest. In the mathematical model of this paper, the calculation equation of permeability is obtained by the cubic law: the larger the fracture aperture and the higher the fracture permeability, the conductivity of fracture will become stronger. This means that the fracture has a stronger impact on the formation of the wormhole. As a result, the breakthrough volume of injected acid will become smaller and the consumption of acid will become less. In the same way, the smaller the fracture aperture, the lower the fracture permeability. The efforts of fracture on the formation of the wormhole are weaker and the breakthrough volume and consumption of injected acid are less than before. When the fracture aperture equals 5e-5m, the breakthrough volume of injected acid is the smallest in the study.
The Effect of Heterogeneity of Initial Fracture Aperture and Dynamic Variation of Fracture Aperture The numerical simulation results are shown in Figure 9. With the increase of injected acid and the Damköhler number, there are five different dissolution patterns in fractured medium: face dissolution, conical dissolution, wormhole, branch dissolution, uniform dissolution.
The breakthrough volume curve of injected acid is plotted in Figure 10. It depicts the effect of homogeneity of the initial fracture aperture, the heterogeneity of initial fracture aperture, and the dynamic change of fracture aperture to the breakthrough volume of injected acid.
As shown in Figure 10, The dissolution patterns are unchangeable regardless of whether the initial fracture aperture is homogeneous, heterogeneous or fracture aperture is dynamic change with the acid etching fracture wall, the optimal injection rate of acid solution will not change when Da 250. FIGURE 10 | Influence of homogeneity of initial fracture aperture, of the heterogeneity of initial fracture aperture, and the changing of fracture aperture on the breakthrough volume of injected acid.