Upscaling of grid properties in reservoir simulation

Upscaling, or homogenization, is substituting a heterogeneous property region consisting of fine grid cells with an equivalent homogeneous region made up of a single coarse-grid cell with an effective property value. (Equivalent in this case means either volume or flux vice, depending on the type of property that is to be upscaled.)

Upscaling is performed for each of the cells in the coarse grid and for each of the grid properties needed in the reservoir flow-simulation model. Therefore, the upscaling process is essentially an averaging procedure in which the static and dynamic characteristics of a fine-scale model are to be approximated by that of a coarse-scale model. A conceptual illustration of the upscaling process is shown in Fig. 1.

• 

  Fig. 1 – Concept of upscaling.

Can upscaling be avoided?

Typically, 3D geological models contain detailed descriptions of the reservoir that can be hard to capture properly with a significantly coarser model. Therefore, it would be preferable if upscaling could be avoided. Currently, an average-sized flow simulation model consists of approximately 100,000 active grid cells. This is to ensure that the CPU consumption of a simulation run will be reasonable (i.e., within practical limits).

Because a typical 3D geological model may consist of approximately 10 million active grid cells, it is obviously infeasible to run fluid-flow simulations directly on the geological model. Hence, upscaling is a required part of current reservoir modeling workflows.

Seen through the eyes of the geologist, the upscaling task may be a painful experience because all the geological details that were put into the model seem to be lost in the process. For a reservoir engineer, on the other hand, effective properties might be all that matter.

For volumetric (additive) properties such as porosity and saturation, the effective flow-cell value is simply given by the bulk and pore volume weighted arithmetic average, respectively, of the geo cells inside it. For the permeability, which is intrinsic (nonadditive) by nature, no such simple averaging method exists. The complexity one needs to take into account when upscaling permeability is considerable; therefore, all current techniques provide only an approximation of the true effective cell permeability. This approximation may range from very good to very poor, depending on the complexity of the fine-scale permeability distribution as well as the upscaling method used.

Upscaling techniques for absolute permeability

Homogenization of absolute permeability does not have an exact analytical solution, except for in a few idealized cases. The challenge of computing an accurate effective permeability has resulted in a large number of upscaling techniques. These techniques range from simple statistical averages to advanced numerical methods.

Tensor methods are the most accurate techniques available for computing the effective cell permeability. These are based on solving a second-order elliptic differential equation describing single-phase, incompressible, steady-state fluid flow in a region without sources and sinks (i.e., wells). Some flow-based methods may provide a full permeability tensor. However, because most multiphase flow simulators can only handle a diagonal permeability tensor because of the use of a seven-point stencil in 3D, diagonal tensor methods are most frequently used whether directly or indirectly (through a diagonalization of a full tensor). For a diagonal tensor, only the effective permeability in the principal directions of flow (x, y, and z) will be nonzero.

The flow equation is usually discretized with a finite-difference scheme, although finite-element methods are also applied occasionally. To compute all the directional components of the permeability tensor, the discretization and solution of the flow equation must be performed for each of the principal flow directions (i.e., three separate single-phase simulations need to be performed). Each simulation involves the iterative solution of a linear equation system (typically, the linear solver is a conjugate gradient method, preconditioned by incomplete Cholesky or LU factorization). The unknowns in this equation system are the geo-cell pressures inside the flow cell, whereas known quantities are the geo-cell dimensions and permeabilities, as well as the pressure conditions along the faces of the flow cell. When the numerical solution of the fine-scale pressure distribution has converged, an effective permeability is computed by equating expressions for the flux through the heterogeneous geo cells with the flux through the equivalent homogeneous flow cell using some form of Darcy’s law.

The pressure field is usually solved locally—that is, for one flow cell at a time. However, the size of the computational region may not necessarily be limited to that of the upscaling region (i.e., the flow cell).

Upscaling schemes for absolute permeability

Based on the size of the computational region, the single-phase upscaling process may either be described as local, regional, or global. With local upscaling techniques, the computational region is identical to the upscaling region (i.e., only geo cells inside the flow cell are considered in the upscaling computations). For regional upscaling, the computational region is expanded beyond that of the flow cell to include a buffer region of neighboring geo cells. In the case of global upscaling, the computational region is that of the entire geo model. Fig. 2 provides a schematic drawing of how the computational region varies with the different upscaling schemes.

