Numerical investigation of fluid phase momentum transfer in carbonate acidizing

Cun-Qi Ji ，Kmy Sepehrnoori ，Hi-Yng Zhng ，Yong-Fei Yng ，Jun Yo ，*

a School of Petroleum Engineering，China University of Petroleum (East China)，Qingdao，266580，Shandong，China

b Hildebrand Department of Petroleum and Geosystems Engineering，The University of Texas at Austin，Austin，TX，78705，USA

c Petroleum Department，Khalifa University，127788，Abu Dhabi，United Arab Emirates

Keywords:Fluid phase momentum transfer Carbonate acidizing Continuum model Reactive flow simulation

ABSTRACT This work mainly studies the effect of fluid phase momentum transfer mechanisms on the acidizing results，including the retardation effect of the porous structure and the interaction between the fluid phase，such as viscous dissipation and inertial effect.The results show that the acid fluid momentum transfer is influenced by the complex porous structure and fluid viscous dissipation.Eventually，the Stokes-Darcy equation is recommended to be adopted to describe the fluid phase momentum transfer in the following numerical simulation studies of the carbonate acidizing process.Based on this model，a parametric research is carried out to investigate the impact of acid on rock physical characteristics in the stimulation process.Increasing the acid concentration appears to minimize the quantity of acid consumed for the breakthrough.The acid surface reaction rate has a considerable impact on the pore volume to breakthrough and the optimum acid injection rate.The influence of permeability on the acidizing results basically shows a negative correlation with the injection rate.The difference between the acidizing curves of different permeability gradually becomes insignificant with the decrease in injection rate.The existence of isolated fracture and vug significantly reduces acid consumption for the breakthrough.

1.Introduction

Reactive flow has many practical applications in porous media，including transport and reaction within catalysts in the chemical industry，storage and sequestration of carbon dioxide in the underground environment，and acidizing stimulation in carbonate reservoirs.This work is mainly inspired by carbonate acidizing.Carbonate reserves play a significant role in the world's fossil energy supply，up to 60% (Xu et al.，2020).Nevertheless，carbonate rocks are more complicated identified with sandstones because of complicated geological functions，which result in some carbonate reservoirs less efficiently developed due to low permeability matrix(Ghommem et al.，2020;Yuan et al.，2019).Acidizing is one effective stimulation technique for these carbonate reservoirs，in which acid solution is used to dissolve the solid matrix for reducing the fluid flow resistance and ultimately improving the recovery (Wei et al.，2019;Zhang et al.2019，2020).

Several studies have been carried out in the past to study the acid stimulation process.The experimental research is mostly performed through core scale displacement experiments (Dong et al.，2014).During the experiments，the acid fluid is injected into the target experimental core samples under certain conditions，and the pressure change between the target stimulated core samples is observed(Furui et al.，2012).When the inlet pressure of the core sample drops significantly，it is said that the carbonate core sample has achieved a breakthrough.The amount of acid consumption at this time is referred to as the breakthrough volume，which is further used to quantitatively compare the influence of different parameters on the acidizing results，including acid type，rock type，acid injection velocity，acid concentration，viscosity，temperature，and core size (Fattah and Nasr-El-Din，2010;Fredd and Fogler，1999;Kumar et al.，2020;Sarmah et al.，2020).Through comparing the breakthrough volume of the acid under different experimental conditions，it is widely discovered that an optimum acid injection rate exists，which corresponds to the least quantity of acid consumption for breakthrough (Dong et al.，2016;Qiu et al.，2018;Wang et al.，1993).When the injection rate deviates from the optimum injection value，the breakthrough volume of the acid solution changes and increases accordingly.Besides，the core breakthrough's dissolution image is also different from the injection rate change (Aidagulov et al.，2018).Usually，the dissolution image formed at the optimum acid injection rate is vividly called a wormhole，because the shape is very similar to the movement trajectory of an earthworm under the ground.When the injection rate of the acid fluid is less than the optimum threshold value，more fresh acid solution is expended to extend the width of the wormhole，and a conical dissolution eventually forms-ramified dissolution image forms when the corrosive acid injection rate is higher than the optimum threshold value since more wormholes are activated at this time (Fredd and Fogler，1999;Maheshwari et al.，2013;Zhang et al.，2020).

