# A New Model for Predicting Flow in Fractured-Vuggy Carbonate Reservoirs under Pseudosteady State Condition

Authors: Nai Cao; Gang Lei; Yuan Su; Hong Li; Jiali Zhang; Ci Zhang
Journal: Geofluids
 (2021)

---

## Abstract

For constant rate production of a single well in a closed boundary fractured-vuggy carbonate reservoir, when the pressure wave propagates to the reservoir boundaries, boundary-dominated flow occurs, and the transient flow period ends, the reservoir will be in a state of pseudoequilibrium (i.e., pseudosteady state flow occurs, and the pressure at any point in the reservoir declines at the same constant rate over time). The characterization of fluid flow in the fractured-vuggy carbonate reservoir under pseudosteady state condition is of significance in describing the productivity index (PI) of the well. However, due to complex mechanisms (e.g., stress dependency of reservoir properties and crossflow between different systems of the reservoirs) during flow in fractured-vuggy carbonate reservoirs, researches on productivity prediction in fractured-vuggy carbonate reservoirs under pseudosteady state are very limited. The present work is aimed at developing a new analytical model for predicting flow in fractured-vuggy carbonate reservoirs under pseudosteady state condition. In the derived model, the crossflow between different systems (i.e., matrix, fracture, and vug) of the reservoirs was taken into account. In addition, based on Hooke’s law, a quantitative model was proposed to study stress-dependent permeability of the fracture system, which connects the reserve spaces. Moreover, the roughness morphology characteristics of fracture surface were taken into account with fractal theory. Finally, with our derived model, influences caused by various related factors on productivity were analyzed. The results show that well productivity during pseudosteady flow will be significantly affected by the morphology of fracture surface (e.g., fracture microstructure parameters) and effective stress. Specifically, due to effective stress, the fracture system in the fractured-vuggy reservoirs will be deformed, and the corresponding properties (e.g., permeability, porosity, and conductivity) will change, leading to the change of well productivity. In addition, there exists a negative relationship between the elastic storage capacity ratio of vugs and well productivity during pseudosteady flow. Moreover, a larger value of matrix-fracture interporosity flow coefficient or vug-fracture interporosity flow coefficient corresponds to a larger value of well productivity during the pseudosteady period. The new derived model is beneficial to improve the productivity prediction accuracy and reduce uncertainty. What is more, the findings of this study can help for providing theoretical reference for the design of efficient development of fractured-vuggy carbonate reservoirs.


---

## Body

## 1. Introduction

Huge reserves are explored in carbonate rocks worldwide, and the known economic reserves of it increase year by year [1, 2]. For example, more than half of the world’s largest crude oil and natural gas reserves are found in carbonate rocks. In addition, carbonate rocks produce a significant portion of the world’s oil and gas. Thus, how to efficiently exploit these resources attracts attention all over of the world [1–4]. Taking China for example, according to the statistics of the Ministry of Natural Resources of China, in 2018, 7176×104t of proved oil and gas reserves of Ordovician carbonate rocks (fractured-vuggy carbonate reservoirs) in Shunbei was reported [2]. As stated in the literature [3–8], in the near future, the resources of carbonate reservoirs will be one of the main battlefields for “increasing reserves and productivity” in China.In general, influenced by a variety of geological processes, the carbonate reservoirs are heterogeneous, and the reservoir spaces are typically dual or triple porosity systems, which are different in size and complex in distribution. For fractured-vuggy carbonate reservoirs, the reservoirs can be classified as matrix system, fracture system, and vug system, in which the vug system is the main reservoir space, and fractures with larger permeability work as the main flow channels [5–7]. Physically speaking, during the developing process, formation permeability of each system in the reservoirs will decrease with the increasing effective stress [8]. For example, affected by effective stress, matrix compression (or vug deformation) will result in permeability decease of matrix system (or vug systems), and fracture closure will lead to permeability decrease of fracture system. Yan et al. [9] conducted stress sensitivity experiments on 12 core samples (e.g., natural fractured samples, artificial fractured samples, single cavity samples, and double cavity samples) from Tahe oilfield and found that the stress sensitive of vug system was weak. In addition, they suggested that the vug system would weak the stress sensitive degree of fractured-vuggy reservoirs [9]. That is, compared with permeability decrease due to matrix compression or deformation of vug systems, fracture closure plays a more crucial role in causing the decrease in formation permeability [8–12]. Thus, the matrix system stress sensitivity (or vug system stress sensitivity) is much less evident in the fractured-vuggy formation. Moreover, when the fractures are highly penetrated, formation stress sensitivity will be dominated by fracture system stress sensitivity [8–13]. In this paper, to simplify the model, the fracture system is assumed to be highly penetrated; thus, the formation stress sensitivity is mostly determined by the stress sensitivity of the fracture system [8–13]. It is well known that, with the decrease of formation permeability due to effective stress, well productivity will decrease, and the effective development becomes more difficult [8–13]. Therefore, it is of theoretical and scientific significance to quantify the stress sensitivity of fracture system permeability and its impact on productivity of fractured-vuggy carbonate reservoirs [14–18].In essence, the surface of the fracture system is rough and has complex morphological characteristics, making the effective stress on the fracture surface unevenly distributed, leading to the extremely complex stress sensitivity of the fracture system [19–23]. However, to the best of our knowledge, many former studies seldom coupled the microstructure characteristics of the fracture surface with effective stress, and these models contained empirical parameters, whose physical meaning was unclear [18–21]. In 1978, Gangi proposed a theoretical model to quantify the stress dependent permeability of fractures [19]. In this derived model, the roughness of fractures was represented by small cylinders unevenly distributed on smooth fracture surfaces. After that, a number of scholars continued to deepen the understanding of stress sensitivity of fractured media [10–13, 15, 20–22]. For example, by combining the Carman-Kozeny model and a derived stress dependent porosity model, Mckee et al. studied the stress dependent permeability of fractures considering the fracture geometry [15]. Based on fractal theory and Gangi’s model, Lei et al. proposed theoretical models to study stress dependent permeability and relative permeability of rough fractures [11, 13]. Moreover, some scholars used small cylinders with different elastic modulus to characterize fracture fillers and studied the closure law of filled fractures under effective stress with the Hertz elastic contact model [20, 21]. Recently, Cao et al. took the fractures as distinct objects and derived a theoretical model to study stress dependent permeability of fractures [22]. In this derived model, they also studied the influence of the extent of fracture penetration on stress-dependent permeability of fractures. Based on the discussions above, it is clear that, though many related work has been done, the relative theoretical research is yet to be further improved for accurate description of fracture system stress sensitivity. Thus, to precisely analyze the stress sensitivity of the rough fracture system, the influence of fracture microstructure on stress sensitivity of the fracture system, which further influence the well productivity, remains to be clarified and analyzed.It is well known that, to better understand the flow in porous media, it is of great significance to conduct productivity predictions. To the best of our knowledge, for productivity prediction of fractured-vuggy reservoirs, applicable mathematical models have been developed from original homogeneous model to double and triple medium models. However, until now, related literature devoted to develop inflow performance relationship curve (IPR curve, i.e., production versus drawdown pressure) of vertical well in fractured-vuggy reservoirs is scarce. In addition, though stress sensitivity has been coupled with dual medium model, the influence of fracture stress sensitivity on well productivity of fractured-vuggy reservoirs has rarely been reported. However, physically speaking, the stress sensitivity of fracture systems directly restricts the communication ability of fractured-vuggy reservoirs under overburden pressure and the fluid supply efficiency of reservoirs to wellbore [24, 25]. Therefore, the influence mechanism of stress sensitivity of rough fracture on productivity of fractured-vuggy carbonate reservoirs is await to be further clarified.For a constant production well in a closed fractured-vuggy reservoir, fluid flow will occur in early-time, middle-time, and late-time regions, respectively [26]. In general, in the early-time region, fluid flow will be affected by wellbore. When the fluid flow progresses to the middle-time region, the flow will be driven by outward diffusion of the pressure drawdown (i.e., the propagation of pressure wave). As production progresses, when the pressure wave propagates to the closed boundary, fluid flow occurs in the late-time region (i.e., pseudo-steady state flow region). Physically speaking, during pseudosteady state flow region, fluid flow is driven by the volumetric compression of the reservoir, and the rate of change of pressure will be constant over time [27–29]. In particular, for fluid flow in multiple media under pseudosteady state flow region, the fluid transfer between different systems (e.g., fracture-matrix system and fracture-vug system) ais proportional to the pressure changes [30]. To describe the well productivity index (PI), many scholars conducted theoretical works to characterize fluid flow under pseudosteady state condition [31–36]. For example, Lei et al. proposed analytical models to study productivity of wells in tight sandstone reservoirs under pseudosteady state [31, 32]. In the derived models, the stress sensitivity of tight matrix was taken into account. However, in these models, the crossflow between fracture system and matrix system was ignored. With material balance constraint, Morgan proposed a production forecasting model for tight gas well [33]. However, the model never took effective stress into account. Moreover, the crossflow between different systems was ignored. Based on the assumption of linear flow from reservoir to fracture and linear flow in the fracture, Zhang et al. established an analytical productivity model of multifractured shale gas wells [34]. They concluded that the prediction from their model was in good agreement with actual data of gas wells in North American and China. However, their model did not take effective stress into account. Al-Rbeawi studied multiphase flow with different wellbore conditions under pseudosteady state [35]. For triple media, Youssef and Alnumaim proposed an analytical model to study well productivity during pseudosteady state flow region. In their model, the crossflow between different systems was taken into account. However, the derived model also failed to take effective stress into account. Thus, developing a theoretical well productivity model for the pseudosteady state flow is needed to deepen the understanding of fluid flow in fractured-vuggy reservoirs.To sum up, the objectives of this work are (1) to develop an analytical model for predicting well productivity of fractured-vuggy carbonate reservoirs under pseudosteady state, (2) to study the stress dependent permeability of rough fractures and its impact on well productivity, and (3) to explore the influence of various factors on well productivity of fractured-vuggy carbonate reservoirs, thereby providing supports for developing designs. Motivated by this, in this work, a productivity prediction model for vertical wells in fractured-vuggy carbonate reservoirs under pseudosteady state is established, in which morphology surfaces of rough fractures are taken into account. In the derived model, the stress sensitivity of the fracture system is quantitatively presented. Moreover, the triple media (matrix, fracture and vug) model is applied to characterize flow in the fractured-vuggy carbonate reservoirs. With the proposed model, the influence brought by the stress sensitivity of rough fracture system on the well productivity is discussed. Based on the general sketch of this work (Figure1), the outline of this work can be summarized as follows: firstly, the mathematical model and the workflow of solution determination will be described in Section 2. Then, model results and the influence of related parameters on well productivity will be discussed in Section 3. Finally, the conclusions will be summarized and presented in Section 4.Figure 1
General sketch of this work.


