Numerical Studies and Analyses on the Acidizing Process in Vug Carbonate Rocks

A vug porosity system, in addition to a matrix, is also the target rocks for the acidizing. In this work, the acidizing process in two typical core-scale separate-vug porosity systems is studied in detail. Numerous cases are conducted to discuss a parametric study on the acidizing process and hydraulic behavior. Results indicate that the presence of vug reduces the pore volume of acid solution consumed to achieve a breakthrough, which is consistent with experimental observations. Increasing the vug diameter and porosity decreases pore volume to breakthrough both for a vugular carbonate rock and isolated vug carbonate rock. In comparison, the acid mass does not change a lot. Typical dissolution patterns can also be observed in the acidizing process when a vug exists. Compared to matrix dissolution patterns, the presence of vug induces wormhole to pass through the vug region.


INTRODUCTION
Reactive flow in porous media plays an essential role in various physical, chemical, and biological processes in nature and has many applications (Fredd and Scott Fogler, 1998;Acharya et al., 2007;Li et al., 2008;Bijeljic et al., 2013;Nick et al., 2013;Menke et al., 2015;Li et al., 2018). This work is conducted to study the carbonate acidizing process in vug carbonate rocks.
Before, lots of studies have been conducted to study the carbonate matrix acidizing process. Both experimental (Bernadiner et al., 1992;Wang et al., 1993;Frick et al., 1994;Fredd and Fogler, 1999;Izgec et al., 2008;Furui et al., 2012a;Furui et al., 2012b;Dong, 2012;Aidagulov et al., 2018;Qiu et al., 2018) and numerical (Kanaka and Panga, 2003;Panga et al., 2005;Kalia and Balakotaiah, 2007;Kalia and Balakotaiah, 2009;Maheshwari et al., 2013;Maheshwari et al., 2016b;Liu et al., 2019) results show that an optimum acid injection rate widely exists, at which the minimum amount of acid solution is consumed to achieve a breakthrough (Jia et al., 2021). Meanwhile, different dissolution patterns (Hoefner and Fogler, 1989;Fredd and Fogler, 1999;Huang et al., 2000;Izgec et al., 2008;Aidagulov et al., 2018) can also be observed through visualization techniques at different acid injection velocities. Apart from the matrix, vug is also a widespread presence in carbonate rocks, which is a void space in the reservoir rock and provides the main storage space for reserves (Kossack and Gurpinar, 2001;Wang et al., 2016;Yue et al., 2018). According to the vug interconnection, a vuggy pore space can be divided into the separate-vug porosity system and toughing-vug porosity system (Nair et al., 2008). Because of complex porosity systems, a flow behavior is more complicated in vug carbonate rocks (Sadeghnejad and Gostick, 2020). In the toughing-vug porosity system, a vug is usually associated with fabric dissolution and karst process (Lucia, 2007), in which an interconnected pore system forms and has a much larger permeability (Zhang et al., 2005) of average 10 3 -10 6 times more than that of the disconnected vug system (Kossack and Gurpinar, 2001). For a separate-vug porosity system, the addition of vug only increases the total porosity but has no significant contribution for the core permeability (Zhang et al., 2005). Such rocks are also the target reservoirs for the acidizing stimulation. However, a few studies have been performed to study the acidizing process in the separate-vug porosity system. For an experimental study, Izgec et al. (2008) have shown that the presence of vug does affect the acidizing process for calcite rocks. The pore volumes to breakthrough decreases with the increase of the vuggy fraction. Furthermore, Izgec et al. (2009) also developed a two-scale continuum model to simulate the acidizing process in the separate-vug porosity system, in which vuggy carbonates are divided into two distinct regions of the matrix region and the vuggy region. The Darcy equation is used in the matrix region, and the Stokes equation is used to account for the fluid flow in the vug (free flow) region. However, it is worth noting that when most of the matrix is dissolved, compared to the damping effect of porous media, the fluid phase viscous dissipation becomes dominant and cannot be ignored anymore. In such cases, the Darcy equation is not applicable any longer. Also, Izgec's model neglects mass exchange between solid and liquid phases. Hence, a more accurate mathematical model is still needed to study the acidizing process in the separate-vug porosity system.
In this work, the Stokes-Darcy equation is used to describe the fluid flow instead of the Darcy equation, and a source term is added to the fluid mass conservation equation for considering the dissolution mass term from the solid phase to the fluid phase. A detailed explanation of the two-scale continuum model is given in the next section. In the third section, numerous cases are conducted to discuss a comprehensive parametric study on the acidizing process, including acid injection velocity, acid injection concentration, vug filling degrees, vug diameter, and matrix porosity heterogeneity.

