Pressure transient analysis for a fractured well in a stress-sensitive tight multi-medium oil reservoir

Abstract

Tight multi-medium oil reservoirs are the main source of hydrocarbon resources around the world. Acid fracturing is the most effective technology to improve productivity in such reservoirs. As carbonates are primarily composed of dolomite and calcite, which are easily dissolved by hydrochloric acid, high-permeability region will be formed near the well along with the main artificial fracture when acid fracturing is implemented in tight multi-medium oil reservoirs. In this study, a comprehensive composite linear flow model was developed to simulate the transient pressure behavior of an acid fracturing vertical well in a naturally fractured vuggy carbonate reservoir. By utilizing Pedrosa’s substitution, perturbation, Laplace transformation and Stehfest numerical inversion technology, the pressure behavior results were obtained in real time domain. Furthermore, the result of this model was validated by comparing with those of previous literature. Additionally, the influences of some prevailing parameters on the type curves were analyzed. Moreover, the proposed model was applied to an acid fracturing well to evaluate the effectiveness of acid fracturing measures, to demonstrate the practicability of the proposed model.

Graphical abstract

Keywords

tight multi-medium oil reservoir / acid fracturing / stress-sensitive permeability / composite linear flow

1 Introduction

Tight multi-medium oil reservoirs are widely distributed around the world, and numerous studies have been conducted on this type of reservoirs (Kossack and Gurpinar, 2001; Kang et al., 2006; Wang et al., 2018a; Xing et al., 2018; Xu et al., 2019). In the last few decades, many scholars noticed the distinctive pore structure and fluid flow mechanism in carbonate reservoirs, where matrix, natural fractures and vugs coexist in naturally fractured vuggy reservoirs. This renders it complicated and intractable to characterize reservoirs accurately. Generally, matrix and vuggy pores act as storage spaces for hydrocarbon fluids, whereas fractures are usually considered as a pathway for fluid flow. Furthermore, vuggy pores can be subdivided as connected and disconnected with natural fractures. This means that fluid in the vug can flow into fracture directly or indirectly via the bridge of matrix.

Due to complex pore types, it is challenging to model fluid flow through tight multi-medium oil reservoirs. After decades of researches, scholars have put forward some effective methods to tackle this problem. Abdassah and Ershaghi (1986) first proposed triple-porosity/single-permeability model. In their model, and unsteady-state interporosity flow between fracture and other systems was considered in their model. Later, Liu et al. (2003) proposed a tri-continuum medium concept considering pseudo-steady interporosity flow. Wu et al. (2006 and 2007) developed an analytical method for transient flow analysis in tight multi-medium oil reservoirs, taking the flow between vug and matrix into consideration. Camacho-Velázquez et al. (2002) established a triple-porosity/dual-permeability model to consider primary flow not only through fracture system to wellbore, but also vugs system to wellbore. Subsequently, Fuentes-Cruz et al. (2004) extended the model to partially penetrated well. Yao et al. (2010) established the discrete fracture-vug network model and provided to describe fluid flow in the fractured-vuggy porous media. Wu et al (2019) applied the discrete fracture-vug model proposed by Yao et al. (2010) to model macrovugs, while microfractures and microvugs are modeled with the triple-continuum concept. Guo et al. (2012) established a test analysis model of a horizontal well, and triple porosity and dual permeability flow behavior were analyzed. Their results showed that type curves were dominated by external boundary conditions as well as the permeability ration of fracture system to the sum of fracture and matrix systems.

In order to reduce flow resistance meanwhile extract hydrocarbon resources in carbonate formations effectively, acid fracturing has been the predominant technology (Abass et al., 2006). Thus, it is significant and attractive to understand the flow behavior of acid fracturing wells in tight multi-medium oil reservoirs. Generally speaking, carbonate minerals are composed of dolomite and calcite which are easy to be dissolved by hydrochloric acid. Therefore, acid fracturing can not only act as “hydraulic fracturing”, but also act as acidification to create different types of wormholes (Fredd, 2000; Dora, 2008; Liu et al., 2013). As for the pressure dynamics analysis for acid fracturing wells in tight multi-medium oil reservoirs, great efforts have been made in the last several years. Wang et al. (2014) developed a theoretical wormhole seepage model for the first time. In this model, wormholes are simplified as multi-branched fractures with an infinite conductivity. Later, a semi-analytical model that derives from the previous one to investigate the pressure behavior for finite conductivity multi-branched fractures in fractured-vuggy reservoirs in detail was also suggested by Wang et al. (2018b). Wang and Yi (2017) studied the flow behavior of a well with a finite-conductivity fracture in a tight multi-medium oil reservoir. In their study, the vugs were are conceptualized as spherical shapes and seven flow regimes were are observed. Recently, Wang et al. (2018b) investigated the flow behavior of acid fracturing wells in a composite fractured-vuggy carbonate reservoirs. Lei et al. (2018) presented an analytical solution considering the heterogeneity of the pore networks in acidized region by the utility of the fractal geometry theory.

It is reported that some tight multi-medium oil reservoirs may exhibit strong stress-sensitive characteristic, which has a significant effect on transient pressure behavior (Zhang et al., 2017; Yang et al., 2017). Yet, the aforementioned researches haven’t taken this factor into consideration. Moreover, all of their wok were based on the assumption that the fracture is completely penetrated in the vertical direction. Unfortunately, to be our best knowledge, fracture may be partially penetrated in actual acid fracturing implementations. With regard to partially penetrating fractured well, many researchers (Raghavan et al., 1978; Rodriguez et al., 1984; Igbokoyi and Tiab, 2008; Zhang et al., 2015; Yuan et al., 2018) developed point/slab source functions and numerical methods to deal with it respectively. However, numerical methods are complex and time-consuming. Compared with numerical simulation, composite linear flow model can solve the problem of single vertically fractured well or multi-stage fractured horizontal well conveniently. Simultaneously, it also avoids the time-consuming process in numerical simulation and the complexity of Green function methods (Brown et al., 2009; Stalgorova and Matter, 2012; Tao et al., 2018; Zeng et al., 2018; Zeng et al., 2017).