## 2. Mathematical Model

### 2.1. Model Assumptions

In this work, the fractured-vuggy carbonate reservoir is simplified as the triple medium shown in Figure2(a). It is assumed that the reservoir is isotropic, and it is consist of evenly distributed matrix, fracture systems, and vug systems. The single-phase fluid in the reservoir is microcompressible, and the flow of which obeys the Darcy’s law. The flow is isothermal during the whole process, and it has reached a pseudosteady state. The fluid supply to the wellbore is stable, and only radial flow in horizontal direction is considered. Figure 2(b) depicts the whole flow process: the fluid flows into the wellbore only through the fracture systems, and there exists crossflow between matrix fracture and vug fracture. That is, both matrix and vug systems will supply fluid to fracture systems.Figure 2
Sketch of the physical model.
(a)
Simplified triple media model(b)
The connection sketch of triple media


### 2.2. Flow Equations and the Solutions

Based on the assumptions stated above, the equations for describing fluid flow in fractured-vuggy reservoir are as follows [36–38]:Fracture system:(1)Kfμ∂2pf∂r2+1r∂pf∂r=αmKmμpf−pm+αvKvμpf−pv+φfCf∂pf∂t.Matrix system:(2)φmCm∂pm∂t=αmKmμpf−pm.Vug system:(3)φvCv∂pv∂t=αvKvμpf−pv,where μ is the fluid viscosity, mPa·s; p is the formation pressure, MPa; K is the permeability, μm2; α is the shape factor, m-2; t is the production time, d; C is the compressibility coefficient, MPa-1; φ is porosity, %; subscripts f, v, and m, respectively, stand for fractures, vugs, and matrix; r is the distance, m.As aforementioned, during pseudosteady state flow region, the pressure drop are stable and consistent throughout the reservoir, and the volumetric compression of the reservoir is the dominated driving force for fluid flow [31–36]. Thus, for fluid flow in fractured-vuggy reservoirs, we have [30, 36]
(4)∂pf∂t=∂pm∂t=∂pv∂t=−qscBVBφfCf+φmCm+φvCv,where B is the fluid bulk coefficient, m3/m3; VB is the fluid volume under formation condition, m3; qsc is the productivity at ground standard condition, m3/d.Substituting Eq. (4) into Eq. (1), we have
(5)Kfμ∂2pf∂r2+1r∂pf∂r=qscBVB.The boundary conditions can be described as(6)pfr=rw=pwf,(7)r∂pf∂rr=re=0,where pwf is the bottom hole pressure, MPa; re is the radius of the reservoir, m; rw is the wellbore radius, m.By solving Eqs. (5)–(7), the pressure at the fracture system is [36]
(8)pfr,t=pwft+qscμB2KfVBre2lnrrw−0.5r2−rw2,VB=πHre2−rw2,where H is the reservoir effective sickness, m.Mathematically speaking, for triple media with multiple systems, the average reservoir pressure can be regarded as the mean of average pressure in each system. Therefore, it is crucial to determine the weighting coefficient of each system. In general, the weighting coefficient is related to pore volume of each system. However, as the fluid flow during pseudosteady state is driven by the volumetric compression of the reservoir, the compressibility coefficient of each system will also greatly affect the weighting coefficient. Since the compressibility coefficient varies with different systems, only pore volume parameter is not enough to be weighting parameters for accurately determining the average reservoir pressure. As stated in the literature, elastic storage capacity ratio of each system is more reasonable to be the weighting parameters [36, 39, 40]. That is, the average reservoir pressure can be obtained as
(9)p¯=p¯fωf+p¯mωm+p¯vωv,where p¯f is the average pressure in the fracture system, p¯m is the average pressure in the matrix system, p¯v is the average pressure in vug system, and the weighting parameter ω is assigned as the elastic storage capacity ratio, which can be defined as
(10)ωi=φiCiφfCf+φvCv+φmCmi=f,v,m.With the assumption that parameterre is much larger than rw, based on Eq. (8), the average pressure in fracture system can be written as
(11)p¯f=∰VBpfdVBπHre2−rw2=pwft+qscμB2πKfHlnrerw−34.Then, based on Eqs. (1) and (3), the average pressure in the matrix system and vug system can be, respectively, obtained as
(12)p¯m=p¯f+qscωmμBrw2λmVBKf,(13)p¯v=p¯f+qscωvμBrw2λvVBKf,where λm is the matrix-fracture interporosity flow coefficient; λv is the vug-fracture interporosity flow coefficient.According to the definition of well production index (PI), from Eqs. (9) to (13), we have [36]
(14)J=qscp¯−pwf=2πKfHμBlnre/rw−3/4⋅D,where
(15)D=11+2rw/re2⋅ωm2/λm+ωv2/λv/lnre/rw−3/4.Eq. (14) demonstrates that, compared with the normal well production index (PI) model for homogeneous single continuum media, the derived well production index (PI) can be treated as two parts. The first part is Darcy’s PI part, which is identical to the normal well production index (PI) model. Physically speaking, this part takes into account how the flow rate is related to pressure drop for the fluid flow in the continuum media. The second part, parameter D, represents how easy the fluid can move from one system to the others, which is dependent on the crossflow parameters.


