SPE 106679 Permeability Upscaling Techniques for Reservoir Simulation J. R. Villa, SPE, PDVSA Intevep; M. O. Salazar, SPE, Universidad Central de Venezuela

Copyright 2007, Society of Petroleum Engineers, Inc. This paper was prepared for presentation at the 2007 SPE Latin American and Caribbean Petroleum Engineering Conference held in Buenos Aires, Argentina, 15-18 April 2007. This paper was selected for presentation by an SPE Program Committee following review of information contained in an abstract submitted by the authors(s). Contents of the paper, as presented, have not been reviewed by the Society of Petroleum Engineers and are subject to correction by the author(s). The material, as presented, does not necessarily reflect any position of the Society of Petroleum Engineers, its officers, or members. Papers presented at SPE meetings are subject to publication review by Editorial Committees of the Society of Petroleum Engineers. Permission to copy is restricted to an abstract of not more than 300 words. Illustrations may not be copied. The abstract should contain conspicuous acknowledgement of where and by whom the paper was presented. Write Librarian, SPE, P.O. Box 833836, Richardson, TX 75083-3836, U.S.A., Fax 1-214-952-9435.

Abstract Upscaling reservoir properties for reservoir simulation is one of the most important steps in the workflow for building reservoir models. Upscaling allows taking high-resolution geostatistical models (107 -108 grid blocks) to coarse scale models (104 -105 grid blocks), manageable for reservoir simulation, while retaining the geological realism and thus effectively representing fluid transport in the reservoir 1,2. This work presents a study of the effectiveness of different available techniques for permeability upscaling and the implementation of a new technique for upscaling of relative permeability curves based on the numerical solution of a two-phase system and the Kyte and Berry method3. The reference fine scale model considered in this study is a conceptual fluvial reservoir based on the Stanford V model4. The reference fine scale isotropic and locally heterogeneous permeability distribution was upscaled to different upscaling ratios by means of analytical (static) and numerical single-phase (pressure solver, dynamic) techniques. Two-phase flow simulations were performed on the reference fine grid and upscaled models using a comercial black-oil simulator. Arithmetic, harmonic, and geometric averages were defined for static upscaling of the permeability distribution. The dynamic upscaling process considered one-phase and two-phase upscaling. One-phase upscaling considered upscaling of the permeability distribution and two-phase upscaling considered upscaling of the permeability distribution and relative permeability curves. Flow simulation results for waterflooding in the coarse scale model indicated relevant discrepancies with the fine grid re-

sults. Compared to fine-scale, flow results of the single-phase upscaling process indicated that the coarsest upscaled models did not match the water breakthrough times, water cut values, or well pressures from the reference model. The finer upscaled models reproduced the reference results more accurately than the coarser models. The two-phase dynamic upscaling technique implemented in this work resulted in the best match with the flow simulation results of the fine grid model. Results show that the most accurate upscaling scheme should be defined using the two-phase dynamic upscaling technique on the model with the smallest upscaling ratio.

Introduction Reservoir models generated by geostatistical techniques, highresolution fine scale models (107 -108 grid blocks), are capable of with great precision reservoir characterization as for compartmentalization, heterogeneity, connectivity and structure. However, the main drawback of high-resolution models is the significant computational cost when performing reservoir simulation. Upscaling reservoir properties allows taking high-resolution models to coarse scale models (104 -105 grid blocks) reducing computational costs during flow simulation for history matching and forecast. Permeability upscaling plays and important role in reservoir characterization5, as shown in Figure 1. The importance of using an appropriate upscaling technique consists in preserving the geological realism of highresolution, fine-scale models, thus preserving the flow response in reservoirs1,2. Permeability upscaling is an active research topic, numerous studies on upscaling have been conducted by university researchers6,7 and industry2,8. In this work, the effectiveness of different permeability upscaling techniques is evaluated using reservoir simulation. Analytical and numerical single-phase upscaling techniques were used with different upscaling ratios for a conceptual fluvial reservoir. Public available software was used for this purpose. In addition, an alternative technique based on two-phase numerical upscaling was developed and implemented. Using these techniques, flow simulation results of upscaled models were compared with the reference fine scale model in terms of flow

2

Permeability Upscaling Techniques for Reservoir Simulation

