A new upscaling technique based on Dykstra–Parsons coefficient: evaluation with streamline reservoir simulation

Journal of Petroleum Science and Engineering 40 (2003) 27 – 36 www.elsevier.com/locate/petrol

A new upscaling technique based on Dykstra–Parsons coefficient: evaluation with streamline reservoir simulation Ce´lio Maschio*, Denis Jose Schiozer Departamento de Engenharia de Petro´leo, Faculdade de Engenharia Mecaˆnica, UNICAMP, Caixa Postal 6122, 13.083-970 Campinas, Sa˜o Paulo, Brazil Received 30 May 2002; accepted 2 March 2003

Abstract The loss of information is inevitable in any upscaling technique. The efficiency of a given method must take in account two basic aspects: the first is the agreement of the results obtained with the coarse grid when compared to the results obtained with the fine grid and the second is the upscaling computational performance. Upscaling is really justifiable for very fine-scale reservoir models, usually with more than 1 million blocks or when many runs are necessary and the computational performance is very important. Due to the size of such fine-scale models, two problems exist. The first problem is that application of a pressure solver technique (or numerical method, which is generally more appropriate) on the upscaling is very time-consuming. The second problem is that the flux simulation of the fine-scale model, in order to validate a given upscaling technique, is very difficult and sometimes impossible through traditional simulators. In this work, an upscaling technique based on a heterogeneity coefficient (Dykstra – Parsons), which is as efficient and faster than numerical methods, is proposed. All tested fine grid reservoir models were modeled by a streamline simulator. The proposed technique was applied to three case studies obtaining good agreement between coarse and fine grids, taking in account several production parameters. D 2003 Elsevier Science B.V. All rights reserved. Keywords: Upscaling; Heterogeneity; Dykstra – Parsons; Reservoir simulation; Streamline simulation

1. Introduction The advances in reservoir characterization technologies, such as the modern seismic data acquisition facilities, seismic inversion, well logs, etc., as well as the use of modern geostatistical methods can provide models of field-scale porous media in which very

* Corresponding author. E-mail addresses: [email protected] (C. Maschio), [email protected] (D.J. Schiozer).

detailed data are incorporated so that complex and large geological models can be generated. Fine-scale geological models have a very detailed structure, often involving several million grid blocks. Such details are very important in the flow behavior, but even considering the advances in the computer and software technologies, there is a great distance between fine geological models and the size supported by traditional simulators. Generally, the number of simulation blocks for traditional simulators is in the range of 10,000 – 100,000. These numbers can increase for more complex models.

0920-4105/03/$ - see front matter D 2003 Elsevier Science B.V. All rights reserved. doi:10.1016/S0920-4105(03)00060-3

28

C. Maschio, D.J. Schiozer / Journal of Petroleum Science and Engineering 40 (2003) 27–36

Therefore, the use of upscaling techniques to adjust the fine grid model to the software and hardware capacities is imperative. The main objective of the upscaling methods is to coarsen the highly resolved geological models to levels of detail suitable for simulation. According to Scheibe and Yabusaki (1998), the works related to upscaling methods used to define equivalent properties can be classified into three main categories: analytical methods, numerical methods and empirical power-averaging approaches to estimate block effective or equivalent parameters. Ideally, such a method should provide the same response as the fine grid model, preserving the structure of porous medium and reproducing the flow behavior in simulation process. However, the loss of information is inevitable in any upscaling method. A complete equivalence between the heterogeneous blocks and the coarse homogenous block is, in practice, impossible. With this in mind, an upscaling method must have two desirable characteristics. First, the response obtained with the coarse grid must approximate the response obtained with the fine grid and, second, the method must be computationally efficient. The logical sense is that the total cost (in terms of time of computational resources consuming) involved in the realization of the upscaling (coarse grid obtaining) plus the cost of flow simulation with the coarse grid must be less than the cost of flow simulation study using the fine grid geological model, if such simulation is possible to be realized (Mukhopadhyay and Sahimi, 2000). There are, in the literature, several descriptions of upscaling methods. Important review papers give a good feedback on most of the methods until the date of the publications. The papers written by Renard and Marsily (1997) and Wen and Herna´ndez (1996) are good examples. In the context of upscaling procedure, an inherent problem is how to validate the techniques, and point out those more suitable, if the simulation of the fine grid model is sometimes prohibitive or impracticable. Using traditional simulators, such as those based on finite difference, is generally impractical to simulate a very large model (greater than 1 million blocks). Due to this problem, there are in the literature a scarcity of works comparing