### 2.3. Stress Sensitivity of the Fracture System

Eq. (14) demonstrates that the permeability of the fracture system plays an important role in the productivity prediction of the fractured-vuggy reservoir in the pseudosteady state stage. However, as stated previously, effective stress will greatly affect the permeability of the fracture system. Thus, it is significant to quantify stress-dependent permeability of the fracture system and couple it with Eq. (14) to study well productivity more accurately. Due to complex morphology of the fracture system, it is difficult to characterize fracture microstructure and establish the quantitative relation between effective stress and fracture permeability. Fortunately, through a series research work, many researchers have found that the fractal geometry theory was applicable in describing the microstructure of the fracture system [11, 13, 41–45]. Specifically, fractal theory has been proved to be applicable and reliable in describing the complex fracture system. For example, many scholars suggested that rough surfaces of fractures had fractal characteristics [41–45]. In this work, for the sake of simplification, fractal theory will be applied to characterize fracture morphology, and then, an analytical model of stress-dependent permeability of the fracture system based on fractal theory and Gangi’s model will be introduced [11, 19, 44].Based on Gangi’s model, roughness of fracture surface can be regarded as “nails” embedded on the smooth fracture walls [11, 19]. Then, for fracture walls distributed with “nails” with various lengths, we assume that the length of “nails” follows fractal theory. The fractal dimension of the fracture system roughness can be expressed as [11]
(16)DR=2−lnφnlnlmin/lmax,where DR is the fractal dimension for fracture surface roughness (0<DR<2), dimensionless. Area ratio of fracture surface φn is the ratio of total roughness area to total surface area of rough fracture, dimensionless; lmin and lmax, respectively, stand for the minimum and maximum length of the roughness, μm. It should be noted that Eq. (16) means that lmin/lmax<0.01.Based on Hooke’s law, the deformation of the “nails” due to the effective stress can be derived. With the detailed derivations stated in AppendixA for a given displacement x of the “nails,” the effective stress exerting on the fracture surfaces can be expressed as follows [11]:
(17)p=FA=EπDRω0−βminDRAε2⋅x−ω0ω01−DR−ω0−x1−DR1−DR+ω02−DR−ω0−x2−DR2−DR=Eφn2−DRlmax2−DR−lmin2−DR⋅x−ω0ω01−DR−ω0−x1−DR1−DR+ω02−DR−ω0−x2−DR2−DR,where p is the effective stress, MPa; A is the area of fracture surface, μm2; E is the rock elastic modulus, GPa; ω0 is the initial fracture aperture, μm; βmin is the minimum shortness of the roughness, μm; F is the force exerting on the fracture, MPa; ε is a constant ratio of the length to the radius of the elements.Based on cubic law, initial permeabilityKf0 and stress-dependent permeability Kf can be obtained, respectively
(18)Kf=112ω0−xR2,(19)Kf0=ω0212.


### 2.4. The Process of Model Determination

Considering the influence brought by the fracture system, the procedure of productivity prediction process of fractured-vuggy reservoir can be summarized as follows:(1)
Based on the initial permeabilityKf0, the initial fracture aperture can be determined by Eq. (19)(2)
With the selected microstructure parameters of fracture, includingφn, lmax, and lmin, the fractal dimension DR can be determined by Eq. (16)(3)
By assigning values toE and xR, the effective stress is determined by Eq. (17), and the permeability of rough fracture system affected by effective stress can be obtained with Eqs. (17) and (18)(4)
By comparing the experimental data or numerical simulation data, steps 2 and 3 are repeated to finally determine the fracture surface microstructure parameters, rock lithology parameters, and fracture system permeability under effective stress(5)
With the determined fracture system permeability under effective stress, the productivity of reservoir in pseudosteady state can be predicted with Eq. (14)


## 2.1. Model Assumptions

In this work, the fractured-vuggy carbonate reservoir is simplified as the triple medium shown in Figure2(a). It is assumed that the reservoir is isotropic, and it is consist of evenly distributed matrix, fracture systems, and vug systems. The single-phase fluid in the reservoir is microcompressible, and the flow of which obeys the Darcy’s law. The flow is isothermal during the whole process, and it has reached a pseudosteady state. The fluid supply to the wellbore is stable, and only radial flow in horizontal direction is considered. Figure 2(b) depicts the whole flow process: the fluid flows into the wellbore only through the fracture systems, and there exists crossflow between matrix fracture and vug fracture. That is, both matrix and vug systems will supply fluid to fracture systems.Figure 2
Sketch of the physical model.
(a)
Simplified triple media model(b)
The connection sketch of triple media


## 2.2. Flow Equations and the Solutions