production, well pressure and saturation distribution. The sensitivity study of different permeability upscaling techniques is referred to a conceptual fluvial reservoir based on the Stanford V model4. As much for the fine scale model as upscaled models the numerical simulation of an oil-water inmiscible displacement process was effected, specifically a waterflooding process. This work includes the evaluation of static upscaling techniques (arithmetic, harmonic, and geometric) and dynamic upscaling techniques (one-phase and two-phase). The fine scale reference model was upscaled to four different upscaling ratios. This work allows to understand the impact of using different upscaling techniques on the simulation results. Investigation on gridding techniques was not considered in this study. Theory Darcy’s Law Darcy’s law express the relationship between fluid velocity and pressure gradient in a porous media. The Darcy velocity can be written in matrix notation for a Cartesian system (neglecting gravity) as: 1 u = − k · ∇p µ

(1)

In Equation 1, u is the velocity vector, k permeability tensor and ∇p the pressure gradient.

Permeability Tensor The permeability of the porous media is a property that can vary at any point and on any direction in the three-dimensional space 5. It is mathematically represented by the full permeability tensor k (Equation 2). Each component of the permeability tensor represents the directional permeability at one point in space. The permeability tensor is usually taken to be locally symmetric (kij = kji ).   kxx kxy kxz k =  kyx kyy kyz  (2) kzx kzy kzz If there is an orientation such that u and ∇p are parallel for a full tensor k, the principal orientation of permeability is obtained. This leads to the principal values and directions of the permeability tensor (eigenvalues and eigenvectors):  ∗  kxx 0 0 ∗ 0  k∗ =  0 kyy (3) ∗ 0 0 kzz

In two-phase incompressible flow, Darcy’s velocity is written as: kr uj = − j k · ∇p (6) µj where, krj is the relative permeability of phase j, and µj is the viscosity of phase j. The conservation equation is written as: ∇ · ut = 0

∇ (λt (s) k · ∇p) = 0

using the Darcy’s law, Equation 4 becomes: (5)

(8)

where λt is the total mobility, defined as: λt (s) =

krw kr + o µw µo

(9)

Classification of Techniques Upscaling techniques can be classified in terms of the parameters to be upscaled. In one-phase parameter upscaling technique, the fine-scale permeability tensor (k) is upscaled to a coarse-scale effective permeability tensor (k∗ ) while retaining the fine-scale relative permeabilities. Analytic and numerical methods are used for this purpose. On the other hand, in two-phase parameter upscaling technique, fine-scale relative permeabilities (krj ) are also upscaled to curves of different shapes (kr∗j ). These curves are usually referred as effective or pseudo-r elative permeability curves and their generation is accomplished by numerical methods. Figure 2 illustrates these upscaling techniques. Analytic methods Analytic methods for computing onephase parameter upscaling involve the solution of Equation 5 with no-flow boundary conditions in the non-communicative layered system shown in Figure 3, with n horizontal layers of permeability ki and dimensionless thickness hi . Flow in x direction is referred as parallel flow and upscaled permeability results in the arithmetic mean of permeability values in each layer: n X ∗ kxx = ki hi (10) i=1

Flow in the z direction is referred as series flow and results in the harmonic mean of permeability values in each layer: n X hi i=1

Governing Equations In single-phase incompressible flow, neglecting gravity and capillary effects, the conservation Equation is written as: ∇·u=0 (4)

(7)

subsituting Equation 6 in Equation 7 yields:

∗ kzz =

The tensor k∗ is the diagonal permeability tensor.

∇ · (k∇p) = 0

SPE 106679

!−1

ki

For this 2D system, the effective permeability tensor is:  ∗  kxx 0 ∗ k = ∗ 0 kyy

(11)

(12)

SPE 106679

J. R. Villa, M. O. Salazar

∗ ∗ where kxx and kyy are the principal values of the permeability tensor. Although permeability values were locally isotropic, effective permeability values in Equation 12 are constant and anisotropic, meaning that when upscaling fine-scale permeability values, coarse-scale permeabilities are anisotropic even with isotropic fine-scale permeability.

An estimate for systems with spatially random permeability can be obtained using the geometric average of permeability values where there is no particular assumption of flow direction: n

kg = exp

1X log ki n i=1

In systems with constant hi , Equation 13 pressed as: ! n1 n Y kg = ki

(13) 9

can be also ex-

(14)

i=1

A generalization of these averages techniques is known as the Power Law10. This empirical relation is written as: n

kω =

1X ω k n i=1 i

!1/ω (15)

