A new analytical model for flow in acidized fractured-vuggy porous media.

Acidizing is one of the widely used technologies that makes the development of naturally fractured-vuggy reservoirs effective. During the process of acidizing, carbonate minerals are dissolved by hydrochloric acid, which can create high conductivity channels and wormholes to connect fractures and pores. In this work, a new analytical model, incorporating the heterogeneity of the pore networks into acidizing region, is proposed to study the flow characteristics in acidized fractured-vuggy reservoirs. The model is coupled by an acidized inner region and a conceptualized outer region of common triple medium. The porosity and permeability of inner region, which are rather heterogeneous and disordered when observed at different length scales, can be well addressed by fractal theory. The properties of the outer region can be described with three basic parameters: the matrix block size LM, the space interval of fracture LF and the radius of the vug Lv. Results show that the flow characteristic curves can be characterized by six flow stages (i.e. wellbore storage stage, radial flow stage in the interior region, fracture-vug inter-porosity flow stage, transition flow stage, fracture-matrix inter-porosity flow stage and external boundary response stage). It can be applied to estimate reservoir parameters for uncertainty reduction using inverse modeling.

Introduction

In the past few decades, characterizing the fluid flow mechanism in fractured-vuggy porous media has become a hot topic in petroleum engineering, as increasing carbonate reservoirs are discovered worldwide[1-6]. It is reported that fractured-vuggy reservoirs have complex pore systems including matrix, fractures and vugs which cause the fluid flow characteristics within these reservoirs more complicated than that in conventional sandstones[7-11].

At present, acidizing or acid fracturing stimulation technology is essential for the development of fractured-vuggy carbonate reservoirs[12,13]. As acid is injected, carbonate minerals are dissolved and highly conductive channels or wormholes are usually formed[14]. Wormhole forming mechanism in carbonate porous media has been widely investigated using experimental and theoretical methods[15-21]. Many scholars have suggested acid flow in the pore networks selectively and create different types of wormhole structures[17,20]. Mostofizadeh and Economides[15] carried out radial acidizing tests on limestone rocks and suggested there exist an optimum injection rate. Frick et al.[16] conducted radial core experiments to investigate the interaction between hydrochloric acid and limestones and identified the optimum injection rate. Fredd[17] presented a dynamic mathematical model for wormhole formation determine, and suggested the derived model can be used to determining optimal acidizing strategies. Tremblay[18] proposed a sand transport model to calculate the wormhole growth. Liu et al.[19] derived a coupled model to study the wormhole propagation behavior under reservoir condition and they concluded that, due to the effect of the compressed zone, there existed a maximum wormhole length for both constant rate condition and constant pressure condition. Zakaria[20] conducted various laboratory tests to study the acid stimulation treatments in carbonate rocks and suggested that acid stimulation treatments are significantly affected by the pore heterogeneity and the acid propagate through the fractured-vuggy carbonate rocks in a selective manner. Li et al.[21] studied the moving interface of acid fingering using both the experimental tests and numerical method. Dejam et al.[22-27] formulated mathematical models to study the solute transport and dispersion in porous media with coupled systems. In addition, a theoretical model was derived by Ugursal[28] to evaluate acid fracturing performances in carbonate rocks.

Recently, the fluid flow and transport behavior in fractured-vuggy porous media have been analyzed with theoretical methods[3,9,29-34]. Wu et al.[29-31] investigated multiphase flow in fractured-vuggy reservoirs using triple-continuum models. Yao et al.[32] presented discrete fracture-vug network model (DFVN) modeling fluid flow in carbonate rocks. However, to the best of my knowledge, reliable characterization of actual pore networks in acidized region is severely limited, and the corresponding analytical models for acidizing wells with wormholes are scarce. In 2014, Wang et al.[33] derived a theoretical model considering acid wormholes to recognize the flow characteristics. And in their model, the pore networks in acidized region were conceptualized as multi-branched wormholes (i.e., the pore networks in acidized region were simplified as multi-branched high-permeability channels), which were used to bridge the gap between the outer region of triple medium and the wellbore. In 2018, Wang et al.[34] also investigated flow characteristics in fractured-vuggy carbonate reservoirs using semi-analytical model considering wormholes as multi-branched fractures. Other scholars proposed composite model to study the pressure transient characteristics in fractured-vuggy reservoirs within which the pore networks in acidized region are equivalent to homogeneous medium of Euclidean geometry[35,36]. Generally speaking, due to the interaction and connectivity between the wormholes and the matrix-fracture-vug systems, the acidized region with multiple property scales and non-Euclidean pore network[37] was highlighted by the characteristics of heterogeneity, which increases the research complexity. As a result, developing an applicable analytical model solving this challenge is more desirable.