results of upscaling techniques in the reservoir forecasting (or evaluation of production parameters) obtained with the fine grid and with the coarse grid systems derived from the studied methods. An alternative for this problem is streamline-based simulation, whose basic concepts are described in the next section.

2. Principles of streamlines simulation Streamline simulation method solves a threedimensional problem by decoupling it into a series of one-dimensional problems. Mathematical details can be found in Batycky et al. (1997), Thiele (1996) and Baker et al. (2002). Modern streamline simulators are fully 3D and account for gravity and fluid mobility change effects. Some other aspects such as changing well conditions, due to rate changes, infill drilling and producer– injector conversions are also accounted for. Streamline simulation uses IMPES numerical technique. The strategy consists of solving the pressure equation implicitly to compute a set of streamlines that represents flow in the reservoir. Each streamline represents a volumetric rate and acts as a one-dimensional grid running perpendicular to pressure contour in the reservoir. In this type of simulation, fluids move along streamlines, in contrast to conventional simulators, by which fluids are confined to the grid cells and move in an orthogonal direction to grid cell faces (Grinestaff, 1999). Streamline simulation is very suitable for big, incompressible or low compressibility models. Another aspect is that this technique is more suitable for displacement than for depletion processes, and situations not dominated by gravity segregation and/or phase transitions and strong change in mobility. It is also very appropriate for waterflood optimization.

3. Power average upscaling methods Power average is a general analytical upscaling method for absolute permeability. The power average is based on the fact that the block permeability must lie in between the arithmetic and harmonic averages of the

C. Maschio, D.J. Schiozer / Journal of Petroleum Science and Engineering 40 (2003) 27–36

cells within the blocks. The power average is written according to Eq. (1), as follows: Keq ¼

N 1X Kp N i¼1 i

! 1p

4. Upper and lower bounds Single averages, such as arithmetic, harmonic and geometric average, for example, are fast methods but do not consider any directional aspects of heterogeneities. An interesting set of combination of these simple averages (basically arithmetic and harmonic), proposed by Cardwell and Parsons and by Le Loc’h (Renard et al., 2000), takes in account some directional aspects of heterogeneity which are similar to numerical methods but are not time-consuming. They showed that the equivalent permeability in a given direction is bounded by two limits: the upper limit (Cmax) is calculated through the harmonic mean of the arithmetic means of the local permeability, and the lower limit (Cmin) is calculated through the arithmetic mean of the harmonic means, in combined directions, as shown in Fig. 1.

Fig. 1. Upper and lower limits.

For a three-dimensional case and considering the x direction, the upper and lower limits are calculated as follows:

ð1Þ

were Ki is the permeability of the i-numbered element of the subdomain, N is the number of elements in subdomain and Keq is the averaged (equivalent) permeability of the subdomain. In Eq. (1), for p =  1, one obtains a harmonic average for the computation of equivalent permeability (Keq), and for p = 1, one obtains an arithmetic average. Geometric average corresponds to p ! 0. In conventional form, this method is used selecting the average exponent from the range  1.0 V p V 1.0.

29

x Cmax

2 !1 31 ny X nx nz X nx 4X 5 ¼ lx ðly ðlz ÞÞ ¼ k xx h h a ny nz i¼1 j¼1 k¼1 i;j;k ð2Þ

x Cmin

" #1 ny X nx nz X nx X xx 1 ¼ ðk Þ ¼ lza ðlya ðlxh ÞÞ ny nz k¼1 j¼1 i¼1 i;j;k ð3Þ

Similarly, for y and z directions: y Cmax ¼ lyh ðlxa ðlza ÞÞ;

y Cmin ¼ lza ðlxa ðlyh ÞÞ

ð4Þ

z ¼ lzh ðlya ðlxa ÞÞ; Cmax

