Effect of produced carbon dioxide on multiphase fluid flow modeling of carbonate acidizing

A two-scale continuum (TSC) numerical model is used at pore and Darcy scales to model and optimize the matrix acidizing dissolution process. During matrix acidizing process in carbonate reservoir, the solubility of carbon dioxide, as one of the hydrochloric acid/calcite reaction products, is limited, which leads to the formation of a separate gas phase in the reservoir based on pressure and temperature. The presence of this free gas has nonlinear effect on the fluid flow and phase distribution in porous media due to the alteration in the relative permeability of the injected/spent acid. Density and viscosity of the spent acid will also be changed. Ignoring these effects on TSC model can reduce the accuracy of the final results compared to similar laboratory studies. The significance of our study is to examine and apply the nonlinear effects of produced carbon dioxide on wormhole propagation, dissolution patterns, and breakthrough curves at different injection rates of the acid. The innovation presented in this study is the evaluation of the effect of the presence of the gaseous carbon dioxide produced during the hydrochloric acid and calcite reaction on the relative permeability of the injected acid. This effect has been numerically investigated and solved along with the governing equations of the matrix acidizing at different time steps. In our work, The TSC model is coupled with a nonlinear relative permeability model to consider the effect of free carbon dioxide. Density, viscosity, and solubility correlations are also coupled with TSC model to update fluid properties during acid propagation and acid/rock interactions. The final results show that the neglecting of free gas affects the accuracy of wormhole propagation estimation. The model shows better agreement with the experimental data to predict the rock dissolution patterns and wormhole breakthrough. By considering free gas, the model can track thicker wormholes due to the free gas blockage and the slower acid/rock reaction rate. Our approach showed that isothermal models are not accurate enough to predict the change in carbon dioxide solubility at higher temperatures. Reduction in solubility affects the acid's relative permeability and increases the volume of acid required for the breakthrough. Hence, the developed model improves our knowledge of wormhole propagation during carbonate acidizing.


Abbreviations η (Etta)
Diameter to length ratio of the core (dimensionless) ε (Porosity) Ratio of pore volume of the rock to bulk volume (dimensionless) ε 0 (Initial porosity) -(Dimensionless) α 0 (Aspect ratio) Height to length ratio of the core (dimensionless) ρ s (Rock density) -(Dimensionless) μ (Viscosity) -(Dimensionless) α (Acid dissolving power) Ratio of dissolved rock mass to gram of consumed acid (dimensionless) ϕ 2 (Thiele modulus) The ratio of reaction rate to diffusion rate (dimensionless) λ X , λ T (constants) Numerical constants related to dispersion coefficient (dimensionless) α os (constants) Numerical constant depends on the respective carbonate rock structure related to tortuosity (dimensionless) v (Kinematic viscosity) The ratio of fluid dynamic viscosity to fluid density (cm 2 /sec) a 0 (Initial specific surface area) The amount of rock surface available to the acid for reaction (cm 2 / cm 3 ) a v (Specific surface area) The amount of rock surface available to the acid for reaction (cm 2 /cm 3 ) A v (Dimensionless specific surface area) -(Dimensionless) Ratio of reaction rate to diffusion rate (dimensionless) r′ (Dimensionless pore mean diameter) -(Dimensionless) Re p (Pore-scale Reynold number) Ratio of inertial force to viscous force(dimensionless) r p (Pore mean diameter) -(Cm) r p0 (Initial pore mean diameter) -(Cm) Sc (Schmidt number) The ratio of kinematic viscosity to molecular diffusivity (dimensionless) Sh (Sherwood number) Ratio of mass transfer by convection to mass transfer by diffusion (dimensionless) Sh ∞ (Asymptotic Sherwood number) -(Dimensionless) t (Dimensionless time) -(Dimensionless) t' (Time) -(Sec) T (Temperature) -(°F) U (Darcy velocity in longitudinal direction) -(Cm/sec) u (Dimensionless Darcy velocity in longitudinal direction) -(Dimensionless) u 0 (Acid injection rate)