In this study, we investigate an acid fracturing vertical well in a rectangular fractured-vuggy carbonate reservoir, and the acidized region is also rectangular. Besides, the artificial fracture can be either fully or partially penetrated. The Pedrosa’s perturbation (Pedrosa, 1986) is utilized to linearize the non-linear equations caused by stress-sensitive permeability. The physical model and relevant assumptions would be elaborated in Section 2. In Section 3, the mathematical model is established and the corresponding solution is given by Laplace transformation and Stehfest inversion technology (Stehfest, 1970). The model verification and parameters sensitivity analysis are presented in Section 4. Finally, some remarkable conclusions are drawn in Section 5.

2 Physical model

As shown in Fig. 1, a single fractured vertical well is located in the center of a tight multi-medium oil reservoir. The reservoir is modeled as triple-porosity medium while the main artificial fracture is modeled as single-porosity medium. As for the triple-porosity medium, natural fractures act as the main pathways, while vug and matrix provide the storage space. The inter-porosity flow is assumed to be in pseudo-steady state, and the process is shown in Fig. 2. The main artificial fracture formed by acid fracturing is symmetric and it partially penetrates the formation in both vertical and horizontal directions. There is an acidized region with lots of wormholes near the main artificial fracture. Slightly compressible fluid with constant viscosity is produced at a fixed rate through the well. For the sake of simplification, some of other tenable assumptions are elaborated as following:

1) The reservoir is horizontal with uniform thickness of h, and the reservoir pressure is pi at the begining time of production and the gadient is uniformly distributed;

2) The upper and bottom boundaries of reservoir are impermeable, and external boundary is also closed;

3) Main artificial fracture and natural fractures in regions 1 to 6 possess stress-sensitive permeability and it can be described as the following formula:

(1)k_f=k_{fi}{e}^{-\alpha (p_i-p_f)},

where k_f denotes fracture permeability at current pressure; k_fidenotes fracture permeability at initial pressure; α refers to permeability modulus; p_iand p_frefer to initial pressure and fracture pressure.

4) The continuity conditions of pressure and flux at the interfaces are used to connect the adjacent flow regions.

5) Isothermal and Darcy flow process is assumed while the gravity and capillary pressure are ignored.

3 Mathematical model and solution

Based on the aforementioned assumption that the hydraulic fracture is symmetrical both in horizontal and vertical direction, one-quarter of the rectangular carbonate reservoir is implemented to simplify the problem. The system is divided into seven regions as shown in Fig. 3, and the shape of each region presents approximately rectangular.

3.1 Flow in regions without acidification (Region 2+ Region 3+ Region 4+ Region 5+ Region 6)

Region 6

In Region 6, considering the effect of stress sensitivity, the governing equations in dimensionless form can be written as:

(2)e^{-{\alpha }_D{p}_{\mathrm{6}fD}}\left[{\displaystyle \frac{{\partial }^{\mathrm{2}}{p}_{\mathrm{6}fD}}{\partial z_D^{\mathrm{2}}}}-{\alpha }_D{\left({\displaystyle \frac{\partial p_{\mathrm{6}fD}}{\partial z_D}}\right)}^{\mathrm{2}}\right]={\displaystyle \frac{\mathrm{1}}{{\eta }_{\mathrm{6}D}}}\left[{\omega }_{\mathrm{6}f}{\displaystyle \frac{\partial p_{\mathrm{6}fD}}{\partial t_D}}+{\omega }_{\mathrm{6}m}{\displaystyle \frac{\partial p_{\mathrm{6}mD}}{\partial t_D}}+{\omega }_{\mathrm{6}v}{\displaystyle \frac{\partial p_{\mathrm{6}vD}}{\partial t_D}}\right],

(3){\omega }_{\mathrm{6}m}{\displaystyle \frac{\mathrm{1}}{{\eta }_{\mathrm{6}D}}}{\displaystyle \frac{\partial p_{\mathrm{6}mD}}{\partial t_D}}+{\lambda }_{\mathrm{6}fm}(p_{\mathrm{6}mD}-p_{\mathrm{6}fD})-{\lambda }_{\mathrm{6}vm}(p_{\mathrm{6}vD}-p_{\mathrm{6}mD})=\mathrm{0},

(4){\omega }_{\mathrm{6}v}{\displaystyle \frac{\mathrm{1}}{{\eta }_{\mathrm{6}D}}}{\displaystyle \frac{\partial p_{\mathrm{6}mD}}{\partial t_D}}+{\lambda }_{\mathrm{6}fv}(p_{\mathrm{6}vD}-p_{\mathrm{6}fD})-{\lambda }_{\mathrm{6}vm}(p_{\mathrm{6}vD}-p_{\mathrm{6}mD})=\mathrm{0},

where, α is permeability modulus, MPa^-1; ω is storativity ratio, dimensionless;ηis diffusivity coefficient, MPa^-1; λ is interporosity flow coefficient, dimensionless; fis natural fracture property; mis matrix property; v is vug property; p is reservoir pressure, MPa; t is time;e is Euler number; D is dimensionless.

Dimensionless parameters in Eqs. (2)‒(4) are defined in Appendix A.

The initial condition is:

(5)p_{\mathrm{6}fD}=p_{\mathrm{6}mD}=p_{\mathrm{6}vD}=\mathrm{0.}

The outer boundary condition (no-flow) is:

(6){{\displaystyle \frac{\partial p_{\mathrm{6}D}}{\partial z_D}}|}_{z_D=z_{\mathrm{2}D}}=\mathrm{0},

where z is distance

Based on pressure continuity, the inner boundary condition is given as:

(7){p_{\mathrm{6}D}|}_{z_D=z_{\mathrm{1}D}}={p_{\mathrm{2}D}|}_{z_D=z_{\mathrm{1}D}}={p_{\mathrm{4}D}|}_{z_D=z_{\mathrm{1}D}}.

With Pedrosa-substitution (Pedrosa, 1986), perturbation and Laplace transformation methods (detailed derivations are in Appendix B), we can obtain:

(8)\begin{array}{l}{{\displaystyle \frac{\partial {\tilde{\xi }}_{\mathrm{6}D}}{\partial z_D}}|}_{z_{\mathrm{1}D}}={\tilde{\xi }}_{\mathrm{2}D}(z_{\mathrm{1}D})\sqrt{{\displaystyle \frac{s{g}_{\mathrm{6}}(s)}{{\eta }_{\mathrm{6}D}}}}\tan h\left[\sqrt{{\displaystyle \frac{s{g}_{\mathrm{6}}(s)}{{\eta }_{\mathrm{6}D}}}}(z_{\mathrm{1}D}-z_{\mathrm{2}D})\right]\\ \mathrm{}={\tilde{\xi }}_{\mathrm{4}D}(z_{\mathrm{1}D})\sqrt{{\displaystyle \frac{s{g}_{\mathrm{6}}(s)}{{\eta }_{\mathrm{6}D}}}}\tan h\left[\sqrt{{\displaystyle \frac{s{g}_{\mathrm{6}}(s)}{{\eta }_{\mathrm{6}D}}}}(z_{\mathrm{1}D}-z_{\mathrm{2}D})\right].\end{array}

where s is dimensionless time variable in Laplace domain, dimensionless; ξ is a variable after Pedrosa’s substitution ;~ is laplace transform; h is formation height, m.

Region 5

Analogously, the governing equations of Region 5 in dimensionless form can be written as:

(9)e^{-{\alpha }_D{p}_{\mathrm{5}fD}}\left[{\displaystyle \frac{{\partial }^{\mathrm{2}}{p}_{\mathrm{5}fD}}{\partial z_D^{\mathrm{2}}}}-{\alpha }_D{\left({\displaystyle \frac{\partial p_{\mathrm{5}fD}}{\partial z_D}}\right)}^{\mathrm{2}}\right]={\displaystyle \frac{\mathrm{1}}{{\eta }_{\mathrm{5}D}}}\left[{\omega }_{\mathrm{5}f}{\displaystyle \frac{\partial p_{\mathrm{5}fD}}{\partial t_D}}+{\omega }_{\mathrm{5}m}{\displaystyle \frac{\partial p_{\mathrm{5}mD}}{\partial t_D}}+{\omega }_{\mathrm{5}v}{\displaystyle \frac{\partial p_{\mathrm{5}vD}}{\partial t_D}}\right],

(10){\omega }_{\mathrm{5}m}{\displaystyle \frac{\mathrm{1}}{{\eta }_{\mathrm{5}D}}}{\displaystyle \frac{\partial p_{\mathrm{5}mD}}{\partial t_D}}+{\lambda }_{\mathrm{5}fm}(p_{\mathrm{5}mD}-p_{\mathrm{5}fD})-{\lambda }_{\mathrm{5}vm}(p_{\mathrm{5}vD}-p_{\mathrm{5}mD})=\mathrm{0},

(11){\omega }_{\mathrm{5}v}{\displaystyle \frac{\mathrm{1}}{{\eta }_{\mathrm{5}D}}}{\displaystyle \frac{\partial p_{\mathrm{5}mD}}{\partial t_D}}+{\lambda }_{\mathrm{5}fv}(p_{\mathrm{5}vD}-p_{\mathrm{5}fD})-{\lambda }_{\mathrm{5}vm}(p_{\mathrm{5}vD}-p_{\mathrm{5}mD})=\mathrm{0.}

The initial condition is:

(12)p_{\mathrm{5}fD}=p_{\mathrm{5}mD}=p_{\mathrm{5}vD}=\mathrm{0.}

The outer boundary condition (no-flow) is:

(13){{\displaystyle \frac{\partial p_{\mathrm{5}fD}}{\partial z_D}}|}_{z_D=z_{\mathrm{2}D}}=\mathrm{0.}

Based on pressure continuity, the inner boundary condition is given as:

(14){p_{\mathrm{5}fD}|}_{z_D=z_{\mathrm{1}D}}={p_{\mathrm{1}fD}|}_{z_D=z_{\mathrm{1}D}}={p_{\mathrm{3}fD}|}_{z_D=z_{\mathrm{1}D}}.

After eliminating the nonlinearity by applying the Pedrosa substitution, the solution of Region 5 can be obtained as following:

(15)\begin{array}{l}{{\displaystyle \frac{\partial {\tilde{\xi }}_{\mathrm{5}D}}{\partial z_D}}|}_{z_{\mathrm{1}D}}={\tilde{\xi }}_{\mathrm{1}D}(z_{\mathrm{1}D})\sqrt{{\displaystyle \frac{s{g}_{\mathrm{5}}(s)}{{\eta }_{\mathrm{5}D}}}}\tan h\left[\sqrt{{\displaystyle \frac{s{g}_{\mathrm{5}}(s)}{{\eta }_{\mathrm{5}D}}}}(z_{\mathrm{1}D}-z_{\mathrm{2}D})\right]\\ \mathrm{}={\tilde{\xi }}_{\mathrm{3}D}(z_{\mathrm{1}D})\sqrt{{\displaystyle \frac{s{g}_{\mathrm{5}}(s)}{{\eta }_{\mathrm{5}D}}}}\tan h\left[\sqrt{{\displaystyle \frac{s{g}_{\mathrm{5}}(s)}{{\eta }_{\mathrm{5}D}}}}(z_{\mathrm{1}D}-z_{\mathrm{2}D})\right].\end{array}

Region 4

In Region 4, the flow is in x and z directions, and the governing equations considering the effect of stress sensitivity in dimensionless form can be written as:

(16)e^{-{\alpha }_D{p}_{\mathrm{4}fD}}\left[{\displaystyle \frac{{\partial }^{\mathrm{2}}{p}_{\mathrm{4}fD}}{\partial x_D^{\mathrm{2}}}}+{\displaystyle \frac{{\partial }^{\mathrm{2}}{p}_{\mathrm{4}fD}}{\partial z_D^{\mathrm{2}}}}-{\alpha }_D{\left({\displaystyle \frac{\partial p_{\mathrm{4}fD}}{\partial x_D}}\right)}^{\mathrm{2}}-{\alpha }_D{\left({\displaystyle \frac{\partial p_{\mathrm{4}fD}}{\partial z_D}}\right)}^{\mathrm{2}}\right]={\displaystyle \frac{\mathrm{1}}{{\eta }_{\mathrm{4}D}}}\left[{\omega }_{\mathrm{4}f}{\displaystyle \frac{\partial p_{\mathrm{4}fD}}{\partial t_D}}+{\omega }_{\mathrm{4}m}{\displaystyle \frac{\partial p_{\mathrm{4}mD}}{\partial t_D}}+{\omega }_{\mathrm{4}v}{\displaystyle \frac{\partial p_{\mathrm{4}vD}}{\partial t_D}}\right],

