Three‐dimensional simulation of the acidizing process under different influencing factors in fractured carbonate reservoirs

Matrix acidizing is widely used to enhance oil/gas production in the exploitation of carbonate reservoirs. In this work, a three‐dimensional (3D) hydro‐chemical‐thermal (H‐C‐T)‐coupled model was presented to improve the understanding of the acidizing process. The influence of different influencing factors was analyzed, especially the coupling effect of natural fractures and in situ stress. With the increase in acid injection concentration, the minimum pore volume of acid required for breakthrough (PVBT) decreases. The optimal injection rate and the minimum PVBT increase with increasing initial reservoir temperature. With the increasing initial reservoir permeability, the minimum PVBT increases. With the increasing initial reservoir pore diameter and specific surface area, the minimum PVBT and the optimal acid injection rate increase. When the fracture direction is perpendicular to the direction of the maximum principal stress, the fracture apertures decrease with the increase of the maximum principal stress, which leads to an increase in PVBT and wider paths of wormholes. Lastly, the present H‐C‐T‐coupled model was applied in the context of Tahe reservoir exploitation, which shows that optimizing the acid injection rate is able to enhance the connection between wellbores and natural caves.


| INTRODUCTION
Carbonate reservoirs contain about 40% of the world's conventional oil/gas reservoirs, such as Tahe reservoir in China. Acidizing treatment is widely used to exploit these carbonate reservoirs. 1,2 Acidizing treatment is achieved by injecting acid into the formation to dissolve the rock at a pressure lower than that of the rock fracture. 3 During the acid-rock reaction, wormholes are created due to the heterogeneity of the rock. 4 These wormholes are highly conductive flow channels that provide favorable paths for fluids to enter the wellbore, thereby enhancing the oil production. 5 There are several theoretical models proposed to simulate the acidizing process in carbonate reservoirs. [6][7][8][9][10][11] The two-scale continuous model proposed by Panga et al. 8 has gradually become the most widely used theoretical model of matrix acidizing because of its strong compatibility. 12,13 It is known that the acidizing process is influenced by many factors. Over the decades, numerous studies have focused on the acidizing process. The acid injection rate has a significant effect on wormhole propagation. [14][15][16][17] Five different dissolution patterns can be observed at different acid injection rates. 18 The acid injection rate that yields dominant dissolution is the optimal acid injection rate because the amount of acid required to break through the core at this rate is minimum. 19 Bazin 17 investigated the effect of acid injection concentration experimentally, which showed that the increase in acid injection concentration leads to a decrease in minimal PV BT . Different acid types correspond to different optimal injection rates, 20 When emulsified acid is used, wormholes form faster in rocks with a high calcite content. 21 Additionally, temperature is a factor that cannot be ignored, [22][23][24][25] as the acid surface reaction rate and diffusion coefficient are sensitive to temperature. 26 With the increase in temperature, there is an increase in the optimum injection rate and corresponding PV BT in limestone acidizing. 17 Ma et al. 27 arrived at the same conclusion as Bazin 17 through simulation. They also concluded that acid temperature is more important than formation temperature in the acidizing process. Furthermore, the heat released from the acid-rock reaction has a minor effect on the acidizing process. 1 According to the authors, the acidizing efficiency of limestone formations is higher at low acid temperatures (below 323 K), but it is higher for dolomite formations at high acid temperatures. Bazin 17 observed that cores with higher permeability require more acid to break through, and the associated optimal injection rates become higher. Using three-dimensional (3D) modelling, Maheshwari et al. 18 concluded that the higher initial permeability of reservoirs leads to an increase in PV BT at the optimal injection rate. By considering the acid-rock reaction in numerical modelling, Jia et al. 28 conducted that the higher initial permeability of reservoirs leads to an increase in PV BT with high injection rates. Jia et al. 12 pointed out that the initial pore diameter of reservoirs has a significant impact on the acidizing process at low injection rates; however, the effect of initial pore diameter can be ignored with high injection rates.
Natural fractures are widespread in carbonate reservoirs, [29][30][31][32] which has a significant effect on the acidizing process. The flow of acid in the matrix is accelerated by the natural fractures. 33 Aidagulov et al. 34 pointed out that most of the injected fluids are consumed by the natural fractures, leading to understimulation in other areas. By comparing reservoirs with and without natural fractures, it was shown that the existence of natural fractures affects PV BT and the dissolution pattern. 35 The influence of natural fracture density and length on the acidizing process was quantitatively analyzed by Zhou et al. 13 Different magnitudes of in situ stress typically exist in reservoirs. The effect of in situ stress on hydraulic fractures has been extensively studied. [36][37][38] Some scholars have tried to take the effect of in situ stress into account in the acidizing treatment, but they have only analyzed the stability of wormholes under stress or the rock's mechanical properties after acidizing. 2,39,40 To date, no numerical studies have systematically analyzed the effect of in situ stress on the acidizing process.
Although several studies have been conducted on the acidizing process, most of them have focused on two-dimensional models. 27 Neither two-dimensional nor linear models can adequately reflect the real reservoir conditions. Furthermore, there is still a lack of systematic understanding of the influences of initial reservoir pore diameter and initial reservoirspecific surface area, especially the coupling effect of natural fractures and in situ stress. In this work, a 3D radial hydro-chemical-thermal (H-C-T)-coupled model was presented to model acidizing process in fractured carbonate reservoirs. In Section 2, the H-C-T coupling model was introduced. In Section 3, the present H-C-T-coupled model was validated by comparing it with pre-existing references. In Section 4, the effect of different influencing factors was analyzed. In Section 5, the effect of coupling effect of natural fractures and in situ stress was analyzed, and the present method was applied under the background of Tahe reservoir exploitation. Finally, the key results are summarized.