z Cmin ¼ lxa ðlya ðlzh ÞÞ

ð5Þ

5. Dykstra –Parsons coefficient Dykstra – Parsons coefficient (DP) is a measurement of heterogeneity of a permeability data set. The main steps for calculating DP are as follows: (i) Arrange the permeability data in order of decreasing values; (ii) Calculate for each value the percent of values which have a greater permeability and express this number as percent greater than (cumulative probability), so that probability of the X value is P(x V X); (iii) Plot the data from previous step on log probability paper. Plot k on the log scale axis (coordinate) and the cumulative probability in abscissa; (iv) Fit the data to a straight line. Select from the fitted line the values of permeability that correspond to probabilities of 0.5 and 0.16; these values are named K50 and K16, respectively. If the data is not logarithmically distributed, prefer to fit better the data between probabilities of 0.16 and 0.84 (i.e.

30

(v) DP ¼

C. Maschio, D.J. Schiozer / Journal of Petroleum Science and Engineering 40 (2003) 27–36

one standard deviation to each side from the media); Finally, DP is calculated according to: K50  K16 K50

ð6Þ

K84  K50 K84

ð7Þ

or DP ¼

if the data is not logarithmically distributed. Dykstra – Parsons coefficient varies between 0 (homogeneous reservoir) and 1 (strongly heterogeneous reservoir). In this paper, Dykstra – Parsons coefficient in conjunction with Cmin and Cmax is used to develop a new method for upscaling absolute permeability. In this way, the ability of Cmin and Cmax to capture directional heterogeneities and DP to measure the heterogeneities level are used. Some authors have suggested the use of Cmin and Cmax to decide about the application of a pressure solver technique (so-called numerical method). In this context, Li et al. (1999), for example, concluded that if the upper (Cmax) and the lower (Cmin) bounds are considerably different, then a numerical method would be more appropriate. However, depending on the spatial distribution of the permeability values, the difference between Cmin and Cmax can be small, therefore not capturing the heterogeneity. Consider the example in Fig. 2, which represents the permeability in y direction for cases A and B, for which Cmin and Cmax are shown below: y ðCmin ÞA ¼ 36:9 mD

y ðCmin ÞB ¼ 6:2 mD

and y ðCmax ÞA ¼ 37:7 mD

y ðCmax ÞB ¼ 37:7 mD

The level of heterogeneity is the same for A and B (DPA = DPB = 0.767), but in case A, this heterogeneity would not be captured by the difference between Cmin

Fig. 2. Illustration of level of heterogeneity and difference in upper and lower bounds.

and Cmax. Therefore, a new method combining Cmin, Cmax and DP is proposed.

6. Proposed method According to Renard and Marsily (1997), Matheron defines a relation between two fundamental bounds by suggesting a formula where the effective permeability is a weighted average of Cmin and Cmax, as follows: Keq ¼ ðCmax Þa ðCmin Þ1a ; with aa½0; 1

ð8Þ

so that a=(D  1)/D, where D is the space dimension. According to Renard et al. (2000), Le Loc’h proposed to take the geometric mean of Cmin and Cmax to calculate the equivalent permeability (Keq): Keq ¼

pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Cmax Cmin

ð9Þ

More recently, Lunati et al. (2001) also have used this approach in the calculation of block effective permeability.

C. Maschio, D.J. Schiozer / Journal of Petroleum Science and Engineering 40 (2003) 27–36

31

ular-shaped models and is applicable for absolute permeability.

7. Applications The proposed method was applied in three models. Model A is a relatively small model. Model B is part of the Tenth SPE Comparative Solution Project—SPE 10 (Model 2) (Christie and Blunt, 2001) and Model C was built with some particular characteristics to show the application of the proposed method. 7.1. Model A Fig. 3. Obtaining equivalent permeability from heterogeneous blocks using DP, Cmax and Cmin.

Model A is composed by 36  36  6 blocks with a regular shape. The six layers are presented in Fig. 4.

In the present paper, an equation similar to Eq. (8) is proposed. The fundamental difference is the introduction of the DP coefficient as a measurement of heterogeneity in place of the factor a. The proposed formula for the equivalent permeability is: Keq ¼ ðCmax ÞDP ðCmin Þð1DPÞ