Reaction and Kinetics
In this study, a solid mineral is supposed to only containing calcite (CaCO 3 ) (Alkhaldi et al., 2010;Maheshwari et al., 2016a;Mahrous et al., 2017). For simplicity, an overall chemical reaction equation is adapted to represent the acid consumption and calcite dissolution, which is given below: Calcite (CaCO 3 ) and hydrochloric acid (HCl) react to form calcium chloride (CaCl 2 ) and carbonic acid (H 2 CO 3 ). The surface reaction rate is dependent on the concentration by the power law expression: R HCl k s c n HCl,s , where R HCl is the acid reaction rate, k s is the acid reaction rate at the solid surface, n is the reaction order, and c HCl,s is the surface acid concentration. The reaction rate of other components can be further determined according to Eq. 1.

Governing Equations
Modeling and numerical simulation of vuggy carbonate reservoirs is still challenging because of vug's presence (Popov et al., 2009). The modeling of vuggy media is traditionally done based on the multiple-continuum concept model. In this model, the matrix and vug are assumed to have uniform and homogeneous properties, and the intercommunication flow between the matrix region and the vug region is at a pseudosteady state (Kang et al., 2006;Wu et al., 2007). The continuum approaches behave pretty well in vug well-developed reservoirs, but the challenging part of this mode is that it is hard to determine accurate rock properties such as porosity and permeability (Yao et al., 2010). Another popular approach is based on the coupled Stokes and Darcy equations, in which the Darcy equation is only used in the porous region, and the Stokes equation is used in the free-flow region. Moreover, at the interface between the porous and free-flow regions, the additional boundary condition is employed to guarantee continuity of mass and momentum across the interface (Izgec et al., 2009;Yao et al., 2010). However, it is still very complicated and less practical to apply the coupled Darcy-Stokes approach to vuggy reservoirs. Because, on the one hand, the interface condition also introduces additional parameters such as material property, α BJ , in Beavers-Joseph-Saffman (BJS) condition (Beavers and Joseph, 1967;Saffman, 1971), which needs to be further determined numerically or experimentally. On the other hand, the Stokes equation used in the free-flow region can only handle the vug region free of any obstacles (Popov et al., 2009). The limitation of the coupled Darcy and Stokes approach can be overcome by using the Stokes-Brinkman equation, which uses a single equation rather than coupled equations to describe the fluid flow in entire vuggy reservoirs and avoids explicit modeling of the interface (Popov et al., 2009;De Oliveira et al., 2012;Fadlelmula F et al., 2016). In this work, a new two-scale continuum model is inspired by the Stokes-Brinkman equation. The Stokes-Darcy equation is employed to predict the momentum transport behavior of the fluid phase. The equation is given below: where ϕ is the porosity, ρ fluid is the density of fluid phase, u u v T is the superficial velocity, k is the permeability, μ is the viscosity, and p is the fluid pressure. The equation can represent all regions during the acidizing process. For example, in an undissolved matrix region with a low permeability, the damping effect of the solid structure is relatively significant, which is represented by the last term in Eq. 3. Moreover, in wormhole and vug regions, the influence of solid skeleton is gradually diminishing with a large permeability, and the viscous dissipation of fluid behaves dominantly, which can be considered by the viscous term in Eq. 3. Then we continue to develop the mass conservation equations. According to Equation 1, the fluid phase consists of water (H 2 O), acid (HCl), calcium chloride (CaCl 2 ), and carbonic acid (H 2 CO 3 ). And the acid component (HCl) is chosen as the main component. The matrix only contains calcite (CaCO 3 ). Therefore, the governing equations for the fluid phase, acid solution, and calcite are independently developed by (Jia et al., 2021), which is as follows: where s stands for the specific surface area; M HCL , M CaCl2 , M H2CO3 , and M CaCO3 are the component for HCl, CaCl 2 , H 2 CO 3 , and CaCO 3 , respectively; c HCl is the acid concentration in the bulk phase; D HCL is the acid diffusion coefficient, determined as D HCL ϕD m I; I is the unit tensor; D m is the molecular diffusion rate; and ρ CaCO3 is the matrix density. The main difference between the conventional two-scale continuum model (Panga et al., 2005;Kalia, 2008) and the model used in this study is reflected in the last term in Eq. 4, which represents a source term to the fluid mass conservation equation for considering the dissolution mass term from the solid phase to the fluid phase. Almost all of the previous numerical studies neglect it. The acid reaction rate (R HCl ) is further determined as follows: R HCl k s c n HCl,s k c c HCl − c HCl,s , where k c is the mass transfer rate, and it is quantified by a model derived by Balakotaiah (Balakotaiah and West, 2002) Sh where Sh is the Sherwood number, Sh ∞ is the asymptotic Sherwood number, the value used in this work is 3.66 for a circular pore structure, Re p is defined as Re p 2|u|r/], ] is defined by ] μ/ρ liquid , and Sc is the Schmidt number defined by Sc ]/D m . The modified Carman-Kozeny equations are used to determine the changes of core parameters (Edery et al., 2011;Jia et al., 2021) which is given as follows: where c, α, and β are the parameters describing pore connectivity and broadening structure, and ϕ 0 , k 0 , r 0 , and s 0 are the initial porosity, permeability, pore diameter, and specific surface area, respectively.

