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 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     : the capillary tube model    ; the network model  ; the averaged model   ; the fractal model   . 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  . 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  . Maheshwari et al. (2013) analyzed the 3-D dissolution patterns  . Smith et al. (2013) simulated core CO2 flooding of low permeability carbonates coupled with X-ray computed microtomography (XCMT)  .
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  ; Liu et al. (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  ; Izgec et al. (2010) further addressed the effects of vugs on matrix acidizing based on experiments  ; Cohen et al. (2008) motivated a thorough study of geometry effects on wormholing from pore scale to wellbore scale in 3-D vision  . 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      . 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     .
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  . Based on the work of Kalia et al., Bousri et al. (2012) have studied the similar problem with local non-equilibrium conditions  and the adiabatic boundary and isothermal boundary in wormholing were discussed under radial flow   . 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  . 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.
2. Mathematical Models
2.1. 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  :
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 × μm2). μ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. Kten 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:
2.2. 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  :
where De is the effective dispersion coefficient in the acid phase (m2/s). is surface area per unit volume in the core available for reaction (m2/m3). Cf is the cup-mixing mass concentration of the acid in the fluid phase (kmol/m3). R(Cs) is the rate of the dissolution reaction, for single step irreversible reaction (m・kmol/s m3). And the Mass balance model for calcium carbonate gives:
α 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/m3).
2.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  .
where kc is the mass transfer coefficient (m/s). ks is the reaction rate constant (m/s). is the activity coefficient of hydrogen ion.
When , the overall calcite dissolution rate is kinetically controlled, then . But when , the overall calcite dissolution rate is mass transfer controlled, then . The activity coefficient of hydrogen ion ( ) in the process of calcite and limestone dissolution can be computed from:
The designate parameters fγ, mCl, BHCl, E and CHCl in Equation (8) are given by Alkattan et al. (1998)  .
2.4. 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 a0 is the initial interfacial area of the medium (m2/m3). K0 is the initial permeability of the reservoir (10−3 × μm2). rp is the pore radius of the medium (m). rp0 is the initial pore radius of the medium (m). β is the exponent determined by the experiment (Dimensionless). 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  , kc. 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). Rep 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). aos is the constant depending on pore geometry (aos = 1.0 for packed-bed of spheres). Dm is the effective molecular diffusivity of acid (m2/s). Dex is the effective dispersion coefficient in the acid phase in the x direction (m2/s). Dey is the effective dispersion coefficient in the acid phase in the y direction (m2/s).
2.5. 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 on the temperature is mainly determined by the solid temperature, Ts, since the reaction happened at the surface of the rock.
where Tf is the temperature of the acid phase (K); Ts is the temperature of the solid phase (K); CPf is the heat capacity of acid (J/kg・K); CPs is the heat capacity of solid (J/kg・K); kef is the thermal conductivity of acid (W/m・K); kes is the thermal conductivity of solid (W/m・K); hc is the acid to solid heat transfer coefficient (W/m2 K); ΔHr (Ts, 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  , the standard molar reaction enthalpies of HCl acid and limestone at porous surface is given by:
where ΔHr (298.15 K) is the molar reaction enthalpy of reaction at 298.15 K (J/mol); Cp,m is the constant-pressure molar heat capacity (J/K mol);
The standard molar reaction enthalpy (obtained from the Handbook of Physical Chemistry  ) of each parameter in Equation (15) is listed in Table 1.
Table 1. Parameters of thermodynamic in Equation (15).
The simplified expression of ΔHr (Ts, 1 atm) can be achieved through substituting values of Table 1 into Equation (15).
2.6. 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.
Figure 1. Relationship of concentration and reaction rate@ 100˚C; 5 MPa; 500 r/min.
The temperature is another factor which affects the reaction rate significantly. In the acid-rock reaction process, the reaction and heat transfer between acid and solid make the changes of reaction rate very complicated. By regressing lab tests on the temperature and reaction rate shown in Figure 2 and integrating
Figure 2. Relationship of temperature and reaction rate@ 20% wt HCl acid; 5 MPa; 500 r/min.
Equation (17), an approximation function describing the reaction rate vs. acid concentration and temperature is obtained as
2.7. Initial and Boundary Conditions
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  .
Boundary condition of HCl concentration
Boundary condition of acid temperature
Boundary condition of solid temperature
where Qinj is the injection capacity (m3/min). Qi is the injection capacity at each grid block at the inlet of core (m3/min). Pe is the constant pressure at the outlet of core (MPa). is the mass concentration of the acid at the inlet, (kmol/m3). is the initial solid temperature (K).
3. 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)  ; 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, ΔHr (Ts, 1 atm), in each core segment at time, n, according to Equation (16); 4) substitute ΔHr (Ts, 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, ks, the mass transfer coefficient, kc, and molecular diffusivity, De, 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  ), then end the calculation; else go to step 2).
4. 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     . 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.
4.1. Analysis of Temperature Profiles under Non-Isothermal Conditions
Figure 3 illustrates the wormhole patterns under different injection capacities with 293.15 K temperature of injected acid, and 393.15 K initial temperature of
Table 2. The corresponding values of parameters in simulation.
Figure 3. Wormhole patterns under non-isothermal conditions using different injection capacity@ (a) Q = 0.2 ml/min; (b) Q = 1.5 ml/min; (c) Q = 10 ml/min; Tinj = 293.15 K; .
solid phase. Comparing to the wormhole patterns under the isothermal conditions proposed by Maheshwari et al. (2013)  , the typical wormhole patterns can also be observed under the non-isothermal conditions as shown in Figure 3. The conical wormhole, dominant wormhole and uniform dissolution are generated under the corresponding injection capacity 0.2 ml/min, 1.5 ml/min and 10 ml/min in Figures 3(a)-(c).
Figure 4 and Figure 5 illustrate the corresponding temperature profiles of acid phase and solid phase when the acid breaks the core shown in Figure 3. 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 Figures 5(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)  .
Figure 4. Temperature profiles of acid phase when the acid breaks the core@ (a) Q = 0.2 ml/min; (b) Q = 1.5 ml/min; (c) Q = 10 ml/min; Tinj = 293.15 K; .
Figure 5. Temperature profiles of solid phase when the acid breaks the core@ (a) Q = 0.2 ml/min; (b) Q = 1.5 ml/min; (c) Q = 10 ml/min; Tinj = 293.15 K; .
4.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    . 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 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: Tinj = 298 K, , ΔHr = 4.86 kJ/mol  .
Figure 6 illustrates the effluent temperature calculated by constant and dependent reaction heat models under different injection rates: 0.2 m3/min, 1 m3/min and 10 m3/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 m3/min to 1 m3/min and 10 m3/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 6. Effluent temperature of different reaction heat models under varying injection rate@ 0.2 m3/min, 1 m3/min and 10 m3/min.
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 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  . 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.
Figure 7. Effluent temperature of different reaction heat models under varying acid mass concentration@ 5wt%, 15wt% and 30wt%.
According to the discussions above, an obvious difference was observed by the simulation results of effluent temperature through both treating models, especially in high injection rate and strong acid concentration. This is mainly caused by the phenomenon that the high injection rate will stimulate the molar heat transfer rate, and the strong acid concentration will make the reaction heat increase. The difference of effluent temperature indirectly proves that the temperature profiles in core are almost different, therefore affecting reaction rate, wormhole propagating velocity, wormhole patterns in the matrix acidizing of carbonate reservoir. Moreover, the dependent reaction heat model takes more comprehensive factors into account and more conforms to reality, therefore, it is of great importance to consider to apply the dependent reaction heat model in the future work referring to the matrix acidizing simulation.
4.3. 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. Figure 8 plots the acid volume required to break the core for various acid injection capacities under 293.15 K, 338.15 K and 393.15 K. As the acid efficiency curves show, no matter what the isothermal conditions are, the PVBT decreases as 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.15 K 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. Volume of acid required to break the core for various acid injection capacity under isothermal conditions@ 293.15 K, 338.15 K and 393.15 K.
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, although 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 PVBT under any injection capacity, especially under low injection capacity. As we can see from the curves, the maximum difference of PVBT is 46.67 PV which occurs at 0.2 ml/min, and the minimum value is 0.31 PV at 1.5 ml/min.
Figure 10 shows the pressure drop as a function of time for non-isothermal cases with different injection temperatures using the optimum injection capacity obtained from Figure 9. 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  . On the other hand, the pressure drop decreases more slowly with increasing injection temperature under optimum injection
Figure 9. Volumes of acid required to break the core for various acid injection capacity under non-isothermal conditions@ Tinj = 293.15 K, Tinj = 338.15 K and Tinj = 383.15 K.
Figure 10. Pressure drop for non-isothermal conditions using optimum injection capacity.
Figure 11. Volumes of acid breakthrough for non-isothermal conditions using optimum injection capacity.
capacity by comparing the curves in Figure 10. However, it is observed that the trend of PVBT does not agree with that for the pressure drop velocity. The PVBT experiences a decrease with the increasing injection temperature, and it is followed by an increase with the increasing injection temperature as shown in Figure 11, so that an optimum injection temperature also exists.
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 CO2 has been ignored, which will affect the simulation results to some extent. When the supercritical CO2 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 CO2 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 CO2; the fracture-vug-matrix medium and the non-Newtonian acid. In the future studies, some of the extension work will be pursued.
The work is supported by grants from the National Natural Science Foundation of China, No. 51474182 and International Science and Technology Cooperation Research Program of Sichuan Province (China) under Grant, No. 2017HH0014.