# Analysis of Flow Behavior for Acid Fracturing Wells in Fractured-Vuggy Carbonate Reservoirs

Authors: Mingxian Wang; Zifei Fan; Xuyang Dong; Heng Song; Wenqi Zhao; Guoqing Xu
Journal: Mathematical Problems in Engineering
 (2018)

---

## Abstract

This study develops a mathematical model for transient flow analysis of acid fracturing wells in fractured-vuggy carbonate reservoirs. This model considers a composite system with the inner region containing finite number of artificial fractures and wormholes and the outer region showing a triple-porosity medium. Both analytical and numerical solutions are derived in this work, and the comparison between two solutions verifies the model accurately. Flow behavior is analyzed thoroughly by examining the standard log-log type curves. Flow in this composite system can be divided into six or eight main flow regimes comprehensively. Three or two characteristic V-shaped segments can be observed on pressure derivative curves. Each V-shaped segment corresponds to a specific flow regime. One or two of the V-shaped segments may be absent in particular cases. Effects of interregional diffusivity ratio and interregional conductivity ratio on transient responses are strong in the early-flow period. The shape and position of type curves are also influenced by interporosity coefficients, storativity ratios, and reservoir radius significantly. Finally, we show the differences between our model and the similar model with single fracture or without acid fracturing and further investigate the pseudo-skin factor caused by acid fracturing.


---

## Body

## 1. Introduction

Acid fracturing is an important stimulation technique that has been widely used in developing fractured-vuggy carbonate reservoirs (FVCRs); thus, evaluating the flow behavior of acid fracturing wells in these reservoirs becomes necessary for further improving the wells’ performance.In general, fractured-vuggy carbonate rock can be considered as a triple-porosity medium, in which natural fractures, matrices, and vugs constitute the flow system [1–3]. Recently, many scholars have made efforts to simulate the fluid transport behavior in such reservoirs through numerical and analytical methods. Wu et al. [2, 3] proposed an analytical approach for transient flow analysis in the naturally fractured-vuggy reservoirs, focusing on single-phase flow and multiphase flow, respectively. Cobett et al. [4] developed a numerical model of fractured carbonate rocks and showed that numerical well testing had its limitations, especially when simulating disperse vugs. Gulbransen et al. [5] presented a multiscale mixed finite element (MsMFE) method for modeling of vuggy and naturally fractured reservoirs. Yao et al. [6] obtained the solution of the discrete fracture-vug network (DFVN) model through a standard Galerkin finite element method and provided a way of modeling realistic fluid flow in the fractured-vuggy porous media. Guo et al. [7] established a model for a horizontal well in a triple-porosity and dual-permeability system, and their results showed that type curves were dominated by external boundary as well as the permeability ratio of fracture system to the sum of fracture and matrix systems. These models provide a series of effective methods for treating the triple-porosity medium, and thus we can gain a comprehensive insight into the flow behavior in FVCRs.However, dealing with acid fracturing wells in these reservoirs is not an easy work. Generally, carbonate minerals are usually composed of dolomite and calcite which are easy to be dissolved by hydrochloric acid [8, 9]. As one of the fundamental ways to stimulate carbonate reservoirs, acid fracturing alters the flow pattern in the reservoirs from a radial one to a linear one [10] and has dual characteristics of hydraulic fracturing and acidification. On the one hand, like hydraulic fracturing, multiple artificial fractures may occur near the wellbore when injecting fluids into a formation at a high pressure [11–14]. Wang and Yi [15] investigated the flow behavior of a well with a finite-conductivity fracture in a carbonate reservoir. However, several field and laboratory observations have shown very complex fracture paths and presence of multiple fractures and multisegmented fractures near the wellbore after hydraulic fracturing [16–20]. On the other hand, the injected acid etches the face of the crack and these acid-etched patterns, namely, wormholes, remain to provide flow channels when hydraulic pressure is released [8, 10, 21–23]. Therefore, wormholes, together with artificial fractures, form the main channels of the fluid flowing into the wellbore. At present, acid fracturing has been widely implemented to enhance productivity of wells with damaged zones or produce from tight carbonate reservoirs [23–25]. Unfortunately, complex seepage systems created by this stimulation pose challenges to study the transport processes in underground natural-resources recovery.According to the concept of composite reservoir, when there is a discontinuity in the material properties in the radial direction from a well, the reservoir can be regarded as a composite reservoir [26–28]. Considering that acid fracturing creates a high permeability zone near the wellbore, the case of a well with wormholes and artificial fractures in FVCRs can be represented by a composite reservoir. Therefore, the idea of composite reservoir can be used to deal with the case in this work. To the best of our knowledge, in the literature there is no similar composite model to study the flow behavior of vertical wells within wormholes and artificial fractures in FVCRs.The purpose of this work is to conduct flow analysis of acid fracturing wells in FVCRs. We consider a composite reservoir with the inner region containing finite number of artificial fractures and wormholes, flow in both of which is linear flow, and the outer region showing a triple-porosity medium. Firstly, analytical solution of the proposed model is derived in the Laplace space, and then numerical solution of the same model is obtained with a finite difference method. The comparison results between the analytical and numerical solutions at different parameters verify the model accurately. Then, using the analytical solution, we plot the type curves and identify the main flow regimes. Effects of interregional diffusivity ratio, interregional conductivity ratio, fluid storativity ratio and interporosity coefficient of each medium, and reservoir radius on the transient response are further analyzed. Finally, we compare the differences between our model and the similar model with single fracture or without acid fracturing and also investigate the pseudo-skin factor caused by acid fracturing in detail.


## 2. Mathematical Model

FVCRs with multiple artificial fractures and wormholes are naturally structured by fracture system, vug system, and matrix system. Figure1(a) shows a physical model of an acid fracturing well in a FVCR. The repetitive matrices and vugs are separated by a set of natural fractures. Using the concept of composite reservoir, the proposed model can be divided into two concentric regions: the inner region (Region I) and the outer region (Region II) (seen in Figure 1(b)). Figure 1(c) shows a sketch of the flow in the composite system. Similar to the dual-porosity model [29], a triple-porosity model is introduced into the outer region [3, 4]. Natural fractures are considered as the main flow pathways and connected with artificial fractures or wormholes directly and ultimately fluid flow from artificial fractures or wormholes into the wellbore. Wormholes are doped with artificial fractures near the wellbore, and it is difficult yet meaningless to determine their exact numbers. In practice, wormholes can be regarded as artificial fractures with larger aperture. By assuming Darcy linear flow in both kinds of flow channels, wormholes can be approximated to artificial fractures.Figure 1
Schematic diagram: (a) an acid fracturing well in a fractured-vuggy carbonate reservoir; (b) a composite reservoir model with two concentric regions; and (c) a sectional view of fluid transport in the composite system.
(a) (b) (c)To simplify the mathematical problem, the basic assumptions are made:( 1 ) In the inner region, there are a finite number, n, of artificial fractures and wormholes with the same aperture, b, permeability, k1, and storage capacity, (ϕc)1, which are conceptualized as vertical plates extending from the top to the bottom. This region’s radius is rf, equal to the length of artificial fractures or wormholes.( 2 ) In the outer region, the permeability and storage capacity for natural fractures are k2f and (ϕc)2f, for vugs are k2v and (ϕc)2v, and for matrices are k2m and (ϕc)2m. Natural fractures, vugs, and matrices are uniformly distributed in this region.( 3 ) All rock properties, such as permeability and porosity, are constant. The flow in the reservoir is considered to be a single-phase flow with slightly compressible fluid of constant viscosity μ and constant formation volume factor B.( 4 ) The reservoir is of uniform thickness with impermeable upper and lower boundaries and closed with a radius, equal to re.( 5 ) A vertical well is located in the center of the inner region and intercepted by these artificial fractures and wormholes. The well is produced at a constant rate q and no wellbore storage is considered.( 6 ) Isothermal and Darcy flow are assumed in two regions.


### 2.1. Establishment of Mathematical Model

#### 2.1.1. Flow in Region I

In the inner region, flow in the artificial fractures and wormholes follows a linear flow pattern; thus, the governing equation for this region can be written as(1)k1μ∂2p1∂r2=ϕct1∂p1∂t.Initial condition is(2)p1r,0=pi0≤r≤rf.Well production condition for constant rate is(3)nbhk1μ∂p1∂rr→0=qB.


#### 2.1.2. Flow in Region II

Based on the triple-porosity model, any infinitesimal volume in the outer region contains a large proportion of three constitutive media. Each point in the model is assigned with three pressures. Each system interacts independently with the other two systems. The matrix and vug systems provide the main storage space but have no direct contributions to flow and transport.The governing equations for natural fracture, vug, and matrix system are(4)k2fμ1r∂∂rr∂p2f∂r=ϕct2f∂p2f∂t-q2mf-q2vf,(5)0=ϕct2v∂p2v∂t+q2mv+q2vf(6)0=ϕct2m∂p2m∂t+q2mf-q2mv,respectively. The subscripts 2f, 2v, and 2m refer to natural fracture, vug, and matrix system, respectively; and q2mf, q2mv, and q2vf are the volumetric flux from matrix to natural fracture, from vug to matrix, and from vug to natural fracture, respectively. Assuming that fluid flow between two systems follows the pseudo-steady state flow [29], q2mf, q2mv, and q2vf can be given as(7)q2mf=γmfk2mp2m-p2fμ,(8)q2mv=γmvk2vp2v-p2mμ,(9)q2vf=γvfk2vp2v-p2fμ,where γmf, γmv, and γvf are the interporosity flow shape factors, which depend on the geometry of the interporosity flow and have dimensions of reciprocal area.By substituting (7)–(9) into (4)–(6), the previous governing equations can be simplified as(10)k2fμ1r∂∂rr∂p2f∂r=ϕct2f∂p2f∂t+ϕct2v∂p2v∂t+ϕct2m∂p2m∂t,(11)0=ϕct2v∂p2v∂t+γmvk2vp2v-p2mμ+γvfk2vp2v-p2fμ,(12)0=ϕct2m∂p2m∂t+γmfk2mp2m-p2fμ-γmvk2vp2v-p2mμ,respectively.The initial and closed external boundary conditions in the outer region are(13)p2fr,0=p2vr,0=p2mr,0=pirf≤r≤re∂p2f∂rr=re=0.


#### 2.1.3. Interface Conditions

Natural fractures in the outer region are connected with artificial fractures and wormholes in the inner region directly. Because of the continuity of pressure and flux at the interface, the following conditions must hold along the surface:(14)p1r=rf=p2fr=rf,nbk1μ∂p1∂rr=rf=2πrk2fμ∂p2f∂rr=rf.


### 2.2. Dimensionless Mathematical Model

#### 2.2.1. Dimensionless Definitions

Dimensionless pressure of artificial fracture and wormhole system in the inner region is(15)p1D=2πk2fhpi-p1qμB.Dimensionless pressure of natural fracture system is(16)p2fD=2πk2fhpi-p2fqμB.Dimensionless pressure of vug system is(17)p2vD=2πk2fhpi-p2vqμB.Dimensionless pressure of matrix system is(18)p2mD=2πk2fhpi-p2mqμB.Dimensionless time is(19)tD=k2ftμϕct2frf2.Dimensionless radius is(20)rD=rrf.Dimensionless radius of external boundary is(21)reD=rerf.Dimensionless radius of wellbore is(22)rwD=rwrf.Interregional diffusivity ratio is(23)α=η2fη1=k2fϕct1k1ϕct2f.Interregional conductivity ratio is(24)β=2πk2frfnbk1.Fluid storativity ratio of fracture system is(25)ω2f=ϕct2fϕct2f+ϕct2m+ϕct2v.Fluid storativity ratio of vug system is(26)ω2v=ϕct2vϕct2f+ϕct2m+ϕct2v.Fluid storativity ratio of matrix system is(27)ω2m=ϕct2mϕct2f+ϕct2m+ϕct2v.Interporosity coefficient of matrix system to fracture system is(28)λmf=γmfk2mrf2k2f.Interporosity coefficient of vug system to fracture system is(29)λvf=γvfk2vrf2k2f.Interporosity coefficient of vug system to matrix system is(30)λmv=γmvk2vrf2k2f.It is worth noting that there is a certain relationship among three storativity ratios:ω2f+ω2m+ω2v=1. In a FVCR, the interporosity coefficient of vug system to fracture system λvf should be larger than that of matrix system to fracture system λmf. The differences in these interporosity coefficients lead to different response time of each system on type curves [1].


#### 2.2.2. Dimensionless Mathematical Model in Real Space

Introducing the dimensionless variables into (1)–(3) and (10)–(14), we obtain the following dimensionless equations.(  1) Region I. The governing differential equation, initial condition, and inner boundary condition become(31)∂2p1D∂rD2=α∂p1D∂tD,(32)p1DrD,0=00≤rD≤rfD,(33)∂p1D∂rDrD→0=-β,respectively.(  2) Region II. The dimensionless governing differential equations of natural fracture, vug, and matrix system are(34)∂2p2fD∂rD2+1rD∂p2fD∂rD=∂p2fD∂tD+ω2mω2f∂p2mD∂tD+ω2vω2f∂p2vD∂tD,(35)ω2vω2f∂p2vD∂tD+λmvp2vD-p2mD+λvfp2vD-p2fD=0,(36)ω2mω2f∂p2mD∂tD+λmfp2mD-p2fD-λmvp2vD-p2mD=0,respectively.The initial and external boundary conditions become(37)p2fDrD,0=p2mDrD,0=p2vDrD,0=0rfD≤rD≤reD,(38)∂p2fD∂rDrD=reD=0.(  3) Interface Conditions. The dimensionless equations of pressure and flux at the interface are given by(39)p1DrD=1=p2fDrD=1,(40)∂p1D∂rDrD=1=βrD∂p2fD∂rDrD=1.