Based on the assumptions stated above, the equations for describing fluid flow in fractured-vuggy reservoir are as follows [36–38]:Fracture system:(1)Kfμ∂2pf∂r2+1r∂pf∂r=αmKmμpf−pm+αvKvμpf−pv+φfCf∂pf∂t.Matrix system:(2)φmCm∂pm∂t=αmKmμpf−pm.Vug system:(3)φvCv∂pv∂t=αvKvμpf−pv,where μ is the fluid viscosity, mPa·s; p is the formation pressure, MPa; K is the permeability, μm2; α is the shape factor, m-2; t is the production time, d; C is the compressibility coefficient, MPa-1; φ is porosity, %; subscripts f, v, and m, respectively, stand for fractures, vugs, and matrix; r is the distance, m.As aforementioned, during pseudosteady state flow region, the pressure drop are stable and consistent throughout the reservoir, and the volumetric compression of the reservoir is the dominated driving force for fluid flow [31–36]. Thus, for fluid flow in fractured-vuggy reservoirs, we have [30, 36]
(4)∂pf∂t=∂pm∂t=∂pv∂t=−qscBVBφfCf+φmCm+φvCv,where B is the fluid bulk coefficient, m3/m3; VB is the fluid volume under formation condition, m3; qsc is the productivity at ground standard condition, m3/d.Substituting Eq. (4) into Eq. (1), we have
(5)Kfμ∂2pf∂r2+1r∂pf∂r=qscBVB.The boundary conditions can be described as(6)pfr=rw=pwf,(7)r∂pf∂rr=re=0,where pwf is the bottom hole pressure, MPa; re is the radius of the reservoir, m; rw is the wellbore radius, m.By solving Eqs. (5)–(7), the pressure at the fracture system is [36]
(8)pfr,t=pwft+qscμB2KfVBre2lnrrw−0.5r2−rw2,VB=πHre2−rw2,where H is the reservoir effective sickness, m.Mathematically speaking, for triple media with multiple systems, the average reservoir pressure can be regarded as the mean of average pressure in each system. Therefore, it is crucial to determine the weighting coefficient of each system. In general, the weighting coefficient is related to pore volume of each system. However, as the fluid flow during pseudosteady state is driven by the volumetric compression of the reservoir, the compressibility coefficient of each system will also greatly affect the weighting coefficient. Since the compressibility coefficient varies with different systems, only pore volume parameter is not enough to be weighting parameters for accurately determining the average reservoir pressure. As stated in the literature, elastic storage capacity ratio of each system is more reasonable to be the weighting parameters [36, 39, 40]. That is, the average reservoir pressure can be obtained as
(9)p¯=p¯fωf+p¯mωm+p¯vωv,where p¯f is the average pressure in the fracture system, p¯m is the average pressure in the matrix system, p¯v is the average pressure in vug system, and the weighting parameter ω is assigned as the elastic storage capacity ratio, which can be defined as
(10)ωi=φiCiφfCf+φvCv+φmCmi=f,v,m.With the assumption that parameterre is much larger than rw, based on Eq. (8), the average pressure in fracture system can be written as
(11)p¯f=∰VBpfdVBπHre2−rw2=pwft+qscμB2πKfHlnrerw−34.Then, based on Eqs. (1) and (3), the average pressure in the matrix system and vug system can be, respectively, obtained as
(12)p¯m=p¯f+qscωmμBrw2λmVBKf,(13)p¯v=p¯f+qscωvμBrw2λvVBKf,where λm is the matrix-fracture interporosity flow coefficient; λv is the vug-fracture interporosity flow coefficient.According to the definition of well production index (PI), from Eqs. (9) to (13), we have [36]
(14)J=qscp¯−pwf=2πKfHμBlnre/rw−3/4⋅D,where
(15)D=11+2rw/re2⋅ωm2/λm+ωv2/λv/lnre/rw−3/4.Eq. (14) demonstrates that, compared with the normal well production index (PI) model for homogeneous single continuum media, the derived well production index (PI) can be treated as two parts. The first part is Darcy’s PI part, which is identical to the normal well production index (PI) model. Physically speaking, this part takes into account how the flow rate is related to pressure drop for the fluid flow in the continuum media. The second part, parameter D, represents how easy the fluid can move from one system to the others, which is dependent on the crossflow parameters.


## 2.3. Stress Sensitivity of the Fracture System

Eq. (14) demonstrates that the permeability of the fracture system plays an important role in the productivity prediction of the fractured-vuggy reservoir in the pseudosteady state stage. However, as stated previously, effective stress will greatly affect the permeability of the fracture system. Thus, it is significant to quantify stress-dependent permeability of the fracture system and couple it with Eq. (14) to study well productivity more accurately. Due to complex morphology of the fracture system, it is difficult to characterize fracture microstructure and establish the quantitative relation between effective stress and fracture permeability. Fortunately, through a series research work, many researchers have found that the fractal geometry theory was applicable in describing the microstructure of the fracture system [11, 13, 41–45]. Specifically, fractal theory has been proved to be applicable and reliable in describing the complex fracture system. For example, many scholars suggested that rough surfaces of fractures had fractal characteristics [41–45]. In this work, for the sake of simplification, fractal theory will be applied to characterize fracture morphology, and then, an analytical model of stress-dependent permeability of the fracture system based on fractal theory and Gangi’s model will be introduced [11, 19, 44].Based on Gangi’s model, roughness of fracture surface can be regarded as “nails” embedded on the smooth fracture walls [11, 19]. Then, for fracture walls distributed with “nails” with various lengths, we assume that the length of “nails” follows fractal theory. The fractal dimension of the fracture system roughness can be expressed as [11]
(16)DR=2−lnφnlnlmin/lmax,where DR is the fractal dimension for fracture surface roughness (0<DR<2), dimensionless. Area ratio of fracture surface φn is the ratio of total roughness area to total surface area of rough fracture, dimensionless; lmin and lmax, respectively, stand for the minimum and maximum length of the roughness, μm. It should be noted that Eq. (16) means that lmin/lmax<0.01.Based on Hooke’s law, the deformation of the “nails” due to the effective stress can be derived. With the detailed derivations stated in AppendixA for a given displacement x of the “nails,” the effective stress exerting on the fracture surfaces can be expressed as follows [11]:
(17)p=FA=EπDRω0−βminDRAε2⋅x−ω0ω01−DR−ω0−x1−DR1−DR+ω02−DR−ω0−x2−DR2−DR=Eφn2−DRlmax2−DR−lmin2−DR⋅x−ω0ω01−DR−ω0−x1−DR1−DR+ω02−DR−ω0−x2−DR2−DR,where p is the effective stress, MPa; A is the area of fracture surface, μm2; E is the rock elastic modulus, GPa; ω0 is the initial fracture aperture, μm; βmin is the minimum shortness of the roughness, μm; F is the force exerting on the fracture, MPa; ε is a constant ratio of the length to the radius of the elements.Based on cubic law, initial permeabilityKf0 and stress-dependent permeability Kf can be obtained, respectively
(18)Kf=112ω0−xR2,(19)Kf0=ω0212.


## 2.4. The Process of Model Determination

Considering the influence brought by the fracture system, the procedure of productivity prediction process of fractured-vuggy reservoir can be summarized as follows:(1)
Based on the initial permeabilityKf0, the initial fracture aperture can be determined by Eq. (19)(2)
With the selected microstructure parameters of fracture, includingφn, lmax, and lmin, the fractal dimension DR can be determined by Eq. (16)(3)
By assigning values toE and xR, the effective stress is determined by Eq. (17), and the permeability of rough fracture system affected by effective stress can be obtained with Eqs. (17) and (18)(4)
By comparing the experimental data or numerical simulation data, steps 2 and 3 are repeated to finally determine the fracture surface microstructure parameters, rock lithology parameters, and fracture system permeability under effective stress(5)
With the determined fracture system permeability under effective stress, the productivity of reservoir in pseudosteady state can be predicted with Eq. (14)


## 3. Results and Discussions

### 3.1. Model Validation

