Document Type: Research Paper
Authors
^{1} M.S. Student, Department of Petroleum Engineering, Ahwaz Faculty of Petroleum, Petroleum University of Technology, Ahwaz, Iran
^{2} M.S., Student, Department of Chemical and Petroleum Engineering, Sharif University of Technology, Tehran, Iran
^{3} Professor, Da Prat & Associates (Consulting Company), Buenos Aires, C1420, Argentina
^{4} Professor, Department of Petroleum Engineering, Montanuniversität Leoben, Austria
Abstract
Keywords
Main Subjects
1. Introduction
Naturally fractured reservoirs consist of matrix blocks and fractures and are known as heterogeneous porous media. Matrix blocks have low permeability but contain most of the fluid. On the other hand, the fractures have high permeability but do not store much. Most of the fluid in the reservoir flows through the fractures into the wellbore. Therefore, the matrixfracture transport capacity, namely interporosity flow, governs the producing capacity of these kinds of reservoirs.
The behavior of naturally fractured reservoirs is similar to that of a homogeneous reservoir in the pseudosteady state condition. However, the transient behaviors of these two reservoirs are different due to the double porosity nature of fractured reservoirs. Therefore, we need a powerful method to estimate the properties of naturally fractured reservoirs as accurately as possible. Also, taking the matrix flow into account is essential for estimating different parameters of naturally fractured reservoirs.
Foundations of different models that have been proposed for describing the flow in double porosity systems were laid by Warren and Root (1963) and Barenblatt (1960). In the model of Warren and Root, a cubic matrix with three orthogonal faces with a flow into the fractures under the pseudosteady state condition was assumed. Different solutions were determined for both infinite and finite systems. For the well test analysis of double porosity systems, the slab model exhibiting the matrix transient flow has mostly been used. DeSwaan (1986) developed a mathematical model for the matrix transient flow and used different dimensions of rectangular parallelepipeds for modeling the matrix outflow, but the obtained results never matched the fracture flow for the well test analysis.
Several methods have been presented for well testing of naturally fractured reservoirs, and they are almost based on the pressure transient. These methods can be classified into two main categories; conventional methods, such as the type curve matching and the Horner plot, and the Tiab direct synthesis (TDS) method. In the TDS technique, there is no need to use the type curve matching or the Horner plot, but the loglog plot of pressure and the log derivative data are critical for evaluating reservoir properties (Tiab, 1995; Tiab et al., 2007). The TDS method is applicable to different areas of application, including the analysis of vertically fractured wells in the closed system, naturally fractured reservoirs, and horizontal well tests in naturally fractured reservoirs without the typecurve matching. However, although these methods require timeconsuming calculations, and different plots of pressure and pressure derivatives must be analyzed using different kinds of software, they are nearly accurate techniques. Therefore, proposing a simple and accurate method which can nearly estimate all the parameters of a naturally fractured reservoir and be applicable to reservoirs with different geometries is crucial.
In this study, the rate transient analysis is presented for constant pressure production from closed and bounded naturally fractured reservoirs, which can estimate the properties of these reservoirs accurately without needing to use complicated calculations. To develop this technique, equations which are based on the model of Warren and Root (1963) and the pseudosteady state and constantrate solution given by Mavor and CincoLey (1979) were used for fracture depletion, and the pseudosteady state and constantpressure solution by Da Prat et al (1981) was employed for the matrix depletion.This method employs the ratetime plot utilized in the decline curve analysis. To use this method, measurements of flow rates versus time are needed, so it requires production at constant bottomhole pressure and under the pseudosteady state flow condition. This technique is based on the exponential solution, and by taking the natural logarithm of both sides of the fracture depletion equation and the total system depletion equation straight lines are found for both equations; then, the slope and intercept of the straight lines can be used for estimating storage capacity ratio (ω), fracture storage capacity (ØC_{t})_{f}, reservoir drainage area (A), interporosity flow parameter (λ), block shape factor (α), and fracture permeability (k_{f}) of reservoirs.
2. Method description
Naturally fractured reservoirs are heterogeneous systems with matrix blocks and fractures that are randomly distributed (see Figure 1). These types of systems are modeled by assuming that the orthogonal system of uniform and continuous fractures separates the matrix elements. In this assumption, systems of fractures are oriented parallel to the permeability principal axes.
Figure 1
Naturally fractured reservoirs as heterogeneous systems: a) real fractured rock (after Warren and Root, 1963); b) a schematic of fractured reservoir rock (after Barenblatt and Zheltov, 1960).
The solutions for the bounded, closed system are the main objective of this study. The behavior of ahomogeneous, closedouterboundary reservoirs has been studied by many authors. Fundamental partial differential equations and solutions for rate transient analysis in reservoirs with different geometries can be found elsewhere (Daryasafar et al., 2017).
Da Prat (1981) shows that for closed reservoirs and in rate transient tests, the flow rate shows a rapid decline at first and then levels off for a long period, after which a final rate decline happens again. In this study, we use this fact to propose a simple method to estimate nearly all the properties of naturally fractured reservoirs. In other words, the reason behind the proposed method is that the fracture depletion dominates for shorter times, while, for longer times, the matrix depletion dominates; consequently, each period can be characterized by a different equation.
A closed naturally fractured reservoir shows a rapid decline in the flow rate at first, which is because of the fracture depletion. The constantrate and pseudosteady state solution for the depletion of fracture is defined by (Mavor and CincoLey, 1979):
(1) 
where
(2) 