## 2.1. Establishment of Mathematical Model

### 2.1.1. Flow in Region I

In the inner region, flow in the artificial fractures and wormholes follows a linear flow pattern; thus, the governing equation for this region can be written as(1)k1μ∂2p1∂r2=ϕct1∂p1∂t.Initial condition is(2)p1r,0=pi0≤r≤rf.Well production condition for constant rate is(3)nbhk1μ∂p1∂rr→0=qB.


### 2.1.2. Flow in Region II

Based on the triple-porosity model, any infinitesimal volume in the outer region contains a large proportion of three constitutive media. Each point in the model is assigned with three pressures. Each system interacts independently with the other two systems. The matrix and vug systems provide the main storage space but have no direct contributions to flow and transport.The governing equations for natural fracture, vug, and matrix system are(4)k2fμ1r∂∂rr∂p2f∂r=ϕct2f∂p2f∂t-q2mf-q2vf,(5)0=ϕct2v∂p2v∂t+q2mv+q2vf(6)0=ϕct2m∂p2m∂t+q2mf-q2mv,respectively. The subscripts 2f, 2v, and 2m refer to natural fracture, vug, and matrix system, respectively; and q2mf, q2mv, and q2vf are the volumetric flux from matrix to natural fracture, from vug to matrix, and from vug to natural fracture, respectively. Assuming that fluid flow between two systems follows the pseudo-steady state flow [29], q2mf, q2mv, and q2vf can be given as(7)q2mf=γmfk2mp2m-p2fμ,(8)q2mv=γmvk2vp2v-p2mμ,(9)q2vf=γvfk2vp2v-p2fμ,where γmf, γmv, and γvf are the interporosity flow shape factors, which depend on the geometry of the interporosity flow and have dimensions of reciprocal area.By substituting (7)–(9) into (4)–(6), the previous governing equations can be simplified as(10)k2fμ1r∂∂rr∂p2f∂r=ϕct2f∂p2f∂t+ϕct2v∂p2v∂t+ϕct2m∂p2m∂t,(11)0=ϕct2v∂p2v∂t+γmvk2vp2v-p2mμ+γvfk2vp2v-p2fμ,(12)0=ϕct2m∂p2m∂t+γmfk2mp2m-p2fμ-γmvk2vp2v-p2mμ,respectively.The initial and closed external boundary conditions in the outer region are(13)p2fr,0=p2vr,0=p2mr,0=pirf≤r≤re∂p2f∂rr=re=0.


### 2.1.3. Interface Conditions

Natural fractures in the outer region are connected with artificial fractures and wormholes in the inner region directly. Because of the continuity of pressure and flux at the interface, the following conditions must hold along the surface:(14)p1r=rf=p2fr=rf,nbk1μ∂p1∂rr=rf=2πrk2fμ∂p2f∂rr=rf.


## 2.1.1. Flow in Region I

In the inner region, flow in the artificial fractures and wormholes follows a linear flow pattern; thus, the governing equation for this region can be written as(1)k1μ∂2p1∂r2=ϕct1∂p1∂t.Initial condition is(2)p1r,0=pi0≤r≤rf.Well production condition for constant rate is(3)nbhk1μ∂p1∂rr→0=qB.


## 2.1.2. Flow in Region II

Based on the triple-porosity model, any infinitesimal volume in the outer region contains a large proportion of three constitutive media. Each point in the model is assigned with three pressures. Each system interacts independently with the other two systems. The matrix and vug systems provide the main storage space but have no direct contributions to flow and transport.The governing equations for natural fracture, vug, and matrix system are(4)k2fμ1r∂∂rr∂p2f∂r=ϕct2f∂p2f∂t-q2mf-q2vf,(5)0=ϕct2v∂p2v∂t+q2mv+q2vf(6)0=ϕct2m∂p2m∂t+q2mf-q2mv,respectively. The subscripts 2f, 2v, and 2m refer to natural fracture, vug, and matrix system, respectively; and q2mf, q2mv, and q2vf are the volumetric flux from matrix to natural fracture, from vug to matrix, and from vug to natural fracture, respectively. Assuming that fluid flow between two systems follows the pseudo-steady state flow [29], q2mf, q2mv, and q2vf can be given as(7)q2mf=γmfk2mp2m-p2fμ,(8)q2mv=γmvk2vp2v-p2mμ,(9)q2vf=γvfk2vp2v-p2fμ,where γmf, γmv, and γvf are the interporosity flow shape factors, which depend on the geometry of the interporosity flow and have dimensions of reciprocal area.By substituting (7)–(9) into (4)–(6), the previous governing equations can be simplified as(10)k2fμ1r∂∂rr∂p2f∂r=ϕct2f∂p2f∂t+ϕct2v∂p2v∂t+ϕct2m∂p2m∂t,(11)0=ϕct2v∂p2v∂t+γmvk2vp2v-p2mμ+γvfk2vp2v-p2fμ,(12)0=ϕct2m∂p2m∂t+γmfk2mp2m-p2fμ-γmvk2vp2v-p2mμ,respectively.The initial and closed external boundary conditions in the outer region are(13)p2fr,0=p2vr,0=p2mr,0=pirf≤r≤re∂p2f∂rr=re=0.


## 2.1.3. Interface Conditions

Natural fractures in the outer region are connected with artificial fractures and wormholes in the inner region directly. Because of the continuity of pressure and flux at the interface, the following conditions must hold along the surface:(14)p1r=rf=p2fr=rf,nbk1μ∂p1∂rr=rf=2πrk2fμ∂p2f∂rr=rf.


## 2.2. Dimensionless Mathematical Model

### 2.2.1. Dimensionless Definitions

Dimensionless pressure of artificial fracture and wormhole system in the inner region is(15)p1D=2πk2fhpi-p1qμB.Dimensionless pressure of natural fracture system is(16)p2fD=2πk2fhpi-p2fqμB.Dimensionless pressure of vug system is(17)p2vD=2πk2fhpi-p2vqμB.Dimensionless pressure of matrix system is(18)p2mD=2πk2fhpi-p2mqμB.Dimensionless time is(19)tD=k2ftμϕct2frf2.Dimensionless radius is(20)rD=rrf.Dimensionless radius of external boundary is(21)reD=rerf.Dimensionless radius of wellbore is(22)rwD=rwrf.Interregional diffusivity ratio is(23)α=η2fη1=k2fϕct1k1ϕct2f.Interregional conductivity ratio is(24)β=2πk2frfnbk1.Fluid storativity ratio of fracture system is(25)ω2f=ϕct2fϕct2f+ϕct2m+ϕct2v.Fluid storativity ratio of vug system is(26)ω2v=ϕct2vϕct2f+ϕct2m+ϕct2v.Fluid storativity ratio of matrix system is(27)ω2m=ϕct2mϕct2f+ϕct2m+ϕct2v.Interporosity coefficient of matrix system to fracture system is(28)λmf=γmfk2mrf2k2f.Interporosity coefficient of vug system to fracture system is(29)λvf=γvfk2vrf2k2f.Interporosity coefficient of vug system to matrix system is(30)λmv=γmvk2vrf2k2f.It is worth noting that there is a certain relationship among three storativity ratios:ω2f+ω2m+ω2v=1. In a FVCR, the interporosity coefficient of vug system to fracture system λvf should be larger than that of matrix system to fracture system λmf. The differences in these interporosity coefficients lead to different response time of each system on type curves [1].


### 2.2.2. Dimensionless Mathematical Model in Real Space

Introducing the dimensionless variables into (1)–(3) and (10)–(14), we obtain the following dimensionless equations.(  1) Region I. The governing differential equation, initial condition, and inner boundary condition become(31)∂2p1D∂rD2=α∂p1D∂tD,(32)p1DrD,0=00≤rD≤rfD,(33)∂p1D∂rDrD→0=-β,respectively.(  2) Region II. The dimensionless governing differential equations of natural fracture, vug, and matrix system are(34)∂2p2fD∂rD2+1rD∂p2fD∂rD=∂p2fD∂tD+ω2mω2f∂p2mD∂tD+ω2vω2f∂p2vD∂tD,(35)ω2vω2f∂p2vD∂tD+λmvp2vD-p2mD+λvfp2vD-p2fD=0,(36)ω2mω2f∂p2mD∂tD+λmfp2mD-p2fD-λmvp2vD-p2mD=0,respectively.The initial and external boundary conditions become(37)p2fDrD,0=p2mDrD,0=p2vDrD,0=0rfD≤rD≤reD,(38)∂p2fD∂rDrD=reD=0.(  3) Interface Conditions. The dimensionless equations of pressure and flux at the interface are given by(39)p1DrD=1=p2fDrD=1,(40)∂p1D∂rDrD=1=βrD∂p2fD∂rDrD=1.


## 2.2.1. Dimensionless Definitions

Dimensionless pressure of artificial fracture and wormhole system in the inner region is(15)p1D=2πk2fhpi-p1qμB.Dimensionless pressure of natural fracture system is(16)p2fD=2πk2fhpi-p2fqμB.Dimensionless pressure of vug system is(17)p2vD=2πk2fhpi-p2vqμB.Dimensionless pressure of matrix system is(18)p2mD=2πk2fhpi-p2mqμB.Dimensionless time is(19)tD=k2ftμϕct2frf2.Dimensionless radius is(20)rD=rrf.Dimensionless radius of external boundary is(21)reD=rerf.Dimensionless radius of wellbore is(22)rwD=rwrf.Interregional diffusivity ratio is(23)α=η2fη1=k2fϕct1k1ϕct2f.Interregional conductivity ratio is(24)β=2πk2frfnbk1.Fluid storativity ratio of fracture system is(25)ω2f=ϕct2fϕct2f+ϕct2m+ϕct2v.Fluid storativity ratio of vug system is(26)ω2v=ϕct2vϕct2f+ϕct2m+ϕct2v.Fluid storativity ratio of matrix system is(27)ω2m=ϕct2mϕct2f+ϕct2m+ϕct2v.Interporosity coefficient of matrix system to fracture system is(28)λmf=γmfk2mrf2k2f.Interporosity coefficient of vug system to fracture system is(29)λvf=γvfk2vrf2k2f.Interporosity coefficient of vug system to matrix system is(30)λmv=γmvk2vrf2k2f.It is worth noting that there is a certain relationship among three storativity ratios:ω2f+ω2m+ω2v=1. In a FVCR, the interporosity coefficient of vug system to fracture system λvf should be larger than that of matrix system to fracture system λmf. The differences in these interporosity coefficients lead to different response time of each system on type curves [1].


## 2.2.2. Dimensionless Mathematical Model in Real Space

Introducing the dimensionless variables into (1)–(3) and (10)–(14), we obtain the following dimensionless equations.(  1) Region I. The governing differential equation, initial condition, and inner boundary condition become(31)∂2p1D∂rD2=α∂p1D∂tD,(32)p1DrD,0=00≤rD≤rfD,(33)∂p1D∂rDrD→0=-β,respectively.(  2) Region II. The dimensionless governing differential equations of natural fracture, vug, and matrix system are(34)∂2p2fD∂rD2+1rD∂p2fD∂rD=∂p2fD∂tD+ω2mω2f∂p2mD∂tD+ω2vω2f∂p2vD∂tD,(35)ω2vω2f∂p2vD∂tD+λmvp2vD-p2mD+λvfp2vD-p2fD=0,(36)ω2mω2f∂p2mD∂tD+λmfp2mD-p2fD-λmvp2vD-p2mD=0,respectively.The initial and external boundary conditions become(37)p2fDrD,0=p2mDrD,0=p2vDrD,0=0rfD≤rD≤reD,(38)∂p2fD∂rDrD=reD=0.(  3) Interface Conditions. The dimensionless equations of pressure and flux at the interface are given by(39)p1DrD=1=p2fDrD=1,(40)∂p1D∂rDrD=1=βrD∂p2fD∂rDrD=1.


## 3. Analytical Solution

### 3.1. Laplace Transform of Dimensionless Mathematical Model

The Laplace transform with respect totD is based on the following function:(41)p~DrD,s=∫0∞pDrD,tDe-stDdtD, where s is the Laplace transform variable.For the inner region, applying the Laplace transform to (31)–(33), we have(42)∂2p~1D∂rD2=αsp~1D,(43)∂p~1D∂rDrD→0=-βs.For the outer region, by employing the Laplace transform, (34)–(37) can be simplified as(44)∂2p~2fD∂rD2+1rD∂p~2fD∂rD-sgsp~2fD=0,where(45)gs=1+λmf+λvfs+w2f1-w2f/w2vw2mλmfλmv+λmvλvf+λmfλvfs2+λmf/w2m+λvf/w2v+1/w2m+1/w2vλmvw2fs+w2f2/w2mw2vλmfλmv+λmvλvf+λmfλvf.The external boundary condition in Laplace space is(46)∂p~2fD∂rDrD=reD=0.Similarly, taking Laplace transform for the dimensionless interface boundary conditions, we have(47)p~1DrD=1=p~2fDrD=1,∂p~1D∂rDrD=1=β∂p~2fD∂rDrD=1.