ð10Þ

According Eq. (10) above, if a set of fine blocks is homogeneous (or approximately homogeneous), DPi0, CminiCmax and, therefore, KeqiCminiCmax. The algorithm for upscaling using Eq. (10) consists basically of the following steps: (i) For a given fine blocks set, the DP coefficient is calculated; (ii) For Kx (permeability in x direction), upper and lower limits (Cmax and Cmin) are calculated, as Eqs. (2) and (3). The same is applied for Ky and Kz; (iii) With DP, Cmax and Cmin, the equivalent permeability (Keq) is then calculated according Eq. (10). Steps (i), (ii) and (iii) are applied for each fine grid blocks to be homogenized (Fig. 3). The proposed method can be used for three-dimensional and irreg-

Fig. 4. Layers of Model A.

32

C. Maschio, D.J. Schiozer / Journal of Petroleum Science and Engineering 40 (2003) 27–36

kPa bottom hole pressure and the injector has a injection rate of 600 m3/day. A coarse grid of size 9  9  2 (scale factor of  4 4  3) was created using all methods analyzed in this work. 7.2. Model B: Tenth SPE Comparative Solution Project (SPE 10) The fine-scale model size is 60  220  85. The x direction permeability is equal to y direction. Kz is different from x and y. A porosity map is also used, varying from practically 0 to 0.35. Four producer wells are used, each one in the corners of the models and an injector well in the center. More details for this model can be found in Christie and Blunt (2001). 7.3. Model C This fine-scale model has 240  120  96 cells (with 1,424,290 active simulation grid blocks). The model has an irregular shape and the x direction permeability (Kx) map is equal to y direction (Ky). The range of Kx and Ky varies from 0.5 to 800 mD. Permeability map in z direction (Kz) is different and varies from 0.5 to 80 mD. Fig. 5 shows the threedimensional (a) permeability distribution (Kx) for the model and three layers chosen at positions 10 (b), 45 (c) and 80 (d), taken from up to down on the model. These layers represent the overall characteristics of the reservoir. The two-dimensional (layers) distribution permeability illustrates the presence of irregular barrier (permeability lower than 6 mD). The porosity map varies from 0.05 to 0.3. The drainage strategy is presented in Fig. 6. It chooses a depletion strategy scheme with 22 producer

Fig. 5. Aspects of Model C.

The permeability is in the range of 5 – 2000 mD (channels characterized by light regions in Fig. 4). A constant porosity of 20% was considered. A producer and an injector, both vertical wells, are used. The production strategy is a quarter five-spot, with the producer on the lower left corner and the injector diagonally opposed. The producer operates at 4900

Fig. 6. Well positions for Model C.

C. Maschio, D.J. Schiozer / Journal of Petroleum Science and Engineering 40 (2003) 27–36

Fig. 7. Comparison of field oil rate curves for Model A.

wells and 12 injector wells (all vertical). We will use a well-by-well analysis for producers; this is labeled in the figure. The producer wells operate at 17,200 kPa bottom hole pressure and the injectors operate at an injection rate of 795 m3/day and a maximum bottom hole pressure of 35,000 kPa. The upscaling methods were applied to generate a coarse grid of size 40  20  16, that is, a reduction factor of 216 (6  6  6 in each direction). For Models B and C, the porosity was upscaled using an arithmetic average.

8. Results and discussion To evaluate the proposed method, a pressure solver and Le Loc’h were chosen as comparative methods. To benchmark such evaluation, the results from fine grid simulation using a fully 3D commercial streamline

Fig. 8. Comparison of field water cut for Model A.

33

Fig. 9. Comparison of field average pressure for Model A.

simulator were used. The well productivity indexes for coarse grid simulations were considered the same as for the fine grid simulations. The main results are presented in the next sections. In order to facilitate the comments, the proposed method will be referred to as ‘‘DP method’’ or ‘‘DP’’. 8.1. Model A Oil rate curves for Model A, obtained for the fine grid simulation and for the upscaling methods, are presented in Fig. 7. All methods exhibit similar efficiency compared to oil rate obtained with fine grid simulation. The three methods present a little underestimation of oil production. The field water cut for 2000 days of production is presented in Fig. 8. There is an overestimation of water production for all methods, but DP method shows a slightly better agreement with the fine grid solution compared to other methods.