It is reported that reservoir heterogeneity provides uncertainty in fluid flow through porous media, so, in order to provide a more realistic representation of flow in acidized fractured-vuggy reservoirs, the heterogeneity of the pore networks in acidized region should be taking into consideration. However, previous models neglect the heterogeneity in acidized region. In this paper, we derive an analytical solution by considering the heterogeneity of the pore networks in acidized region to analyze flow characteristic in fractured-vuggy carbonate reservoirs. In this presented model, pore networks properties (porosity and permeability) in acidized region recognized as an inner zone is modeled based on the fractal geometry theory. The schematic diagram of the physical model is shown in Fig. 1. Compared with the previous models, our model takes the heterogeneity of the pore networks in acidized region into account. The derived model can help to reduce the uncertainty in flow through porous media and obtain data with high accuracy.

Figure 1

Schematic diagram of the physical model.

As the outline of this work, the theoretical models and the corresponding analytical solution in Laplace space are provided in the following section. Then, the results will be analyzed, followed by the discussions of the derived model. Finally, the conclusions are presented.

Methods

Physical concept

As observed in the carbonate formation of the Tahe and Tarim Oilfield in western China, fractured-vuggy reservoirs typically consist of matrix, fracture and vug systems[5]. It is reported that acidizing or acid fracturing stimulation technology has become an important stimulation measure, which is essential to improve or recover productivity for fractured-vuggy carbonate reservoirs[13,38,39]. During the process of acidizing, wormhole dissolution patterns will be affected by a series of complicated physical chemical transportations and reactions[40,41]. As shown in Fig. 2 (i.e. the schematic of the fractured-vuggy model), the physical model was conceptualized as a composite model, which could be divided into the acidized region (interior region) and outer region. In the acidized region, the complicated wormholes are usually connected with matrix-vug-fracture system which is treated as a fractal medium. And the outer region is treated as a common triple medium[29,31,42]. Using the conceptual model, the new composite model is derived under the following assumptions:

Figure 2

Schematic of the fractured-vuggy physical model.

As fractal geometry theory is a natural candidate for the representation of heterogeneity of the acidized region[43-45], fractal geometry theory is applied in this article to account for the pore network in the acidized region. However, the outer region is assumed as a common triple medium, in which matrix and vug are regarded as the main sources of hydrocarbon, and fracture is treated as the flow path[35,46]. The formation effective thickness is constant.

During the flow, the property parameters, such as permeability, porosity and fluid viscosity, are constant. Gravity is neglected. Moreover, the isothermal and Darcy flow process is assumed.

Mathematical model

As depicted in Fig. 3, the composite model can be divided into the acidized region (interior region) and outer region. Based on fractal geometry theory, the permeability and porosity in the interior region can be expressed as[47-49]where k1 is the permeability at the edge of the wellbore, φ1 is porosity at the edge of the wellbore, rw is wellbore radius, df is pore spaces fractal dimension, θ is the conductivity index[49,50].

Figure 3

Comparison of PwD and dPwD obtained from the proposed model and numerical simulation: (a) case 1; (b) case 2; (c) case 3; (d) case 4.

Based on Eq. (1), the corresponding governing equation in interior region isin which p1 is the pressure in interior region, μ is the fluid viscosity, Ct is the rock compressibility in interior region, rc is the outer boundary of the acidized region, and t is the time.

In outer region, the corresponding governing equation can be given bywhere re is the outer boundary of the reservoir, the subscript F, M and v represent the fracture system, matrix system and vug system. kτ (τ = F, M and v) is the permeability in outer region. pτ (τ = F, M and v) is the pressure in outer region. φM is matrix porosity, φF is fracture porosity and φv is vug porosity. CM is matrix compressibility, CF is fracture compressibility and Cv is vug compressibility. αFM and αFv are shape factors for fracture-vug and fracture-matrix, respectively. More details about the calculation equations for physical parameters in triple medium, such as kF, kv, φF, φv, αFM and αFv, can be found in Appendix A.

At the initial state, the following pressure equation can be written aswhere p0 is the initial pressure.

Three types of external boundary conditions (constant pressure boundary, infinite boundary or closed boundary) can be expressed respectively as[5]

The inner boundary condition can be expressed aswhere h is effective thickness, B is fluid volume factor, and q is the flow rate at wellhead, which is positive for producers and negative for injectors.

The interface continuity conditions are governed by

By using dimensionless variables listed in Table 1, Eqs (2–7) can be rewritten as following:

Table 1

Summary of the variables in the model.