### 3.2. Analytical Solution in Laplace Space

The general solutions for (42) and (44) are given by(48)p~1D=AcoshsαrD+BsinhsαrD,p~2fD=CI0sgsrD+DK0sgsrD,where I0 and K0 are modified Bessel functions of zero order of the first and second kind, respectively.Combining the general solutions with all the boundary conditions, we can get the analytical solution and then the wellbore pressure in Laplace domain can be obtained:(49)p~wDs=βssαx2sinh⁡sα1-rwD-x1cosh⁡sα1-rwDx2cosh⁡sα-x1sinh⁡sα, where(50)x1=K0sgs+K1sgsreDI1sgsreDI0sgs,x2=-βgsαK1sgs-K1sgsreDI1sgsreDI1sgs.In real space, the dimensionless bottomhole pressure (BHP) can be obtained by using Stehfest numerical inversion [30] to convert p~wDs to pwDtD.


## 3.1. Laplace Transform of Dimensionless Mathematical Model

The Laplace transform with respect totD is based on the following function:(41)p~DrD,s=∫0∞pDrD,tDe-stDdtD, where s is the Laplace transform variable.For the inner region, applying the Laplace transform to (31)–(33), we have(42)∂2p~1D∂rD2=αsp~1D,(43)∂p~1D∂rDrD→0=-βs.For the outer region, by employing the Laplace transform, (34)–(37) can be simplified as(44)∂2p~2fD∂rD2+1rD∂p~2fD∂rD-sgsp~2fD=0,where(45)gs=1+λmf+λvfs+w2f1-w2f/w2vw2mλmfλmv+λmvλvf+λmfλvfs2+λmf/w2m+λvf/w2v+1/w2m+1/w2vλmvw2fs+w2f2/w2mw2vλmfλmv+λmvλvf+λmfλvf.The external boundary condition in Laplace space is(46)∂p~2fD∂rDrD=reD=0.Similarly, taking Laplace transform for the dimensionless interface boundary conditions, we have(47)p~1DrD=1=p~2fDrD=1,∂p~1D∂rDrD=1=β∂p~2fD∂rDrD=1.


## 3.2. Analytical Solution in Laplace Space

The general solutions for (42) and (44) are given by(48)p~1D=AcoshsαrD+BsinhsαrD,p~2fD=CI0sgsrD+DK0sgsrD,where I0 and K0 are modified Bessel functions of zero order of the first and second kind, respectively.Combining the general solutions with all the boundary conditions, we can get the analytical solution and then the wellbore pressure in Laplace domain can be obtained:(49)p~wDs=βssαx2sinh⁡sα1-rwD-x1cosh⁡sα1-rwDx2cosh⁡sα-x1sinh⁡sα, where(50)x1=K0sgs+K1sgsreDI1sgsreDI0sgs,x2=-βgsαK1sgs-K1sgsreDI1sgsreDI1sgs.In real space, the dimensionless bottomhole pressure (BHP) can be obtained by using Stehfest numerical inversion [30] to convert p~wDs to pwDtD.


## 4. Numerical Solution

### 4.1. Finite Difference Form of Dimensionless Mathematical Model

Though the analytical solution has been obtained in the previous section, we still seek to get the numerical solution of this model with a finite difference method [31], which can help us to verify the accuracy of the analytical solution.As shown in Figure2, we divide the reservoir into I segments in the radial direction, in which the inner region contains J segments of equal length and the outer region contains I-J segments of equal length. To ensure the computational efficiency and accuracy, the segment length of the inner region is shorter than that of the outer region. Meanwhile, in real space, time is discretized to obtain transient bottomhole pressure (BHP).Figure 2
Composite reservoir model divided intoI segments.


#### 4.1.1. Finite Difference Model of Region I

According to the uniform block-centered grid system shown in Figure2, we can differentiate the governing equation (31) at segment i at a timestep n as(51)pi-1n+1-2pin+1+pi+1n+1Δ2rD=pin+1-pinΔtDwhich can be rearranged as(52)pi-1n+1-a1pin+1+pi+1n+1=-a1-2pin,where a1=α(Δ2rD/ΔtD)+2 and pin represents the dimensionless pressure at the ith segment at a timestep n in the inner region. Equation (52) applies at all segments in the inner region, i=2,3,…,J-1, except the first and last segments of the inner region.Douglas and Peaceman [32] suggest treating no-flow boundaries by using a mirror image method. In the discretized system, taking the segment located at the wellbore (i=1) as a mirror, a pseudo segment (i=0) corresponding to the second segment (i=2) can be obtained on the other side of the wellbore (shown in Figure 2). p0n, p1n, and p2n represent the dimensionless pressure at the pseudo segment (i=0), the first segment (i=1), and the second segment (i=2) at a timestep n, respectively. Then, the inner boundary condition (see (32)) can be given as(53)∂p1D∂rDrD→0=p2n+1-p0n+12ΔrD=-β.Rearranging (53) gives(54)p2n+1-p0n+1=-2βΔrD.Applying (52) at the inner boundary of segment i=1, we obtain(55)p0n+1-a1p1n+1+p2n+1=-a1-2p1n.Then, adding (54) and (55) together gives(56)-a1p1n+1+2p2n+1=-b1+a1-2p1n,where b1=2βΔrD. Equation (56) is the differentiated form of flow equation for the first segment (i=1), which also corresponds to the inner boundary condition.


#### 4.1.2. Finite Difference Model of Region II

Similarly, the dimensionless governing equations of natural fracture system, vug system, and matrix system in the outer region can be differentiated as(57)fi-1n+1-2fin+1+fi+1n+1ΔrD12+1rDifi+1n+1-fi-1n+12ΔrD1=fin+1-finΔtD+ω2mω2fmin+1-minΔtD+ω2vω2fvin+1-vinΔtD,(58)-λmvvin+1-min+1-λvfvin+1-fin+1=ω2vω2fvin+1-vinΔtD,(59)λmffin+1-min+1+λmvvin+1-min+1=ω2mω2fmin+1-minΔtD, respectively. fin, vin, and min represent the dimensionless pressure of fracture, vug, and matrix at the ith segment at a timestep n, respectively.Combining (58) and (59), the following equations can be derived:(60)vin+1=Bλvf+λmvλmfAB-λmv2fin+1+Bω2vω2fΔtDAB-λmv2vin+λmvω2mω2fΔtDAB-λmv2min,(61)min+1=Aλmf+λmvλvfAB-λmv2fin+1+λmvω2vω2fΔtDAB-λmv2vin+Aω2mω2fΔtDAB-λmv2min,where(62)A=λmv+λvf+ω2vω2fΔtD;B=λmf+λmv+ω2mω2fΔtD.Then, substituting (60) and (61) into (57) gives(63)a2ifi-1n+1-a3fin+1+a4ifi+1n+1=-a5fin+a6min+a7vin,where(64)a2i=1-ΔrD121+i-J-1ΔrD1i=J+1,J+2,…,I;a3=2+ΔrD12ΔtD+ΔrD12ω2fΔtDω2mAλmf+λmvλvf+ω2vBλvf+λmvλmfAB-λmv2;a4i=1+ΔrD121+i-J-1ΔrD1i=J+1,J+2,…,I,rDi=1+i-J-1ΔrD1;a5=ΔrD12ΔtD;a6=ω2mΔrD12ω2fΔtDAω2m+λmvω2vω2fΔtDAB-λmv2-1;a7=ω2vΔrD12ω2fΔtDBω2v+λmvω2mω2fΔtDAB-λmv2-1.Equation (63) applies at all segments in the outer region, i=J+2,J+3,…,I-1, except the first and last segments of the outer region.Similarly, using the mirror image method at the external boundary, a pseudo segment (i=I+1) corresponding to the segment (i=I-1) can be obtained at the outside of the outer region (shown in Figure 2). fI-1n, fIn, and fI+1n represent the dimensionless fracture pressure at segment I-1, I, and I+1 at a timestep n, respectively. Then, the external boundary condition (see (38)) can be given as(65)∂p2D∂rDrD=reD=fI+1n-fI-1n2ΔrD1=0which can be rearranged as(66)fI+1n-fI-1n=0.Equation (63) applies at the external boundary of segment i=I, and the following equation is set up:(67)a2IfI-1n+1-a3fIn+1+a4IfI+1n+1=-a5fIn+a6mIn+a7vIn.Combining (66) and (67), the following equation can be obtained:(68)a2I+a4IfI-1n+1-a3fIn+1=-a5fIn+a6mIn+a7vIn.Equation (68) is the differentiated form of flow equation for the last segment (i=I) in the outer region.


#### 4.1.3. Finite Difference Treatment of Interface Condition

Applying (52) for the segment (i=J) and (63) for the segment (i=J+1) gives(69)pJ-1n+1-a1pJn+1+pJ+1n+1=-a1-2pJn,(70)a2J+1fJn+1-a3fJ+1n+1+a4J+1fJ+2n+1=-a5fJ+1n+a6mJ+1n+a7vJ+1n.The differentiated form of the continuity of pressure at the interface (see (39)) can be easily got:(71)pJ+1n+1=fJ+1n+1.Combining (71) and (69) will yield(72)pJ-1n+1-a1pJn+1+fJ+1n+1=-a1-2pJn.Equation (72) is the differentiated form of flow equation for the segment (i=J) at the inside of the interface, which just corresponds to the last segment of the inner region.The equation of the continuity of flux at the interface (see (40)) can be given as(73)pJ+1n+1-pJ-1n+12ΔrD=βfJ+2n+1-fJn+12ΔrD1which can be rearranged as(74)pJ+1n+1-pJ-1n+1=βΔrDΔrD1fJ+2n+1-fJn+1.Combining (70), (72), and (74) gives(75)-a1pJn+1+a9fJ+1n+1-a10fJ+2n+1=a11fJ+1n-a12mJ+1n-a13vJ+1n-a1-2pJn,where(76)a8=βΔrDΔrD1;a9=2+a3a8a2J+1;a10=a8a4J+1a2J+1+1;a11=a5a8a2J+1;a12=a6a8a2J+1;a13=a7a8a2J+1.Equation (75) is the differentiated form of flow equation for the segment (i=J+1) at the outside of the interface, corresponding to the first segment of the outer region.


### 4.2. Numerical Solution of Dimensionless Mathematical Model

Using the differentiated form for all the segments (i=1,2,3,…,J,J+1,…,I), a tridiagonal matrix with the rank of I×I is formed, which is a linear system:(77)-a121-a11⋯⋯⋯1-a11-a1a9-a10a2J+2-a3a4J+2⋯⋯⋯a2I+a4I-a3p1n+1p2n+1⋯pJn+1fJ+1n+1fJ+2n+1⋯fIn+1=-a1-2p1n-b1-a1-2p2n⋯-a1-2pJna11fJ+1n-a12mJ+1n-a13vJ+1n-a1-2pJn-a5fJ+2n+a6mJ+2n+a7vJ+2n⋯-a5fIn+a6mIn+a7vIn.The initial conditions in the differentiated form for the inner and outer region are given as(78)pi0=0i=1,2,…,J(79)fi0=vi0=mi0=0i=J+1,J+2,…,I.p 1 n represents the wellbore pressure (BHP) at a timestep n. An iterative method is used to get the wellbore pressure at any timestep; thus, transient responses corresponding to different time can be obtained. In this paper, Thomas chasing method, one of the classical iterative methods, is used to solve the tridiagonal matrix. It is found that the solution does not change appreciably when the reservoir with a dimensionless radius of 100 is divided into 4000 segments. In all cases studied, 10 time intervals are taken in each log cycle of dimensionless time.


## 4.1. Finite Difference Form of Dimensionless Mathematical Model

Though the analytical solution has been obtained in the previous section, we still seek to get the numerical solution of this model with a finite difference method [31], which can help us to verify the accuracy of the analytical solution.As shown in Figure2, we divide the reservoir into I segments in the radial direction, in which the inner region contains J segments of equal length and the outer region contains I-J segments of equal length. To ensure the computational efficiency and accuracy, the segment length of the inner region is shorter than that of the outer region. Meanwhile, in real space, time is discretized to obtain transient bottomhole pressure (BHP).Figure 2
Composite reservoir model divided intoI segments.


### 4.1.1. Finite Difference Model of Region I

According to the uniform block-centered grid system shown in Figure2, we can differentiate the governing equation (31) at segment i at a timestep n as(51)pi-1n+1-2pin+1+pi+1n+1Δ2rD=pin+1-pinΔtDwhich can be rearranged as(52)pi-1n+1-a1pin+1+pi+1n+1=-a1-2pin,where a1=α(Δ2rD/ΔtD)+2 and pin represents the dimensionless pressure at the ith segment at a timestep n in the inner region. Equation (52) applies at all segments in the inner region, i=2,3,…,J-1, except the first and last segments of the inner region.Douglas and Peaceman [32] suggest treating no-flow boundaries by using a mirror image method. In the discretized system, taking the segment located at the wellbore (i=1) as a mirror, a pseudo segment (i=0) corresponding to the second segment (i=2) can be obtained on the other side of the wellbore (shown in Figure 2). p0n, p1n, and p2n represent the dimensionless pressure at the pseudo segment (i=0), the first segment (i=1), and the second segment (i=2) at a timestep n, respectively. Then, the inner boundary condition (see (32)) can be given as(53)∂p1D∂rDrD→0=p2n+1-p0n+12ΔrD=-β.Rearranging (53) gives(54)p2n+1-p0n+1=-2βΔrD.Applying (52) at the inner boundary of segment i=1, we obtain(55)p0n+1-a1p1n+1+p2n+1=-a1-2p1n.Then, adding (54) and (55) together gives(56)-a1p1n+1+2p2n+1=-b1+a1-2p1n,where b1=2βΔrD. Equation (56) is the differentiated form of flow equation for the first segment (i=1), which also corresponds to the inner boundary condition.