Boundary and Initial Conditions
The calculation domain is selected as a 2D cross section of a cylindrical core with 4 cm diameter and 10 cm height. The boundary and initial conditions refer to the core flooding experiments (Izgec et al., 2008). Acid is injected at a constant velocity and concentration from one side of the core sample. The acid injection concentration is identical to the value, 15%, used in the experimental study. Constant pressure is held on the other core sample side. Because in the experiment process, the core sample is always surrounded by a high confining pressure to avoid the acid solution going through the gap between the core sample and the core holder plug. Initially, the core sample is saturated with water, and matrix porosity has a value of 0.15. Matrix permeability is assumed isotropic and has a value of 5 mD. Furthermore, a uniform distribution random function creates a heterogeneous distribution of porosity (Maheshwari et al., 2016a;Wei et al., 2017). Boundary and initial conditions are summarized below: n · v 0, n · (uc HCl + D HCl ∇c HCl ) 0 on transverse boundaries, (16) c HCl 0, ϕ ϕ 0 + Δϕ 0 at t 0, where n is the boundary unit direction vector, u 0 is the injection velocity, c 0 is the injection concentration, p e is the pressure at the core outside, and Δϕ 0 is the initial heterogeneity magnitude of porosity.

RESULTS AND DISCUSSIONS
The separate-vug porosity system covers a wide range of structures from millimeters to meters in the diameter (Arns et al., 2005;Sok et al., 2010;Mousavi et al., 2012;Sadeghnejad and Gostick, 2020). This study is concerned with two typical corescale separate-vug porosity systems: the vugular carbonate rock and the isolated vug carbonate rock, as shown in Figure 1. For numerical studies, vug in vugular carbonate rock is supposed uniformly embedded in the matrix region, and each vug initially has a diameter of 3 mm (Casar-Gonzalez and Suro-Perez, 2000; Izgec et al., 2010). Moreover, for the isolated vug carbonate rock, vug has a centimeter diameter, 1.5 cm (Zhang et al., 2004;Casar-Gonzalez and Suro-Perez, 2000). For a separate-vug porosity system, the addition of vug only contributes to the increasing storage space. The acid breakthrough condition is still consistent with matrix acidizing, determined by a 100 times increase of the core average permeability (De Oliveira et al., 2012). The pore volume at a breakthrough condition is counted and normalized with the sample initial pore volume as an evaluation parameter cited as pore volume to breakthrough (PVBT) (Jia et al., 2021). Besides, acid breakthrough injection mass (AMBT) is also collected in this work because PVBT is simultaneously a function of porosity.

