Three-dimensional simulation of acidizing process in carbonate rocks using the Darcy–Forchheimer framework

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 3-D reactive flow model with non-Darcy framework is developed based on the two-scale 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 non-Darcy flow on dissolution pattern and breakthrough volume is analyzed. It is found that there is no effect of non-Darcy 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 non-Darcy flow or not.


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][3][4][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][9][10]; (2) capillary tube model [11][12][13][14][15]; (3) network model [7,[16][17][18][19]; and (4) continuum model [6,[20][21][22][23]. Reviews of these models, along with summaries of their limitations and assumptions, can be found in [24][25][26][27][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][31][32], completion methods [33], the type of injected acid (hydrochloric acid and self-diverting acid) [34][35][36][37][38], medium heterogeneity (the presence of vugs and fractures) [20,23,39,40], and reservoir temperature [41][42][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 reactive-transport process. However, their work is limited to the 2-D condition, which cannot accurately represent the actual acidizing operation. The main aim of this paper is to extend the reactivetransport model to 3-D non-Darcy 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 non-Darcy flow on dissolution patterns and breakthrough volume.

Mathematical model
The two-scale continuum model developed by Panga et al. [26] consists of a Darcy-scale model and a pore-scale model. The Darcy-scale model describes the reactive-transport 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 Darcy-scale model can be expressed as follows: ð3Þ where, P is the pressure of the pore-fluid; l is the viscosity of the fluid phase; K is the permeability tensor, which is calculated from a pore-scale model or determined by lab measurement; u is the velocity vector; b is referred as the Forchheimer number that can be obtained experimentally; q is the density of the fluid; t is the time; / is the porosity of the rock; C f is the cup-mixing 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; a is the dissolving power of the acid, defined as grams of dissolved solid per mole of acid reacted; and q 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 pore-scale model that captures the evolution of these variables with the variation in porosity is needed to complete the Darcy-scale model. We adopt here a fractal theory based formulation to relate the permeability with the porosity as: with the definition of the fractal dimension of pore spaces and tortuosity, D f and D T , respectively as, where d E is the Euclidean dimension; k min and k max are the smallest pore diameter and the largest diameter, respectively. The typical value of k min =k 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: The mass transfer coefficient, and the effective dispersion coefficients, inside the pores are expressed by the relations: 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 =m, m is the kinematic viscosity; Sc is the Schmidt number defined as Sc ¼ m=D m ; D eX and D eT are the longitudinal and transverse dispersion coefficient, respectively. a os , k X , and k T are constants that depend upon the pore structure and their typical values are 0.5, 0.5, 0.1 for a packed-bed of spheres, respectively. A 3-D 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: 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;f 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 ; D/ 0 .

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. g is the pore-to-domain length ratio; a LH and a LW are the aspect ratios; j is the dimensionless permeability; B is the dimensionless non-Darcy coefficient. The parameters obtained from the non-dimensionalization are the pore scale Thiele modulus h 2 T , 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 U 2 is core-scale equivalent of the pore-scale 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 non-dimensionalized form are as follows: and the resulting boundary and initial conditions are,

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 reactive-transport equation (Eq. (21)) to calculate the acid concentration. The operator splitting method is adopted here for solving the reactive-transport 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: and the following equations in the second step: 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 pore-scale 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][47][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.

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 3-D 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.

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][51][52][53][54], it will not be repeated here. This work only focuses on the influence of non-Darcy flow on the dissolution structure and breakthrough volume. 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 non-Darcy 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 high-speed non-Darcy resistance term in the Darcy-Forchheimer framework. Wherein, this non-Darcy resistance term is equivalent to reducing the mobility of the high permeable area to make the injected acid flow into other low permeable areas. 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:

Breakthrough curve
where A is the area of the injection end; V is the bulk volume of the rock; t BT is the breakthrough time.
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.

Conclusion
In this paper, the effect of non-Darcy flow on dissolution structure and breakthrough volume is studied by numerically solving the 3-D two-scale continuum model with the Darcy-Forchheimer framework. The main results can be summarized as following: 1. 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. 2. 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. 3. There is no effect on optimal injection rate whether or not to consider non-Darcy flow influence. The optimal injection rate calculated using the Forchheimer equation and the Darcy equation is the same.