### 4.1.2. Finite Difference Model of Region II

Similarly, the dimensionless governing equations of natural fracture system, vug system, and matrix system in the outer region can be differentiated as(57)fi-1n+1-2fin+1+fi+1n+1ΔrD12+1rDifi+1n+1-fi-1n+12ΔrD1=fin+1-finΔtD+ω2mω2fmin+1-minΔtD+ω2vω2fvin+1-vinΔtD,(58)-λmvvin+1-min+1-λvfvin+1-fin+1=ω2vω2fvin+1-vinΔtD,(59)λmffin+1-min+1+λmvvin+1-min+1=ω2mω2fmin+1-minΔtD, respectively. fin, vin, and min represent the dimensionless pressure of fracture, vug, and matrix at the ith segment at a timestep n, respectively.Combining (58) and (59), the following equations can be derived:(60)vin+1=Bλvf+λmvλmfAB-λmv2fin+1+Bω2vω2fΔtDAB-λmv2vin+λmvω2mω2fΔtDAB-λmv2min,(61)min+1=Aλmf+λmvλvfAB-λmv2fin+1+λmvω2vω2fΔtDAB-λmv2vin+Aω2mω2fΔtDAB-λmv2min,where(62)A=λmv+λvf+ω2vω2fΔtD;B=λmf+λmv+ω2mω2fΔtD.Then, substituting (60) and (61) into (57) gives(63)a2ifi-1n+1-a3fin+1+a4ifi+1n+1=-a5fin+a6min+a7vin,where(64)a2i=1-ΔrD121+i-J-1ΔrD1i=J+1,J+2,…,I;a3=2+ΔrD12ΔtD+ΔrD12ω2fΔtDω2mAλmf+λmvλvf+ω2vBλvf+λmvλmfAB-λmv2;a4i=1+ΔrD121+i-J-1ΔrD1i=J+1,J+2,…,I,rDi=1+i-J-1ΔrD1;a5=ΔrD12ΔtD;a6=ω2mΔrD12ω2fΔtDAω2m+λmvω2vω2fΔtDAB-λmv2-1;a7=ω2vΔrD12ω2fΔtDBω2v+λmvω2mω2fΔtDAB-λmv2-1.Equation (63) applies at all segments in the outer region, i=J+2,J+3,…,I-1, except the first and last segments of the outer region.Similarly, using the mirror image method at the external boundary, a pseudo segment (i=I+1) corresponding to the segment (i=I-1) can be obtained at the outside of the outer region (shown in Figure 2). fI-1n, fIn, and fI+1n represent the dimensionless fracture pressure at segment I-1, I, and I+1 at a timestep n, respectively. Then, the external boundary condition (see (38)) can be given as(65)∂p2D∂rDrD=reD=fI+1n-fI-1n2ΔrD1=0which can be rearranged as(66)fI+1n-fI-1n=0.Equation (63) applies at the external boundary of segment i=I, and the following equation is set up:(67)a2IfI-1n+1-a3fIn+1+a4IfI+1n+1=-a5fIn+a6mIn+a7vIn.Combining (66) and (67), the following equation can be obtained:(68)a2I+a4IfI-1n+1-a3fIn+1=-a5fIn+a6mIn+a7vIn.Equation (68) is the differentiated form of flow equation for the last segment (i=I) in the outer region.


### 4.1.3. Finite Difference Treatment of Interface Condition

Applying (52) for the segment (i=J) and (63) for the segment (i=J+1) gives(69)pJ-1n+1-a1pJn+1+pJ+1n+1=-a1-2pJn,(70)a2J+1fJn+1-a3fJ+1n+1+a4J+1fJ+2n+1=-a5fJ+1n+a6mJ+1n+a7vJ+1n.The differentiated form of the continuity of pressure at the interface (see (39)) can be easily got:(71)pJ+1n+1=fJ+1n+1.Combining (71) and (69) will yield(72)pJ-1n+1-a1pJn+1+fJ+1n+1=-a1-2pJn.Equation (72) is the differentiated form of flow equation for the segment (i=J) at the inside of the interface, which just corresponds to the last segment of the inner region.The equation of the continuity of flux at the interface (see (40)) can be given as(73)pJ+1n+1-pJ-1n+12ΔrD=βfJ+2n+1-fJn+12ΔrD1which can be rearranged as(74)pJ+1n+1-pJ-1n+1=βΔrDΔrD1fJ+2n+1-fJn+1.Combining (70), (72), and (74) gives(75)-a1pJn+1+a9fJ+1n+1-a10fJ+2n+1=a11fJ+1n-a12mJ+1n-a13vJ+1n-a1-2pJn,where(76)a8=βΔrDΔrD1;a9=2+a3a8a2J+1;a10=a8a4J+1a2J+1+1;a11=a5a8a2J+1;a12=a6a8a2J+1;a13=a7a8a2J+1.Equation (75) is the differentiated form of flow equation for the segment (i=J+1) at the outside of the interface, corresponding to the first segment of the outer region.


## 4.1.1. Finite Difference Model of Region I

According to the uniform block-centered grid system shown in Figure2, we can differentiate the governing equation (31) at segment i at a timestep n as(51)pi-1n+1-2pin+1+pi+1n+1Δ2rD=pin+1-pinΔtDwhich can be rearranged as(52)pi-1n+1-a1pin+1+pi+1n+1=-a1-2pin,where a1=α(Δ2rD/ΔtD)+2 and pin represents the dimensionless pressure at the ith segment at a timestep n in the inner region. Equation (52) applies at all segments in the inner region, i=2,3,…,J-1, except the first and last segments of the inner region.Douglas and Peaceman [32] suggest treating no-flow boundaries by using a mirror image method. In the discretized system, taking the segment located at the wellbore (i=1) as a mirror, a pseudo segment (i=0) corresponding to the second segment (i=2) can be obtained on the other side of the wellbore (shown in Figure 2). p0n, p1n, and p2n represent the dimensionless pressure at the pseudo segment (i=0), the first segment (i=1), and the second segment (i=2) at a timestep n, respectively. Then, the inner boundary condition (see (32)) can be given as(53)∂p1D∂rDrD→0=p2n+1-p0n+12ΔrD=-β.Rearranging (53) gives(54)p2n+1-p0n+1=-2βΔrD.Applying (52) at the inner boundary of segment i=1, we obtain(55)p0n+1-a1p1n+1+p2n+1=-a1-2p1n.Then, adding (54) and (55) together gives(56)-a1p1n+1+2p2n+1=-b1+a1-2p1n,where b1=2βΔrD. Equation (56) is the differentiated form of flow equation for the first segment (i=1), which also corresponds to the inner boundary condition.


## 4.1.2. Finite Difference Model of Region II

Similarly, the dimensionless governing equations of natural fracture system, vug system, and matrix system in the outer region can be differentiated as(57)fi-1n+1-2fin+1+fi+1n+1ΔrD12+1rDifi+1n+1-fi-1n+12ΔrD1=fin+1-finΔtD+ω2mω2fmin+1-minΔtD+ω2vω2fvin+1-vinΔtD,(58)-λmvvin+1-min+1-λvfvin+1-fin+1=ω2vω2fvin+1-vinΔtD,(59)λmffin+1-min+1+λmvvin+1-min+1=ω2mω2fmin+1-minΔtD, respectively. fin, vin, and min represent the dimensionless pressure of fracture, vug, and matrix at the ith segment at a timestep n, respectively.Combining (58) and (59), the following equations can be derived:(60)vin+1=Bλvf+λmvλmfAB-λmv2fin+1+Bω2vω2fΔtDAB-λmv2vin+λmvω2mω2fΔtDAB-λmv2min,(61)min+1=Aλmf+λmvλvfAB-λmv2fin+1+λmvω2vω2fΔtDAB-λmv2vin+Aω2mω2fΔtDAB-λmv2min,where(62)A=λmv+λvf+ω2vω2fΔtD;B=λmf+λmv+ω2mω2fΔtD.Then, substituting (60) and (61) into (57) gives(63)a2ifi-1n+1-a3fin+1+a4ifi+1n+1=-a5fin+a6min+a7vin,where(64)a2i=1-ΔrD121+i-J-1ΔrD1i=J+1,J+2,…,I;a3=2+ΔrD12ΔtD+ΔrD12ω2fΔtDω2mAλmf+λmvλvf+ω2vBλvf+λmvλmfAB-λmv2;a4i=1+ΔrD121+i-J-1ΔrD1i=J+1,J+2,…,I,rDi=1+i-J-1ΔrD1;a5=ΔrD12ΔtD;a6=ω2mΔrD12ω2fΔtDAω2m+λmvω2vω2fΔtDAB-λmv2-1;a7=ω2vΔrD12ω2fΔtDBω2v+λmvω2mω2fΔtDAB-λmv2-1.Equation (63) applies at all segments in the outer region, i=J+2,J+3,…,I-1, except the first and last segments of the outer region.Similarly, using the mirror image method at the external boundary, a pseudo segment (i=I+1) corresponding to the segment (i=I-1) can be obtained at the outside of the outer region (shown in Figure 2). fI-1n, fIn, and fI+1n represent the dimensionless fracture pressure at segment I-1, I, and I+1 at a timestep n, respectively. Then, the external boundary condition (see (38)) can be given as(65)∂p2D∂rDrD=reD=fI+1n-fI-1n2ΔrD1=0which can be rearranged as(66)fI+1n-fI-1n=0.Equation (63) applies at the external boundary of segment i=I, and the following equation is set up:(67)a2IfI-1n+1-a3fIn+1+a4IfI+1n+1=-a5fIn+a6mIn+a7vIn.Combining (66) and (67), the following equation can be obtained:(68)a2I+a4IfI-1n+1-a3fIn+1=-a5fIn+a6mIn+a7vIn.Equation (68) is the differentiated form of flow equation for the last segment (i=I) in the outer region.


## 4.1.3. Finite Difference Treatment of Interface Condition

Applying (52) for the segment (i=J) and (63) for the segment (i=J+1) gives(69)pJ-1n+1-a1pJn+1+pJ+1n+1=-a1-2pJn,(70)a2J+1fJn+1-a3fJ+1n+1+a4J+1fJ+2n+1=-a5fJ+1n+a6mJ+1n+a7vJ+1n.The differentiated form of the continuity of pressure at the interface (see (39)) can be easily got:(71)pJ+1n+1=fJ+1n+1.Combining (71) and (69) will yield(72)pJ-1n+1-a1pJn+1+fJ+1n+1=-a1-2pJn.Equation (72) is the differentiated form of flow equation for the segment (i=J) at the inside of the interface, which just corresponds to the last segment of the inner region.The equation of the continuity of flux at the interface (see (40)) can be given as(73)pJ+1n+1-pJ-1n+12ΔrD=βfJ+2n+1-fJn+12ΔrD1which can be rearranged as(74)pJ+1n+1-pJ-1n+1=βΔrDΔrD1fJ+2n+1-fJn+1.Combining (70), (72), and (74) gives(75)-a1pJn+1+a9fJ+1n+1-a10fJ+2n+1=a11fJ+1n-a12mJ+1n-a13vJ+1n-a1-2pJn,where(76)a8=βΔrDΔrD1;a9=2+a3a8a2J+1;a10=a8a4J+1a2J+1+1;a11=a5a8a2J+1;a12=a6a8a2J+1;a13=a7a8a2J+1.Equation (75) is the differentiated form of flow equation for the segment (i=J+1) at the outside of the interface, corresponding to the first segment of the outer region.


## 4.2. Numerical Solution of Dimensionless Mathematical Model

Using the differentiated form for all the segments (i=1,2,3,…,J,J+1,…,I), a tridiagonal matrix with the rank of I×I is formed, which is a linear system:(77)-a121-a11⋯⋯⋯1-a11-a1a9-a10a2J+2-a3a4J+2⋯⋯⋯a2I+a4I-a3p1n+1p2n+1⋯pJn+1fJ+1n+1fJ+2n+1⋯fIn+1=-a1-2p1n-b1-a1-2p2n⋯-a1-2pJna11fJ+1n-a12mJ+1n-a13vJ+1n-a1-2pJn-a5fJ+2n+a6mJ+2n+a7vJ+2n⋯-a5fIn+a6mIn+a7vIn.The initial conditions in the differentiated form for the inner and outer region are given as(78)pi0=0i=1,2,…,J(79)fi0=vi0=mi0=0i=J+1,J+2,…,I.p 1 n represents the wellbore pressure (BHP) at a timestep n. An iterative method is used to get the wellbore pressure at any timestep; thus, transient responses corresponding to different time can be obtained. In this paper, Thomas chasing method, one of the classical iterative methods, is used to solve the tridiagonal matrix. It is found that the solution does not change appreciably when the reservoir with a dimensionless radius of 100 is divided into 4000 segments. In all cases studied, 10 time intervals are taken in each log cycle of dimensionless time.