Model Comparison and Validation
The two-scale continuum model is widely accepted and used to describe the core scale stimulation process (Kalia and Balakotaiah, 2007;Liu et al., 2017;Maheshwari et al., 2013). The Darcy equation is usually employed to determine the superficial velocity of the fluid phase. In this study, the Stokes-Darcy equation, as shown in Eq. 3, is employed to predict the fluid phase momentum transport behavior. When the acid solution acquires a very high injection rate, the core sample is uniformly dissolved by acid and still keeps the characteristic of porous media because of no existence of a wormhole. At this time, the Darcy equation can still be applicable. The acid breakthrough consumption volume calculated based on the Darcy equation has little difference from the value of PVBT calculated based on the Stokes-Darcy equation (Jia et al., 2021). Thus, the model can be validated by comparing PVBT values with the previous two-scale continuum model with high acid injection velocities. Table 1 shows the values of PVBT separately calculated by the current and previous two-scale continuum models with high acid injection velocities.
The results show that the difference of PVBT values between the current and previous two-scale continuum model is not noticeable when acid acquires a very high injection rate. Figure 2 further compares dissolution images of current and previous two-scale continuum models with very high acid injection velocities. The current two-scale continuum model is well validated.

Effect of Acid Injection Velocity
This section studies the effect of the acid injection rate on the acidizing process of the vugular carbonate rock and the isolated vug carbonate rock. The porosity of the vug region is initially set as a high value of 0.99. The whole sample porosity is calculated as 0.3875 for vugular carbonate rock and 0.1871 for isolated vug carbonate rock. The whole sample mean permeability of the vugular carbonate rock and the isolated vug carbonate rock is measured as 10 and 6.6 mD, respectively. In comparison, the permeability of the matrix zone is 5 mD. The results show that the existence of vug in the vugular carbonate rock and the isolated vug carbonate rock does provide more storage space but has a little contribution to hydraulic conductivity for the matrix core sample. Figure 3 shows the acidizing results of vugular carbonate rock. It can be observed from Figure 3A that the pore volume to breakthrough (PVBT) of vugular carbonate rock is less than the PVBT of matrix carbonate rock. The same change trend of the PVBT is also observed in the experimental study (Izgec et al., 2008). The acid mass to breakthrough (AMBT) change of the vugular carbonate rock and the matrix carbonate rock is not apparent compared to the change of the PVBT. The AMBT value   Figure 3B. Figure 4 compares the dissolution patterns of the vugular carbonate rock and the matrix carbonate rock with different acid injection velocities. Typical dissolution patterns, including uniform, ramified wormhole, wormhole, and conical wormhole dissolution patterns (Jia et al., 2021), are also observed in the acidizing process of vugular carbonate rock. Compared to dissolution patterns of matrix carbonate rock, vug presence induces wormhole to pass through the vug region. Besides, more wormholes can be stimulated due to vug existence with a high acid injection rate, as shown in Figure 4B. Also, because of the existence of vug, less acid is used to dissolve the solid matrix. The apparent difference in pore volume to breakthrough is mainly caused by the vug increasing initial pore volume. Furthermore, according to the definition of PVBT, a higher initial pore volume leads to lower PVBT. Figure 5 and Figure 6 compare acidizing curves and dissolution patterns of the isolated vug carbonate rock and matrix carbonate rock, respectively. Similar conclusions can also be observed with the vugular carbonate rock. Because the presence of isolated vug is increasing the initial pore volume, the PVBT of isolated vug carbonate rock is higher than the PVBT of matrix carbonate rock, which is also consistent with the experimental observation (Izgec et al., 2008). Four typical dissolution patterns can also be observed in the acidizing process of the isolated vug carbonate rock, as shown in Figure 6. The induction effect of isolated vug on wormhole growth is more intuitive for the isolated vug     Frontiers in Earth Science | www.frontiersin.org July 2021 | Volume 9 | Article 712566 6