Defined parameters	Expressions
Dimensionless time, tD	t_D=k_Ft/[(ϕ_MC_M+ϕ_FC_F+ϕ_vC_v)μrw2]
Dimensionless pressure, pD	p_D=2πk_Fh(p_0−p)/(qBμ)
Dimensionless radius, rD	r_D=r/r_w
Diffusivity ratio between the two regions, η1D	η_1D=C_t1ϕ_1/(ϕ_MC_M+ϕ_FC_F+ϕ_vC_v)
Fracture storage capacitance coefficient, ηFD	η_FD=C_Fϕ_F/(ϕ_MC_M+ϕ_FC_F+ϕ_vC_v)
Matrix storage capacitance coefficient, ηMD	η_MD=C_Mϕ_M/(ϕ_MC_M+ϕ_FC_F+ϕ_vC_v)
Vug storage capacitance coefficient, ηvD	η_vD=C_vϕ_v/(ϕ_MC_M+ϕ_FC_F+ϕ_vC_v)
Permeability ratio, M1	M_1=k_1/k_F
Fracture-matrix inter-porosity coefficient, λFM	λ_FM=α_FMrw2k_M/k_F
Fracture-vug inter-porosity coefficient, λFv	λ_Fv=α_Fvrw2k_v/k_F

Based on the derived process stated in Appendix B, the analytical solution for Eq. (8) in Laplace space can be obtained. The wellbore pressure in Laplace space considering the wellbore storage and the skin effects iswhere S is the skin factor, CD is the dimensionless wellbore storage coefficient. Then, with the Stehfest inversion algorithm[51], the wellbore pressure in real space can be obtained.

Results

Validation of the proposed model

To validate the derived analytical solution, the dimensionless pressure and derivative (PwD and dPwD) of the new derived model are compared with numerical solutions. In the calculation, the outer boundary condition is set as constant pressure, and the function g(z) in Eq. (11) is set to be 1 (the outer region of the model can be simplified as the homogeneous porous media). Table 2 summarizes the data used for the comparison. Result (Fig. 3) illustrates that calculations from the proposed model agree well with those from numerical simulation.

Table 2

Synthetic data used for validation.

Cases	Parameters
	r eD	r cD	C D	S	d f	θ	η 1D	M 1
1	7460	260	100	0	2	0	1	1
2	7460	260	100	0	2	0	1	0.05
3	7460	260	100	0	1.75	0	1	0.05
4	7460	260	100	0	1.75	0.25	3	10

Comparison with former models

When parameters df, θ and g(z) in our derived model are assigned as 2, 0 and 1, respectively, the acidized region (inner region) can be considered homogeneous, and the reservoir can be regarded as two different but uniform regions. In this case, our derived model can be validated with radial composite model. Figure 4a presents the comparison of the results (PwD versus tD) from the derived model and that from fractal radial composite model proposed by Razminia et al.[52] When parameters df and θ in fractal radial composite model are set to be 2 and 0, respectively, fractal radial composite model can be simplified as the conventional radial composite model. In the calculation, the outer boundary condition is set as infinite boundary, other parameters applied in the model assigned were: CD = 200, S = 2 and rcD = 500. Results (Fig. 4a) suggest that the derived model is validated. When parameters df, θ, M1, η1D, and g(z) in our derived model are assigned as 2, 0, 1, 1 and 1, respectively, the reservoir can be regarded as a homogeneous porous media. In this case, the derived model can be compared with the classic radial flow model derived by Agarwal et al.[53] Fig. 4b presents the comparison of the results (PwD versus tD) from the derived model and that from the classic analytical solution[53]. In the calculation, the outer boundary condition is set as infinite boundary, other parameters applied in the model assigned were: CD = 200 and S = 2. Results (Fig. 4b) also suggest that the derived model is validated. Figure 4a,b show the former classic models are special cases of the new derived model, and the new derived analytical model can be used to analyze more complex flow behaviors.

Figure 4

Comparison of PwD obtained from the proposed model and other models: (a) for radial composite model; (b) for radial model. Results from former models are shown with points, and calculations from the derived model are solid lines.

Flow characteristic curves