## 5. Model Verification

The model in this study can be verified by comparing the results of its analytical and numerical solution. Figure3 shows the comparison results of the analytical and numerical solutions at different parameters. All the analytical solutions agree well with the corresponding numerical solutions within the most production life. At early time and some of intermediate time there is a slight discrepancy between the two kinds of solutions on pressure derivative curves, which may be caused by the long length of discrete segments, but the error is within an acceptable limit. These validations indicate that the analytical model in this work is reliable. In addition to acting as an effective way to verify the accuracy of the analytical solution, the numerical solution in this study has some advantages over the analytical solution in computing speed. However, the analytical solution also has its unique advantages in some respects: for example, wellbore storage and skin effect can be easily added to the analytical solution in Laplace domain [33–35], which may be extremely difficult for the numerical solution.Figure 3
Comparison results of analytical and numerical solutions at different parameters.
(a)
α = 0.005, 0.5 (b)
β = 0.005, 0.5 (c)
ω 2 f = 0.002, 0.5 (d)
ω 2 m = 0.6, 0.9 (e)
ω 2 v = 0.3, 0.7 (f)
λ m f = 10 - 4, 10−2  (g)
λ m v = 0.0001, 0.05 (h)
λ v f = 0.05, 0.1


## 6. Results and Discussion

One objective of this research is to calculate type curves for various reservoirs and observe the flow behavior from these curves. By type curves matching, some reservoir parameters, such as permeability, skin factor, and reservoir drainage area, can be obtained [36, 37]. In this section, we discuss the characteristics of transient response of an acid fracturing well in a FVCR and also analyze the differences between our model and the similar model with single fracture or without acid fracturing and further investigate the pseudo-skin factor caused by acid fracturing.


### 6.1. Flow Regime Recognition

In the carbonate reservoirs, either matrix system or vug system provides the major storage space; thus, these reservoirs can be subdivided into two categories: reserves dominated by matrix system and reserves dominated by vug system. Figures4 and 5 show the complete transient flow processes in these two categories. The values of the parameters used in Figures 4 and 5 are listed in Tables 1 and 2, respectively.Table 1
The data used in Figure4.
Model parameter Value Model parameter Value α 0.005 r w D 0.001 β 0.5 r e D 70 ω2f 0.001 λ m f 0.001 ω2m 0.96 λ m v 0.01 ω2v 0.039 λ v f 0.1Table 2
The data used in Figure5.
Model parameter Value Model parameter Value α 0.005 r w D 0.001 β 0.5 r e D 70 ω 2 f 0.001 λ m f 0.001 ω 2 m 0.249 λ m v 0.01 ω 2 v 0.75 λ v f 0.1Figure 4
Flow regimes on dimensionless pressure and pressure derivative curves (reserves dominated by matrix system).Figure 5
Flow regimes on dimensionless pressure and pressure derivative curves (reserves dominated by vug system).


#### 6.1.1. Case 1: Hydrocarbon Reserves Dominated by Matrix System

As is described in Figure4, the flow process in this reservoir can be divided into eight flow regimes. Initially, there is a linear flow period characterized by a half-slope straight line on the pressure derivative curve. After a transition flow regime, the first interporosity flow develops and is followed with the second transition flow. Then, the second interporosity flow occurs and immediately is followed by the third transition flow. Eventually, a short radial flow regime may be observed before the system reaches a pseudo-steady flow. Interestingly, three characteristic V-shaped segments can be observed on the pressure derivative curve.(  1) Linear Flow Regime. In this regime, the production is dominated by the fluid stored in artificial fractures and wormholes. It occurs at very small values of dimensionless time. The pressure and pressure derivative curves are an upward straight line with a slope of 0.5. This regime is typical of fluid flow for acid fracturing wells. Before the interface is felt, the wellbore pressure solution in this study can be approximate to(80)pwDtD≈2βαtDπ.Corresponding to the square root of the dimensionless time, this equation implies that a log-log plot of function results a straight line with the slope of one-half, which is the typical characteristic of linear flow. This equation also indicates that the wellbore pressure in this regime depends on the values of interregional diffusivity ratio and interregional conductivity ratio (seen in Figures 6 and 7).Figure 6
Effect of interregional diffusivity ratio on transient flow behavior.Figure 7
Effect of interregional conductivity ratio on transient flow behavior.(  2) Transition Flow Regime (a). The first transition flow takes place when the interface takes effect. The derivative curve in this regime shows the first V-shaped segment, which represents fluid flow from natural fractures into artificial fractures or wormholes. It can be interpreted by the fact that as the pressure of the inner region depletes gradually, natural fracture system responds first owing to the direct connection to the inner region. The shape of this segment is also controlled by interregional diffusivity ratio and interregional conductivity ratio, and this V-shaped segment may not be observed at particular parameters (seen in Figures 6 and 7).(  3) Interporosity Flow Regime (a). The first interporosity flow regime corresponds to the second V-shaped segment on the derivative curve, which represents the crossflow between vug system and natural fracture system. This interporosity flow takes place first because vug permeability is larger than matrix permeability. This segment is predominantly determined by storativity ratios of natural fracture system and matrix system, and interporosity coefficient between vugs and natural fractures, and also may be absent from the curves at particular parameters (seen in Figures 8–10, and 13).Figure 8
Effect of fracture storativity ratio on transient flow behavior.Figure 9
Effect of matrix storativity ratio on transient flow behavior.Figure 10
Effect of vug storativity ratio on transient flow behavior.(  4) Transition Flow Regime (b). With the depletion of vug system, the matrix system begins to supply fluid to natural fracture system gradually and a short transition flow occurs. The duration of this regime depends on the difference of the ability supplying fluid between matrix system and vug system; thus, this regime may be so short that it can hardly be identified from the derivative curve.(  5) Interporosity Flow Regime (b). In this interporosity flow regime, the production is mainly dominated by the fluid flow from matrix system to fracture system. This regime is characterized by the third V-shaped segment on the derivative curve, which is governed by two physical processes: one is interporosity flow from matrix system to fracture system and the other is interporosity flow from vug system to matrix system. The former has greater influence on the shape of this regime than the latter (seen in Figures 11 and 12).Figure 11
Effect of matrix-fracture interporosity coefficient on transient flow behavior.Figure 12
Effect of matrix-vug interporosity coefficient on transient flow behavior.Figure 13
Effect of vug-fracture interporosity coefficient on transient flow behavior.(  6) Transition Flow Regime (c). With the depletion of matrix system, fluids far away from the wellbore begin to flow. Like the second transition flow, another short transition flow regime may develop before the dynamic balance reaches the whole reservoir. It also may not be observed on the pressure derivative curve.(  7) Radial Flow Regime. When the matrix system depletes fast and the reservoir radius is large enough, fluid flow in the whole reservoir is likely to reach a dynamic balance; thus, a radial flow occurs after the third transition flow. In this regime, the slope of pressure derivative curve is zero, and the pressure derivative converges to “0.5 lines.” As shown in Figure 14, reservoir radius has a great influence on the duration of radial flow.Figure 14
Effect of dimensionless reservoir radius on transient flow behavior.(  8) Pseudo-Steady Flow Regime. For closed boundary, the long-time behavior corresponds to a pseudo-steady flow. Both the pressure and derivative curves show a unit slope straight line, which is the typical characteristic of closed reservoirs. By means of asymptotic analysis, the approximate solution in this regime can be yielded:(81)pwDtD≈2ω2ftDreD2-1+reD22reD2-1ln⁡ω2ftD+reD2ln⁡4-0.5772+12reD2-1.This solution shows that wellbore pressure in pseudo-steady regime only depends on natural fracture storativity ratio and reservoir radius (seen in Figures6–14).


#### 6.1.2. Case 2: Hydrocarbon Reserves Dominated by Vug System

As is shown in Figure5, flow in this reservoir can be divided into six flow regimes: linear flow, transition flow (a), interporosity flow, transition flow (b), radial flow, and pseudo-steady flow. The major difference of flow characteristics between two kinds of carbonate reservoirs lies in the interporosity flow period. The characteristics of the other flow regimes are the same with the corresponding regimes of the previous case, including the approximate solutions of linear flow and pseudo-steady flow.Because the ability of vug system supplying fluid to natural fracture system is much stronger than that of matrix system supplying fluid to natural fracture system, the last two V-shaped segments are connected to each other and form one V-shaped segment. Consequently, only two V-shaped segments can be observed on the pressure derivative curve and the second V-shaped segment merges with the radial flow.


### 6.2. Parameter Sensitivity Analysis

The shape rendered from real data may be distorted by noises in well testing, which makes it necessary to establish the stylized shapes under different parameter conditions [36]. Figures 6–14 contain the type curves of pressure-transient response for an acid fracturing well in FVCRs. Later we will show how the shape of the type curves may be affected by changing the parameters.


#### 6.2.1. Effect of Interregional Diffusivity Ratio

According to the dimensionless definition, interregional diffusivity ratioα represents the ratio of fluid transmissivity ability between natural fractures in the outer region and artificial fractures or wormholes in the inner region. A larger α means a weaker diffusion ability of wormholes or artificial fractures. Figure 6 shows the effect of this parameter on pressure-transient response (α=0.001, 0.005, 0.05, and 0.5). We can find that interregional diffusivity ratio mainly affects early-flow period. Before the first interporosity flow regime, a larger α leads to lower pressure depletion at the same production rate. Meanwhile, as this parameter increases, due to a weaker diffusion ability in the inner region, duration of the linear flow is prolonged, while duration of the first interporosity flow is shortened. When interregional diffusivity ratio gradually increases to 1, the first V-shaped segment becomes shallower and even disappears. All the curves for different interregional diffusivity ratios begin to normalize at the end of the first interporosity flow regime.


#### 6.2.2. Effect of Interregional Conductivity Ratio

According to the dimensionless definition, interregional conductivity ratioβ represents a measure of the conductivity between wormholes or artificial fractures and natural fractures. In practice, a larger β means a lower conductivity of the inner region. Transient responses on the log-log plots for varied interregional conductivity ratio are shown in Figure 7 (β=0.2, 0.4, 0.6, and 0.8). It can be concluded that this parameter mainly affects early- and middle-flow periods. A larger interregional conductivity ratio leads to more pressure depletion at the same production. Meanwhile, with the increase of this parameter, the hump corresponding to the first transition flow becomes more significant. Besides, the pressure curves for varied interregional conductivity ratio normalize in the late-flow period, while the pressure derivative curves begin to normalize from the first interporosity flow regime.


#### 6.2.3. Effect of Storativity Ratios

There are three storativity ratios in a FVCR. Each ratio represents the relative capacity of fluid stored in a corresponding system. Generally fluid stored in the natural fracture system is much less than that stored in the vug or matrix system.Figure8 presents the results of sensitivity analysis of fracture storativity ratio (ω2f = 0.001, 0.005, 0.025, and 0.125). It indicates that a larger ω2f leads to more pressure depletion in the middle- and late-flow periods. Fracture storativity ratio begins to affect the transient response at the start of the first interporosity flow, which extends to the pseudo-state flow period; thus, the curves do not normalize, but they are parallel to each other in the late-flow period. Due to a higher conductivity and diffusion capacity for natural fracture system provided by a larger fracture storativity ratio, when this ratio increases from 0.001 to 0.125, the last two V-shaped segments become shallower and even the second V-shaped segment may disappear. Moreover, the pseudo-steady flow begins earlier with the increase of this parameter.When investigating the effect of another two storativity ratios, we keepω2f constant; thus, matrix storativity ratio has a reciprocal relationship with vug storativity ratio. Therefore, it is not necessary to discuss all the values in their own range.Figure9 shows the response of varying matrix storativity ratio, ω2m = 0.6, 0.7, 0.8, and 0.9, which indicates that matrix system provides the majority of storage space. As is described in Figure 9, matrix storativity ratio mainly affects the middle-flow period, and the pressure and derivative curves for different matrix storativity ratios overlap completely in the early- and late-flow periods. But the shape of the last two V-shaped segments is affected by this storativity ratio significantly. With the increase of matrix storativity ratio, the second V-shaped segment becomes shallower, while the third V-shaped segment becomes deeper. In fact, matrix storativity ratio only affects the third V-shaped segment and vug storativity ratio only affects the second V-shaped segment. However, due to the reciprocal relationship between two storativity ratios, matrix storativity ratio affects the second V-shaped segment indirectly.Figure10 presents the results of sensitivity analysis of vug storativity ratio, ω2v = 0.6, 0.7, 0.8, and 0.9, which indicates that the reserves are dominated by vug system. As shown in Figure 10, vug storativity ratio only has a slight effect on the middle-flow period and the related curves normalize in the late-flow period. With the increase of vug storativity ratio, the concave of the second V-shaped segment deepens slightly.


#### 6.2.4. Effect of Interporosity Coefficients