| MATHEMATICAL MODEL
In this paper, the present H-C-T coupled model is accomplished by introducing temperature effects into the two-scale continuum model. The two-scale model consists of Darcy-scale model (see Section 2.1) and pore-scale model (see Section 2.2). In this work, the parameters in Darcy-scale model are achieved by solving the pore-scale model. The initial and boundary conditions of the current model are given in Section 2.3. The H-C-T-coupled model is generated based on the following assumptions: (1) the fluid flow obeying Darcy's law occurs in both rock matrix and natural fractures; (2) fluid is incompressible; (3) the influences of gravity and capillary pressure are ignored.

| Fluid flow
The fluid flow in both the fractures and matrix can be described by Darcy's law, and the governing equation is as follows 13 : where φ s and φ j are the porosities of the rock matrix and the fractures, respectively; u s and u j are Darcy's velocity vectors of the rock matrix and the fractures, respectively (m/s); k s and k j are the permeability tensors of the rock matrix and the fracture, respectively (m 2 ); k j = b 2 /12, where b is the fracture aperture, 27,37 and it is an implicit aperture 38,41 ; μ is the acid viscosity of the acid solution (Pa·s); and P is the fluid pressure (Pa). The viscosity of the acid is temperature dependent and can be calculated as 5,42 μ T ( ) = 0.0214 × 10 , where T is the temperature (K).

| Reactive mass transport
The governing equation of the solute transport in both the matrix and fractures is presented as follows: where c f is the acid concentration (mol/m 3 ); D e s and D e j are the effective dispersion tensors of the rock matrix and the fractures, respectively (m 2 /s); and a v is the specific surface area (m −1 ). The acid-rock reaction rate R(c s , T) is defined as a function of temperature T and the acid concentration at the fluid-solid interface c s (mol/(m 3 ·s)). The acid-rock reaction rate can be represented as 26 where k s is the surface dissolution reaction rate (m/s); k c is the local mass transfer coefficient (m/s); E a is the reaction activation energy (J/mol); c s is the concentration of the acid at the fluid-solid interface (mol/m 3 ); and R is the universal gas constant (J/(mol·K)).

| Heat transport
Based on energy balance, the temperature field governing equation is described as follows: where the effective thermal capacity (ρC p ) eff can be represented as Equation (12); ρ s and ρ f are the densities of the rock and acid, respectively (kg/m 3 ); C p,s and C p,f are the thermal capacities of the rock and acid, respectively (J/(kg·K)); the effective thermal conductivity q eff , can be represented as Equation (13); q s T and q f T are the thermal conductivities of the rock and acid, respectively (J/(kg·K)); and δ(i) is the heat generated by the acid-rock reaction (W/m 3 ), as described in Equation (14). ΔH r (T) is the heat created by the reaction when a unit mole of acid is consumed (J/mol), and can be described as follows 26 :

| Rock porosity evolution
For the acid-rock reaction, the porosity evolution process is described as follows 18 : where α is the dissolving power of the acid, defined as the mass of solid dissolved per unit mole of acid reacted (kg/mol).