The Power Law allows to determine the effective permeability kω of a set of volume elements for different values of ω ranging between -1 and 1. For values of ω = 1, 0, −1, Equation 15 results in arithmetic, geometric and harmonic averages, respectively. In general, arithmetic average provides an upper bound to k∗ , and the harmonic average provides a lower bound. Figure 4 shows the effective permeability computed for differents values of ω in a heterogeneous system10. Numerical methods Numerical methods involve the solution of Equation 5 or Equation 8 for the pressure distribution using finite-difference methods. In one-phase dynamic upscaling, each component of the diagonal permeability tensor is calculated separately depending of the flow direction. To solve the pressure distribution, arbitrary boundary conditions are assigned in the target coarse grid block (Figure 5) and the pressure of each fine grid block inside the target coarse grid block is computed by the solution of the single phase incompressible flow steady state Equation 5. This equation can be expressed in finite difference form and written as a matrix equation as: Tp = b

(16)

where, T is the transmissibility matrix, b a vector representing the source / sink term and p the unknown pressure vector. Once the pressure vector is calculated, the effective permeability of ∗ the target coarse grid block in the x direction (kxx ) is computed as follows: n

∗ kxx

ny

z X X nx = k1jk (p1jk − pin ) ny nz (pin − pout ) j=1

k=1

(17)

3

where, nx , ny , and nz are the number of fine grid blocks in the x, y, and z direction respectively, pin is the pressure in the inlet of the coarse gridblock, pout is the pressure in the outlet of the coarse grid block, k1jk is the absolute permeability in each fine grid block next to the inlet, and p1jk is the pressure in each ∗ ∗ fine grid block next to the inlet. The terms, kyy and kzz can be obtained in a similar fashion. In two-phase dynamic upscaling, it is recognized that it is not enough to upscale the absolute permeability to characterize transport in porous media under inmiscible displacement processes2. Therefore, the fine-grid relative permeability curves ∗ (krj ) are upscaled to different curves (krj ) and thus the fluidrock interaction in the coarse-scale model is considered. This kind of upscaling performed by generating pseudo-relative permeability curves allows to better represent fluid flow when the fine-grid relative permeability curves function are not able to represent. The main methods for generating pseudo-relative permeability functions 11 are Kyte and Berry method, Stone method, weighted porous volume method, weighted relative permeabilities method, and the Kirchoff’s Law method. All these methods use numerical flow simulation results on the high-resolution fine-scale model to generate the relative permeability curves of the coarse-scale model. Upscaling ratio The upscaling ratio is defined as: r=

n N

(18)

where, r is the upscaling ratio, n the number of fine grid blocks, and N the number of coarse grid blocks. The upscaling ratio represents a measure of how coarse is the coarse model. The larger upscaling ratio the coarser the upscaled model. Figure 6 shows two coarse models at different upscaling ratios.

Methodology The methodology propossed for this work involves the use of analytic and numerical methods to upscale a fine-grid model to four different upscaling ratios. Figure 7 illustrates the methodology followed in this work. Flow simulations are performed on the reference model and on the coarse-scale models and comparisons are made based on the calculated error. Analytic upscaling and single-phase dynamic upscaling of the finegrid model was performed using the Fortran-based program flowsim12. For two-phase dynamic upscaling, a new Fortranbased program flowsim2p was coded for this porpuse, following a similar structure to GSLIB suite of programs 12. Porosity upscaling was performed using program upscaler12 and program gsl2ecl 12 was used for output compatibility to the black-oil reservoir simulator ECLIPSE13. An automated workflow and post-processing of results were implemented using the application MATLAB14,15.

4

Permeability Upscaling Techniques for Reservoir Simulation

flowsim2p The computational tool developed in this work and termed flowsim2p was coded in Fortran9016 and implements the Kyte and Berry method for generating pseudo-relative permeability curves11,17. In the Kyte and Berry method, the fine-scale properties are weighted to obtain coarse-scale properties and later used in the Darcy’s law to obtain pseudo-relative permeabilities curves corresponding for each phase (j) and for each coarse-grid block: kr∗j = −

µ∗j qj∗ T ∗ ∆p∗j − ∆ρ∗j g∆D∗




SPE 106679

Cases Four cases were defined in this work. Each case corresponds to a different upscaling ratio. Different upscaling techniques were used for all cases. Table 1 shows the different cases proposed in this study. Case 1 is the finest upscaled model with 50x65x10 (32,500 grid blocks), Case 2 has 50x26x5 (6,500 grid blocks), Case 3 has 10x13x5 (650 grid blocks) and Case 4 is the coarsest upscaled model with 5x5x5 (125 grid blocks). Upscaling ratios r are 4, 20, 200, and 1,040 respectively.

(19)

