An Integrated Model for Acid Fracturing without Prepad Considering Wormhole Growth

Abstract

Acid fracturing is an effective stimulation technology that is widely applied in carbonate reservoirs. An integrated model for acid fracturing without prepad treatment has been established. Compared with the previous models which use prepad for generating hydraulic fractures, this model can simultaneously simulate the fracture propagation and the acid etching of fracture surfaces, as well as the wormhole growth during acid fracturing. The influences of some essential factors have been studied through a series of numerical simulations, and the main conclusions are as follows. First, increasing the injected acid volume can expand the size of the formed hydraulic fractures and extend the propagation distance of the wormhole. Increasing the injected acid volume can also expand the etched width and extend the effective distance of the injected acid. Second, a high injection rate impels more acid to flow into the depth of a fracture before infiltration and reaction, resulting in the augmentation of a hydraulic fracture’s geometric size and the extension of the effective distance. But the maximum etched width decreases as the injection rate rises. A high injection rate can also enable wormholes to grow in the natural fracture area farther away from the hydraulic fracture inlet, but shorten the length of the original wormhole near the hydraulic fracture inlet. Third, an increase in acid viscosity can enlarge the geometric size of the hydraulic fracture and reduce the propagation distance of wormholes. In addition, an increase in the acid viscosity blocks the acid flow from fracture inlet to tip, reducing the effective distance of acid fracturing. Fourth, the natural fracture is the vital inducement of wormhole growth, and wormholes are apt to grow in the natural fracture area. Moreover, the geometric size of the hydraulic fracture and the effective distance of acid fracturing decrease with an increasing number of natural fractures. This research can provide a reference for field applications of acid fracturing without prepad.

1. Introduction

Acid treatment technology is an essential stimulation measure and has been used to improve well performance for a long time, especially in carbonate reservoirs [1,2]. Acid treatment technology can be further classified into different technologies according to the construction method.

When the stimulation purpose is to increase or restore permeability near the wellbore, matrix acidizing is the preferred stimulation method [3,4]. Acid is injected into the reservoir at pressures below the formation breakdown pressure and dissolves the carbonate rock non-uniformly, forming several branched channels. These channels, referred to as wormholes, can penetrate the damage zone and improve flow conditions. However, the stimulation result of matrix acidizing is often unsatisfactory in low permeability reservoirs, and it is required to adopt a more aggressive acid treatment, acid fracturing. In the construction of acid fracturing, acid is injected into the reservoir at pressures higher than the formation breakdown pressure in order to create an artificial fracture or open natural fractures. Unlike the artificial fracture, propped by proppants in conventional hydraulic fracturing [5,6,7], the artificial fracture in acid fracturing is etched by the acid in non-uniform patterns and generates a certain conductivity. These etched fractures can enhance fluid flow from the reservoir into the wellbore and increase the well’s production significantly.

A successful acid fracturing treatment requires an appropriate acid fracturing design. However, because acid fracturing involves several physical and chemical processes, the mechanism of acid fracturing is more complex than traditional hydraulic fracturing [8,9,10]. This undoubtedly brings a great challenge to the acid fracturing design. To clear the acid fracturing mechanism and optimize the acid fracturing design, many scholars have conducted relevant research and established numerous acid fracturing models. Because the acid fracturing process can be divided into hydraulic fracture propagation, the etching of the hydraulic fracture and wormhole growth, previous studies mainly focused on these aspects. By taking the fracture geometry calculated using the hydraulic fracturing simulator as the initial geometry of the acid-etched fracture, Oeth et al. [11] realized the numerical stimulation of acid etching of fractures at an engineering scale. To study the influence of fracture surface morphology on the flow and reaction of acid fluids, Mou et al. [12] established the initial fracture geometry using spatial random numbers to characterize the heterogeneity of the carbonate reservoir. Williams et al. [13] built an early 1D acid etching fracture model. They calculated the acid concentration in the fracture using an analytical solution of the heat transfer model for parallel plates. Lo and Dean [14] introduced a 1D acid solution coupled with a multiple-fluid geometry model. Roberts et al. [15] developed a 2D steady-state acid etching fracture model by considering the convection effect in the fracture length direction and the diffusion effect in the fracture width direction. They also introduced the concept of the mass transfer coefficient in order to facilitate the solution. Settari et al. [16] built a transient 2D acid fracturing model, where the acid concentration was resolved along the fracture width. On their basis, Ugursal et al. [17] established a 2D acid fracturing model, which considers natural fractures and enables the simulation of the acid fracturing treatment in naturally fractured limestone reservoirs. Romero et al. [18] had earlier proposed a pseudo-3D acid fracturing model by assuming a 3D analytical fluid flow velocity profile. Aljawad et al. [19] developed a pseudo-3D acid fracturing model by updating the velocity profile with the Berman flow function equation. Their model can additionally consider the influence of formation heterogeneity and temperature fields on the reaction. Mou et al. [20] proposed a mesoscale steady-state 3D acid fracturing model, which constructs the 3D hydraulic fracture geometry using a geostatistical method to characterize the mineral’s heterogeneity. Oeth et al. [21] improved Mou et al.’s research to establish a fully 3D numerical acid fracturing model, which adopts the semi-implicit method for pressure-linked equations (SIMPLE) algorithm to resolve the acid velocity profile. Alhubail et al. [22] proposed a transient 3D acid fracturing model where the temperature field is considered. They realized the solution of the model using the finite element method.