There are also three interporosity coefficients in the triple-porosity reservoir. A sensitivity analysis is also performed by varying these three coefficients.Figure11 demonstrates the effect of matrix-fracture interporosity coefficient on the transient behavior (λmf = 5 × 10−5, 5 × 10−4, 5 × 10−3, and 10−2). It can be seen that this parameter mainly affects the middle-flow period and type curves for different interporosity coefficients normalize in the late-flow period. As this coefficient increases, the second V-shaped segment deepens, the third V-shaped segment moves to the left side, and the duration of radial flow regime is prolonged. This is mainly because a larger matrix-fracture interporosity coefficient means a stronger ability of matrix system supplying fluid for natural fracture system, which makes the arrival time of matrix-fracture interporosity flow regime earlier.Figure12 presents the results of sensitivity analysis of matrix-vug interporosity coefficient (λmv = 0.001, 0.005, 0.025, and 0.05). This coefficient also mainly affects middle-flow period, especially the matrix-fracture interporosity flow regime, corresponding to the third V-shaped segment. It can be interpreted by that this coefficient only affects the ability of vug system to supply fluid for matrix system. With the increase of this coefficient, the starting time of the second interporosity flow moves back and the corresponding V-shaped segment deepens greatly. Similarly, the pressure and derivative curves will normalize in the late-flow period.Figure13 shows the response of varying vug-fracture interporosity coefficient (λvf = 0.01, 0.05, 0.25, and 0.5). This parameter mainly affects the middle-flow period, and the normalization can also be observed in the late-flow period. As this coefficient increases, the starting time of the first interporosity flow moves back and the third V-shaped segment becomes shallower. It can be interpreted by that the ability of vug system to supply fluid for fracture system becomes much stronger than that of matrix system supplying fluid for fracture system with the increase of this coefficient.


#### 6.2.5. Effect of Dimensionless Reservoir Radius

Figure14 presents the results of sensitivity analysis of dimensionless reservoir radius (reD = 30, 50, 70, and 90). There is no difference in the shape of type curves in the early- and middle-flow periods under different reservoir radius; however, the pressure and pressure derivative curves show a unit slope straight line and go up rapidly in the late-flow period. In addition, with the increase of the dimensionless reservoir radius, the duration of radial flow regime is prolonged significantly.


### 6.3. Comparison with Other Similar Models

Gringarten et al. [38] investigated unsteady-state pressure distributions created by a well with a single infinite-conductivity vertical fracture in a homogeneous reservoir. Further, we can easily obtain the transient-pressure solution for the same well type in a FVCR. There are a finite number of artificial fractures and wormholes in our model’s inner region; thus, it is necessary to compare these two models’ difference in flow behavior. Meanwhile, it is also significant to make a comparison of our model with the model not having acid fracturing to know the effect of acidification. In order to have a convenient comparison, we firstly make a comparison to the parameters used in three models and then compare their type curves by fixing a set of the same parameter values. Table 3 shows the differences between three models’ parameters. It can be seen that the same parameters are ω2f, ω2m, ω2v, λmf, λmv, λvf, rwD, and reD, but the parameters α and β only belong to the model with acid fracturing. To make an effective comparison, we set the same values for the same parameters, which are shown in Table 3.Table 3
The comparison of different parameters between different models.
Model Model parameters Acid fracturing α β ω 2 f ω 2 m ω 2 v λ m f λ m v λ v f r w D r e D 0.001 0.05 0.001 0.96 0.039 0.001 0.01 0.1 0.001 40 Single-vertical fracture / / ω 2 f ω 2 m ω 2 v λ m f λ m v λ v f r w D r e D / / 0.001 0.96 0.039 0.001 0.01 0.1 0.001 40 No acid fracturing / / ω 2 f ω 2 m ω 2 v λ m f λ m v λ v f r w D r e D / / 0.001 0.96 0.039 0.001 0.01 0.1 0.001 40Figure15 shows the comparison result of type curves between three models, and we can observe the differences clearly: (i) in the early-flow period, the pressure and derivative curves for the model with acid fracturing and the model with single fracture coincide perfectly, respectively. It demonstrates that these two models’ flow behavior is the same in this period. However, because drainage area of the well without acid fracturing is much smaller than that of the well with acid fracturing or with single fracture, the pressure depletion of the former case is faster than that of the latter two cases at the same group of production parameters; (ii) in the middle-flow period, the pressure depletion of the well without single fracture becomes faster than that of the well with acid fracturing, which also can be interpreted by the difference of drainage area; (iii) in the models with single fracture and without acid fracturing, only two V-shaped segments can be observed on the derivative curves, and the first V-shaped segment in the model with acid fracturing is replaced by a unit slope straight line or a horizontal line with the value of 0.5, which are the typical characteristics of linear flow and radial flow, respectively.Figure 15
The comparison of type curves between three different models.


### 6.4. Parameter Influence on Pseudo-Skin Factor

During the pseudo-steady flow period, for the same position, there is a difference of the dimensionless pressure between the models with acid fracturing and without acid fracturing. This difference remains constant at large value of dimensionless time and can be handled as a pseudo-skin factor [39, 40] that can be defined as(82)SrD,tD;α,β,λ,ω,reD=p1DrD,tD;α,β,λvf,λmv,λmf,ω2f,ω2v,ω2m,reD-pNAFDrD,tD;λvf,λmv,λmf,ω2f,ω2v,ω2m,reD, where p1D and pNAFD represent the dimensionless pressure near the wellbore for the model with acid fracturing and without the stimulation, respectively.Applying the Laplace transform to (82), we have(83)S~rD,s;α,β,λ,ω,reD=p~1DrD,s;α,β,λvf,λmv,λmf,ω2f,ω2v,ω2m,reD-p~NAFDrD,s;λvf,λmv,λmf,ω2f,ω2v,ω2m,reD,where(84)p~1D=βssαx2sinh⁡sα1-rD-x1cosh⁡sα1-rDx2cosh⁡sα-x1sinh⁡sα,p~NAFD=1sK0rDsfs+K1reDsfsI1reDsfsI0rDsfs.Note that x1 and x2 have been given in (49).Because of higher connectivity between the wellbore and the reservoir caused by acid fracturing, the pseudo-skin factor defined here is always negative and indicates an increase in well productivity. The larger the absolute value of the pseudo-skin factor, the better the stimulated effect of acid fracturing.As is shown in Figure16, under different interregional diffusivity ratios, the pseudo-skin factors at the same position near the wellbore are the same, which indicates that this parameter has no influence on the ability of enhancing the productivity in the late-flow period. In the pseudo-steady flow regime, the diffusion capacity of the inner region well satisfies the demand of fluid flow feeding the wellbore and the flow in the model reaches a dynamic balance; thus, the change of interregional diffusivity ratio does not affect the productivity. It also can be concluded that the closer the position to the wellbore, the better the stimulation effect.Figure 16
Pseudo-skin factor versus dimensionless formation radius for an acid fracturing well at different interregional diffusivity ratios.Figure17 presents different pseudo-skin factors for the vertical well of varied interregional conductivity ratio. With the increase of interregional conductivity ratio, the absolute value of the pseudo-skin factor at the same position decreases, which implies that the ability of enhancing the well production decreases. As we know, a larger interregional conductivity ratio corresponds to a weaker conductivity of the inner region and also means the increase of the flow resistance, which explains the result above. Similarly, in this group, the closer the position to the wellbore, the better the stimulation effect.Figure 17
Pseudo-skin factor versus dimensionless reservoir radius for an acid fracturing well at different interregional conductivity ratios.


## 6.1. Flow Regime Recognition

In the carbonate reservoirs, either matrix system or vug system provides the major storage space; thus, these reservoirs can be subdivided into two categories: reserves dominated by matrix system and reserves dominated by vug system. Figures4 and 5 show the complete transient flow processes in these two categories. The values of the parameters used in Figures 4 and 5 are listed in Tables 1 and 2, respectively.Table 1
The data used in Figure4.
Model parameter Value Model parameter Value α 0.005 r w D 0.001 β 0.5 r e D 70 ω2f 0.001 λ m f 0.001 ω2m 0.96 λ m v 0.01 ω2v 0.039 λ v f 0.1Table 2
The data used in Figure5.
Model parameter Value Model parameter Value α 0.005 r w D 0.001 β 0.5 r e D 70 ω 2 f 0.001 λ m f 0.001 ω 2 m 0.249 λ m v 0.01 ω 2 v 0.75 λ v f 0.1Figure 4
Flow regimes on dimensionless pressure and pressure derivative curves (reserves dominated by matrix system).Figure 5
Flow regimes on dimensionless pressure and pressure derivative curves (reserves dominated by vug system).


### 6.1.1. Case 1: Hydrocarbon Reserves Dominated by Matrix System

As is described in Figure4, the flow process in this reservoir can be divided into eight flow regimes. Initially, there is a linear flow period characterized by a half-slope straight line on the pressure derivative curve. After a transition flow regime, the first interporosity flow develops and is followed with the second transition flow. Then, the second interporosity flow occurs and immediately is followed by the third transition flow. Eventually, a short radial flow regime may be observed before the system reaches a pseudo-steady flow. Interestingly, three characteristic V-shaped segments can be observed on the pressure derivative curve.(  1) Linear Flow Regime. In this regime, the production is dominated by the fluid stored in artificial fractures and wormholes. It occurs at very small values of dimensionless time. The pressure and pressure derivative curves are an upward straight line with a slope of 0.5. This regime is typical of fluid flow for acid fracturing wells. Before the interface is felt, the wellbore pressure solution in this study can be approximate to(80)pwDtD≈2βαtDπ.Corresponding to the square root of the dimensionless time, this equation implies that a log-log plot of function results a straight line with the slope of one-half, which is the typical characteristic of linear flow. This equation also indicates that the wellbore pressure in this regime depends on the values of interregional diffusivity ratio and interregional conductivity ratio (seen in Figures 6 and 7).Figure 6
Effect of interregional diffusivity ratio on transient flow behavior.Figure 7
Effect of interregional conductivity ratio on transient flow behavior.(  2) Transition Flow Regime (a). The first transition flow takes place when the interface takes effect. The derivative curve in this regime shows the first V-shaped segment, which represents fluid flow from natural fractures into artificial fractures or wormholes. It can be interpreted by the fact that as the pressure of the inner region depletes gradually, natural fracture system responds first owing to the direct connection to the inner region. The shape of this segment is also controlled by interregional diffusivity ratio and interregional conductivity ratio, and this V-shaped segment may not be observed at particular parameters (seen in Figures 6 and 7).(  3) Interporosity Flow Regime (a). The first interporosity flow regime corresponds to the second V-shaped segment on the derivative curve, which represents the crossflow between vug system and natural fracture system. This interporosity flow takes place first because vug permeability is larger than matrix permeability. This segment is predominantly determined by storativity ratios of natural fracture system and matrix system, and interporosity coefficient between vugs and natural fractures, and also may be absent from the curves at particular parameters (seen in Figures 8–10, and 13).Figure 8
Effect of fracture storativity ratio on transient flow behavior.Figure 9
Effect of matrix storativity ratio on transient flow behavior.Figure 10
Effect of vug storativity ratio on transient flow behavior.(  4) Transition Flow Regime (b). With the depletion of vug system, the matrix system begins to supply fluid to natural fracture system gradually and a short transition flow occurs. The duration of this regime depends on the difference of the ability supplying fluid between matrix system and vug system; thus, this regime may be so short that it can hardly be identified from the derivative curve.(  5) Interporosity Flow Regime (b). In this interporosity flow regime, the production is mainly dominated by the fluid flow from matrix system to fracture system. This regime is characterized by the third V-shaped segment on the derivative curve, which is governed by two physical processes: one is interporosity flow from matrix system to fracture system and the other is interporosity flow from vug system to matrix system. The former has greater influence on the shape of this regime than the latter (seen in Figures 11 and 12).Figure 11
Effect of matrix-fracture interporosity coefficient on transient flow behavior.Figure 12
Effect of matrix-vug interporosity coefficient on transient flow behavior.Figure 13
Effect of vug-fracture interporosity coefficient on transient flow behavior.(  6) Transition Flow Regime (c). With the depletion of matrix system, fluids far away from the wellbore begin to flow. Like the second transition flow, another short transition flow regime may develop before the dynamic balance reaches the whole reservoir. It also may not be observed on the pressure derivative curve.(  7) Radial Flow Regime. When the matrix system depletes fast and the reservoir radius is large enough, fluid flow in the whole reservoir is likely to reach a dynamic balance; thus, a radial flow occurs after the third transition flow. In this regime, the slope of pressure derivative curve is zero, and the pressure derivative converges to “0.5 lines.” As shown in Figure 14, reservoir radius has a great influence on the duration of radial flow.Figure 14
Effect of dimensionless reservoir radius on transient flow behavior.(  8) Pseudo-Steady Flow Regime. For closed boundary, the long-time behavior corresponds to a pseudo-steady flow. Both the pressure and derivative curves show a unit slope straight line, which is the typical characteristic of closed reservoirs. By means of asymptotic analysis, the approximate solution in this regime can be yielded:(81)pwDtD≈2ω2ftDreD2-1+reD22reD2-1ln⁡ω2ftD+reD2ln⁡4-0.5772+12reD2-1.This solution shows that wellbore pressure in pseudo-steady regime only depends on natural fracture storativity ratio and reservoir radius (seen in Figures6–14).


### 6.1.2. Case 2: Hydrocarbon Reserves Dominated by Vug System