• 

  Fig. 2 – Upscaling schemes and the size of the computational region (geo grid in white, flow grid in black).

It should be noted that the different upscaling schemes are only relevant when considering flow-based (tensor) methods. It is also important to realize that even though the computational region may vary according to the scheme used, the upscaling region remains unchanged and is of course defined by the flow cell, as in the case of the simple, analytical upscaling techniques.

Local upscaling

Because it used to be too time-consuming to compute the fine-scale pressure field for the complete geo grid in a single operation, the flow-based methods have traditionally been restricted to solving the pressure field locally—that is, for a single flow cell at a time. Hence, the effective cell permeability is computed separately and independently of the other flow cells, which may or may not be correct depending on how representative the imposed pressure conditions along the faces of the flow cell are.

Different types of artificial boundary conditions for the flow cell have been suggested over the years, all with the objective of providing as good an approximation of the real boundary conditions as possible. An important design criterion for the artificial boundary conditions is the conservation of flux in and out of the flow cell.

The first type of boundary conditions proposed for the local solution of the pressure equation was published by Warren and Price in 1961. Their approach is to impose a constant pressure gradient in a selected direction of flow by specifying a pressure of 1 on the inflow face and a pressure of 0 on the outflow face. By allowing no flow to pass through the sides of the cell, all fluxes are forced to go in the principal direction of flow. Therefore, this type of boundary conditions is often referred to as the no-flow or sealed-sides boundary conditions. The sealed-sides boundary conditions are graphically illustrated in Fig. 3 for flow in the vertical direction (here in the case of a flow cell containing a barrier).

• 

  Fig. 3 – The sealed sides boundary conditions (in z-direction).

The choice of boundary conditions emulates the way core permeability is measured in the lab. This is hardly a coincidence. As in the coreflood experiment, the local numerical flow simulation is in effect 1D because the cell faces parallel to the main flow direction are sealed. This implies that the estimated effective permeability will be scalar. Hence, the maximum number of directional permeability components that can be obtained with this type of boundary conditions is three, one for each of the principal directions of flow. In practice, the diagonal permeability tensor is derived by setting up the boundary conditions for x, y, and z directions, respectively, in three independent single-phase simulations.

As documented in two different sources, a tensor technique based on the sealed-sides boundary conditions tends to bias the estimated effective permeability toward a low value. The physical implication of this is most clearly seen in the case of a bimodal permeability system of sand and shale. This is because the sealed-sides method consistently underestimates the reservoir flow characteristics by thickening shale barriers and narrowing sand channels. The latter effect also has a tendency of disconnecting stacked sand channels.

Take, for example, the flow illustrated by Fig. 3. Because the barrier extends across the entire length of the local upscaling region, the resulting effective permeability (in the z-direction) will be zero. For vertical flow, the result, therefore, is a thickening of the shale in the flow model equal to the thickness of the flow cell. Depending on which factors that affect fluid flow in the region of the cell, this may or may not be a representative value for that particular flow cell.

Strictly speaking, the sealed-sides boundary conditions are only valid if no wells are present and the flow cells are symmetric in each direction of the grid as illustrated in 2D by Fig. 4. Hence, the sealed-sides boundary conditions assume that the flow cell is surrounded by mirror images of itself.

• 

  Fig. 4 – The fine-scale permeabilities are symmetric on the scale of the flow cell (geo grid in white, flow grid in black).

By the end of the 1980s, 3D geological models had started to appear more regularly on the modeling scene. This resulted in a new demand for advanced upscaling. In this renewed effort, two alternative boundary conditions for solving the local pressure solution in a flow-based method were suggested more or less at the same time. One was based on linear boundary conditions, the other on periodic boundary conditions.

The use of linear boundary conditions in flow-based upscaling was suggested by Guerillot et al. in 1989 and Samier in 1990 to enable the computation of a full-permeability tensor. Instead of setting the flow through the sides of the cell to zero, the pressure along the sides is allowed to vary in a linear fashion that matches the constant pressure on the two cell faces perpendicular to the flow. Hence, the imposed pressure gradient is still constant, but the flow is allowed to enter and leave the cell at any point along the sides parallel to the main flow direction. Therefore, this type of boundary conditions is also referred to as the open-sides boundary conditions. The situation is graphically illustrated in Fig. 5 for flow in the vertical direction (here in the case of a flow cell containing a barrier).