Introduction
From the initiation of drilling a well, various operations may lead to formation damage in the reservoir and production system (Civan 2015). Increased pressure drop and reduced formation permeability are the most common aspects of formation damage phenomenon (Radwan et al. 2019). The source of formation damage could be mud and cement invasion during drilling and cementing, deposition of asphaltene and wax at different temperatures and pressures, migration of fine particles, clay swelling, or incompatibility of fluids, etc. . A well-known operation to remedy formation damage and increase productivity is matrix acidizing (Radwan et al. 2019). Dissolution in pores during carbonate acidizing occurs in three steps. Initially, there is a mass transfer of acid, as the reactive component takes place from the bulk fluid to the rock surface. At the next step, the reaction between the acid and the rock occurs, and finally, the products of the acid-rock reaction are transferred to the bulk fluid, which is called the mass transfer step of the reaction components (Panga et al. 2005;Kalia and Balakotaiah 2007). The wormhole propagation strongly depends on the acid flow rate in the formation (Golfier et al. 2002). Different dissolution patterns may form, which affects the depth of acid invasion and bypasses the damaged zone. Figure 1 shows a typical acid response curve with different acid flow patterns and the volume required for acid to breakthrough during an acid flooding test .
Numerous experiments have been conducted, and Various models have been developed to predict the formation and propagation of wormholes at different acidizing patterns, which will be discussed briefly in the following.
One of the first experimental studies to investigate the wormhole propagation in carbonate rocks was conducted by Williams et al. (1979), who suggested that acid should be injected at the highest possible rate to prevent the formation of face dissolution patterns. Liu et al. (1997) modeled fluid flow equations in porous media and considered the rock dissolution effect, porosity/permeability changes, mass transfer, and acid/rock chemical reactions. In this model, the overall dissolution rate was controlled by the reaction. The model was validated by experimental results. In an empirical study, Fred and Fogler (1998) explored the dissolution process in limestone, using acetic acid, hydrochloric acid, and retardedgelled acids, and concluded that, for different types of acid/ rock systems, the volume of acid required for breakthrough strongly depends on the dimensionless Damkohler number. Panga et al. (2004) proposed a two-scale continuum model to describe the wormhole formation in carbonate rocks by considering the reaction and mass transfer mechanisms in linear flow. In 2007, Kalia and Balakotaiah developed the model proposed by Panga by applying it to radial flow and studied the effect of convection and molecular diffusion parameters on dissolution structures (Kalia and Balakotaiah 2007). In 2009, the same researchers investigated the effect of radial flow heterogeneity on dissolution patterns using the continuous two-scale model (Kalia et al. 2009).
By conducting experiments at different temperatures, Shukla et al. (2006) studied the effect of the presence of immiscible gas and oil phases on the development of wormholes during acid injection. They found that the reduction of the acid's relative permeability (due to the presence of an external phase) decreases the fluid loss and forms longer wormholes. The effect of the presence of gas phase before acidizing was confirmed as an effective factor in reducing the required volume of acid for breakthrough of the rock. Bastami et al. (2014) investigated the thermodynamic behavior of the H 2 O-CaCl 2 -CO 2 system as products of acid-rock reaction, to estimate the solubility of carbon dioxide and free gas saturation. Bastami and Pourafshary (2016) developed a continuum model considering the acid consumption and the effect of reaction products, such as free gas, on the acid flow and the wormhole propagation. Bastami et al. (2018) extended this work by considering a thermodynamic model to calculate the density of the H 2 O-CaCl 2 -CO 2 system. Cheng et al. (2017) performed core flooding experiments and CT scan observations to investigate the effect of carbon dioxide at different temperatures and pressures on wormhole behavior. They concluded that at low acid injection rates, the produced gas phase slows down the development of the wormhole and increases its diameter. At injection rates higher than the optimum rate, as well as at very high temperatures, the effect of carbon dioxide is usually negligible. Mahmoodi et al. (2018) developed a new two-phase twoscale continuum model to simulate the dissolution process during the carbonate acidizing by considering the secondary phase of oil. The effect of flow constraint parameters, due to relative permeability changes and viscous fingering on the volume of acid to breakthrough the rock and dissolution patterns in the synthetic model, were investigated. These researchers concluded that the presence of two-phase flow and the phenomenon of fingering lead to the development of narrower and longer wormholes. Muhemmed et al. (2021) presented a two-scale, three-phase model, considering the effect of soluble CO2. By modeling CO2 as a separate phase, their model correlated better with experiments. Sajjaat Muhemmed et al. (2021) study investigated the effects of evolved CO2 caused by carbonate-mineral dissolution, and its ensuing activity during the preflush stages in matrix acidizing of sandstone reservoirs The potential of evolved CO2, a byproduct of the sandstone-acidizing preflush stage, toward its contribution in swelling the surrounding oil, lowering its viscosity, and thus mobilizing the trapped oil has been depicted in this study.
In this study, we considered the effect of CO 2 and CaCl 2 as reaction products of acid/rock on the rock/fluid properties such as relative permeability, density, and viscosity during the acidizing, to enhance the accuracy of acid flooding models. The assumption used in this study is being first order of the reaction between acid and rock. CO 2 forms a separate gas phase at volumes higher than its solubility in reservoir fluid, which affects the flowing fluid permeability and dissolution patterns as well as the acidizing efficiency. In this study, the main focus is on CO 2 and its different behavior on the matrix acidizing efficiency. In a numerical study, Bastami and Pourafshary (2016) presented a model, which the effect of CO 2 (in the form of free gas) was considered to be a linear function of the relative permeability of the injected   fluid saturation, which is not an accurate assumption. This study develops a new nonlinear relative permeability model to consider the presence of acid/rock reaction products, their effects on the matrix acidizing optimal conditions and changes in the relative permeability during the injection/ flowing of acid. We believed that existence of CO 2 in the form of free gas prevents acid from entering some pores by affecting relative permeability of the injected acid, which increases the pore volume to breakthrough (PVTB). Also, the effects of temperature changes on acidizing efficiency curve in both quantitative and qualitative are presented in this work. According to experimental results, existence of carbon dioxide at high temperatures increases the volume of required injected acid for achieving breakthrough.