Effect of Acid Injection Concentration
Acid injection concentration is also significant for the practical acidizing application. Initially, Fredd and Fogler (1999) adopted the concentration of 2% in core-flooding experiments. The concentration of 15% is widely used in the subsequent experimental studies (Furui et al., 2012a;Dong et al., 2014;Aidagulov et al., 2018;Kumar et al., 2020). Wang et al. (1993) ever used different concentrations such as 0.5, 3.4, and 15% to experimentally study the influence of acid injection concentration on the acidizing results. The results indicate that increasing acid injection concentration decreases the PVBT, while the optimum injection rate change is not apparent. Figure 7 numerically calculates the acidizing curves with acid injection concentration changing from 5 to 25%. It can be found that when the acid injection concentration increases from 5 to 25%, PVBT values shows a pronounced decrease both for the vugular carbonate rock and isolated vug carbonate rock, which is consistent with experimental studies conducted by Wang et al. (1993) and numerical studies of matrix acidizing presented by Jia et al. (2020a). Besides, acid mass to breakthrough (AMBT) is also compared with different acid injection concentrations. Because it is very important to note that when we are studying the influence of acid injection concentration on acidizing results, the amount of acid solute and acid solvent are both changing. The change of PVTB reflects more information on the amount of acid solute (Fredd and Fogler, 1999;Panga et al., 2005). Few studies before noticed and discussed the change of acid solvent mass. Figure 7 shows that compared to the apparent change of PVBT, the amount of acid mass required to breakthrough does not change apparently with increasing acid injection concentration. Hence, the amount of solute in the acid system has an apparent relationship with the acid injection concentration, but the solvent required to achieve a breakthrough has a less obvious influence.

Effect of Vug Porosity
The filling phenomenon is widespread in vug carbonate rocks. It is affected by multiple factors, such as mechanical deposition, collapsed deposition, chemical deposition, and so on, which eventually leads to the vug with different porosities (Li, 2017). This section studies the influence of vug porosity on the acidizing process of vugular carbonate rock and isolated vug carbonate rock. Figure 8 summarizes the change of core porosity and permeability, with the vug porosity varying from 15 to 99% considering different degrees of filling. The results show that both core porosity and permeability increase as the vug porosity increases. For the vugular carbonate rock, compared to the change of core porosity increasing from 15 to 40%, the change of core permeability is limited, increasing from 5 to 10 mD. The increase of core porosity and permeability in isolated vug carbonate rock is not apparent. Core porosity increases from 15 to 19%, and core permeability increases from 5 to 5.6 mD. Figure 9 further compares the acidizing results with different vug porosities from the low filling degree with vug porosity of 15% to the high filling degree with vug porosity of 99%. The results show that vug porosity mainly has a noticeable influence on the pore volume to breakthrough. Higher vug porosity results in a lower value of the PVBT. However, the difference of acid injection mass consumed to achieve breakthrough between different vug porosities is still not apparent.

Effect of Vug Diameter
This section studies the influence of vug diameter on the acidizing performance of vugular carbonate rock and isolated vug carbonate rock. Here, we suppose that for the vugular carbonate rock, the diameter of vug has the magnitude of a millimeter (Casar-Gonzalez and Suro-Perez, 2000; Ramakrishnan et al., 2001;Sok et al., 2010;Mousavi et al., 2012). Figure 10 shows core porosity and permeability with a vug diameter varying from 0.5 to 3.5 mm. It can be observed that when the vug diameter reaches 3.5 mm, core permeability is different from core permeability of other core samples with a vug diameter less than 3.5 mm. Besides, core porosity with a vug diameter of 3.5 mm reaches 47.33%. Hence, the final vug diameter below 3.5 mm is determined as the threshold diameter for the acidizing stimulation for the vugular carbonate rock. Then, we study the effect of isolated vug diameter on the acidizing process. The diameter of isolated vug on the core scale is supposed to have a centimeter magnitude (Casar-Gonzalez and Suro-Perez, 2000;Zhang et al., 2004;Arns et al., 2005). Unlike the vugular carbonate rock, the core porosity and permeability slowly increase with an increase of the isolated vug diameter. With the isolated vug diameter varying from 0.5 to 3 cm, the core porosity increases from 15.41 to 29.84% and core permeability increases from 5.1 to 6.9 mD. Finally, two varying vug diameters from 0.5 to 3 mm and 0.5 to 2.5 cm are separately chosen to study the influence of vug diameter on acidizing curves, as shown in Figure 11. The results indicate that the PVBT decreases as the vug diameter increases. The difference of the AMBT between different vug diameters is also not noticeable. For the isolated vug carbonate rock, the results indicate that the PVBT decreases as the isolated vug diameter increases. Acid mass does not change a lot with different isolated vug diameters. The noticeable difference of  Frontiers in Earth Science | www.frontiersin.org July 2021 | Volume 9 | Article 712566 10 the AMBT in different isolated vug diameters is only shown when acid solution is pumped at an intermediate injection velocity of 0.05 cm/s. Compared to a core sample with a big isolated vug diameter, less acid is needed to achieve the breakthrough with a low isolated vug diameter. In acid solution with intermediated rates, a wormhole-shaped dissolution image is formed, as shown in Figure 12. The isolated vug with a higher diameter induces a wormhole to pass through the isolated vug, which leads to more acid consumption. Lower isolated vug diameter has little effect on the development of wormhole.

