Wormholing Influenced by Injection Temperature in Carbonate Rocks

This paper presented empirical models of describing reaction rate vs. hydrochloric acid temperature and concentration by regressing experimental data. And this paper also introduced the dependent reaction heat model into the thermal non-equilibrium models and coupled with two-scale continuum model to obtain governing equations for describing wormholing under non-isothermal conditions. The governing equations were discretized by implicit difference method and solved by programing. The effects of temperature on wormholing have been investigated based on the simulation results on 2-D vision. A significant difference of the effluent temperature between the dependent reaction heat model in present work and the constant reaction heat model in available literatures was observed, especially in high injection rate and strong acid concentration. In addition, the tendencies of optimum injection capacity vs. injection temperature under isothermal and non-isothermal conditions were almost different. Finally, an optimum injection temperature was found by changing the injection temperature under non-isothermal conditions.


Introduction
Matrix acidizing is a most important method to stimulate carbonate reservoir through bypassing the damage zone with wormholes using a small volume of acid.Therefore, nearly all of the studies focusing on recovering or stimulating How to cite this paper: Xue, H., Huang, Z.X., Zhao, L.Q., Wang the production of carbonate reservoirs were based on investigations of the initiation and propagation of wormholes or unstable dissolution fronts.Many mathematical models were developed to describe the wormhole patterns observed in experiment results [1] [2] [3] [4]: the capillary tube model [5] [6] [7]; the network model [8]; the averaged model [4] [9]; the fractal model [1]  (XCMT) [13].
The wormholing process is very complicated since there are many known factors influencing convection, dispersion and reaction properties, such as acid system, concentration of acid, porosity or permeability distribution, fracture or vug medium, injection conditions, acid flow types and temperature.In the physical heterogeneities, Kalia and Balakotaiah (2009) used a unique heterogeneity parameter to characterize porosity heterogeneity, and found that the dissolution depended on the orientation of fractures [14]; Liu et al. (2012) have studied the similar problem with local non-equilibrium conditions Open Journal of Yangtze Gas and Oil [25] and the adiabatic boundary and isothermal boundary in wormholing were discussed under radial flow [23] [26].It demonstrated that the temperature also played a leading role in governing fluid flow as dispersion, convection and reaction did.The reaction heat is a constant as treated in the previous work while it is confirmed that it will be affected by temperature and pressure in matrix and fracture in carbonate reservoirs [27].In order to simulate temperature profile of core more accurately, the dependent reaction heat model was introduced into two-scale continuum model to simulate wormhole initiation and propagation under different injection conditions.The simulation results calculated from both reaction heat models were compared and discussed.To understand effects of injection temperatures on wormholing, we have studied acid volume of breakthrough, optimum injection capacity and overall pressure drops under different acid injection temperatures under isothermal and non-isothermal conditions using this simulation method.

Darcy Flow Model
The fluid flow in porous medium conforms to the equivalent Darcy's law given by Ratnakar et al. ( 2013) for linear flow [22]: where .represents the norm of a vector/matrix.U is the velocity (m/s).P is the formation pressure (MPa).K is the permeability of the reservoir (10 −3 × μm 2 ).μ eff is the effective viscosity (mPa•s).
Equation (2) shows the effective viscosity model of acid fluid: where φ is the porosity (Dimensionless).n is the power index.K ten is the permeability tensor.μ 0 is the apparent viscosity (mPa•s).It should be noted that Darcy law is adapted to the flow of Newtonian acid when the parameter, n, in the effective viscosity term equals to 1. While, it will be adapted to non-Newtonian acid flow when n is less than 1.In this paper, we assume that HCl acid system is the Newtonian fluid.And the Mass balance or continuity model of acid phase is given by:

Darcy-Scale Model
In order to keep track of the concentration of hydrogen ion, the model to describe the hydrogen ion using species balance is shown as follows [10]: where D e is the effective dispersion coefficient in the acid phase (m 2 /s).v a is surface area per unit volume in the core available for reaction (m 2 /m 3 ).C f is the cup-mixing mass concentration of the acid in the fluid phase (kmol/m 3 ).R(C s ) is the rate of the dissolution reaction, for single step irreversible reaction ( ) ( ) α is the dissolving power of the acid, defined as mass of solid dissolved by unit mass of reacted acid (kg/kg).ρ s is the density of the solid phase (kg/m 3 ).