As stated in the literatures[5,33], the flow characteristic curves can be used for recognizing flow behavior in the reservoir. Taking the closed outer boundary condition as example, the flow characteristic curve will be discussed in detail. Figure 5a presents the log-log curves (i.e., tD versus PwD and dPwD) for closed outer boundary condition. In the calculation, fractal dimension of the pore spaces df is 2, the conductivity index θ is 0, the dimensionless parameters reD and rcD are 30000 and 2000, respectively, skin factor S assigned is 0, CD assigned is 200, diffusivity ratio between the two regions η1D is 0.127, fracture storage capacitance coefficient ηFD is 0.022, matrix storage capacitance coefficient ηMD is 0.106, vug storage capacitance coefficient ηvD is 0.872, permeability ratio M1 is 1.251, λFM assigned is 2.01 × 10−7, and λFv assigned is 9.70 × 10−3. With the same parameters, Fig. 5b compares the differences of three well testing curves to study the effect of outer boundary on flow characteristics. Results (Fig. 5a,b) illustrate that the flow characteristic curves can be characterized by six flow stages:

Figure 5

Flow characteristic curves: (a) pressure & pressure derivative curves for closed boundary condition; (b) pressure derivative curves for three different outer boundary conditions.

The first stage is well storage stage, in which the fluid flow is dominated by the parameters S and CD.

The second stage is radial flow in the inner region. The dPwD curve shows zero slope horizontal line. During this period, the fluid in the interior region flows into the wellbore. The flow characteristic is mainly controlled by the parameters rcD, pore spaces fractal dimension and the conductivity index. The larger the range of the inner region is, the longer the stage 2 will be.

The third stage is fracture-vug inter-porosity flow stage, in which the dPwD curve shows the first concave. The concave shape represents the supplement between the fracture system and vug system. In this stage, the pressure wave has propagated to the coupling interface, and the flow characteristic in this period is mainly controlled by fracture-vug inter-porosity coefficient.

The fourth stage is transition flow stage, which presents a straight line with the slope 0.5 on the log-log curve of tD versus dPwD.

The fifth stage is fracture-matrix inter-porosity flow stage, in which the dPwD curve shows the second concave shape. The second concave shape segment reflects supplement from matrix to fracture. During this stage, flow characteristic is mainly controlled by fracture-matrix inter-porosity coefficient.

The sixth stage is external boundary response stage, in which flow characteristic in the reservoir is mainly controlled by the parameter reD.

Sensitivity analysis

In this part, we will study the influences of parameters on well testing curves in detail. The outer boundary is set as closed pressure, the parameter values is listed in Table 3. The presented model demonstrates that df, θ, rc, LM, LF and Lv are the basic parameters, which can determine the flow characteristic in acidized fractured-vuggy reservoirs. The basic parameters and their ranges used for sensitivity analysis are summarized in Table 4.

Table 3

Parameter values.

Parameter	Parameter value	Parameter	Parameter value
φ M	0.1	r eD	30000
rw/m	0.1	r cD	2000
LM/m	2	η 1D	0.127
LF/m	8 × 10−4	η FD	0.022
Lv/m	1.08	η MD	0.106
k1/m2	1.0 × 10−10	η vD	0.872
KM/m2	1.6 × 10−14	λ FM	2.01 × 10−7
d f	2	λ Fv	9.66 × 10−3
θ	0	M 1	1.56

Table 4

Ranges of parameters for sensitivity analysis.

Parameter	Range of the parameter	Parameter	Range of the parameter
d f	1.65 ≤ df ≤ 2	LM/m	1.5 ≤ LM ≤ 2.5
θ	−0.25 ≤ θ ≤ 0.25	LF/m	3 × 10−4 ≤ LF ≤ 8 × 10−4
rc/m	100 ≤ rc ≤ 300	Lv/m	0.8 ≤ Lv ≤ 1.2

Figure 6 shows flow characteristic curves with different inner region parameters (i.e., inner region fractal dimension, inner region connectivity index and inner region radius). Figure 6a presents the effect of parameter df on the flow characteristic curves (for 1.65 ≤ df ≤ 2). As shown in Fig. 6a, PwD and dPwD in stages 1, 2 and 3 decrease with the increased inner region fractal dimension df. This can be interpreted as that the larger inner region fractal dimension df implies inner region permeability increasing, leading to the decrease of flow resistance and the dimensionless pressure. Figure 6b demonstrates the effect of inner region connectivity index θ on the flow characteristic curves (for −0.25 ≤ θ ≤ 0.25). Figure 6b demonstrates PwD and dPwD in stages 1, 2 and 3 increase with the increased inner region connectivity index θ. Figure 6c illustrates the effect of inner region radius on the flow characteristic curves (for 100 m ≤ rc ≤ 300 m). Results (Fig. 6c) show PwD and dPwD in stages 3 to 5 decrease with the increased inner region radius. The main reason for this is the permeability at the edge of the wellbore k1 is larger than fracture system permeability kF (M1 > 1). Figure 6c also shows a longer period of radial flow (stage 2) corresponds to a larger inner region radius, which is expected. Figure 6 demonstrates that forming larger acidized region can help improve well production and reduce fluid pressure loss.