As is shown in Figure5, flow in this reservoir can be divided into six flow regimes: linear flow, transition flow (a), interporosity flow, transition flow (b), radial flow, and pseudo-steady flow. The major difference of flow characteristics between two kinds of carbonate reservoirs lies in the interporosity flow period. The characteristics of the other flow regimes are the same with the corresponding regimes of the previous case, including the approximate solutions of linear flow and pseudo-steady flow.Because the ability of vug system supplying fluid to natural fracture system is much stronger than that of matrix system supplying fluid to natural fracture system, the last two V-shaped segments are connected to each other and form one V-shaped segment. Consequently, only two V-shaped segments can be observed on the pressure derivative curve and the second V-shaped segment merges with the radial flow.


## 6.1.1. Case 1: Hydrocarbon Reserves Dominated by Matrix System

As is described in Figure4, the flow process in this reservoir can be divided into eight flow regimes. Initially, there is a linear flow period characterized by a half-slope straight line on the pressure derivative curve. After a transition flow regime, the first interporosity flow develops and is followed with the second transition flow. Then, the second interporosity flow occurs and immediately is followed by the third transition flow. Eventually, a short radial flow regime may be observed before the system reaches a pseudo-steady flow. Interestingly, three characteristic V-shaped segments can be observed on the pressure derivative curve.(  1) Linear Flow Regime. In this regime, the production is dominated by the fluid stored in artificial fractures and wormholes. It occurs at very small values of dimensionless time. The pressure and pressure derivative curves are an upward straight line with a slope of 0.5. This regime is typical of fluid flow for acid fracturing wells. Before the interface is felt, the wellbore pressure solution in this study can be approximate to(80)pwDtD≈2βαtDπ.Corresponding to the square root of the dimensionless time, this equation implies that a log-log plot of function results a straight line with the slope of one-half, which is the typical characteristic of linear flow. This equation also indicates that the wellbore pressure in this regime depends on the values of interregional diffusivity ratio and interregional conductivity ratio (seen in Figures 6 and 7).Figure 6
Effect of interregional diffusivity ratio on transient flow behavior.Figure 7
Effect of interregional conductivity ratio on transient flow behavior.(  2) Transition Flow Regime (a). The first transition flow takes place when the interface takes effect. The derivative curve in this regime shows the first V-shaped segment, which represents fluid flow from natural fractures into artificial fractures or wormholes. It can be interpreted by the fact that as the pressure of the inner region depletes gradually, natural fracture system responds first owing to the direct connection to the inner region. The shape of this segment is also controlled by interregional diffusivity ratio and interregional conductivity ratio, and this V-shaped segment may not be observed at particular parameters (seen in Figures 6 and 7).(  3) Interporosity Flow Regime (a). The first interporosity flow regime corresponds to the second V-shaped segment on the derivative curve, which represents the crossflow between vug system and natural fracture system. This interporosity flow takes place first because vug permeability is larger than matrix permeability. This segment is predominantly determined by storativity ratios of natural fracture system and matrix system, and interporosity coefficient between vugs and natural fractures, and also may be absent from the curves at particular parameters (seen in Figures 8–10, and 13).Figure 8
Effect of fracture storativity ratio on transient flow behavior.Figure 9
Effect of matrix storativity ratio on transient flow behavior.Figure 10
Effect of vug storativity ratio on transient flow behavior.(  4) Transition Flow Regime (b). With the depletion of vug system, the matrix system begins to supply fluid to natural fracture system gradually and a short transition flow occurs. The duration of this regime depends on the difference of the ability supplying fluid between matrix system and vug system; thus, this regime may be so short that it can hardly be identified from the derivative curve.(  5) Interporosity Flow Regime (b). In this interporosity flow regime, the production is mainly dominated by the fluid flow from matrix system to fracture system. This regime is characterized by the third V-shaped segment on the derivative curve, which is governed by two physical processes: one is interporosity flow from matrix system to fracture system and the other is interporosity flow from vug system to matrix system. The former has greater influence on the shape of this regime than the latter (seen in Figures 11 and 12).Figure 11
Effect of matrix-fracture interporosity coefficient on transient flow behavior.Figure 12
Effect of matrix-vug interporosity coefficient on transient flow behavior.Figure 13
Effect of vug-fracture interporosity coefficient on transient flow behavior.(  6) Transition Flow Regime (c). With the depletion of matrix system, fluids far away from the wellbore begin to flow. Like the second transition flow, another short transition flow regime may develop before the dynamic balance reaches the whole reservoir. It also may not be observed on the pressure derivative curve.(  7) Radial Flow Regime. When the matrix system depletes fast and the reservoir radius is large enough, fluid flow in the whole reservoir is likely to reach a dynamic balance; thus, a radial flow occurs after the third transition flow. In this regime, the slope of pressure derivative curve is zero, and the pressure derivative converges to “0.5 lines.” As shown in Figure 14, reservoir radius has a great influence on the duration of radial flow.Figure 14
Effect of dimensionless reservoir radius on transient flow behavior.(  8) Pseudo-Steady Flow Regime. For closed boundary, the long-time behavior corresponds to a pseudo-steady flow. Both the pressure and derivative curves show a unit slope straight line, which is the typical characteristic of closed reservoirs. By means of asymptotic analysis, the approximate solution in this regime can be yielded:(81)pwDtD≈2ω2ftDreD2-1+reD22reD2-1ln⁡ω2ftD+reD2ln⁡4-0.5772+12reD2-1.This solution shows that wellbore pressure in pseudo-steady regime only depends on natural fracture storativity ratio and reservoir radius (seen in Figures6–14).


## 6.1.2. Case 2: Hydrocarbon Reserves Dominated by Vug System

As is shown in Figure5, flow in this reservoir can be divided into six flow regimes: linear flow, transition flow (a), interporosity flow, transition flow (b), radial flow, and pseudo-steady flow. The major difference of flow characteristics between two kinds of carbonate reservoirs lies in the interporosity flow period. The characteristics of the other flow regimes are the same with the corresponding regimes of the previous case, including the approximate solutions of linear flow and pseudo-steady flow.Because the ability of vug system supplying fluid to natural fracture system is much stronger than that of matrix system supplying fluid to natural fracture system, the last two V-shaped segments are connected to each other and form one V-shaped segment. Consequently, only two V-shaped segments can be observed on the pressure derivative curve and the second V-shaped segment merges with the radial flow.


## 6.2. Parameter Sensitivity Analysis

The shape rendered from real data may be distorted by noises in well testing, which makes it necessary to establish the stylized shapes under different parameter conditions [36]. Figures 6–14 contain the type curves of pressure-transient response for an acid fracturing well in FVCRs. Later we will show how the shape of the type curves may be affected by changing the parameters.


### 6.2.1. Effect of Interregional Diffusivity Ratio

According to the dimensionless definition, interregional diffusivity ratioα represents the ratio of fluid transmissivity ability between natural fractures in the outer region and artificial fractures or wormholes in the inner region. A larger α means a weaker diffusion ability of wormholes or artificial fractures. Figure 6 shows the effect of this parameter on pressure-transient response (α=0.001, 0.005, 0.05, and 0.5). We can find that interregional diffusivity ratio mainly affects early-flow period. Before the first interporosity flow regime, a larger α leads to lower pressure depletion at the same production rate. Meanwhile, as this parameter increases, due to a weaker diffusion ability in the inner region, duration of the linear flow is prolonged, while duration of the first interporosity flow is shortened. When interregional diffusivity ratio gradually increases to 1, the first V-shaped segment becomes shallower and even disappears. All the curves for different interregional diffusivity ratios begin to normalize at the end of the first interporosity flow regime.


### 6.2.2. Effect of Interregional Conductivity Ratio

According to the dimensionless definition, interregional conductivity ratioβ represents a measure of the conductivity between wormholes or artificial fractures and natural fractures. In practice, a larger β means a lower conductivity of the inner region. Transient responses on the log-log plots for varied interregional conductivity ratio are shown in Figure 7 (β=0.2, 0.4, 0.6, and 0.8). It can be concluded that this parameter mainly affects early- and middle-flow periods. A larger interregional conductivity ratio leads to more pressure depletion at the same production. Meanwhile, with the increase of this parameter, the hump corresponding to the first transition flow becomes more significant. Besides, the pressure curves for varied interregional conductivity ratio normalize in the late-flow period, while the pressure derivative curves begin to normalize from the first interporosity flow regime.


### 6.2.3. Effect of Storativity Ratios

There are three storativity ratios in a FVCR. Each ratio represents the relative capacity of fluid stored in a corresponding system. Generally fluid stored in the natural fracture system is much less than that stored in the vug or matrix system.Figure8 presents the results of sensitivity analysis of fracture storativity ratio (ω2f = 0.001, 0.005, 0.025, and 0.125). It indicates that a larger ω2f leads to more pressure depletion in the middle- and late-flow periods. Fracture storativity ratio begins to affect the transient response at the start of the first interporosity flow, which extends to the pseudo-state flow period; thus, the curves do not normalize, but they are parallel to each other in the late-flow period. Due to a higher conductivity and diffusion capacity for natural fracture system provided by a larger fracture storativity ratio, when this ratio increases from 0.001 to 0.125, the last two V-shaped segments become shallower and even the second V-shaped segment may disappear. Moreover, the pseudo-steady flow begins earlier with the increase of this parameter.When investigating the effect of another two storativity ratios, we keepω2f constant; thus, matrix storativity ratio has a reciprocal relationship with vug storativity ratio. Therefore, it is not necessary to discuss all the values in their own range.Figure9 shows the response of varying matrix storativity ratio, ω2m = 0.6, 0.7, 0.8, and 0.9, which indicates that matrix system provides the majority of storage space. As is described in Figure 9, matrix storativity ratio mainly affects the middle-flow period, and the pressure and derivative curves for different matrix storativity ratios overlap completely in the early- and late-flow periods. But the shape of the last two V-shaped segments is affected by this storativity ratio significantly. With the increase of matrix storativity ratio, the second V-shaped segment becomes shallower, while the third V-shaped segment becomes deeper. In fact, matrix storativity ratio only affects the third V-shaped segment and vug storativity ratio only affects the second V-shaped segment. However, due to the reciprocal relationship between two storativity ratios, matrix storativity ratio affects the second V-shaped segment indirectly.Figure10 presents the results of sensitivity analysis of vug storativity ratio, ω2v = 0.6, 0.7, 0.8, and 0.9, which indicates that the reserves are dominated by vug system. As shown in Figure 10, vug storativity ratio only has a slight effect on the middle-flow period and the related curves normalize in the late-flow period. With the increase of vug storativity ratio, the concave of the second V-shaped segment deepens slightly.


### 6.2.4. Effect of Interporosity Coefficients

There are also three interporosity coefficients in the triple-porosity reservoir. A sensitivity analysis is also performed by varying these three coefficients.Figure11 demonstrates the effect of matrix-fracture interporosity coefficient on the transient behavior (λmf = 5 × 10−5, 5 × 10−4, 5 × 10−3, and 10−2). It can be seen that this parameter mainly affects the middle-flow period and type curves for different interporosity coefficients normalize in the late-flow period. As this coefficient increases, the second V-shaped segment deepens, the third V-shaped segment moves to the left side, and the duration of radial flow regime is prolonged. This is mainly because a larger matrix-fracture interporosity coefficient means a stronger ability of matrix system supplying fluid for natural fracture system, which makes the arrival time of matrix-fracture interporosity flow regime earlier.Figure12 presents the results of sensitivity analysis of matrix-vug interporosity coefficient (λmv = 0.001, 0.005, 0.025, and 0.05). This coefficient also mainly affects middle-flow period, especially the matrix-fracture interporosity flow regime, corresponding to the third V-shaped segment. It can be interpreted by that this coefficient only affects the ability of vug system to supply fluid for matrix system. With the increase of this coefficient, the starting time of the second interporosity flow moves back and the corresponding V-shaped segment deepens greatly. Similarly, the pressure and derivative curves will normalize in the late-flow period.Figure13 shows the response of varying vug-fracture interporosity coefficient (λvf = 0.01, 0.05, 0.25, and 0.5). This parameter mainly affects the middle-flow period, and the normalization can also be observed in the late-flow period. As this coefficient increases, the starting time of the first interporosity flow moves back and the third V-shaped segment becomes shallower. It can be interpreted by that the ability of vug system to supply fluid for fracture system becomes much stronger than that of matrix system supplying fluid for fracture system with the increase of this coefficient.


### 6.2.5. Effect of Dimensionless Reservoir Radius

Figure14 presents the results of sensitivity analysis of dimensionless reservoir radius (reD = 30, 50, 70, and 90). There is no difference in the shape of type curves in the early- and middle-flow periods under different reservoir radius; however, the pressure and pressure derivative curves show a unit slope straight line and go up rapidly in the late-flow period. In addition, with the increase of the dimensionless reservoir radius, the duration of radial flow regime is prolonged significantly.


## 6.2.1. Effect of Interregional Diffusivity Ratio