Modeling methodology
Equation 1 shows the reaction of hydrochloric acid (HCl) with limestone (CaCO 3 ).
It should be noted that in Eq. 1, the HCl liquid phase is reacted with solid rock of CaCO 3 , and solid salt (CaCl 2 or calcium chloride) and liquid phases of H 2 O and CO 2 are produced. Based on Eq. 1, during the carbonate acidizing with HCL, a ternary system of calcium carbonate, carbon dioxide, and water is formed as products. Depending on the reservoir temperature and pressure and the soluble salts composition of the reservoir fluid, carbon dioxide can remain in solution or form a separate gas phase (Duan and Sun 2003). The presence of free gas changes the viscosity and density, which in this research is taken into account to update the fluid properties at each time step. Also, a novel correlation between acid/gas relative permeability is presented. By using this correlation, the free gas saturation at each time step updates the relative permeability of the injected acid in the porous media. The complete process for solving the developed model equations is shown in Fig. 2. The details and sequences are shown in Fig. 2. The pressure, concentration, porosity, permeability, and velocity profiles should be updated by considering the effect of spent acid at every time step, to achieve breakthrough in the rock. The criterion for achieving breakthrough in this model is that the porosity of the last grid block reaches twice the average porosity. The acid volume required for the breakthrough and the dissolution pattern of the porous medium are the main results of the developed model. Different models and correlations are required to estimate the parameters shown in gray in Fig. 2, which will be presented in the appendix sections.
During numerical simulation, first, the magnitudes of pressure, velocities vector in different directions, acid concentration, and porosity distribution are calculated. Then, Fig. 2 Flowchart of solving the development continuum model section using product concentration equation, the molality of products is calculated. The maximum solubility of carbon dioxide is calculated by Duan and Sun (2003). Finally, the effect of free gas phase on the relative permeability of the flowing acid is estimated by a novel relative permeability-volume model.