In order to show it is reasonable to use the proposed fractal model to quantitatively present the stress sensitivity of the fracture system, and to show how to obtain fracture parameters using the model, core slices of 2 cores sampled from two production wells in Jianghan Basin, China, were used for verification and illustration [11]. Experiments were conducted under formation condition using the proppants with a diameter range of 30-60 mesh. It can be seen from Figure 3. that the experimental data of both the two cores showed a friendly match with the calculated data using the fractal model, which means the quantitative model for stress sensitivity of the fracture system is reasonable. Besides, with the stress-dependent permeability model and the experimental data, we successfully obtained the fracture parameters of these two core samples. For example, for core sample 1, the initial fracture aperture was 5.2×10−3cm; the minimum and maximum length of the roughness lmin and lmax were 2.5×10−5cm and 5.2×10−3cm, respectively; the area ratio φn was 0.5%; the fractal dimension for fracture surface roughness DR was 1.01; and the elastic modulus E was 45 GPa. However, for core sample 2, the initial fracture aperture was 5.2×10−3cm; the minimum and maximum length of the roughness lmin and lmax were 2.2×10−5cm and 5.2×10−3cm, respectively; the area ratio φn was 1.2%; the fractal dimension for fracture surface roughness DR was 1.18; and the elastic modulus E was 45 GPa. It should be noted that the value of elastic modulus in the model is identical to that determined from the test results, which validates the rationality of the parameters used in the model.Figure 3
Comparison of stress-dependent permeability of fractures from experimental data and results from the derived model.Lu and Ghedan [46] stated that, for a fully penetrating vertical well located at the center of a close homogeneous, isotropic circular reservoir, the production during pseudosteady state could be written as
(20)qsc=2πK¯Hp¯−pwfμBlnre/rw−3/4,where parameter K¯ is the average permeability of the reservoir. To further validate the proposed model Eq. (13), we assume that the reservoir is just composed of fracture system. In other words, the matrix system and vug system can be neglected in the reservoir (i.e., ωv=ωm=0); then, the model Eq. (14) can be simplified as
(21)qsc=2πKfHp¯−pwfμBlnre/rw−3/4.As the reservoir only contains the fracture system, the average permeability of the reservoir is the permeability of the fracture system. That is, the average permeability of the reservoir isKf. Thus, Eq. (21) is identical to Eq. (20), which validates the rationality of Eq. (14).Recently, Youssef [47] derived an analytical model to study IPR curve of fractured-vuggy carbonate reservoirs, which was
(22)J=qscp¯−pwf=2πKf+KmHE23μBlnre/rw−3/4;1E23=1+2rw/re2lnre/rw−3/4.ωv2λv+λvm+ϑ2γ1+rwre2κ1−κγ⋅K1βreI1βrw−K1βrwI1βreK1βreI0βrw+K0βrwI1βre,ϑ=ωm+ωvλvmλv+λvm+κ−1;κ=KfKf+Km,γ=λm+λv⋅λvmλv+λvm;β=1rwγκ1−κ.It should be noted that the derived model Eq. (22) is based on the assumption that both the fracture system and matrix system are flow paths for fluid flow in the carbonate reservoirs, and the vug system has direct connections with other two systems. If we assume Kf is much larger than Km, ignore the crossflow between the matrix system and vug system (i.e., ωvm=0), then Eq. (22) can be rewritten as
(23)qsc=2πKf+KmHE23p¯−pwfμBlnre/rw−3/4≈2πKfHE23p¯−pwfμBlnre/rw−3/4,1E23=1+2rw/re2lnre/rw−3/4⋅ωv2λv+ωm2λm.By comparing Eqs. (14), (15), and (23), we can see that parameter qsc in Eq. (14) is identical to that in Eq. (23). Besides, parameter D in Eq. (15) is the same with parameter E23 in Eq. (23). Thus, our model Eq. (14) can be further validated by Eqs. (22) and (23). By comparing the predictions from the Eq. (22) and numerical solution and analytical solutions of the fluid flow equations of fractured-vuggy media, Youssef [47] validated the model Eq. (22). Thus, to make our paper more concise, we will not compare the results from our model Eq. (14) with results from numerical simulation and other analytical solutions. However, as the predictions from Eq. (22) are agreement with the results from other models, the predicted results from our model Eq. (14) are surely consistent with those from other models.


### 3.2. Influence of Effective Stress on Well Productivity

Based on Eqs. (18) and (19), the initial fracture permeability and fracture permeability under effective stress will be predicted. Then, the determined fracture permeabilities from these two equations will be, respectively, substituted into Eq. (14) to contrastively analyze the IPR (inflow performance relationship) curves under different stress conditions. The predicted results (IPR curves) and the corresponding parameters used in the models are shown in Figure 4. It can be seen from Figure 4 that, under the effective stress condition, the well production is significantly lower than that without the effective stress, indicating that stress sensitivity of fracture system will greatly affect the well productivity in the pseudosteady state stage. As the bottom hole pressure decreases (or draw-down pressure increases), the difference between these two curves increases. This phenomenon indicates that, as the development progresses, the influence of stress sensitivity of fracture system on productivity gradually increases. Thus, to make the prediction more accurate, it is recommended to take the stress sensitivity of the fracture system into account.Figure 4
The impact of effective stress on IPR curves.


### 3.3. Parameter Sensitivity Analysis

#### 3.3.1. Elastic Storage Capacity Ratio of Matrix

The IPR curves with different matrix elastic storage capacity ratios are shown in Figure5. It shows that the production decreases slightly with the increase of matrix elastic storage capacity ratio. A possible explanation for this is that the matrix system elastic storage capacity ratio reflects its relative storage capacity compared with those of fracture system and vug system. Mathematically speaking, under a given vug elastic storage capacity, a higher matrix elastic storage capacity means a relative lower fracture storage capacity. Since the fracture is the main flow channel for fluid flow, lower fracture storage capacity will cause a decrease in the efficiency of fluid supply to the wellbore.Figure 5
IPR curves under differentωm.


#### 3.3.2. Elastic Storage Capacity Ratio of Vug

The IPR curves with different elastic storage capacity ratios of vug are shown in Figure6. It shows that the productivity decreases largely with the increase of the elastic storage capacity ratio of vug. Generally speaking, the elastic storage capacity ratio of vug reflects its larger relative storage capacity compared with those of the fracture system and matrix system. That is, under a given matrix elastic storage capacity, a higher value of elastic storage capacity ratio of the vug system is always accompanied by a lower value of fracture storage capacity. Thus, there exists a negative relationship between the productivity and the elastic storage capacity ratio of the vug system.Figure 6
IPR curves under differentωv.


#### 3.3.3. Matrix-Fracture Interporosity Flow Coefficient

The IPR curves with different matrix-fracture interporosity flow coefficients are shown in Figure7. It shows that the productivity increases with the increase of the matrix-fracture interporosity flow coefficient. This can be explained by the fact that a larger matrix-fracture interporosity flow coefficient means a stronger flow capacity from the matrix into the fracture, which can improve the fluid supply capacity of fracture system to the wellbore, thereby increasing the productivity.Figure 7
IPR curves under differentλm.


#### 3.3.4. Vug-Fracture Interporosity Flow Coefficient