where, kr∗j is the pseudo-relative permeability of a coarse-grid block for phase j, µ∗j the phase viscosity of a coarse-grid block, qj∗ the phase rate in one specific direction associated to a coarse grid block, T ∗ the weighted transmissibility between a coarse grid block and its neighboring grid block, ∆p∗j the pressure difference between the target coarse grid block and its neighboring grid block, ρ∗j the phase density associated to a coarse grid block, g the gravity term, and ∆D∗ the thickness difference between the target coarse grid block and its neighboring grid block. As shown in Equation 19, time-dependent phase pressure and saturation is computed in each fine-grid block inside the target coarse-grid block. Some drawbacks associated with the Kyte and Berry method can be found18,19. These include: problematic gridblocks due to flow restrictions, unability to generate flow direction-dependent curves, and high computational cost for global boundary conditions. In flowsim2p, an inmiscible oil-water displacement in a target coarse grid block is formulated and specific boundary conditions are imposed in the target coarse grid block (Figure 8). The displacement process is focused locally and the global boundary conditions are ignored, the oil phase has no dissolved gas, the gravitational and capilar effects are negligible, the rock and fluid compressibilities are dismissed, series flow is assumed between neighboring fine grid blocks, the fluid viscosities are taken constant, and the duration time of the process is variable depending of the coarse grid block size. Figure 9 shows the workflow of program flowsim2p. First, the input data is stored, then the program computes the pressure and saturation distribution of each phase inside the target coarse grid block, and finally the pseudo-relative permeability curves are generated. This process is repeated a certain time step number depending of the coarse-grid block size.

Reference fine-scale model In this work, the fine-scale model used is a conceptual fluvial reservoir based on the Stanford V model4. This is a Cartesian model with 130,000 grid blocks (100x130x10) as shown in Figure 10. The heterogeneous porosity and locally isotropic permeability distributions were generated using nonconditional sequential Gaussian simulation 20,21. The reservoir has four oil production wells and one water injector well. Permeability and porosity distributions are shown in Figures 11 and 12, respectively.

Results and Discussion A waterflooding process was performed in the fine-grid model and in each coarse-scale model using the black-oil reservoir simulator ECLIPSE13. This program is a fully-implicit, three dimensional, three phase, general purpose black-oil reservoir simulation. There are four vertical producing wells located at the reservoir boundaries. Water is injected through one vertical wells located in the middle on the reservoir. All wells are open to flow over the entire thickness. All wells exhibit no formation damage. Well constraints include maximum oil production of 20,000 STB/d and maximum bottom-hole pressure of 1050 psia at the injection well. No restriction were defined to the field production or well water cut. Simulation results were post-processed for analysis and visualization using MATLAB14. First, discussion of numerical simulation results is presented and later, an evaluation of the upscaling errors is considered. Simulation results Figures 13 to 17 shows the permeability field and water saturation distribution for different upscaling ratios (r) and at the lower layer of each corresponding model. In the case of two-phase dynamic upscaling (Figure 17), flow simulation of the coarse-scale models was performed with the correspondent pseudo-relative permeability curves. It is important to emphasize that the end-points of the pseudo-relative permeability curves are the same as the end-points of the relative permeability curves of the fine-scale model, since displacement efficiency of the reference fine-scale model is preserved3. The increase in the upscaling ratio, the difference in shape between the pseudo-relative permeability curves and the fine-scale relative permeability curves. When upscaling absolute permeability, the capacity of capturing heterogeneities is lost as the upscaling ratio increases, and the saturation profile loses representativety with respect to the reference fine-scale model. It is precisely the shape of the pseudo-relative permeability curves that compensates this lost of representativity. Also, this difference between the pseudo-curves and the original curves does not indicate that the pseudo-curves are better or worst, this fact only indicate the upscaling ratio used. For the cases of arithmetic, harmonic, geometric and one-phase dynamic upscaling (Figures 13 to 16) it is maintained the tendency of losing the representativity on the permeability field and saturation profiles as the upscaling ratio is increased.

SPE 106679

J. R. Villa, M. O. Salazar

Figures 18 to 22 show the bottomhole pressure and water cut at the four production wells for the differents techniques and upscaling ratios. Each well have early breakthrough times respect to the reference fine-scale model. Performance of the bottomhole pressure response is also. This difference less marked in the case of two-phase dynamic upscaling (Figure 22), which is the technique that presents the best fit of the flow response by well respect to reference fine-scale model. On the contrary, the harmonic static upscaling (Figure 19) is the technique that generates the flow response less representative. Figures 23 and 24 show the responses of field water cut and cumulative oil production, corresponding to different techniques and upscaling ratios. Once again, it can be observed how as it increases the upscaling ratio the flow responses become less representative of what is happening in the reference fine-scale model. For instance, in the cumulative oil production graphics for the models of greater upscaling ratios the technique of harmonic upscaling begins to fail remarkably on representing the flow performance. Figure 25 presents the CPU time required to perform the flow simulation on the studied models. Here, one motivation of performing upscaling is demonstrated by reducing the simulation time decreasing the number of grid blocks of the model. In fact, it can be observed that increasing on the magnitude order of the number of grid blocks generate increasings on the magnitude order of the simulation time. However, it should be established a balance between the decreasing of the simulation time and the representativity maintenance of the flow performance that happens on the reference fine-scale model.