Relative permeability model
By using the reaction product concentration equations and the equation of state (Duan and Sun 2003), the product's molar volume is calculated and, subsequently, the volume of carbon dioxide V CO 2 in the free gas phase is obtained. By using Eq. 2, the free gas CO 2 volume is converted to S CO 2 , and the saturation of gaseous CO 2 in each block is: where dx, dy, and dz represent the length elements and represent the porosity. Bastami and Pourafshary (2016) developed a linear model to consider the free gas relative permeability on injection acid, as shown in Eq. (3).
In this study, by considering matching coefficients (ω and β) in Eq. 3 a nonlinear model is developed, as shown in Eq. (4).
where K rel is the relative permeability of flowing acid, which plays an important role in the permeability distribution of the whole system of the two-scale continuum model. S CO 2 is the saturation of gaseous CO 2 , which is calculated by Eq. 2, S ir is the saturation of the liquid phase (injection fluid) at each time step, and ω and β are experimental matching coefficients and are calculated based on experimental results. in this study, the magnitudes of 1.3 and 1.8 have been selected for omega and beta, respectively. S ir refers to the injected/ spent acid saturation, which decreases at every time step, due to the increase in gaseous carbon dioxide in the environment and the effect of CO 2 on the relative permeability of the flowing acid. By increasing the saturation of insoluble carbon dioxide, the relative permeability of the injected acid is reduced and the acid is prevented from entering some of the pores, which can deepen the penetration of structural patterns and slow down the reaction. On the other hand, by decreasing S ir , K rel decreases and the permeability distribution in the basic equations is also decreased. It should be noted that the relative permeability model is a function (4) k rel = S ir − S CO 2 S ir of gas phase saturation (SCO2) and liquid phase saturation (S ir ).

Model validation
In the first step, the base model is constructed through input data used by the numerical study of Kalia et al. (2009), and subsequently validated with the result of that study, as shown in Fig. 3. The base model shows an acceptable match with the numerical result presented by Kalia.
At low injection rates, as the acid cannot invade into the formation, the rock surface is dissolved and a face dissolution pattern is developed by consuming a high volume of acid. If the injection rate increases, the acid penetrates the formation and develops conical channels. By further increasing the injection rate, the optimal dominant wormhole is developed. In this condition, the volume of acid required for breakthrough is the minimum, while the acidizing process is at the optimum condition. At higher injection rates, a branching wormhole pattern is formed. Finally, at very high injection rates, the acid penetrates all pores and dissolves uniformly (Panga et al. 2005). In  addition to the qualitative comparison in Fig. 3, quantitative results show the good match between this study and Kalia et al. (2009) as shown in Table 1: The dissolution patterns predicted by the model are compared with the CT scan observations conducted by Golfier et al. (2002), as shown in Fig. 4. The current model is capable to recognize a dominant wormhole at a rate of 50 cc/ hr as Golfier et al. captured the similar single wormhole at this rate. Figure 4 shows that the current model is capable of predicting the other dissolution patterns at similar rates of Golfier et al. study.
To validate the developed model predictions, the acid flooding experiment conducted by  was modeled. The model was run to estimate the injection volume required to breakthrough for different acid injection rates. The parameters applied in the model, based on the experiments, are shown in Table 1.
The model output is shown in Fig. 5. The developed model shows an acceptable match with the experimental results . It should be noted that detailed information of the rock properties, especially porosity distribution, used in the  is not available. Therefore, an exact match is impossible (Table 2).
To validate the model in the case of the presence of free gas during acid flooding, the results of the developed model were compared with the observations of a two-phase experimental study by Cheng et al. (2017). Figure 6 shows the acid response curves measured in the laboratory and those predicted by the present model for different temperatures. Experimental studies show that the required acid volume to breakthrough increases as the temperature increases, which is similar to what was estimated by the developed model. At higher temperatures, the solubility of carbon dioxide in the liquid phase decreases and a higher volume of gaseous CO 2 is developed, which affects the relative permeability of the injected acid and, subsequently, increases the pore volume to     Table 3. The structural patterns of two-phase output at the two desired temperatures are presented in Fig. 7, which confirms the last point. As shown in Fig. 7, at higher temperatures, due to the presence of more CO 2 in the gas phase and slower dissolution, the wormhole pattern becomes thicker. Hence, the validated model can be used to analyze the effect of governing parameters, such as acid concentration, heterogeneity, and temperature, on the propagation of wormholes in the presence of spent acid and acid/rock products. Also, the pressure parameter shows its effect in solubility, density, and viscosity relations, which discussed before, It should be noted that the effect of pressure and temperature are not considered in most of the previous acidizing models. However, in the current model, the effect of pressure and temperature on parameters such as solubility, viscosity, and density are considered, which alters the optimum injection rate. The good agreement between this study and Cheng et al. (2017) is shown in Table 4.

