Effect of Isolated Fracture on the Carbonate Acidizing Process

Carbonate reservoirs are one of the most important fossil fuel sources, and the acidizing stimulation is a practical technique for improving the recovery of carbonate reservoirs. In this study, the improved two-scale continuum model, including the representative elementary volume (REV) scale model and the upscaling model, is used to study the acidizing process with an isolated fracture. Based on this model, a comprehensive discussion is presented to study the effect of the physical parameters of the isolated fracture on the acidizing results and dissolution images, including the isolated fracture geometry, location, and morphology. Results show that the isolated fracture system is still the target system for the acidizing stimulation. The isolated fracture provides a limited contribution to the core porosity. The permeability of the core sample with fracture can be obviously increased only when the fracture penetrates through the whole sample. The existence of the isolated fracture reduces the consumption of acid solution to achieve a breakthrough. The acidizing curve is sensitive to the change of the length, aperture, and position of the isolated fracture. The acidizing curve difference corresponding to different rotation angles has not changed significantly for clockwise rotation and anticlockwise rotation groups.


INTRODUCTION
Carbonate reservoirs are one of the significant parts of fossil fuel sources (Zhao et al., 2018;Liu et al., 2020). Compared with sandstone sedimentary rocks, carbonate reservoirs have strong heterogeneity due to complex geological processes, which result in more formation energy consumption in some carbonate reservoirs with low permeability (Zhao et al., 2020;Liu et al., 2021). Carbonate acidizing is a widely used stimulation technology to effectively improve the oil recovery of carbonate reservoirs using the acid solution to dissolve the low permeability matrix (Wei et al., 2019).
Large numbers of articles have been published to study the carbonate matrix acidizing process. Experimental studies are mainly conducted through core flooding experiments, in which constant concentration acid solution is pumped at the inlet boundary of the core sample (Fredd and Fogler, 1999;Dong, 2012). During the experimental process, the pressure at the inlet boundary of the core sample is always recorded as an indirect indicator variable for monitoring the acidizing implementation. As the pore structure inside the core sample is constantly dissolved by the acid solution, the migration resistance of the acid fluid inside the core sample also constantly changes, causing the pressure at the core entrance to vary from time to time. Usually, the core displacement test is stopped when the inlet boundary pressure of the core sample is observed to decrease significantly (Zakaria and Nasr-El-Din, 2016). The obvious decrease of pressure means that the fluid migration resistance inside the core sample has dropped significantly, which is usually defined as the breakthrough moment during the acidizing process. The amount of acid consumed at the breakthrough time has been widely determined as a quantitative indicator to study the effectiveness of the acidizing process and evaluate the impact of different influencing factors on the acidizing process, including the acid solution types, concentration, injection rate, and core sample geometry (Dong, 2012;Furui et al., 2012;Wadekar and Pandya, 2014;Sarmah et al., 2020). One of the most important findings of these experimental studies is that there exists an optimum value for the acid breakthrough consumption curve, at which the injected acid solution consumes the minimum amount of the acid solution to achieve a breakthrough (Fredd and Fogler, 1999;Furui et al., 2012). Another important finding is that when acid solution is injected into the core sample at different rates, the final dissolution structure inside the core sample is also different when achieving the breakthrough (Sarmah et al., 2020). The high acid injection rate makes the acid solution to be transported to most of the core sample, causing the solid matrix to be uniformly dissolved. When the acid injection rate is reduced at the optimum value, the migration of the acid solution is controlled by the convection, diffusion, and core heterogeneity together. Finally, a major dissolution channel is formed inside the core sample named as the wormhole. When the acid solution injection rate further decreases, the width of the main wormhole is enlarged because the acid solution obtains sufficient time to dissolve the solid matrix.
Several mathematical models have also been developed to capture and reproduce the important results observed in the experimental studies, including the pore network model, empirical model, and continuum model (Maheshwari et al., 2013). The pore network model is the most realistic model, in which the initial calculation domain uses the pore network obtained through the true target stimulated core sample. The acidizing process is simulated inside each pore and throat of the real core sample (Tansey, 2014). However, huge calculation expense decides that the pore network model can only be conducted within the limited domain (Tansey and Balhoff, 2016). The empirical model has the smallest calculation expense, but the model accuracy depends on the assumptions for deriving the empirical model and the benchmark data obtained through the experimental studies (Maheshwari et al., 2016;Palharini Schwalbert et al., 2019). Finally, the continuum model has been extensively developed recently due to the high model accuracy compared to the empirical model and less computational consumption compared to the pore network model. The continuum model was first proposed by Panga et al. (2005) including the Darcy-scale model and pore-scale model. Hence, the continuum model is also referred to as the two-scale continuum model (Panga et al., 2005;Maheshwari et al., 2013;Liu et al., 2017). Recently, the twoscale continuum model was further improved by Jia et al. (2021a) to more accurately simulate the stimulation process, in which the additional mass term is added in the fluid phase continuity equation to consider mass exchange term between the solid matrix and acid solution. Also, the Stoke-Darcy equation is used instead of the Darcy equation to describe acid solution flow (Jia et al., 2021b).
Apart from the solid matrix, fractures are also very common in carbonate reservoirs. And, the carbonate rocks with isolated fractures are also the target reservoirs for the acidizing stimulation (Li and Voskov, 2021). However, few studies are conducted to study the acidizing process in carbonate rocks with isolated fractures. One of the reasons is that the efficient and accurate characterization of fractures is still a challenge in numerical simulation. Usually, the fracture models can be divided into the discrete fracture model and the duel-medium model according to the fractures' descriptions. The discrete fracture model explicitly distinguishes the location of each fracture, and the fluid flow in each single fracture can be simulated (Neuzil and Tracy, 1981;Tsang and Witherspoon, 1981). The simplest discrete fracture model uses two parallel, infinite planes to represent single fracture (Snow, 1969). The cubic law is used to determine the fracture permeability used in the Darcy equation (Muralidhar, 1990). In actual application, the fracture aperture is generated through a random function to consider the effect of the roughness of the fracture surface (Zimmerman et al., 1991;Zimmerman et al., 1992). Different from the discrete fracture model, the duel-medium model assumes that the simulated domain contains both the solid matrix and the fracture (Moreno et al., 1988;Tsang and Tsang, 1989). The governing equations are developed both for the solid matrix region and the fracture region (Meyer and Bazan, 2011). The flow exchange between the solid matrix region and the fracture region is a function of the fluid pressure difference in the solid matrix region and the fracture region (Wu, 2015). It is worth noting that the fluid flow in the above fracture models is all described by the Darcy equation. However, the matrix acidizing results have proven that the applicability of the Darcy equation is only validated when the acid injection rate is relatively high, and the solid matrix is uniformly dissolved. When wormholes appear, the Darcy equation is not able to accurately describe the fluid flow because it only considers the effect of the structure of porous media (Jia et al., 2021b). However, in the wormhole region, the complex porous structure is dissolved into free-flow region. And, for the fracture region, it can also be seen as a special wormhole region, in which the applicability of the Darcy equation is still challenged and needs to be further validated. Hence, a more general mathematical model is still needed to simulate the acidizing process with fractures.
In this study, the improved two-scale continuum model is used to study the acidizing process with an isolated fracture. The detailed model descriptions are given in the next section. In the discussions and results section, both the classical and the improved two-scale models are used to simulate the acidizing process with an isolated fracture. Several cases are designed to investigate the effect of the fracture physical parameters on the stimulation results and compare the differences between the classical and the improved two-scale models. The main conclusions are summarized in the last section.