According to the dimensionless definition, interregional diffusivity ratioα represents the ratio of fluid transmissivity ability between natural fractures in the outer region and artificial fractures or wormholes in the inner region. A larger α means a weaker diffusion ability of wormholes or artificial fractures. Figure 6 shows the effect of this parameter on pressure-transient response (α=0.001, 0.005, 0.05, and 0.5). We can find that interregional diffusivity ratio mainly affects early-flow period. Before the first interporosity flow regime, a larger α leads to lower pressure depletion at the same production rate. Meanwhile, as this parameter increases, due to a weaker diffusion ability in the inner region, duration of the linear flow is prolonged, while duration of the first interporosity flow is shortened. When interregional diffusivity ratio gradually increases to 1, the first V-shaped segment becomes shallower and even disappears. All the curves for different interregional diffusivity ratios begin to normalize at the end of the first interporosity flow regime.


## 6.2.2. Effect of Interregional Conductivity Ratio

According to the dimensionless definition, interregional conductivity ratioβ represents a measure of the conductivity between wormholes or artificial fractures and natural fractures. In practice, a larger β means a lower conductivity of the inner region. Transient responses on the log-log plots for varied interregional conductivity ratio are shown in Figure 7 (β=0.2, 0.4, 0.6, and 0.8). It can be concluded that this parameter mainly affects early- and middle-flow periods. A larger interregional conductivity ratio leads to more pressure depletion at the same production. Meanwhile, with the increase of this parameter, the hump corresponding to the first transition flow becomes more significant. Besides, the pressure curves for varied interregional conductivity ratio normalize in the late-flow period, while the pressure derivative curves begin to normalize from the first interporosity flow regime.


## 6.2.3. Effect of Storativity Ratios

There are three storativity ratios in a FVCR. Each ratio represents the relative capacity of fluid stored in a corresponding system. Generally fluid stored in the natural fracture system is much less than that stored in the vug or matrix system.Figure8 presents the results of sensitivity analysis of fracture storativity ratio (ω2f = 0.001, 0.005, 0.025, and 0.125). It indicates that a larger ω2f leads to more pressure depletion in the middle- and late-flow periods. Fracture storativity ratio begins to affect the transient response at the start of the first interporosity flow, which extends to the pseudo-state flow period; thus, the curves do not normalize, but they are parallel to each other in the late-flow period. Due to a higher conductivity and diffusion capacity for natural fracture system provided by a larger fracture storativity ratio, when this ratio increases from 0.001 to 0.125, the last two V-shaped segments become shallower and even the second V-shaped segment may disappear. Moreover, the pseudo-steady flow begins earlier with the increase of this parameter.When investigating the effect of another two storativity ratios, we keepω2f constant; thus, matrix storativity ratio has a reciprocal relationship with vug storativity ratio. Therefore, it is not necessary to discuss all the values in their own range.Figure9 shows the response of varying matrix storativity ratio, ω2m = 0.6, 0.7, 0.8, and 0.9, which indicates that matrix system provides the majority of storage space. As is described in Figure 9, matrix storativity ratio mainly affects the middle-flow period, and the pressure and derivative curves for different matrix storativity ratios overlap completely in the early- and late-flow periods. But the shape of the last two V-shaped segments is affected by this storativity ratio significantly. With the increase of matrix storativity ratio, the second V-shaped segment becomes shallower, while the third V-shaped segment becomes deeper. In fact, matrix storativity ratio only affects the third V-shaped segment and vug storativity ratio only affects the second V-shaped segment. However, due to the reciprocal relationship between two storativity ratios, matrix storativity ratio affects the second V-shaped segment indirectly.Figure10 presents the results of sensitivity analysis of vug storativity ratio, ω2v = 0.6, 0.7, 0.8, and 0.9, which indicates that the reserves are dominated by vug system. As shown in Figure 10, vug storativity ratio only has a slight effect on the middle-flow period and the related curves normalize in the late-flow period. With the increase of vug storativity ratio, the concave of the second V-shaped segment deepens slightly.


## 6.2.4. Effect of Interporosity Coefficients

There are also three interporosity coefficients in the triple-porosity reservoir. A sensitivity analysis is also performed by varying these three coefficients.Figure11 demonstrates the effect of matrix-fracture interporosity coefficient on the transient behavior (λmf = 5 × 10−5, 5 × 10−4, 5 × 10−3, and 10−2). It can be seen that this parameter mainly affects the middle-flow period and type curves for different interporosity coefficients normalize in the late-flow period. As this coefficient increases, the second V-shaped segment deepens, the third V-shaped segment moves to the left side, and the duration of radial flow regime is prolonged. This is mainly because a larger matrix-fracture interporosity coefficient means a stronger ability of matrix system supplying fluid for natural fracture system, which makes the arrival time of matrix-fracture interporosity flow regime earlier.Figure12 presents the results of sensitivity analysis of matrix-vug interporosity coefficient (λmv = 0.001, 0.005, 0.025, and 0.05). This coefficient also mainly affects middle-flow period, especially the matrix-fracture interporosity flow regime, corresponding to the third V-shaped segment. It can be interpreted by that this coefficient only affects the ability of vug system to supply fluid for matrix system. With the increase of this coefficient, the starting time of the second interporosity flow moves back and the corresponding V-shaped segment deepens greatly. Similarly, the pressure and derivative curves will normalize in the late-flow period.Figure13 shows the response of varying vug-fracture interporosity coefficient (λvf = 0.01, 0.05, 0.25, and 0.5). This parameter mainly affects the middle-flow period, and the normalization can also be observed in the late-flow period. As this coefficient increases, the starting time of the first interporosity flow moves back and the third V-shaped segment becomes shallower. It can be interpreted by that the ability of vug system to supply fluid for fracture system becomes much stronger than that of matrix system supplying fluid for fracture system with the increase of this coefficient.


## 6.2.5. Effect of Dimensionless Reservoir Radius

Figure14 presents the results of sensitivity analysis of dimensionless reservoir radius (reD = 30, 50, 70, and 90). There is no difference in the shape of type curves in the early- and middle-flow periods under different reservoir radius; however, the pressure and pressure derivative curves show a unit slope straight line and go up rapidly in the late-flow period. In addition, with the increase of the dimensionless reservoir radius, the duration of radial flow regime is prolonged significantly.


## 6.3. Comparison with Other Similar Models

Gringarten et al. [38] investigated unsteady-state pressure distributions created by a well with a single infinite-conductivity vertical fracture in a homogeneous reservoir. Further, we can easily obtain the transient-pressure solution for the same well type in a FVCR. There are a finite number of artificial fractures and wormholes in our model’s inner region; thus, it is necessary to compare these two models’ difference in flow behavior. Meanwhile, it is also significant to make a comparison of our model with the model not having acid fracturing to know the effect of acidification. In order to have a convenient comparison, we firstly make a comparison to the parameters used in three models and then compare their type curves by fixing a set of the same parameter values. Table 3 shows the differences between three models’ parameters. It can be seen that the same parameters are ω2f, ω2m, ω2v, λmf, λmv, λvf, rwD, and reD, but the parameters α and β only belong to the model with acid fracturing. To make an effective comparison, we set the same values for the same parameters, which are shown in Table 3.Table 3
The comparison of different parameters between different models.
Model Model parameters Acid fracturing α β ω 2 f ω 2 m ω 2 v λ m f λ m v λ v f r w D r e D 0.001 0.05 0.001 0.96 0.039 0.001 0.01 0.1 0.001 40 Single-vertical fracture / / ω 2 f ω 2 m ω 2 v λ m f λ m v λ v f r w D r e D / / 0.001 0.96 0.039 0.001 0.01 0.1 0.001 40 No acid fracturing / / ω 2 f ω 2 m ω 2 v λ m f λ m v λ v f r w D r e D / / 0.001 0.96 0.039 0.001 0.01 0.1 0.001 40Figure15 shows the comparison result of type curves between three models, and we can observe the differences clearly: (i) in the early-flow period, the pressure and derivative curves for the model with acid fracturing and the model with single fracture coincide perfectly, respectively. It demonstrates that these two models’ flow behavior is the same in this period. However, because drainage area of the well without acid fracturing is much smaller than that of the well with acid fracturing or with single fracture, the pressure depletion of the former case is faster than that of the latter two cases at the same group of production parameters; (ii) in the middle-flow period, the pressure depletion of the well without single fracture becomes faster than that of the well with acid fracturing, which also can be interpreted by the difference of drainage area; (iii) in the models with single fracture and without acid fracturing, only two V-shaped segments can be observed on the derivative curves, and the first V-shaped segment in the model with acid fracturing is replaced by a unit slope straight line or a horizontal line with the value of 0.5, which are the typical characteristics of linear flow and radial flow, respectively.Figure 15
The comparison of type curves between three different models.


## 6.4. Parameter Influence on Pseudo-Skin Factor

During the pseudo-steady flow period, for the same position, there is a difference of the dimensionless pressure between the models with acid fracturing and without acid fracturing. This difference remains constant at large value of dimensionless time and can be handled as a pseudo-skin factor [39, 40] that can be defined as(82)SrD,tD;α,β,λ,ω,reD=p1DrD,tD;α,β,λvf,λmv,λmf,ω2f,ω2v,ω2m,reD-pNAFDrD,tD;λvf,λmv,λmf,ω2f,ω2v,ω2m,reD, where p1D and pNAFD represent the dimensionless pressure near the wellbore for the model with acid fracturing and without the stimulation, respectively.Applying the Laplace transform to (82), we have(83)S~rD,s;α,β,λ,ω,reD=p~1DrD,s;α,β,λvf,λmv,λmf,ω2f,ω2v,ω2m,reD-p~NAFDrD,s;λvf,λmv,λmf,ω2f,ω2v,ω2m,reD,where(84)p~1D=βssαx2sinh⁡sα1-rD-x1cosh⁡sα1-rDx2cosh⁡sα-x1sinh⁡sα,p~NAFD=1sK0rDsfs+K1reDsfsI1reDsfsI0rDsfs.Note that x1 and x2 have been given in (49).Because of higher connectivity between the wellbore and the reservoir caused by acid fracturing, the pseudo-skin factor defined here is always negative and indicates an increase in well productivity. The larger the absolute value of the pseudo-skin factor, the better the stimulated effect of acid fracturing.As is shown in Figure16, under different interregional diffusivity ratios, the pseudo-skin factors at the same position near the wellbore are the same, which indicates that this parameter has no influence on the ability of enhancing the productivity in the late-flow period. In the pseudo-steady flow regime, the diffusion capacity of the inner region well satisfies the demand of fluid flow feeding the wellbore and the flow in the model reaches a dynamic balance; thus, the change of interregional diffusivity ratio does not affect the productivity. It also can be concluded that the closer the position to the wellbore, the better the stimulation effect.Figure 16
Pseudo-skin factor versus dimensionless formation radius for an acid fracturing well at different interregional diffusivity ratios.Figure17 presents different pseudo-skin factors for the vertical well of varied interregional conductivity ratio. With the increase of interregional conductivity ratio, the absolute value of the pseudo-skin factor at the same position decreases, which implies that the ability of enhancing the well production decreases. As we know, a larger interregional conductivity ratio corresponds to a weaker conductivity of the inner region and also means the increase of the flow resistance, which explains the result above. Similarly, in this group, the closer the position to the wellbore, the better the stimulation effect.Figure 17
Pseudo-skin factor versus dimensionless reservoir radius for an acid fracturing well at different interregional conductivity ratios.


## 7. Conclusions

This work presented a composite reservoir model to analyze the flow behavior of an acid fracturing well in a fractured-vuggy carbonate reservoir. Both the analytical and numerical solutions are obtained in this study. Based on this work, several important conclusions are obtained:( 1 ) All the analytical solutions show well agreement with the corresponding numerical solutions within the most production life for different parameters. These verifications indicate that the analytical model in this work is reliable. The numerical solution has some advantages over the analytical solution in computing speed, but the analytical solution is convenient to handle the condition with wellbore storage or skin effect.( 2 ) Due to that reserves in the fractured-vuggy carbonate reservoir are dominated by either matrix system or vug system; for an acid fracturing well, flow in this composite model can be divided into six or eight flow regimes comprehensively. Linear flow regime is typical of fluid flow for acid fracturing wells in such reservoirs. For the two scenarios above, three and two characteristic V-shaped segments can be observed on their own pressure derivative curves, respectively.( 3 ) Effects of interregional diffusivity ratio and interregional conductivity ratio on transient responses are strong in the early-flow period. The existence and duration of linear flow and the first transition flow are mainly determined by these two parameters. As interregional diffusivity ratio increases, duration of linear flow is prolonged while that of the first transition flow is shortened. Besides, a lower interregional diffusivity ratio or a larger interregional conductivity ratio leads to more pressure depletion in the early-flow period at the same production.( 4 ) The shape and position of type curves are also dominated by storativity ratios and interporosity coefficients. In the late-flow period, for different interregional diffusivity ratio, interregional conductivity ratio, and interporosity coefficient, the pressure and pressure derivative curves will normalize. However, type curves for varied fracture storativity ratio and dimensionless reservoir radius are parallel to each other.( 5 ) Compared with the models with single fracture or without acid fracturing, the pressure depletion of the well with acid fracturing is slower than that of another two models at the same group of production parameters in the middle- and late-flow period. Due to a higher connectivity caused by acid fracturing, the pseudo-skin factor defined in this work is always negative and indicates an increase in well productivity. Interregional diffusivity ratio has no influence on the ability of enhancing the productivity in the late-flow period. However, the ability of enhancing the well production decreases with the increase of interregional conductivity ratio.



---

*Source: 6431910-2018-03-27.xml*