• 

  Fig. 5 – The open sides boundary conditions (in z-direction).

As with the sealed-sides boundary conditions, three independent single-phase simulations, with the main flow direction in x, y, and z, respectively, are needed to yield all of the components of the permeability tensor. With open-sides boundary conditions, however, also the off-diagonal components will generally be nonzero. Hence, unlike the sealed-sides boundary conditions where the effective permeability is limited to that of a diagonal tensor, the open-sides boundary conditions, as previously mentioned, give a full permeability tensor. The resulting full tensor may be either symmetric or nonsymmetric depending on the properties of the method under consideration.

As documented in two different sources, a tensor technique based on the open-sides boundary conditions tends to bias the estimated effective permeability toward a high value. The physical implication of this is most clearly seen in the case of a bimodal permeability system of sand and shale. This is because the open-sides method consistently overestimates the reservoir flow characteristics by narrowing shale barriers and thickening sand channels. The latter effect also has a tendency of connecting isolated sand channels.

Take, for example, the situation illustrated by Fig. 5. Even though the barrier extends across the entire length of the local upscaling region, the resulting effective permeability (in the z-direction) will be significantly larger than zero. For vertical flow, the result is therefore a narrowing of the shale in the flow model equal to the horizontal dimensions of the flow cell. Depending on which factors affect fluid flow in the region of the cell, this may or may not be a representative value for that particular flow cell.

The use of linear boundary conditions has its origin in the effective medium theory, which states that any region of permeability behaves as if embedded within the average medium. Strictly speaking, these boundary conditions are therefore only valid if the neighboring flow cells are of a uniform, nonzero permeability. This is illustrated in 2D by Fig. 6.

• 

  Fig. 6 – The fine-scale properties are uniform and nonzero outside the upscaling region (the geo cells outside the upscaling region are drawn in gray to reflect that the average value need not be known).

The use of periodic boundary conditions originates from the volume averaging theory, and its use in flow-based upscaling was first introduced by Durlofsky and Chung in 1990 and by Durlofsky in 1991. Durlofsky used periodic boundary conditions, together with Darcy’s law and the classic requirement of flux conservation, to derive a full permeability tensor. A somewhat different approach, which also uses a periodic pressure field around the flow cell, was proposed by Øistein Bøe et al. in 1994. This uses a weak form of Darcy’s law to prove that periodic boundary conditions result in a full permeability tensor that is both symmetric and positive definite. The Norsk Hydro tensor method is based on the conservation of dissipation (mechanical energy per unit weight of fluid), although it turns out that fluxes are conserved as well.

A simplistic illustration of the periodic boundary conditions is given in Fig. 7 for flow in the vertical direction (here in the case of a flow cell containing a barrier).

• 

  Fig. 7 – The periodic boundary conditions in z-direction (P_1k = P_nx,k, P_i1 = P_i,nz).

Although the periodic boundary conditions generally result in an effective permeability that is higher than that computed with the sealed-sides boundary conditions, the effective vertical permeability for the upscaling region illustrated in Fig. 7 will also be zero.

Strictly speaking, the periodic boundary conditions are only valid if no wells are present and the fine-scale medium is periodic on the scale of the flow cells (i.e., the fine-scale property distribution inside each flow cell must be identical). This is illustrated in 2D by Fig. 8. Please note that if a medium is symmetric on the scale of ΔL, then it will be periodic on the scale of 2ΔL.

• 

  Fig. 8 – The fine-scale properties are periodic on the scale of the flow cell (geo grid in white, flow grid in black).

The relative performance of the tensor methods that is caused by the various boundary conditions has proven to be of considerable interest. As it happens, the sealed-sides method provides a lower bound and the open-sides method an upper bound of the effective permeability. The periodic-based method turns out to give an effective permeability estimate that generally lies in between the two previous methods.