Apart from experimental studies，few numerical models have also been proposed to reproduce key conclusions observed in experimental research and guide the practical application，including the pore network model，empirical model，and continuous model (Jia et al.2021a，2021b;Li et al.，2008;Maheshwari et al.，2013;Palharini Schwalbert et al.，2020;Tansey，2014;Zhang et al.，2021).The pore network model is simulated at the porescale level，and its accuracy is higher than the empirical and continuous models (Tansey，2014).The pore network model accurately distinguishes the solid skeleton and pores to simulate the entire process of solid skeleton dissolution during the acidizing process.However，accurate characterization of the pore structure also brings a huge amount of calculation expense，which greatly limits applying the high-precision pore network model (Liu and Mostaghimi，2017a，b;Tansey and Balhoff，2016).The empirical model has been widely used due to its simplicity，which is mainly based on the principle of mass conservation.The empirical model correlates the main factors affecting the growth of wormholes by summarizing and analyzing experimental and numerical simulation findings，and finally obtains the empirical relationship for predicting wormhole characteristics (Akanni and Nasr-El-Din，2015;Li et al.，2008).Therefore，the accuracy of the empirical model is largely determined by the number of experimental studies and the accuracy of numerical simulation studies (Li et al.，2015;Palharini Schwalbert et al.，2019).The continuum model is developed as a compromise between the pore network model and the empirical model.The continuum model adopts the continuum hypothesis to focus on the mass，momentum，and energy transfer of the acid solution occurring inside the representative elementary volume (Jia et al.，2021a).Hence，the continuum model can avoid the huge calculation expense caused by distinguishing the pore structure.Meanwhile，to improve the simulation accuracy，several pore-scale models are supplemented as constitutive equations to describe the acid transfer and reaction process occurring at the pore-scale and update the constantly changing core physical parameters during the acidizing process (Maheshwari et al.，2013;Panga et al.，2005).Therefore，the continuum model is also called the two-scale continuum model and has been extensively utilized to investigate the acidifying process of different applications (2-dimensional and 3-dimensional) and explore the influence of different parameters including rock heterogeneity，oil saturation，and fracture existence (Liu et al.，2017c;Liu et al.，2021;Maheshwari et al.，2013;Mou et al.，2019;Qi et al.，2018).More detailed mathematical model descriptions can be found elsewhere(Jia et al.，2021a).

The two-scale continuum model has been broadly expanded due to the effective reduction of calculation expense and ensuring high simulation accuracy.The researchers also find that the continuum model still has some hidden shortcomings that need further improvement.For example，the fluid continuity equation of acid solution is modified by Jia et al.(2021a) to add a mass source term to characterize the mass change caused by the dissolution of the solid matrix into the fluid phase during the acidizing process.Meanwhile，Jia et al.(2021a) also point out that the wormhole region also influences the determination of the velocity of acid fluid migration.The Darcy equation was widely used in previous numerical studies to determine the apparent velocity of the reactive fluid system in porous media，such as the reactive flow involved in the catalyst and CO2geological capture process.Unlike the above reactive processes，the chemical reaction involved in the acidizing process is more intense (Fredd and Fogler，1999).The solid porous matrix can be completely dissolved and degenerated into an almost free flow domain within the wormhole region is an implicit challenge to the applicability of the Darcy equation.The Darcy equation is mainly used to describe the retardation effect of the permeable structure for fluid flow(Panga et al.，2005;Qi et al.，2018).The fluid movement within the wormholes is influenced by the permeable structure and the interaction between fluids，such as viscous dissipation and inertial effect，gradually increasing as the porous structure degrades.Recently，a few works have begun to pay attention to the influence of fluid viscous dissipation and inertial effect on the acidizing process，but there is still no definitive conclusion.Zhang et al.(2014)and Mou et al.(2019)think that the effect of viscous dissipation is dominant compared to inertial terms.Akanni et al.(2017) hold the opposite view.Moreover，both terms are important in Ali and Nasr-El-Din models (Ali and Nasr-El-Din，2019).Some studies still think that the transport of acid solution is still dominated by the retardation effect of the porous structure，and the interaction between fluids can be ignored(Liu et al.，2021;Wei et al.，2019).This work is mainly to conduct a comprehensive study on fluid phase momentum transfer in carbonate acidizing.

In this work，the different mechanisms affecting the momentum transfer of the acid system，such as the retardation effect of the porous structure，fluid viscous dissipation，and inertial effect，are studied and compared separately to find out the main control mechanism and obtain a general fluid phase momentum transfer equation.Based on this model，a parametric study is conducted to compare the influence of acid and rock physical parameters on the acidizing results.A detailed description of the model is given in the next section.Finally，the conclusions are displayed within the final segment.

2.Model descriptions

The governing equations are mainly developed based on the continuum hypothesis to describe the changes of mass and momentum in the acidizing system.In this study，the acid solution is assumed to be a hydrochloric acid (HCl) system，and the mineral composition of the solid matrix is composed of calcite(CaCO3)(Liu and Mostaghimi，2017b;Wei et al.，2017).The reaction relationship between the solid component calcite and the hydrochloric acid(HCl) solution is given by