(3) 
We use Laplace transform to obtain Equation 4:
(4) 
This expression is valid for .
As indicated by Van Everdingen and Hurst (1949), there exists a relationship between the Laplace transformed solutions for the constant pressure and constant rate problems, so if the relation of one of them is known, the other one can be simply derived. Van Everdingen and Hurst showed that q_{D} and p_{D} can be related by the below equation in the Laplace space:
(5) 

(6) 
From Equation 4 and Equation 5, one may obtain:
(7) 
The inverse Laplace transforms of Equation 7 yields:
(8) 
The above equation is used for fracture depletion under the pseudosteady state and constant bottomhole pressure conditions.
Moreover, Da Prat et al. (1981) suggested a solution for the matrix depletion at a constant bottomhole pressure as follows (Da Prat et al., 1981):
(9) 
Taking the natural logarithm of both sides of Equation 8 gives a straight line with the slope of m_{Df} and the intercept of b_{Df} as follows:
(10) 

(11) 
By taking the natural logarithm of both sides of Equation 9, a straight line with the slope of m_{Dm}and the intercept of b_{Dm}can be obtained as follows:
(12) 
and
(13) 
One may obtain the below equations from Equations 10 and 11:
(14) 

(15) 
From Equations 12 and 13 and according to the assumption that , we have:
(16) 
Letting Equation 15 equal Equation 16 gives:
(17) 
The below dimensionless variables are defined for determining different properties of naturally fractured reservoirs. In this method, since the slope and the intercept of the curve of ln(rate) versus flow rate are used, the dimensionless slope and the dimensionless intercept parameters must be defined for simplifying the relations and correlating the fractured parameters; for the sake of simplicity, the evaluation of these parameters is presented elsewhere in detail (Daryasafar et al., 2017).
(18) 

(19) 

(20) 
Now, the obtained slopes and intercepts can be utilized to estimate the parameters of a naturally fractured system.
(21) 
The storage capacity ratio (ω) can be determined by combining Equation 16 into Equation 21:
The dimensionless parameter (ω) defines the ratio of fractures storage to the storage of the total system. Mathematically, it is given by:
(22) 
By using the above equations and given (ϕc_{t})_{m}, the fracture storage capacity (ϕc_{t})_{f} and the total storage capacity (ϕc_{t})_{t} can be estimated:
Using Equations 18 and 19 for the matrixes:
(23) 
Using Equations 18 and 19 for fractures:
(24) 
By combining Equations 17, 23, and 24, the reservoir drainage area can be estimated using the below equation:
(25) 
From Equations 11 and 19, the fracture permeability is defined as:
(26) 
Interporosity flow coefficient (λ) defines the ability of the fluid to flow from the matrix blocks to the fractures system. This parameter relates the block shape factor (α) to the matrix and fracture permeabilities as follows:
(27) 
By combining Equations 13, 19, 27, and 28, the interporosity flow parameter and the block shape factor can de estimated as follows:
(28) 