With regard to the outer bounds of effective permeability, it is well known that the harmonic and arithmetic means provide the absolute lower and upper limit of the effective permeability, respectively. It is less known that the uncertainty range in the effective permeability may be narrowed using the composite averages. In fact, it may be mathematically proven that the harmonic-arithmetic average provides a closer lower limit than the pure harmonic mean, whereas the arithmetic-harmonic average provides a closer upper limit than the pure arithmetic mean (truly valid only for regular grids). In this context, it is important to realize that the two flow-based methods (sealed and open sides) provide an even narrower uncertainty band for the effective permeability, but at the expense of increased CPU time.

The relative performance of the most important local upscaling techniques is shown in Fig. 9. Fig. 9 also indicates the "inner" uncertainty range of the true effective permeability for the sake of comparison.

• 

  Fig. 9 – The relative performance of local upscaling techniques.

Regional upscaling

Regional upscaling is applied to reduce the influence of the artificial boundary conditions on the effective permeability estimate by moving the boundary of the computational region away from the flow cell. This implies that the influence of neighboring geo cells is taken into account in addition to the geo cells inside the flow cell. In other words, regional upscaling represents an expansion of the local computational region outside the volume of the flow cell. The size of the so-called buffer or skin around each flow cell is usually given in number of neighboring geo cells to either side of the flow cell and must be specified by the modeler for each of the three coordinate directions.

The permeability estimate of a regional upscaling method will improve as the size of the buffer region increases, and it will ultimately be equal to the "true" effective permeability when the buffer size has reached the boundaries of the geo model for all three directions. The gain in accuracy is largest in the beginning (i.e., for small buffer values). This is illustrated by Fig. 10, showing the behavior of the lower- and upper-bound tensor techniques in the case of increasing buffer size. Please note that in Fig. 10, the outer bounds are shown to be symmetric around the "true" effective permeability. Generally, this is not the case.

• 

  Fig. 10 – The outer bounds of effective permeability and their behavior with increasing buffer size.

Global upscaling

Strictly, the fine-scale pressure field must be determined for the entire geo grid simultaneously to compute "exact" effective permeabilities for the flow cells. In the past, however, this has been too CPU-intensive to be performed in practice. With the introduction of new and promising solution algorithms such as the Output Least Squares (OSL) method, global upscaling schemes can now be realized. In the paper by Holden and Nielsen, the OSL method is used to minimize the difference in pressure, as well as velocity, between the geo and flow grids in an iterative process. Because the CPU consumption of the applied equation solver is proportional to the number of geo cells, a global solution will use approximately the same amount of computational time as the sum of all the local computations. Therefore, the new global upscaling scheme is just as fast as any local method.

An obvious advantage with the global upscaling approach is that one avoids using artificial boundary conditions around the upscaling region (i.e., instead of guessing what the boundary conditions for the flow cells might be, the pressure conditions surrounding the cells are explicitly known). Another important benefit is that a poor separation of scales in the upscaling will no longer occur because the size of the computational region is the same as the geo model.

Although still in its research stage, global upscaling has much potential for improving today’s permeability estimators, especially for models containing a complex facies architecture with large permeability contrasts between facies. In fact, according to Holden and Nielsen, preliminary results show an improvement factor of 10 in some cases.

Still, as discussed by Holden and Nielsen, the global upscaling approach is not enough to ensure maximum accuracy in the modeling of the effective permeability. Because the value of the effective permeability is influenced by changes in the pressure field, the flow-cell permeabilities should strictly be recomputed by the global method for every timestep taken by the multiphase flow simulator. In practice, though, it might be good enough to update the effective permeability field whenever a significant change occurs as a result of altering the well configuration or production/injection rates and so on. Hence, the ultimate upscaling scheme for the absolute permeability might be the one that is coupled with the multiphase flow simulator and automatically updates the absolute effective permeability field for each timestep. With the current computer power and the lack of proper integration between the geological model and the simulation model, this is hardly achievable yet.

Best practice guidelines

As may be understood from the previous sections, there exists no single upscaling method for absolute permeability that is superior to all other methods in all situations, at least not until it has been fully established that the global upscaling scheme represents the ultimate method of choice. Selecting the proper upscaling method from the many available choices can be quite a challenge. The choice of sophistication in the upscaling method generally depends on one or several of the following factors:

• Complexity of the fine-scale permeability distribution (i.e., the geo model)
• Degree of upscaling that needs to be performed (i.e., the coarsening factor)
• Number of permeability realizations that need to be upscaled
• Time available to the project for performing upscaling
• Intended use of the flow model

Because an exact validation of the upscaling process cannot be performed unless a multiphase flow simulation is carried out on the geo model itself, two alternative upscaling approaches for identifying the proper homogenization method are presented here.

Absolute upscaling approach

This approach assumes that there exists a way to properly validate the absolute performance of an upscaling method without resorting to an extremely time-consuming (if at all feasible) finite-volume simulation of the geo model. As documented in Samier et al., streamline simulation offers a very efficient way of validating the performance of upscaling methods. The validation process is carried out by first running a streamline simulation on the geo model to compute the reference solution. Then a streamline simulation is run on the flow model for each of the upscaling methods that are to be evaluated. The simulated performance of the various upscaling methods is then compared to that of the geo model. The validation of upscaling methods is best done under single-phase flow conditions to avoid introducing other model parameters (e.g., relative permeabilities and associated rock types) that also need to be upscaled in one way or the other. A higher confidence may also be obtained for the validation process if the actual well pattern is used in the streamline simulations.

Using the previously described validation scheme, the modeler may choose to evaluate any upscaling method until one with a satisfactory performance is found. Still, a more systematic way of identifying the optimum upscaling method is desirable. With the absolute upscaling approach that is presented here, the modeler is offered a multistep procedure that is to be terminated as soon as a satisfactory upscaling method has been identified. The recommended procedure involves the following steps:

1. Compute the upper and lower bounds of the effective permeability using the arithmetic-harmonic (or pure arithmetic) and harmonic-arithmetic (or pure harmonic) average techniques, respectively. Being of the analytical type, these methods are very fast and will provide a first quantification of the upscaling uncertainty. Validate the performance of the two composite methods against that of the geo model using a single-phase streamline simulator with the actual well pattern.
2. If the performance of any of the two methods in Step 1 is within an acceptable range of the geo model, then terminate the procedure and choose the appropriate method. If, on the other hand, the performance of both methods is unsatisfactory because of the complexity of the geo model, then use the upper (open-sides) and lower (sealed-sides) bound diagonal tensor methods to narrow the uncertainty in the flow-model performance. Validate the performance of the two tensor methods against that of the reference solution.
3. If the performance of any of the methods in Step 2 is within an acceptable range of the reference solution, then terminate the procedure and choose the appropriate method. If, on the other hand, the performance of both methods is unsatisfactory, then the following alternatives may be worth considering:

a. If time allows, refine or coarsen the flow grid (whatever is best) to achieve a better separation of the length scales. Then repeat the upscaling of the outer bounds (in Step 2) and redo the validation to check if the performance of either method has improved.

b. Apply a tensor method with periodic or semi-open-boundary conditions [the semi-open boundary conditions alternative is available in some applications using a multiplier between 0 (sealed) and 1 (open) to the side faces of the computational region]. As previously mentioned, this should result in an intermediate estimate of the effective permeability tensor and therefore provide a flow model performance that lies somewhere in between the two methods in Step 2.

c. Select the best of the two tensor methods in Step 2 and convert the local method to a regional method using a buffer region of modest size. Validate its performance. If necessary, repeat this step using an increasingly larger buffer region until a satisfactory performance of the flow model is obtained.

d. If for some reason none of the previous alternatives are an option, then one needs to apply the method that best satisfies the wanted flow behavior of a given cell or cells in a given region. In other words, a combination of the outer bound techniques within one and the same model may be a fourth alternative.

As previously mentioned, the method using open-sides boundary conditions is a good estimator of sand continuity and quality, whereas the sealed-sides boundary conditions method is better at detecting the presence and effect of barriers.

Consider a long horizontal oil producer in the Troll West Gas Province that is protected against coning from the overlying gas cap by a calcite barrier just above the well. If the vertical grid resolution in the Troll full-field model was such that one could apply the open-sides boundary conditions technique on the cells containing the well, and the sealed-sides boundary conditions technique on the cells containing the barrier, then this would be the optimum local upscaling approach.