As illustrated in Figure 1, the acid leak-off in acid fracturing will generate acid wormholes in reservoirs and shorten the effective distance of the acid. On the other hand, these wormholes can accelerate the acid leak-off significantly, and the leak-off acid will further promote the growth of wormholes. However, these acid fracturing models mentioned above are unable to simulate the wormhole growth on the fracture walls. These previous acid fracturing models mostly adopt a modified leak-off coefficient to signify the influence of wormholes on the acid leak-off. This simplified treatment method will undoubtedly reduce the accuracy of acid fracturing simulations. Additionally, most previous acid fracturing models are only applicable to prepad acid fracturing. These models assume that prepad is used to create an artificial fracture first, then the acid is injected into the fracture and etches the fracture wall. Nevertheless, in small acid fracturing constructions, acid is often injected directly at pressures higher than the formation breakdown pressure instead of injecting prepad first. Acid takes both roles of creating hydraulic fracture and etching the fracture wall. Thus, the previous acid fracturing models are inapplicable to the acid fracturing without prepad.

In this paper, we established a mathematical model for acid fracturing without prepad treatment, which also incorporates the wormhole growth module. In contrast with the sequential simulation of the previous models, this model can simultaneously simulate the fracture propagation and acid etching of the fracture surface, making it more applicable to acid fracturing without prepad treatment. Then, based on the established model, a series of sensitivity analyses were conducted to investigate the influences of various factors on acid fracturing. The analysed factors include the injection rate of acid, the injected acid volume, the number of natural fractures and the acid’s viscosity.

2. Numerical Model

This integrated acid fracturing model consists of three sub-models: a fracture propagation model, an acid etching model and a wormhole growth model. These sub-models and the solution method for the integrated acid fracturing model are successively described in the following section.

2.1. Fracture Propagation Model

Acid flow, reaction and leak-off are all carried out in the fracture. At each time step, the fracture propagation model provides the solution domain for the acid etching model and the wormhole growth model. Thus, the accuracy of the acid fracturing simulation is directly determined by the accurate calculation of the artificial fracture geometry. To balance the simulation accuracy and calculation time, the typical pseudo-3D fracture propagation model under symmetric stresses is adopted in this integrated model. The geometry of the hydraulic fracture is shown in Figure 2, and the adopted pseudo-3D fracture propagation model can be described as follows [23].

The continuity equation of the fracture propagation model is written as Equation (1):

$$-\frac{\partial q_{\mathrm{f}}(x,t)}{\partial x}=v_{\mathrm{l}}\left(x,t\right)+\frac{\partial A_{\mathrm{f}}(x,t)}{\partial t}$$

(1)

where x is the distance in the fracture length direction, m; t is the construction time, min; q_f (x, t) is the rate of flow at position x at time t, m³/min; v_l (x, t) is the filtration rate at position x at time t, m³/min; and A_f (x, t) is the cross-sectional area of the hydraulic fracture at position x at time t, m².

The pressure drop equation of fluid in the fracture can be simplified as Equation (2):

$$\frac{\partial p(x,0,t)}{\partial x}=-\frac{64{\mu }_{\mathrm{a}}{q}_{\mathrm{f}}(x,t)}{\pi H(x,t)W_p{(x,0,t)}^3}$$

(2)

where p (x, 0, t) is the pressure at the fracture centre at position x at time t, MPa; μ_a is the acid viscosity, mPa·s; H (x, t) is the fracture height at position x at time t, m; and W_p (x, 0, t) is the fracture width at the fracture centre at position x at time t, m.

The fracture width equation is shown as Equation (3):

$$W_p(x,y,t)=\frac{8\left(1-{v_{\mathrm{p}}}^2\right)H(x,t)}{\pi E}{\displaystyle {\int }_{{\eta }_{\mathrm{f}}}^1\frac{\tau \mathrm{d}\tau }{\sqrt{{\tau }^2-{\eta }_{\mathrm{f}}^2}}}{\displaystyle {\int }_0^{\tau }\frac{p(y)\mathrm{d}y}{\sqrt{{\tau }^2-y^2}}}$$

(3)

of which:

$${\eta }_{\mathrm{f}}=\frac{y}{H(x,t)/2}$$

where y is the distance in the fracture height direction, m; W(x, y, t) is the fracture width at position (x, y) at time t, m; v_p is Poisson’s ratio of the rock; E is the Young’s modulus of the rock, MPa; and p (y) is the pressure at the fracture height y, MPa.

The fracture height equation is written as Equation (4):

$$-\frac{\partial p}{\partial x}=\frac{1}{H}\left(\frac{K_I}{2\sqrt{\pi H_{\mathrm{half} }}}-\frac{2}{\pi }\Delta s\frac{f}{\sqrt{1-f^2}}\right)\frac{\mathrm{d}{H}_{\mathrm{half} }}{\mathrm{d}x}$$

(4)

where K_I is the stress intensity factor, MPa·m^0.5; H_half is the half-height of the fracture, m; ∆s is the stress difference between the reservoir and the interbed, MPa; and f is the ratio of reservoir thickness to fracture height.

By solving Equations (1)–(4) with the Runge–Kutta method and iterative calculation, the fracture geometry at each time can be obtained. Since the pseudo-3D fracture model has been introduced in numerous articles, the specific solution steps of the fracture propagation model are no longer discussed in detail.

2.2. Acid Etching Model