Partial Liquid Equilibrium
The reaction between carbonate rock and aqueous of HCl acid system involves the procedures, the transport of the acid molecules from the bulk fluid to the rock surface, and the chemical reaction at the rock surface.And the overall reaction is controlled by the slower procedure.Almost all the previous work of simulating wormhole propagation in carbonate rocks uses hydrogen ion concentration particularly to describe chemical reaction from the equilibrium given by: ( ) ( ) While in this high TDS (Total Dissolved Solids) system, we have to use hydrogen ion activity in equilibrium conditions.Alkattan et al. (1998) gave the developed function of overall calcite dissolution rate of hydrogen ion activity as follows [28].
( ) where k c is the mass transfer coefficient (m/s).k s is the reaction rate constant (m/s).

H ,s
γ + is the activity coefficient of hydrogen ion.
, the overall calcite dissolution rate is kinetically controlled, then ( ) rate is mass transfer controlled, then ( ) . The activity coefficient of hydrogen ion ( H ,s γ + ) in the process of calcite and limestone dissolution can be computed from: The designate parameters f γ , m Cl , B HCl , E and C HCl in Equation ( 8) are given by Alkattan et al. (1998) [28].

Pore Scale Model
The dissolution process determined local quantities in porous media, therefore the modified Garman-Kozeny correlation models are introduced to describe the relationship between porosity, permeability, interfacial area and pore radius.
where a 0 is the initial interfacial area of the medium (m 2 /m 3 ).K 0 is the initial permeability of the reservoir (10 −3 × μm 2 ).r p is the pore radius of the medium (m).r p0 is the initial pore radius of the medium (m).β is the exponent determined by the experiment (Dimensionless).0 φ is the porosity of the reservoir (Dimensionless).
Panga et al. (2005) proposed the formulas to describe the relationship between Sherwood number and local mass transfer coefficient [10], k c .and the diffusion coefficient determined by molecular diffusion and convection diffusion in x direction and y direction respectively as are listed as Equation (11) and Equation (12).
where Sh is the Sherwood number, which is defined as the ratio of convective to diffusive mass transport (Dimensionless).Sh ∞ is the asymptotic Sherwood number (Dimensionless).m is the ratio of pore length to pore diameter (Dimensionless).Re p is the pore scale Reynold's number (Dimensionless).Sc is the Schmidt number (Dimensionless).|U| is the magnitude of v and u (m/s).λ y is the constant depending on pore geometry (λ y = 0.1 for packed-bed of spheres).
λ x is the constant depending on pore geometry (λ x = 0.5 for packed-bed of spheres).a os is the constant depending on pore geometry (a os = 1.0 for packed-bed of spheres).D m is the effective molecular diffusivity of acid (m 2 /s).D ex is the effective dispersion coefficient in the acid phase in the x direction (m 2 /s).D ey is the effective dispersion coefficient in the acid phase in the y direction (m 2 /s).