During the stimulation procedure，the solid calcite (CaCO3) is continuously dissolved by the reactive acid solution (HCl) and converted into calcium chloride(CaCl2)and carbonic acid(H2CO3).The carbon dioxide always keeps soluble in the produced products under high pressure of 1000 psi and above(Dong，2012;Dong et al.，2014).The reaction rate is determined by

where ksis the rate of acid surface reaction，cHCl，sis the concentration of reactive acid on the solid surface，n is the reaction order，usually valued as 1 for the hydrochloric acid and calcite system.Here，hydrochloric acid is selected as the main component，and Eq.(1) can be used to calculate the chemical reaction rate for other components.Then，the governing equations are developed to investigate the change of mass and momentum of the acid fluid phase，main acid component，and the solid matrix during the carbonate acidizing process.The continuity equation gives the mass conservation of the acid fluid phase:

where ρfis the fluid phase density，φ is the porosity，u is the velocity of fluid phase，s is the specific surface area，MCaCO3is the calcite molar mass.Eq.(3) indicates that the fluid phase mass change in the controlling volume is not only caused by the inflow and outflow but also the mass supplement caused by the continuous dissolving of a solid matrix into a fluid phase.Most previous numerical simulations almost ignore the mass exchange term caused by the solid matrix dissolution (Maheshwari et al.，2013;Panga et al.，2005).It has been demonstrated incorrect to ignore the mass exchange between the solid matrix and the fluid phase，which eventually leads to a larger breakthrough amount of acid consumption (Jia et al.，2021a).The mass conservation equation of the main component in the acid system is subsequently given by

where DHClis the acid diffusion rate，cHClis the acid bulk concentration.Eq.(4)shows that the mass change of the main component in the acid system is mainly caused by the acid transport process and reaction process.The transport of acid solution mainly includes convection and diffusion processes.Unlike the diffusion process occurring in the solution system，the complex pore structure of porous media also has an impact on the diffusion process of the acid solution，and its influence decreases with the simplification of the pore structure，which is given by(Panga et al.，2005)

where r is the pore radius，｜u｜is the norm of the fluid phase velocity，u and v are the fluid phase velocity components，Dmis the molecular diffusion rate，I is the unit vector.Besides，it is important to note that the acid concentration related to the chemical reaction rate is the concentration on the solid surface.The concentration solved by Eq.(4)is the concentration in the bulk phase.These two acid concentrations are linked with the acid mass transfer formula:

where kcis the rate of acid mass transfer，determined by(Balakotaiah and West，2002)

where Sh，Sh∞，Rep，and Sc are the dimensionless numbers，the asymptotic Sherwood number，Sh∞，commonly regarded as 3.66 for circular pore structure，pore Reynold number，Rep，and Schmidt number，Sc，is defined by Rep=2r｜u｜/ν and Sc=ν/Dm，respectively，ν is the dynamic viscosity.The fluid phase velocity is determined by the momentum conservation equation of the acid solution.However，the fluid phase momentum transfer equation has no unique expression in numerical studies on carbonate acidizing studies，as shown in Table 1.Hence，it is necessary to conduct a comprehensive study on fluid phase momentum transfer to obtain a general fluid phase momentum transfer equation for modeling the carbonate acidizing process.

Table 1Fluid phase momentum transfer equations used in the continuum model.

The momentum transfer of the fluid phase in porous media is theoretically determined by the retardation effect of the porous structure，fluid viscous dissipation，and inertial effect.However，the momentum transfer equation can have different expressions based on the above three different mechanisms (Le Bars and Worster，2006;Nield and Bejan 2006;Whitaker，1999).The most extensive expression is derived from the up-scaling method and referred to as the Naiver-Stokes-Darcy equation，which can simultaneously describe the above three mechanisms within one expression:

where k is he rock permeability，μ is the viscosity of the fluid phase，p is the fluid phase pressure.Furthermore，the most widely used expression is the Darcy equation in the actual application process:

The Darcy equation is also the most used in acidizing numerical studies(Panga et al.，2005;Liu et al.，2017c).It can be regarded as a special form of Eq.(8)，which only considers the retardation impact of the porous structure for fluid flowing and ignores the viscous dissipation and inertial effect of the fluid phase.When considering the fluid phase viscous dissipation and the retardation effect of porous media，Eq.(8) can be degenerated into

Eq.(10)is referred to as the Stokes-Darcy equation.In this work，Eqs.(8)-(10) are separately used as the fluid phase momentum transfer equation to calculate the fluid phase velocity and study the effect of the different momentum transfer mechanisms on the acidizing results.