The acid etching model is used to simulate the acid flow and reaction in the fracture and can be divided into three parts: a mass conservation model, an acid transport model and a fracture width change model. In the simulation process, firstly, the mass conservation model is applied to solve the pressure and velocity fields. Then, the acid concentration distribution is calculated using the acid transport model. Finally, the acid etched width distribution is obtained using the fracture width change model.

2.2.1. Mass Conservation Model

As shown in Figure 3, a control volume is taken in the hydraulic fracture and parallelepiped with dimensions of Δx, Δy and W. By conducting a mass balance analysis with the control volume, the mass conservation equation of the acid flow in the fracture can be written as Equation (5).

$$-\frac{\partial \left(v_{\mathrm{f}x}W\right)}{\partial x}-\frac{\partial \left(v_{\mathrm{f}y}W\right)}{\partial y}-2v_{\mathrm{l}}=\frac{\partial W}{\partial t}$$

(5)

where v_fx is the flow rate of acid in the fracture length direction, m/s; and v_fy is the flow rate of acid in the fracture height direction, m/s.

Based on Darcy’s law and the theory of the flow between parallel plates, the velocities of acid in the hydraulic fracture can be expressed as Equations (6) and (7).

$$v_{\mathrm{f}x}=-\frac{W^2}{12{\mu }_{\mathrm{a}}}\frac{\partial p}{\partial x}$$

(6)

$$v_{\mathrm{f}y}=-\frac{W^2}{12{\mu }_{\mathrm{a}}}\frac{\partial p}{\partial y}$$

(7)

Based on Equations (5)–(7), the mass conservation equation of acid can be written as Equation (8).

$$\frac{1}{12{\mu }_{\mathrm{a}}}\frac{\partial }{\partial x}\left(W^3\frac{\partial p}{\partial x}\right)+\frac{1}{12{\mu }_{\mathrm{a}}}\frac{\partial }{\partial y}\left(W^3\frac{\partial p}{\partial y}\right)-2v_{\mathrm{l}}=\frac{\partial W}{\partial t}$$

(8)

Assuming the acid flow is steady and the fracture width is constant during a time step, Δt, Equation (8) can be reformulated as Equation (9).

$$\frac{1}{12{\mu }_{\mathrm{a}}}\frac{\partial }{\partial x}\left(W^3\frac{\partial p}{\partial x}\right)+\frac{1}{12{\mu }_{\mathrm{a}}}\frac{\partial }{\partial y}\left(W^3\frac{\partial p}{\partial y}\right)-2v_{\mathrm{l}}=0$$

(9)

2.2.2. Acid Transport Model

As long as the flow velocity of acid is obtained using the mass conservation model, the acid concentration profile can be acquired using the acid transport model. Because the acid accumulation in the control volume is equal to the amount of acid flowing into the control volume minus the amount of acid flowing out of the control volume, the expression for acid transport in the fracture can be obtained as Equation (10).

$$\frac{1}{12{\mu }_{\mathrm{a}}}\frac{\partial }{\partial x}\left(\overline{C}{W}^3\frac{\partial p}{\partial x}\right)+\frac{1}{12{\mu }_{\mathrm{a}}}\frac{\partial }{\partial y}\left(\overline{C}{W}^3\frac{\partial p}{\partial y}\right)-2\overline{C}{v}_{\mathrm{l}}-2k_g(\overline{C}-C_{\mathrm{w}})=\frac{\partial \left(\overline{C}W\right)}{\partial t}$$

(10)

where $\overline{C}$ is the acid concentration averaged across the fracture width, kmol/m³; k_g is the mass transfer coefficient in the hydraulic fracture, m/s; and C_w is the acid concentration on the fracture surface, kmol/m³.

The equation of the local reaction on the fracture surface is given by Equation (11).

$$k_{\mathrm{g}}(\overline{C}-C_{\mathrm{w}})=R\left(C_{\mathrm{w}}\right)$$

(11)

where R(C_w) is the dissolution rate of one-step irreversible reaction, kmol/(s·m²).

The rate of the acid–rock reaction varies with different rock mineral compositions. In this research, the model adopts the acid–rock reaction rate formula of limestone in the Iraqi Missan oil field, which is obtained from a dissolution experiment of rock samples in this area. According to a previous study [24], the dissolution rate of limestone by hydrochloric acid can be defined by Equation (12).

$$R\left(C_{\mathrm{w}}\right)=0.015{\mathrm{e}}^{-18616/RT}{\overline{C}}^{1.1997}$$

(12)

where R is the gas constant, 8.314 J/(mol·K); and T is the temperature, K.

Assuming the acid flow is steady and the fracture width is constant during a time step, Δt, Equation (10) can be rewritten as Equation (13).

$$\frac{1}{12{\mu }_{\mathrm{a}}}\frac{\partial }{\partial x}\left(\overline{C}{W}^3\frac{\partial p}{\partial x}\right)+\frac{1}{12{\mu }_{\mathrm{a}}}\frac{\partial }{\partial y}\left(\overline{C}{W}^3\frac{\partial p}{\partial y}\right)-2\overline{C}{v}_{\mathrm{l}}-2k_g(\overline{C}-C_{\mathrm{w}})=\frac{W\partial \left(\overline{C}\right)}{\partial t}$$

(13)

2.2.3. Fracture Width Change Model

When the acid concentration profile in the fracture is obtained using the acid transport model, the fracture width change caused by the acid etching can be calculated based on the acid dissolving power and the amount of acid available for the reaction.

In the control volume, the rock volume dissolved by the filtration and diffusion acid is shown in Equation (14).

