- Email: [email protected]
New upgridding method to capture the dynamic performance of the fine scale heterogeneous reservoir
Journal of Petroleum Science and Engineering 86-87 (2012) 225–236
Contents lists available at SciVerse ScienceDirect
Journal of Petroleum Science and Engineering journal homepage: www.elsevier.com/locate/petrol
New upgridding method to capture the dynamic performance of the fine scale heterogeneous reservoir Mohammad Sharifi ⁎, Mohan Kelkar McDougall School of Petroleum Engineering, The University of Tulsa, 800 S Tucker Drive, Tulsa, OK 74104, United States
a r t i c l e
i n f o
Article history: Received 29 June 2011 Accepted 14 March 2012 Available online 27 March 2012 Keywords: upgridding upscaling vertical communication pressure equilibrium tight gas reservoir
a b s t r a c t Many tight gas reservoirs are characterized by large gross interval with significant presence of shale and interspersed sand. The geo-cellular models of these reservoirs, which honor the resolution of log data, contain large number of layers. Without significant upgridding of vertical layers, flow simulation of such reservoirs is very difficult. Ideally, the upgridding process should be conducted such that the dynamic performance of geo-cellular model is preserved in the upgridded model. One way to ensure efficient upgridding is to preserve fine scale pressure profile in the coarse scale model. In this study we observed that, in heterogeneous reservoirs, the vertical communication across layers has a big impact on the flow behavior in a fine scale model. In the presence of vertical communication, fine scale heterogeneities are homogenized. That is, layers with different permeabilities have similar pressure profiles during the depletion phase. Hence by combining these layers into a single layer would still predict a similar dynamic performance. In contrast, when crossflow across layers does not exist, combining those layers would result in much different dynamic behavior than fine scale behavior. It was observed that upgridding performance is not only influenced by the extent of sand connectivity across the layers (which tells us about vertical communication) but also how areally dispersed the connections are present. In the case of existence of dispersed shale within reservoir it was found that the vertical flow is more helped by smaller sized vertical connections between layers, which are areally dispersed, than large sized connections concentrated within a certain region. Analytical solution is provided to quantitatively prove the effect of size as well as dispersed nature of connection across the layers. Using the underpinnings of analytical solution, a new methodology is developed which combines layers based on sand connectivity (vertical permeability) across adjacent layers as well as their respective petrophysical properties (horizontal permeability). Our approach is analytical. However, it was shown that the analytical upgridding error we estimate closely mirrors the simulation error which represents the difference in the dynamic performance of fine versus coarse layers. We validate our method by applying it to both synthetic as well as field cases. Through these examples, we demonstrate that newly proposed method is superior to proportional as well as variance and flux based upgridding methods. © 2012 Elsevier B.V. All rights reserved.
1. Introduction Conducting computationally efficient and an accurate simulation of a reservoir model is the first priority of any simulation engineer. Over the last decade, the development of detailed three dimensional geologic descriptions has become common-place within asset and project teams worldwide. As the computational speed has improved, both geo-cellular and simulation models have become bigger. However computational requirements for geo-cellular models are always smaller than that of the corresponding flow simulation models. Hence, there is a need to upscale the geo-cellular models.
⁎ Corresponding author. E-mail address: mohammad-sharifi@utulsa.edu (M. Sharifi). 0920-4105/$ – see front matter © 2012 Elsevier B.V. All rights reserved. doi:10.1016/j.petrol.2012.03.016
A typical reservoir simulator can handle 10 6 to 10 7 grid blocks. The number of grid blocks is a function of the production mechanism, heterogeneity of reservoir and even post-processing of the simulation data for history matching purposes. It is increasingly important that we not only predict the history, but also, uncertainties in the future performance. This requires generating multiple reservoir models for history matching purposes. Creation of multiple models requires arriving at a compromise between a very detailed, history matched, single model versus multiple coarser models which are reasonably history matched and can predict the future performance uncertainties (Ates et al., 2005; Taware et al., 2011). That is, as long as coarse models are able to properly resolve production data, it is better to use coarser models which can capture the future uncertainties. Therefore, reducing the number of grid blocks in a geo-cellular model is very important. The key to upscaling the geological model is to retain sufficient amount of heterogeneities in the model so that
226
M. Sharifi, M. Kelkar / Journal of Petroleum Science and Engineering 86-87 (2012) 225–236
the dynamic performance of the fine scale model is captured. Determining the optimum number of grid blocks in the coarse scale model to capture the dynamic performance is an active area of research. In the literature, the process of reducing the number of reservoir layers is called upgridding and determining the properties of the generated grid is called upscaling. Different upscaling and upgridding procedures are appropriate in different situations. The ideal procedure to use on a particular problem depends on the simulation question being addressed, the production mechanism, and the level of detail that can be accommodated in the coarse model. For a problem involving primary production with only oil or gas being produced, the coarse model should correctly capture the effects of near well heterogeneity as well as the general large scale flow response of the reservoir. For scenarios involving displacement of oil by water or gas, it may be important to accurately capture the effects of main flow paths between injection and production wells. This may require the use of specialized gridding procedures (Durlofsky, 2005). Many upgridding and upscaling techniques are available in the literature (Fincham and Christensen, 2004; Stern, 2005). These methods vary from simple statistical methods to using advanced flow based simulation techniques (Durlofsky, 1996; Durlofsky et al., 1994; King, 2006). These methods can be divided in two main categories. Those methods that are based on static analysis of the reservoir properties are called static based methods and those that involve using flow simulation are called dynamic methods. 1.1. Static based methods Proportional and variance based methods are the two common static based methods used in the industry. The simplest approach is a proportional method. In this method, the numbers of grid blocks in vertical direction are combined in a fixed proportion without considering their underlying properties. As an example, for reducing 100 vertical layers to 20 vertical layers, we will simply combine five fine scale layers into one upscaled layer. Although this method is currently used in the industry, the dynamic performance of the upscaled model may not reproduce the fine scale performance (Hosseini and Kelkar, 2010). For overcoming this deficiency, some researchers used the variance based method. Testerman (1962) was the first one who implemented such a method. Li et al. (1995), Li and Beckner (2000, 2001), Stern and Dawson (1999)and more recently King et al. (2006) used variance based method for upgridding purposes. In all of these methods, layers with similar static characteristics are combined into a coarse layer. The logic of static based method is to preserve the variance or heterogeneity of the fine scale model by isolating layers which are dissimilar. The key assumption behind these methods is that key static properties influence the dynamic characteristics; therefore, preservation of heterogeneities would result in better reproduction of dynamic performance. However in many instances, this assumption does not hold true. As an example, in the case of a three layer model with permeabilities of 1, 10 and 30 md, the variance based method would combine 1 and 10 md layers first since the difference is smaller; however, simulation of a gas production indicates that combining 10 and 30 md layers would preserve the fine scale behavior better (Hosseini and Kelkar, 2010). 1.2. Dynamic based methods Dynamic based methods rely on simulating the fine scale model under realistic conditions and then determining the upgridding of this model so that dynamic performance is preserved. One of such methods is a flux based method that was suggested by Durlofsky (1996). In that method first the single phase flow in fine scale model is run and then the fluxes are output from each layer. Once the basic structure of coarse grid is specified, fine grid blocks with
similar fluxes are combined to reach an appropriate number of coarse grid blocks (Durlofsky et al., 1997). In general, the procedure tries to maintain the resolution in an area which has high fluxes and lump the layers with low fluxes. Flow based methods may provide a better upgridding procedure for a given well configuration, but the methods require using fine scale simulations which may be difficult when the fine scale model contains millions of cells. Further, the method may be optimal for a given well configuration, but would be sub-optimal under other well configuration scenarios. 2. Problem statement To summarize the current state of the art in upgridding, the static methods are computationally efficient but may not reproduce dynamic characteristics of the fine scale model well; whereas, the dynamic methods are computationally demanding and may not provide optimal upgridding without precisely knowing the well configuration as well as the flow processes under actual field conditions. Our goal in this research is to address the deficiencies in the current static and dynamic models. We want to develop a method which accounts for the dynamic characteristics of the fine scale model, yet is efficient and robust enough that it can be used for various well configurations as well as different flow processes. In this paper, we only concentrate on primary depletion process. In this research, our goal is to develop an optimum algorithm for upgridding to reduce the number of layers of fine scale geological models for building a simulation model. Instead of focusing on static properties, our focus is on the dynamic characteristic of the fine scale model (i.e. production profile) so that the upgridded model will capture the dynamic behavior of the fine scale model. The method has to be efficient so that it can be applied quickly for a large number of grid blocks. In addition, the method has to be dependent on the flow process. To summarize, we want to develop an upgridding process which has the following characteristics: (1) It can capture most important dynamic characteristics of the fine scale model; (2) It is robust and efficient so that it can be applied quickly and for various well configurations; (3) It would capture primary depletion process in actual field conditions. 3. Proposed method/hypothesis The main hypothesis in the current methodology is that by preserving the pressure distribution in the fine scale model, dynamic characteristics of the fine scale model can be preserved (i.e. fluxes). The potential or the pressure distribution determines the flow behavior; therefore, by developing an upgridding scheme which captures the pressure distribution, it can also capture the dynamic characteristics of the coarse scale model. To determine the pressure profile in the fine scale model, we intend to use an analytical method. By nature, any analytical method requires some model simplification. Therefore, second hypothesis in our research is that the assumptions we make for deriving the analytical model are valid for actual field conditions. An actual fine scale model consists of three dimensional distributions of grid blocks of various sizes. To determine the combination of the layers, we will need to determine the errors between the layers. The error is defined later. However, imagine a column of grid blocks. We take two adjacent grid blocks and calculate the error. The error is the difference between the dynamic behavior of an upscaled grid block (which combines the two fine scale grid blocks) versus dynamic behavior in the two fine scaled grid blocks. Since we are not conducting any simulation, we are representing the dynamic behavior by average pressure within fine scale versus coarse scale grid blocks. Once we
M. Sharifi, M. Kelkar / Journal of Petroleum Science and Engineering 86-87 (2012) 225–236
227
Fig. 1. Simple 3 layer model.
calculate the error for every pair of fine scale grid blocks within two adjacent layers, we add the error for all the fine scale grid blocks within two layers to calculate error between two layers. 3.1. Approach Reservoir simulation solves a system of nonlinear partial differential equations for pressure and saturation distributions in the reservoir. The driving force in any reservoir is pressure; this can be supplied by external forces such as injection wells or by internal energy of the reservoir fluids (e.g., expansion of fluids) itself during the primary depletion. In our derivation, the focus is on a simple, single phase flow. The single phase can be either an oil or gas. 3.2. No crossflow scenario In this section, we propose a method based on a single phase flow with two additional assumptions. We assume no crossflow between the layers (Kz = 0) and a pseudo steady state condition in the reservoir. We start our model based on a simple three layer reservoir with different petrophysical properties in each layer and a well at the center of the model. This is shown in Fig. 1. We define the error as the difference between the average pressures in fine layers during depletion compared to average pressure in coarse layer during the depletion. That is, if we combine layers 1 and 2 the error term (Error12) will be the difference between the average pressures of layers 1 and 2 in fine scale model with the average pressure that will be generated in a coarse model. It is worth mentioning that after combining any two layers, the horizontal and vertical permeabilities of a coarse layer are calculated by arithmetic and harmonic averages of corresponding fine layer model permeabilities respectively. After combining layers 1 and 2 or 2 and 3 we can define these errors mathematically as: ∞ ∞ error1; 2 ¼ ∫0 ðp1 −p12 Þdt þ ∫0 ðp2 −p12 Þdt
ð1Þ
∞ ∞ error2; 3 ¼ ∫0 ðp2 −p23 Þdt þ ∫0 ðp3 −p23 Þdt :
ð2Þ
All the integrals are calculated between time zero and ∞ that is our final simulation time. The reason we used infinity is because we wanted our analytical expression to be utilized in upgridding process. By considering only a limiting case, our solution becomes independent of time and hence easy to use. The subscripts 1, 2 and 3 represent the first, second and third layers respectively and subscripts 12 and 23 represent the combined layers 1 and 2, and 2 and 3 respectively. For pseudo-steady state flow and single phase liquid, Eqs. (1) and (2) can be solved analytically. We used the upper limit as infinity so that we can get an expression which is independent of time. The detailed derivation is shown in Appendix A. The final expression for the error can be written in a simple form: i.e.,: φ φ EIk;kþ1 ¼ k − kþ1 : K k K kþ1
ð3Þ
That is, the best combination of layers is the one which minimizes error in Eq. (3). The details are provided in Section 3.6. This result is in agreement with the result that was observed by Kim et al. (2010). Please note that in real 3-D reservoir models, there are no homogeneous layers. Therefore, the error will have to be calculated grid block by grid block in vertical direction and the errors for all the grid blocks within two adjacent layers have to be combined to calculate the total error between the two layers. 3.3. Analytical solution for partial communication scenario As mentioned before, Eq. (3) is based on no crossflow. In reality, reservoir rarely contains layers with no crossflow. Many tight gas reservoirs contain combination of shale and sand. Sand grid blocks may communicate vertically, whereas, shale will impose a barrier. Depending on the percentage of shale and sand, partial blocking in vertical direction is likely. Higher the shale percentage, less is the vertical communication. The two layer or three layer problem with partial vertical communication has no analytical solution without simplifying assumptions. In this section, we extend our expression to account for the cross flow effect. Here the effect of a barrier between two porous media in the presence of two small openings is examined. In Fig. 2 well_A is producing at a constant flow rate q = 1000 STBD and is located at
Fig. 2. Effect of connection distribution on flow across the barrier.
228
M. Sharifi, M. Kelkar / Journal of Petroleum Science and Engineering 86-87 (2012) 225–236
the distance of 400 ft from the barrier (D0 = 400 ft). The distances between openings 1 and 2 from the well_A are D1 and D2 respectively. The distance between points 1 and 2 is called DL. For solving this system analytically we simulate the connections as injection wells (Zhao and Thompson, 2002). Both flux and pressure at the connections are changing with time. Then pressure equation for the system that consists of five unknowns is solved analytically (pressures and flow rates at points 1 and 2 and bottom hole pressure at well _A). The detailed calculation procedure is given in Appendix B. We do not use the equation derived in Appendix B directly in our solution procedure. However, it is shown to illustrate the importance of opening in a barrier on the flow behavior of a well. This provides us with an analog when we are considering vertical communication through discontinuous shale barriers. A production result (from other side of barrier) from analytical solution is shown in Fig. 2. It can be seen that only by having two small connections along the barrier, at late time, about 40% of the production of the well_A is coming from the other side of the barrier through points 1 and 2. Also, when the distance between the points 1 and 2 is bigger (DL= 300 ft) in Fig. 2b, the amount of supply from other side of the barrier is greater than the case where the connection points are closer to each other (DL= 10 ft) in Fig. 2a. 3.4. Effect of crossflow in 3 layer model In the Section 3.3 we showed that existence of opening within the barrier can cause more homogenization. In this section we are going to see effect of the existence of connections between two layers that are separated by an impermeable layer in the middle. In Fig. 3 we have a three layer model including a well that operates under constant bottom hole pressure. This is similar to Fig. 1 except the middle layer is a shale barrier with two small openings which connect the two layers. So, the middle layer is not contributing to the flow except providing communication across top and bottom layers. The well is fully perforated. The top layer has a higher permeability than the bottom layer. The middle layer is shale except two connections through which gas can flow from the bottom layer into top layer. Intuitively, the top layer would be depleted faster and due to pressure difference between the top and the bottom layers, the gas will flow from the bottom layer into top layer. Fig. 3 shows the case where two connections are close to each other; whereas in Fig. 4 the distance between these connections is larger. Finite difference simulation results for both cases are shown in Fig. 5. We compare our results with coarse scale model when the two productive layers are combined into a single layer. It is clear that the error from combining layers in the presence of connection points is much smaller compared to the no vertical communication scenario. More than that, it can be seen that the result from upscaled model is closer to fine scale model when the distance between two connection points is larger. These results are in agreement to the analytical solution that was obtained in Appendix B.
Fig. 4. 3 layer model containing shale barrier and two distant connections.
behavior of combined layers (upgridded) would be similar to the fine scale model when significant vertical communication exists. This also sheds light on the correct upgridding process. If significant vertical communication exists, heterogeneous layers can be combined under the depletion process; whereas, two relatively similar layers may not be combined if vertical communication does not exist across those layers. Another important consideration in incorporating the vertical communication between the layers is how well distributed the vertical communication is across the layer. That is, if the two adjacent layers are communicating only in one part of the area between the two layers, the vertical communication is not as effective. On the other hand, if there are many discontinuous regions where vertical communication exists, it would be lot more effective in “homogenizing” the two layers. This is what we observed in our analytical solution as well as simple three layer model. We want to validate this observation by considering more complex situations. For simplicity, we built a three layer model with the top layer with 200 md and bottom layer with 10 md. The middle layer is a shale barrier except that we created few openings in that layer by assigning non-zero permeability and very small porosity in that barrier. Fig. 6 shows the permeable paths in the shale barrier. As can be seen, various configurations of openings in shale barrier are considered. In the first figure (top, left) no opening exists; in the second figure (top, middle) only one grid block within shale barrier has nonzero vertical permeability, in the third (top, right), a total of 36 grid blocks have nonzero vertical permeability. In the bottom left 16 grid blocks have nonzero vertical permeability for communication and are evenly distributed. In the bottom middle and bottom right there are 36 grid blocks with nonzero vertical permeability for vertical communication, except that the locations of those openings are at different places. Fig. 7 shows the effect of the variability in configuration on the error (difference between three layer fine scale and one layer coarse scale production). It is clear from the graph that both the extent and distribution of connections are important. For quantifying the effect of these two factors (shale extent and distribution) to the upscaling procedure, we used the RL concept that was originally proposed by Zapata and Lake (1989). They defined qffiffiffiffi ffi
RL ¼ HL
Kz Kh
and proposed that for RL greater than 10, the layers will
reach the vertical equilibrium. L in the RL equation represents the
3.5. Heterogeneous reservoir upgridding (partial crossflow) The analytical solution as well as the results from a simplified three layer model illustrates the importance of vertical communication. In effect, vertical communication “homogenizes” the static heterogeneity. Because the pressure equilibrates in vertical direction, the flow
Fig. 3. 3 layer model containing shale barrier and two close connections.
Fig. 5. Simulation result for upscaled model with fine scale model in the presence and absence of vertical communication.
M. Sharifi, M. Kelkar / Journal of Petroleum Science and Engineering 86-87 (2012) 225–236
229
Fig. 6. Schematic diagram of different connection configuration (kz ≠ 0).
characteristic length in horizontal direction. For layered model, it is the total length in horizontal direction. H is the thickness of the two layers combined, Kz is the harmonic average of the layer vertical permeabilities and Kh is the arithmetic average of the layer horizontal permeabilities. High value of RL signifies higher amount of vertical communication and vice versa. Note that L is not a grid block dimension. Therefore, for determining the vertical communication, we need to consider overall dimension of the layer. If ‘x’ and ‘y’ dimensions of the two layers are different, then we can consider square root of area as a representative length. We used the same concept with some modification. For applying this formula to each of the two vertically adjacent grid blocks, “RL” was calculated locally. In the formula, we used harmonic average for vertical permeability and arithmetic average for horizontal permeability of top and bottom grid blocks. The height also is the summation of height of both layers. Instead of ‘L’, we used the effective length concept. Effective length is calculated globally as a square root of the effective area. Effective area is the summation of areas of all of the grid blocks in that layer multiplied by an effective fraction. In the limit, if all the grid blocks in vertical direction are communicating, the fraction is one. If there is no communication, the fraction is zero. The procedure for calculating the effective fraction is discussed in the next step. The reason for using effective fraction is to ensure that we consider the distribution of vertical crossflow areas in the layers as well as their extent. For calculating the effective fraction, we assume that if the vertical connections between two layers are well distributed, they provide better
Fig. 7. Difference in oil production between coarse scale (different connection configuration) and the fine scale models.
communication. On the other hand, if the vertical communication is concentrated in one part of the layer, it does not provide an effective communication. To determine the effective vertical communication, we superpose a coarse grid on the fine grid and calculate the new fraction. For having a more realistic fraction number that represents both the extent and the distribution of connected sands within the layers, the procedure was applied for at least four different superimposed coarse configurations, i.e., different sizes of superimposed grids. In each step, the fraction is calculated and finally the minimum of those four calculated fractions was considered as an effective fraction. Formula for calculating effective fraction for layer ‘k’ is shown in Eq. (4).
f ef f ;k ¼
∑Nt i¼1 mi wi ∑Nt i¼1 wi
ð4Þ
In this equation ‘Nt’ is the total number of grid blocks in each layer and the weight (wi) assigned to a grid in a coarse model is inversely proportional to the number of fine scale connections in the coarse grid. The method is illustrated using a simple example. Assume our fine scale model has Nt grids whereas our coarse superposed model is 2 by 2. If within one of coarse grids that consists of Nt/4 fine scale grids, ‘F’ connected grids exist, then wi for all of the fine grids within the coarse block will be 4F/Nt. In the case of no connected fine grids in the coarse grid, we assign wi = 1 for all of fine grids inside that coarse grid. In addition to the weight factor, for accounting for the extent of connection we assign mi = 0 for non-connected grids and mi = 1 for connected grids. For example, let us consider two different connection distributions in Fig. 8a and c. In both cases the number of connected grid blocks is the same (real fraction is 0.11). Here we superimpose once two coarse grid blocks (Fig. 4a and c) and once 4 coarse grid blocks (Fig. 4b and d) on the top of the fine model. From Eqs. (5-a) to (5-d), for different configurations we get different fraction and based on our methodology we choose the minimum fraction as the effective fraction of the system. Let us explain the procedure for calculating of effective fraction of Fig. 8b. In this case, the fine scale model was divided to four coarse grids. In the top-left superimposed grid, four small openings exist out of nine fine grids. So for all the 9 grids wi = 4/9. In other three coarse grids due to nonexistence of
230
M. Sharifi, M. Kelkar / Journal of Petroleum Science and Engineering 86-87 (2012) 225–236
Fig. 8. Different effective fractions for different connection configuration.
connection wi = 1. Finally, as it was mentioned mi = 0 for all nonconnected fine grids. The result is shown in Eq. (5-b). The interesting point is that in the case of uniform distribution of connected grids (Fig. 8c and d), the effective fraction is higher than the case of concentrated distribution (Fig. 8a and b) and the fractions are equal to the real fraction. (Fig. 8a)
fi ¼
4 4 ∑i¼1 18 18 4 ∑i¼1 18 þ ∑18 i¼1 1
¼ 0:04
ð5 aÞ
(Fig. 8b)
fi ¼
∑9i¼1 49 þ ∑27 i¼1 1
2 ∑36 i¼1 18
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Layer area ðAk Þ f k Arithmatic average of KHk ; KH kþ1 Harmonic average of KZ k ; KZ kþ1 Hk þ Hkþ1 ð6Þ
Fig. 9 shows the relation between RL that was obtained from our calculation and the dynamic error for all three layer model configurations that were shown in Fig. 7. As it was mentioned before, error is the difference between cumulative production of coarse and fine scale models: t¼ts
∑4i¼1 49
2 ∑4i¼1 18
RLk;kþ1 ¼
Dynamic Error ¼ ∑t¼0 ðQ –C ðt Þ−Q F ðt ÞÞ ¼ 0:057
ð5 bÞ
(Fig. 8c)
fi ¼
In this step we can calculate the RL by using the effective fraction from previous step in Eq. (6).
¼ 0:11
ð5 cÞ
¼ 0:11:
ð5 dÞ
ð7Þ
where Q_C and Q_F are coarse and fine scale cumulative production, and ts is final simulation time. We would expect that better vertical communication would result in a smaller error. That means, the larger the value of RL (better vertical communication), the smaller would be the simulation error. We use the effective fraction for all the configurations in Fig. 6 using our analytical approach and plot that value against the simulation error. Fig. 9 demonstrates that the dynamic error in flow simulation is well presented by the static error based on our analytical approach. In other words we are able
(Fig. 8d) 4
fi ¼
∑i¼1 19 1 ∑36 i¼1 9
By this procedure, smaller fraction is calculated for the layers in which connections are not evenly distributed. As an actual example, the comparison between real and effective fractions of the real reservoir is discussed in the result section (Fig. 13). What this implies is that, more evenly distributed the communication between the two layers, closer the pressure will be within the layers.
Fig. 9. Effect of RL on production error on 6-different 3-layer model configuration.
M. Sharifi, M. Kelkar / Journal of Petroleum Science and Engineering 86-87 (2012) 225–236
231
Start Upgridding Number of initial Layers=Nz
Column by column error term for all the grid blocks ϕ ϕ E I i , j ,k ,k +1 = i , j , k − i , j , k + 1 K i, j,k K i, j ,k +1 Calculating the effective fraction fi and RL for each two adjacent layers using Equations 4 and 6
EF i , j ,k ,k +1 = EI i , j ,k ,k +1 (0.9 exp(−0.5 RLk , k +1 ) + 0.1)
ny
nx
EL k ,k +1 = ∑∑ EF i , j ,k ,k +1 i=1 j=1
Find Minimum Layer Error ) ( EL m −1 , m
Combine that layer and update error
EL
U
EL
m −1 , m
U m , m +1
= EL
= EL
m −1 , m
m , m +1
+ EL
m , m +1
+ E L m +1,m + 2
Nz=Nz-1
If Nz> Desired Coarse Layer Number
No End Fig. 10. Upgridding algorithm flowchart.
Yes
232
M. Sharifi, M. Kelkar / Journal of Petroleum Science and Engineering 86-87 (2012) 225–236
to correctly predict what type of dynamic error will result from upgridding fine scale layers. Finally we need to define the final error term (E_F) based on both preserving the pressure profile criteria (Eq. (3)) and crossflow effect (Eq. (6)). By running different synthetic models we found that the following formula can be considered as the representative of the final error term. EF i; j;k;kþ1 ¼ EIi; j;k;kþ1 0:9Exp −0:5RLk;kþ1 þ 0:1
ð8Þ
In this way, when RL = 0 final error is equal to the initial static error (Eq. (3)). The error gets smaller as RL increases; however, never reaches below 10% of the static error. 3.6. Upgridding algorithm In this section we summarize the workflow for the proposed upgridding procedure. (1) Basic error term is calculated between any two adjacent grid blocks in vertical direction. This error represents the maximum error possible if no vertical communication exists. φ i; j;k φi; j;kþ1 EIi; j;k;kþ1 ¼ − K i; j;k K i;j;kþ1
ð9Þ
(5) After we have calculated Nz-1 errors corresponding Nz vertical layers, we arrange them in ascending order. We select the first two layers combined as the one with the smallest error. The only errors updated are the ones which are below and above the combined layers. Assume we combine layers m and m + 1. After combining m and m + 1 layers, only the error terms adjusted are: u
ð11Þ
u
ð12Þ
EL m−1;m ¼ ELm−1;m þ ELm;mþ1 EL m;mþ1 ¼ ELm;mþ1 þ ELmþ1;mþ2
and the remaining error term will not change. This is done to ensure that prior errors are accounted for while making the decision about which next two layers to combine. (6) We continue the process regressively so that in each step on more vertical layer is eliminated by combining with another layer. This continues till we reduce our fine scale model enough for running reservoir simulation. The summary of upgridding procedure based on both initial error and crossflow effect is shown in Fig. 10. 4. Results 4.1. Synthetic examples
(2) Degree of vertical communication between all adjacent layers is calculated using the procedure that was discussed in Section 3.4. So for layer k and k + 1, RLk,k + 1 is calculated for every pair of vertically adjacent layers. (3) We update the error term (EFi,j,k,k + 1) using Eq. (8) for every pair of vertically adjacent grid blocks. (4) Since our upgridding scheme is based on combining layers (and not combining vertical grid blocks), we calculate the total error between layers k and k + 1 by using the following equation. It is a summation of all the adjacent vertical block errors over grid blocks in x and y directions. nx
ny
ELi; j;k;kþ1 ¼ ∑i¼1 ∑j¼1 EF i; j;k;kþ1
ð10Þ
Fig. 11 shows the results from different methods for upscaling from 50 layers model to a 10 layer coarse model. Here we assume that Kz = 0. We compare the fine scale cumulative production (50 layers) with 10 layer models generated using various upgridding methods. These methods include proportional layering, flux based upgridding and the proposed, new, method. It can be seen that the methodology based on Eq. (3) does better job in terms of having a production profile closer to the fine scale model. For this case the result from both proportional and flux based methods is almost the same and optimistic. Eq. (3) is based on no crossflow assumption. In Fig. 12 it can be seen that in the case of a perfect crossflow between the layers (Kz = Kh), all the methods will give almost the same result and it is
Fig. 11. Resulting of upgridding from 50 to 10 layers.
M. Sharifi, M. Kelkar / Journal of Petroleum Science and Engineering 86-87 (2012) 225–236
233
Fig. 12. Oil production comparison in the presence and absence of crossflow.
very close to the fine scale result. By comparing Fig. 12a and b it is clear that crossflow will help the reservoir appear as a homogeneous reservoir and the errors from all methods are very small. Using the formula for calculating the effective fraction for perfect crossflow Eq. (4) always gives the same answer and the effective fraction is equal to one. Because of that, the value of RL will be maximum (Eq. (6)) and the calculated final error will be minimum (Eq. (8)). In a real reservoir we have regions where crossflow exists, and regions where the crossflow is negligible. We test proposed upgridding methodology on the two real cases in the next section. 4.2. Real case example_1 In this section new methodology is applied for upgridding Vible reservoir. Vible is located in Green River Basin, Wyoming, USA. Because of the depositional character of the field, the average rock permeabilities are around 10− 3 md. The model has fifteen producing wells. The reservoir contains single gas phase and includes 53 × 69 × 228 grid blocks in x,y and vertical directions respectively. In this case the reservoir was upgridded from 228 to 28 layers. Before showing the production profile result, the result of calculating effective fraction and the comparison with real fraction is shown in Fig. 13. It shows that the fraction number obtained from our method is always less than or equal to real fraction of connected
Fig. 13. Comparison between real and effective fraction in Vible case.
grids according to their extent and distribution. This indicates an importance of distribution of sand connectivity across the layers. The result from different upgridding techniques can be seen in Fig. 14. It can be seen that the result from the new method is better than the result from other methods and it could capture the dynamic response of fine scale model better than the other methods. 4.3. Real case example_2 Second real case is a 33 × 33 × 575 model from one heterogeneous shaly reservoir (reservoir_A). The initial reservoir pressure was 4000 psi and the reservoir fluid was dead oil. All the wells are producing at 1200 psi bottom hole pressure control. Due to high heterogeneous nature of reservoir the number of layers was reduced from 575 layers to 68 layers. The result of different upgridding scheme is shown in Fig. 15. It can be seen here that again new method performance in prediction of cumulative oil production is better than the proportional and flux based method. It is worth mentioning that proposed method is analytical and hence does not require any flow simulation. This method is very efficient and takes just about 3 min for both real cases using a standard personal computer.
Fig. 14. Gas production of Vible reservoir.
234
M. Sharifi, M. Kelkar / Journal of Petroleum Science and Engineering 86-87 (2012) 225–236
dpave −J p ¼ ∫t0 ∫pinit dt: C t vp pave −pwf
ð6 aÞ
By solving Eq. (6-a), we can express pressure changes in the each layer. i.e., in the first layer " # −J 1 t p1 −pwf ¼ pinit −pwf exp ct V p1
ð7 aÞ
"
p1 −pwf
# −J 1 ∝ exp t : ct V p1
ð8 aÞ
Now we define the simulation error as the difference between the average pressure in fine layers that are to be combined, and the coarse grid combined layer average pressure between time zero and end of simulation time ts. Fig. 15. Oil cumulative production of reservoir_A.
5. Summary and conclusions In this study relatively simple but computationally efficient (no need for simulation) method was proposed for upgridding. We used analytical solution to determine the best layering solution so that pressures in fine scale will be closer to the coarse scale model during the depletion process. Further, we also considered the impact of discontinuous sand shale distribution on vertical communication between the layers. We observed that the vertical communication between two adjacent layers is a function of the total number of connected grid blocks as well as how they are distributed. We developed a method to calculate effective communication between two adjacent layers. We tested our methodology by comparing the new method with proportional as well as flux based methods. For both the synthetic as well as field cases, we observed that our method is superior to other methods. Appendix A For simple three layer model with different permeability and porosity and similar height, we can develop a simple solution for pressure difference between fine and coarse scale models in the following steps: Flow rate and productivity index of a pseudo steady state system can be expressed as Eqs. (1-a) to (3-a). (Dake, 1978) q¼ " J¼
ln
kh pave −pwf re 3 rw − 4 þ s
# kh re 3 −4 þ s ln rw
q ¼ J pave −pwf
ð1 aÞ
ð2 aÞ
1 ∂vp vp ∂p
ð3 aÞ
ð4 aÞ
Using Eqs. 3-a and 4-a pressure changes in the system as a function of time: −ct vp
dpave dvp ¼ ¼ q ¼ J pave −pwf dt dt
ð9 aÞ
ts ts error2; 3 ¼ ∫0 ðp2 −p23 Þdt þ ∫0 ðp3 −p23 Þdt
ð10 aÞ
Eq. (9-a) can be expanded in the following line: ts error1; 2 ¼ ∫0 p1 −pwf − p12 −pwf dt
ð11 aÞ
ts þ∫0 p2 −pwf − p12 −pwf dt : The reason that we were able to put absolute value outside of integral in Eqs. (9) and (10) is that the sign of each term does not change with time. The reason that sign would not change is because the average pressure in fine scale layers would either be smaller or larger than average pressure in the combined layer throughout the depletion process. For example we want to prove that {exp[−At]exp[− Bt]} always is negative or positive where A = J1/(ctVp1) and B = J12/(ctVp12). Regarding the fact that A and B are constants in time, above statement is always true. If A > B then the expression is always negative and is B > A then the expression is always positive. To ensure that we can obtain a solution which is independent of time, we change the upper limit of integral to infinity. This would allow us to get a solution independent of time (that can be easily implemented in our upgridding scheme) Substituting, (
" #
" #)∞ cV ct V p12 −J 1 −J 12 t p1 error1; 2 ¼ dt − dt exp exp −J 1 ct V p1 −J 12 ct V p12 0
(
" #
" #)∞ cV ct V p12 −J 2 −J 12 t p2 dt − dt exp exp þ −J 2 ct V p2 −J 12 ct V p12 0
We also know that we can express wellbore pressure as a function of productivity index by following Eqs. (4-a) to (7-a): ct ¼ −
ts ts error1; 2 ¼ ∫0 ðp1 −p12 Þdt þ ∫0 ðp2 −p12 Þdt
ð5 aÞ
ct V p1 ct V p12 ct V p2 ct V p12 φ1 φ1 þ φ2 þ ∝ − − − ¼ J K J J J K þK 1
12
2
12
1
1
2
φ φ þ φ2 φ1 K 2 −φ2 K 1 φ2 K 1 −φ1 K 2 ¼ þ :ð12 aÞ þ 2 − 1 K2 K1 þ K2 ðK 1 þ K 2 ÞK 1 ðK 1 þ K 2 ÞK 2 Two final terms always have opposite sign, so we can simplify φ them more, i.e.: φK kk − K kþ1 . kþ1 Generally speaking, the final result can be expressed in a very simple form: φ φ errork;kþ1 ¼ k − kþ1 K k K kþ1
ð13 aÞ
M. Sharifi, M. Kelkar / Journal of Petroleum Science and Engineering 86-87 (2012) 225–236
These flow rates can be calculated as:
Appendix B Here we determine the amount of contribution from the other side of a barrier in total production of well_A in Fig. 2. For simplicity and without loss of generality, we assume that reservoir properties on both part of reservoir are the same. In the absence of any connection point across the barrier it is very easy to solve this system analytically by using image method (Dake, 1978). However, having a connection through the barrier, makes the problem complicated. For this purpose we simulate each connection point by an imaginary well. Then instead of solving a very complicated geometry we consider our system as two regions (Zhao and Thompson, 2002). Region_I consists of well_A, its image and two injection well_1 and well_2. Region_II consists of two imaginary producing wells, well_1 and well_2. The amount of production from well_1 and 2 is equal to their injection rate. For better understanding let us explain the pressure drop at well_A that is located in Region_I (Eq. (1-b)). Pressure change at well_A is affected by producing from well_A itself, its image well and points 1 and 2. The reason is that we know that in absence of any connection (perfect barrier) pressure drop at well_A is summation of pressure drop caused by well_A that is working in infinite reservoir, plus the pressure drop caused by its image well with the same rate that is located at distance 2d1 from well_A (Dake, 1978). But we know that because of existence of two connections (imaginary wells), the real pressure drop should be less than summation of well_A and its image. The effect of these two imaginary wells is shown in third and fourth terms in Eq. (1-b). Third and fourth terms are pressure support from connection numbers one and two respectively. It is worth mentioning that neither pressure nor flow rate is not constant at the imaginary well_1 and well_2. Regarding the fact that both pressure and flow rate are changing at the wells_1 and 2, we write these equations by applying the Duhamel principle (Thompson and Reynolds, 1986). Here qc is flow rate at well_A that is constant. qv1 and qv2 are variable rates at points 1 and 2. Δpcu, I and Δp ' cu, I are pressure change and derivative of pressure change with time, assuming that wells are in Region_I. Similar to this, we can write Eqs. (2)–(5-b). Eq. (2-b) which show the pressure drop at point_1 assuming it is an injection well in Region_I. Eq. (3-b) shows the pressure drop at point_1 considering it as a production well in a second region. Eq. (4-b) shows the pressure drop at point_2 assuming it is an injection well in Region_I. Eq. (5-b) shows the pressure drop at point_2 considering it as a production well in the second region. ′
t
ΔpA ¼ qc Δpcu;I ðr w Þþqc Δpcu;I ð2d0 Þ−2∫0 qv1 Δp cu;I ðd1 ; t−τÞdτ−2 ð1 bÞ ′
t
×∫0 qv2 Δp cu;I ðd2 ; t−τÞdτ ′
t
′
t
Δpo1 ¼ −2∫0 qv1 Δp cu;I ðr w ; t−τ Þdτ−2∫0 qv2 Δp cu;I ðdL ; t−τ Þdτ
ð2 bÞ
þ2qc Δpcu;I ðd1 ; t Þ ′
t
′
t
Δpo1 ¼ 2∫0 qv1 Δp cu;II ðr w ; t−τ Þdτ þ 2∫0 qv2 Δp cu;II ðdL ; t−τÞdτ ′
t
′
t
Δpo2 ¼ −2∫0 qv2 Δp cu;I ðr w ; t−τ Þdτ−2∫0 qv1 Δp cu;I ðdL ; t−τ Þdτ
ð3 bÞ ð4 bÞ
þ2qc Δpcu;I ðd2 ; t Þ t
235
′
t
′
Δpo2 ¼ 2∫0 qv2 Δp cu;II ðr w ; t−τ Þdτ þ 2∫0 qv1 Δp cu;II ðdL ; t−τÞdτ
ð5 bÞ
We can simplify these five equations for calculating the five unknowns. After solving Eqs. (1-b)–(5-b) for the flow rate and applying the Laplace transform, we get the expression for flow through points_1 and 2 in Fig. 2 (Eqs. (6-b) and (7-b)). Here we are interested to see the amount of flow rate that is passing through these points.
n o qc Δpcu;I ðd2 Þ−Uqv1 Δpcu;II ðdL Þ þ Δpcu;I ðdL Þ n o ¼ U Δpcu;II ðr w Þ þ Δpcu;I ðr w Þ
qv2
qv1 ¼
qc Δpcu;I ðd1 Þ−Δpcu;I ðd2 Þ
Δpcu;II ðdL ÞþΔpcu;I ðdL Þ Δpcu;II ðrwÞþΔpcu;I ðrwÞ
fΔpcu;II ðdL ÞþΔpcu;I ðdL Þg 2 fΔpcu;II ðrw ÞþΔpcu;I ðrw Þg 2
1−
ð6 bÞ
ð7 bÞ
Where: Δpcu;I ðrw; U Þ ¼
sffiffiffiffi!! 5:615 U K rw hφct 4πηU 0 η
ð8 bÞ
I
Δpcu;II ðrw; U Þ ¼
sffiffiffiffi!! 5:615 U K rw hφct 4πηU 0 η
ð9 bÞ
II
η¼
0:006328 k : μct φ
ð10 bÞ
K0 is zero order modified Bessel function and U is Laplace variable. Finally by using Stephest algorithm, we can solve Laplace Eqs. (6b)-(10-b) in time domain for obtaining flow rate. The final result is shown in Fig. 2c for different distances between points_1 and 2. Nomenclature
Symbol
Definition
Ct D0 D1 D2 DL EF EI EL feff H k kh kz L m Nt Nx Ny Nz Pwf Pinit q Q_C Q_F qc qv rw ts U Vp w J η φ
Total compressibility Distance between well and barrier Distance between well and point 1 Distance between well and point 2 Distance between point 1 and point 2 Grid block Final error Grid block Static error Layer error Effective fraction Height permeability Horizontal permeability Vertical permeability Length Connection indicator Total number of grids in each layer Number of grid in x direction Number of grid in y direction Number of layer in z direction Bottom hole pressure Reservoir initial pressure flow rate Coarse scale cumulative production Fine scale cumulative production Constant flow rate Variable flow rate Wellbore radius Simulation time Laplace variable Pore volume Weighting factor Productivity index Diffusivity coefficient Porosity
236
M. Sharifi, M. Kelkar / Journal of Petroleum Science and Engineering 86-87 (2012) 225–236
References Ates, H., Bahar, A., El-Abd, S., Charfeddine, M., Kelkar, M., Datta-Gupta, A., 2005. Ranking and upscaling of geostatistical reservoir models using streamline simulation: a field case study. SPE Reserv. Eval. Eng. 8 (1), 22–32. Dake, L.P., 1978. Fundamental of reservoir engineering. Elsevier, Amsterdam. Durlofsky, L.J., 1996. Scale up of heterogeneous three dimensional reservoir descriptions. SPE J. 313–326. Durlofsky, L.J., 2005. Upscaling and gridding of fine scale geological models for flow simulation. Presented at the 8th International Forum on Reservoir Simulation (Stresa,Italy). Durlofsky, L.J., Behrens, R.A., Jones, R.C., Bernath, A., 1994. Application of a new scale up methodology to the simulation of displacement processes in heterogeneous reservoirs. SPE paper 28704, presented at the SPE International Petroleum Conference and Exhibition of Mexico, Veracruz (Oct 10–13). Durlofsky, L.J., Jones, R.C., Milliken, W.J., 1997. A nonuniform coarsening approach for the scale up of displacement processes in heterogeneous porous media. Adv. Water Resour. 335–347. Fincham, A., Christensen, F.A., 2004. Up-gridding from geological model to simulation model: review, applications and limitations. SPE paper 90921,presented at SPE Annual Technical Conference and Exhibition (September 26–29). Hosseini, S.A., Kelkar, M., 2010. Analytical upgridding method to preserve dynamic flow behavior. SPE Reserv. Eval. Eng. 473–484 (June). Kim, J.U., Datta-Gupta, A., Jimenez, E.A., Hohl, D., 2010. A dual scale approach to production data integration into high resolution geologic models. J. Pet. Sci. Eng. 147–159. King, M.J., 2006. Upgridding and upscaling: current trends and future directions. SPE 112810-DL, presented during Distinguished Lecture.
King, M.J., Burn, K.S., Wang, P., Muralidharan, V., Alvardo, F., Ma, X., Datta-Gupta, A., 2006. Optimal coarsening of 3D reservoir models for flow simulation. SPE Reserv. Eval. Eng. 317–334 (October 20–22). Li, D., Beckner, B., 2000. Optimal uplayering for scaleup of multimillion-cell geologic models. SPE Annual Technical Conference and Exhibition (October 1–4). Li, D., Beckner, B., 2001. A new efficient averaging technique for scaleup of multimillion-cell geologic models. SPE Reserv. Eval. Eng. 297–307. Li, D., Cullick, A.S., Lake, L.W., 1995. Global scale-up of reservoir model permeability with local grid refinement. J. Pet. Sci. Technol. 149–157. Stern, D., 2005. Practical aspects of scaleup of simulation models. J. Pet. Technol. 57 (9), 74–81. Stern, D., Dawson, A.G., 1999. A technique for generating reservoir simulation grids to preserve geologic heterogeneity. SPE paper 51942, presented at SPE Reservoir Simulation Symposium held in Houston, Texas (Feb 14–17). Taware, S., Friedel, T., Datta-Gupta, A., 2011. A practical approach for assisted history matching using grid coarsening and streamline-based inversion: experiences in a giant carbonate reservoir. SPE paper 141606, Presented at The Woodlands, Texas (February). Testerman, J.D., 1962. A statistical reservoir-zonation technique. J. Pet. Technol. 889–893. Thompson, L.G., Reynolds, A.C., 1986. Analysis of variable-rate well-test pressure data using Duhamel's principle. SPE formation evaluation, pp. 453–469. Zapata, V.J., Lake, L.W., 1989. A theoretical analysis of viscous crossflow. SPE Annual Technical Conference and Exhibition, San Antonio, Texas (October 4–7). Zhao, G., Thompson, L.G., 2002. Semi analytical modeling of complex-geometry reservoirs. SPE Reserv. Eval. Eng. 437–446.
Contents lists available at SciVerse ScienceDirect
Journal of Petroleum Science and Engineering journal homepage: www.elsevier.com/locate/petrol
New upgridding method to capture the dynamic performance of the fine scale heterogeneous reservoir Mohammad Sharifi ⁎, Mohan Kelkar McDougall School of Petroleum Engineering, The University of Tulsa, 800 S Tucker Drive, Tulsa, OK 74104, United States
a r t i c l e
i n f o
Article history: Received 29 June 2011 Accepted 14 March 2012 Available online 27 March 2012 Keywords: upgridding upscaling vertical communication pressure equilibrium tight gas reservoir
a b s t r a c t Many tight gas reservoirs are characterized by large gross interval with significant presence of shale and interspersed sand. The geo-cellular models of these reservoirs, which honor the resolution of log data, contain large number of layers. Without significant upgridding of vertical layers, flow simulation of such reservoirs is very difficult. Ideally, the upgridding process should be conducted such that the dynamic performance of geo-cellular model is preserved in the upgridded model. One way to ensure efficient upgridding is to preserve fine scale pressure profile in the coarse scale model. In this study we observed that, in heterogeneous reservoirs, the vertical communication across layers has a big impact on the flow behavior in a fine scale model. In the presence of vertical communication, fine scale heterogeneities are homogenized. That is, layers with different permeabilities have similar pressure profiles during the depletion phase. Hence by combining these layers into a single layer would still predict a similar dynamic performance. In contrast, when crossflow across layers does not exist, combining those layers would result in much different dynamic behavior than fine scale behavior. It was observed that upgridding performance is not only influenced by the extent of sand connectivity across the layers (which tells us about vertical communication) but also how areally dispersed the connections are present. In the case of existence of dispersed shale within reservoir it was found that the vertical flow is more helped by smaller sized vertical connections between layers, which are areally dispersed, than large sized connections concentrated within a certain region. Analytical solution is provided to quantitatively prove the effect of size as well as dispersed nature of connection across the layers. Using the underpinnings of analytical solution, a new methodology is developed which combines layers based on sand connectivity (vertical permeability) across adjacent layers as well as their respective petrophysical properties (horizontal permeability). Our approach is analytical. However, it was shown that the analytical upgridding error we estimate closely mirrors the simulation error which represents the difference in the dynamic performance of fine versus coarse layers. We validate our method by applying it to both synthetic as well as field cases. Through these examples, we demonstrate that newly proposed method is superior to proportional as well as variance and flux based upgridding methods. © 2012 Elsevier B.V. All rights reserved.
1. Introduction Conducting computationally efficient and an accurate simulation of a reservoir model is the first priority of any simulation engineer. Over the last decade, the development of detailed three dimensional geologic descriptions has become common-place within asset and project teams worldwide. As the computational speed has improved, both geo-cellular and simulation models have become bigger. However computational requirements for geo-cellular models are always smaller than that of the corresponding flow simulation models. Hence, there is a need to upscale the geo-cellular models.
⁎ Corresponding author. E-mail address: mohammad-sharifi@utulsa.edu (M. Sharifi). 0920-4105/$ – see front matter © 2012 Elsevier B.V. All rights reserved. doi:10.1016/j.petrol.2012.03.016
A typical reservoir simulator can handle 10 6 to 10 7 grid blocks. The number of grid blocks is a function of the production mechanism, heterogeneity of reservoir and even post-processing of the simulation data for history matching purposes. It is increasingly important that we not only predict the history, but also, uncertainties in the future performance. This requires generating multiple reservoir models for history matching purposes. Creation of multiple models requires arriving at a compromise between a very detailed, history matched, single model versus multiple coarser models which are reasonably history matched and can predict the future performance uncertainties (Ates et al., 2005; Taware et al., 2011). That is, as long as coarse models are able to properly resolve production data, it is better to use coarser models which can capture the future uncertainties. Therefore, reducing the number of grid blocks in a geo-cellular model is very important. The key to upscaling the geological model is to retain sufficient amount of heterogeneities in the model so that
226
M. Sharifi, M. Kelkar / Journal of Petroleum Science and Engineering 86-87 (2012) 225–236
the dynamic performance of the fine scale model is captured. Determining the optimum number of grid blocks in the coarse scale model to capture the dynamic performance is an active area of research. In the literature, the process of reducing the number of reservoir layers is called upgridding and determining the properties of the generated grid is called upscaling. Different upscaling and upgridding procedures are appropriate in different situations. The ideal procedure to use on a particular problem depends on the simulation question being addressed, the production mechanism, and the level of detail that can be accommodated in the coarse model. For a problem involving primary production with only oil or gas being produced, the coarse model should correctly capture the effects of near well heterogeneity as well as the general large scale flow response of the reservoir. For scenarios involving displacement of oil by water or gas, it may be important to accurately capture the effects of main flow paths between injection and production wells. This may require the use of specialized gridding procedures (Durlofsky, 2005). Many upgridding and upscaling techniques are available in the literature (Fincham and Christensen, 2004; Stern, 2005). These methods vary from simple statistical methods to using advanced flow based simulation techniques (Durlofsky, 1996; Durlofsky et al., 1994; King, 2006). These methods can be divided in two main categories. Those methods that are based on static analysis of the reservoir properties are called static based methods and those that involve using flow simulation are called dynamic methods. 1.1. Static based methods Proportional and variance based methods are the two common static based methods used in the industry. The simplest approach is a proportional method. In this method, the numbers of grid blocks in vertical direction are combined in a fixed proportion without considering their underlying properties. As an example, for reducing 100 vertical layers to 20 vertical layers, we will simply combine five fine scale layers into one upscaled layer. Although this method is currently used in the industry, the dynamic performance of the upscaled model may not reproduce the fine scale performance (Hosseini and Kelkar, 2010). For overcoming this deficiency, some researchers used the variance based method. Testerman (1962) was the first one who implemented such a method. Li et al. (1995), Li and Beckner (2000, 2001), Stern and Dawson (1999)and more recently King et al. (2006) used variance based method for upgridding purposes. In all of these methods, layers with similar static characteristics are combined into a coarse layer. The logic of static based method is to preserve the variance or heterogeneity of the fine scale model by isolating layers which are dissimilar. The key assumption behind these methods is that key static properties influence the dynamic characteristics; therefore, preservation of heterogeneities would result in better reproduction of dynamic performance. However in many instances, this assumption does not hold true. As an example, in the case of a three layer model with permeabilities of 1, 10 and 30 md, the variance based method would combine 1 and 10 md layers first since the difference is smaller; however, simulation of a gas production indicates that combining 10 and 30 md layers would preserve the fine scale behavior better (Hosseini and Kelkar, 2010). 1.2. Dynamic based methods Dynamic based methods rely on simulating the fine scale model under realistic conditions and then determining the upgridding of this model so that dynamic performance is preserved. One of such methods is a flux based method that was suggested by Durlofsky (1996). In that method first the single phase flow in fine scale model is run and then the fluxes are output from each layer. Once the basic structure of coarse grid is specified, fine grid blocks with
similar fluxes are combined to reach an appropriate number of coarse grid blocks (Durlofsky et al., 1997). In general, the procedure tries to maintain the resolution in an area which has high fluxes and lump the layers with low fluxes. Flow based methods may provide a better upgridding procedure for a given well configuration, but the methods require using fine scale simulations which may be difficult when the fine scale model contains millions of cells. Further, the method may be optimal for a given well configuration, but would be sub-optimal under other well configuration scenarios. 2. Problem statement To summarize the current state of the art in upgridding, the static methods are computationally efficient but may not reproduce dynamic characteristics of the fine scale model well; whereas, the dynamic methods are computationally demanding and may not provide optimal upgridding without precisely knowing the well configuration as well as the flow processes under actual field conditions. Our goal in this research is to address the deficiencies in the current static and dynamic models. We want to develop a method which accounts for the dynamic characteristics of the fine scale model, yet is efficient and robust enough that it can be used for various well configurations as well as different flow processes. In this paper, we only concentrate on primary depletion process. In this research, our goal is to develop an optimum algorithm for upgridding to reduce the number of layers of fine scale geological models for building a simulation model. Instead of focusing on static properties, our focus is on the dynamic characteristic of the fine scale model (i.e. production profile) so that the upgridded model will capture the dynamic behavior of the fine scale model. The method has to be efficient so that it can be applied quickly for a large number of grid blocks. In addition, the method has to be dependent on the flow process. To summarize, we want to develop an upgridding process which has the following characteristics: (1) It can capture most important dynamic characteristics of the fine scale model; (2) It is robust and efficient so that it can be applied quickly and for various well configurations; (3) It would capture primary depletion process in actual field conditions. 3. Proposed method/hypothesis The main hypothesis in the current methodology is that by preserving the pressure distribution in the fine scale model, dynamic characteristics of the fine scale model can be preserved (i.e. fluxes). The potential or the pressure distribution determines the flow behavior; therefore, by developing an upgridding scheme which captures the pressure distribution, it can also capture the dynamic characteristics of the coarse scale model. To determine the pressure profile in the fine scale model, we intend to use an analytical method. By nature, any analytical method requires some model simplification. Therefore, second hypothesis in our research is that the assumptions we make for deriving the analytical model are valid for actual field conditions. An actual fine scale model consists of three dimensional distributions of grid blocks of various sizes. To determine the combination of the layers, we will need to determine the errors between the layers. The error is defined later. However, imagine a column of grid blocks. We take two adjacent grid blocks and calculate the error. The error is the difference between the dynamic behavior of an upscaled grid block (which combines the two fine scale grid blocks) versus dynamic behavior in the two fine scaled grid blocks. Since we are not conducting any simulation, we are representing the dynamic behavior by average pressure within fine scale versus coarse scale grid blocks. Once we
M. Sharifi, M. Kelkar / Journal of Petroleum Science and Engineering 86-87 (2012) 225–236
227
Fig. 1. Simple 3 layer model.
calculate the error for every pair of fine scale grid blocks within two adjacent layers, we add the error for all the fine scale grid blocks within two layers to calculate error between two layers. 3.1. Approach Reservoir simulation solves a system of nonlinear partial differential equations for pressure and saturation distributions in the reservoir. The driving force in any reservoir is pressure; this can be supplied by external forces such as injection wells or by internal energy of the reservoir fluids (e.g., expansion of fluids) itself during the primary depletion. In our derivation, the focus is on a simple, single phase flow. The single phase can be either an oil or gas. 3.2. No crossflow scenario In this section, we propose a method based on a single phase flow with two additional assumptions. We assume no crossflow between the layers (Kz = 0) and a pseudo steady state condition in the reservoir. We start our model based on a simple three layer reservoir with different petrophysical properties in each layer and a well at the center of the model. This is shown in Fig. 1. We define the error as the difference between the average pressures in fine layers during depletion compared to average pressure in coarse layer during the depletion. That is, if we combine layers 1 and 2 the error term (Error12) will be the difference between the average pressures of layers 1 and 2 in fine scale model with the average pressure that will be generated in a coarse model. It is worth mentioning that after combining any two layers, the horizontal and vertical permeabilities of a coarse layer are calculated by arithmetic and harmonic averages of corresponding fine layer model permeabilities respectively. After combining layers 1 and 2 or 2 and 3 we can define these errors mathematically as: ∞ ∞ error1; 2 ¼ ∫0 ðp1 −p12 Þdt þ ∫0 ðp2 −p12 Þdt
ð1Þ
∞ ∞ error2; 3 ¼ ∫0 ðp2 −p23 Þdt þ ∫0 ðp3 −p23 Þdt :
ð2Þ
All the integrals are calculated between time zero and ∞ that is our final simulation time. The reason we used infinity is because we wanted our analytical expression to be utilized in upgridding process. By considering only a limiting case, our solution becomes independent of time and hence easy to use. The subscripts 1, 2 and 3 represent the first, second and third layers respectively and subscripts 12 and 23 represent the combined layers 1 and 2, and 2 and 3 respectively. For pseudo-steady state flow and single phase liquid, Eqs. (1) and (2) can be solved analytically. We used the upper limit as infinity so that we can get an expression which is independent of time. The detailed derivation is shown in Appendix A. The final expression for the error can be written in a simple form: i.e.,: φ φ EIk;kþ1 ¼ k − kþ1 : K k K kþ1
ð3Þ
That is, the best combination of layers is the one which minimizes error in Eq. (3). The details are provided in Section 3.6. This result is in agreement with the result that was observed by Kim et al. (2010). Please note that in real 3-D reservoir models, there are no homogeneous layers. Therefore, the error will have to be calculated grid block by grid block in vertical direction and the errors for all the grid blocks within two adjacent layers have to be combined to calculate the total error between the two layers. 3.3. Analytical solution for partial communication scenario As mentioned before, Eq. (3) is based on no crossflow. In reality, reservoir rarely contains layers with no crossflow. Many tight gas reservoirs contain combination of shale and sand. Sand grid blocks may communicate vertically, whereas, shale will impose a barrier. Depending on the percentage of shale and sand, partial blocking in vertical direction is likely. Higher the shale percentage, less is the vertical communication. The two layer or three layer problem with partial vertical communication has no analytical solution without simplifying assumptions. In this section, we extend our expression to account for the cross flow effect. Here the effect of a barrier between two porous media in the presence of two small openings is examined. In Fig. 2 well_A is producing at a constant flow rate q = 1000 STBD and is located at
Fig. 2. Effect of connection distribution on flow across the barrier.
228
M. Sharifi, M. Kelkar / Journal of Petroleum Science and Engineering 86-87 (2012) 225–236
the distance of 400 ft from the barrier (D0 = 400 ft). The distances between openings 1 and 2 from the well_A are D1 and D2 respectively. The distance between points 1 and 2 is called DL. For solving this system analytically we simulate the connections as injection wells (Zhao and Thompson, 2002). Both flux and pressure at the connections are changing with time. Then pressure equation for the system that consists of five unknowns is solved analytically (pressures and flow rates at points 1 and 2 and bottom hole pressure at well _A). The detailed calculation procedure is given in Appendix B. We do not use the equation derived in Appendix B directly in our solution procedure. However, it is shown to illustrate the importance of opening in a barrier on the flow behavior of a well. This provides us with an analog when we are considering vertical communication through discontinuous shale barriers. A production result (from other side of barrier) from analytical solution is shown in Fig. 2. It can be seen that only by having two small connections along the barrier, at late time, about 40% of the production of the well_A is coming from the other side of the barrier through points 1 and 2. Also, when the distance between the points 1 and 2 is bigger (DL= 300 ft) in Fig. 2b, the amount of supply from other side of the barrier is greater than the case where the connection points are closer to each other (DL= 10 ft) in Fig. 2a. 3.4. Effect of crossflow in 3 layer model In the Section 3.3 we showed that existence of opening within the barrier can cause more homogenization. In this section we are going to see effect of the existence of connections between two layers that are separated by an impermeable layer in the middle. In Fig. 3 we have a three layer model including a well that operates under constant bottom hole pressure. This is similar to Fig. 1 except the middle layer is a shale barrier with two small openings which connect the two layers. So, the middle layer is not contributing to the flow except providing communication across top and bottom layers. The well is fully perforated. The top layer has a higher permeability than the bottom layer. The middle layer is shale except two connections through which gas can flow from the bottom layer into top layer. Intuitively, the top layer would be depleted faster and due to pressure difference between the top and the bottom layers, the gas will flow from the bottom layer into top layer. Fig. 3 shows the case where two connections are close to each other; whereas in Fig. 4 the distance between these connections is larger. Finite difference simulation results for both cases are shown in Fig. 5. We compare our results with coarse scale model when the two productive layers are combined into a single layer. It is clear that the error from combining layers in the presence of connection points is much smaller compared to the no vertical communication scenario. More than that, it can be seen that the result from upscaled model is closer to fine scale model when the distance between two connection points is larger. These results are in agreement to the analytical solution that was obtained in Appendix B.
Fig. 4. 3 layer model containing shale barrier and two distant connections.
behavior of combined layers (upgridded) would be similar to the fine scale model when significant vertical communication exists. This also sheds light on the correct upgridding process. If significant vertical communication exists, heterogeneous layers can be combined under the depletion process; whereas, two relatively similar layers may not be combined if vertical communication does not exist across those layers. Another important consideration in incorporating the vertical communication between the layers is how well distributed the vertical communication is across the layer. That is, if the two adjacent layers are communicating only in one part of the area between the two layers, the vertical communication is not as effective. On the other hand, if there are many discontinuous regions where vertical communication exists, it would be lot more effective in “homogenizing” the two layers. This is what we observed in our analytical solution as well as simple three layer model. We want to validate this observation by considering more complex situations. For simplicity, we built a three layer model with the top layer with 200 md and bottom layer with 10 md. The middle layer is a shale barrier except that we created few openings in that layer by assigning non-zero permeability and very small porosity in that barrier. Fig. 6 shows the permeable paths in the shale barrier. As can be seen, various configurations of openings in shale barrier are considered. In the first figure (top, left) no opening exists; in the second figure (top, middle) only one grid block within shale barrier has nonzero vertical permeability, in the third (top, right), a total of 36 grid blocks have nonzero vertical permeability. In the bottom left 16 grid blocks have nonzero vertical permeability for communication and are evenly distributed. In the bottom middle and bottom right there are 36 grid blocks with nonzero vertical permeability for vertical communication, except that the locations of those openings are at different places. Fig. 7 shows the effect of the variability in configuration on the error (difference between three layer fine scale and one layer coarse scale production). It is clear from the graph that both the extent and distribution of connections are important. For quantifying the effect of these two factors (shale extent and distribution) to the upscaling procedure, we used the RL concept that was originally proposed by Zapata and Lake (1989). They defined qffiffiffiffi ffi
RL ¼ HL
Kz Kh
and proposed that for RL greater than 10, the layers will
reach the vertical equilibrium. L in the RL equation represents the
3.5. Heterogeneous reservoir upgridding (partial crossflow) The analytical solution as well as the results from a simplified three layer model illustrates the importance of vertical communication. In effect, vertical communication “homogenizes” the static heterogeneity. Because the pressure equilibrates in vertical direction, the flow
Fig. 3. 3 layer model containing shale barrier and two close connections.
Fig. 5. Simulation result for upscaled model with fine scale model in the presence and absence of vertical communication.
M. Sharifi, M. Kelkar / Journal of Petroleum Science and Engineering 86-87 (2012) 225–236
229
Fig. 6. Schematic diagram of different connection configuration (kz ≠ 0).
characteristic length in horizontal direction. For layered model, it is the total length in horizontal direction. H is the thickness of the two layers combined, Kz is the harmonic average of the layer vertical permeabilities and Kh is the arithmetic average of the layer horizontal permeabilities. High value of RL signifies higher amount of vertical communication and vice versa. Note that L is not a grid block dimension. Therefore, for determining the vertical communication, we need to consider overall dimension of the layer. If ‘x’ and ‘y’ dimensions of the two layers are different, then we can consider square root of area as a representative length. We used the same concept with some modification. For applying this formula to each of the two vertically adjacent grid blocks, “RL” was calculated locally. In the formula, we used harmonic average for vertical permeability and arithmetic average for horizontal permeability of top and bottom grid blocks. The height also is the summation of height of both layers. Instead of ‘L’, we used the effective length concept. Effective length is calculated globally as a square root of the effective area. Effective area is the summation of areas of all of the grid blocks in that layer multiplied by an effective fraction. In the limit, if all the grid blocks in vertical direction are communicating, the fraction is one. If there is no communication, the fraction is zero. The procedure for calculating the effective fraction is discussed in the next step. The reason for using effective fraction is to ensure that we consider the distribution of vertical crossflow areas in the layers as well as their extent. For calculating the effective fraction, we assume that if the vertical connections between two layers are well distributed, they provide better
Fig. 7. Difference in oil production between coarse scale (different connection configuration) and the fine scale models.
communication. On the other hand, if the vertical communication is concentrated in one part of the layer, it does not provide an effective communication. To determine the effective vertical communication, we superpose a coarse grid on the fine grid and calculate the new fraction. For having a more realistic fraction number that represents both the extent and the distribution of connected sands within the layers, the procedure was applied for at least four different superimposed coarse configurations, i.e., different sizes of superimposed grids. In each step, the fraction is calculated and finally the minimum of those four calculated fractions was considered as an effective fraction. Formula for calculating effective fraction for layer ‘k’ is shown in Eq. (4).
f ef f ;k ¼
∑Nt i¼1 mi wi ∑Nt i¼1 wi
ð4Þ
In this equation ‘Nt’ is the total number of grid blocks in each layer and the weight (wi) assigned to a grid in a coarse model is inversely proportional to the number of fine scale connections in the coarse grid. The method is illustrated using a simple example. Assume our fine scale model has Nt grids whereas our coarse superposed model is 2 by 2. If within one of coarse grids that consists of Nt/4 fine scale grids, ‘F’ connected grids exist, then wi for all of the fine grids within the coarse block will be 4F/Nt. In the case of no connected fine grids in the coarse grid, we assign wi = 1 for all of fine grids inside that coarse grid. In addition to the weight factor, for accounting for the extent of connection we assign mi = 0 for non-connected grids and mi = 1 for connected grids. For example, let us consider two different connection distributions in Fig. 8a and c. In both cases the number of connected grid blocks is the same (real fraction is 0.11). Here we superimpose once two coarse grid blocks (Fig. 4a and c) and once 4 coarse grid blocks (Fig. 4b and d) on the top of the fine model. From Eqs. (5-a) to (5-d), for different configurations we get different fraction and based on our methodology we choose the minimum fraction as the effective fraction of the system. Let us explain the procedure for calculating of effective fraction of Fig. 8b. In this case, the fine scale model was divided to four coarse grids. In the top-left superimposed grid, four small openings exist out of nine fine grids. So for all the 9 grids wi = 4/9. In other three coarse grids due to nonexistence of
230
M. Sharifi, M. Kelkar / Journal of Petroleum Science and Engineering 86-87 (2012) 225–236
Fig. 8. Different effective fractions for different connection configuration.
connection wi = 1. Finally, as it was mentioned mi = 0 for all nonconnected fine grids. The result is shown in Eq. (5-b). The interesting point is that in the case of uniform distribution of connected grids (Fig. 8c and d), the effective fraction is higher than the case of concentrated distribution (Fig. 8a and b) and the fractions are equal to the real fraction. (Fig. 8a)
fi ¼
4 4 ∑i¼1 18 18 4 ∑i¼1 18 þ ∑18 i¼1 1
¼ 0:04
ð5 aÞ
(Fig. 8b)
fi ¼
∑9i¼1 49 þ ∑27 i¼1 1
2 ∑36 i¼1 18
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi sffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Layer area ðAk Þ f k Arithmatic average of KHk ; KH kþ1 Harmonic average of KZ k ; KZ kþ1 Hk þ Hkþ1 ð6Þ
Fig. 9 shows the relation between RL that was obtained from our calculation and the dynamic error for all three layer model configurations that were shown in Fig. 7. As it was mentioned before, error is the difference between cumulative production of coarse and fine scale models: t¼ts
∑4i¼1 49
2 ∑4i¼1 18
RLk;kþ1 ¼
Dynamic Error ¼ ∑t¼0 ðQ –C ðt Þ−Q F ðt ÞÞ ¼ 0:057
ð5 bÞ
(Fig. 8c)
fi ¼
In this step we can calculate the RL by using the effective fraction from previous step in Eq. (6).
¼ 0:11
ð5 cÞ
¼ 0:11:
ð5 dÞ
ð7Þ
where Q_C and Q_F are coarse and fine scale cumulative production, and ts is final simulation time. We would expect that better vertical communication would result in a smaller error. That means, the larger the value of RL (better vertical communication), the smaller would be the simulation error. We use the effective fraction for all the configurations in Fig. 6 using our analytical approach and plot that value against the simulation error. Fig. 9 demonstrates that the dynamic error in flow simulation is well presented by the static error based on our analytical approach. In other words we are able
(Fig. 8d) 4
fi ¼
∑i¼1 19 1 ∑36 i¼1 9
By this procedure, smaller fraction is calculated for the layers in which connections are not evenly distributed. As an actual example, the comparison between real and effective fractions of the real reservoir is discussed in the result section (Fig. 13). What this implies is that, more evenly distributed the communication between the two layers, closer the pressure will be within the layers.
Fig. 9. Effect of RL on production error on 6-different 3-layer model configuration.
M. Sharifi, M. Kelkar / Journal of Petroleum Science and Engineering 86-87 (2012) 225–236
231
Start Upgridding Number of initial Layers=Nz
Column by column error term for all the grid blocks ϕ ϕ E I i , j ,k ,k +1 = i , j , k − i , j , k + 1 K i, j,k K i, j ,k +1 Calculating the effective fraction fi and RL for each two adjacent layers using Equations 4 and 6
EF i , j ,k ,k +1 = EI i , j ,k ,k +1 (0.9 exp(−0.5 RLk , k +1 ) + 0.1)
ny
nx
EL k ,k +1 = ∑∑ EF i , j ,k ,k +1 i=1 j=1
Find Minimum Layer Error ) ( EL m −1 , m
Combine that layer and update error
EL
U
EL
m −1 , m
U m , m +1
= EL
= EL
m −1 , m
m , m +1
+ EL
m , m +1
+ E L m +1,m + 2
Nz=Nz-1
If Nz> Desired Coarse Layer Number
No End Fig. 10. Upgridding algorithm flowchart.
Yes
232
M. Sharifi, M. Kelkar / Journal of Petroleum Science and Engineering 86-87 (2012) 225–236
to correctly predict what type of dynamic error will result from upgridding fine scale layers. Finally we need to define the final error term (E_F) based on both preserving the pressure profile criteria (Eq. (3)) and crossflow effect (Eq. (6)). By running different synthetic models we found that the following formula can be considered as the representative of the final error term. EF i; j;k;kþ1 ¼ EIi; j;k;kþ1 0:9Exp −0:5RLk;kþ1 þ 0:1
ð8Þ
In this way, when RL = 0 final error is equal to the initial static error (Eq. (3)). The error gets smaller as RL increases; however, never reaches below 10% of the static error. 3.6. Upgridding algorithm In this section we summarize the workflow for the proposed upgridding procedure. (1) Basic error term is calculated between any two adjacent grid blocks in vertical direction. This error represents the maximum error possible if no vertical communication exists. φ i; j;k φi; j;kþ1 EIi; j;k;kþ1 ¼ − K i; j;k K i;j;kþ1
ð9Þ
(5) After we have calculated Nz-1 errors corresponding Nz vertical layers, we arrange them in ascending order. We select the first two layers combined as the one with the smallest error. The only errors updated are the ones which are below and above the combined layers. Assume we combine layers m and m + 1. After combining m and m + 1 layers, only the error terms adjusted are: u
ð11Þ
u
ð12Þ
EL m−1;m ¼ ELm−1;m þ ELm;mþ1 EL m;mþ1 ¼ ELm;mþ1 þ ELmþ1;mþ2
and the remaining error term will not change. This is done to ensure that prior errors are accounted for while making the decision about which next two layers to combine. (6) We continue the process regressively so that in each step on more vertical layer is eliminated by combining with another layer. This continues till we reduce our fine scale model enough for running reservoir simulation. The summary of upgridding procedure based on both initial error and crossflow effect is shown in Fig. 10. 4. Results 4.1. Synthetic examples
(2) Degree of vertical communication between all adjacent layers is calculated using the procedure that was discussed in Section 3.4. So for layer k and k + 1, RLk,k + 1 is calculated for every pair of vertically adjacent layers. (3) We update the error term (EFi,j,k,k + 1) using Eq. (8) for every pair of vertically adjacent grid blocks. (4) Since our upgridding scheme is based on combining layers (and not combining vertical grid blocks), we calculate the total error between layers k and k + 1 by using the following equation. It is a summation of all the adjacent vertical block errors over grid blocks in x and y directions. nx
ny
ELi; j;k;kþ1 ¼ ∑i¼1 ∑j¼1 EF i; j;k;kþ1
ð10Þ
Fig. 11 shows the results from different methods for upscaling from 50 layers model to a 10 layer coarse model. Here we assume that Kz = 0. We compare the fine scale cumulative production (50 layers) with 10 layer models generated using various upgridding methods. These methods include proportional layering, flux based upgridding and the proposed, new, method. It can be seen that the methodology based on Eq. (3) does better job in terms of having a production profile closer to the fine scale model. For this case the result from both proportional and flux based methods is almost the same and optimistic. Eq. (3) is based on no crossflow assumption. In Fig. 12 it can be seen that in the case of a perfect crossflow between the layers (Kz = Kh), all the methods will give almost the same result and it is
Fig. 11. Resulting of upgridding from 50 to 10 layers.
M. Sharifi, M. Kelkar / Journal of Petroleum Science and Engineering 86-87 (2012) 225–236
233
Fig. 12. Oil production comparison in the presence and absence of crossflow.
very close to the fine scale result. By comparing Fig. 12a and b it is clear that crossflow will help the reservoir appear as a homogeneous reservoir and the errors from all methods are very small. Using the formula for calculating the effective fraction for perfect crossflow Eq. (4) always gives the same answer and the effective fraction is equal to one. Because of that, the value of RL will be maximum (Eq. (6)) and the calculated final error will be minimum (Eq. (8)). In a real reservoir we have regions where crossflow exists, and regions where the crossflow is negligible. We test proposed upgridding methodology on the two real cases in the next section. 4.2. Real case example_1 In this section new methodology is applied for upgridding Vible reservoir. Vible is located in Green River Basin, Wyoming, USA. Because of the depositional character of the field, the average rock permeabilities are around 10− 3 md. The model has fifteen producing wells. The reservoir contains single gas phase and includes 53 × 69 × 228 grid blocks in x,y and vertical directions respectively. In this case the reservoir was upgridded from 228 to 28 layers. Before showing the production profile result, the result of calculating effective fraction and the comparison with real fraction is shown in Fig. 13. It shows that the fraction number obtained from our method is always less than or equal to real fraction of connected
Fig. 13. Comparison between real and effective fraction in Vible case.
grids according to their extent and distribution. This indicates an importance of distribution of sand connectivity across the layers. The result from different upgridding techniques can be seen in Fig. 14. It can be seen that the result from the new method is better than the result from other methods and it could capture the dynamic response of fine scale model better than the other methods. 4.3. Real case example_2 Second real case is a 33 × 33 × 575 model from one heterogeneous shaly reservoir (reservoir_A). The initial reservoir pressure was 4000 psi and the reservoir fluid was dead oil. All the wells are producing at 1200 psi bottom hole pressure control. Due to high heterogeneous nature of reservoir the number of layers was reduced from 575 layers to 68 layers. The result of different upgridding scheme is shown in Fig. 15. It can be seen here that again new method performance in prediction of cumulative oil production is better than the proportional and flux based method. It is worth mentioning that proposed method is analytical and hence does not require any flow simulation. This method is very efficient and takes just about 3 min for both real cases using a standard personal computer.
Fig. 14. Gas production of Vible reservoir.
234
M. Sharifi, M. Kelkar / Journal of Petroleum Science and Engineering 86-87 (2012) 225–236
dpave −J p ¼ ∫t0 ∫pinit dt: C t vp pave −pwf
ð6 aÞ
By solving Eq. (6-a), we can express pressure changes in the each layer. i.e., in the first layer " # −J 1 t p1 −pwf ¼ pinit −pwf exp ct V p1
ð7 aÞ
"
p1 −pwf
# −J 1 ∝ exp t : ct V p1
ð8 aÞ
Now we define the simulation error as the difference between the average pressure in fine layers that are to be combined, and the coarse grid combined layer average pressure between time zero and end of simulation time ts. Fig. 15. Oil cumulative production of reservoir_A.
5. Summary and conclusions In this study relatively simple but computationally efficient (no need for simulation) method was proposed for upgridding. We used analytical solution to determine the best layering solution so that pressures in fine scale will be closer to the coarse scale model during the depletion process. Further, we also considered the impact of discontinuous sand shale distribution on vertical communication between the layers. We observed that the vertical communication between two adjacent layers is a function of the total number of connected grid blocks as well as how they are distributed. We developed a method to calculate effective communication between two adjacent layers. We tested our methodology by comparing the new method with proportional as well as flux based methods. For both the synthetic as well as field cases, we observed that our method is superior to other methods. Appendix A For simple three layer model with different permeability and porosity and similar height, we can develop a simple solution for pressure difference between fine and coarse scale models in the following steps: Flow rate and productivity index of a pseudo steady state system can be expressed as Eqs. (1-a) to (3-a). (Dake, 1978) q¼ " J¼
ln
kh pave −pwf re 3 rw − 4 þ s
# kh re 3 −4 þ s ln rw
q ¼ J pave −pwf
ð1 aÞ
ð2 aÞ
1 ∂vp vp ∂p
ð3 aÞ
ð4 aÞ
Using Eqs. 3-a and 4-a pressure changes in the system as a function of time: −ct vp
dpave dvp ¼ ¼ q ¼ J pave −pwf dt dt
ð9 aÞ
ts ts error2; 3 ¼ ∫0 ðp2 −p23 Þdt þ ∫0 ðp3 −p23 Þdt
ð10 aÞ
Eq. (9-a) can be expanded in the following line: ts error1; 2 ¼ ∫0 p1 −pwf − p12 −pwf dt
ð11 aÞ
ts þ∫0 p2 −pwf − p12 −pwf dt : The reason that we were able to put absolute value outside of integral in Eqs. (9) and (10) is that the sign of each term does not change with time. The reason that sign would not change is because the average pressure in fine scale layers would either be smaller or larger than average pressure in the combined layer throughout the depletion process. For example we want to prove that {exp[−At]exp[− Bt]} always is negative or positive where A = J1/(ctVp1) and B = J12/(ctVp12). Regarding the fact that A and B are constants in time, above statement is always true. If A > B then the expression is always negative and is B > A then the expression is always positive. To ensure that we can obtain a solution which is independent of time, we change the upper limit of integral to infinity. This would allow us to get a solution independent of time (that can be easily implemented in our upgridding scheme) Substituting, (
" #
" #)∞ cV ct V p12 −J 1 −J 12 t p1 error1; 2 ¼ dt − dt exp exp −J 1 ct V p1 −J 12 ct V p12 0
(
" #
" #)∞ cV ct V p12 −J 2 −J 12 t p2 dt − dt exp exp þ −J 2 ct V p2 −J 12 ct V p12 0
We also know that we can express wellbore pressure as a function of productivity index by following Eqs. (4-a) to (7-a): ct ¼ −
ts ts error1; 2 ¼ ∫0 ðp1 −p12 Þdt þ ∫0 ðp2 −p12 Þdt
ð5 aÞ
ct V p1 ct V p12 ct V p2 ct V p12 φ1 φ1 þ φ2 þ ∝ − − − ¼ J K J J J K þK 1
12
2
12
1
1
2
φ φ þ φ2 φ1 K 2 −φ2 K 1 φ2 K 1 −φ1 K 2 ¼ þ :ð12 aÞ þ 2 − 1 K2 K1 þ K2 ðK 1 þ K 2 ÞK 1 ðK 1 þ K 2 ÞK 2 Two final terms always have opposite sign, so we can simplify φ them more, i.e.: φK kk − K kþ1 . kþ1 Generally speaking, the final result can be expressed in a very simple form: φ φ errork;kþ1 ¼ k − kþ1 K k K kþ1
ð13 aÞ
M. Sharifi, M. Kelkar / Journal of Petroleum Science and Engineering 86-87 (2012) 225–236
These flow rates can be calculated as:
Appendix B Here we determine the amount of contribution from the other side of a barrier in total production of well_A in Fig. 2. For simplicity and without loss of generality, we assume that reservoir properties on both part of reservoir are the same. In the absence of any connection point across the barrier it is very easy to solve this system analytically by using image method (Dake, 1978). However, having a connection through the barrier, makes the problem complicated. For this purpose we simulate each connection point by an imaginary well. Then instead of solving a very complicated geometry we consider our system as two regions (Zhao and Thompson, 2002). Region_I consists of well_A, its image and two injection well_1 and well_2. Region_II consists of two imaginary producing wells, well_1 and well_2. The amount of production from well_1 and 2 is equal to their injection rate. For better understanding let us explain the pressure drop at well_A that is located in Region_I (Eq. (1-b)). Pressure change at well_A is affected by producing from well_A itself, its image well and points 1 and 2. The reason is that we know that in absence of any connection (perfect barrier) pressure drop at well_A is summation of pressure drop caused by well_A that is working in infinite reservoir, plus the pressure drop caused by its image well with the same rate that is located at distance 2d1 from well_A (Dake, 1978). But we know that because of existence of two connections (imaginary wells), the real pressure drop should be less than summation of well_A and its image. The effect of these two imaginary wells is shown in third and fourth terms in Eq. (1-b). Third and fourth terms are pressure support from connection numbers one and two respectively. It is worth mentioning that neither pressure nor flow rate is not constant at the imaginary well_1 and well_2. Regarding the fact that both pressure and flow rate are changing at the wells_1 and 2, we write these equations by applying the Duhamel principle (Thompson and Reynolds, 1986). Here qc is flow rate at well_A that is constant. qv1 and qv2 are variable rates at points 1 and 2. Δpcu, I and Δp ' cu, I are pressure change and derivative of pressure change with time, assuming that wells are in Region_I. Similar to this, we can write Eqs. (2)–(5-b). Eq. (2-b) which show the pressure drop at point_1 assuming it is an injection well in Region_I. Eq. (3-b) shows the pressure drop at point_1 considering it as a production well in a second region. Eq. (4-b) shows the pressure drop at point_2 assuming it is an injection well in Region_I. Eq. (5-b) shows the pressure drop at point_2 considering it as a production well in the second region. ′
t
ΔpA ¼ qc Δpcu;I ðr w Þþqc Δpcu;I ð2d0 Þ−2∫0 qv1 Δp cu;I ðd1 ; t−τÞdτ−2 ð1 bÞ ′
t
×∫0 qv2 Δp cu;I ðd2 ; t−τÞdτ ′
t
′
t
Δpo1 ¼ −2∫0 qv1 Δp cu;I ðr w ; t−τ Þdτ−2∫0 qv2 Δp cu;I ðdL ; t−τ Þdτ
ð2 bÞ
þ2qc Δpcu;I ðd1 ; t Þ ′
t
′
t
Δpo1 ¼ 2∫0 qv1 Δp cu;II ðr w ; t−τ Þdτ þ 2∫0 qv2 Δp cu;II ðdL ; t−τÞdτ ′
t
′
t
Δpo2 ¼ −2∫0 qv2 Δp cu;I ðr w ; t−τ Þdτ−2∫0 qv1 Δp cu;I ðdL ; t−τ Þdτ
ð3 bÞ ð4 bÞ
þ2qc Δpcu;I ðd2 ; t Þ t
235
′
t
′
Δpo2 ¼ 2∫0 qv2 Δp cu;II ðr w ; t−τ Þdτ þ 2∫0 qv1 Δp cu;II ðdL ; t−τÞdτ
ð5 bÞ
We can simplify these five equations for calculating the five unknowns. After solving Eqs. (1-b)–(5-b) for the flow rate and applying the Laplace transform, we get the expression for flow through points_1 and 2 in Fig. 2 (Eqs. (6-b) and (7-b)). Here we are interested to see the amount of flow rate that is passing through these points.
n o qc Δpcu;I ðd2 Þ−Uqv1 Δpcu;II ðdL Þ þ Δpcu;I ðdL Þ n o ¼ U Δpcu;II ðr w Þ þ Δpcu;I ðr w Þ
qv2
qv1 ¼
qc Δpcu;I ðd1 Þ−Δpcu;I ðd2 Þ
Δpcu;II ðdL ÞþΔpcu;I ðdL Þ Δpcu;II ðrwÞþΔpcu;I ðrwÞ
fΔpcu;II ðdL ÞþΔpcu;I ðdL Þg 2 fΔpcu;II ðrw ÞþΔpcu;I ðrw Þg 2
1−
ð6 bÞ
ð7 bÞ
Where: Δpcu;I ðrw; U Þ ¼
sffiffiffiffi!! 5:615 U K rw hφct 4πηU 0 η
ð8 bÞ
I
Δpcu;II ðrw; U Þ ¼
sffiffiffiffi!! 5:615 U K rw hφct 4πηU 0 η
ð9 bÞ
II
η¼
0:006328 k : μct φ
ð10 bÞ
K0 is zero order modified Bessel function and U is Laplace variable. Finally by using Stephest algorithm, we can solve Laplace Eqs. (6b)-(10-b) in time domain for obtaining flow rate. The final result is shown in Fig. 2c for different distances between points_1 and 2. Nomenclature
Symbol
Definition
Ct D0 D1 D2 DL EF EI EL feff H k kh kz L m Nt Nx Ny Nz Pwf Pinit q Q_C Q_F qc qv rw ts U Vp w J η φ
Total compressibility Distance between well and barrier Distance between well and point 1 Distance between well and point 2 Distance between point 1 and point 2 Grid block Final error Grid block Static error Layer error Effective fraction Height permeability Horizontal permeability Vertical permeability Length Connection indicator Total number of grids in each layer Number of grid in x direction Number of grid in y direction Number of layer in z direction Bottom hole pressure Reservoir initial pressure flow rate Coarse scale cumulative production Fine scale cumulative production Constant flow rate Variable flow rate Wellbore radius Simulation time Laplace variable Pore volume Weighting factor Productivity index Diffusivity coefficient Porosity
236
M. Sharifi, M. Kelkar / Journal of Petroleum Science and Engineering 86-87 (2012) 225–236
References Ates, H., Bahar, A., El-Abd, S., Charfeddine, M., Kelkar, M., Datta-Gupta, A., 2005. Ranking and upscaling of geostatistical reservoir models using streamline simulation: a field case study. SPE Reserv. Eval. Eng. 8 (1), 22–32. Dake, L.P., 1978. Fundamental of reservoir engineering. Elsevier, Amsterdam. Durlofsky, L.J., 1996. Scale up of heterogeneous three dimensional reservoir descriptions. SPE J. 313–326. Durlofsky, L.J., 2005. Upscaling and gridding of fine scale geological models for flow simulation. Presented at the 8th International Forum on Reservoir Simulation (Stresa,Italy). Durlofsky, L.J., Behrens, R.A., Jones, R.C., Bernath, A., 1994. Application of a new scale up methodology to the simulation of displacement processes in heterogeneous reservoirs. SPE paper 28704, presented at the SPE International Petroleum Conference and Exhibition of Mexico, Veracruz (Oct 10–13). Durlofsky, L.J., Jones, R.C., Milliken, W.J., 1997. A nonuniform coarsening approach for the scale up of displacement processes in heterogeneous porous media. Adv. Water Resour. 335–347. Fincham, A., Christensen, F.A., 2004. Up-gridding from geological model to simulation model: review, applications and limitations. SPE paper 90921,presented at SPE Annual Technical Conference and Exhibition (September 26–29). Hosseini, S.A., Kelkar, M., 2010. Analytical upgridding method to preserve dynamic flow behavior. SPE Reserv. Eval. Eng. 473–484 (June). Kim, J.U., Datta-Gupta, A., Jimenez, E.A., Hohl, D., 2010. A dual scale approach to production data integration into high resolution geologic models. J. Pet. Sci. Eng. 147–159. King, M.J., 2006. Upgridding and upscaling: current trends and future directions. SPE 112810-DL, presented during Distinguished Lecture.
King, M.J., Burn, K.S., Wang, P., Muralidharan, V., Alvardo, F., Ma, X., Datta-Gupta, A., 2006. Optimal coarsening of 3D reservoir models for flow simulation. SPE Reserv. Eval. Eng. 317–334 (October 20–22). Li, D., Beckner, B., 2000. Optimal uplayering for scaleup of multimillion-cell geologic models. SPE Annual Technical Conference and Exhibition (October 1–4). Li, D., Beckner, B., 2001. A new efficient averaging technique for scaleup of multimillion-cell geologic models. SPE Reserv. Eval. Eng. 297–307. Li, D., Cullick, A.S., Lake, L.W., 1995. Global scale-up of reservoir model permeability with local grid refinement. J. Pet. Sci. Technol. 149–157. Stern, D., 2005. Practical aspects of scaleup of simulation models. J. Pet. Technol. 57 (9), 74–81. Stern, D., Dawson, A.G., 1999. A technique for generating reservoir simulation grids to preserve geologic heterogeneity. SPE paper 51942, presented at SPE Reservoir Simulation Symposium held in Houston, Texas (Feb 14–17). Taware, S., Friedel, T., Datta-Gupta, A., 2011. A practical approach for assisted history matching using grid coarsening and streamline-based inversion: experiences in a giant carbonate reservoir. SPE paper 141606, Presented at The Woodlands, Texas (February). Testerman, J.D., 1962. A statistical reservoir-zonation technique. J. Pet. Technol. 889–893. Thompson, L.G., Reynolds, A.C., 1986. Analysis of variable-rate well-test pressure data using Duhamel's principle. SPE formation evaluation, pp. 453–469. Zapata, V.J., Lake, L.W., 1989. A theoretical analysis of viscous crossflow. SPE Annual Technical Conference and Exhibition, San Antonio, Texas (October 4–7). Zhao, G., Thompson, L.G., 2002. Semi analytical modeling of complex-geometry reservoirs. SPE Reserv. Eval. Eng. 437–446.