Model Descriptions
The mathematical models used in the acidizing process simulation with an isolated fracture is referred as the improved two-scale continuum model (Panga et al., 2005), which contains the representative elementary volume scale model and the upscaling model.

Representative Elementary Volume Scale Model
The representative elementary volume scale model is developed based on the mass and momentum conservation principles. The chemical reaction between the solid matrix and the acid solution proceeds according to the following equation: (1) The specific reaction rate is expressed by the concentration of hydrochloric acid in the acid solution system: where k s is the acid surface reaction rate and c HCl,s is the acid surface concentration. Then, the mass conservation equations of the fluid phase, acid solution, and solid matrix are given in turn as follows: z ϕc HCl zt −∇ · (D HCl · ∇c HCl + uc HCl ) − R HCl s, where ϕ is the porosity, ρ f is the fluid phase density, u is the fluid phase velocity, M CaCO 3 is the molar mass of calcite, s is the specific surface area, c HCl is the acid bulk concentration, D HCl is the acid diffusion coefficient, and ρ s is the solid phase density. As for the determination of fluid phase velocity, the Darcy equation is most commonly used in both matrix and fracture regions (Wu, 2015): where μ is the fluid viscosity, k is the permeability, and p is the fluid phase pressure. Because it is usually accepted that fluid flow inside the porous media is mainly influenced by the porous media structure, however, our recent study results (Jia et al., 2021b) have proved that acid solution migration in porous media is affected by the porous structure for acidizing processes. The influence of viscous dissipation between acid solutions becomes more significant in the wormhole region and should also be considered. Instead, the Stokes-Darcy equation is adopted in the improved two-scale continuum model to determine acid solution convection velocity (Jia et al., 2021b): In this work, both the Darcy equation and the Stokes-Darcy equation are used to simulate the acidizing process with an isolated fracture to further compare the effect of differences between the classical and the improved two-scale models on the acidizing results.