However, if both a segment of the well trajectory and a segment of the calcite are present inside the same flow cell, then the open-sides boundary conditions technique will give a good estimate of the well’s pressure indicator (PI) but result in a much too early gas breakthrough, whereas the sealed-sides boundary conditions method will better capture the effect of the calcite but give a too low estimate of the well’s PI. If this is the case, one needs to consider applying one of the alternatives A or C (alternative B will, in this particular example, give the same result as the sealed-sides method).

The multistep procedure of the absolute upscaling approach is graphically illustrated in Fig. 11.

• 

  Fig. 11 – Flow chart for the absolute upscaling approach.

Relative upscaling approach

This approach acknowledges the fact that an exact validation of the upscaling results cannot be achieved in practice. Therefore, instead of trying to validate the absolute performance of an upscaling method, the approach diagnoses the relative performance of outer bound methods using the actual multiphase finite-volume simulator. This implies that a full black-oil simulation is run on the flow model for each of the upscaling methods that are to be evaluated. The deviation in the simulated performance between outer bound methods will then reflect the part of the model uncertainty that originates from the upscaling process itself. To ensure a high degree of relevance in the diagnostics, it is important that the test simulations contain a representative description of the actual flow model.

If the project is pressed for time, the simulations may be skipped altogether in favor of a faster, although less robust, way of performing the diagnostics. Instead of analyzing simulation profiles, a normalized difference parameter may be computed on a cell by cell basis using the formula

....................(1)

The relative upscaling approach that is presented here utilizes a multistep procedure that applies outer bound methods of increasing accuracy until the best possible upscaling method can be identified. The recommended procedure involves the following steps:

1. Compute the upper and lower bounds of the effective permeability using the arithmetic-harmonic (or pure arithmetic) and harmonic-arithmetic (or pure harmonic) average techniques, respectively. Being of the analytical type, these methods are very fast and will provide a first quantification of the upscaling uncertainty. Run the finite-volume simulator for each of the two composite methods (or compute a grid-based difference parameter) and perform the diagnostics.
2. If the performance gap between the two methods in Step 1 is acceptable (small), then terminate the procedure and choose either of the two. If the deviation in performance is unsatisfactory because of the complexity of the geo model, then use the upper (open-sides) and lower (sealed-sides) bound diagonal tensor methods to narrow the uncertainty in the flow model performance. Run the finite-volume simulator for each of the two tensor methods (or compute a grid-based difference parameter) and perform the diagnostics.
3. If this reduces the upscaling uncertainty to within acceptable limits, then either of the two diagonal tensor methods may be used to provide the final permeability field for the flow model. If the deviation in the performance is unsatisfactory (large), then the following alternatives may be worth considering:
   1. If time allows, refine or coarsen the flow grid (whatever is the best) to achieve a better separation of the length scales. Then repeat the upscaling of the outer bounds (in Step 2) and rerun the simulations (or recompute the difference parameter) to check if the performance gap (upscaling uncertainty) has narrowed.
   2. Apply a tensor method with semi-open or periodic boundary conditions. As previously mentioned, this should result in an intermediate estimate of the effective permeability tensor and hence provide a flow model performance that lies somewhere in between the two methods in Step 2.
   3. Convert the two local tensor methods in Step 2 into regional methods using a buffer region of modest size. Rerun the finite-volume simulator (or recompute the difference parameter) and check the performance gap. If necessary, repeat this step using an increasingly larger buffer region until the upscaling uncertainty reaches acceptable limits, at least as far as practically possible.
   4. If for some reason none of the previous alternatives are an option, then one needs to apply the method that best satisfies the wanted flow behavior of a given cell or cells in a given region. In other words, a combination of the outer bound techniques within one and the same model may be a fourth alternative.

The multistep procedure of the relative upscaling approach is graphically illustrated in Fig. 12.

• 

  Fig. 12 – Flow chart for the relative upscaling approach.

Nomenclature

∂

=

normalized difference parameter

Noteworthy papers in OnePetro

Use this section to list papers in OnePetro that a reader who wants to learn more should definitely read

External links

Use this section to provide links to relevant material on websites other than PetroWiki and OnePetro

See also

Reservoir simulation

Gridding in reservoir simulation

Streamline simulation

Reservoir simulation applications

Geostatistical reservoir modeling

PEH:Reservoir_Simulation