The IPR curves with different vug-fracture interporosity flow coefficients are shown in Figure8. It can be seen from Figure 8 that, when other parameters are fixed, the productivity markedly increases with the increase of the vug-fracture interporosity flow coefficient. The main reason is that the vug-fracture interporosity flow coefficient indicates the capacity of fluid in the vug system flowing into the fracture system. Physically speaking, under the same condition, the total volume of fluid in the fracture system will increase with increasing vug-fracture interporosity flow coefficient. Since the fluid flows into the wellbore only through the fracture systems, a larger volume of fluid in the fracture system can directly cause an increase in the fluid supply to the wellbore, which also means an increase in the productivity.Figure 8
IPR curves under differentλv.


#### 3.3.5. Fractal Dimension of Fracture Surface Roughness

The IPR curves with different fractal dimension of fracture surface roughness are shown in Figure9. It shows that the productivity decreases with the increase of the fractal dimension of fracture surface roughness. The main reason is that the larger the fractal dimension of fracture surface roughness is, the larger the proportion of small sized roughness is. As a result, with a certain effective pressure, a larger value of the fractal dimension corresponds a larger change in fracture aperture and a larger extent of the fracture closure. Consequently, the IPR curve with a larger value of fractal dimension of fracture surface roughness is lower than that with a smaller value of fractal dimension of fracture surface roughness.Figure 9
IPR curves under differentDR.


#### 3.3.6. Rock Elastic Modulus

The IPR curves with different rock elastic modulus are shown in Figure10. Results (Figure 10) suggest that there exists a positive relationship between the well productivity and rock elastic modulus. That is, the productivity shows a noticeable increase with the increase of the elastic modulus. In general, under a certain effective pressure, the reservoir rock is less prone to deform with lager value of elastic modulus. In other words, under a certain effective pressure, a lager value of elastic modulus corresponds to a smaller value of change in fracture aperture and a larger permeability of fracture. Thus, to some extent, the stress sensitivity of fractured-vuggy reservoir can be ignored, if rock elastic modulus of the reservoir is large enough. However, since there are a number of clay minerals with small elastic modulus exist in the fractured-vuggy carbonate reservoirs, rock elastic modulus of the reservoir cannot be large enough, and the stress sensitivity of these type reservoirs cannot be ignored.Figure 10
IPR curves under differentE.


#### 3.3.7. Area Ratio of Fracture Surface

As the definition of area ratio of fracture surface (i.e., the ratio of total roughness area to total surface area of rough fracture) stated before, this parameter can be used to characterize the microstructure of fracture surface. Figure11 shows the influence of area ratio of fracture surface on IPR curves. It can be seen that the productivity increases with the increasing area ratio of fracture surface. Physically speaking, a larger value of area ratio of fracture surface corresponds to a smaller the proportion of small sized roughness. Thus, under a given overburden pressure (or effective pressure), with the increase of area ratio of fracture surface, the move distance of fracture surfaces decreases, resulting in a smaller extent of the fracture closure and a larger permeability of fracture system. As a result, IPR curve with a larger value of area ratio of fracture surface is higher than that with a smaller value of area ratio of fracture surface.Figure 11
IPR curves under differentφn.


## 3.1. Model Validation

In order to show it is reasonable to use the proposed fractal model to quantitatively present the stress sensitivity of the fracture system, and to show how to obtain fracture parameters using the model, core slices of 2 cores sampled from two production wells in Jianghan Basin, China, were used for verification and illustration [11]. Experiments were conducted under formation condition using the proppants with a diameter range of 30-60 mesh. It can be seen from Figure 3. that the experimental data of both the two cores showed a friendly match with the calculated data using the fractal model, which means the quantitative model for stress sensitivity of the fracture system is reasonable. Besides, with the stress-dependent permeability model and the experimental data, we successfully obtained the fracture parameters of these two core samples. For example, for core sample 1, the initial fracture aperture was 5.2×10−3cm; the minimum and maximum length of the roughness lmin and lmax were 2.5×10−5cm and 5.2×10−3cm, respectively; the area ratio φn was 0.5%; the fractal dimension for fracture surface roughness DR was 1.01; and the elastic modulus E was 45 GPa. However, for core sample 2, the initial fracture aperture was 5.2×10−3cm; the minimum and maximum length of the roughness lmin and lmax were 2.2×10−5cm and 5.2×10−3cm, respectively; the area ratio φn was 1.2%; the fractal dimension for fracture surface roughness DR was 1.18; and the elastic modulus E was 45 GPa. It should be noted that the value of elastic modulus in the model is identical to that determined from the test results, which validates the rationality of the parameters used in the model.Figure 3
Comparison of stress-dependent permeability of fractures from experimental data and results from the derived model.Lu and Ghedan [46] stated that, for a fully penetrating vertical well located at the center of a close homogeneous, isotropic circular reservoir, the production during pseudosteady state could be written as
(20)qsc=2πK¯Hp¯−pwfμBlnre/rw−3/4,where parameter K¯ is the average permeability of the reservoir. To further validate the proposed model Eq. (13), we assume that the reservoir is just composed of fracture system. In other words, the matrix system and vug system can be neglected in the reservoir (i.e., ωv=ωm=0); then, the model Eq. (14) can be simplified as
(21)qsc=2πKfHp¯−pwfμBlnre/rw−3/4.As the reservoir only contains the fracture system, the average permeability of the reservoir is the permeability of the fracture system. That is, the average permeability of the reservoir isKf. Thus, Eq. (21) is identical to Eq. (20), which validates the rationality of Eq. (14).Recently, Youssef [47] derived an analytical model to study IPR curve of fractured-vuggy carbonate reservoirs, which was
(22)J=qscp¯−pwf=2πKf+KmHE23μBlnre/rw−3/4;1E23=1+2rw/re2lnre/rw−3/4.ωv2λv+λvm+ϑ2γ1+rwre2κ1−κγ⋅K1βreI1βrw−K1βrwI1βreK1βreI0βrw+K0βrwI1βre,ϑ=ωm+ωvλvmλv+λvm+κ−1;κ=KfKf+Km,γ=λm+λv⋅λvmλv+λvm;β=1rwγκ1−κ.It should be noted that the derived model Eq. (22) is based on the assumption that both the fracture system and matrix system are flow paths for fluid flow in the carbonate reservoirs, and the vug system has direct connections with other two systems. If we assume Kf is much larger than Km, ignore the crossflow between the matrix system and vug system (i.e., ωvm=0), then Eq. (22) can be rewritten as
(23)qsc=2πKf+KmHE23p¯−pwfμBlnre/rw−3/4≈2πKfHE23p¯−pwfμBlnre/rw−3/4,1E23=1+2rw/re2lnre/rw−3/4⋅ωv2λv+ωm2λm.By comparing Eqs. (14), (15), and (23), we can see that parameter qsc in Eq. (14) is identical to that in Eq. (23). Besides, parameter D in Eq. (15) is the same with parameter E23 in Eq. (23). Thus, our model Eq. (14) can be further validated by Eqs. (22) and (23). By comparing the predictions from the Eq. (22) and numerical solution and analytical solutions of the fluid flow equations of fractured-vuggy media, Youssef [47] validated the model Eq. (22). Thus, to make our paper more concise, we will not compare the results from our model Eq. (14) with results from numerical simulation and other analytical solutions. However, as the predictions from Eq. (22) are agreement with the results from other models, the predicted results from our model Eq. (14) are surely consistent with those from other models.