Influence of acid concentration and heterogeneity
In order to investigate the effect of acid concentration, pore volume to breakthrough at different injection rate values of three different concentrations was simulated. As shown in Fig. 8, by increasing the acid concentration, more solid is dissolved at the same time, which decreases the optimum required acid volume. It means that at high acid concentration cases, lower injection rate should apply to reach the PVTB. However, the optimum acid concentration should be calculated based on economic criteria.

Effect of free CO 2 on wormhole propagation pattern
The effect of CO 2 in the form of free gas on fluid flow parameters, such as relative permeability, affects the wormhole shape and size during acid flooding. To show this effect, in one case, the effect of free gas was neglected and, in the other case, it was considered. The difference in the results at different injection rates is shown in Fig. 9. At all injection rates, the estimated pore volume to breakthrough was higher    Fig. 8 Investigation of the effect of acid concentration on two-phase model when the effect of free gas was active in the model. The presence of gas decreases the movement of the injected fluid and slows down the acid flow in pores, which enhances the reaction mechanism. Hence, the effect of spent acid should be considered in models to improve the accuracy of the wormhole propagation prediction. The positive effect of the presence of free CO 2 on dissolution patterns is shown in Fig. 10. In Fig. 10, case (a) shows the estimated dissolution pattern of the single-phase model in which the carbon dioxide is completely soluble and the free gas is neglected. At the same injection rate, by considering the presence of free gas, the wormhole propagation pattern is altered, which can be observed in case (b). When the free CO 2 exists in the model, the dissolution patterns become thinner and less branched. The reason for this phenomenon at one certain flow rate, can be attributed to the change in the chemistry of the acid/rock reaction by insoluble carbon dioxide, which slows down the reaction and prevents acid from entering some pores and makes the wormhole narrower. It should be noted that, in this case, the volume required for breakthrough for single-phase and two-phase cases was 3.1 and 2.8, respectively, which shows a noticeable effect of the presence of the secondary phase, which cannot be neglected in the modeling.
In addition, it should be mentioned that the acceptable pressure range in this study is between 500 and 3000 psia and acceptable temperature range is between 70 and 150 degrees Fahrenheit. In the following, the plots presented in Figs. 4, 7, and 10 are presented in Fig. 11.