Upscaling Model
The upscaling model is developed to link the actual acidizing process with the REV scale model. In the actual stimulation process, the chemical reaction between the acid solution and the solid matrix occurs on the pore surface inside the solid matrix. However, the reaction terms involved in the model established on the REV scale are based on the control volume. The connection between the actual chemical reaction rate at the pore scale and the equivalent, chemical reaction term on the REV scale is achieved through the mass transfer equation: where k c is the acid mass transfer rate, which is determined by an analytical expression (Balakotaiah and West, 2002): where Sh, Sh ∞ , Re p , and Sc are the dimensionless numbers, the asymptotic Sherwood number, Sh ∞ , is usually valued as 3.66 for circular pore structure, pore Reynold number, Re p , and Schmidt number, and Sc is defined as Re p 2r|u|/] and Sc ]/D m , ] is the dynamic viscosity. During the acidizing process, the continuous dissolution of the solid matrix by the acid solution causes the porous structure to change from time to time. The modified Carman-Kozeny equations update the change of the physical parameters of the solid matrix (Panga et al., 2005;Maheshwari et al., 2013): where k 0 , r 0 , and s 0 are the initial rock parameters of permeability, pore radius, and specific surface area, respectively, and γ and β are the parameters related to pore connection and cementation. At the same time, the porous structure of the solid matrix also affects the diffusion process of the acid solution: where D m is the molecular diffusion coefficient, I is the unit vector, r is the pore radius, |u| is the fluid velocity norm, and u and v are the components of fluid phase velocity.

Initial and Boundary Conditions
The application of the initial and boundary conditions is consistent with the core displacement experiment (Fredd and Fogler, 1999;Dong, 2012). The core sample is saturated with water at the beginning. When the acidizing stimulation begins, the acid solution is continuously injected from the core inlet boundary at a constant injection rate and concentration, and the right boundary of the core sample maintains a constant pressure. The surrounding of the core sample is kept closed (Dong et al., 2014).

RESULTS AND DISCUSSION
The study domain is selected as a common core sample size adopted in the experimental study with a dimension of 4 cm width and 10 cm length. The initial porosity and permeability of the core sample are measured as 0.15 and 5 m D (Dong et al., 2014). The isolated fracture is located in the center of the core sample with a length of 5 cm and an aperture of 0.15 mm. The acid system is selected as a hydrochloric acid solution with an injection concentration of 15% (4.42 mol/L). The molecular diffusion coefficient of the acid solution is determined as 3.5 × 10 −5 cm 2 /s and surface reaction rate is determined as 1.1 × 10 −3 cm/s (Panga et al., 2005;Jia et al., 2021b). The values of other physical parameters used in numerical models are summarized in Table 1. The governing equations are discrete and solved based on the finite element method. During the calculation process, the core inlet boundary pressure is detected to monitor the acidizing process. The calculation stops when the core apparent permeability increases 100 times (Dong et al., 2014). The apparent permeability is determined based on the Darcy equation using the pressure difference between the inlet and the outlet boundaries. The volume of acid consumed at the calculation stops is normalized by the core initial pore volume to obtain the evaluation variable for the acidizing result, pore volume to breakthrough (PVBT) (Panga et al., 2005).