In summary，the changes of mass and momentum in the acid system are governed by Eq.(1)to Eq.(10).The boundary conditions and initial conditions are determined based on the experimental applications (Dong et al.2014，2016;Fredd and Fogler，1999).During the core acidizing flooding experiment，the carbonate core sample is initially saturated with water.The acid is then pumped into the core sample at a constant rate and concentration，and constant pressure is also applied to the outlet boundary of the core sample to eliminate the influence of gaseous carbon dioxide.The transverse boundaries of the core sample always remain closed.The equations are given by

At the inlet of the core sample，

At the outlet of the core sample，

As for the governing equation for the solid matrix，it is only controlled by the mass conservation equation since the velocity of the solid matrix is usually zero，which is

where ρsis the solid phase density.Eq.(15) shows that the acid dissolution causes the mass change of the solid phase.Because the solid phase density usually remains constant，Eq.(15) is mainly used to update porosity in the simulation of the acidizing process.The initial value of the core porosity is obtained from the experimental observations.The rock heterogeneity is then represented through a random function，where the initial porosity distribution is produced by a random function fulfilling a uniform distribution(Wei et al.，2017).The mean value of the uniformly distributed random function is set as the initial core porosity.As the acid solution continuously dissolves the porous medium，the change of the other rock physical parameters，including permeability，pore radius，and specific surface area，is described by the modified Carman-Kozeny equations (Edery et al.，2011;Maheshwari et al.，2013;Panga et al.，2005):

where γ and β are the pore connectivity and cementation characteristics，k0，r0and s0are the initial permeability，pore radius，and specific surface area of the rock，respectively.

3.Results and discussion

In this section，several numerical cases are performed to study the effect of different momentum transfer mechanisms on the acidizing results and eventually obtain a general fluid phase momentum transfer equation.In addition，a parametric study is conducted to discuss the influence of acid and rock physical parameters on the stimulation process.The study domain is performed in a 2-dimensional rectangular area with a length of 10 cm and a width of 4 cm.The fluid system is assumed to contain the hydrochloric acid solution，and the solid matrix consists of calcite.The precise acid and rock physical parameter values utilized in this work are consistent with the values used already conducted numerical studies (Jia et al.，2021a;Panga et al.，2005) and experimental measurements (Dong et al.，2014)，as summarized in Table 2.

In the specific calculation process，the governing equations are discretely solved by the finite element method.The segregated method decouples the unknown variables instead of a fully coupled method to improve the numerical model's stability and reduce the calculation expense.The fluid velocity is initially calculated by combining the fluid phase's mass and momentum conservation equations.The acid transport and reaction equation is then solved in order to get the acid concentration distribution.Finally，the solid phase continuity equation is solved for updating the solid matrix porosity.The updated porosity is substituted into the constitutive equations to prepare the other physical parameters for the next step.During the core displacement test，the experiment is stopped when the core inlet boundary pressure is observed as an obvious drop.Correspondingly，the core inlet pressure is constantly detected in the calculation process.The numerical simulation also stops when the core inlet pressure decreases significantly.The moment is generally determined as the pressure difference between core samples reduces to one percent of the initial pressure difference.The consumption amount of the acid is usually recorded as the pore volume to breakthrough (PVBT).

3.1.Model validation

This section is mainly to verify the numerical simulation method used in this study.The previous classical numerical simulation study conducted by Panga et al.(2005)is selected as the referencesolution，in which the calculation model is performed on a rectangular porous media of dimensions 2 cm×5 cm.The acid solution is injected at a constant rate at the core sample inlet boundary with a concentration of 5 mol/L.Acid surface reaction rate is valued as 0.2 cm/s.Moreover，the fluid continuity equation's mass exchange term temporarily is not considered to keep consistent with the reference model.Fig.1 first compares the influence of different grid sizes on the calculation results，in which degrees of freedom change from 6023 to 73493.The degrees of freedom are set at 68171 in the following calculation to eliminate the influence of the grid size.Table 3 summarizes the pore volume to breakthrough determined by the current work and the reference model.The results show that the amount of acid consumption at breakthrough calculated by current work fits well with the values in reference numerical simulations.

Table 2The acid and rock physical parameter values utilized in the current work.

Table 3Comparison of the acidizing results between this work and reference solution.

3.2.Comparison of fluid phase momentum transfer equations