(29) 
Based on the above formulation and calculations, if a ratetransient well test analysis is performed in a naturally fractured reservoir, the important parameters of the reservoir can be calculated directly. It is worth mentioning that r_{w} in all the above calculations is the effective wellbore radius, which can be determined by r_{w} = r_{wa} e^{–S}, where S is the skin factor and should be estimated by other tests such as the buildup test and r_{wa} is the actual wellbore radius. Needing a short buildup test for the evaluation of the skin factor is the main limitation of the proposed well test approach. The accuracy and the effectiveness of the proposed method are discussed through the following simulated and field test examples.
3. Simulated and field test examples
3.1. Simulated example
A simulated example was carried out for the rate transient well test in the vertical well of a cylindrical naturally fractured reservoir. The data on the flow rate versus time are listed in Table 1. Figure 2 illustrates the variation of flow rate with time, and the change in ln(q) as a function of time is presented in Figure 3. Rock and fluid parameters used in this example are given below.
p_{i} = 11500 psi p_{wf} =5000 psi r'_{w }= 0.25 ft s = –4.09 k_{fi }= 0.147 mD
h = 480 ft Ø_{m }= 10.96% μ = 1 cP c_{tm} = 2.54×10^{–8} psi^{–1} B = 1 RB/STB
k_{m} = 0.1 mD
Solution
As can be seen in Figure 3, data points at early times follow the fracture depletion behavior, and as the time increases, the interaction between fracture and matrix blocks becomes the same as that of homogeneous systems, thereby reflecting the total system depletion. The following data are read from Figure 3:
Storage capacity ratio is estimated using Equation 21:
(30) 
Fracture and total storage capacities are given by:
(31) 
(32) 
From Equation 25, the reservoir drainage area is defined by:
(33) 
The fracture permeability is equal to:
(34) 
From Equation 28, the interporosity flow parameter is expressed by:
(35) 
From Equation 29, the block shape factor is calculated as:
(36) 
The type curve matching method is also used for estimating parameters for this example. The results are presented in Table 2; the method is presented in Appendix A. It should be noticed that in the type curve matching method, ω and λ need to be first obtained by other methods like the buildup analysis, but in the proposed method, just the skin factor needs to be known, which is one of the advantages of our method.
Table 1
Transient rate data for the simulated example.
Q (bbl/day) 
Time (h) 
740 
0.48 
680 
0.72 
620 
0.96 
540 
1.44 
470 
1.92 
360 
2.88 
250 
4.08 
210 
5.04 
160 
6.48 
125 
8.16 
105 
9.6 
88 
12 
79 
16.8 
76 
21.6 
76 
28.8 
74 
48 
Table 2
Comparison between the results of the simulated example obtained by different methods.
Parameter 
Type curve method 
This study 
Error (%) 
ω 
0.01 
0.01 
0 
Drainage area (ft^{2}) 
7068583.47 
7198718.501 
1.8 
(psi^{–1}) 
2.8×10^{–11} 
2.8×10^{–11} 
0 
(psi^{–1}) 
2.8×10^{–9} 
2.8×10^{–9} 
0 
λ 
5×10^{–6} 
4.8×10^{–6} 
–4 
α (ft^{–2}) 
3.33×10^{–8} 
3.12×10^{–8} 
–6.3 
(mD) 
0.15 
0.146 
–2.6 
Figure 2
Variation in oil flow rate versus time for the simulated example.
Figure 3
Variation of ln(flow rate) as a function of time for the simulated example.
3.2. Field test example
Chen (1985) reported the data on a field rate transient test for a fractured system in the Austin Chalk formation as listed in Table 3. Other rock and fluid properties of this formation are also as follows:
p_{i} = 3800 psi p_{wf} =0 psi r’_{w }= 0.25 ft s = –4.0
k_{fi }= 0.82 mD h = 40 ft Ø_{m }= 10% μ = 0.26 cP
C_{tm }= 3.041×10^{–6} psi^{–1} B = 1.58 RB/STB k_{m }= 0.28 mD
Solution
For a closed boundary reservoir, the plot of ln(flow rate) versus time shows a rapid decline at first (fracture depletion), then a curved shape which shows the transient between the fracture depletion and the total system depletion, and finally another straight line for a long period of time, which shows the total (fracture + matrix) depletion. As a result, regions related to the fracture system in pseudosteady state (PSS) flow and the total system in PSS flow can be indicated (see Figure 4). These regions can also be identified using the diagnostic plot as discussed in the previous examples, and then the data points related to these regions must be plotted on the graph of ln(rate) versus time to find the slopes and the intercepts of the curves.
Using Figure 4, the below data can be determined:
Storage capacity ratio is estimated by using Equation 21:
(37) 
Fracture and total storage capacities are given by:
(38) 
(39) 
The reservoir drainage area is obtained from Equation 25:
(40) 
The fracture permeability is equal to:
(41) 
Using Equation 28, the interporosity flow parameter is expressed by:
(42) 
Using Equation 29, the block shape factor is defined as:
(43) 
Figure 4
Variation of ln(flow rate) as a function of time for the field example.
The type curve matching method was used for estimating parameters in this example. The results are presented in Table 4; the method is presented in Appendix B. It should be mentioned that, in the type curve matching method, parameters such as ω and λ must first be obtained by other methods like the buildup analysis and cannot be obtained using type curves; also, the published type curves can be used just for circular reservoirs; nevertheless, in the proposed method, only the skin factor needs to be known, which is one of the advantages of our method.
Table 3
Transient rate data on the field example.
Q (bbl/day) 
Time (h) 
238.86 
1306.4 
248.3 
1941.4 
234.3 
2783.1 
189.37 
3459.3 
153.06 
4241.1 
119.02 
4968.4 
110.1 
5756.4 
89.04 
6327.1 
82.4 
7220.7 
72 
7953.2 
59.3 
8630.4 
57 
9420.2 
54.9 
9998.9 
50.8 
10787 
46.1 
11468.5 
46.1 
12207.2 
42.7 
12995.3 
38.7 
13676.8 
38 
14362 
34.5 
15149 
33.2 
15886 
32.5 
16518 
28.4 
17304 
27.6 
18044 
25.81 
18724 
23.4 
19353 
Table 4
Comparison between the results of the field example obtained by different methods.
Parameter 
Type curve method 
This study 
Error (%) 
ω 
0.5527 (From Chen, 1985) 
0.5527 
0 
Drainage area (ft^{2}) 
1742258.45 
1751764.71 
0.54 
(psi^{–1}) 
3.757×10^{–7} 
3.757×10^{–7} 
0 
(psi^{–1}) 
6.8×10^{–7} 
6.8×10^{–7} 
0 
λ 
4.50×10^{–5 }(From Chen, 1985) 
4.45×10^{–5} 
–1.1 
(mD) 
0.64 
0.72 
12.1 
α (ft^{–2}) 
5.6×10^{–7} 
6.194×10^{–7} 
10.6 
4. Conclusions
A new method was proposed for the accurate estimation of the characteristics of naturally fractured reservoirs using the analysis of the data on transient rates. The analysis combines the fracture and the total system depletion equations to accurately calculate the storage capacity ratio, the drainage area, the fracture storage capacity, the interporosity flow parameter, the block shape factor, the reservoir shape factor, the fracture permeability, and all the other parameters by figuring out the slope and the intercept of the matrix and the fracture depletion straight lines obtained when we plot ln(rate) versus time. It should be pointed out that the slope for the fracture system and that of the total system (matrix and fracture) should be taken in the time interval in which pressure for the constant rate solution exhibits a pseudosteady state flow regime. In this method, no predesigned charts are used, and by using the PSS flow regime, the impact of wellbore storage on the value of different parameters found by this method is negligible. The proposed method was tested using different fields and simulated examples, and its accuracy and efficiency were confirmed by comparing the results of this method with those of the type curve matching approach.
Nomenclature
A 
Drainage area (ft^{2}) 
B 
Oil formation volume factor (RB/STB) 
b 
Intercept of semilog plot (h.^{–1}) 
C 
System compressibility (psi^{–1}) 
C_{A} 
Reservoir shape factor (dimensionless) 
h 
Net pay thickness (ft) 
k 
Permeability (mD) 
m 
Slope of semilog plot (h.^{–1}) 
p 
Pressure (psi) 
q 
Oil flow rate (STB/D) 
r 
Radius (ft) 
r_{D } 
r/r_{w} (dimensionless) 
r_{w} 
Effective well bore radius (ft) 
s 
Laplace variable^{ } 
S 
Skin factor (dimensionless) 
t 
Time (h.) 
Greek 