Upscaling errors This sections presents an evaluation of the upscaling error computed for each technique used. Error is is defined by Equations 20 and 21:    ¯w P    S k,r     1 − k S wk,r    e= (20) nr and, P s¯wk,r =

i swi , ∀i  k nr

(21)

where r is the upscaling ratio, swi the water saturation of a fine grid block i inside a grid block k, swk,r the water saturation of a coarse grid block k at upscaling ratio r, s¯wk,r the weighted water saturation of fine grid blocks respect to a coarse grid block k, nr the number of coarse grid blocks at upscaling ratio r, and e the upscaling error. Equations 20 and 21 basically establish the difference between the block saturation on the upscaled model and the average fine-grid block saturations in the corresponding coarse-grid block. Figure 26 illustrates the procedure for the calculation of the upscaling error. Figures 27 to 31 show the error maps in gray scale for each one of the upscaling methods. Larger errors are located in grid blocks on the water-oil interfase of the waterflooding process,

5

indicating that upscaling errors increment in blocks with partial water and oil saturation. When two-phase dynamic upscaling method is used, errors are significantly reduced on the oilwater interfase (Figure 31), indicating that upscaling errors are reduced in those grid blocks. Table 2 shows the total upscaling errors in each case for different techniques and upscaling ratios. When two-phase dynamic upscaling method is used in Case 1, the smallest upscaling error is obtained. The harmonic static upscaling in the Case 4 generates larger errors. Figure 32 shows the upscaling errors for each case and method evaluated. In general, analytic and numerical upscaling techniques can be evaluated in terms of a total error. Figure 33 shows that numerical techniques generate smaller upscaling error that analytic techniques. The difference in the upscaling error can reach up to 10%, for the studied model, when upscaling ratio is increased. Conclusions Based on the oil and water production, well pressures, and saturation results obtained from the simulation results using dynamic and static upscaling, it can be see that by upscaling, earlier water breakthrough times result. According to the reproduction of the reference water cut, the most accurate dynamically upscaled result came from 50 × 65 × 10 model. Numerical upscaling techniques provides better prediction than does analytic upscaling. Dynamic upscaling will provide more accurate results relative to static for a given set of boundary conditions because different flow regimes can be accounted for in a single model. Lastly, reliable future predictions can only be obtained when geologic models depict the fine-scale case accurately and precisely. Regardless of the upscaling method, the averaging process will alter the original permeability field. The degree of difference between the resulting effective permeability and the true reference controls how different the grid block pressure and water saturation, well water cuts and bottom-hole pressures will be in the upscaled models relative to the fine scale. When analytic techniques were used, geometric static upscaling generated the best flow representativity. Upscaling relative permeabilities showed an improvement in upscaling results. Two-phase dynamic upscaling generated better flow response than the one-phase dynamic upscaling. An important aspect is the upscaling ratio, which determines the accuracy of production predictability. With a high upscaling ratio, the accuracy of the production prediction decreases. There is a limit on how coarse a model can be without introducing significant errors. This limit is important particularly when extreme permeability features are present in the reservoir, which directly affect the fluid flow. Results show that there is no an optimal upscaling technique suitable for any fine-scale model. For each fine-scale model,

6

Permeability Upscaling Techniques for Reservoir Simulation