Conclusions
In this study, the effect of considering CO 2 was investigated as a separate gas phase in spent acid on the matrix acidizing modeling of carbonate formations. Evaluation of the effect of the presence of the gaseous carbon dioxide produced during the hydrochloric acid and calcite reaction on the relative permeability of the injected acid is the innovation of this study. This effect had been numerically investigated and along with the governing equations of the matrix acidizing (TSC models) in different time steps solved. Acidizing efficiency curve for both single-phase and two-phase mode was simulated and compared with the acid flooding experiments in validation sections. The model shows acceptable agreement with the experiments result in both single-and two-phase modes. The following points are obtained from our study: • Applying the effect of free gas to model enables it to capture mechanisms such as preventing the entry of injection fluid into some pores and branching. • Considering the presence of free carbon dioxide in the spent acid in the modeling affects the fluid flow parameters such as viscosity and density of the fluid flow. • By using the developed model, the effect of temperature on wormhole pattern and the optimum acid injection rate was analyzed. At high temperatures, the solubility of the gas in the flowing fluid decreases. Therefore, the CO 2 saturation in the gas phase is increased, which has a negative effect on the relative permeability of the injected acid and leads to a higher acid volume for breakthrough. • By evaluating the structural patterns at different temperatures, it was found that the presence of carbon dioxide thickens the flow paths by slowing down the movement of the acid and delaying the reaction. On the other hand, at the optimum condition, presence of CO 2 makes wormholes thinner and longer and decreases PVBT, in contrast to the single-phase flow. • The effect of CO2 on dissolution pattern in both quantitative and qualitative sights are investigated. The volume required for breakthrough for single-phase and twophase cases was 3.1 and 2.8, respectively, which shows a noticeable effect of the presence of the secondary phase, which cannot be neglected during the matrix acidizing modeling.

Appendix Appendix A
The following equations are the governing relations of TSC model: The following dimensionless numbers were used to develop dimensionless governing equations.
Hence, the dimensionless TSC model is: To express the pore structure parameters (permeability, pore radius, and specific surface area distribution) in acidizing modeling equations, semi-empirical property-structure models by Kalia and Balakotaiah (2007) are used, as shown in Eqs. (12)(13)(14).
The transfer of acid species from the reservoir of bulk fluid to the pore-fluid interface is quantified by the Sherwood number. The local pore structure, the reaction rate, and the local velocity of the fluid affect this number. In this study, the model suggested by Kalia and Balakotaiah (2007) was used to estimate the Sherwood number as:

Carbon Dioxide Solubility Model
The thermodynamic model proposed by Duan and Sun (2003) is used to estimate the solubility of carbon dioxide in the electrolytic solution of pure water-calcium chloride, as The method is based on the balance between the chemical potential of carbon dioxide in the liquid and gas phases. m CO 2 , y CO 2 and CO 2 express the soluble mole fraction, mole fraction in vapor phase, and fugacity coefficient related to carbon dioxide, respectively. P is the total pressure, R is the universal gas constant, T is the absolute temperature, and m shows the molality of different salts in solution. Three interaction parameters ( CO 2 , CO 2 −Na and CO 2 −Na−Cl ) are also used in the model.

Density and Viscosity Models
Adding CO 2 to the electrolyte solution changes the density and viscosity of the solution. To calculate the density, one must estimate the density of pure water, water containing calcium chloride, and the density of the ternary solution. Batzle and Wang (1992) presented a simple model to calculate the density of pure water w as a function of temperature and pressure, as shown in Eq. 24. Kemp et al. (1989) proposed the model presented in Eq. (4) to estimate the density of water containing a combination of several salts b .
A relation proposed by Bastami and Pourafshary (2018) was used to estimate the density of the aqueous solution of calcium chloride containing CO 2 aq , as shown in Eq. 26.
To predict the viscosity of such solutions, first the viscosity of pure water w is estimated by the Islam and Carlson (2012), see Eq. (27).
Then, according to the model presented by Mao and Duan (2009), the viscosity of water containing salt sol is estimated by: where m and r represent the salts molality and the ratio of the brine viscosity to the pure water viscosity, respectively. The A, B, and C constants are temperature-dependent parameters. Finally, Eq. 30 is used to calculate the viscosity of a ternary solution containing CO 2 (Islam and Carlson 2012).
x CO 2 and H 2 O−NaCl represent the molar fraction of carbon dioxide in the aqueous phase and viscosity of pure water along with the electrolytic components, respectively.

Products concentration equation
In Eq. 31, C p is the concentration of the product and ν is the stoichiometric ratio of the product. This equation is solved for each one of the products separately after calculating the concentration and porosity of media in every time step.
Funding There is no funding for this manuscript.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes (26) aq = b + 2.271 * 10 −1 x CO 2 + 1.6129 * 10 2 x 2 CO 2 . (31) were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.