(17){\omega }_{\mathrm{4}m}{\displaystyle \frac{\mathrm{1}}{{\eta }_{\mathrm{4}D}}}{\displaystyle \frac{\partial p_{\mathrm{4}mD}}{\partial t_D}}+{\lambda }_{\mathrm{4}fm}(p_{\mathrm{4}mD}-p_{\mathrm{4}fD})-{\lambda }_{\mathrm{4}vm}(p_{\mathrm{4}vD}-p_{\mathrm{4}mD})=\mathrm{0},

(18){\omega }_{\mathrm{4}v}{\displaystyle \frac{\mathrm{1}}{{\eta }_{\mathrm{4}D}}}{\displaystyle \frac{\partial p_{\mathrm{4}mD}}{\partial t_D}}+{\lambda }_{\mathrm{4}fv}(p_{\mathrm{4}vD}-p_{\mathrm{4}fD})-{\lambda }_{\mathrm{4}vm}(p_{\mathrm{4}vD}-p_{\mathrm{4}mD})=\mathrm{0.}

The no-flow outer boundary in x direction is:

(19){{\displaystyle \frac{\partial p_{\mathrm{4}fD}}{\partial x_D}}|}_{x_{\mathrm{2}D}}=\mathrm{0.}

The pressure continuity condition between Region 4 and Region 2 on the interface is:

(20){p_{\mathrm{4}fD}|}_{x_D=x_{\mathrm{1}D}}={p_{\mathrm{2}fD}|}_{x_D=x_{\mathrm{1}D}}.

Based on the definite conditions, the solution of Region 4 can be written as:

(21){{\displaystyle \frac{\partial {\tilde{\xi }}_{\mathrm{4}D}}{\partial x_D}}|}_{x_{\mathrm{1}D}}={\tilde{\xi }}_{\mathrm{2}D}(x_{\mathrm{1}D})\sqrt{{\alpha }_{\mathrm{4}}}\tan h[\sqrt{{\alpha }_{\mathrm{4}}}(x_{\mathrm{1}D}-x_{\mathrm{2}D})].

Region 3

In Region 3, the flow is in x and z directions, and the governing equations considering the effect of stress sensitivity in dimensionless form can be written as:

(22)e^{-{\alpha }_D{p}_{\mathrm{3}fD}}\left[{\displaystyle \frac{{\partial }^{\mathrm{2}}{p}_{\mathrm{3}fD}}{\partial x_D^{\mathrm{2}}}}+{\displaystyle \frac{{\partial }^{\mathrm{2}}{p}_{\mathrm{3}fD}}{\partial z_D^{\mathrm{2}}}}-{\alpha }_D{\left({\displaystyle \frac{\partial p_{\mathrm{3}fD}}{\partial x_D}}\right)}^{\mathrm{2}}-{\alpha }_D{\left({\displaystyle \frac{\partial p_{\mathrm{3}fD}}{\partial z_D}}\right)}^{\mathrm{2}}\right]={\displaystyle \frac{\mathrm{1}}{{\eta }_{\mathrm{3}D}}}\left[{\omega }_{\mathrm{3}f}{\displaystyle \frac{\partial p_{\mathrm{3}fD}}{\partial t_D}}+{\omega }_{\mathrm{3}m}{\displaystyle \frac{\partial p_{\mathrm{3}mD}}{\partial t_D}}+{\omega }_{\mathrm{3}v}{\displaystyle \frac{\partial p_{\mathrm{3}vD}}{\partial t_D}}\right],

(23){\omega }_{\mathrm{3}m}{\displaystyle \frac{\mathrm{1}}{{\eta }_{\mathrm{3}D}}}{\displaystyle \frac{\partial p_{\mathrm{3}mD}}{\partial t_D}}+{\lambda }_{\mathrm{3}fm}(p_{\mathrm{3}mD}-p_{\mathrm{3}fD})-{\lambda }_{\mathrm{3}vm}(p_{\mathrm{3}vD}-p_{\mathrm{3}mD})=\mathrm{0},

(24){\omega }_{\mathrm{3}v}{\displaystyle \frac{\mathrm{1}}{{\eta }_{\mathrm{3}D}}}{\displaystyle \frac{\partial p_{\mathrm{3}mD}}{\partial t_D}}+{\lambda }_{\mathrm{3}fv}(p_{\mathrm{3}vD}-p_{\mathrm{3}fD})-{\lambda }_{\mathrm{3}vm}(p_{\mathrm{3}vD}-p_{\mathrm{3}mD})=\mathrm{0.}

The no-flow outer boundary in x direction is:

(25){{\displaystyle \frac{\partial p_{\mathrm{3}fD}}{\partial x_D}}|}_{x_{\mathrm{2}D}}=\mathrm{0.}

The pressure continuity condition between Region 3 and Region 1 on the interface is:

(26){p_{\mathrm{3}fD}|}_{x_D=x_{\mathrm{1}D}}={p_{\mathrm{1}fD}|}_{x_D=x_{\mathrm{1}D}}.

Based on the definite conditions, the solution of Region 3 can be written as:

(27){{\displaystyle \frac{\partial {\tilde{\xi }}_{\mathrm{3}D}}{\partial x_D}}|}_{x_{\mathrm{1}D}}={\tilde{\xi }}_{\mathrm{1}D}(x_{\mathrm{1}D})\sqrt{{\alpha }_{\mathrm{3}}}\tan h[\sqrt{{\alpha }_{\mathrm{3}}}(x_{\mathrm{1}D}-x_{\mathrm{2}D})].

Region 2

In Region 2, the flow is in x, y and z directions, and the governing equations considering the effect of stress sensitivity in dimensionless form can be written as:

(28)\begin{array}{l}{e}^{-{\alpha }_D{p}_{\mathrm{2}fD}}\left[{\displaystyle \frac{{\partial }^{\mathrm{2}}{p}_{\mathrm{2}fD}}{\partial y_D^{\mathrm{2}}}}+{\displaystyle \frac{{\partial }^{\mathrm{2}}{p}_{\mathrm{2}fD}}{\partial x_D^{\mathrm{2}}}}+{\displaystyle \frac{{\partial }^{\mathrm{2}}{p}_{\mathrm{2}fD}}{\partial z_D^{\mathrm{2}}}}-{\alpha }_D{\left({\displaystyle \frac{\partial p_{\mathrm{2}fD}}{\partial y_D}}\right)}^{\mathrm{2}}-{\alpha }_D{\left({\displaystyle \frac{\partial p_{\mathrm{2}fD}}{\partial x_D}}\right)}^{\mathrm{2}}-{\alpha }_D{\left({\displaystyle \frac{\partial p_{\mathrm{2}fD}}{\partial z_D}}\right)}^{\mathrm{2}}\right]\\ ={\displaystyle \frac{\mathrm{1}}{{\eta }_{\mathrm{2}D}}}\left[{\omega }_{\mathrm{2}f}{\displaystyle \frac{\partial p_{\mathrm{2}fD}}{\partial t_D}}+{\omega }_{\mathrm{2}m}{\displaystyle \frac{\partial p_{\mathrm{2}mD}}{\partial t_D}}+{\omega }_{\mathrm{2}v}{\displaystyle \frac{\partial p_{\mathrm{2}vD}}{\partial t_D}}\right],\end{array}

(29){\omega }_{\mathrm{2}m}{\displaystyle \frac{\mathrm{1}}{{\eta }_{\mathrm{2}D}}}{\displaystyle \frac{\partial p_{\mathrm{2}mD}}{\partial t_D}}+{\lambda }_{\mathrm{2}fm}(p_{\mathrm{2}mD}-p_{\mathrm{2}fD})-{\lambda }_{\mathrm{2}vm}(p_{\mathrm{2}vD}-p_{\mathrm{2}mD})=\mathrm{0},

(30){\omega }_{\mathrm{2}v}{\displaystyle \frac{\mathrm{1}}{{\eta }_{\mathrm{2}D}}}{\displaystyle \frac{\partial p_{\mathrm{2}mD}}{\partial t_D}}+{\lambda }_{\mathrm{2}fv}(p_{\mathrm{2}vD}-p_{\mathrm{2}fD})-{\lambda }_{\mathrm{2}vm}(p_{\mathrm{2}vD}-p_{\mathrm{2}mD})=\mathrm{0.}

The no-flow outer boundary in y direction is:

(31){{\displaystyle \frac{\partial p_{\mathrm{2}fD}}{\partial y_D}}|}_{y_{\mathrm{2}D}}=\mathrm{0.}

The pressure continuity condition between Region 2 and Region 1 on the interface is:

(32)p_{\mathrm{2}fD}(y_{\mathrm{1}D})=p_{\mathrm{1}fD}(y_{\mathrm{1}D}).

(33){{\displaystyle \frac{\partial {\tilde{\xi }}_{\mathrm{2}fD}}{\partial y_D}}|}_{y_{\mathrm{1}D}}={\tilde{\xi }}_{\mathrm{1}fD}(y_{\mathrm{1}D})\sqrt{{\alpha }_{\mathrm{2}}}\tan h[\sqrt{{\alpha }_{\mathrm{2}}}(y_{\mathrm{1}D}-y_{\mathrm{2}D})].

3.2 Flow in acidized region (Region 1)

In Region 1, the flow is also in x, y and z directions, and the governing equations considering the effect of stress sensitivity in dimensionless form can be written as:

(34)\begin{array}{l}{e}^{-{\alpha }_D{p}_{\mathrm{6}fD}}\left[{\displaystyle \frac{{\partial }^{\mathrm{2}}{p}_{\mathrm{6}fD}}{\partial y_D^{\mathrm{2}}}}+{\displaystyle \frac{{\partial }^{\mathrm{2}}{p}_{\mathrm{6}fD}}{\partial x_D^{\mathrm{2}}}}+{\displaystyle \frac{{\partial }^{\mathrm{2}}{p}_{\mathrm{6}fD}}{\partial z_D^{\mathrm{2}}}}-{\alpha }_D{\left({\displaystyle \frac{\partial p_{\mathrm{6}fD}}{\partial y_D}}\right)}^{\mathrm{2}}-{\alpha }_D{\left({\displaystyle \frac{\partial p_{\mathrm{6}fD}}{\partial x_D}}\right)}^{\mathrm{2}}-{\alpha }_D{\left({\displaystyle \frac{\partial p_{\mathrm{6}fD}}{\partial z_D}}\right)}^{\mathrm{2}}\right]\\ ={\displaystyle \frac{\mathrm{1}}{{\eta }_{\mathrm{6}D}}}\left[{\omega }_{\mathrm{6}f}{\displaystyle \frac{\partial p_{\mathrm{6}fD}}{\partial t_D}}+{\omega }_{\mathrm{6}m}{\displaystyle \frac{\partial p_{\mathrm{6}mD}}{\partial t_D}}+{\omega }_{\mathrm{6}v}{\displaystyle \frac{\partial p_{\mathrm{6}vD}}{\partial t_D}}\right].\end{array}

The solution of region 1 can be obtained as:

(35){{\displaystyle \frac{\partial {\tilde{\xi }}_{\mathrm{1}D}}{\partial y_D}}|}_{\frac{w_D}{\mathrm{2}}}=-{\beta }_{\mathrm{2}}{\tilde{\xi }}_{FD}\left({\displaystyle \frac{w_D}{\mathrm{2}}}\right),

where