a sensitivity study of different upscaling techniques with different upscaling ratios should be performed to determine the least upscaling error technique. Acknowledgements The authors would like to thank Universidad Central de Venezuela and PDVSA Intevep for the use of its computational facilities, and the Department of Petroleum Engineering at Stanford University for the use of its computational tools. Nomenclature u = velocity vector p = pressure vector s = saturation vector k = permeability tensor k∗ = diagonal permeability tensor T = transmissibility matrix b = source/sink vector ∇p = pressure gradient n = total number of fine grid blocks N = total number of coarse grid blocks ip = identificator number of fine grid blocks k = absolute permeability of a fine grid block h = thickness of a fine grid block λt = total mobility ω = power parameter kω = effective permeability obtained by Power Law kg = geometric average of permeability values pin = pressure at the inlet of a coarse grid block pout = pressure at the outlet of a coarse grid block p1,j,k = pressure on fine grid blocks next to the inlet k1,j,k = permeability on fine grid blocks next to the inlet vi,j,k = volume of a fine grid block ∆t = timestep size γ, m = Corey’s terms A = cross-sectional area to the flow References 1. Durlofsky, L. Advanced Reservoir Engineering PE222. Stanford University, California, USA, 2002. 2. Christie, M. Upscaling for Reservoir Simulation. Journal of Petroleum Technology, November 1996. 3. Westhead, A. Upscaling for Two-Phase Flow in Porous Media. California Institute of Technology, 2005. 4. Mao, S. and Journel, A. Generation of a Reference Petrophysical-Seismic Data Set: The Stanford V Reservoir. Technical report, Stanford University, California, USA, 1999. 5. Villa, J. R. Simulaci´on de Yacimientos. Universidad Central de Venezuela, Caracas, Venezuela, 2005. 6. Durlofsky, L. “Upscaling of Geocellular Models for Reservoir Flow Simulations: A Review of Recent Progress”. Paper presented at 7 International Forum on Reservoir Simulation, Germany, June 2003.

SPE 106679

7. Holden, L. and Nielsen, B. F. Global Upscaling of Permeability. 8. Stern, D. “Practical Aspects of Scaleup of Simulation Models”. Paper SPE 89032, 2005. 9. Ahmed, T. and McKinney, P. D. Reservoir Engineering Handbook. Gulf Professional Publishing, Houston, TX, USA, 2001. 10. Journel, A., Deutsch, C., and Desbarats, A. “Power Averaging for Block Effective Permeability”. Paper SPE 15128 presented at the 56 California Regional Meeting, California, USA, September 1986. 11. Cao, H. Evaluation of Pseudo Functions. Master’s thesis, Stanford University, 1988. 12. Deutsch, C. and Journel, A. GSLIB - Geostatistical Software Lybrary and User’s Guide. Oxford University Press, New York, USA, 1998. 13. Schlumberger. ECLIPSE Reference Manual, 2003. 14. The Mathworks, Inc. MATLAB, The Language of Technical Computing, 2005. 15. Hanselmanz, D. and Littlefield, B. Mastering Matlab 6: A Comprenhensive Tutorial and Reference. Prentice Hall, New Jersey, USA, 2001. 16. Chapman, J. FORTRAN 90/95 for Scientists and Engineers. McGraw-Hill, New York, USA, 2005. 17. Inanc, O. A Sensitivity Study on the Effectiveness of the Pseudo Relative Permeability Concept, PE224 Class Project. Technical report, Stanford University, California, USA, 2000. 18. Pickup, G. and Stephen, K. An Assessment of SteadyState Scale-Up for Small-Scale Geological Models. Technical report, Heriot-Watt University, 2000. 19. Barker, J. and Thibeau, S. “A Critical Review of the Use of Pseudo Relative Permeabilities for Upscaling”. Paper SPE 35491 presented at European 3-D Reservoir Modelling Conference, Stavanger, Norway, April 1996. 20. Deutsch, C. V. Geostatistical Reservoir Modeling. Oxford University Press, New York, USA, 2002. 21. Stanford University. Geostatistical Earth Modeling Software, 2004. 22. Kleppe, J. Reservoir Simulation. Technical report, Norwegian University of Science and Tecnology, January 2006. 23. Vaca, P. Simulaci´on de Yacimientos. Universidad Central de Venezuela, Caracas, Venezuela, 2003. 24. de la Garza, F. R. Simulaci´on Num´erica de Yacimientos. Technical report, PEMEX, M´exico, 2000. 25. Mattax, C. and Dalton, R. Reservoir Simulation. SPE, Richardson, Texas, USA, 1990. 26. Fanchi, J. Principles of Applied Reservoir Simulation. Butterworth-Heinemann, Texas, USA, 2001. 27. Crotti, M. and Cobenas, R. “Scaling Up of Laboratory Relative Permeability Curves. An Advantageous Approach Based on Realistic Average Water Saturations”. Paper SPE 69394 presented at the SPE Latin

SPE 106679

J. R. Villa, M. O. Salazar

7

American and Caribbean Petroleum Engineering Conference, Buenos Aires, Argentina, March 2001. Appendix A: flowsim2p Algorithm The first step flowsim2p executes is to read and store the input data defined in the parameter file, including porosity and permeability data (Figure A-1). Considering the fine-grid model showed in Figure A-2, an index number is assigned to each fine grid block: ip = (k − 1)nx ny + (j − 1)nx + i

