- Email: [email protected]
Comparison of streamline-based and grid-based dual porosity simulation
Journal of Petroleum Science and Engineering 43 (2004) 129 – 137 www.elsevier.com/locate/petrol
Comparison of streamline-based and grid-based dual porosity simulation Wenfen Huang a,b, Ginevra Di Donato a, Martin J. Blunt a,* a
Department of Earth Science and Engineering, Imperial College, London SW7 2AZ, UK b Shengli Geological Scientific Research Institute, Sinopec Inc., China Received 21 May 2003; accepted 22 January 2004
Abstract In field-scale dual porosity models, fractured reservoirs are represented as two intercommunicating domains: a flowing fraction composed of the connected fracture network and high permeability matrix; and a stagnant fraction composed of lower permeability matrix containing most of the original oil. The transfer of fluids between fracture and matrix is modeled by an empirical function of saturation. Two transfer functions are used in a streamline-based dual porosity code: (1) a conventional transfer function used in traditional grid-based simulators; and (2) a new linear transfer function that reproduces the results of core-scale imbibition experiments. The results of the simulations are compared with those from a conventional grid-based finite difference code. The two codes give very similar results, but the streamline code requires a much shorter run time and has a lower memory requirement. Using the linear transfer function, the streamline-based code is orders of magnitude faster and requires many times less memory than the conventional grid-based code. The larger the number of grid blocks, the greater the relative speed of the streamline-based approach. It is possible to run a 4.5 million grid block dual porosity model on a standard PC using a streamline code. D 2004 Elsevier B.V. All rights reserved. Keywords: Streamlines; Dual porosity; Fractured reservoir; Transfer function
1. Introduction Traditionally, to model a fractured reservoir at the field scale, a flowing region with high permeability represents the fracture network, while a stagnant region with low permeability represents the matrix (Barenblatt et al., 1960; Warren and Root, 1963). A dual porosity (no flow between matrix blocks, as implemented in this paper) or dual permeability (there
* Corresponding author. Tel.: +44-20-7594-6500; fax: +44-207594-7444. E-mail address: [email protected] (M.J. Blunt). 0920-4105/$ - see front matter D 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.petrol.2004.01.002
is flow between matrix blocks) approach is applicable where the majority of the flow is through a wellconnected fracture network (Gilman and Kazemi, 1983; Nakashima et al., 2001). The fluid flow between fracture and matrix is modeled using empirical transfer functions. In most studies a conventional steady-state transfer function is used that varies non-linearly with fracture and matrix saturations. The nonlinearity of the transfer function makes an iterative approach to solving the transport equations necessary, resulting in longer run times than equivalent single porosity (non-fractured) simulations. Streamline-based simulation is frequently used for modeling recovery by displacement because it is
130
W. Huang et al. / Journal of Petroleum Science and Engineering 43 (2004) 129–137
generally faster, suffers from less numerical dispersion and can handle larger models than grid-based codes (Batycky et al., 1997). Streamline methods are ideal for simulating incompressible, or nearly incompressible, systems in highly heterogeneous reservoirs. When the system becomes more compressible and the voidage replacement ratio is much less than one, the method becomes less efficient (Baker et al., 2002). Streamline methods have been successfully applied to several field studies (Samier et al., 2002; Baker et al., 2002; Grinestaff and Caffrey, 2000). Recently, a dual porosity streamline-based code to simulate waterflooding was published (Di Donato et al., 2003). This code will be used to simulate flow in a synthetic fractured reservoir model using both the conventional transfer function, and a linear function that reproduces the results of imbibition experiments. The results will be compared to a grid-based dual porosity simulator with comparison of the oil recovery, run time and memory requirements as the number of grid blocks in the simulation is increased.
2. Streamline-based dual porosity simulator In a dual porosity simulator the flowing (fracture) and stagnant (matrix) domains are separately defined. It is assumed that there is no significant fluid flow between matrix grid blocks. The main flow in the reservoir occurs through exchange of fluid from the matrix to the fractures and in the fracture network itself. In the streamline-based dual porosity model, a high permeability matrix is combined with the fractures conceptually by considering its contribution to the effective permeability of the flowing fraction. The streamline code assumes incompressible two-phase flow of oil and water. While gravitational forces are included in the transport of fluid in the flowing regions, the only mechanism for fracture/matrix transfer allowed is capillary-controlled imbibition. The volume conservation equations are: /f
/m
BSwf þ vt jfwf þ j gG ¼ T Bt BSwm ¼T Bt
ð1Þ
ð2Þ
where T is the transfer function. The time-of-flight s(s) is defined as the time taken for neutral tracer to move a distance s along a streamline (Batycky et al., 1997; King and Datta-Gupta, 1998): Z s /f ds ð3Þ sðsÞ ¼ 0 Avt A Using Eq. (3) the following conservation equations are solved for both flowing and stagnant fractions along each streamline (ignoring gravity): BSwf Bfwf T þ ¼ /f Bt Bs
ð4Þ
BSwm T ¼ /m Bt
ð5Þ
The term involving gravity in Eq. (1) is handled using operator splitting. The details of how these equations are solved are presented by Di Donato et al. (2003). The key feature is that fluid is transported along periodically recomputed streamlines that define the flow field in the fractures. The fracture/matrix transfer T is accommodated through source/sink terms in the streamline transport equations.
3. Transfer functions and conventional simulation An industry-standard code, Eclipse (E100, 2001A Version) was selected as a comparative simulator. In Eclipse, the number of layers in the simulation is double the number of actual layers. The first half of the set is associated with the matrix blocks, and the second half with the fractures (Schlumberger, 2001). The porosity, permeability, depth, etc., are defined independently and a matrix-fracture coupling transmissibility is calculated automatically to simulate flow between the two systems. The transfer function assumes steady-state Darcy flow with upstream weighting (Gilman and Kazemi, 1983; Schechter, 1997): Tw uT ¼ rKm kwf ðPwf Pwm Þ
ð6Þ
To ¼ rKm kom ðPof Pom Þ
ð7Þ
where k is the mobility = kr/l and r, with the dimensions of 1/length2, is the shape factor and represents the
W. Huang et al. / Journal of Petroleum Science and Engineering 43 (2004) 129–137
inverse of the fracture spacing squared. r can be calculated by (Gilman and Kazemi, 1983): ! 1 1 1 r¼4 2þ 2þ 2 ; Lx Ly Lz
ð8Þ
where Lx, Ly and Lz are the dimensions of the matrix blocks in the x-, y- and z-directions, respectively, and represent the fracture spacing. Since the distribution of fractures is generally unknown, r is often treated as a history matching parameter. Because there is no viscous mediated flow in the stagnant regions, the matrix capillary pressure Pcm = Pom Pwm. For incompressible flow To + Tw = 0. Assuming that the capillary pressure in the matrix is much higher than in the fracture ( Pof = Pwf, PcfbPcm): T ¼ rKm
kwf kom Pcm kwf þ kom
ð9Þ
Note that the transfer function in Eq. (9) is a nonlinear function of both the fracture and matrix saturations through the mobilitites. This conventional transfer function, Eq. (9), was implemented in the streamline code. Another transfer function that varies linearly with the matrix saturation was also incorporated in the streamline code. Its functional form was chosen so that the simulator matched recoveries measured experimentally in water-wet blocks surrounded by water. This linear transfer function is given by (Di Donato et al., 2003): T ¼ b/m ð1 Sorm Swm Þ ¼0
Swf > 0 : Swf ¼ 0
ð10Þ
Values for the transfer rate b, and shape factor, r, were chosen so that Eqs. (9) and (10) gave the same value of T for a matrix water saturation of 0.5 and a fracture saturation of 1. The details of how this was done are given in Di Donato et al. (2003).
4. Reservoir model The reservoir model used was derived from the 10th SPE Comparative Solution Project, SPE-10 (Christie
131
and Blunt, 2001). It is a non-fractured sandstone, fluvial reservoir, based on a North Sea oilfield with high permeability meandering channels surrounded by low permeability shale. The permeability varies over more than four orders of magnitude. The dimensions of the reservoir in the x, y and z directions are 366, 670 and 52 m, respectively. The original data includes a description of the heterogeneous permeability and porosity distributions, compressible two-phase PVT data and relative permeabilities. To satisfy the need for testing the new streamline-based code, it was changed to a quasi-incompressible model by reducing the compressibility of the rock and fluid. This adjustment has a negligible impact on predictions of oil recovery, because the field is operated above the bubble point. In addition, for the dual porosity models: (1) a constant porosity for both matrix and fracture was assumed; (2) the original heterogeneous permeability data was treated as defining the fracture, with the exception that a minimum fracture permeability of 1 mD was set, while a single uniform permeability of 1 mD was assigned to the matrix; and (3) the original relative permeability curves were assigned to the matrix. For the fracture, a linear fractional flow for water was assumed, even if the oil and water viscosities were not equal. From this assumption, the following expressions were derived for the relative permeabilities in the fracture:
krwf
1 Swf lo ¼ 1þ Swf lw
krof
Swf lw ¼ 1þ 1 Swf lo
1 ;
ð11Þ
:
ð12Þ
1
The fracture capillary pressure was assumed to be zero. The matrix capillary pressure was: Pcm ¼
Pcmax
a Swm ð1 Sorm Þa ; a Swmi ð1 Sorm Þa
ð13Þ
where a = 1, S orm = 0.2, S wmi = 0.2 and P cmax = 7.4 105 Pa. To assess the influence of the number of grid blocks on the results of the simulations, in addition to the fine model (60, 220 and 85 grid blocks in the x, y and z directions, respectively) from SPE-10,
132
W. Huang et al. / Journal of Petroleum Science and Engineering 43 (2004) 129–137
Table 1 Parameters used in the simulations /m = 0.2 /f = 0.02 Km = 1.0 mD Kf = 1 – 20,000 mD Swim = Sorm = 0.2 Swif = Sorf = 0.0 krof(Swi) = krwf(Sor) = 1.0 krom(Swi) = krwm(Sor) = 1.0 lw = 0.3 10 3 Pa s lo = 3 10 3 Pa s qw = 1026 kg m 3 qo = 850 kg m 3 Bw = 1.01 cw = 4.35 10 10 Pa 1 cr = 1.45 10 10 Pa 1 b = 5 10 9 s 1 r (shape factor) = 0.05 m 2 (for dual porosity models) Q = 795 m3 day 1 (for single porosity models) Q = 300 m3 day 1 (for dual-porosity models) Bottom hole pressure = 2.758 107 Pa
another three coarser models, namely 60 55 85, 20 55 85, and 20 55 17, were generated by upscaling the original permeability data. Upscaling was performed using geometric averaging for Kx and Kz. Ky is equal to Kx. Geometric averaging was chosen since it was quickly and easily performed. This study was not designed to test different upscaling techniques. A finer, 4.5 million cell model (120 440 85) was created by downscaling the absolute permeability. The method used to do this is described in Appendix A. The same rock curves were used for all simulations (i.e. there is no upscaling of relative permeability). Each model was described on a regular Cartesian grid with a constant grid block size. The reference pressure was 4.13 105 Pa at a depth of 3657.6 m. For all models, one injector was located at the center with four producers at the corners. All wells were vertical and perforated throughout the formation. In addition, a constant wellbore rate was used for the injector and a constant wellbore pressure for the producers as the inner boundary conditions in addition to no flow at the outer boundary. The parameters used in the simulations are listed in Table 1.
The computed oil rate and average field pressure are essentially identical to the published results (Christie and Blunt, 2001). The same models were then run using a streamline-based single porosity code. As expected, for the same number of grid blocks, the results from the streamline code are almost identical to those from the grid-based method, Fig. 1. However, due to limitations on memory and computer time, the finest (1 million and 4.5 million grid block) single porosity models were not run using Eclipse. The oil rate decreases slightly as the grid is refined. This is because the finer grid can resolve small-scale channels of high permeability that result in earlier breakthrough. This effect of grid size has been previously discussed for this dataset by Christie and Blunt (2001).
6. Comparisons for dual porosity simulation Dual porosity simulation was then conducted using the two codes. For the streamline-based code, the two transfer functions were implemented in the conservation equation, as previously discussed. The oil and water rates for the two simulators with the same conventional transfer function and the same sized grid give almost identical results, Figs. 2 and 3. Again, due to memory and computer time limitations, the finest dual porosity Eclipse models were not run. As the grid is refined the oil rate decreases slightly since fine, high permeability channels are resolved. This leads to lower sweep efficiency with more oil being left behind in the matrix and lower recoveries.
5. Comparisons for single porosity simulation To validate the models, single porosity simulations were performed using Eclipse for different grid block sizes with the same parameters as the ones in SPE-10.
Fig. 1. Comparisons of oil rate for grid-based (Eclipse) and streamline-based (SL) single porosity models for different grid sizes.
W. Huang et al. / Journal of Petroleum Science and Engineering 43 (2004) 129–137
Fig. 2. Comparisons of oil rate for grid-based (Eclipse) and streamline-based (SL) dual porosity models using the conventional transfer function.
Fig. 4 shows comparisons of the oil rate using the streamline-based code with a linear transfer function. The reasons for using this transfer function are twofold: (1) for recovery controlled by counter-current imbibition, the transfer function reproduces corescale experimental results and hence is likely to be more physically valid than the conventional steadystate treatment of the problem; and (2) since the model is linear, it avoids an iterative approach to solve the transport equations along streamlines, leading to much faster run times, as discussed in the next section. At early times the linear model gives lower oil rates. This is because, at low matrix saturations corresponding to early time, the conventional model gives much higher transfer rates, thus boosting recovery. At late times, corresponding to high matrix saturations, the conventional transfer function gives lower transfer rates, meaning that recovery of matrix
Fig. 3. Comparisons of water production rate for grid-based (Eclipse) and streamline-based (SL) dual porosity models using the conventional transfer function.
133
Fig. 4. Comparisons of oil rate for streamline-based models with a linear transfer function (SL) and grid-based models with the conventional transfer function (Eclipse).
oil is slower and as a consequence lower oil rates are observed. This is discussed in more detail by Di Donato et al. (2003). As the models become finer, the disparities between the two methods become smaller. This point can be seen by comparing the predictions of oil recovery, Fig. 5. For the pair of 20 55 17 models, the disparity of two methods at 9000 days is a recovery factor of 0.05; for the pair of 60 55 85 models, the disparity is just 0.03. Therefore it is reasonable to expect that the results from a million grid block model will be closer to each other. This is because the finer models resolve the influence of heterogeneity more accurately, and in comparison, the effect of the matrix/fracture transfer rate on recovery becomes less significant. Fig. 6 shows saturation profiles after 2000 days for the 60th layer of the 60 55 85 model from the streamline code (middle and bottom) and the gridbased code (top). The areas with lighter color, corresponding to low saturation, are where the fracture permeability is small; the areas with darker color are regions with relatively high permeability fracture through which the injected water has channeled. The streamline code gives very similar predictions for the distribution of the fluid in both matrix (left) and fracture (right) as the grid-based code, since the flow is dominated by the extreme permeability contrast in the fractures. However, there has been less transfer into the matrix in the streamline code using the linear transfer function. This is consistent with the lower recoveries at 2000 days shown in Fig. 5. The linear model predicts less transfer into the matrix at early
134
W. Huang et al. / Journal of Petroleum Science and Engineering 43 (2004) 129–137
7. Comparisons of CPU time and memory requirements
Fig. 5. Comparisons of recovery for grid-based (Eclipse) and linear streamline-based dual porosity (SL) models for different grid sizes. The results are sensitive to the grid size. The larger grid size, the smaller the difference between the two methods.
and intermediate times than an equivalent conventional model (Di Donato et al., 2003). When the conventional transfer function is used in the streamline code (bottom figure), the saturation profiles are virtually identical to those obtained using the gridbased approach.
Both single porosity and dual porosity models were run with different grid sizes for a simulation time of 2000 days. All models were run on a 2 GB RAM, 2.4 GHz dual-processor Pentium Xeon workstation operated under Windows XP. Fig. 7 gives the CPU time as a function of the number of grid blocks for each model. Table 2 gives the run times, time steps and the memory requirements. Table 2 also shows the speed-up factors for the streamline-based code, which is defined as the ratio of the run time for a grid-based model to the run time for a corresponding streamline-based model. The grid-based simulations were all run in fully implicit mode and the time step size was chosen to give optimal run times. Smaller time steps lead to longer run times without affecting the results, whereas longer time steps lead to convergence problems and instabilities in the solutions. Fig. 8 plots the speed-up factors as a function of the number of grid blocks. Fig. 9 shows
Sw
Fig. 6. Saturation profiles of the 60th layer in matrix (left) and fracture (right) after 2000 days for dual porosity simulation of the 605585 model for the grid-based code (top), the linear streamline code (middle) and the conventional streamline code (bottom).
W. Huang et al. / Journal of Petroleum Science and Engineering 43 (2004) 129–137
Fig. 7. CPU times as a function of the number of grid blocks for the streamline-based code and the grid-based code (Eclipse) for both single porosity and dual porosity models.
the memory requirements as a function of the number of grid blocks. The following observations are evident from Figs. 7 –9 and Table 2: (1) For the streamline-based code, in both single porosity and dual porosity mode, the run time scales approximately linearly with the number of grid blocks. In contrast, the grid-based code gives run times that are approximately proportional to the number of grid blocks squared. (2) For both simulators, dual porosity models are slower than corresponding single porosity models. This is because with the conventional transfer
135
Fig. 8. Speed-up factors for the streamline-based code for different grid sizes for both single porosity and dual porosity models.
function, a more complex one-dimensional transport solver is required that requires an iterative solution. The linear transfer function gives faster run times than the conventional transfer function since this iterative solution is avoided. (3) For the streamline-based code, fewer than 20 time steps were used to simulate 2000 days of production. Using more time steps did not affect the results. In the grid-based code, however, the number of time steps increased rapidly with the increase in the number of grid blocks, particularly for the dual porosity models. Convergence and stability constraints mean that the grid-based simulator requires many more time steps than
Table 2 Run times, memory requirements and speed-up factors for each model Number of grid blocks
18,700 Single porosity 18,700 Dual porosity 93,500 Single porosity 93,500 Dual porosity 280,500 Single porosity 280,500 Dual porosity 1,122,000 Single porosity 1,122,000 Dual porosity 4,488,000 Single porosity 4,488,000 Dual porosity
Streamline-based code
Linear Conventional Linear Conventional Linear Conventional Linear Conventional Linear Conventional
Time (min)
RAM (MB)
0.2 0.3 1.5 2.4 8 15 9 30 58 64 188 383 416 1389 3825
109 109 147 147 202 202 507 507 1799 1799
Grid-based code Time steps 12 13 15 12 14 16 12 14 17 16 14 15 13 15 17
Time (min) 1.7 12
RAM (MB) 52 72
Time steps 79 851
253 526
522 582
417 1,073
1739 4268
1110 1312
1029 10,009
/
>2GB, 64bit
/
/
>2GB, 64bit
/
Speed-up factor
7 36 8 109 65 36 185 143 74 / / / / / /
136
W. Huang et al. / Journal of Petroleum Science and Engineering 43 (2004) 129–137
Fig. 9. Memory requirements for the streamline-based and gridbased codes for both single porosity and dual porosity models.
the streamline model that has no restriction on the size of time step. This explains why the grid-based code takes more time. (4) For dual porosity models, the speed-up factors for the linear transfer function are larger than for the conventional transfer function. As mentioned above, this is attributed to simplified saturation updating along the streamlines for the linear model. (5) As the number of grid blocks increases, the speedup factor increases for all models (Fig. 9). This is consistent with the scaling of run time with the number of grid blocks mentioned in point 1. (6) The memory requirements for streamline-based models are much less than for grid-based models. In the grid-based method, since the simulator generates one set of grid blocks for matrix parameters and the other for the fracture parameters, the number of grid blocks in a dual porosity model is double that of single porosity models. Hence, the memory requirements are larger in dual porosity models. In contrast, in the streamlinebased method, the memory requirements are the same for both single and dual porosity models since the same number of grid blocks is used.
appropriate—incompressible displacement in a highly heterogeneous reservoir. For the situations studied, streamline-based simulation gave almost identical results to commercial grid-based simulation for single and dual porosity models. In the streamline code a linear transfer function was introduced that gave much faster run times and reproduces the results of core-scale imbibition experiments. With an appropriate choice of the transfer rate, the results using the linear transfer function were similar to those obtained using the conventional steady-state transfer function. The grid-based code is orders of magnitude slower than streamline-based simulation for both single and dual porosity models. The linear streamline-based models are particularly fast since they avoid an iterative solution in the transport equations along each streamline. The larger the number of grid blocks, the larger is the speed-up. Therefore, for the model with one million grid blocks it is expected that the streamline approach will be at least one hundred times faster than using a grid-based code. It was possible to run a 4.5 million grid block dual porosity model using a standard PC. The formulation presented assumes incompressible flow. For compressible systems, such as when gas is present, or where there is significant pressure depletion, the performance of the streamline code is unlikely to be so favorable.
8. Conclusions
Nomenclature f fractional flow, dimensionless g acceleration due to gravity, Lt 2, ms 2 G gravity fractional flow, dimensionless K permeability, L2, m2 or D kr relative permeability Pc capillary pressure, ml 1 t 2, Pa Q injection rate, L3 t 1, m3 s 1 R recovery, dimensionless S saturation, dimensionless t time, T, s T transfer function, T 1, s 1 vt total velocity, LT 1, ms 1
The results of streamline-based simulations with different grid sizes for both single porosity and dual porosity models were compared with equivalent predictions from a grid-based simulator. The cases were those where streamline-based simulation is particularly
Greek Symbols b rate constant in linear transfer function, T 1, s 1 / porosity, dimensionless k mobility, M 1Lt, Pa 1.s 1
W. Huang et al. / Journal of Petroleum Science and Engineering 43 (2004) 129–137
l q r s
viscosity, ML 1t 1, Pa.s density, ML 3, kg.m 3 shape factor, L 2, m 2 time of flight, T, s
for each new part. a is a parameter to be determined. Here it is set to 3.5. For grid blocks with a permeability K = 1 mD (the minimum allowed fracture permeability), Ki = 1 mD. This method ensures that the geometric average permeability is the same as the original field. The parameter a is a measure of the small-scale randomness in the permeability.
Subscripts f flowing or fracture i initial m matrix or stagnant o oil r residual t total w water
References
Acknowledgements We would like to thank Shell and the sponsors of the ITF project on ‘Improved Simulation of Flow in Fractured and Faulted Reservoirs’ for supporting this research. Wenfen Huang thanks the China Scholarship Council for financial support. We are also grateful to E I Obi for his help in running the streamline simulations.
Appendix A . Downscaling permeability data To generate a finer fracture permeability dataset (with 4,488,000 grid blocks) from SPE-10, a downscaling technique that preserves the geometric average of the permeability was employed. Each block was divided into four parts, by dividing the block in half in the x and y directions. Four random numbers (x1, x2, x3, x4) were generated uniformly in the range 0 – 1. The permeabilities of the downscaled grid blocks were then computed as follows: Ki ¼ K eaðxi ¯xÞ ; x¯ ¼
i ¼ 1; 2; 3; 4
1 ðx1 þ x2 þ x3 þ x4 Þ: 4
137
ð14Þ ð15Þ
where K is the original permeability of the block before downscaling and Ki is the new permeability
Baker, R.O., Kuppe, F., Chugh, S., Bora, R., Stojanovic, S., Batycky, R., 2002. Full-field modeling using streamline-based simulation: four case studies. SPE REE 5 (2), 126 – 134. Barenblatt, G.I., Zheltov, I.P., Kochina, I.N., 1960. Basic concepts in the theory of seepage of homogeneous liquids in fissured rocks. J. Appl. Math. Mech., Engl. Transl. 24, 1286 – 1303. Batycky, R.P., Blunt, M.J., Thiele, M.R., 1997. 3D field-scale streamline-based reservoir simulator. SPE RE 11, 246 – 254. Christie, M.A., Blunt, M.J., 2001. Tenth SPE comparative solution project: a comparison of upscaling techniques. SPE REE 4 (4), 308 – 317 (http://www.spe.org/csp). Di Donato, G., Huang, W., Blunt, M.J., 2003. Streamline-based dual porosity simulation of fractured reservoirs. SPE paper 84036. Proceedings of the SPE Annual Technical Meeting and Exhibition, Denver, Colorado, October. Gilman, J.R., Kazemi, H., 1983. Improvement in simulation of naturally fractured reservoirs. SPE J. 23, 695 – 707. Grinestaff, G.H., Caffrey, D.J., 2000. Waterflood management: a case study of the Northwest fault block area of Prudhoe Bay, Alaska, using streamline simulation and traditional waterflood analysis. SPE paper 63152. Proceedings of the SPE Annual Technical Conference and Exhibition, Dallas, 1 – 4 October. King, M.J., Datta-Gupta, A., 1998. Streamline simulation: a current perspective. In Situ 22 (1), 91 – 140. Nakashima, T., Arihara, N., Sutopo, N., Sato, K., 2001. Effective permeability estimation for modeling naturally fractured reservoirs. SPE paper 68124. Proceedings of the SPE Middle East Oil Show held in Bahrain, 17 – 20 March. Samier, P., Quettier, L., Thiele, M., 2002. Applications of streamline simulations to reservoir studies. SPE REE 5 (4), 324 – 332. Schechter D.S., 1997. Numerical simulation of a water flood pilot in the naturally fractured Sprabeery trend. 2nd Annual Technical Report to DOE Class III Field Demonstration Project in the Naturally Fractured Spraberry Trend. Schlumberger, 2001. Eclipse technical description. Manuals 2001A, 119 – 136. Warren, J.E., Root, P.J., 1963. The behavior of naturally fractured reservoirs. SPE J. 3, 245 – 255; Trans. AIME 228, 245 – 255.
Comparison of streamline-based and grid-based dual porosity simulation Wenfen Huang a,b, Ginevra Di Donato a, Martin J. Blunt a,* a
Department of Earth Science and Engineering, Imperial College, London SW7 2AZ, UK b Shengli Geological Scientific Research Institute, Sinopec Inc., China Received 21 May 2003; accepted 22 January 2004
Abstract In field-scale dual porosity models, fractured reservoirs are represented as two intercommunicating domains: a flowing fraction composed of the connected fracture network and high permeability matrix; and a stagnant fraction composed of lower permeability matrix containing most of the original oil. The transfer of fluids between fracture and matrix is modeled by an empirical function of saturation. Two transfer functions are used in a streamline-based dual porosity code: (1) a conventional transfer function used in traditional grid-based simulators; and (2) a new linear transfer function that reproduces the results of core-scale imbibition experiments. The results of the simulations are compared with those from a conventional grid-based finite difference code. The two codes give very similar results, but the streamline code requires a much shorter run time and has a lower memory requirement. Using the linear transfer function, the streamline-based code is orders of magnitude faster and requires many times less memory than the conventional grid-based code. The larger the number of grid blocks, the greater the relative speed of the streamline-based approach. It is possible to run a 4.5 million grid block dual porosity model on a standard PC using a streamline code. D 2004 Elsevier B.V. All rights reserved. Keywords: Streamlines; Dual porosity; Fractured reservoir; Transfer function
1. Introduction Traditionally, to model a fractured reservoir at the field scale, a flowing region with high permeability represents the fracture network, while a stagnant region with low permeability represents the matrix (Barenblatt et al., 1960; Warren and Root, 1963). A dual porosity (no flow between matrix blocks, as implemented in this paper) or dual permeability (there
* Corresponding author. Tel.: +44-20-7594-6500; fax: +44-207594-7444. E-mail address: [email protected] (M.J. Blunt). 0920-4105/$ - see front matter D 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.petrol.2004.01.002
is flow between matrix blocks) approach is applicable where the majority of the flow is through a wellconnected fracture network (Gilman and Kazemi, 1983; Nakashima et al., 2001). The fluid flow between fracture and matrix is modeled using empirical transfer functions. In most studies a conventional steady-state transfer function is used that varies non-linearly with fracture and matrix saturations. The nonlinearity of the transfer function makes an iterative approach to solving the transport equations necessary, resulting in longer run times than equivalent single porosity (non-fractured) simulations. Streamline-based simulation is frequently used for modeling recovery by displacement because it is
130
W. Huang et al. / Journal of Petroleum Science and Engineering 43 (2004) 129–137
generally faster, suffers from less numerical dispersion and can handle larger models than grid-based codes (Batycky et al., 1997). Streamline methods are ideal for simulating incompressible, or nearly incompressible, systems in highly heterogeneous reservoirs. When the system becomes more compressible and the voidage replacement ratio is much less than one, the method becomes less efficient (Baker et al., 2002). Streamline methods have been successfully applied to several field studies (Samier et al., 2002; Baker et al., 2002; Grinestaff and Caffrey, 2000). Recently, a dual porosity streamline-based code to simulate waterflooding was published (Di Donato et al., 2003). This code will be used to simulate flow in a synthetic fractured reservoir model using both the conventional transfer function, and a linear function that reproduces the results of imbibition experiments. The results will be compared to a grid-based dual porosity simulator with comparison of the oil recovery, run time and memory requirements as the number of grid blocks in the simulation is increased.
2. Streamline-based dual porosity simulator In a dual porosity simulator the flowing (fracture) and stagnant (matrix) domains are separately defined. It is assumed that there is no significant fluid flow between matrix grid blocks. The main flow in the reservoir occurs through exchange of fluid from the matrix to the fractures and in the fracture network itself. In the streamline-based dual porosity model, a high permeability matrix is combined with the fractures conceptually by considering its contribution to the effective permeability of the flowing fraction. The streamline code assumes incompressible two-phase flow of oil and water. While gravitational forces are included in the transport of fluid in the flowing regions, the only mechanism for fracture/matrix transfer allowed is capillary-controlled imbibition. The volume conservation equations are: /f
/m
BSwf þ vt jfwf þ j gG ¼ T Bt BSwm ¼T Bt
ð1Þ
ð2Þ
where T is the transfer function. The time-of-flight s(s) is defined as the time taken for neutral tracer to move a distance s along a streamline (Batycky et al., 1997; King and Datta-Gupta, 1998): Z s /f ds ð3Þ sðsÞ ¼ 0 Avt A Using Eq. (3) the following conservation equations are solved for both flowing and stagnant fractions along each streamline (ignoring gravity): BSwf Bfwf T þ ¼ /f Bt Bs
ð4Þ
BSwm T ¼ /m Bt
ð5Þ
The term involving gravity in Eq. (1) is handled using operator splitting. The details of how these equations are solved are presented by Di Donato et al. (2003). The key feature is that fluid is transported along periodically recomputed streamlines that define the flow field in the fractures. The fracture/matrix transfer T is accommodated through source/sink terms in the streamline transport equations.
3. Transfer functions and conventional simulation An industry-standard code, Eclipse (E100, 2001A Version) was selected as a comparative simulator. In Eclipse, the number of layers in the simulation is double the number of actual layers. The first half of the set is associated with the matrix blocks, and the second half with the fractures (Schlumberger, 2001). The porosity, permeability, depth, etc., are defined independently and a matrix-fracture coupling transmissibility is calculated automatically to simulate flow between the two systems. The transfer function assumes steady-state Darcy flow with upstream weighting (Gilman and Kazemi, 1983; Schechter, 1997): Tw uT ¼ rKm kwf ðPwf Pwm Þ
ð6Þ
To ¼ rKm kom ðPof Pom Þ
ð7Þ
where k is the mobility = kr/l and r, with the dimensions of 1/length2, is the shape factor and represents the
W. Huang et al. / Journal of Petroleum Science and Engineering 43 (2004) 129–137
inverse of the fracture spacing squared. r can be calculated by (Gilman and Kazemi, 1983): ! 1 1 1 r¼4 2þ 2þ 2 ; Lx Ly Lz
ð8Þ
where Lx, Ly and Lz are the dimensions of the matrix blocks in the x-, y- and z-directions, respectively, and represent the fracture spacing. Since the distribution of fractures is generally unknown, r is often treated as a history matching parameter. Because there is no viscous mediated flow in the stagnant regions, the matrix capillary pressure Pcm = Pom Pwm. For incompressible flow To + Tw = 0. Assuming that the capillary pressure in the matrix is much higher than in the fracture ( Pof = Pwf, PcfbPcm): T ¼ rKm
kwf kom Pcm kwf þ kom
ð9Þ
Note that the transfer function in Eq. (9) is a nonlinear function of both the fracture and matrix saturations through the mobilitites. This conventional transfer function, Eq. (9), was implemented in the streamline code. Another transfer function that varies linearly with the matrix saturation was also incorporated in the streamline code. Its functional form was chosen so that the simulator matched recoveries measured experimentally in water-wet blocks surrounded by water. This linear transfer function is given by (Di Donato et al., 2003): T ¼ b/m ð1 Sorm Swm Þ ¼0
Swf > 0 : Swf ¼ 0
ð10Þ
Values for the transfer rate b, and shape factor, r, were chosen so that Eqs. (9) and (10) gave the same value of T for a matrix water saturation of 0.5 and a fracture saturation of 1. The details of how this was done are given in Di Donato et al. (2003).
4. Reservoir model The reservoir model used was derived from the 10th SPE Comparative Solution Project, SPE-10 (Christie
131
and Blunt, 2001). It is a non-fractured sandstone, fluvial reservoir, based on a North Sea oilfield with high permeability meandering channels surrounded by low permeability shale. The permeability varies over more than four orders of magnitude. The dimensions of the reservoir in the x, y and z directions are 366, 670 and 52 m, respectively. The original data includes a description of the heterogeneous permeability and porosity distributions, compressible two-phase PVT data and relative permeabilities. To satisfy the need for testing the new streamline-based code, it was changed to a quasi-incompressible model by reducing the compressibility of the rock and fluid. This adjustment has a negligible impact on predictions of oil recovery, because the field is operated above the bubble point. In addition, for the dual porosity models: (1) a constant porosity for both matrix and fracture was assumed; (2) the original heterogeneous permeability data was treated as defining the fracture, with the exception that a minimum fracture permeability of 1 mD was set, while a single uniform permeability of 1 mD was assigned to the matrix; and (3) the original relative permeability curves were assigned to the matrix. For the fracture, a linear fractional flow for water was assumed, even if the oil and water viscosities were not equal. From this assumption, the following expressions were derived for the relative permeabilities in the fracture:
krwf
1 Swf lo ¼ 1þ Swf lw
krof
Swf lw ¼ 1þ 1 Swf lo
1 ;
ð11Þ
:
ð12Þ
1
The fracture capillary pressure was assumed to be zero. The matrix capillary pressure was: Pcm ¼
Pcmax
a Swm ð1 Sorm Þa ; a Swmi ð1 Sorm Þa
ð13Þ
where a = 1, S orm = 0.2, S wmi = 0.2 and P cmax = 7.4 105 Pa. To assess the influence of the number of grid blocks on the results of the simulations, in addition to the fine model (60, 220 and 85 grid blocks in the x, y and z directions, respectively) from SPE-10,
132
W. Huang et al. / Journal of Petroleum Science and Engineering 43 (2004) 129–137
Table 1 Parameters used in the simulations /m = 0.2 /f = 0.02 Km = 1.0 mD Kf = 1 – 20,000 mD Swim = Sorm = 0.2 Swif = Sorf = 0.0 krof(Swi) = krwf(Sor) = 1.0 krom(Swi) = krwm(Sor) = 1.0 lw = 0.3 10 3 Pa s lo = 3 10 3 Pa s qw = 1026 kg m 3 qo = 850 kg m 3 Bw = 1.01 cw = 4.35 10 10 Pa 1 cr = 1.45 10 10 Pa 1 b = 5 10 9 s 1 r (shape factor) = 0.05 m 2 (for dual porosity models) Q = 795 m3 day 1 (for single porosity models) Q = 300 m3 day 1 (for dual-porosity models) Bottom hole pressure = 2.758 107 Pa
another three coarser models, namely 60 55 85, 20 55 85, and 20 55 17, were generated by upscaling the original permeability data. Upscaling was performed using geometric averaging for Kx and Kz. Ky is equal to Kx. Geometric averaging was chosen since it was quickly and easily performed. This study was not designed to test different upscaling techniques. A finer, 4.5 million cell model (120 440 85) was created by downscaling the absolute permeability. The method used to do this is described in Appendix A. The same rock curves were used for all simulations (i.e. there is no upscaling of relative permeability). Each model was described on a regular Cartesian grid with a constant grid block size. The reference pressure was 4.13 105 Pa at a depth of 3657.6 m. For all models, one injector was located at the center with four producers at the corners. All wells were vertical and perforated throughout the formation. In addition, a constant wellbore rate was used for the injector and a constant wellbore pressure for the producers as the inner boundary conditions in addition to no flow at the outer boundary. The parameters used in the simulations are listed in Table 1.
The computed oil rate and average field pressure are essentially identical to the published results (Christie and Blunt, 2001). The same models were then run using a streamline-based single porosity code. As expected, for the same number of grid blocks, the results from the streamline code are almost identical to those from the grid-based method, Fig. 1. However, due to limitations on memory and computer time, the finest (1 million and 4.5 million grid block) single porosity models were not run using Eclipse. The oil rate decreases slightly as the grid is refined. This is because the finer grid can resolve small-scale channels of high permeability that result in earlier breakthrough. This effect of grid size has been previously discussed for this dataset by Christie and Blunt (2001).
6. Comparisons for dual porosity simulation Dual porosity simulation was then conducted using the two codes. For the streamline-based code, the two transfer functions were implemented in the conservation equation, as previously discussed. The oil and water rates for the two simulators with the same conventional transfer function and the same sized grid give almost identical results, Figs. 2 and 3. Again, due to memory and computer time limitations, the finest dual porosity Eclipse models were not run. As the grid is refined the oil rate decreases slightly since fine, high permeability channels are resolved. This leads to lower sweep efficiency with more oil being left behind in the matrix and lower recoveries.
5. Comparisons for single porosity simulation To validate the models, single porosity simulations were performed using Eclipse for different grid block sizes with the same parameters as the ones in SPE-10.
Fig. 1. Comparisons of oil rate for grid-based (Eclipse) and streamline-based (SL) single porosity models for different grid sizes.
W. Huang et al. / Journal of Petroleum Science and Engineering 43 (2004) 129–137
Fig. 2. Comparisons of oil rate for grid-based (Eclipse) and streamline-based (SL) dual porosity models using the conventional transfer function.
Fig. 4 shows comparisons of the oil rate using the streamline-based code with a linear transfer function. The reasons for using this transfer function are twofold: (1) for recovery controlled by counter-current imbibition, the transfer function reproduces corescale experimental results and hence is likely to be more physically valid than the conventional steadystate treatment of the problem; and (2) since the model is linear, it avoids an iterative approach to solve the transport equations along streamlines, leading to much faster run times, as discussed in the next section. At early times the linear model gives lower oil rates. This is because, at low matrix saturations corresponding to early time, the conventional model gives much higher transfer rates, thus boosting recovery. At late times, corresponding to high matrix saturations, the conventional transfer function gives lower transfer rates, meaning that recovery of matrix
Fig. 3. Comparisons of water production rate for grid-based (Eclipse) and streamline-based (SL) dual porosity models using the conventional transfer function.
133
Fig. 4. Comparisons of oil rate for streamline-based models with a linear transfer function (SL) and grid-based models with the conventional transfer function (Eclipse).
oil is slower and as a consequence lower oil rates are observed. This is discussed in more detail by Di Donato et al. (2003). As the models become finer, the disparities between the two methods become smaller. This point can be seen by comparing the predictions of oil recovery, Fig. 5. For the pair of 20 55 17 models, the disparity of two methods at 9000 days is a recovery factor of 0.05; for the pair of 60 55 85 models, the disparity is just 0.03. Therefore it is reasonable to expect that the results from a million grid block model will be closer to each other. This is because the finer models resolve the influence of heterogeneity more accurately, and in comparison, the effect of the matrix/fracture transfer rate on recovery becomes less significant. Fig. 6 shows saturation profiles after 2000 days for the 60th layer of the 60 55 85 model from the streamline code (middle and bottom) and the gridbased code (top). The areas with lighter color, corresponding to low saturation, are where the fracture permeability is small; the areas with darker color are regions with relatively high permeability fracture through which the injected water has channeled. The streamline code gives very similar predictions for the distribution of the fluid in both matrix (left) and fracture (right) as the grid-based code, since the flow is dominated by the extreme permeability contrast in the fractures. However, there has been less transfer into the matrix in the streamline code using the linear transfer function. This is consistent with the lower recoveries at 2000 days shown in Fig. 5. The linear model predicts less transfer into the matrix at early
134
W. Huang et al. / Journal of Petroleum Science and Engineering 43 (2004) 129–137
7. Comparisons of CPU time and memory requirements
Fig. 5. Comparisons of recovery for grid-based (Eclipse) and linear streamline-based dual porosity (SL) models for different grid sizes. The results are sensitive to the grid size. The larger grid size, the smaller the difference between the two methods.
and intermediate times than an equivalent conventional model (Di Donato et al., 2003). When the conventional transfer function is used in the streamline code (bottom figure), the saturation profiles are virtually identical to those obtained using the gridbased approach.
Both single porosity and dual porosity models were run with different grid sizes for a simulation time of 2000 days. All models were run on a 2 GB RAM, 2.4 GHz dual-processor Pentium Xeon workstation operated under Windows XP. Fig. 7 gives the CPU time as a function of the number of grid blocks for each model. Table 2 gives the run times, time steps and the memory requirements. Table 2 also shows the speed-up factors for the streamline-based code, which is defined as the ratio of the run time for a grid-based model to the run time for a corresponding streamline-based model. The grid-based simulations were all run in fully implicit mode and the time step size was chosen to give optimal run times. Smaller time steps lead to longer run times without affecting the results, whereas longer time steps lead to convergence problems and instabilities in the solutions. Fig. 8 plots the speed-up factors as a function of the number of grid blocks. Fig. 9 shows
Sw
Fig. 6. Saturation profiles of the 60th layer in matrix (left) and fracture (right) after 2000 days for dual porosity simulation of the 605585 model for the grid-based code (top), the linear streamline code (middle) and the conventional streamline code (bottom).
W. Huang et al. / Journal of Petroleum Science and Engineering 43 (2004) 129–137
Fig. 7. CPU times as a function of the number of grid blocks for the streamline-based code and the grid-based code (Eclipse) for both single porosity and dual porosity models.
the memory requirements as a function of the number of grid blocks. The following observations are evident from Figs. 7 –9 and Table 2: (1) For the streamline-based code, in both single porosity and dual porosity mode, the run time scales approximately linearly with the number of grid blocks. In contrast, the grid-based code gives run times that are approximately proportional to the number of grid blocks squared. (2) For both simulators, dual porosity models are slower than corresponding single porosity models. This is because with the conventional transfer
135
Fig. 8. Speed-up factors for the streamline-based code for different grid sizes for both single porosity and dual porosity models.
function, a more complex one-dimensional transport solver is required that requires an iterative solution. The linear transfer function gives faster run times than the conventional transfer function since this iterative solution is avoided. (3) For the streamline-based code, fewer than 20 time steps were used to simulate 2000 days of production. Using more time steps did not affect the results. In the grid-based code, however, the number of time steps increased rapidly with the increase in the number of grid blocks, particularly for the dual porosity models. Convergence and stability constraints mean that the grid-based simulator requires many more time steps than
Table 2 Run times, memory requirements and speed-up factors for each model Number of grid blocks
18,700 Single porosity 18,700 Dual porosity 93,500 Single porosity 93,500 Dual porosity 280,500 Single porosity 280,500 Dual porosity 1,122,000 Single porosity 1,122,000 Dual porosity 4,488,000 Single porosity 4,488,000 Dual porosity
Streamline-based code
Linear Conventional Linear Conventional Linear Conventional Linear Conventional Linear Conventional
Time (min)
RAM (MB)
0.2 0.3 1.5 2.4 8 15 9 30 58 64 188 383 416 1389 3825
109 109 147 147 202 202 507 507 1799 1799
Grid-based code Time steps 12 13 15 12 14 16 12 14 17 16 14 15 13 15 17
Time (min) 1.7 12
RAM (MB) 52 72
Time steps 79 851
253 526
522 582
417 1,073
1739 4268
1110 1312
1029 10,009
/
>2GB, 64bit
/
/
>2GB, 64bit
/
Speed-up factor
7 36 8 109 65 36 185 143 74 / / / / / /
136
W. Huang et al. / Journal of Petroleum Science and Engineering 43 (2004) 129–137
Fig. 9. Memory requirements for the streamline-based and gridbased codes for both single porosity and dual porosity models.
the streamline model that has no restriction on the size of time step. This explains why the grid-based code takes more time. (4) For dual porosity models, the speed-up factors for the linear transfer function are larger than for the conventional transfer function. As mentioned above, this is attributed to simplified saturation updating along the streamlines for the linear model. (5) As the number of grid blocks increases, the speedup factor increases for all models (Fig. 9). This is consistent with the scaling of run time with the number of grid blocks mentioned in point 1. (6) The memory requirements for streamline-based models are much less than for grid-based models. In the grid-based method, since the simulator generates one set of grid blocks for matrix parameters and the other for the fracture parameters, the number of grid blocks in a dual porosity model is double that of single porosity models. Hence, the memory requirements are larger in dual porosity models. In contrast, in the streamlinebased method, the memory requirements are the same for both single and dual porosity models since the same number of grid blocks is used.
appropriate—incompressible displacement in a highly heterogeneous reservoir. For the situations studied, streamline-based simulation gave almost identical results to commercial grid-based simulation for single and dual porosity models. In the streamline code a linear transfer function was introduced that gave much faster run times and reproduces the results of core-scale imbibition experiments. With an appropriate choice of the transfer rate, the results using the linear transfer function were similar to those obtained using the conventional steady-state transfer function. The grid-based code is orders of magnitude slower than streamline-based simulation for both single and dual porosity models. The linear streamline-based models are particularly fast since they avoid an iterative solution in the transport equations along each streamline. The larger the number of grid blocks, the larger is the speed-up. Therefore, for the model with one million grid blocks it is expected that the streamline approach will be at least one hundred times faster than using a grid-based code. It was possible to run a 4.5 million grid block dual porosity model using a standard PC. The formulation presented assumes incompressible flow. For compressible systems, such as when gas is present, or where there is significant pressure depletion, the performance of the streamline code is unlikely to be so favorable.
8. Conclusions
Nomenclature f fractional flow, dimensionless g acceleration due to gravity, Lt 2, ms 2 G gravity fractional flow, dimensionless K permeability, L2, m2 or D kr relative permeability Pc capillary pressure, ml 1 t 2, Pa Q injection rate, L3 t 1, m3 s 1 R recovery, dimensionless S saturation, dimensionless t time, T, s T transfer function, T 1, s 1 vt total velocity, LT 1, ms 1
The results of streamline-based simulations with different grid sizes for both single porosity and dual porosity models were compared with equivalent predictions from a grid-based simulator. The cases were those where streamline-based simulation is particularly
Greek Symbols b rate constant in linear transfer function, T 1, s 1 / porosity, dimensionless k mobility, M 1Lt, Pa 1.s 1
W. Huang et al. / Journal of Petroleum Science and Engineering 43 (2004) 129–137
l q r s
viscosity, ML 1t 1, Pa.s density, ML 3, kg.m 3 shape factor, L 2, m 2 time of flight, T, s
for each new part. a is a parameter to be determined. Here it is set to 3.5. For grid blocks with a permeability K = 1 mD (the minimum allowed fracture permeability), Ki = 1 mD. This method ensures that the geometric average permeability is the same as the original field. The parameter a is a measure of the small-scale randomness in the permeability.
Subscripts f flowing or fracture i initial m matrix or stagnant o oil r residual t total w water
References
Acknowledgements We would like to thank Shell and the sponsors of the ITF project on ‘Improved Simulation of Flow in Fractured and Faulted Reservoirs’ for supporting this research. Wenfen Huang thanks the China Scholarship Council for financial support. We are also grateful to E I Obi for his help in running the streamline simulations.
Appendix A . Downscaling permeability data To generate a finer fracture permeability dataset (with 4,488,000 grid blocks) from SPE-10, a downscaling technique that preserves the geometric average of the permeability was employed. Each block was divided into four parts, by dividing the block in half in the x and y directions. Four random numbers (x1, x2, x3, x4) were generated uniformly in the range 0 – 1. The permeabilities of the downscaled grid blocks were then computed as follows: Ki ¼ K eaðxi ¯xÞ ; x¯ ¼
i ¼ 1; 2; 3; 4
1 ðx1 þ x2 þ x3 þ x4 Þ: 4
137
ð14Þ ð15Þ
where K is the original permeability of the block before downscaling and Ki is the new permeability
Baker, R.O., Kuppe, F., Chugh, S., Bora, R., Stojanovic, S., Batycky, R., 2002. Full-field modeling using streamline-based simulation: four case studies. SPE REE 5 (2), 126 – 134. Barenblatt, G.I., Zheltov, I.P., Kochina, I.N., 1960. Basic concepts in the theory of seepage of homogeneous liquids in fissured rocks. J. Appl. Math. Mech., Engl. Transl. 24, 1286 – 1303. Batycky, R.P., Blunt, M.J., Thiele, M.R., 1997. 3D field-scale streamline-based reservoir simulator. SPE RE 11, 246 – 254. Christie, M.A., Blunt, M.J., 2001. Tenth SPE comparative solution project: a comparison of upscaling techniques. SPE REE 4 (4), 308 – 317 (http://www.spe.org/csp). Di Donato, G., Huang, W., Blunt, M.J., 2003. Streamline-based dual porosity simulation of fractured reservoirs. SPE paper 84036. Proceedings of the SPE Annual Technical Meeting and Exhibition, Denver, Colorado, October. Gilman, J.R., Kazemi, H., 1983. Improvement in simulation of naturally fractured reservoirs. SPE J. 23, 695 – 707. Grinestaff, G.H., Caffrey, D.J., 2000. Waterflood management: a case study of the Northwest fault block area of Prudhoe Bay, Alaska, using streamline simulation and traditional waterflood analysis. SPE paper 63152. Proceedings of the SPE Annual Technical Conference and Exhibition, Dallas, 1 – 4 October. King, M.J., Datta-Gupta, A., 1998. Streamline simulation: a current perspective. In Situ 22 (1), 91 – 140. Nakashima, T., Arihara, N., Sutopo, N., Sato, K., 2001. Effective permeability estimation for modeling naturally fractured reservoirs. SPE paper 68124. Proceedings of the SPE Middle East Oil Show held in Bahrain, 17 – 20 March. Samier, P., Quettier, L., Thiele, M., 2002. Applications of streamline simulations to reservoir studies. SPE REE 5 (4), 324 – 332. Schechter D.S., 1997. Numerical simulation of a water flood pilot in the naturally fractured Sprabeery trend. 2nd Annual Technical Report to DOE Class III Field Demonstration Project in the Naturally Fractured Spraberry Trend. Schlumberger, 2001. Eclipse technical description. Manuals 2001A, 119 – 136. Warren, J.E., Root, P.J., 1963. The behavior of naturally fractured reservoirs. SPE J. 3, 245 – 255; Trans. AIME 228, 245 – 255.