Fig. 10. Comparison of field oil rate for Model B.

34

C. Maschio, D.J. Schiozer / Journal of Petroleum Science and Engineering 40 (2003) 27–36

Fig. 11. Comparison of field water cut for Model B.

Fig. 13. Comparison of field oil rate for a lower coarse grid (Model B).

Analyzing the average reservoir pressure (Fig. 9), DP method clearly presents better results. It was better, inclusive, than the pressure solver. Le Loc’h method provided the higher values for pressure.

method shows a little improvement compared to Le Loc’h. Fig. 11 presents the comparison of field water cut for Model B. Again, Le Loc’h and DP methods provided better results, DP being slightly closer to the fine grid solution. For average field pressure curves (Fig. 12), there is a lower level of agreement for all methods, but in this case, Le Loc’h and pressure solver methods show a lower difference from the curve related to fine grid solution. The effect of scale factor in upscaling numerical method was also analyzed. A lower coarse model (20  73  28) was generated using a scale factor of 27 (3  3  3). Field oil rate for this situation is presented in Fig. 13. It can be seen that there was only a little improvement on the oil rate prediction compared to the 12  44  17 (scale factor 5  5  5 = 125) model.

8.2. Model B The field oil rate curves for Model B (SPE 10— Model 2) are presented in Fig. 10. The largest discrepancy in this case is related to pressure solver method. This is in agreement with the results presented by Christie and Blunt (2001), in which the worst method was the use of linear pressure gradient boundary conditions for oil rate predictions. This method significantly overestimated the oil flow rate. In our work, Le Loc’h and DP methods presented a good prediction for the field oil rate, but the DP

Fig. 12. Field average pressure for Model B.

Fig. 14. Comparison of field oil rate for Model C.

C. Maschio, D.J. Schiozer / Journal of Petroleum Science and Engineering 40 (2003) 27–36

Fig. 15. Cumulative oil produced for Model C.

8.3. Model C The field oil rate for Model C is presented in Fig. 14. The overall level of agreement is relatively good for all methods. Until approximately 750 days of production, Le Loc’h method seriously underestimates the oil rate compared with other methods, but after this time, oil rate is closest to the fine grid solution. It is also observed that after 750 days, DP and pressure solver exhibit a slight overestimation of oil rate. To better visualize the difference in field oil production, a cumulative curve is shown in Fig. 15. As the Le Loc’h method showed the higher discrepancy in the period of higher production, the cumulative deviation was greater than that for DP and numerical methods. Fig. 16 shows the results for all 22 producer wells for Model C. The total oil produced by each well, as

35

the result of the simulation with the several coarse grids, was compared to the total oil produced by each well obtained from fine grid simulation. The deviation of each method is presented in Fig. 16. The first general observation is that the deviations from numerical method and DP method are very similar. The lower deviations (absolute value equal or below 5%) that can be observed from these two methods are related to most of the wells: P01, P02, P04, P07, P08, P12, P15, P16, P18, P19, P20 and P22. These methods resulted in a deviation greater than 10% only for four wells (P05, P06, P11 and P13). On the other hand, the deviation approximately equal to or greater than 10% associated to Le Loc’h results occurred for the wells P01, P02, P03, P04, P05, P06, P13, P15, P18, P20, P21 and P22, that is, most of producer wells. Another interesting observation is that for Le Loc’h method; the lower deviations are associated to wells P08, P10, P11, P12, P14, P16 and P19. In Fig. 6, it can be observed that these wells are in the central portion (look at top view) of the reservoir, and observing Fig. 5, one can see that in the central portion, the barrier concentration is lower, that is, for this case, Le Loc’h method was more efficient for the relatively more homogeneous regions. 8.4. Time for upscaling The main contribution of this work is that the proposed upscaling method provides good predictions very similar to a pressure solver, but the time required to execute the upscaling is much lower than the numerical method.

Fig. 16. Comparison of the deviation (related to fine grid solution) in oil produced by the 22 wells (Model C).

36