$$\Delta V_{\mathrm{Rock}}=\frac{\alpha }{{\rho }_r\left(1-\varphi \right)}\Delta t\left(\eta 2v_{\mathrm{l}}\Delta x\Delta y\overline{C}+2k_g\Delta x\Delta y\overline{C}\right)$$

(14)

where ρ_r is the density of the rock, kg/ m³; ϕ is the rock’s porosity; and η is the fraction of the acid reacting with the rock at fracture walls before leaking off, and takes the value of 0 in most cases.

On the other hand, the volume change in the control volume caused by the acid–rock reaction can be expressed by Equation (15).

$$\Delta V_{\mathrm{CV}}=\Delta x\Delta y\left\{W(x,y,n+1)-W(x,y,n)\right\}$$

(15)

where W(x, y, n + 1) is the fracture width at position (x, y) at time n + 1, m; and W(x, y, n) is the fracture width at position (x, y) at time n, m.

Because the volume change in the control volume equals the dissolved rock volume, Equation (16) can be acquired directly.

$$\frac{\beta }{{\rho }_r\left(1-\varphi \right)}\Delta t\left(\eta 2v_l\Delta x\Delta y\overline{C}+2k_g\Delta x\Delta y\overline{C}\right)=\Delta x\Delta y\left\{W(x,y,n+1)-W(x,y,n)\right\}$$

(16)

By conducting several transformations and taking the limit of Equation (16) at Δt → 0, the fracture width change due to the acid etching can be estimated by Equation (17).

$$\frac{\partial W}{\partial t}=\frac{\beta }{{\rho }_r\left(1-\varphi \right)}\left(\eta 2v_{\mathrm{l}}\overline{C}+2k_g\overline{C}\right)$$

(17)

2.3. Wormhole Growth Model

The wormholes formed during the acid fracturing can significantly increase the acid filtration and then affect the effect of acid fracturing. In order to simulate the wormhole growth, the two-scale continuum model proposed by Panga et al. [25] is coupled into the acid fracturing model in this work. However, because the 3D simulation of wormhole growth is very complicated and time-consuming, this coupled wormhole growth model only conducted a 2D wormhole growth simulation to facilitate field use. The two-scale continuum model includes two sub-models, Darcy’s scale model and a pore-scale model. Darcy’s scale model describes the flow of acid in rock, and the pore-scale model describes the property changes in microporous media. In the simulation of wormhole growth, the reservoir is meshed, as illustrated in Figure 4, and then Darcy’s scale model is used to obtain the pressure field, flow velocity field and acid concentration field in formation. Finally, the pore-scale model is applied to calculate the changes in microphysical fields such as the permeability field.

2.3.1. Darcy’s Scale Model

Acid flows in Darcy’s scale porous media, and corrodes rock, resulting in changes in physical parameters such as porosity and permeability. For describing the specific process, it is necessary to calculate the distributions of pressure, velocity and acid concentration in the reservoir using Darcy’s scale model. Darcy’s scale model consists of the Darcy equation, a material balance equation, a convection–diffusion equation and a porosity equation.

Darcy’s equation describes the flow velocity of acid and is directly derived from Darcy’s law, as shown in Equation (18).

$$\left(v_{\mathrm{r}x},v_{\mathrm{rz}}\right)=-\frac{K}{{\mu }_{\mathrm{a}}}\cdot \left(\frac{\partial P}{\partial x},\frac{\partial P}{\partial z}\right)$$

(18)

Through material balance analysis, the material balance equation can be given directly in Equation (19).

$$-\frac{\partial {\rho }_{\mathrm{a}}\varphi }{\partial t}=\frac{\partial {\rho }_{\mathrm{a}}{v}_{\mathrm{r}x}}{\partial x}+\frac{\partial {\rho }_{\mathrm{a}}{v}_{\mathrm{rz}}}{\partial z}$$

(19)

Combining Darcy’s Equation (18) with Equation (19), the material balance equation can be rewritten as Equation (20), which is used to acquire the distributions of pressure and velocity in the reservoir.

$$\frac{\partial }{\partial x}(\frac{K}{{\mu }_{\mathrm{a}}}\frac{\partial P}{\partial x})+\frac{\partial }{\partial z}(\frac{K}{{\mu }_{\mathrm{a}}}\frac{\partial P}{\partial z})=C_{\mathrm{p}}\varphi \frac{\partial P}{\partial t}$$

(20)

where ρ_a is the density of the acid, kg/m³; and C_p is the comprehensive compressibility of the reservoir, MPa⁻¹.

After the pressure and velocity fields in the reservoir are obtained using the material balance equation, the concentration distribution of the acid can be calculated using the convection–diffusion equation, shown as Equation (21).

$$\frac{\partial \left(\varphi C_f\right)}{\partial t}+\frac{\partial }{\partial x}\left(v_{\mathrm{r}x}{C}_f\right)+\frac{\partial }{\partial z}\left(v_{\mathrm{rz}}{C}_f\right)=\frac{\partial }{\partial x}\left(\varphi D_{ex}\frac{\partial C_f}{\partial x}\right)+\frac{\partial }{\partial z}\left(\varphi D_{ez}\frac{\partial C_f}{\partial z}\right)-R\left(C_s\right)a_v$$

(21)

where C_f is the cup-mixing concentration of the acid in the liquid phase, kmol/m³; D_e is the effective dispersion coefficient in the acid phase, m²/s; a_v is the interfacial area available for reaction per unit volume of the medium; R(C_s) is the rate of the dissolution reaction, kmol/(s·m²); and C_s is the cup-mixing mass concentration of the acid at the fluid–solid interface, kmol/m³.