## 3.2. Influence of Effective Stress on Well Productivity

Based on Eqs. (18) and (19), the initial fracture permeability and fracture permeability under effective stress will be predicted. Then, the determined fracture permeabilities from these two equations will be, respectively, substituted into Eq. (14) to contrastively analyze the IPR (inflow performance relationship) curves under different stress conditions. The predicted results (IPR curves) and the corresponding parameters used in the models are shown in Figure 4. It can be seen from Figure 4 that, under the effective stress condition, the well production is significantly lower than that without the effective stress, indicating that stress sensitivity of fracture system will greatly affect the well productivity in the pseudosteady state stage. As the bottom hole pressure decreases (or draw-down pressure increases), the difference between these two curves increases. This phenomenon indicates that, as the development progresses, the influence of stress sensitivity of fracture system on productivity gradually increases. Thus, to make the prediction more accurate, it is recommended to take the stress sensitivity of the fracture system into account.Figure 4
The impact of effective stress on IPR curves.


## 3.3. Parameter Sensitivity Analysis

### 3.3.1. Elastic Storage Capacity Ratio of Matrix

The IPR curves with different matrix elastic storage capacity ratios are shown in Figure5. It shows that the production decreases slightly with the increase of matrix elastic storage capacity ratio. A possible explanation for this is that the matrix system elastic storage capacity ratio reflects its relative storage capacity compared with those of fracture system and vug system. Mathematically speaking, under a given vug elastic storage capacity, a higher matrix elastic storage capacity means a relative lower fracture storage capacity. Since the fracture is the main flow channel for fluid flow, lower fracture storage capacity will cause a decrease in the efficiency of fluid supply to the wellbore.Figure 5
IPR curves under differentωm.


### 3.3.2. Elastic Storage Capacity Ratio of Vug

The IPR curves with different elastic storage capacity ratios of vug are shown in Figure6. It shows that the productivity decreases largely with the increase of the elastic storage capacity ratio of vug. Generally speaking, the elastic storage capacity ratio of vug reflects its larger relative storage capacity compared with those of the fracture system and matrix system. That is, under a given matrix elastic storage capacity, a higher value of elastic storage capacity ratio of the vug system is always accompanied by a lower value of fracture storage capacity. Thus, there exists a negative relationship between the productivity and the elastic storage capacity ratio of the vug system.Figure 6
IPR curves under differentωv.


### 3.3.3. Matrix-Fracture Interporosity Flow Coefficient

The IPR curves with different matrix-fracture interporosity flow coefficients are shown in Figure7. It shows that the productivity increases with the increase of the matrix-fracture interporosity flow coefficient. This can be explained by the fact that a larger matrix-fracture interporosity flow coefficient means a stronger flow capacity from the matrix into the fracture, which can improve the fluid supply capacity of fracture system to the wellbore, thereby increasing the productivity.Figure 7
IPR curves under differentλm.


### 3.3.4. Vug-Fracture Interporosity Flow Coefficient

The IPR curves with different vug-fracture interporosity flow coefficients are shown in Figure8. It can be seen from Figure 8 that, when other parameters are fixed, the productivity markedly increases with the increase of the vug-fracture interporosity flow coefficient. The main reason is that the vug-fracture interporosity flow coefficient indicates the capacity of fluid in the vug system flowing into the fracture system. Physically speaking, under the same condition, the total volume of fluid in the fracture system will increase with increasing vug-fracture interporosity flow coefficient. Since the fluid flows into the wellbore only through the fracture systems, a larger volume of fluid in the fracture system can directly cause an increase in the fluid supply to the wellbore, which also means an increase in the productivity.Figure 8
IPR curves under differentλv.


### 3.3.5. Fractal Dimension of Fracture Surface Roughness

The IPR curves with different fractal dimension of fracture surface roughness are shown in Figure9. It shows that the productivity decreases with the increase of the fractal dimension of fracture surface roughness. The main reason is that the larger the fractal dimension of fracture surface roughness is, the larger the proportion of small sized roughness is. As a result, with a certain effective pressure, a larger value of the fractal dimension corresponds a larger change in fracture aperture and a larger extent of the fracture closure. Consequently, the IPR curve with a larger value of fractal dimension of fracture surface roughness is lower than that with a smaller value of fractal dimension of fracture surface roughness.Figure 9
IPR curves under differentDR.


### 3.3.6. Rock Elastic Modulus

The IPR curves with different rock elastic modulus are shown in Figure10. Results (Figure 10) suggest that there exists a positive relationship between the well productivity and rock elastic modulus. That is, the productivity shows a noticeable increase with the increase of the elastic modulus. In general, under a certain effective pressure, the reservoir rock is less prone to deform with lager value of elastic modulus. In other words, under a certain effective pressure, a lager value of elastic modulus corresponds to a smaller value of change in fracture aperture and a larger permeability of fracture. Thus, to some extent, the stress sensitivity of fractured-vuggy reservoir can be ignored, if rock elastic modulus of the reservoir is large enough. However, since there are a number of clay minerals with small elastic modulus exist in the fractured-vuggy carbonate reservoirs, rock elastic modulus of the reservoir cannot be large enough, and the stress sensitivity of these type reservoirs cannot be ignored.Figure 10
IPR curves under differentE.


### 3.3.7. Area Ratio of Fracture Surface

As the definition of area ratio of fracture surface (i.e., the ratio of total roughness area to total surface area of rough fracture) stated before, this parameter can be used to characterize the microstructure of fracture surface. Figure11 shows the influence of area ratio of fracture surface on IPR curves. It can be seen that the productivity increases with the increasing area ratio of fracture surface. Physically speaking, a larger value of area ratio of fracture surface corresponds to a smaller the proportion of small sized roughness. Thus, under a given overburden pressure (or effective pressure), with the increase of area ratio of fracture surface, the move distance of fracture surfaces decreases, resulting in a smaller extent of the fracture closure and a larger permeability of fracture system. As a result, IPR curve with a larger value of area ratio of fracture surface is higher than that with a smaller value of area ratio of fracture surface.Figure 11
IPR curves under differentφn.


## 3.3.1. Elastic Storage Capacity Ratio of Matrix

The IPR curves with different matrix elastic storage capacity ratios are shown in Figure5. It shows that the production decreases slightly with the increase of matrix elastic storage capacity ratio. A possible explanation for this is that the matrix system elastic storage capacity ratio reflects its relative storage capacity compared with those of fracture system and vug system. Mathematically speaking, under a given vug elastic storage capacity, a higher matrix elastic storage capacity means a relative lower fracture storage capacity. Since the fracture is the main flow channel for fluid flow, lower fracture storage capacity will cause a decrease in the efficiency of fluid supply to the wellbore.Figure 5
IPR curves under differentωm.


## 3.3.2. Elastic Storage Capacity Ratio of Vug