The momentum transfer process of the fluid phase has not been adequately discussed in previous studies.In theory，fluid flow in porous media is primarily influenced by the retardation impact of the porous structure and fluid-fluid interaction，which includes viscous dissipation and inertial effect (Le Bars and Worster，2006;Nield and Bejan 2006;Whitaker，1999).However，in the actual applications，the interaction between fluids is often ignored due to the limited pore space，and only the influence of porous structure on fluid migration is considered.However，for the acidizing process，the original complex porous structure is continuously weakened as the acid solution continues to dissolve the solid matrix，which correspondingly aggravates the interaction between the fluids.Therefore，whether the viscous dissipation and inertial effects between fluids should be reconsidered requires further research.

Fig.1.Effect of the mesh size on the calculation results.

Fig.2.Effect of different fluid momentum transfer mechanisms on the acidizing results.

In this section，Eqs.(8)-(10) are separately used as the fluid phase momentum transfer equation to discuss the influence of the fluid inertia term，porous structure，and fluid viscous dissipation on the acidizing results.Fig.2 shows the pore volume values to breakthrough separately calculated using the fluid phase momentum transfer equations.The results indicate that the difference between the values of PVBT calculated by the Stokes-Darcy equation and the Naiver-Stokes-Darcy equation is insignificant with the change of injection rate.The main difference in the amount of acid consumption for core sample breakthrough is reflected in the calculation results of the Darcy equation and the Stokes-Darcy equation.Hence，the acid solution's momentum transfer is mainly affected by the complex permeable structure and fluid viscous dissipation in the acidizing process.The inertial effect of the fluid phase has little contribution to the fluid momentum transportation process.Specifically，F(xiàn)ig.2 shows that the difference between the PVBT calculated by the Darcy equation and the Stokes-Darcy equation is not apparent when the acid solution is injected at a rapid rate.The difference between the acid breakthrough amount becomes obvious when an acid solution is injected at intermediate and low injection rates.The value of PVBT calculated by the Darcy equation is higher than that determined by the Stokes-Darcy equation with intermediate velocity.When the acid injection rate is lowered further，the revised changing trend of the PVBT is observed.The quantity of acid used determined by the Darcy equation is smaller than the amount anticipated by the Stokes-Darcy equation.

Fig.3.The acidizing dissolution patterns at breakthrough with different injection rates:(a) 1 cm/s;(b) 0.1 cm/s;(c) 0.01 cm/s;(d) 0.001 cm/s.

When the rate of acid injection is relatively high，the acid solution's transport process is dominant compared to the reaction process.The acid solution prefers to be transported throughout the entire core sample，which eventually causes the solid matrix to be uniformly eroded，as shown in Fig.3a.The core sample still maintains a porous structure，and the viscous dissipation between fluids is relatively weak.The difference between the acid breakthrough consumption amount determined by the Darcy equation and the Stokes-Darcy equation is small.When the acid injection rate decreases to an intermediate value，the strength of the acid transport process is reduced，and multiple wormholes begin to appear inside the core sample due to the heterogeneity，as shown in Fig.3b and c.The viscous dissipation between fluids becomes prominent with the emergence of wormholes，increasing the amount of time of acid solution to react with the solid matrix.Fig.4a shows that when the same pore volume of the acid solution is injected into the core sample，the acid solution，which takes into account the viscous dissipation of the fluid phase，is found to have a deeper penetration distance inside the core sample compared to the acid solution only considering the effect of the complex porous structure.Conversely，the deeper penetration distance reduces the acid fluid phase's flowing resistance within the core sample，as shown in Fig.4b，in which the lower pressure difference is observed with the implementation of the Stokes-Darcy equation.The lower flowing resistance further attracts more acid solution to be transported through the wormholes and breakthrough the core sample.The increasing viscous dissipation between fluids results in the acid solution being effectively applied to dissolve the solid matrix instead of flowing out directly from the core outlet boundary.The pore volume to breakthrough calculated by the Stokes-Darcy equation is smaller than the value determined by the Darcy equation due to consideration of the effect of viscous dissipation.When the injection rate is very low，the acid reaction process becomes dominant，and more and more acid is consumed to promote the development of wormholes，as shown in Fig.3d.The increase in the wormhole fraction further enhances the effect of viscous dissipation on the momentum transfer of the acid solution.Unlike the intermediate acid injection rate，the effect of the fluid viscous dissipation is not positive anymore.The viscous dissipation between fluids makes the injected acid solution more inefficiently react with the solid matrix instead of moving forward to break through the target core sample.However，because the Darcy equation does not take into account the effect of fluid viscous dissipation，the acid solution is mistakenly used more efficiently，as shown in Fig.5a.The results show that the dissolution profile decided by the Darcy equation is very close to the breakthrough with the same amount of the injected acid solution.While the dissolution profile determined by the Stokes-Darcy equation shows a limited penetration distance，which also brings higher flowing resistance of the acid fluid phase，as shown in Fig.5b.The pressure differential indicated by the Stokes-Darcy equation between the inlet and outlet of the core sample is greater than that determined by the Darcy equation.A lower value of the acid breakthrough volume is finally gained through the Darcy equation.