| Pore-scale model
The parameters (Equations 2, 6, and 16) of a Darcy-scale model depend on the pore-scale phenomenon. 43 To complete the Darcy-scale model, expressions describing the relationships between rock permeability, pore radius, specific surface area, and local porosity are needed. For this, a simple pore-scale model was adopted 44 where k s 0 is the rock initial permeability (m 2 ), β is a dimensionless pore broadening parameter, r p0 is the rock initial pore radius (m), and a v0 is the rock initial specific surface area (m −1 ). The local mass-transfer coefficient and the effective diffusion coefficients are represented as follows 13 : where D ei is the i-directional diffusion coefficients (m 2 /s); i = X denotes the injecting direction, i = T r denotes the transverse direction; D m 0 is the effective molecular diffusivity (m 2 /s) when the temperature is T 0 (K); α os and λ i are numerical constants that are dependent on the pore structure; P ep is the Peclet number, as described in Equation (21); E d is the activation energy of diffusion (J/mol); Sh is the Sherwood number; Sh ∞ is the asymptotic Sherwood number; R ep is the pore Reynolds number, as described in Equation (24); μ v is the kinematic viscosity (m 2 /s), as described in Equation (25); and S c is the Schmidt number, as described by Equation (26).

| Boundary and initial conditions of the H-C-T-coupled model
The computational domain is W. The boundary conditions at the inlet and exit are described as follows: where P e is the fluid pressure at the outlet.
At the upper and lower boundaries, we have the following equations: The initial conditions can be described as follows:

| Computational platform
The numerical models were solved using the finite element method on the COMSOL Multiphysics platform. The simulations were run on a computer workstation with an Intel (R) Xeon (R) Gold 6125 CPU @ 2.10 GHz and with 1 TB of RAM. The fractional step method was used to sequentially solve the fluid pressure, acid concentration, temperature field, and rock porosity. When the total pressure difference in the model was reduced to 1% of its starting value, the computation was complete. 27 The acid required for a breakthrough was normalized with the initial pore volume of the sample, which was called the pore volume of acid required for a breakthrough (PV BT ). 13 All the dissolution patterns are the final porosity field with porosity greater than 0.75.

| VALIDATION OF THE PRESENT H-C-T COUPLING MODEL
A square specimen with a dimension of 2 cm × 2 cm was used, as described by Kalia and Glasbergen. 45 The acid, with an initial temperature of 298 K, was injected into the computational model with an initial temperature of 419 K. The parameters of the computational model are listed in Table 1, which were taken from Kalia and Glasbergen. 45 Figure 1 illustrates the comparison between the breakthrough curves computed using the present H-C-T model and pre-existing reference. The pore volume of acid required for breakthrough (PV BT ) obtained using the present model is in good agreement with the results reported by Kalia and Glasbergen. 45

| Acid injection concentration
As shown in Figure 2, the computational model of fractured reservoirs with one wellbore is established, in which the size of tetrahedral elements changes between 0.02 and 1.8 cm. The dimensions of this cuboid-shaped computational domain are 1 m × 1 m × 0.15 m, and the wellbore radius is r w = 0.02 m. There is one set of natural fractures distributed randomly, which are perpendicular to the xy-plane. The angle of natural fractures is defined as the angle between xaxis and fractures. The length of these natural fractures is 10 ± 2.5 cm, the angle of natural fractures is 180 ± 5°, and the total fracture number is 75. The physical parameters of rock and acid are listed in Table 2. All the values remained the same throughout the study, unless otherwise stated. Three groups of computational models are generated with different acid injection concentrations, including 5%, 15%, and 25%. In each group, seven computational models with different acid injection rates are taken into account. The breakthrough curves of a total of 21 computational models and PV BT are shown in Figure 3. Table 3 shows the minimum PV BT and minimum acid injection mass with different acid injection concentrations at the optimal acid injection rate of 0.1 cm/s. The PV BT represents the amount of acid solution required to break through the formation. Higher acid injection concentration leads to an increase in the minimum PV BT . As the dissolution capacity of the acid increases with increasing acid injection concentration, it leads to a decrease in the minimum PV BT .
To further explore the effect of acid injection concentration on the acidizing process, the mass of solute needed to break through the formation (acid injection mass) was calculated. Figure 4 shows the breakthrough curves for the acid injection mass at different concentrations. It can be seen that acid injection mass increases with the increase in acid injection concentration. The final dissolution patterns with different acid injection concentrations of 5%, 15%, and 25% are shown in Figure 5. Acid injection concentration has no significant effect on the dissolution pattern.