The IPR curves with different elastic storage capacity ratios of vug are shown in Figure6. It shows that the productivity decreases largely with the increase of the elastic storage capacity ratio of vug. Generally speaking, the elastic storage capacity ratio of vug reflects its larger relative storage capacity compared with those of the fracture system and matrix system. That is, under a given matrix elastic storage capacity, a higher value of elastic storage capacity ratio of the vug system is always accompanied by a lower value of fracture storage capacity. Thus, there exists a negative relationship between the productivity and the elastic storage capacity ratio of the vug system.Figure 6
IPR curves under differentωv.


## 3.3.3. Matrix-Fracture Interporosity Flow Coefficient

The IPR curves with different matrix-fracture interporosity flow coefficients are shown in Figure7. It shows that the productivity increases with the increase of the matrix-fracture interporosity flow coefficient. This can be explained by the fact that a larger matrix-fracture interporosity flow coefficient means a stronger flow capacity from the matrix into the fracture, which can improve the fluid supply capacity of fracture system to the wellbore, thereby increasing the productivity.Figure 7
IPR curves under differentλm.


## 3.3.4. Vug-Fracture Interporosity Flow Coefficient

The IPR curves with different vug-fracture interporosity flow coefficients are shown in Figure8. It can be seen from Figure 8 that, when other parameters are fixed, the productivity markedly increases with the increase of the vug-fracture interporosity flow coefficient. The main reason is that the vug-fracture interporosity flow coefficient indicates the capacity of fluid in the vug system flowing into the fracture system. Physically speaking, under the same condition, the total volume of fluid in the fracture system will increase with increasing vug-fracture interporosity flow coefficient. Since the fluid flows into the wellbore only through the fracture systems, a larger volume of fluid in the fracture system can directly cause an increase in the fluid supply to the wellbore, which also means an increase in the productivity.Figure 8
IPR curves under differentλv.


## 3.3.5. Fractal Dimension of Fracture Surface Roughness

The IPR curves with different fractal dimension of fracture surface roughness are shown in Figure9. It shows that the productivity decreases with the increase of the fractal dimension of fracture surface roughness. The main reason is that the larger the fractal dimension of fracture surface roughness is, the larger the proportion of small sized roughness is. As a result, with a certain effective pressure, a larger value of the fractal dimension corresponds a larger change in fracture aperture and a larger extent of the fracture closure. Consequently, the IPR curve with a larger value of fractal dimension of fracture surface roughness is lower than that with a smaller value of fractal dimension of fracture surface roughness.Figure 9
IPR curves under differentDR.


## 3.3.6. Rock Elastic Modulus

The IPR curves with different rock elastic modulus are shown in Figure10. Results (Figure 10) suggest that there exists a positive relationship between the well productivity and rock elastic modulus. That is, the productivity shows a noticeable increase with the increase of the elastic modulus. In general, under a certain effective pressure, the reservoir rock is less prone to deform with lager value of elastic modulus. In other words, under a certain effective pressure, a lager value of elastic modulus corresponds to a smaller value of change in fracture aperture and a larger permeability of fracture. Thus, to some extent, the stress sensitivity of fractured-vuggy reservoir can be ignored, if rock elastic modulus of the reservoir is large enough. However, since there are a number of clay minerals with small elastic modulus exist in the fractured-vuggy carbonate reservoirs, rock elastic modulus of the reservoir cannot be large enough, and the stress sensitivity of these type reservoirs cannot be ignored.Figure 10
IPR curves under differentE.


## 3.3.7. Area Ratio of Fracture Surface

As the definition of area ratio of fracture surface (i.e., the ratio of total roughness area to total surface area of rough fracture) stated before, this parameter can be used to characterize the microstructure of fracture surface. Figure11 shows the influence of area ratio of fracture surface on IPR curves. It can be seen that the productivity increases with the increasing area ratio of fracture surface. Physically speaking, a larger value of area ratio of fracture surface corresponds to a smaller the proportion of small sized roughness. Thus, under a given overburden pressure (or effective pressure), with the increase of area ratio of fracture surface, the move distance of fracture surfaces decreases, resulting in a smaller extent of the fracture closure and a larger permeability of fracture system. As a result, IPR curve with a larger value of area ratio of fracture surface is higher than that with a smaller value of area ratio of fracture surface.Figure 11
IPR curves under differentφn.


## 4. Discussions

Eq. (14) illustrates that, except for the parameter D, the IPR curves of fractured-vuggy carbonate reservoirs are identical to the normal Darcy’s IPR curves generated for single continuum reservoir. Thus, if the parameter D is assigned as unity, the IPR curves of fractured-vuggy carbonate reservoirs can be simplified as the normal Darcy’s IPR curves generated for single continuum reservoir. As a result, parameter D can be regarded as a parameter representing of how easy the fluid can move from one system to the others. In this paper, we ignore the fluid flow within the matrix system and assume that only the fracture system is directly connected with the wellbore. However, if the permeability of the matrix system is not small enough, we can modify the model and take the fluid flow in matrix system into account. In that case, we just need to modify parameter D. Physically speaking, parameter D accounts for the pressure drop that occurs due to the crossflow between different systems in the reservoirs. In the normal straight line IPR model, we ignore the crossflow between different systems and only consider the radial flow pressure drop. In other words, when the value of parameter D approaches unity, the IPR curve will go up due to the reduction in the pressure drop caused by crossflow between different systems.As stated by Youssef [47], if we take the effect of damage zone with small radius around wellbore into account, Eq. (14) can be modified as
(24)J=qscp¯−pwf=2πKfHDμBlnre/rw−3/4+s,where parameter s is the skin factor.However, it should be noted that Eq. (24) is based on the assumption that the damaged zone only affects the radial flow and has nothing to do with the crossflow between different systems. That is, the damaged zone has little influence on parameter D. Physically, parameter D will be affected by the damaged zone. In our future work, we will further study the effect of damaged zone on parameter D.


## 5. Summary and Conclusions

In this paper, based on fractal theory, a theoretical model for stress-dependent permeability of rough fracture is developed. Then, by coupling flow equations in triple media, and the derived fracture permeability model, a new productivity prediction model for fractured-vuggy reservoir under pseudosteady state is established. Compared with former models, the derived models in this work account for complex morphology of rough fracture, fracture closure due to effective stress, and crossflow between different systems in the fractured-vuggy carbonate reservoirs. In addition, every model parameter in this work has exact physical meaning. Followings are the main conclusions:(1)
The influence of effective stress on well productivity is significant, and it will become stronger as the development progresses. Thus, our derived model can help to predict the well productivity more accurately with lower the uncertainty(2)
Compared with matrix elastic storage capacity ratio, the elastic storage capacity ratio of vug will have a greater impact on well productivity. It is found that higher elastic storage capacity ratio of vug will cause a decrease in the well production. However, increases in both matrix-fracture interporosity flow coefficient and vug-fracture interporosity flow coefficient will increase the well production, which was anticipated(3)
It is worth noting that the model is established on the assumption of evenly distributed matrix, fracture system, and vug system, which is largely simplified than real reservoir situation. In addition, the model is limited to single phase flow, which ignores the interaction between multiple phases. To better understand the well productivity in fractured-vuggy reservoirs, more works are needed to be done. In our future work, the derived model will be extended to multiple phases flow in the fractured-vuggy reservoirs



---

*Source: 5613966-2021-10-22.xml*