Figure 6

Effects of parameters of inner region on the well testing curves: (a) for inner region fractal dimension; (b) for inner region connectivity index; (c) for parameter rc.

Figure 7 presents the influences of the outer region parameters, such as the space interval of fracture system LF, matrix block size LM and radius of the vug Lv, on the transient pressure behavior. Figure 7a presents the effect of the space interval of fracture system LF on the flow characteristic curves (for 3 × 10−4 m ≤ LF ≤ 8 × 10−4 m). As shown in Fig. 7a, a larger LF corresponds to a larger pressure depletion. The main reason is that a greater parameter LF value results in a smaller parameter M1 and the smaller parameter η1D. We can see that the effect of parameter LF on the well testing curves is concentrated in stages 1 to 4. Figure 7b demonstrates the effect of matrix block size in outer region on the flow characteristic curves (for 1.5 m ≤ LM ≤ 2.5 m). We can see that PwD in stages 1 to 4 decreases with the increasing matrix block size LM. As a larger LM corresponding to a larger M1, which leads to a smaller pressure depletion in stages 1 and 2. Furthermore, parameter ηvD decreases as parameter LM increases, indicating the supplement of the vug system to fracture system decreases will the increasing parameter LM. As a result, the dimensionless pressure derivation in stage 3 increases with the increasing LM. Figure 7c shows the effect of parameter Lv in outer region on the flow characteristic curves (for 0.5 m ≤ Lv ≤ 1.2 m). It can be seen that the effects of Lv on PwD and dPwD curves are concentrated in stage 3, where a larger Lv will lead to a smaller pressure depletion. The main reason is that, a larger Lv implies larger energy in the vugs, leading to increased supply from the vug system to the fracture system, which makes the pressure derivative values smaller.

Figure 7

Effect of basic parameters in outer region on the dimensionless wellbore pressure: (a) for the space interval of fracture in outer region; (b) for the matrix block size in outer region; (c) for the radius of the vug in outer region.

Case example

In this case, pressure transient data from an acidized well in Tahe oilfield (Tarim basin, China), is selected to interpret. The actual test data of the acidized well are cited from Wang et al.[33]. In the calculation, the parameters measured are PwD and dPwD and the corresponding tD, and the fitting parameters mainly concludes: df; θ; rcD; η1D; ηFD; ηMD; ηvD; M1; λFv; λFM; CD and S. Figure 8 presents the fitting results and the actual test data, and the corresponding results by well test interpretation are as follows: df = 1.95; θ = 0.15; rcD = 100; η1D = 9 × 10−6; ηFD = 9 × 10−5; ηMD = 4.8 × 10−2; ηvD = 0.952; M1 = 8; λFv = 4 × 10−5; λFM = 6.9 × 10−3; CD = 0.0017; S = 0. Figure 8 illustrates that the results agree with each other, which verifies the presented model.

Figure 8

Comparison of the calculated results with actual test data.

Discussions

Advantages and limitations of the derived model

The analytical model lays theoretical foundations for characterizing the flow characteristic curves and dynamic pressure of acidized well in acidized carbonate porous media. With this new analytical solution, it can help to reduce the uncertainty in flow through porous media and obtain data with high accuracy. In addition, it can be applied to field case study and estimate more accurate reservoir parameters reflecting flow characteristic using inverse modeling. However, it should be noted that the derived model is limited to single-phase flow. Multi-phase fluid flow in fractured-vuggy carbonate porous medium is an interesting and challenging topic, and the corresponding work is still under processing.

Conclusions

Based on fractal theory and triple medium concept model, a novel analytical solution to a fractal composite model is derived to study the flow characteristic in acidized carbonate reservoirs. Predicted results deduced from the derived analytical solution agree well with numerical simulation results. For the closed outer boundary condition, the influence of basic parameters, such as df, θ, rc, LM, LF and Lv, on flow characteristic curves are studied in detail. Followings are the main conclusions:

The calculated results (both PwD and d PwD) based on the derived analytical solution match the actual test data very well, which illustrates this proposed model can be applied to field case study.

As the inner region permeability is larger than fracture system permeability, forming larger acidized region can help improve well production and reduce fluid pressure loss.

The new analytical solution to the fractal composite model, considering the heterogeneity of the pore networks in acidized region, is useful for flow characteristic analysis of acidized well in acidized carbonate porous media. With this new analytical solution, it can help to reduce the uncertainty in flow through porous media and obtain data with high accuracy. In addition, it can be applied to estimate more accurate reservoir parameters reflecting flow characteristic using inverse modeling.

Appendix