Fig.4.Comparison of the numerical result profiles with an intermediate acid injection rate of 0.01 cm/s and acid injection volume at 12 PV(pore volume):(a)Acid concentration profile;(b) Pressure profile.

Fig.5.Comparison of the numerical result profiles with a low acid injection rate of 0.001 cm/s and acid injection volume at 24 PV(pore volume):(a)Acid concentration profile;(b)Pressure profile.

To sum up，the acid solution momentum transfer is mainly affected by the complex porous structure and fluid viscous dissipation during the acidizing process.The inertial effect of the fluid phase has little contribution to the fluid momentum transportation process.The variation of the acid breakthrough volume is mainly caused due to the competition between the influence of the fluid viscous dissipation and the retardation effect of the porous structure.The Stokes-Darcy equation is recommended to describe the fluid phase momentum transfer in the following numerical simulation studies on the carbonate acidizing process.It can take into account the effect of permeable structure as well as the effect of fluid viscous dissipation.

3.3.Effect of acid physical parameters

In this section，numerical studies are designed to discuss the influence of acid solution parameters on the acidizing results，including acid concentration，acid surface reaction rate，and viscosity.

The acid concentration is variable in different experimental studies.Fredd and Fogler (1999) always use an acid injection concentration of 0.5 mol/L to confirm the wormhole formation in porous carbonate media.Sidaoui et al.(2018) use the high concentration of 28% acid to measure the stability of chelating agents for a particular reservoir.The acid system with a mass fraction of 15%has been widely used in recent experimental studies(Aldakkan et al.，2018;Kumar et al.，2020).Wang et al.(1993) discussed the influence of acid concentration on the results through experiments，in which a varying acid injection concentration from 0.147 mol/L to 4.4 mol/L is used to stimulate Indiana limestone core sample with 1-inch diameter and 6-inch length.The results show that acid injection concentration obviously affects the breakthrough consumption amount of acid solution，as shown in Fig.6a.The PVBT of low acid injection concentration is higher than that of high acid injection concentration.A similar conclusion can also be observed in the results of numerical simulations，as shown in Fig.6b.Increasing the acid injection concentration apparently reduces the pore volume to breakthrough.Moreover，the difference of PVBT also decreases with the increase in acid injection concentration.The optimum acid injection rate widely exists in the acidizing curves corresponding to different injection concentrations，and it is less affected by the change of acid injection concentrations.The specific value of the optimum acid injection velocity is around 0.01 cm/s，which is close to the value observed in the experimental studies.The values of 0.02 cm/s and 0.025 cm/s are observed in experiments separately performed by Dong (2012) and Furui et al.(2012).

The acid surface reaction rate is mainly determined by the type of acid solution.Nevertheless，the specific value of acid reaction rate is difficult to determine accurately.Because the acid consumption in the acidizing process contains a heterogeneous chemical reaction，which is not only affected by the chemical reaction occurring on the rock surface but also the migration process of reactants and products between the bulk phase and the rock surface，it leads to the different values of acid surface reaction rate used in the numerical studies.The commonly used value varies from 0.2 cm/s to 0.000018 cm/s(Liu and Mostaghimi，2017b;Panga et al.，2005;Wei et al.，2017).Fig.7 shows pore volume to breakthrough with a wide range of acid surface reaction rates from 0.2 cm/s to 0.00002 cm/s.It can be observed that both the acid breakthrough volume and the optimum acid injection rate are intuitively affected by the acid surface reaction rate.The optimum acid injection rate widely exists in the acidizing curves with different acid surface reaction rates，at which the minimum amount of the acid solution is consumed for achieving a breakthrough.Generally，the value of the optimum acid injection rate increases as the acid reaction rate increases.For example，when the acid reaction rate is at a very low value of 0.00002 cm/s，the optimum acid injection rate is also observed at a low value of around 0.001 cm/s.When the acid reaction rate increases to a high value of 0.2 cm/s，the optimum acid injection rate correspondingly increases to 1 cm/s.Moreover，increasing the acid reaction rate decreases the optimum acid breakthrough volume，as shown in Fig.7.Apart from the reaction rate，the viscosity of the acid solution also changes with different acid systems.In addition，acid solution viscosity is also affected by some chemical additives.For example，polymers are usually added to slow down the corrosion rate between the acid solution and the metal tube.Fig.8 studies the effect of viscosity of Newtonian fluid on the stimulation results.The results indicate that the influence of acid viscosity on the acidizing process is limited compared with the influence of acid surface reaction rate.The viscosity mainly affects the acidizing process when a very small rate of injection rate is adopted.High acid viscosity increases the amount of acid consumption for the breakthrough.The low acid injection rate allows the acid solution to gain sufficient time to dissolve the solid matrix，which leads to the full development of the wormhole area.The continuous increase in the wormhole area makes the effect of the acid solution's viscosity dissipation continue to increase.Due to viscous dissipation，more acid is needed to dissolve the solid matrix.