Energy Model
In order to study wormhole propagation in carbonate matrix acidizing affected by temperature, the models considering reaction heat, which are used to describe a thermal non-equilibrium between the solid and acid phases, are introduced.
We assumed that the exothermic heat of reaction between acid and rock was transported from the solid phase to the acid phase, and the reaction rate based Open Journal of Yangtze Gas and Oil on the temperature is mainly determined by the solid temperature, T s , since the reaction happened at the surface of the rock.
where T f is the temperature of the acid phase (K); T s is the temperature of the solid phase (K); C Pf is the heat capacity of acid (J/kg•K); C Ps is the heat capacity of solid (J/kg•K); k ef is the thermal conductivity of acid (W/m•K); k es is the thermal conductivity of solid (W/m•K); h c is the acid to solid heat transfer coefficient (W/m 2 K); ΔH r (T s , 1 atm) is the acid-rock molar reaction enthalpy (J/mol).As for the reaction heat depending on the temperature and pressure in porous media [29], the standard molar reaction enthalpies of HCl acid and limestone at porous surface is given by: ( where ΔH r (298.15K) is the molar reaction enthalpy of reaction at 298.15 K (J/mol); C p,m is the constant-pressure molar heat capacity (J/K mol); The standard molar reaction enthalpy (obtained from the Handbook of Physical Chemistry [30]) of each parameter in Equation ( 15) is listed in Table 1.

Reaction Kinetics
The chemical reaction kinetics of HCl acid with various concentrations and temperatures on MISSAN limestone from Iraq has been studied in this paper with rotating disk apparatus at a constant pressure and shear rate.And the functions are used to calculate the reaction rate with the changing temperature and the consumption of hydrochloric acid, which are obtained by regressing experimental data.
Since the reaction rate decreases with the consumption of HCl acid, the function describing the relationship between the concentration of HCl acid and the reaction rate is obtained by regressing laboratory test data shown in Figure 1.

Pressure boundary condition
It is assumed that the total injection capacity is constant at the inlet of the carbonate core during acidizing.However, the injection capacity for each grid block changes with time because of the permeability improvement.Li (2004) listed out the way to solve the constant flux boundary [31].

Boundary condition of HCl concentration
Boundary condition of acid temperature Boundary condition of solid temperature where Q inj is the injection capacity (m 3 /min).Q i is the injection capacity at each Open Journal of Yangtze Gas and Oil grid block at the inlet of core (m 3 /min).P e is the constant pressure at the outlet of core (MPa). 0f C is the mass concentration of the acid at the inlet, (kmol/m 3 ).

s
T is the initial solid temperature (K).

Computational Methods
The pressure, temperature and HCl concentration for each core segment at each time step are solved numerically by using boundary and initial conditions in Equations ( 20)- (23).The partial differential Equations ((3), ( 4), ( 13) and ( 14)) are performed with spatial and time discretization through implicit finite difference method to obtain the corresponding ODEs.And we programed an algorithm to save the coefficients of ODEs into a matrix by using MATLAB.The results of pressure, temperature and HCl concentration in each time step are calculated using Ax = B, where x is the vector to be calculated.The detailed solving procedures are listed as: 1) initialize the mesh of core and generate primary distribution of porosity through the method developed by Liu et al. (2012) [15]; 2) calculate the pressure profile and velocity distribution in the core at time, n + 1, by Equations ( 1)-(3); 3) calculate the molar reaction enthalpy of reaction, ΔH r (T s , 1 atm), in each core segment at time, n, according to Equation ( 16); 4) substitute ΔH r (T s , 1 atm) into Equation ( 14), and calculate the temperature profile of acid phase and solid phase by solving the integration of Equation ( 13) and Equation (14) using coupled method; 5) calculate the reaction rate constant, k s , the mass transfer coefficient, k c , and molecular diffusivity, D e , for each core segment; 6) calculate the HCl concentration profile at time, n + 1, through Equation (4); 7) solve Equation (5) and Equation ( 9) in sequence, the porosity, permeability, pore radius and interfacial area can be obtained; 8) judge: if the acid breaks the core (when the inlet pressure drops to 1% of its initial value, breakthrough phenomenon is considered [14]), then end the calculation; else go to step 2).

Analysis of Simulation Results
The temperature is one of the important factors which influence the wormhole pattern in matrix acidizing in carbonate reservoirs, thereby affecting the optimum injection capacity and the acid volume to break the damage zone of formation [32] [33] [34] [35].Notably, there are no reports that have simulated the wormhole propagation process under different injection temperatures of acid in carbonate reservoirs.To understand the injection temperatures on wormholing, we focus on studying the acid volume of breakthrough and optimum injection capacity under different temperatures of injected acid using the simulation method with the 2-D vision.The analysis of simulation results is based on the conditions listed in Table 2.

Effects of Acid-Rock Reaction Heat on Effluent Temperature
A series of studies have shown that the reaction heat of acid-rock should be considered when simulating acidizing progress [23] [24] [26].All of them assumed acid-molar reaction heat as a constant value while it will change with the temperature and pressure as illustrated in Equation (18).The difference of reaction

H. Xue et al. Open Journal of Yangtze Gas and Oil
heat treating methods will affect the predicted temperature to some extent, therefore affecting reaction rate.In order to compare the available model (with constant reaction heat) and present theoretical model (with dependent reaction heat), the temperature in the effluent of core changing with time under different injection rates and acid concentrations were analyzed.The values of the constant parameters in calculation are given by: T inj = 298 K, 0 s 365 K T = , ΔH r = 4.86 kJ/mol [23].
Figure 6 illustrates the effluent temperature calculated by constant and dependent reaction heat models under different injection rates: 0.2 m 3 /min, 1 m 3 /min and 10 m 3 /min.The curves indicate that the effluent temperature declines from 365 K to about 298 K under both reaction heat treating methods in different injection rates.This is because the colder acid supplying at the inlet of core persistently replaces the heated acid in the core.However, the temperature decreases faster with the increase of injection rate.It is contributed to two points.
In the first place, the velocity of cold supplying acid is faster in the higher injection rate, leading to a faster replacing rate of heated acid.The second reason is that the high injection rate increases the heat transfer rate between acid and rock.
The trends of temperature are the same under both reaction heat treating models.Nevertheless, the temperatures treated by dependent reaction heat model are much higher than that of the constant reaction heat model by comparing the line-graph and dash-graph in Figure 6.Moreover, raising the injection rate from 0.2 m 3 /min to 1 m 3 /min and 10 m 3 /min, the maximum differences of temperature for both simulation results experience an increase from 1.32˚C to 8.12˚C and 13.57˚C as shown in stripe-graph in Figure 6. Figure 7 shows the variation of effluent temperature declines solely caused by persistently injecting cold acid at the inlet of core under different acid mass Open Journal of Yangtze Gas and Oil concentrations, 5wt%, 15wt% and 30wt%.When the constant reaction heat model is chosen, the effluent temperature under strong acid concentration is a little higher than that of the weak one in the whole acidizing process as shown in the line-graph of Figure 7.It is caused by the phenomenon that the reaction heat will increase with the raising acid concentration.And Guo et al. ( 2014) also achieved the same conclusion in the study of acid fracturing [27].When it comes to the dash-graph drawn by the simulation results from the dependent reaction heat model, the trends and conclusions are almost the same as the constant reaction heat model did before.However, the difference of effluent temperature between strong acid concentration and weak acid concentration is more remarkable than that of the line-graph.Furthermore, with the increase of the acid concentration from 5wt%, to 15wt% and 30wt%, the maximum differences of temperature for both treating models catch an increase from 1.12˚C to 6.17˚C and 12.16˚C as shown in the stripe-graph of Figure 7.