| Initial reservoir temperature
The mesh and boundary conditions of the present models are the same as that in Section 4.1. Due to the heat exchange between acid and wellbores, the acid temperature inside the wellbore is commonly approximately equal to the temperature of nearby reservoirs. 48 Therefore, in this work, the initial temperature of acid is assumed to be the same as the temperature of reservoirs. Three groups of computational models are generated with different initial reservoir temperatures, including 293, 323, and 353 K. In each group, six computational models with different acid injection rates are taken into account. The breakthrough curves of a total of 18 computational models and PV BT are shown in Figure 6. Table 4 shows the minimum PV BT and the optimal acid injection rate with different initial reservoir temperatures. The optimal injection rate and the minimum PV BT increase with increasing initial reservoir temperature. Due to the increase in reservoir temperature, the viscosity, μ, decreases (see Equation 5); meanwhile the surface dissolution reaction rate, k s , (see Equation 9) and the diffusion coefficient, D m , (see Equation 22) both increase, which leads to the fast-reacting acid with low viscosity. Therefore, higher optimal injection rates are required to develop dominant dissolution patterns. Meanwhile, an increase in the optimal injection rate implies that more acid is required to break through the formation, leading to an increase in the minimum PV BT . Additionally, at high injection rates, a larger PV BT is caused by a lower initial acid temperature, and vice versa at high injection rates. These findings support pre-existing experiments reported by Bazin. 17 The evolution of the dissolution pattern and temperature field are shown in Figure 7, with the | 3107 initial reservoir temperature of 323 K and the optimal injection rate of 0.25 cm/s. Initially, the acid first dissolves in the vicinity of the wellbore which leads to an increase in porosity (see Figure 7I(A)) and temperature (see Figure 7II(A)) in the vicinity of the wellbores. Then, as shown in Figure 7I(B), the acid mostly flows into natural fractures and reacts with fracture walls, as the permeability of natural fractures is higher. Simultaneously, the acid-rock reaction and T A B L E 2 Parameters of reservoirs and acid.
T A B L E 3 Minimum PV BT and minimum acid injection mass with different acid injection concentrations at the optimal injection rate of 0.1 cm/s. heat release are reduced in the vicinity of the wellbores. As shown in Figure 7II(B), the hightemperature zone appears along the front of the fracture and wormholes. Finally, as shown in Figure 7I(C), the acid breaks through the formation. The numerical results show that: (1) natural fractures have a significant impact on the dissolution pattern; (2) the high-temperature zones induced by the acidrock reaction are mainly distributed along the front and walls of wormholes (see Figure 7II(C)).

| Initial reservoir permeability
The mesh and boundary conditions of the present computational models are the same as that in Section 4.1. Three groups of computational models are generated with different initial permeability, including 1 × 10 −15 , 1 × 10 −14 , and 1 × 10 −13 m 2 . In each group, seven computational models with different acid injection rates are taken into account. The breakthrough curves of a total of 21 computational models and PV BT are shown in Figure 8. Table 5 shows the influence of initial permeability on the minimum PV BT , which shows that higher initial reservoir permeability leads to an increase in the minimum PV BT . As high initial permeability reduces the resistance of fluid flow, more acid is able to penetrate into the rock matrix. As a result, the minimum PV BT increases with increasing initial reservoir F I G U R E 4 Breakthrough curves of acid injection mass at different acid injection concentrations.
F I G U R E 5 Dissolution patterns with different acid injection concentrations at the optimal injection rate of 0.1 cm/s.

F I G U R E 6
Breakthrough curves with different initial reservoir temperatures.
T A B L E 4 Minimum PV BT and the optimal acid injection rate with different initial reservoir temperatures.
Initial reservoir temperature (K) 293 323 353 Minimum PV BT 11.14 11.49 11.61 Optimal acid injection rate (cm/s) 0.1 0.25 0.5 permeability, which is consistent with the previous experiment. 17 The final dissolution patterns with different initial reservoir permeability are shown in Figure 9. The wormholes become wider as the initial reservoir permeability increases. In the higher initial permeability formation, more acid is consumed when acid breakthroughs the formation, and wider wormholes are generated.
Dissolution patterns at different injection rates are shown in Figure 10. At lower injection rates, most of the acid is consumed before it penetrates deep into the formation, resulting in a face or conical dissolution pattern (see Figure 10A). In contrast, at higher injection rates, the acid can penetrate deeper into the formation, but the insufficient residence time leads to incomplete reactions, producing a ramified or uniform dissolution pattern (see Figure 10C). At intermediate flow rates, that is, the optimal injection rate, a dominant dissolution pattern develops (see Figure 10B). In this dissolution mode, the wormholes are finer and the amount of acid required to break through the formation is minimal. The optimal injection rate plays a key role in achieving efficient and economical acidizing treatments.