Ø 
Porosity (dimensionless) 
Block shape factor (1/L^{2}) 

λ 
Interporosity flow parameter (dimensionless) 
μ 
Viscosity (cP) 
ω 
Storage capacity ratio (dimensionless) 
Subscripts 

D 
Dimensionless 
f 
Fracture 
h 
Homogeneous 
i 
Initial 
m 
Matrix 
t 
Total 
wf 
Well bore flowing 
Appendix A
The type curve matching method for estimating reservoir parameters in the simulated example
Step 1: Utilizing this method requires plotting flow rate (bbl/day) versus time (days) on the loglog scale, as shown in Figure A1. The flow rate was plotted on tracing paper and then placed over the type curve corresponding to ω = 0.01 and also λ = 5×10^{–6} , which were already obtained from the buildup analysis. From the match point, parameters can be estimated as follows:
We choose a match point in Figure A2 as follows:
Step 2: From this match point, parameters can be estimated:
(A1)

(A2) 
(A3) 
(A4) 
Appendix B
The type curve matching method for estimating reservoir parameters in the field example
Since there is no published type curve for ω = 0.5527 and λ = 4.5×10^{–5}, we have to find the reservoir parameters using different type curves and calculate the parameters by interpolating and extrapolating between the values obtained via these type curves. For this example, four type curves with 1) ω = 0.1 and λ = 10^{–5}, 2) ω = 0.1 and λ = 10^{–4}, 3) ω = 0.01 and λ = 10^{–5}, and 4) ω = 0.01 and λ = 10^{–4} were used for estimating different parameters.
Step 1:The flow rate was plotted on tracing paper (Figure B1) and then placed over the type curves corresponding to the above values. The results related to the matching points are obtained as follows:
The obtained match point using Figure B2 is expressed by:
The obtained match point using Figure B3 is described by:
The obtained match point using Figure B4 is defined as:
The obtained match point using Figure B5 is equal to:
Step 2: Using the above match points and equation , k_{f} can be calculated for each type curve.
Step 3: From , we obtain r_{w} which is the effective wellbore radius considering the skin effect. Then, by using the obtained r_{eD} for each type curve, r_{e} and the drainage area (A) can be estimated.
Step 4: For each type curve, can be calculated using equation , and then we have by .
Step 5: The block shape factor can be obtained using equation .
Using the above match points to calculate the reservoir parameters for each curve and then by interpolating and extrapolating between the results, the final results of a naturally fractured reservoir with ω = 0.5527 and λ = 4.5×10^{–5} are given by:
Figure A1
Variation of the flow rate of the simulated example versus time on a loglog scale.
Figure A2
Type curve matching for the simulated example.
Figure B1
Variation of the flow rate of the field example versus time on a loglog scale.
Figure B2
Matching the data of the field example on a type curve with ω = 0.1 and λ = 10^{–5}.
Figure B3
Matching the data of the field example on a type curve with ω = 0.1 and λ = 10^{–4}.
Figure B4
Matching the data of the field example on a type curve with ω = 0.01 and λ = 10^{–5}.
^{ }
Figure B5
Matching the data of the field example on a type curve with ω = 0.01 and λ = 10^{–4}.