In terms of acid–rock reaction, the reaction rate can be calculated using Equation (22).

$$R\left(C_s\right)a_v=\frac{k_c{k}_s}{k_c+k_s}{C}_f{a}_v$$

(22)

where k_s is the reaction rate constant, m/s; and k_c is mass transfer coefficient, m/s.

Thus, once the acid concentration distribution is obtained, the dissolution of rock can be computed using the porosity equation, written as Equation (23).

$$\frac{\partial \varphi }{\partial t}=\frac{k_c{a}_v{C}_f\left(\frac{k_s}{k_s+k_c}\right)\alpha }{{\rho }_r}$$

(23)

where α is the dissolving power of the acid, defined as grams of dissolved solid per mole of reacted acid, kg/mol.

2.3.2. Pore-Scale Model

The dissolution of rock makes the rock properties changes dynamically during the acid fracturing. In order to conduct an accurate simulation, it is necessary to update rock properties in time. The pore-scale model describes the pore volume change and is composed of the calculation of coefficients and structure–porosity relation.

The structure–porosity relation is the relationship between porosity, permeability, average pore radius and pore specific surface. After the change in porosity is obtained using the porosity equation, the change in each structural property is computed using the structure–porosity relation, shown as Equations (24)–(26).

$$\frac{K}{K_0}=\frac{\varphi }{{\varphi }_0}{\left[\frac{\varphi \left(1-{\varphi }_0\right)}{{\varphi }_0\left(1-\varphi \right)}\right]}^{2{\beta }_{\mathrm{w}}}$$

(24)

$$\frac{r_p}{r_{p0}}=\sqrt{\frac{K{\varphi }_0}{K_0\varphi }}$$

(25)

$$\frac{a_v}{a_0}=\frac{\varphi r_{p0}}{{\varphi }_0{r}_P}$$

(26)

where K₀ is the initial permeability of the reservoir, 10⁻³ × μm²; r_p0 is the initial pore radius of the medium, m; β_w is a constant that depends on the structure of the medium; ϕ₀ is the initial porosity of the reservoir; and a₀ is the initial pore-specific surface, m²/m³.

The mass transfer coefficient and the dispersion coefficient are dispensable calculation parameters in Equations (21) and (22), while their values are related to the microstructure of the rock. The calculation methods of the mass transfer coefficient and the dispersion coefficient are given by Equations (27)–(29).

$$Sh=\frac{2k_c{\overline{r}}_p}{D_m}=S{h}_{\infty }+0.7{\mathrm{Re}}_p^{1/2}S{c}^{1/3}$$

(27)

$$D_{ex}={\alpha }_{os}{D}_m+\frac{{\lambda }_x\left|U\right|r_p}{\varphi }$$

(28)

$$D_{ez}={\alpha }_{os}{D}_m+\frac{{\lambda }_z\left|U\right|r_p}{\varphi }$$

(29)

where Sh is the Sherwood number and represents the dimensionless mass transfer coefficient; D_m is the effective molecular diffusivity of acid, m²/s; Sh_∞ is the asymptotic Sherwood number; Re_p is the pore Reynolds number; Sc is the Schmidt number; D_ex and D_ez are the effective dispersion coefficients in x and z directions, m²/s; |U| is the acid flow velocity, m/s; α_os, λ_x and λ_z are constants that depend on the pore geometry and take the values of 0.5, 0.5 and 0.1, respectively, for a packed bed of spheres.

2.4. Solution Method of the Integrated Model

The integrated acid fracturing model established in this work contains three sub-models. Therefore, these sub-models need to be successively solved in the simulation process. At the beginning of a time step, the fracture dimensions are calculated using the fracture propagation model. The filtration rate of acid, v_l, is acquired using Equation (30).

$$v_{\mathrm{l}}(x,t)=\frac{K}{{\mu }_{\mathrm{a}}}\frac{p(x,t)-P(x,\Delta z_1,t)}{\Delta z_1}$$

(30)

where v_l(x,t) is the filtration rate of acid at position x along the fracture length at time t, m/s; p(x,t) is the fracture pressure at position x along the fracture length at time t, MPa; P(x, ∆z₁, t) is the reservoir pressure in the grid row closest to the fracture at time t, MPa; and ∆z₁ is the width of the first grid row, m.

Then, the fracture geometry is provided to the acid etching model as the solution domain. The profiles of pressure, acid concentration and etched width are calculated in turn during the acid etching simulation. The initial conditions and boundary conditions of the acid etching model are shown as Equations (31) and (32).

Initial Conditions:

$$\left\{\begin{array}{l}p(x,y)=p_r,\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\forall x,y,t=0\hfill \\ \overline{C}(x,y)=0,\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\forall x,y,t=0\hfill \\ W(x,y)=W_p(x,y),\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\forall x,y,t=0\hfill \end{array}\right.$$

(31)

Boundary Conditions:

$$\left\{\begin{array}{l}{\left.{\displaystyle {\int }_{-H_{\mathrm{half}}}^{H_{\mathrm{half}}}\frac{W^3}{12{\mu }_{\mathrm{a}}}}\frac{\partial p}{\partial x}\right|}_{x=0} \mathrm{d}y=q_{\mathrm{inj}},\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}t>0\hfill \\ p(L_{\mathrm{f}},y)=p_r+p_{\mathrm{tip}},\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\forall y,t>0\hfill \\ {\left.\frac{\partial p}{\partial y}\right|}_{y=-H_{\mathrm{half}}}=0,\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\forall x,t>0\hfill \\ {\left.\frac{\partial p}{\partial y}\right|}_{y=H_{\mathrm{half}}}=0,\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\forall x,t>0\hfill \\ \overline{C}(0,y)=C_0,\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}\forall y,t>0\hfill \end{array}\right.$$