C. Maschio, D.J. Schiozer / Journal of Petroleum Science and Engineering 40 (2003) 27–36

The time required for upscaling using numerical method is strongly dependent on the scale factor. For example, the upscaling of Model B (SPE 10 case) from 60  220  85 to 20  73  28 (scale factor of 3  3  3 = 27) is about 4 times slower than for a scale factor of 125 (5  5  5). On the other hand, for the same model and using the proposed method, the time for upscaling from 60  220  85 to 20  73  28 is about 80 times lower than the tested numerical method and 200 times lower for upscaling from 60  220  85 to 12  44  17, and the increase on time for upscaling using the proposed method from a scale factor of 125 to a scale factor of 27 is only 25%.

9. Conclusions This work presented a fast and relatively efficient upscaling technique based on the Dykstra – Parsons heterogeneity coefficient, and upper and lower permeability limits. The method was tested in three models with critical heterogeneities (channels and barriers). The method uses the ability of Cmin and Cmax in capturing directional effects of heterogeneities and the ability of DP coefficient in measuring the level of heterogeneity. The efficiency of the method, in terms of reproducing the results from fine grid, was similar to the numerical method, but in terms of computational performance, the proposed method is much faster. In the case of SPE 10, for example, for a scale factor of 27, numerical method took 213 min and DP method took 0.84 min (about 250 times faster). Another important aspect treated in this work was the benchmark of the proposed method using the results obtained from the flux simulation of the fine grid models. This was possible through a fully 3D commercial streamline simulator. Because the models presented low compressibility (situation in which streamline simulator is more compatible with traditional simulators), the results obtained from fine grid solution consist of a reliable form of validating the upscaling technique.

Acknowledgements The authors would to thank Petrobras and FINEP for financial support of this project and the availability of the program UPA by Petrobras.

References Baker, R.O., Kuppe, F., Chugh, S., et al., 2002. Full-field modeling using streamline-based simulation: 4. Case studies. SPE Reservoir Evaluation & Engineering 5 (2), 126 – 134 (April). Batycky, R.P., Blunt, M.J., Thiele, M.R., 1997. A 3D field-scale streamline-based reservoir simulator. SPE Reservoir Engineering 12 (4), 246 – 254 (November). Christie, M.A., Blunt, M.J., 2001. Tenth SPE comparative solution project: a comparison of upscaling techniques. SPE Reservoir Evaluation & Engineering 4 (4), 308 – 317 (August). Grinestaff, G.H., 1999. Waterflood pattern allocations: quantifying the injector to producer relationship with streamline simulation. SPE 54616, SPE Western Regional Meeting, Anchorage, AK, 26 – 28 May. Li, D., Beckner, B., Kumar, A., 1999. A new efficient averaging technique for scaleup of multimillion-cell geologic models. SPE 56554, SPE Annual Technical Conference and Exhibition, Houston, TX, 3 – 6 October. Lunati, I., Bernad, D., et al., 2001. A numerical comparison between upscaling techniques: non-local inverse based scaling and simplified renormalization. Advances in Water Resources 24, 913 – 929. Mukhopadhyay, S., Sahimi, M., 2000. Calculation of the effective permeabilities of field-scale porous media. Chemical Engineering Science 55, 4495 – 4513. Renard, Ph., Marsily, G., 1997. Calculating equivalent permeability: a review. Advances in Water Resources 20 (5 – 6), 253 – 278. Renard, P., Le Loc’h, G., et al., 2000. A fast algorithm for the estimation of the equivalent hydraulic conductivity of heterogenous media. Water Resources Research 36 (12), 3567 – 3580. Scheibe, T., Yabusaki, S., 1998. Scaling of flow and transport behavior in heterogeneous groundwater systems. Advances in Water Resources 22 (3), 223 – 238. Thiele, M.R., Batycky, R.P., Blunt, M.J., Orr, F.M., 1996. Simulating flow in heterogeneous systems using streamtubes and streamlines. SPE Reservoir Engineering 11, 5 – 12 (February). Wen, X.H., Herna´ndez, J.G., 1996. Upscaling hydraulic conductivities in heterogeneous media: an overview. Journal of Hydrology 183, 9 – 32.