(36){\beta }_{\mathrm{2}}=\sqrt{{\alpha }_{\mathrm{1}}}{\displaystyle \frac{\exp \left[-\sqrt{{\alpha }_{\mathrm{1}}}{\displaystyle \frac{w_D}{\mathrm{2}}}\right]-{\beta }_{\mathrm{1}}\exp \left[\sqrt{{\alpha }_{\mathrm{1}}}{\displaystyle \frac{w_D}{\mathrm{2}}}\right]}{\exp \left[-\sqrt{{\alpha }_{\mathrm{1}}}{\displaystyle \frac{w_D}{\mathrm{2}}}\right]+{\beta }_{\mathrm{1}}\exp \left[\sqrt{{\alpha }_{\mathrm{1}}}{\displaystyle \frac{w_D}{\mathrm{2}}}\right]}},

(37){\beta }_{\mathrm{1}}=\exp [-\sqrt{{\alpha }_{\mathrm{1}}}\mathrm{2}{y}_{\mathrm{1}D}]{\displaystyle \frac{(k_{\mathrm{1}f,h}\sqrt{{\alpha }_{\mathrm{1}}}+k_{\mathrm{2}f,h}\sqrt{{\alpha }_{\mathrm{2}}}\tan h[\sqrt{{\alpha }_{\mathrm{2}}}(y_{\mathrm{1}D}-y_{\mathrm{2}D})])}{(k_{\mathrm{1}f,h}\sqrt{{\alpha }_{\mathrm{1}}}-k_{\mathrm{2}f,h}\sqrt{{\alpha }_{\mathrm{2}}}\tan h[\sqrt{{\alpha }_{\mathrm{2}}}(y_{\mathrm{1}D}-y_{\mathrm{2}D})])}},

where w is wellbore property; β is intermediate variable; FD is dimensionless fracture conductivity.

3.3 Flow in artificial main fracture

In the main artificial fracture, the single-porosity model is applied. Thus, the diffusivity equation can be written as:

(38){\displaystyle \frac{{\partial }^{\mathrm{2}}{p}_{FD}}{\partial x_D^{\mathrm{2}}}}+{\displaystyle \frac{\partial p_{FD}}{\partial y_D}}={\displaystyle \frac{\mathrm{1}}{{\eta }_{FD}}}{\displaystyle \frac{\partial p_{FD}}{\partial t_D}},

The boundary conditions in x direction are:

(39){{\displaystyle \frac{\partial p_{FD}}{\partial x_D}}|}_{x_{\mathrm{1}D}}=\mathrm{0},

(40){{\displaystyle \frac{\partial p_{FD}}{\partial x_D}}|}_{x_{\mathrm{1}D}}=\mathrm{0}{{\displaystyle \frac{\partial p_{FD}}{\partial x_D}}|}_{\mathrm{0}}={\displaystyle \frac{\pi }{F_{CD}}},

where F_CD is dimensionless artificial fracture conductivity, dimensionless.

.

Therefore, the pressure solution for fracture region is:

(41){\tilde{\xi }}_{FD}={\displaystyle \frac{\pi \cos h[\sqrt{{\alpha }_F}(x_D-x_{\mathrm{1}D})]}{s{F}_{CD}\sqrt{{\alpha }_F}\sin h[\sqrt{{\alpha }_F}(x_{\mathrm{1}D})]}}.

Set x_D=0, we can obtain the final solution for well bottom-hole pressure in Laplace domain as following:

(42){\tilde{\xi }}_{wD}={\displaystyle \frac{\pi }{s{F}_{CD}\sqrt{{\alpha }_F}\tan h[\sqrt{{\alpha }_F}(x_{\mathrm{1}D})]}}.

According to the Duhamel and superposition theory, the dimensionless well bottomhole pressure responses incorporating wellbore storage and skin effect in Laplace domain can be written as follows:

(43){\tilde{\xi }}_{wD}(s,C_D,S_c)={\displaystyle \frac{S_c+s{\tilde{\xi }}_{wD}}{s+C_D{s}^{\mathrm{2}}\left(s{\tilde{\xi }}_{wD}+S_c\right)}},

where C_D denotes dimensionless wellbore storage coefficient; S_c denotes skin factor.

Finally, the well bottomhole pressure can be obtained by the following equation:

(44)p_{wD}=-{\displaystyle \frac{\mathrm{1}}{{\alpha }_D}}\ln (\mathrm{1}-{\alpha }_D{\xi }_{wD}).

4 Results and discussion

4.1 Model verification

In order to verify the proposed model, the results obtained by this model and by trilinear flow model in the reference (Brown et al., 2009) are compared with. Let g(s) equal to 1, y_1D equal to y_2D and z_1D equal to z_2D to make the two models the same. In addition, skin factor is not considered in this case. Some other relevant parameters are presented in Table 1.

Figure 4 presents a log/log plot of dimensionless pressure and pressure derivatives vs. dimensionless time solved by the two model. As we can see, the results obtained by this model shows good agreement with the trilinear flow model, which demonstrates the accuracy of this model.

4.2 Flow regimes recognition

In order to obtain the transient pressure type curves, a series of parameters are set as listed in Table 2, and the corresponding type curves are shown in Fig. 5.There are eight flow regimes can be observed as following respectively from Fig. 5:

Regime I: bilinear flow in artificial fracture and in the acidized region (Region 1). The pressure derivative curve’s slope is 1/4 in the regime;

Regime II: first linear flow in the acidized region. The pressure derivative curve shows a straight line with a slope of 1/2;

Regime III: interporosity flow from vug system to natural fracture system in acidized region;

Regime IV: interporosity flow from matrix system to natural fracture system and from vug system to matrix system in acidized region;

Regime V: second linear flow from un-acidized region to acidized region. The pressure derivative curve shows a straight line with a slope of 1/2;

Regime VI: interporosity flow from vug system to natural fracture system in un-acidized region;

Regime VII: interporosity flow from matrix system to natural fracture system and from vug system to matrix system in un-acidized region;

Region VIII: pseudo-steady flow (boundary-dominant flow).

4.3 Parameters analysis

4.3.1 Effect of fracture conductivity

Figure 6 shows the effect of fracture conductivity on dimensionless pressure and pressure derivative. It can be seen that fracture conductivity mainly influences the type curves at the early stage. The greater the fracture conductivity is, the smaller the dimensionless pressure is. It means that lower flow resistance in the artificial main fracture. It can also be seen that the bilinear flow stage disappears when fracture conductivity equals to 10, and the “concaves” in the third and fourth flow regimes have an apparent “dip”, which indicates that fracture conductivity has a significant influence on the pressure behavior of acidized region.