(32)

where W_p is the fracture width provided by fracture propagation model, m; L_f is fracture length, m; p_r is original formation pressure, MPa; q_inj is the injection rate of acid, m³/min; p_tip is the net pressure of the fracture tip provided by the fracture propagation model, MPa; and C₀ is the initial concentration of acid, kmol/m³.

After the profiles of pressure and acid concentration are acquired, they are transmitted into the wormhole growth model. In order to simulate the wormhole growth in this time step, the model takes the pressure and acid concentration distributions along the fracture as the indispensable boundary conditions. The specific initial conditions and boundary conditions of the wormhole growth model are shown in Equations (33) and (34).

Initial Conditions:

$$\left\{\begin{array}{l}P(x,z)=p_r,\hspace{0.33em}\hspace{0.33em}\forall x,z,t=0\\ C_f(x,z)=0,\hspace{0.33em}\hspace{0.33em}\forall x,z,t=0\end{array}\right.$$

(33)

Boundary Conditions:

$$\left\{\begin{array}{l}P(x,0)=p(x,0),\hspace{0.33em}\hspace{0.33em}\forall x,t>0\\ C_f(x,0)=\overline{C}(x,0),\hspace{0.33em}\hspace{0.33em}\forall x,t>0\\ {\left.\frac{\partial P}{\partial x}\right|}_{x=0,x=L_{\mathrm{f}}}=0\\ {\left.\frac{\partial C_f}{\partial x}\right|}_{x=0,x=L_{\mathrm{f}}}=0\end{array}\right.$$

(34)

After the rock property fields at the end of this time step are obtained using the wormhole growth model, the acid fracturing simulation of the current time step is completed. Then, the obtained physical properties of the reservoir are transported into the next time step and utilized as essential parameters to calculate the acid filtration rate. This simulation procedure is carried out at each time step until the final injection time is reached.

3. Case Calculation and Sensitivity Analysis

This section may be divided by subheadings. It should provide a concise and precise description of the experimental results, their interpretation, as well as the experimental conclusions that can be drawn.

In this part, the influence of different factors on acid fracturing are analysed through a series of sensitivity analyses. In order to take natural fractures into account, a characterization method of natural fractures is adopted, which is similar to the treatment methods in previous studies. As illustrated in Figure 5, Ren et al. [26,27] simplify these interconnected nature fractures to a single nature fracture perpendicular to the hydraulic fracture. According to the number of natural fractures, we evenly select several grid columns within 30 m from the hydraulic fracture inlet and endow them with higher porosity, 0.6–0.7, to represent natural fractures. Meanwhile, the porosity of the reservoir matrix is set to fluctuate between 0.2 and 0.3 to reflect formation heterogeneity. In the analysis, the effective area of treatment is defined as the area where the etched width is greater than 1 mm and the farthest distance from the fracture inlet to the area is defined as the effective distance of treatment. Meanwhile, 8% acid concentration is defined as the acid concentration boundary, and the area with an acid concentration greater than this standard is considered as the high concentration area. The borders of the effective area of treatment and the high concentration area are marked with red lines in the corresponding figures. These analysed factors include injection rate, acid viscosity, number of natural fractures and the amount of injected acid. The values of basic calculation parameters used in the simulation are shown in

Table 1. Basic calculation parameters.

Parameter	Value	Unit
Injection rate, qinj	5	m3/min
Temperature of injected acid, T	293.15	K
The amount of injected acid, Q	100	m3
Reservoir thickness, Hre	23	m
Young’s modulus of the rock, E	2.59 × 104	MPa
Poisson’s ratio of the rock, vp	0.22	-
Acid viscosity, μa	30	mPa·s
Density of the rock, ρr	2400	kg/m3
Original formation pressure, pr	59	MPa
Dissolving power of the acid, α	1.37	-
Effective dispersion coefficient, De	3.6 × 10−9	m2/s
Initial reservoir porosity, ϕ0	0.2	-
Initial reservoir permeability, K0	20	10−3 × μm2
Initial pore radius of the medium, rp0	1 × 10−5	m
Initial pore specific surface, a0	5000	m2/m3
Length of time step, ∆t	0.5	min
Size of x-direction grinds, ∆x	0.25	m
Size of z-direction grinds, ∆z	0.01	m
Reaction rate constant, ks	0.2 × 10−2	m/s
Acid initial concentration, C0	15	%
Number of natural fractures, Nf	3	-
Comprehensive compressibility, Cp	6 × 10−3	MPa−1

.

3.1. Influence of Injected Acid Amount

In this part, the other parameters take the values of basic calculation parameters, and the amount of injected acid, Q, is changed gradually to study its effect on acid fracturing.

As shown in Figure 6, the hydraulic fracture propagates with the increase in injected acid amount. As the injected acid amount increases from 20 m³ to 100 m³, the length of the hydraulic fracture ascends from 49.1 m to 62.5 m and the max fracture height ascends from 29.5 m to 30.6 m. Thus, the increase in injected acid amount can enlarge the geometric size of the hydraulic fracture, but the rate of increase decreases with the increase in injected acid amount. Meanwhile, the acid concentration distribution is basically unaffected by the injected acid amount. Regardless of the injected acid amount, the acid concentration is highest at the centre of the fracture inlet, which equals the initial concentration, and gradually reduces along the fracture length and height.

From the simulation results in Figure 7, the etched width also descends along the fracture length and height. The area with the largest etched width always appears at the centre of the fracture inlet. By comparing the etched width distribution in different cases, it can be concluded that the increase in injected acid amount can strengthen the etching degree of the fracture surface. At a certain acid injection rate, the increase in injected acid amount will prolong the etching time, and the area of acid–rock reaction will spread due to a larger fracture geometry. With the increase in the injected acid amount from 20 m³ to 100 m³, the max etched width rises from 1.6 mm to 9.4 mm. Besides, it can be seen from Figure 8 that the effective distance of treatment expands significantly from 6.5 m to 49.5 m. Therefore, the enlargement of injected acid amount can increase the etched width and broaden the effective scope of the acid fracturing treatment.

Additionally, the increasing injected acid amount can augment the propagation length of the wormhole and enable the wormhole to grow in the farther natural fracture. As illustrated in Figure 9, almost all leak-off acid enters the natural fracture and induces the wormhole to grow in the natural fracture area. This simulation result coincides with the previous research. When the injected acid amount is 20 m³, the wormhole only grows in the first natural fracture, which is closest to the hydraulic fracture inlet. At this time, the wormhole length is only 2 cm. Then, with the increase in injected acid amount, wormholes slowly grow and also start to grow in the other natural fractures. When the injected acid amount rises up to 100 m³, the wormhole length in three natural fractures reaches 3 cm, 5 cm and 7 cm. Thus, the increase in injected acid amount can enhance wormhole growth.

3.2. Influence of Injection Rate

The injection rate is an essential construction factor, and a reasonable injection rate can improve the effect of the acid fracturing treatment. In this part, the other parameters are set to basic calculation parameters, and the injection rate of acid is gradually elevated to investigate the influence of injection rate on the acid fracturing effect.

As shown in Figure 10, the increase in injection rate can enlarge the hydraulic fracture geometry and boost the acid concentration in the depths of the hydraulic fracture. In the case of a 2 m³/min injection rate, the length and height of the hydraulic fracture are only 27 m and 24.6 m. As the injection rate rises to 10 m³/min, the length and height of the hydraulic fracture substantially increase to 95.2 m and 39 m. Meanwhile, a large injection rate can boost the acid concentration in a wide area. In the case of a 10 m³/min injection rate, the acid concentration at the fracture centre can keep over 8% within 41 m from the fracture inlet. But the acid concentration at the fracture centre only can keep over 8% within 22.5 m when the injection rate is 2 m³/min. This is because with the increase in the injection rate more acid solution can be transmitted to the fracture depths in a shorter time, reducing the leak-off and reaction of acid solution.

It can be concluded from Figure 11 and Figure 12 that the increase in the injection rate reduces the etched width near the fracture inlet, but it can enable the acid to effectively etch the fracture wall in a wider range. When the injection rate is 2 m³/min, the max etched width near the fracture inlet is 24 mm, and the effective distance of treatment is only 26.8 m. When the injection rate rises to 10 m³/min, more acid is transmitted to the fracture depths. Thus, the max etched width near the fracture inlet decreases to 4.5 mm. Additionally, the effective area of treatment greatly extends, and the effective distance increases to 40 m.

As illustrated in Figure 13, a large injection rate can make the wormhole propagate in more natural fractures, but it will also shorten the extension distance of the wormhole in the natural fractures near the hydraulic fracture inlet. When the injection rate is 2 m³/min, there is no wormhole growth in the third fracture. The wormholes extend to 9 cm and 5 cm, respectively, in the first two natural fractures. As the injection rate rises, the wormhole begins to grow in the third natural fracture, and the extension distances of the wormholes in the first two natural fractures decrease. The wormhole in the third natural fracture propagates to 2 cm, and the wormholes in the other natural fractures shorten to 4.5 cm and 3 cm when the injection rate reaches 10 m³/min.

3.3. Influence of Acid Viscosity

Acid viscosity is another important property of acid which can influence the flow and leak-off behaviours of acid. In this section, the other parameters are set as the basic calculation parameters, and the cases with different acid viscosities are simulated to investigate the influence of acid viscosity.

As shown in Figure 14, the length and height significantly increase with the increasing acid viscosity. The increase in the fracture height is mainly due to the fact that high-viscosity fracturing fluid is more likely to generate a short and high hydraulic fracture. Nevertheless, the increase in acid viscosity also decreases the amount of leak-off acid simultaneously, and more acid solution can contribute to the propagation of the hydraulic fracture. Therefore, the length of the hydraulic fracture does not reduce, but further increases. Meanwhile, the increase in acid viscosity hinders the flow of acid to the fracture tip and decreases the acid concentration in the depth of the hydraulic fracture at the end of the injection. When the acid viscosity is 30 mPa·s, the acid concentration at the fracture centre can keep over 8% within 30.5 m from the fracture inlet. But when the acid viscosity rises to 110 mPa·s, the acid concentration at the fracture centre can keep over 8% within only 20 m from the fracture inlet.

According to Figure 15 and Figure 16, the increase in acid viscosity cannot extend the effective distance of acid fracturing treatment. On the contrary, the effective distance decreases slowly with the increase in viscosity. In the case of the 30 mPa·s acid viscosity, the effective distance is 49.5 m. When the acid viscosity increases to 110 mPa·s, the effective distance shortens to 37 m. Unlike the effective distance of acid fracturing treatment, the max etched width is not affected by the acid viscosity and remains at 9.3 mm.

It can be seen directly from Figure 17 that the propagation length of wormholes decreases with the increase in the acid viscosity. This is because the increase in acid viscosity reduces the filtration of acid from the hydraulic fracture into formation, thereby reducing the acid amount that dissolves the formation rock to form wormholes. When the acid viscosity is 30 mPa·s, the wormholes in the three natural fracture areas have the maximum propagation distances, which are 3 cm, 4 cm and 7 cm, respectively. With the increase in the acid viscosity, the propagation distance of the three wormholes decreases gradually. When the acid viscosity reaches 90 mPa·s, there is almost no wormhole growth in the third natural fracture area. The wormhole in the first natural fracture area shortens from 7 cm to 2 cm when the acid viscosity finally rises to 110 mPa·s.

3.4. Influence of Natural Fracture Number

During the acid fracturing construction, the presence of natural fractures significantly increases the magnitude of the acid leak-off and induces the growth of the wormhole. In this part, the influence of the natural fracture number is studied by changing the natural fracture number continuously and keeping the other parameters to the basic calculation parameters. It is assumed that no matter how the number of natural fractures changes, the natural fractures are evenly distributed in the area within 30 m from the inlet of hydraulic fracture.

As illustrated in Figure 18, the geometry of hydraulic fracture gradually shrinks with the increasing number of natural fractures. When there is no natural fracture, the length and height of the hydraulic fracture formed at the end of the acid injection are 109.2 m and 34.2 m, respectively. The increase in natural fracture number causes the augmentation of acid filtration, which leads to the reduction in hydraulic fracture size. When the number of natural fractures reaches four, the length and height of the hydraulic fracture reduce to 53.9 m and 29.8 m separately. Additionally, compared with the other areas of the hydraulic fracture, natural fractures also make the acid concentration more average in the intersection area. Several abnormal acid concentration areas can be seen in Figure 18. The acid concentration in these areas is more average due to the existence of natural fractures.

In addition, as shown in Figure 19 and Figure 20, the effective distance of acid fracturing decreases gently as the number of natural fractures increases. When there is no natural fracture, the fracture surface within 50 m from the hydraulic fracture inlet can be effectively etched. When the number of natural fractures increases to four, the effective distance of the acid fracturing treatment drops to 46 m. However, regardless of the number of natural fractures, the maximum etched width appears at the inlet of the hydraulic fracture, and its value remains 9.4 mm.

It can be concluded from Figure 21 that the natural fracture is the cardinal inducement for the growth of wormholes in the formation. When there is no natural fracture, the acid solution can only penetrate into the formation through the rock matrix, and the amount of leak-off acid is very low. There is almost no wormhole growth in the reservoir under this condition. Natural fractures act as the dominant channels for acid solution filtration. When natural fractures exist, the amount of the leak-off acid greatly increases, and wormholes begin to grow in the natural fracture area.

4. Conclusions and Recommendation

In this paper, we improve previous studies and establish a numerical model for acid fracturing without prepad, which can simultaneously simulate the fracture propagation and acid etching of fracture surfaces. Additionally, the wormhole growth model is also coupled into the numerical model. Therefore, the established model can additionally simulate the wormhole growth during the acid fracturing. Then, a series of sensitivity analyses is conducted to investigate the effects of various factors, including the amount of injected acid, the injection rate of acid, the acid viscosity and the number of natural fractures. The main conclusions are as follows.

(1)

Increasing the amount of injected acid can not only enlarge the geometric size of the formed hydraulic fracture, but also increase the etched width and extend the effective distance of acid fracturing treatment. Additionally, the propagation distance of the wormhole also extends with the increase in the injected acid amount.

(2)

The increase in the injection rate impels more acid solution to flow into the depths of the hydraulic fracture before leak-off and reaction. Thus, the increase in the injection rate can enlarge the hydraulic fracture geometry and boost the acid concentration away from the fracture inlet. The increase in the injection rate can extend the effective distance of acid fracturing treatment, but also shorten the total injection time, which leads to a decrease in maximum etched width. Meanwhile, with the increase in the injection rate, the wormhole starts to grow in the natural fracture area farther away from the hydraulic fracture inlet, but the length of the original wormhole in the natural fracture area near the inlet is shortened.

(3)

The adoption of high viscosity in acid fracturing without prepad can increase the geometric size of the hydraulic fracture. However, because the high viscosity of acid impedes the flow of acid from the fracture inlet to the fracture tip, it will reduce the acid concentration in the fracture depth and shorten the effective distance of acid fracturing. In addition, the increase in acid viscosity reduces the propagation length of wormholes.

(4)

Natural fractures are the cardinal inducement for the growth of wormholes in the formation, and there is almost no wormhole growth in the reservoir without natural fractures. Wormholes are apt to grow in the natural fracture area. Additionally, with the increasing number of natural fractures, the geometry of the hydraulic fracture gradually shrinks, and the effective distance of acid fracturing decreases.

(5)

In contrast with previous models for acid fracturing, this model can simultaneously simulate the fracture propagation and acid etching of fracture surfaces. Additionally, it can additionally simulate wormhole growth during acid fracturing. In the design of acid fracturing without prepad, this model can be taken to predict the geometry of hydraulic fractures, acid concentration profiles, etched width profiles and the growth of wormholes.