Model Comparison with the Existence of an Isolated Fracture
This section mainly investigates the influence of the isolated fracture on the acidizing process. We first study the influence of the single fracture on the core porosity and permeability. Figure 1 shows two typical single fracture patterns inside the core sample, separately named as the connected fracture and the isolated fracture. The porosity values are determined as 0.1815 and 0.1657, respectively, for the connected fracture and isolated fracture core samples. The permeability of the isolated fracture core sample is numerically measured at 7.94 mD, while the  Frontiers in Earth Science | www.frontiersin.org July 2021 | Volume 9 | Article 698086 permeability of the connected fracture core sample is at 50 mD. The matrix porosity and permeability is 0.15 and 5 mD. It can be observed that both the isolated fracture and the connected fracture have a limited influence on porosity. The isolated fracture also has limited contribution for the core permeability. Only the connected fracture obviously increases the core permeability. Hence, the core sample with an isolated fracture is still the target core sample for the acidizing stimulation.
In the classical two-scale continuum model (Panga et al., 2005), the Darcy equation is used to describe the fluid flow during the acidizing process. In our recent work (Jia et al., 2021b), the Stokes-Darcy equation is suggested to determine the fluid phase velocity instead of the Darcy equation. Because the applicability of the Darcy equation is questioned within the wormhole regions, in the wormhole regions, the porosity value is very high, even up to 1. The effect of the porous structure on the fluid flow becomes less and less obvious. The viscous dispersion between the fluid phase becomes more important and should be considered. Figure 2 shows the carbonate matrix acidizing curves separately determined by the Darcy equation and the Stokes-Darcy equation. The results support the application of the Stokes-Darcy equation. When the acid solution is injected at high rates, the difference between the two acidizing curves is not very obvious. Because most of the solid matrix is not totally dissolved into wormholes with high acid injection rates. The effect of porous structure on the fluid flow is still dominated in this case compared to fluid phase viscous dispersion. The acidizing curve difference becomes apparent with low acid injection rates. When the acid injection rate is very low, the acid solution has sufficient time to react with the solid matrix, causing most of the solid matrix to completely dissolve into wormholes. Because the Darcy equation cannot consider the effect of fluid viscosity dissipation, the acid breakthrough volume calculated by the Darcy equation is obviously less than that determined by the Stokes-Darcy equation. The acidizing efficiency of the acid solution is artificially enhanced by using the Darcy equation. Figure 3 further compares the acidizing curves separately determined by the Darcy equation and the Stokes-Darcy equation with an isolated fracture. With contrast to Figure 2, the difference between the acidizing curves with an isolated fracture becomes obvious within the whole injection rate interval, even though the acid injection rate is at high values. The results show that the existence of an isolated fracture also challenges the applicability of the Darcy equation at high injection rates. Figure 4 continues to compare the acidizing curves with an isolated fracture and without an isolated fracture. The results show that the acid breakthrough volume with an isolated fracture is smaller than that of the matrix acidizing result within in both acidizing curves determined by the Darcy equation and the Stokes-Darcy equation. The acid breakthrough volume difference between the acidizing curves with and without an isolated fracture is quite obvious when the acid injection rate is at a high value for both the Darcy equation and the Stokes equation. The acid breakthrough volume difference between the acidizing curves with and without an isolated fracture starts to differ at lower injection rates. The acid breakthrough volume difference determined by the Stokes-Darcy equation gradually becomes less apparent at low injection rates. The acid breakthrough volume difference determined by the Darcy equation is still significant when the injection rate is low. Figure 5 compares the dissolution images of the core sample with and without an isolated fracture. The isolated fracture greatly influences the final dissolution patterns especially when the acid solution is injected into the core sample with high and intermediate rates. The isolated fracture provides an advantageous passage for acid migration in porous media, in which isolated fracture prevents more acid from being ineffectively consumed to dissolve the solid matrix. Table 2 further studies the effect of mass exchange term between solid matrix and acid fluid solution on the stimulation results with an isolated fractured, which is neglected by most previous numerical studies (Panga et al., 2005). The results show that ignoring the mass exchange term still leads to an overestimated acid breakthrough volume, which is consistent with the conclusions obtained from the matrix acidification process (Jia et al., 2021a).