(A-1)

where, ip is the index number, (i, j, k) the Cartesian coordinates of a fine grid block in x, y, and z, and nx , ny , nz the number of grid blocks in each direction. The next step is to calculate the effective permeability of the target coarse grid block by harmonic mean (Equation 11). Next, fine grid block pressures and saturations inside a coarse grid block are computed. The pressure equation is defined and solved by the IMPES method22 as shown in Equation A-2: α

α

[To + Tw ] pα+1 ± [bo + bw ] = 0

(A-2)

where, To and Tw are the transmissibility matrix for oil and water respectively, p the pressure vector, bo and bw the source/sink terms of oil and water respectively, α represents the previous timestep, α + 1 represents the current timestep, and 0 the null vector. In flowsim2p, the capillary effects are neglected and therefore po =pw =p. This equation is solved by the LSOR method 23,24. Once the pressure distribution inside the target coarse grid block is calculated, the phase saturation is calculated by the following equations25,26: To α po α+1 ± bo α = Tw α pw α+1 ± bw α =

 φvi,j,k  α+1 so − so α ∆t

(A-3)

 φvi,j,k  α+1 sw − sw α ∆t

(A-4)

where, φ is the porosity of a fine grid block, vi,j,k the volume of the respective fine grid block, ∆t the timestep size, and so , sw the phase saturation vectors. The phase relative permeability value of each fine grid block inside the target coarse grid block can be determinated by the Corey relations9: m1  sw − swc krw = γ1 (A-5) 1 − swc − sor

 kro = γ2

1 − sw − sor 1 − swc − sor

m2 (A-6)

where, krw and kro are the relative permeabilities of water and oil respectively, sw the water saturation in the current time step, swc the connate water saturation, and sor the residual oil saturation. The terms γ1 , γ2 , m1 , and m2 are real numbers determined from the original relative permeability curves and solving a linear system of equations. The flow rate in each fine grid block is determined using Darcy’s law for two-phase flow 27:     krj ∂p qj = − kA (A-7) µj ∂x where, (j) represents the respective phase, qj the phase flow rate at one specific direction, krj the phase relative permeability, µj the phase viscosity, k the absolute permeability, A the cross-sectional area to the flow, and ∂p/∂x the pressure gradient. The phase pressure difference and phase saturation of the target coarse grid block is determined by the following equations11,17: ∆p∗j = p∗jinlet − p∗joutlet p∗j


P krj k h pj
= P krj k h s∗j =

(A-8) (A-9)

P

v sj v∗

(A-10)

where, (j) represents the respective phase, ∆p∗j the effective phase pressure difference at the target coarse grid block, p∗jinlet and p∗joutlet the effective phase pressures at the inlet and outlet of the target coarse grid block respectively, krj the phase relative permeability of each fine grid block inside the target coarse grid block, k the absolute permeability of each fine grid block inside the target coarse grid block, h the thickness of each fine grid block inside the target coarse grid block, pj the phase pressure of each fine grid block at the respective face (inlet or outlet) inside the target coarse grid block, s∗j the effective phase saturation of the target coarse grid block, sj the phase saturation of each fine grid block inside the target coarse grid block, v the volume of each fine grid block inside the target coarse grid block, and v ∗ the volume of the target coarse grid block. flowsim2p calculates the phase difference pressure considering solely local boundary conditions and therefore only the fine grid blocks at the inlet and outlet of the target coarse grid block are used. If global boundary conditions were considered, it would take into account the central sides of two neighboring coarse blocks. Finally, the phase pseudo-relative permeabilities of the target coarse grid block are determined by Equation 19 adapted to the premises and suppositions of flowsim2p:

8

Permeability Upscaling Techniques for Reservoir Simulation

SPE 106679

µ∗ q ∗ D∗ = ∗ w ∗ w∗ x ∗ Dy Dz k ∆pw

(A-11)

µ∗o qo∗ Dx∗ Dz∗ k ∗ ∆p∗o

(A-12)

effective phase rate of the target coarse grid block calculated by the summation of the phase rates of the fine grid blocks at the central side of the target coarse grid block, and µ∗j the effective phase viscosity of the target coarse grid block calculated by volume weighted mean of the viscosities of the fine grid blocks.

where, kr∗j are the phase pseudo-relative permeabilities, Dx∗ , Dy∗ , Dz∗ , the coarse dimensions at each main direction, k ∗ the effective permeability of the target coarse grid block, ∆p∗j the phase pressure difference at the target coarse grid block, qj∗ the