Fig.6.Comparison of acidizing curves with different concentrations:(a) Experimental study (Wang et al.，1993);(b) Numerical study in current work.

Fig.7.The acidizing curve of pore volume to breakthrough with different acid surface reaction rates.

3.4.Effect of rock physical parameters

The heterogeneity of carbonate rock is more obvious than sandstones due to complex geological processes.The rock physical parameters of carbonate cores vary with the heterogeneity of the pore structure.This section focuses on the impact of rock physical characteristics on acidizing results，such as porosity，permeability，and heterogeneity.

Fig.8.The acidizing curve of pore volume to breakthrough with different viscosities.

Jia et al.(2021c) comprehensively collect and compares carbonate reservoirs' porosity and permeability data in the literature.According to the findings，the porosity and permeability of acidifying core samples are typically between 0.05 and 0.3 and 1 mD to 100 mD，respectively.Because the definition of pore volume to breakthrough is related to porosity，the PVBT is not used to discuss the influence of porosity on the stimulation results (Jia et al.，2021c).Instead，the core volume to breakthrough (CVBT) is employed to discuss the porosity influence，as shown in Fig.9a，defined as the multiple pore volume to breakthrough and porosity.The results show that the influence of the porosity on the acidizing process differs from the change of acid injection rate.When the injection rate is high，the acid breakthrough consumption amount increases with the increase in porosity.When the rate of injection is low，the consumption of acid solution decreases with the increase in porosity.With a high rate of acid injection，the convection transportation process of the acid solution is so strong that most of the acid solution tends to flow through the core sample without sufficient reaction with the solid matrix.Therefore，more acid solution is needed for the core sample with high porosity.Because higher porosity contributes to more storage space.When the rate of injection is extremely low，compact dissolution prefers to form because the acid solution spends enough time dissolving the solid matrix.Hence，the acid breakthrough volume mainly depends on the fraction of the solid matrix.Compared to low-porosity core samples，it is relatively easy for high-porosity core samples to achieve a breakthrough due to a smaller fraction of the solid matrix.Fig.9b further shows the acidizing curve of pore volume to breakthrough.Unlike the change of the core volume to breakthrough，the pore volume to breakthrough always decreases as the rock porosity increases.The results also show that it is necessary to use the definition of core volume to breakthrough to study the effect of the rock porosity on the acidizing results.

Fig.9.The acidizing curve of core volume to breakthrough with different porosity values:(a) Core volume to breakthrough (CVBT);(b) Pore volume to breakthrough (PVBT).

The influence of permeability on the acidizing results basically shows a negative correlation with the acid injection rate.The pore volume to breakthrough decreases as the permeability increases，as shown in Fig.10.The difference between the acidizing curves of different permeability gradually becomes insignificant with the decrease in the injection rate.The consumption of acid solution is mainly controlled by the transport process and reaction process.The high injection rate enhances the influence of the acid transport process.Moreover，higher porosity contributes to more storage space resulting in more acid requirement.The momentum transfer of acid solution is mainly affected by the porous structure because the wormholes are not well developed.Therefore，the amount of acid consumed is very sensitive to the change of permeability.The momentum consumption of the acid solution decreases with the increase in permeability during the acidizing process.When the injection rate is low，the reaction process is highlighted，in which the acid solution is mainly used to dissolve the solid matrix.Low porosity requires more acid solution.The developed wormholes enhance the effect of fluid viscosity dissipation on the momentum transfer of the acid solution and weaken the influence of the porous structure.Therefore，the acidizing results become insignificant in response to changes in permeability.

Fig.10.The acidizing curve of pore volume to breakthrough with different permeability values.

Fig.11.The acidizing curve of pore volume to breakthrough with different heterogeneity magnitudes.

Fig.12.The comparison of pore volume to breakthrough:(a) Isolated fracture and matrix;(b) Isolated vug and matrix.

Fig.13.The acidizing dissolution patterns of isolated fracture and isolated vug with different injection rates:(a) 1 cm/s;(b) 0.1 cm/s;(c) 0.01 cm/s;(d) 0.001 cm/s.