Effect of Isolated Fracture Geometry
This section mainly studies the influence of fracture geometry on the stimulation process, including fracture length and fracture aperture. The fracture length is selected to vary from 3 to 8 cm, in which the fracture remains in the middle of the core sample with an aperture of 1.5 mm. The results show that the porosity of the core increases from 0.16 for a 3 cm fracture length to 0.18 for an 8 cm fracture length. The core permeability increases from 6.1 mD for a 3 cm fracture length to 15.1 mD for an 8 cm fracture length. The permeability and porosity of the matrix core sample correspond to 0.15 and 5 mD, separately. As for the influence of the aperture, the porosity increases from the matrix porosity of 0.15-0.18, corresponding to the fracture aperture of 3 cm. The permeability increases from the matrix  Frontiers in Earth Science | www.frontiersin.org July 2021 | Volume 9 | Article 698086 6 permeability of 5 to 8.2 mD. Hence, different fracture lengths and apertures have a very limited effect on the growth of core porosity and permeability. Figure 6 shows the acidizing curves corresponding to different fracture lengths. The results show that the acid breakthrough volume increases with the increase of the isolated fracture length. And, the existence of isolated fracture can obviously affect the distribution of pressure and velocity inside the core sample, as shown in Figure 7. Inside the matrix core sample, the acid fluid injection rate is parallel to the length direction of the core sample, and the pressure gradient also decreases uniformly along the core length. When an isolated fracture exists, the uniformly changing pressure gradient inside the matrix core sample is intuitively broken. The pressure within an isolated fracture is generally lower than the pressure of the matrix core sample. A pressure convergence zone and a pressure divergence zone separately forms at the inlet and the outlet boundaries of the isolated fracture. Hence, the fluid tends to flow into the convergence zone and the divergence zone and also leads to the formation of wormhole branches at the outlet boundary of the isolated fracture. Figure 8 continues to study the influence of fracture aperture on the acidizing results. The results show that the increase of the fracture opening can reduce the amount of acid consumed at breakthrough. High fracture apertures increase the affected domain generated by the isolated fracture and causes the acid solution migration within the core sample affected more obviously. More acid solution is induced by the high fracture aperture, resulting in less acid solution consumed.

Effect of Isolated Fracture Location
In this section, several cases are designed to study the isolated fracture locations on the acidizing results. The isolated fracture developed in the center of the core sample is selected as the reference comparison case. Here, nine different isolated    center-up, center-center, and center-down), and right group (right-up, right-center, and right-down). These nine types of core samples at different isolated fracture locations have the same porosity of 0.17. Figure 9 compares the acidizing curves corresponding to different groups of fracture positions. It can be found that there is no obvious difference between the acidizing curves with the same group of isolated fractures in the up, center, and down positions. The difference in the acidizing curves becomes obvious between the different groups of left, center, and right positions. The acid consumption volume with the fracture position at the left boundary of the core sample is larger than the breakthrough volume of the acid corresponding to the fracture position at the right boundary of the core sample. When the fracture is located at the entrance boundary, the existence of fracture provides an advantageous channel for the migration of acid solution. During the acidizing process, the acid solution tends to migrate along the main direction of the fracture. Until the acid moves to the right end of the fracture, the acid begins to come into contact with the core sample with low permeability, which causes the acid solution not to be transported as efficiently as it moves in the fracture. The acid solution is forced to accumulate in the fracture and ineffectively consumed. More acid solution is used to increase the width of the fracture, as shown in Figure 10. When the fracture is at the outlet boundary, the acid tends to contact the fracture. Moreover, once the wormhole formed inside the matrix begins to contact the fracture, the core breakthrough can be easily realized. Because the fracture directly communicates with the right boundary, the negative acid accumulation no longer exists. Less acid is needed to break through the core.