| Initial reservoir pore diameter
The mass transfer mechanism from the acid to the surface of the porous medium is mostly influenced by the pore diameter. 12 There are few studies on the influence of initial pore diameter on the acidizing process because the initial pore diameter is complex to determine. Maheshwari et al. 49 experimentally measured that the pore diameter of carbonate rocks is mostly distributed between 1 × 10 −7 and 1 × 10 −5 m. The mesh and boundary conditions of the present computational models are the same as that in Section 4.1. In this section, three groups of computational models are generated with different initial pore diameters, including 1 × 10 −7 , 1 × 10 −6 , and 1 × 10 −5 m. In each group, seven computational models with different acid injection rates are taken into account. The breakthrough curves of a total of 21 computational models and PV BT are shown in Figure 11.  Table 6 shows the influence of initial pore diameter on the minimum PV BT at the same injection rate (0.1 cm/s). The minimum PV BT and the optimal acid injection rate increase with increasing initial reservoir pore diameter. These findings support pre-existing experiments reported by Yoo et al. 50 They believe that there is a high correlation between the optimal acid injection rate and the distribution of pore diameters because the optimal acid injection rate is derived from diffusion coefficients, which are strongly correlated with the distribution of pore diameters.
The final dissolution patterns with different initial reservoir pore diameters at the injection rate of 0.1 cm/s are shown in Figure 12, which shows that a higher initial pore diameter leads to slightly wider wormholes.

| Initial reservoir-specific surface area
The specific surface area is defined as the area available for reaction per unit volume of rocks, and its value can be measured via the gas adsorption (nitrogen) method and petrographic image analysis. 51-53 Jia et al. 12 suggest that the specific surface area mostly varies between 5 × 10 3 and 5 × 10 5 m −1 . The mesh and boundary conditions of the present computational models are the same as that in Section 4.1. Three groups of computational models are generated with different initial specific surface areas, including 5 × 10 3 , 1 × 10 4 , and 2 × 10 4 m −1 . Meanwhile, in each group, seven computational models with different acid injection rates are taken into account. The breakthrough curves of a total of 21 computational models and PV BT are shown in Figure 13. Table 7 shows the influence of the initial specific surface area on the minimum PV BT and optimal injection rate. Higher initial specific surface area leads to an increase in the optimal injection rate and the minimum PV BT . When the initial specific surface area is larger, the acid consumed per unit F I G U R E 8 Breakthrough curves with different initial reservoir permeabilities.
T A B L E 5 Minimum PV BT with different initial reservoir permeability at the optimal injection rate of 0.1 cm/s. volume of the rock matrix increases. Therefore, higher injection rates are required to develop dominant dissolution patterns. Meanwhile, an increase in the optimal injection rate implies that more acid is required to break through the formation, leading to an increase in the minimum PV BT . The final dissolution patterns with different initial specific surface areas at the same injection rate of 0.1 cm/s are shown in Figure 14. When the initial specific surface area increases, the reaction between the acid and the rock is enhanced, resulting in increased acid consumption, and more rock is dissolved. As a result, wider wormholes are formed.

| THE COMBINED AIFRAC AND THE H-C-T CALCULATION METHOD
In this section, The AiFrac is introduced to model the deformation and stress-strain of reservoirs based on the hybrid of finite element and mesh-free method (FEMM). 41,[54][55] The computation process is divided into two stages: first, the reservoir deformation and fracture aperture evolution are solved using AiFrac; then, the acidizing process is solved using the present H-C-T coupling model.

| The coupling effect of natural fractures and in situ stress
In this section, the coupling effect of natural fractures and in situ stress is investigated by combining the AiFrac and the present H-C-T coupling model. As shown in Figure 15, the computational model of fractured reservoirs with one wellbore is established. The dimensions of this cuboid-shaped computational domain are 1 m × 1 m × 0.15 m, and the wellbore radius is r w = 0.02 m. The size of tetrahedral elements varies between 0.02 and 1.8 cm. There are a total of 44 natural fractures, which are parallel to the x-direction and are perpendicular to the xy-plane; meanwhile, the length of these natural  Figure 16. Second, the acidizing process is solved using the present H-C-T coupling model. The dissolution patterns for different in situ stress at the same injection rate of 0.1 cm/s are shown in Figure 17. Because σ x and σ y are parallel and perpendicular to the direction of the natural fractures, respectively, the natural fracture aperture decreases with the increase in situ stress deviation Δσ. The PV BT for Case I, Case II, and Case III are 7.71, 8.54, and 16.38, respectively. When the fracture direction is perpendicular to the direction of the maximum principal stress, an increase in the maximum principal stress leads to a decrease in the fracture aperture and an increase in the PV BT . This is because the decrease in fracture aperture decreases the fracture permeability, which increases the resistance of acid flow in the fracture. Consequently, the time for the acid to break through the formation is increased, resulting in an increase in the PV BT and wider wormholes are formed.

