Theory

The coarsening approach implemented in pycopm is based on industry standard grid formats which is preserved during coarsening. This allows for direct application in standard reservoir simulators. In coarsening accuracy is traded for speed. This is often appropriate, but care must be taken to assure that the coarse solution preserve the key properties of the model. The Lie 2019 textbook gives an excellent introduction to grid coarsening and upscaling geological properties, and the methods implemented in pycopm follow similar techniques.

Grid coarsening

Corner-point grids are commonly use in subsurface simulations, due to its flexibility to handle complex geometries such as faults. The grid is defined by vertical pillars and horizontal lines connecting the pillars, resulting in cells that can have up to six faces (hexahedrons), but also these grids allow to define connections between non-neighboring cells (NNC). A special case of corner-point grids are Cartesian regular grids. These grids also feature inactive cells, which are cells with zero pore volume and no connected to the active cells.

We define the grid coarsening operator \(\Xi_{xyz}^C\), which creates from an input grid \(\Omega\) a coarser grid \(\Omega^*\) (we adopt the notation * for coarser quantities, while \({}^C\) in the coarsening operator emphasize the different approaches for the cell clustering) by removing vertical pillars and horizontal lines:

\[\Xi_{xyz}^C(\Omega)\rightarrow\Omega^*\]

To give flexibility in the coarsening, the cell clustering can be given as x, y, and z arrays to define the pillars/lines to be removed. One natural question is how to handle the inactive cells in a cluster, and for this, three options to define the coarser cells are implemented: min, max, and mode (mode is the default, i.e., the coarser cells is active if the number of active cells is the most common value in the cluster). Figure 1 shows a simple 2D corner-point grid with different coarsening using min, max, and mode.

Figure 1: Example of coarser models from a grid with 18 cells, where the cell #11 is inactive. One example where the max option could be useful is for models where there are a lot of inactive cells in the z direction, while the min option could be useful for models applying coarsening in the xy plane, since using min results in coarser models that do not generate new connections across inactive cells.

Upscaling geophysical properties

Naturally, the pore volume in a coarser cell \(\Phi_{i^*,j^*,k^*}^*\) (\(i^*\), \(j^*\), and \(k^*\) referring to the cell coarse indices in the x, y, and z direction respectively) are equal to the sum of pore volume from the corresponding cells \(\Phi_{i,j,k}\) in the input model, which are part of the cluster \(\mathbb{C}_{i^*,j^*,k^*}\):

\[\Phi_{i^*,j^*,k^*}^* = \sum_{(i,j,k)\in\mathbb{C}_{i^*,j^*,k^*}}\Phi_{i,j,k}\]

For example, in Figure 1 when coarsening into two cells using mode/max resulted in two coarse cells, where:

\[\mathbb{C}_{1^*,1^*,1^*}=\{ (1,1,1),(2,1,1),(3,1,1),(1,1,2),(2,1,2),(3,1,2),(1,1,3),(2,1,3),(3,1,3) \}.\]

From this definition, the porosity in the coarse model \(\phi_{i^*,j^*,k^*}^*\) can be simply computed by:

\[\phi_{i^*,j^*,k^*}^* = \frac{\Phi_{i^*,j^*,k^*}^*}{\mathbb{V}_{i^*,j^*,k^*}^*}\]

where \(\mathbb{V}_{i^*,j^*,k^*}^*\) is the geometric volume of the coarser cell.

For the rock permeability, there are different upscaling methods (e.g., arithmetic or harmonic average) that are case dependent and perform different, see the Lie 2019 textbook for comparison of these methods. In pycopm, by default the permeability in the x and y directions are computed using the arithmetic average, while the permeability in the z direction by the harmonic average. As additional options, the permeabilities in the coarser cells \(\mathbb{K}_{i^*,j^*,k^*}^*\) can be set to equal the min, max, and mean values from the permeabilities in the corresponding cluster \(\mathbb{K}_{i,j,k}\). For example, using the max for permeabilities could be useful for history matching studies, where the parameters to history match are saturation functions (relative permeabilities and capillary pressure).

The above line mentions the initial application to develop pycopm (coarsening to history match saturation functions), as such there are no upscaling methods implemented in pycopm for saturation functions. In a geological model, it is common to define different regions (referred as satnum) to assign different saturation function tables. Then, if a cluster \(\mathbb{C}_{i^*,j^*,k^*}\) involves different values for satnum, the mode (the most frequent value) is used to assign the value in the coarser cell (this is also used to assign additional discrete coarser values such as fluid-in-place regions (fipnum)).

For grids with large number of non-neighboring connections (faults) and inactive cells, then a better approach is to upscale transmisibilities. A drawback of upscaling transmissibilities is that permeabilities cannot be used in history matching, but instead, transmissibilities multipliers, which increases the number of parameters to history match and might break history match workflows where different permeability fields are generated from spatial correlations. To this end, two approaches to upscale transmissibilities are implemented in pycopm. The former computes the coarser transmissibility

\[\scriptsize\mathbb{T}_{i^*,j^*,k^*}^*=\left[ T_{i^*,j^*,k^*\rightarrow i^*+1,j^*,k^*}^*,T_{i^*-1,j^*,k^*\rightarrow i^*,j^*,k^*}^*,T_{i^*,j^*,k^*\rightarrow i^*,j^*+1,k^*}^*,T_{i^*,j^*-1,k^*\rightarrow i^*,j^*,k^*}^*,T_{i^*,j^*,k^*\rightarrow i^*,j^*,k^*+1}^*,T_{i^*,j^*,k^*-1\rightarrow i^*,j^*,k^*}^* \right]\]

using the armonic averaging along the transsmissibility direction and summing over these values over the cell coarser face. For example, for the transmissibility in the z+ direction:

\[T_{i^*,j^*,k^*\rightarrow i^*,j^*,k^*+1}^*=\sum_{(i,j)\in\mathbb{C}_{i^*,j^*,k^*}}\left( \sum_{(i,j,k)\rightarrow(i,j,k+1)\in\mathbb{C}_{i^*,j^*,k^*}} \frac{1}{T_{i,j,k\rightarrow i,j,k+1}} \right)^{-1}\]

For cases where the coarsening is only along one direction (e.g., the z direction), the second method sets the transmissibility on the coarse cell faces in that direction equal to the overlapping cell face values in the corresponding cluster (instead of computing the armonic average). For input models with a large number of inactive cells, this approach has resulted in better results with respect to the input model simulations than using the armonic average. For both approaches, the transmissibilities are scaled with the ratio of the cell face effective areas (input model) to the coarse cell area. For non-neighbouring connections, this approach is also implemented, i.e., the non-neighbouring connections in the coarser model sum the values from the non-neighbouring connections in the input model, which is important in order to honor the pressure connections along open faults communicating different formations.