Effects of Injection Temperature on Pore Volume of Breakthrough
Under the isothermal conditions, we assume that the injected acid temperature equals to the solid temperature, and they are uniform throughout the core medium.Therefore, the concentration of HCl acid is the only driving mechanism for wormhole propagation under the certain initial temperature of acid/solid.The wormhole propagation process is more complicated when the non-isothermal conditions are considered.The reaction rate in each core segment at each time step is influenced not only by the acid concentration, but also by the temperature significantly, which will affect the volume of acid breakthrough and optimum injection capacity. Figure 9 plots the acid efficiency curves with different injection temperatures under non-isothermal conditions.Different from the isothermal conditions shown in Figure 8, there is no such tendency that the optimum injection capacity increases with increasing injection temperature, al-though the optimum injection capacity exists in each group of simulation results.
And they obtain a nearly similar value for the optimum injection capacity, 1.5 ml/min.However, the injection temperature affects PV BT under any injection capacity, especially under low injection capacity.As we can see from the curves, the maximum difference of PV BT is 46.67 PV which occurs at 0.2 ml/min, and the minimum value is 0.31 PV at 1.5 ml/min.

Conclusions
In the present work, the dependent reaction heat model was developed and introduced into the integrated governing models of matrix acidizing in linear core.
And the effluent temperatures were compared on both reaction heat treating methods.Finally, the acid wormholing influenced by temperature under isothermal and non-isothermal conditions was studied by this simulation method.
The conclusions coming from the calculation and analysis could be summarized as: 1) Since there is an important difference of simulation results between the dependent reaction heat model and constant reaction model, and the former takes more comprehensive factors into account and more conforms to reality, therefore, it is suggested to use the dependent reaction heat model when modeling the acidizing process; 2) Under isothermal conditions, the optimum injection capacity decreases with the decreasing temperature; 3) Under non-isothermal conditions, the wormhole patterns, i.e. uniform dissolution, ramified wormhole, dominant wormhole, conical wormhole and face dissolution, can also be generated by changing the injection capacity, and the optimum injection capacity is slightly affected by the injection temperature; 4) Similar to the optimum injection capacity, there also exists an optimum injection temperature for developing dominant wormhole under non-isothermal conditions.
In the present work, the influence of supercritical CO 2 has been ignored, Open Journal of Yangtze Gas and Oil which will affect the simulation results to some extent.When the supercritical CO 2 is considered, it will affect the reaction heat definitely, and thereby affect the wormholing patterns.By the way, the two-phase flow should be considered in the modeling instead of one-phase flow in the present work when the supercritical CO 2 is considered.
The wormhole propagation process in acidizing carbonate is very complex, especially when using the unconventional acid system such as in-situ self-diverting acid and in-situ cross-linked acid etc.It requires further studies to capture the details of wormholing of unconventional acid system by considering the affection of temperature.The work presented in this paper can be improved by including the effects of supercritical fluid CO 2 ; the fracture-vug-matrix medium and the non-Newtonian acid.In the future studies, some of the extension work will be pursued.
, H.H. and Liu, P.L. (2019) Wormholing Influenced by Injection Temperature in Carbonate Rocks.Open Journal of Yangtze Gas and Oil, 4, 12-30.https://doi.org/10.4236/ojogas.2019.41002H. Xue et al.DOI: 10.4236/ojogas.2019.4100213 Open Journal of Yangtze Gas and Oil [2].Panga et al.(2005)  coupled Darcy-scale dissolution model with the pore-scale model to summarize the governing equations which is called the two-scale continuum model[10].In this model, the macroscopic properties such as mass transfer coefficients and dispersion properties were properly considered.Since the results calculated from the two-scale continuum model matched experiment results very well, most of the extension simulation studies were based on them.Kalia and Balakotaiah (2007) studied dissolution patterns in radial flow by extending the linear model to radial model [11].Maheshwari et al. (2013) analyzed the 3-D dissolution patterns [12].Smith et al. (2013) simulated core CO 2 flooding of low permeability carbonates coupled with X-ray computed microtomography (2012) introduced a normal distribution instead of a uniform distribution of porosity to study the wormholing in radial flow, and they also claimed that only a large number of continuous vugs could affect wormholing [15]; Izgec et al. (2010) further addressed the effects of vugs on matrix acidizing based on experiments [16]; Cohen et al. (2008) motivated a thorough study of geometry effects on wormholingfrom pore scale to wellbore scale in 3-D vision[17].And in the acid system, most of the studies focus on diverting acid because of extensive applying in both sandstone and carbonate heterogeneous reservoirs.The simulation results of diverting effects achieved a high agreement with experimental results[18] [19][20] [21][22].The acid types (HAc, CDTA, EDTA and HCl acid system) affecting wormhole patterns, pressure drops, PVBT and wormhole density also achieved an extensively discussion through both laboratory experiments and numerical simulations[15] [18]  [22][23].In addition to the work mentioned above, almost all of them ignored the effects of temperature on wormholing in matrix acidizing simulation.However, temperature is also an important factor influencing wormhole patterns formation.Kalia et al. (2009)  have introduced the temperature model into two-scale continuum model and compared the wormhole patterns under isothermal and non-isothermal conditions by changing the viscosity of acid according to temperature under liner flow[24].Based on the work of Kalia et al.,Bousri et al.