This process is iteratively repeated for each time step. The number of time steps dependes on the size of the coarse grid block. Finally, an output file is generated (Figure A-3), containing the pseudo-relative permeability curves of a coarse-scale model for a defined upscaling ratio.

kr∗w

kr∗o =

Dy∗

SPE 106679

J. R. Villa, M. O. Salazar

Cases Fine: 100x130x10 Case 1: 50x65x10 Case 2: 50x26x5 Case 3: 10x13x5 Case 4: 5x5x5

Table 1: Cases Number Upscaling of blocks ratio 130,000 1 32,500 4 6,500 20 650 200 125 1,040

9

fine/coarse ni /Ni [1 1 1] [2 2 1] [2 5 2] [10 10 2] [20 26 2]

Figure 2: One-phase and two-phase parameter upscaling Table 2: Upscaling errors Upscaling method Arithmetic Harmonic Geometric Single-phase Two-phase

Case 1

Case 2

Case 3

Case 4

3.5% 3.9% 3.2% 2.9% 2.4%

12.3% 17.6% 9.2% 6.9% 5.2%

14.1% 19.9% 11.6% 10.9% 10.5%

18.1% 23.9% 17.2% 16.6% 15.8%

Figure 3: a) Parallel flow, b) Series flow

Figure 1: Reservoir modeling workflow

Figure 4: Power Law10

10

Permeability Upscaling Techniques for Reservoir Simulation

Figure 5: Flow scheme considered in one-phase dynamic upscaling

SPE 106679

Figure 8: Boundary conditions and flow scheme at the target coarse grid block considered on flowsim2p

Figure 6: Two reservoir models at different upscaling ratios

Figure 9: Workflow for flowsim2p

Figure 7: Methodology for upscaling the reference model

Figure 10: Permeability distribution and well locations in the reference fine-scale model

SPE 106679

J. R. Villa, M. O. Salazar

11

Figure 11: Histogram of permeability distribution in the reference fine-scale model

Figure 13: Permeability and water saturation distribution for different upscaling ratios - arithmetic static upscaling

Figure 12: Histogram of porosity distribution in the reference fine-scale model

Figure 14: Permeability and water saturation distribution for different upscaling ratios - harmonic static upscaling

12

Permeability Upscaling Techniques for Reservoir Simulation

SPE 106679

Figure 15: Permeability and water saturation distribution for different upscaling ratios - geometric static upscaling

Figure 17: Permeability and water saturation distribution for different upscaling ratios - two-phase dynamic upscaling

Figure 16: Permeability and water saturation distribution for different upscaling ratios - single-phase dynamic upscaling

Figure 18: Bottom-hole pressure and well water cut for different upscaling ratios using arithmetic upscaling

SPE 106679

J. R. Villa, M. O. Salazar

13

Figure 19: Bottom-hole pressure and well water cut for different upscaling ratios using harmonic upscaling

Figure 21: Bottom-hole pressure and well water cut for different upscaling ratios using single-phase upscaling

Figure 20: Bottom-hole pressure and well water cut for different upscaling ratios using geometric upscaling

Figure 22: Bottom-hole pressure and well water cut for different upscaling ratios using two-phase upscaling

14

Permeability Upscaling Techniques for Reservoir Simulation

SPE 106679

Figure 25: CPU time vs. number of blocks

Figure 23: Field water cut for different upscaling ratios and upscaling methods

Figure 24: Field cumulative oil production for different upscaling ratios and upscaling methods

Figure 26: Upscaling error calculation

SPE 106679

J. R. Villa, M. O. Salazar

15

Figure 27: Upscaling error distribution for different upscaling ratios using arithmetic upscaling

Figure 29: Upscaling error distribution for different upscaling ratios using geometric upscaling

Figure 28: Upscaling error distribution for different upscaling ratios using harmonic upscaling

Figure 30: Upscaling error distribution for different upscaling ratios using single-phase upscaling

16

Permeability Upscaling Techniques for Reservoir Simulation

SPE 106679

Figure 33: Upscaling error vs. upscaling ratio for analytic and numerical methods

Figure 31: Upscaling error distribution for different upscaling ratios using two-phase upscaling

Figure 32: Upscaling error vs. upscaling ratio for different upscaling methods

Figure A-1: Parameter file for flowsim2p

SPE 106679

J. R. Villa, M. O. Salazar

Figure A-2: Fine grid blocks arrange inside a coarse grid block - flowsim2p

Figure A-3: Output file of flowsim2p

17