Regular Article
Threedimensional simulation of acidizing process in carbonate rocks using the Darcy–Forchheimer framework
^{1}
Qingdao Key Laboratory for Geomechanics and Offshore Underground Engineering, School of Science, Qingdao University of Technology, 266033 Qingdao, PR China
^{2}
School of Petroleum Engineering, China University of Petroleum (Huadong), 266555 Qingdao, PR China
^{*} Corresponding author: peiyang988@gmail.com
Received:
10
January
2020
Accepted:
12
May
2020
Acidizing is an economical and effective practice to remove the near wellbore damage, which is performed by injecting acid into the formation through the wellbore. The injected acid dissolves the rock, by which the permeability nearby the wellbore can be improved. For a carbonate reservoir, the injected acid dissolves some of the minerals and some narrow and long channels, named wormholes, are formed then. These wormholes can bypass the damaged zone and hence improve the productivity of the well. The process for acid dissolving rocks involves complex physicochemical change, including the chemical reactions at the pore scale and the fluid flow at Darcy scale. In this paper, a 3D reactive flow model with nonDarcy framework is developed based on the twoscale continuum model, and is solved by using the finite volume method. Five types of dissolution patterns, named face dissolution, conical wormhole, wormhole, ramified wormhole, and uniform dissolution, are obtained as the injection velocity increases. The effect of nonDarcy flow on dissolution pattern and breakthrough volume is analyzed. It is found that there is no effect of nonDarcy on dissolution structure and breakthrough volume when the injection velocity is very low. However, when the injection velocity is very high, the generated wormhole has more branches when using the Forchheimer equation than using the Darcy equation. Moreover, the optimal injection velocity is found to be the same whether considering the nonDarcy flow or not.
© P. Liu et al., published by IFP Energies nouvelles, 2020
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
1 Introduction
Acidizing is a commonly performed stimulation treatment for carbonate reservoirs. Acid is injected into the formation through the wellbore under the fracture pressure. The injected acid reacts with the rock and, if successful, creates some highly permeable channels. The channels are known as wormholes, which can bypass the damaged zone and hence yield a productivity enhancement [1]. However, the wormhole has not always been formed after injecting acid into the formation. Experimental observations demonstrate that different types of dissolution patterns are formed when injecting acid into carbonate rocks. The dissolution pattern depends on the injection rate, type and concentration of the acid, type of the minerals within the rock, heterogeneity of the rock, as well as the temperature, and so on [2–5]. For example, when the injection rate is quite low, the injected acid is completely consumed before it penetrates deeply into the rock. The rock near the inlet end is dissolved completely, and as a result, a face dissolution pattern is formed. On the contrary, when the acid is injected into the rock at a very high rate, live acid will invade almost all of the pores of the rock. The porosity is increased uniformly which leads to a uniform dissolution pattern. As the injection rate increases, between these two extremes, conical dissolution, wormhole, and ramified dissolution patterns are formed subsequently [6]. By comparing the volume of acid required to break through the core, it is found that the breakthrough volume is a minimum when the wormhole dissolution pattern is formed [7]. Consequently, there is a practical motivation for finding an optimal injection condition, at which the wormhole dissolution is formed.
To investigate the dissolution process, numerous experiments have been carried out, which provide a direct observation of the dissolution dynamic in rocks as well as a fundament for mathematical modelling to predict the wormhole formation. The numerical models that have been proposed over the last few decades, can be broadly classified into four types: (1) dimensionless model [8–10]; (2) capillary tube model [11–15]; (3) network model [7, 16–19]; and (4) continuum model [6, 20–23]. Reviews of these models, along with summaries of their limitations and assumptions, can be found in [24–28]. Because the continuum model is better at forecasting the dissolution patterns observed in experiments, it can estimate the breakthrough volumes accurately. Till now, this model has been extensively used over the last few years. The continuum model was first presented by Liu et al. [29], which, unfortunately, not valid in the mass transfer controlled regime. Then Golfier et al. [25] developed a model considering the effect of mass transfer but it is not valid in the kinetic regime. Based on previous works, Panga et al. [26] presented a continuum model that can capture both the kinetic and mass transfer controlled regime simultaneously. After that, the continuum model has been widely used to investigate the effect of rock or fluid properties, including flow geometry (linear and radial) [30–32], completion methods [33], the type of injected acid (hydrochloric acid and selfdiverting acid) [34–38], medium heterogeneity (the presence of vugs and fractures) [20, 23, 39, 40], and reservoir temperature [41–43] on wormhole propagation and optimal injection volume.
Although the above studies bring important insights to the dissolution process, the equations that refer to the fluid flow in the rock still remain further improved. In most of the published studies, Darcy’s law is used to relate the flow velocity in the pressured drop, rock permeability, and fluid viscosity. Nevertheless, Darcy’s law is valid only for laminar flows. During the dissolution process, the porosity in the dissolved area will increase as the acid continues being injected. Highly permeable channels may also be generated in some areas, where the rock is completely dissolved. Additionally, the flow in these areas is nonlinear and Darcy’s law is no longer applicable. In order to accurately describe the fluid flow in dissolved area, some investigators [44, 45] use Darcy–Brinkman–Forchheimer framework instead of Darcy’s equation to simulate the reactivetransport process. However, their work is limited to the 2D condition, which cannot accurately represent the actual acidizing operation. The main aim of this paper is to extend the reactivetransport model to 3D nonDarcy flow condition by basing on the Darcy–Forchheimer framework. The simulation results are compared with that calculated using the Darcy framework in order to investigate the effect of nonDarcy flow on dissolution patterns and breakthrough volume.
2 Mathematical model
The twoscale continuum model developed by Panga et al. [26] consists of a Darcyscale model and a porescale model. The Darcyscale model describes the reactivetransport process of acid and the evolution of the rock. By replacing the Darcy framework with the Darcy–Forchheimer framework to describe the fluid flow field, the Darcyscale model can be expressed as follows:(1)(2)(3)(4)(5)where, P is the pressure of the porefluid; μ is the viscosity of the fluid phase; K is the permeability tensor, which is calculated from a porescale model or determined by lab measurement; u is the velocity vector; β is referred as the Forchheimer number that can be obtained experimentally; ρ is the density of the fluid; t is the time; ϕ is the porosity of the rock; C_{f} is the cupmixing concentration of the solute in the fluid phase; C_{s} is the concentration of the acid at the solid–fluid interface; D_{e} is the effective diffusion tensor; a_{v} is the interfacial area available for reaction per unit volume of the medium; k_{c} is the mass transfer coefficient; R(C_{s}) is the dissolution rate; α is the dissolving power of the acid, defined as grams of dissolved solid per mole of acid reacted; and ρ_{s} is the density of the solid.
The porosity increases as the rock is dissolved, which results in the changes in permeability, specific area a_{v} and pore radius r_{p}. Because these variables are pore structure dependent, a porescale model that captures the evolution of these variables with the variation in porosity is needed to complete the Darcyscale model. We adopt here a fractal theory based formulation to relate the permeability with the porosity as:(6)with the definition of the fractal dimension of pore spaces and tortuosity, D_{f} and D_{T}, respectively as,(7)(8)where d_{E} is the Euclidean dimension; λ_{min} and λ_{max} are the smallest pore diameter and the largest diameter, respectively. The typical value of λ_{min}/λ_{max} is 0.01 as analyzed in Liu et al. [6]. In equation (6), the parameter with subscript 0 indicates the value of variable when porosity is equal to ϕ_{0}, and ϕ_{0} represents the initial porosity.
The pore radius and specific area are calculated by:(9)(10)
The mass transfer coefficient, and the effective dispersion coefficients, inside the pores are expressed by the relations:(11)(12)(13)where Sh is the Sherwood number, which represents dimensionless mass transfer coefficient; D_{m} is the molecular diffusivity; Re_{p} is the pore Reynolds number defined as Re_{p} = 2ur_{p}/ν, ν is the kinematic viscosity; Sc is the Schmidt number defined as Sc = ν/D_{m}; D_{eX} and D_{eT} are the longitudinal and transverse dispersion coefficient, respectively. α_{os}, λ_{X}, and λ_{T} are constants that depend upon the pore structure and their typical values are 0.5, 0.5, 0.1 for a packedbed of spheres, respectively.
A 3D cubic domain is considered in our simulation to mimic the real core acidizing experiment. The acid is injected into the carbonate core at a constant injection rate and the pressure at the outlet boundary is fixed. Acid cannot flow out through the transverse boundaries. Therefore, the boundary and initial conditions are given by:(14)(15)(16)(17)(18)where v_{0} is the injection velocity; C_{0} is the concentration of the injected acid; P_{e} is the pressure at the exit boundary of the domain; is a random fluctuation in the initial porosity field that is introduced to represent heterogeneity in the medium, and is assumed to be uniformly distributed in the interval [−Δϕ_{0}, Δϕ_{0}].
3 Dimensionless model
The dimensionless variables and parameters are defined as follow,where the subscript D represents dimensionless variable; L is the length of the core; H and W are the height and width of the domain, respectively; U is the dimensionless velocities vector. η is the poretodomain length ratio; α_{LH} and α_{LW} are the aspect ratios; κ is the dimensionless permeability; B is the dimensionless nonDarcy coefficient. The parameters obtained from the nondimensionalization are the pore scale Thiele modulus , defined as the ratio of diffusion time to reaction time, Damkӧhler number Da, defined as the ratio of convection time to reaction time; the axial Peclet number Pe_{L}, defined as the ratio of advective transport rate to diffusive transport rate; the acid capacity number N_{ac}, defined as the volume of solid dissolved per unit volume of acid consumed; the macroscopic Thiele modulus Φ^{2} is corescale equivalent of the porescale Thiele modulus. More details on the definition of these dimensionless groups can be found in Panga et al. [26]. The PDEs after conversion into the nondimensionalized form are as follows:(19)(20)(21)(22)and the resulting boundary and initial conditions are,(23)(24)
4 Numerical schemes
The mathematical model presented in the above section is solved sequentially. Firstly, the momentum and continuity equations are solved to obtain the velocity and pressure fields. Because a strong coupling exists between velocity and pressure, equations (19) and (20) are solved simultaneously using a staggered grid. The pressure field is stored at cell centroids, while the velocity field is stored at cell faces. The coupling between velocity and pressure is enforced in a staggered grid because no interpolation is needed to calculate the velocity field in the continuity equation. After the velocity field is known, it is substituted into the reactivetransport equation (Eq. (21)) to calculate the acid concentration. The operator splitting method is adopted here for solving the reactivetransport coupling problem as used in Maheshwari et al. [27] and Liu et al. [43]. In a time step, the diffustion–convection operator is solved firstly to calculate the acid concentration at time t*, at which the transport process has completed, but the reaction process has not yet started. Then, the reaction operator is solved by taking the concentration value at time t* as the initial value. Mathematically, we solve the following equations in the first step:(25)and the following equations in the second step:(26)
By solving equations (25) and (26), the distribution of acid concentration and porosity field at the new time step can be obtained. Substituting the porosity value into the porescale model, the permeability field involved in the momentum equation can be updated, and then the algorithm marches forward.
The procedure described above is repeated until the acid breaks through the rock. The breakthrough is defined when the pressure difference between the outlet and the inlet is negligible. Numerically, the breakthrough time is taken as the time when the overall pressure drop across the medium drops to 1/100th of its initial value.
The governing equations are discretized using the finite volume method [32]. In order to avoid the instability issues, the model is solved implicitly in time, i.e. the transient term is discretized using the backward Euler scheme [46–48]. For spatial term, the central difference and upwind schemes are used to discretize the diffusion and convection terms, respectively. Therefore, the numerical scheme is stable and has a first order accuracy. The simulations presented in this work are performed with the number of grid fixed at 100 × 40 × 40 (100 cells in x direction, 40 cells in y and z directions). The grid size is established by comparing the simulation accuracy as described in Liu et al. [43]. To put it briefly, the grid is refined until the simulation result is insensitive to any grid size changes.
5 Model testing
Before performing the numerical simulation, the validity of the proposed model is tested. This is achieved by comparing the numerical results with the experimental observations reported by Seagraves et al. [49]. In their acidizing experiments, 15 wt% HCl was injected into Indiana limestone at a constant rate of 25 mL/min. The limestone cores are 17.8 cm in diameter and 45.7 cm in length. After acidizing, the core samples were CT scanned and 3D reconstruction of the internal void space was created to show the wormhole patterns generated by acid dissolution. We use our model to replicate the experimental conditions of Test 7 in Seagraves et al. [49]. The core used in Test 7 has a permeability of 22.2 mD and a porosity of 0.133. Because some parameters used in our numerical model, such as the interface area and average pore diameter, are not given in the experimental report, we use the values reported in previous similar acidizing numerical simulation studies to estimate the values of these unspecified or unknown parameters. The values of these parameters are given in Table 1. Figure 1 displays the comparison of the simulated wormhole pattern with the experimental result reported in Seagraves et al. [49]. It can be seen that a repeatedly branched wormhole pattern observed in Seagraves’s acidizing experiment is also predicted in our numerical simulation. The slight difference in wormhole structure can be ascribed to the difference in initial porosity field of rock used in the experiment and porous medium generated for numerical simulation.
Fig. 1 Comparison of (a) the wormhole structures simulated using the presented model with and (b) the experimental observation reported by Seagraves et al. [49]. 
List of parameters and dimensionless numbers used in simulation.
6 Results and analysis
In this section, the numerical simulations are performed to analyze the sensitivity of the dissolution dynamics. The values of dimensionless numbers and parameters used in the calculations are listed in Table 1. All the values remain fixed unless otherwise stated. Because the effect of injection rate, acid properties and rock heterogeneous on dissolution patterns and breakthrough volume has been studied in detail in previous literature [6, 27, 32, 44, 50–54], it will not be repeated here. This work only focuses on the influence of nonDarcy flow on the dissolution structure and breakthrough volume.
6.1 Dissolution patterns
Figure 2 shows the dissolution patterns obtained from the simulation by using the Darcy framework and Darcy–Forchheimer framework, respectively. It can be seen that, five types of dissolution patterns, named face dissolution, conical wormhole, dominate wormhole, ramified wormhole, and uniform dissolution patterns are observed for both frameworks as the injection rate varies. The formation of different dissolution patterns is determined by the combined effects of the solute transport rate and the reaction rate, and is independent of the type of flow equation (whether nonDarcy or not). Another conclusion obtained from Figure 2 is that the dissolution structures are almost same for using Darcy framework and Darcy–Forchheimer framework when the injection rate is lower (Figs. 2a–2c) by comparing each volume. However, when the injection rate is higher, which corresponds to the formation of the ramified wormhole and uniform dissolution patterns, the wormhole structure calculated by using the Darcy–Forchheimer framework is more branching than that using Darcy framework. This is caused by the effect of inertial forces and the formation of highspeed nonDarcy resistance term in the Darcy–Forchheimer framework. Wherein, this nonDarcy resistance term is equivalent to reducing the mobility of the high permeable area to make the injected acid flow into other low permeable areas.
Fig. 2 Dissolution patterns obtained from simulations by using of Darcy and Forchheimer framework. 
6.2 Breakthrough curve
Figure 3 shows the breakthrough curves for Darcy framework or Darcy–Forchheimer framework. Referring to the previous works [27, 39], the breakthrough curves are expressed by plotting the pore volumes to breakthrough (PV_{BT}) against the injection velocity (v). The PV_{BT} is defined as the acid volume required for breakthrough, per unit initial pore volume of the media. Mathematically, it can be expressed as:(27)where A is the area of the injection end; V is the bulk volume of the rock; t_{BT} is the breakthrough time.
Fig. 3 Breakthrough curves obtained from simulations with of Darcy and Forchheimer framework. 
It can be seen from Figure 3 that the values of the breakthrough volume calculated using the Darcy equation and the Forchheimer equation are the same when the injection rate is low. However, as the injection rate increases, the breakthrough volume calculated with Darcy–Forchheimer equation is larger than that with the Darcy equation, and this difference becomes larger as the injection rate increases. This trend in breakthrough curves is consistent with the changing trend of the dissolution structure, because the more branches of the wormhole are dissolved, the more acid will be needed. The curves in Figure 3 show a minimum at intermediate Da, which is corresponding to the wormhole formation and called the optimum injection rate. The optimum injection rate is practically important, since it is usually extrapolated to define the optimum flow rate to be used in the field in order to have a successful treatment at low cost. It can be seen from Figure 3 that the optimum injection rate for the using of Darcy framework and Darcy–Forchheimer framework is the same. The reason is that the optimal injection velocity is only related to the property of the rock and acid itself, irrespective of the equation that describes the fluid flow in the rock.
7 Conclusion
In this paper, the effect of nonDarcy flow on dissolution structure and breakthrough volume is studied by numerically solving the 3D twoscale continuum model with the Darcy–Forchheimer framework. The main results can be summarized as following:
When the injection rate is high, the wormhole structure calculated by using the Darcy–Forchheimer framework is more branching than that using Darcy framework. When the injection rate is low, the dissolution structure is insensitive to the type of flow equation.
When the injection rate is high, the breakthrough volume calculated using the Forchheimer equation is larger than that calculated by using the Darcy equation, and the difference increases as the injection rate increases.
There is no effect on optimal injection rate whether or not to consider nonDarcy flow influence. The optimal injection rate calculated using the Forchheimer equation and the Darcy equation is the same.
Acknowledgments
The authors gratefully acknowledge the support from the National Natural Science Foundation of China (No. 51804325), the China Postdoctoral Science Foundation (No. 2019M652508), and the National Natural Science Foundation of China (No. 51874262).
References
 Economides M.J., Nolte K.G. (2000) Reservoir stimulation, Wiley, Chichester, pp. 131–1715. [Google Scholar]
 Daccord G. (1987) Chemical dissolution of a porous medium by a reactive fluid, Phys. Rev. Lett. 58, 5, 479. [CrossRef] [PubMed] [Google Scholar]
 Bazin B., Abdulahad G. (1999) Experimental investigation of some properties of emulsified acid systems for stimulation of carbonate formations, in: Proceedings of the SPE Middle East Oil Show & Conference, pp. 347–356. [Google Scholar]
 Bazin B. (2001) From matrix acidizing to acid fracturing: A laboratory evaluation of acid/rock interactions, SPE Prod. Facil. 16, 1, 22–29. [CrossRef] [Google Scholar]
 Fredd C.N., Fogler H.S. (1999) Optimum conditions for wormhole formation in carbonate porous media: Influence of transport and reaction, SPE J. 4, 3, 196–205. [CrossRef] [Google Scholar]
 Liu P., Yao J., Couples G.D., Ma J., Iliev O. (2017) 3D modelling and experimental comparison of reactive flow in carbonates under radial flow conditions, Sci. Rep. 7, 1, 17711. [CrossRef] [PubMed] [Google Scholar]
 Fredd C.N., Fogler H.S. (1998) Influence of transport and reaction on wormhole formation in porous media, AlChE J. 44, 9, 1933–1949. [CrossRef] [Google Scholar]
 Daccord G., Touboul E., Lenormand R. (1989) Carbonate acidizing: toward a quantitative model of the wormholing phenomenon, SPE Prod. Eng. 4, 1, 63–68. [CrossRef] [Google Scholar]
 Daccord G., Lenormand R., Liétard O. (1993) Chemical dissolution of a porous medium by a reactive fluid – I. Model for the “wormholing” phenomenon, Chem. Eng. Sci. 48, 1, 169–178. [Google Scholar]
 Daccord G., Lietard O., Lenormand R. (1993) Chemical dissolution of a porous medium by a reactive fluid – II. Convection vs. reaction, behavior diagram, Chem. Eng. Sci. 48, 1, 179–186. [Google Scholar]
 Schechter R.S., Gidley J.L. (1969) The change in pore size distribution from surface reactions in porous media, AlChE J. 15, 3, 339–350. [CrossRef] [Google Scholar]
 Hung K.M., Hill A.D., Sepehrnoori K. (1989) A mechanistic model of wormhole growth in carbonate matrix acidizing and acid fracturing, J. Petrol. Technol. 41, 1, 59–66. [CrossRef] [Google Scholar]
 Gdanski R. (1999) A fundamentally new model of acid wormholing in carbonates, in: Proceedings of the SPE European Formation Damage Conference, Society of Petroleum Engineers, The Hague, Netherlands, pp. 1–10. [Google Scholar]
 Gong M., ElRabaa A.M. (1999) Quantitative model of wormholing process in carbonate acidizing, in: Proceedings of Society of Petroleum Engineers, pp. 1–11. [Google Scholar]
 Buijse M.A. (2000) Understanding wormholing mechanisms can improve acid treatments in carbonate formations, SPE Prod. Facil. 15, 3, 168–175. [CrossRef] [Google Scholar]
 Hoefner M.L., Fogler H.S. (1988) Pore evolution and channel formation during flow and reaction in porous media, AlChE J. 34, 1, 45–54. [CrossRef] [Google Scholar]
 Kang Q., Lichtner P.C., Viswanathan H.S., AbdelFattah A.I. (2010) Pore scale modeling of reactive transport involved in geologic CO_{2} sequestration, Transp. Porous Media 82, 1, 197–213. [Google Scholar]
 Budek A., Szymczak P. (2012) Network models of dissolution of porous media, Phys. Rev. E 86, 5, 056318. [Google Scholar]
 Tansey J. (2014) Porenetwork modeling of carbonate acidization, in: Proceedings of the SPE Annual Technical Conference and Exhibition, Society of Petroleum Engineers, pp. 1–10. [Google Scholar]
 Liu P., Yao J., Couples G.D., Ma J., Huang Z., Sun H. (2017) Modelling and simulation of wormhole formation during acidization of fractured carbonate rocks, J. Petrol. Sci. Eng. 154, 284–301. [CrossRef] [Google Scholar]
 Yao J., Liu P., Huang Z., Wang Y., Yan X. (2017) Analysis of influencing factors on the optimum stimulation conditions of the acidizing treatment in carbonate reservoirs, Sci. Sin. Tech. 47, 1, 1–16. [CrossRef] [Google Scholar]
 Yao J., Liu P., Huang Z., Wang Y., Yan X., Zeng Q. (2017) Status and progress of reactive flow simulations for carbonate reservoirs, Earth Sci. 42, 8, 1263–1272. [Google Scholar]
 Liu P., Yao J., Couples G.D., Huang Z., Hai S., Ma J. (2017) Numerical modelling and analysis of reactive flow and wormhole formation in fractured carbonate rocks, Chem. Eng. Sci. 172, 143–157. [Google Scholar]
 Fredd C.N., Miller M.J. (2000) Validation of carbonate matrix stimulation models, in: Proceedings of the SPE International Symposium on Formation Damage Control, Society of Petroleum Engineers, pp. 1–14 [Google Scholar]
 Golfier F., Cesar Z., Bazin B., Lenormand R., Lasseux D., Quintard M. (2002) On the ability of a Darcyscale model to capture wormhole formation during the dissolution of a porous medium, J. Fluid. Mech. 457, 213–254. [Google Scholar]
 Panga M.K.R., Ziauddin M., Balakotaiah V. (2005) Twoscale continuum model for simulation of wormholes in carbonate acidization, AlChE J. 51, 12, 3231–3248. [CrossRef] [Google Scholar]
 Maheshwari P., Ratnakar R.R., Kalia N., Balakotaiah V. (2013) 3D simulation and analysis of reactive dissolution and wormhole formation in carbonate rocks, Chem. Eng. Sci. 90, 258–274. [Google Scholar]
 Ghommem M., Zhao W., Dyer S., Qiu X., Brady D. (2015) Carbonate acidizing: Modeling, analysis, and characterization of wormhole formation and propagation, J. Petrol. Sci. Eng. 131, 18–33. [CrossRef] [Google Scholar]
 Liu X., Ormond A., Bartko K., Li Y., Ortoleva P. (1997) A geochemical reactiontransport simulator for matrix acidizing analysis and design, J. Petrol. Sci. Eng. 17, 1, 181–196. [CrossRef] [Google Scholar]
 Kalia N., Balakotaiah V. (2007) Modeling and analysis of wormhole formation in reactive dissolution of carbonate rocks, Chem. Eng. Sci. 62, 4, 919–928. [Google Scholar]
 Cohen C.E., Ding D., Quintard M., Bazin B. (2008) From pore scale to wellbore scale: Impact of geometry on wormhole growth in carbonate acidization, Chem. Eng. Sci. 63, 12, 3088–3099. [Google Scholar]
 Liu P., Couples G.D., Yao J., Huang Z., Song W., Ma J. (2018) A general method for simulating reactive dissolution in carbonate rocks with arbitrary geometry, Comput. Geosci. 22, 5, 1187–1201. [Google Scholar]
 Kalia N., Balakotaiah V. (2010) Wormholing in perforated completions, in: Proceedings of the SPE International Symposium and Exhibiton on Formation Damage Control, Society of Petroleum Engineers, Lafayette, pp. 1–17. [Google Scholar]
 Ratnakar R., Kalia N., Balakotaiah V. (2012) Carbonate matrix acidizing with gelled acids: An experimentbased modeling study, in: Proceedings of the SPE International Production and Operations Conference & Exhibition, Society of Petroleum Engineers, pp. 1–16. [Google Scholar]
 Maheshwari P., Balakotaiah V. (2013) 3D simulation of carbonate acidization with HCl: Comparison with experiments, in: Proceedings of the SPE Production and Operations Symposium, Society of Petroleum Engineers, Oklahoma, pp. 1–17. [Google Scholar]
 Ratnakar R.R., Kalia N., Balakotaiah V. (2013) Modeling, analysis and simulation of wormhole formation in carbonate rocks with in situ crosslinked acids, Chem. Eng. Sci. 90, 179–199. [Google Scholar]
 Maheshwari P., Maxey J., Balakotaiah V. (2015) Reactivedissolution modeling and experimental comparison of wormhole formation in carbonates with gelled and emulsified acids, SPE Prod. Oper. 31, 2, 103–119. [Google Scholar]
 Maheshwari P., Maxey J.E., Balakotaiah V. (2014) Simulation and analysis of carbonate acidization with gelled and emulsified acids, in: Proceedings of the Abu Dhabi International Petroleum Exhibition and Conference, Society of Petroleum Engineers, Abu Dhabi, pp. 1–28. [Google Scholar]
 Kalia N., Balakotaiah V. (2009) Effect of medium heterogeneities on reactive dissolution of carbonates, Chem. Eng. Sci. 64, 2, 376–390. [Google Scholar]
 Izgec O., Zhu D., Hill A.D. (2010) Numerical and experimental investigation of acid wormholing during acidization of vuggy carbonate rocks, J. Petrol. Sci. Eng. 74, 1, 51–66. [CrossRef] [Google Scholar]
 Kalia N., Glasbergen G. (2009) Wormhole formation in carbonates under varying temperature conditions, in: Proceedings of the European Formation Damage Conference, Society of Petroleum Engineers, Scheveningen, pp. 1–19. [Google Scholar]
 Kalia N., Glasbergen G. (2010) Fluid temperature as a design parameter in carbonate matrix acidizing, in: Proceedings of the SPE Production and Operations Conference and Exhibition, Society of Petroleum Engineers, pp. 1–21. [Google Scholar]
 Liu P., Yan X., Yao J., Sun S. (2019) Modeling and analysis of the acidizing process in carbonate rocks using a twophase thermalhydrologicchemical coupled model, Chem. Eng. Sci. 207, 215–234. [Google Scholar]
 Wu Y., Salama A., Sun S. (2015) Parallel simulation of wormhole propagation with the Darcy–Brinkman–Forchheimer framework, Comput. Geotech. 69, 564–577. [Google Scholar]
 Yuan T., Yang N., Qin G. (2016) Numerical modeling and simulation of coupled processes of mineral dissolution and fluid flow in fractured carbonate formations, Transp. Porous Media 114, 3, 747–775. [Google Scholar]
 Zeng Q., Yao J., Shao J. (2020) An extended finite element solution for hydraulic fracturing with thermohydroelastic–plastic coupling, Comput. Methods Appl. Mech. Eng. 364, 112967. [Google Scholar]
 Yan X., Huang Z., Yao J., Li Y., Fan D., Sun H., Zhang K. (2018) An efficient numerical hybrid model for multiphase flow in deformable fracturedshale reservoirs, SPE J. 23, 4, 1412–1437. [CrossRef] [Google Scholar]
 Li J., Zhang T., Sun S., Yu B. (2019) Numerical investigation of the POD reducedorder model for fast predictions of twophase flows in porous media, Int. J. Numer. Meth. Heat Fluid Flow 29, 11, 4167–4204. [Google Scholar]
 Seagraves A.N., Smart M.E., Ziauddin M.E. (2018) Fundamental wormhole characteristics in acid stimulation of perforated carbonates, in: Proceedings of the SPE International Conference and Exhibition on Formation Damage Control, Society of Petroleum Engineers, Lafayette, Louisiana, USA, 50 p. [Google Scholar]
 Li Y., Liao Y., Zhao J., Peng Y., Pu X. (2017) Simulation and analysis of wormhole formation in carbonate rocks considering heat transmission process, J. Nat. Gas. Sci. Eng. 42, 120–132. [Google Scholar]
 Liu P., Xue H., Zhao L., Zhao X., Cui M. (2016) Simulation of 3D multiscale wormhole propagation in carbonates considering correlation spatial distribution of petrophysical properties, J. Nat. Gas Sci. Eng. 32, 81–94. [Google Scholar]
 Liu N., Liu M. (2016) Simulation and analysis of wormhole propagation by VES acid in carbonate acidizing, J. Petrol. Sci. Eng. 138, 57–65. [CrossRef] [Google Scholar]
 Maheshwari P., Balakotaiah V. (2013) Comparison of carbonate HCl acidizing experiments with 3D simulations, SPE J. 28, 4, 402–413. [Google Scholar]
 Liu M., Zhang S., Mou J., Zhou F. (2013) Wormhole propagation behavior under reservoir condition in carbonate acidizing, Transp. Porous Media 96, 1, 203–220. [Google Scholar]
All Tables
All Figures
Fig. 1 Comparison of (a) the wormhole structures simulated using the presented model with and (b) the experimental observation reported by Seagraves et al. [49]. 

In the text 
Fig. 2 Dissolution patterns obtained from simulations by using of Darcy and Forchheimer framework. 

In the text 
Fig. 3 Breakthrough curves obtained from simulations with of Darcy and Forchheimer framework. 

In the text 