Effect of Isolated Fracture Morphology
This section studies the influence of fracture morphology on acidizing process by rotating the isolated fracture at a certain angle. Two rotation directions, clockwise and anticlockwise rotation, are adopted. Each type is rotated by 10, 20, and 30°, respectively. Finally, six different fracture morphologies are obtained. The isolated fracture results with 0°rotation is selected for reference comparison. Although the isolated fracture has different rotation angles inside the core sample, the change of permeability is not obvious. The permeability values are estimated around 7 mD, and the permeability of the matrix core corresponds to 5 mD. Figure 11 studies the acidizing curves corresponding to different rotation angles. The acidizing curves corresponding to different rotation angles are divided into the clockwise rotation and anticlockwise rotation groups. The results show that the change of the acidizing curve corresponding to different rotation angles of the same rotation sequence is very similar. The difference between them has not changed significantly both for clockwise rotation and anticlockwise rotation groups.
Then, we continue to consider an extreme case, in which the rotation angle of the isolated fracture reaches 90°, and the fracture is vertically distributed inside the core sample. Both the vertically distributed and the horizontally distributed isolated fractured cores have a porosity of 0.16. The permeability of vertical and horizontal isolated fracture cores is measured as 5.2 mD and 6.1 mD, respectively. The increase in porosity and permeability is still not obvious. However, the difference between the acidizing curves corresponding to the horizontally and the vertically distributed fracture begins to become obvious, as shown in Figure 12. Compared with the horizontal distribution of fracture, the vertical distribution of fracture simultaneously reduces the breakthrough volume of acid solution and the optimal injection rate. When the fracture is distributed horizontally, the horizontal tips of the fracture form a convergence zone and a divergence zone to attract the acid solution. The acid solution tends to penetrate the fracture through the inlet boundary and form a branch structure at the outlet boundary of the fracture, as shown in Figure 13. When the fracture becomes vertically distributed inside the core sample, the convergence zone and divergence zone are also formed around the fracture. Nevertheless, unlike the divergence zone formed by horizontal fracture, the divergence direction of the velocity corresponding to the vertical fracture divergence zone is consistent with the direction of core breakthrough flow. In contrast, the divergence zone formed by horizontal fracture maintains a certain angle with the direction of core breakthrough. Hence, the formation of the horizontal fracture divergence zone shows an impediment to the breakthrough of acid solution in contrast to vertically developed fracture.

CONCLUSION
Apart from the low-permeability matrix, carbonate rocks with isolated fractures are also the target reservoirs for the acidizing stimulation. In this study, the improved two-scale continuum model is used to study the acidizing process with an isolated Frontiers in Earth Science | www.frontiersin.org July 2021 | Volume 9 | Article 698086 9 fracture. Several cases are designed to present a comprehensive discussion on the isolated fracture parameters, including the isolated fracture geometry, location, and morphology. The main conclusions are summarized below: 1) Both the isolated fracture and the connected fracture have a limited influence on the core sample porosity. Only the connected fracture obviously increases the core permeability. The core sample with an isolated fracture is still the target core sample for the acidizing stimulation. 2) The isolated fracture provides an advantageous passage for the acid migration, in which the isolated fracture prevents more acid from being ineffectively consumed to dissolve the solid matrix reducing the amount of acid breakthrough. 3) Ignoring the mass exchange term still leads to an overestimated acid breakthrough volume, which is consistent with the conclusions obtained from the matrix acidizing process. 4) The existence of an isolated fracture also challenges the applicability of the Darcy equation at high injection rates.  Frontiers in Earth Science | www.frontiersin.org July 2021 | Volume 9 | Article 698086 5) The isolated fracture length, aperture, and position have an obvious influence on the acidizing process. The difference of the acidizing curve corresponding to different rotation angles has not changed significantly both for clockwise rotation and anticlockwise rotation groups.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material; further inquiries can be directed to the corresponding author.

AUTHOR CONTRIBUTIONS
CJ conceptualized the study, developed the methodology, and wrote the original draft; JY supervised and reviewed and edited the study; HX and HZ helped with software support; TH: supervised the study and carried out funding acquisition.