| Application in Tahe reservoir exploitation
In fractured-cavity carbonate reservoirs, oil/gas resources are mainly distributed in natural cavities. 56 Therefore, in Tahe reservoir, as a typical fractured-cavity carbonate F I G U R E 12 Dissolution patterns with different initial pore diameters at the optimal injection rate of 0.1 cm/s. (A) Dissolution pattern when the initial pore diameter is 1 × 10 -7 m. (B) Dissolution pattern when the initial pore diameter is 1 × 10 -6 m. (C) Dissolution pattern when the initial pore diameter is 1 × 10 -5 m.
F I G U R E 13 Breakthrough curves with different initial specific surface areas.
T A B L E 7 Minimum PV BT and the optimal acid injection rate with different initial reservoir surface areas.
Initial reservoir-specific surface area (m −1 ) 5 × 10 3 1 × 10 4 2 × 10 4 Minimum PV BT 11.14 11.25 11.56 The optimal acid injection rate (cm/s) reservoir, it is critical to generate high-permeability channels between wellbores and natural cavities. The natural caves might be filled with sediments (see Figure 18A). 56 Seimic surveys are widely used to determine the location and shape of natural caves (see Figure 18B). 55,58 Then, the treatment parameters are optimized to control the geometry of dissolution patterns, to connect wellbores and natural caves. As shown in Figure 19, the computational model of fractured reservoirs with one wellbore is established. The dimensions of this cuboid-shaped computational domain are 1 m × 1 m × 0.15 m and the wellbore radius is r w = 0.02 m. The wellbore center is located at the coordinate of (0, 0.77, 0.075). The size of tetrahedral elements varies between 0.02 and 1.8 cm. There are a total of 60 natural fractures. These fractures are parallel to the x-direction and are perpendicular to xy-plane; meanwhile, the length of these natural fractures is 12 cm. The natural cave shape is set according to the field observation ( Figure 18(B)). The caves are naturally pre-existing, whose porosity is initially 0.9. 13 The in situ stress is set as σ x = 60 MPa and σ y = 30 MPa. Young's module of the natural cavity is 5 GPa. 56 The mechanical parameters of the rock matrix are the same as those in Section 5.1.
The final dissolution patterns with different injection rates of 0.05 and 0.1 cm/s are shown in Figure 20. When the acid injection rate is lower (0.05 cm/s), the achieved wormholes are able to connect with the cavity (see Figure 20A); however, at a higher injection rate (0.1 cm/s), it fails to connect with the cavity (see Figure 20B). At the optimal injection rate, the wormholes formed are thin and long, and the area dissolved in the formation is not large. When the injection rate is less than the optimal injection rate, the diffusion of acid in porous media is more pronounced, and the acid can react completely with the rock in the flow zone. This leads to an increase in the dissolved area in the formation, resulting in wider wormholes, which means that a lower injection rate is helpful to connect more natural caves in the vicinity of wellbores (see Figure 20B).

| CONCLUSIONS
The H-C-T-coupled model was presented to investigate the acidizing process. The effects of different influencing factors were analyzed numerically. The application of the H-C-T-coupled model in Tahe reservoir exploitation is demonstrated. The main conclusions are summarized as follows: 1. The higher acid injection concentration led to a decrease in the minimal PV BT . The higher initial reservoir temperature, pore diameter, and specific surface area leads to an increase in the optimal acid injection rate and the minimum PV BT . The higher initial reservoir permeability leads to an increase in the minimal PV BT . At the optimal injection rate, the wormholes are slightly wider with high reservoir permeability. 2. When the fracture direction is perpendicular to the direction of the maximum principal stress, an increase in the maximum principal stress leads to an increase in the PV BT and wider wormholes. 3. By modelling the acidizing process under the application background of Tahe reservoir exploitation, it is shown that a lower injection rate is helpful to connect more natural caves in the vicinity of wellbores.