Effect of Matrix Porosity
Unlike sandstones, the pore structure of carbonate rock is more complex due to complicated geology processes, which makes the carbonate rock with obvious porosity heterogeneity (Sok et al., 2010;Kim et al., 2011;Milad et al., 2020). This section studies the effect of matrix porosity on the acidizing process of vugular carbonate rock and isolated vug carbonate rock. In a recent work by Jia et al. (2020b), lots of rock sample porosity values measured by experiments are collected and compared. The results show that matrix acidizing target carbonate rock samples usually have rock porosity varying between 0.05 and 0.3. Figure 13 first studies the effect of matrix porosity on carbonate core porosity. It can be founded that the matrix porosity has a more evident effect on the vugular carbonate rock porosity. The core porosity of vugular carbonate rock increases from 0.31 to 0.46, with an increase of the matrix porosity from 0.05 to 0.25. While for isolated vug carbonate rock, the core porosity only increases from 0.09 to 0.28. Figure 14 further compares the effect of matrix porosity on acidizing results. It can be observed from Figure 14 that matrix porosity mainly has an apparent influence on the stimulation process when acid is injected into the vugular carbonate rock and isolated vug carbonate rock core sample with high rates. Acid consumption mass increases as matrix porosity increases. Compared to the high acid injection rates, the acidizing curve difference becomes less obvious. It is because when the acid injection rate is high, the acid transportation process is more governed by the convection process, and acid is preferred to being uniformly transported into the core sample. Core samples with higher matrix porosity can provide more storage space, which results in more acid being stored in the solid matrix. When acid injection velocity is low, the PVBT slightly decreases as the matrix porosity increases.  Frontiers in Earth Science | www.frontiersin.org July 2021 | Volume 9 | Article 712566 CONCLUSION Apart from the matrix region, the presence of vug is very common in the carbonate rock, which provides the main storage space for reserves. In this work, the acidizing process in two typical core-scale separate-vug porosity systems, the vugular carbonate rock and the isolated vug carbonate rock, is studied in detail. The Stokes-Darcy equation is used to describe the fluid flow to avoid using different equations in the free-flow region and porous media region. Parameter sensitivity analysis is preformed to ascertain the effect of acid injection velocity, vug filling degrees, and vug diameter on acidizing results. The main conclusions are presented below: 1) Acid injection velocity and concentration have a noticeable influence on pore volume to breakthrough (PVBT) both for the vugular carbonate rock and isolated vug carbonate rock. The existence of vug reduces the PVBT values of both vugular carbonate rock and isolated vug carbonate rock, which is consistent with the experimental observations. In contrast, acid injection velocity has little influence on acid mass to breakthrough (AMBT). 2) Typical dissolution patterns, including uniform, ramified wormhole, wormhole, and conical wormhole dissolutions, can also be observed in the acidizing process when vug exists. Compared to dissolution patterns of matrix carbonate rock, the wormhole is induced to pass through the vug region.
3) The presence of vug in a separate-vug porosity system contributes to more storage space but has little contribution to hydraulic conductivity. Core porosity and permeability both increase as filling degrees of vug and vug diameter increase. However, the increasing core permeability is limited. 4) Increasing the vug diameter and vug porosity decreases the pore volume to breakthrough both for the vugular carbonate rock and isolated vug carbonate rock.
To sum up, this study provides a powerful two-scale continuum model for the 2D vug acidizing simulation. The improved two-scale continuum model can avoid to use different governing equations in the vug region and the matrix region. The fluid phase continuity equation is also modified to consider the solid matrix dissolution mass term into the fluid phase. Further, the 3D vug simulation will be performed in our following studies.

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. Frontiers in Earth Science | www.frontiersin.org July 2021 | Volume 9 | Article 712566