4.3.2 Effect of permeability modulus

Figure 7 illustrates the effect of permeability modulus on transient pressure behavior. As shown in the picture, stress permeability has an apparent impact on the type curves in the late time period. The greater the permeability modulus is, the faster the permeability decreases and it leads to larger pressure consumption at the same production. As the permeability modulus increases from 0.00009 to 0.00019, the boundary dominated flow shifts to an early time.

4.3.3 Effect of the coefficient of interporosity flow between matrix and vug system

Figure 8 depicts the effect of coefficient of cross-flow from vug to fracture on transient pressure behavior. As is known, the coefficient of interporosity flow from vug to fracture in acidized region represents permeability ratio of the vug system in the fracture system. Therefore, a large value of λ_fv1 denotes a high permeability of vug and leads to an early occurrence of first interporosity flow stage. As Fig. 8 shows, the first interporosity flow stage advances with the increase of the parameter λ_fv1. It can also be seen that a bigger will pose a smaller pressure depletion because dimensionless pressure curve descends as parameter λ_fv1 increases from 0.1 to 1.

4.3.4 Effect of storativity ratio of natural fracture in acidized region

Figure 9 shows the effect of storage ratio of natural fracture system in Region 1 on the flow regimes. The natural fracture storativity ratio refers to the ratio of the fluid storing volume in natural fracture system to the total capacity of fluid storing in the reservoir. Furthermore, it reflects the fracture growth level as well, so it can be called fracture intensity as well. Similar to the coefficient of interporosity flow from vug to fracture, storativity ratio of natural fracture in acidized region mainly takes effect during the first interporosity flow stage as well. The smaller the parameter ω_f1 is, the lower the first V-segment becomes.

4.3.5 Effect of width of acidized region

To investigate the width of acidized region separately, the length of acidized region is considered as a constant. Figure 10 indicates the effect of the width of acidized region on the transient pressure behavior. As shown in Fig. 10, the interporosity flow in acidized region happens later with the increasing of the width of acidized region, while the interporosity flow in unacidized region happens earlier. Meanwhile, the larger the size of acidized region is, the lower the dimensionless pressure curve falls, and the fact implies that engineer should try to enlarge the volume of acidized region to reduce the flow obstacle and improve well productivity.

4.3.6 Effect of reservoir size

Figure 11 displays the effect of the reservoir size on the transient pressure behavior. As is shown in the figure, the dimensionless pressure decreases especially in the late time with the increase of reservoir size. It is mainly because that larger reservoir size can provide more fluid to slow down the pressure depletion when well produces at the same rate. Thus, a smaller reservoir size leads to an early boundary-dominant flow.

4.4 Real case application

This section provides an application of the presented model. Well A is located in Tarim oil field in Northwest China. The effective height of well A is 28m. Figure 12 presents the buildup data fitting with the proposed model through an algorithm of auto history matching. As presented in Fig. 12, the model is able to match the real testing data perfectly. The interpretation parameters are shown in Table 3. As we can see in Fig. 12, the permeability of Region 1 is 325mD, which demonstrates that the acidized fracturing makes an effective work.

5 Conclusions

In this paper, an efficient alternative for the analysis of acid fracturing well in carbonate oil reservoir is developed. Acidized regions in acid fracturing carbonate reservoir are considered as the same as SRV concept in hydraulic fracturing tight and shale reservoir. The fracture, matrix and vug are conceptualized as multiple-continuum medium and the composite linear model is used to solve the problem. The stress-sensitive permeability and partially penetrated fracture are both taken into consideration. Based on the investigation of this paper, the following conclusions can be drawn:

1) As for an acidized fracturing well in a carbonate oil reservoir, the flow regimes can be subdivided into the following stages, e.g., bilinear flow stage; linear flow stage; cross flow in acidized region; second linear flow from un-acidized region to acidized region; cross flow in un-acidized region; boundary-dominant flow.

2) Parameters sensitivity analysis demonstrates that fracture conductivity mainly influences the type curves in the early stage while permeability modulus influences the late stage; storativity ratio and interporosity flow mainly influences the cross flow stage; large acidized region size can reduce the pressure loss and enhance well productivity.

3) The presented solution combining with the algorithm of auto history matching can be used to obtain reservoir parameters and evaluate the effectiveness of acid fracturing measures.

References

[1] Abass H H, Al-Mulhem A A, Alqam M H, Khan M R (2006). Acid fracturing or proppant fracturing in carbonate formation? A rock mechanics view. In: SPE Annual Technical Conference and Exhibition, San Antonio, USA

[2] Abdassah D, Ershaghi I (1986). Triple-porosity system for representing naturally fractured reservoirs. SPE Form Eval, 1(02): 113–127

[3] Brown M, Ozkan E, Raghavan R, Kazemi H (2011). Practical solutions for pressure transient responses of fractured horizontal wells in unconventional reservoirs. SPE Reservoir Eval Eng, 14(06): 663–676

[4] Camacho-Velázquez R, Vasquez-Cruz M, Castrejon-Aivar R, Arana-Ortiz V (2002). Pressure transient and decline curve behaviors in naturally fractured vuggy carbonate reservoirs. In: SPE Annual Technical Conference and Exhibition, San Antonio, USA

[5] Dora P R (2008). Pressure behavior of a system containing multiple vertical fractures, Dissertation for the Doctoral Degree. Norman: University of Oklahoma

[6] Fredd C N (2000). Dynamic model of wormhole formation demonstrates conditions for effective skin reduction during carbonate matrix acidizing. In: SPE Permian Basin Oil and Gas Recovery Conference, Midland, USA

[7] Fuentes-Cruz G, Camacho-Velazquez R, Vasquez-Cruz M (2004). Pressure transient and decline curve behaviors for partially penetrating wells completed in naturally fractured-vuggy reservoirs. In: SPE International Petroleum Conference in Mexico, Puebla Pue., Mexico

[8] Guo J, Nie R, Jia Y (2012). Dual permeability fow behavior for modeling horizontal well production in fractured vuggy carbonate reservoirs. J Hydrol (Amst), 464-465: 281–293