Figure 3 Figure 3 .
Figure 3 illustrates the wormhole patterns under different injection capacities with 293.15K temperature of injected acid, and 393.15K initial temperature of

Figure 4 and
Figure4and Figure5illustrate the corresponding temperature profiles of acid phase and solid phase when the acid breaks the core shown in Figure3.With the colder acid entering the core, the solid phase is cooled down through the thermal conduction between acid phase and solid phase.Under the low injection capacity, the temperature of solid phase decreases more significantly than that of the moderate and high injection capacity by comparing Figures5(a)-(c).It might be because the breakthrough volume of injected acid is much more under conical wormhole condition or face dissolution condition than that under dominant wormhole condition or uniform dissolution condition.By comparing the temperature in the wormholes for acid phase and solid phase, the temperature along the wormholes near the inlet of core nearly equals to the injection temperature, 293.15 K, and the temperature along the wormholes in the acid phase is slightly lower than the solid temperature because of the reaction on the solid surface.It is interesting to note that the maximum temperature is found at the tips of wormholes since the majority of the reaction occurs here, which agrees with the modeling results of Kalia et al. (2009)[24].

Figure 8
Figure 8 plots the acid volume required to break the core for various acid injection capacities under 293.15K, 338.15K and 393.15 K.As the acid efficiency curves show, no matter what the isothermal conditions are, the PV BT decreasesas the injection capacity increases when it is smaller than the optimum injection capacity.However, it experiences a relatively slow growth with the increase of injection capacity when it exceeds the optimum injection capacity.On the other hand, the optimum injection capacity increases with the increase of initial temperature of acid/solid, and the corresponding optimum injection capacities are 0.7 ml/min, 1.2 ml/min and 1.7 ml/min for the temperatures 293.15 K, 338.15K and 393.15 K, respectively.This is because that reaction rate decreases with the decreasing temperature from Equation (18), and the injection capacity should be lowered to make sure the convection and reaction become comparable to generate the dominant wormhole.

Figure 8 .
Figure 8. Volume of acid required to break the core for various acid injection capacity under isothermal conditions@ 293.15 K, 338.15K and 393.15 K.

Figure 10 Figure 9 .
Figure10shows the pressure drop as a function of time for non-isothermal cases with different injection temperatures using the optimum injection capacity obtained from Figure9.All of the pressure drops decrease almost linearly with time until the acid breaks the core, which means the dominant wormholes are generated in those cases[36].On the other hand, the pressure drop decreases more slowly with increasing injection temperature under optimum injection

Figure 10 .
Figure 10.Pressure drop for non-isothermal conditions using optimum injection capacity.

Figure 11 .
Figure 11.Volumes of acid breakthrough for non-isothermal conditions using optimum injection capacity.