Next，we continue to investigate the impact of heterogeneity on the acidizing process.Fig.11 shows the acidizing curves of the PVBT with heterogeneity magnitudes of 0.01，0.05，0.1，and 0.15.The results show that the PVBT does not change significantly when the heterogeneity magnitude is less than 0.1.A similar conclusion is also observed in a numerical study performed by Kalia and Balakotaiah (2009).In addition to matrix heterogeneity，fracture and vug are also important parts of carbonate heterogeneity(Mohammed et al.，2021).Isolated fracture and vug carbonate reserves are also potential targets for acidizing treatment.The contribution of the isolated fracture and vug to the rock permeability is limited.Here，we mainly study the acidizing process in two types of typical isolated fracture and vug carbonate core samples，in which the isolated fracture has a length of 5 cm and an aperture of 2 mm (Aghli et al.，2020)，and the isolated vug is a square with a length of 1.5 cm (Zhang et al.，2005).Both isolated fracture and vug are located in the middle of the core sample.Fig.12 shows the acidizing curves of the isolated fracture and isolated vug samples.The results show that the appearance of isolated fracture and vug obviously affects the pore volume to breakthrough.The consumption volume of the acid solution decreases significantly within the isolated fracture and vug core samples.However，the isolated fracture and vug has an unobvious impact on core porosity and permeability.The porosity of fracture and vug core samples increases from 15% of the matrix core to 16.875% and 20.22%，respectively.The permeability separately increases from 5 mD of the matrix core to 8.83 mD for an isolated fracture core sample and 5.55 mD for an isolated vug core sample.Simultaneously，the isolated fracture and vug affect acidizing dissolution patterns，as shown in Fig.13.Isolated fracture and vug in the core sample offer a favorable route for acid solution migration，which makes the acid solution more effectively stimulated the core sample for the breakthrough.For such cores，the breakthrough can be relatively easy to achieve.Because compared with matrix core breakthrough，acid is mainly used to communicate the advantageous channel of isolated fracture and isolate vug with the inlet and outlet boundaries of the core sample.

4.Conclusions

The continuum model has undergone significant development to study the acidizing process due to the effective reduction of the calculation expense and high simulation accuracy.In this work，the Darcy equation，the Stokes-Darcy equation，and the Naiver-Stokes-Darcy equation are separately used as the fluid phase momentum transfer equation to study the effect of the different momentum transfer mechanisms on the acidizing results，including the retardation effect of the permeable structure，fluid viscous dissipation，and inertial effect.The results indicate that the momentum transfer of the acid solution is not only affected by the complex permeable structure of the core sample，but with the development of the wormholes regions，the interactions between the acid solution molecules are gradually promoted and further affect the acid solution transportation process.Finally，the Stokes-Darcy equation is recommended as the general fluid phase momentum transfer equation to calculate fluid phase velocity in the numerical simulation studies of the carbonate acidizing process because it can consider both the effect of porous media and the effect of fluid viscous dissipation.In addition，a parametric sensitivity analysis is conducted to discuss the influence of the acid and rock physical property parameters on the acidizing process，including acid concentration，acid surface reaction rate，viscosity，core porosity，permeability，and heterogeneity.The conclusions are presented below:

1.The optimum acid injection rate is generally present in the acidizing curves corresponding to different injection concentrations，and it is less affected by the varying acid injection concentrations.

2.The acid injection concentration affects the breakthrough consumption amount of acid solution.Increasing the acid injection concentration apparently reduces the pore volume to breakthrough(PVBT).The difference in PVBT also decreases with the increase in acid concentration.

3.The influence of acid viscosity on the acidizing process is limited compared with the acid surface reaction rate.The acid breakthrough volume increases as the acid reaction rate decreases.The optimum acid injection rate decreases with the reduction of the acid surface reaction rate.

4.When the injection rate is high，the acid breakthrough consumption amount increases with the increase in porosity.Conversely，when the rate of injection is low，the consumption of acid solution decreases with increased porosity.

5.The influence of permeability on the acidizing results basically shows a negative correlation with the injection rate.The difference between the acidizing curves of different permeability gradually becomes insignificant with the decrease in the injection rate.

6.Compared with the matrix heterogeneity，the appearance of isolated fracture and vug obviously affects the acid breakthrough volume.The consumption volume of the acid solution decreases significantly within the isolated fracture and vug core samples.

Acknowledgments

We appreciate the financial support from the Key Project of the National Natural Science Foundation of China (No.52034010) and the China Scholarship Council(201906450038).The authors would also like to thank Dr.Chowdhury Mamun for evaluating the paper and making recommendations to enhance it.