[9] Igbokoyi A O, Tiab D (2008). Pressure transient analysis in partially penetrating infinite conductivity hydraulic fractures in naturally fractured reservoirs. In: SPE Annual Technical Conference and Exhibition, Denver, USA

[10] Kossack C A, Gurpinar O A (2001). Methodology for simulation of vuggy and fractured reservoirs. In: SPE Reservoir Simulation Symposium, Houston, USA.

[11] Kang Z, Wu Y S, Li J, Wu Y, Zhang J, Wang G (2006). Modeling multiphase flow in naturally fractured vuggy petroleum reservoirs. In: SPE Annual Technical Conference and Exhibition, San Antonio, USA

[12] Liu J, Bodvarsson G S, Wu Y S (2003). Analysis of flow behavior in fractured lithophysal reservoirs. J ContamHydrol, 62–63(1): 189–211

[13] 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

[14] Lei G, Liao Q, Zhang D A (2018). New analytical model for flow in acidized fractured-vuggy porous media. Scientifc Reports, 9(1): 8293

[15] Pedrosa O A (1986) Pressure transient response in stress-sensitive formations. In: SPE California Regional Meeting, Oakland, USA

[16] Raghavan R, Uraiet A, Tomas G W (1978).Vertical fracture height: effect on transient flow behavior. Soc Petrol Eng J, 18(4): 265–277

[17] Rodriguez F, Horne R N, Cinco-Ley H (1984). Partially penetrating fractures: pressure transient analysis of an infinite conductivity fracture. In: SPE Annual Technical Conference and Exhibition, Houston, USA

[18] Stalgorova E, Mattar L (2012).Practical analytical model to simulate production of horizontal wells with branch fractures. In: SPE Canadian Unconventional Resource Conference, Calgary, Canada

[19] Stehfest H (1970). Numerical inversion of Laplace transforms. Commun ACM, 13(1): 47–49

[20] Tao H, Zhang L, Liu Q, Deng Q, Luo M, Zhao Y (2018). An analytical flow model for heterogeneous multi-fractured systems in shale gas reservoirs. Energies, 3422(12):

[21] Wu Y S, Qin G, Ewing R E, Efendiev Y, Kang Z, Ren Y (2006). A multiple-continuum approach for modeling multiphase flow in naturally fractured vuggy petroleum reservoirs. In: International Oil & Gas Conference and Exhibition in China, Beijing, China

[22] Wu Y S, Ehlig-Economides C A, Qin G, Kang Z, Zhang W, Ajayi B T, Tao Q (2007). A triple-continuum pressure-transient model for a naturally fractured vuggy reservoir. In: SPE Annual Technical Conference and Exhibition, Anaheim, USA

[23] Wang L, Wang X, Luo E, Wang J (2014). Analytical modeling of fow behavior for wormholes in naturally fractured-vuggy porous media. Transp Porous Media, 105(3): 539–558

[24] Wang L, Chen X, Xia Z (2018a). A novel semi-analytical model for multi-branched fractures in naturally fractured-vuggy reservoirs. Sci Rep, 8(1): 11586

[25] Wang M X, Fan Z F, Dong X Y, Song H, Zhao W, Xu G (2018b). Analysis of flow behavior for acid fracturing wells in fractured-vuggy carbonate reservoirs. Mathematical Problems in Engineering, 2018(PT.3): 6431910.1–6431910.20

[26] Wang Y, Yi X (2017). Transient pressure behavior of a fractured vertical well with a finite-conductivity fracture in triple media carbonate reservoir. J Porous Media, 20(8): 707–722

[27] Wu Y, Cheng L, Huang S, Fang S, Jia P, Wang S (2019). A semianalytical model for simulating fluid flow in naturally fractured reservoirs with nonhomogeneous vugs and fractures. SPE J, 24(01): 334–348

[28] Wang S, Cheng L, Xue Y, Huang S, Wu Y, Jia P, Sun Z (2018). A semi-analytical method for simulating two-phase flow performance of horizontal volatile oil wells in fractured carbonate reservoirs. Energies, 11(10): 2700

[29] Xing C, Yin H, Liu K, Li X, Fu J (2018). Well test analysis for fractured and vuggy carbonate reservoirs of well drilling in large scale cave. Energies, 11(1): 80

[30] Xu Y, Li P, Du X, Lu Z, Li D, Lu D (2019). A novel analytical well test model for fractured vuggy carbonate reservoirs considering the coupling between oil flow and wave propagation. J Petrol Sci Eng, 173: 447–461

[31] Yang Y, Liu Z, Sun Z, An S, Zhang W, Liu P, Yao J, Ma J (2017). Research on Stress Sensitivity of Fractured Carbonate Reservoirs Based on CT Technology. Energies, 10(11): 1833

[32] Yuan J, Jiang R, Zhang W (2018). The workflow to analyze hydraulic fracture effect on hydraulic fractured horizontal well production in composite formation system. Advances in Geo-Energy Research, 2(3): 319–342

[33] Yao J, Huang Z, Li Y, Wang C, Lv X ( 2010). Discrete fracture-vug network model for modeling fluid flow in fractured vuggy porous media. In: International Oil and Gas Conference and Exhibition in China, Beijing, China

[34] Zhang F, An M, Yan B, Wang Y (2017). Modeling the depletion of fractured vuggy carbonate reservoir by coupling geomechanics with reservoir flow. In: SPE Reservoir Characterisation and Simulation Conference and Exhibition

[35] Zhang H,Wang X,Wang L (2015). An analytical solution of partially penetrating hydraulic fractures in a box-shaped reservoir. Math Probl in Eng, 2015:11 (in Chinese)

[36] Zeng Y, Wang Q, Ning Z, Sun H (2018). A mathematical pressure transient analysis model for multiple fractured horizontal wells in shale gas reservoirs. Geofluids, 2018: 8065949

[37] Zeng J, Wang X, Guo J, Zeng F (2017). Composite linear flow model for multi-fractured horizontal wells in heterogeneous shale reservoir. J Nat Gas Sci Eng, 38: 527–548

RIGHTS & PERMISSIONS

Higher Education Press