Transmissibility Upscaling on Unstructured Grids for Highly Heterogeneous Reservoirs

One critical point of modeling of flow in porous media is the capacity to consider parameters that are highly variable in space. It is then very challenging to simulate numerically fluid flow on such heterogeneous porous media. The continuous increase in computing power makes it possible to integrate smaller and smaller heterogeneities into geological models of up to tens of millions of cells. On such meshes, despite computer performance, multi-phase flow equations cannot be solved in an acceptable time for hydrogeologists and reservoir engineers, especially when the modeling considers several components in each fluid and when taking into account rock-fluid interactions. Taking average reservoir properties is a common approach to reducing mesh size. During the last decades, many authors studied the upscaling topic. Two different ways have been investigated to upscale the absolute permeability: (1) an average of the permeability for each cell, which is then used for standard transmissibility calculation, or (2) computing directly the upscaled transmissibility values using the high-resolution permeability values. This paper is related to the second approach. The proposed method uses the half-block approach and combines the finite volume principles with algebraic methods to provide an upper and a lower bound of the upscaled transmissibility values. An application on an extracted map of the SPE10 model shows that this approach is more accurate and faster than the classical transmissibility upscaling method based on flow simulation. This approach keeps the contrast of transmissibility values observed at the high-resolution geological scale and improves the accuracy of field-scale flow simulation for highly heterogeneous reservoirs. Moreover, the upper and lower bounds delivered by the algebraic method allow checking the quality of the upscaling and the gridding.


Introduction
One of the essential points for the hydrogeologists and the reservoir engineers is the transition from the geological mesh to the reservoir mesh. The properties of the reservoir mesh, being composed of cells supposed to represent an average of the properties of the geological model, affect the results of the simulations essential to solving different problems such as history matching, the optimization of the production, the production forecasting and their associated uncertainties [1]. This upscaling step is even more critical for geological reservoirs with high-heterogeneous properties, particularly for the management of aquifers (e.g., seawater intrusion in coastal aquifer [2]), the management of hydrocarbon reservoirs [3], or the geological sequestration of CO 2 [4,5].
Usually, geostatistical approaches are used to model the geological model and stochastically generate the distribution of the reservoir heterogeneity on meshes of tens of millions of cells. However, the numerical fluid flow simulations with commercial reservoir simulators are not conceivable with such mesh size. Despite the increase in computational speed, an averaging step of the heterogeneous properties is mandatory to reduce the mesh size and make it possible to perform multi-phase flow reservoir simulations. Many authors have investigated several averaging methods of reservoir properties in the last decades [6][7][8][9]. Regarding the averaging of absolute permeability, these methods can be split into two categories: (1) Approach #1: first, the absolute permeability is averaged for each reservoir cell, and then the discretization of the flow equations lead computing the transmissibility values [10][11][12][13][14]; (2) Approach #2: the average transmissibility values are calculated directly from the high-resolution scale permeability distributions [15][16][17][18][19][20].
Regardless of the approach chosen, the upscaling step leads to a loss of information. Different criteria have been proposed to quantify the error introduced by the upscaling step [21][22][23][24][25].
Generally, when the main flow direction is not parallel to the main axes of the grid, the average permeability should be represented by a full tensor for adequately taking into account the characteristics of the flow. But, the permeability obtained cannot be used directly in the discretized equation. A harmonic average between the permeability of two neighboring cells is then usually used to obtain the transmissibility values. Samier [15] proposes a finite element volume technique with schemes varying from 9 to 27 points, applied on non-uniform grid spacing. The framework developed by Ding and Urgelli [16] and Urgelli [17], deliver a technic to calculate the large-scale reservoir transmissibility values on a shifted block approach with a finite volume discretization of flow equation for Cartesian grids and corner point geometry (CPG) grids. In addition to the numerical scheme, the choice of the boundary condition differs according to the methods. The shifted block that has been proved in one dimension in space (1D) by Urgelli [17] to be more accurate than the use of full blocks (Approach #1).
Panfili et al. [19] have used upscaled transmissibility to a huge carbonate reservoir, the Kashagan field, located below the North Caspian Sea, and they have found this approach more accurate by numerical comparisons with simulations on high-resolution geological cells.
This paper follows the second approach. The half-cell volumes on both sides of the boundaries are taken into account to preserve the contrast of the heterogeneity better. Compared to previous works, the shifted block technique is considered, and the finite volume principles are combined with algebraic approaches to provide the two bounds (upper and lower) of the upscaled transmissibility values. The presented method uses a combined approach, including the advantages of algebraic methods and the flow-based method, i.e., the speed and accuracy. Moreover, the two bounds of the upscaled transmissibility values allow controlling the quality of upscaling and the gridding.
In the next section, the proposed approach is detailed, and the implementation steps are given. Then, to show the benefits of this method, a comparison of transmissibility estimation that we propose is done with the results obtains with the "classical" approach (Approach #1), on an extracted map of the SPE10 model.

Local Transmissibility Upscaling Method
The proposed transmissibility upscaling approach extends the algebraic upscaling method initially introduced by Guérillot et al. [11], called Cardwell and Parsons. We suggest a transmissibility upscaling version of this approach. This methodology, primarily based on the permeability upscaling techniques, enables the spatial distribution of local transmissibility values to be taken into account. The precision of the estimation is improved by considering the neighboring volumes on each side of the edge of the cells in a way leading to better preserves the transmissibility contrast introduced by the high-resolution geological cell-scale heterogeneities.
One of the benefits of the algebraic approach is that the upper and lower bounds are given: T + and T − . The accuracy of the upscaling and the quality of the gridding process can be controlled using these boundaries. In addition, if T + /T − , the ratio of T + versus T − is lower than a specified given threshold, then the averaged transmissibility values proposed by Guérillot et al. [11] would be the geometric mean of these two bounds. This geometric mean degenerates well in the case where the arithmetic average or the harmonic average is the correct solution. A ratio of the bounds above this given threshold indicates that the griding must be redesigned.
Combining the half-block approach and the finite volume principles with algebraic methods, the introduced "local" method consists of: (1) Remeshing each reservoir cell with a high-resolution geological mesh finer than the geological mesh; (2) Decomposing each reservoir cell in local sub-domain containing geological cells; (3) Calculate an upper and lower admissible half-transmissibility values on each sub-domain of each reservoir cell, by defining for each shifted block, 1D local problems with linear pressure boundary conditions [11] and apply a finite-volume numerical scheme to discretize these local problems; (4) Calculate upscaled transmissibility values between each reservoir cells using half-transmissibility values of corresponding cells weighted by the surface exchange between these cells.
Let us assume that the geological mesh refers to the ensemble of three-dimensions in space (3D) cells with their associated porosities and permeability tensor used to represent the geological model, that the simulator mesh can be built independently of the geological mesh, without constraint between both, and that the geological and reservoir mesh could be both built using cells represented by eight points in space without restrictions on their alignments. Throughout this paper, the lowercase refers to variables associated with the geological mesh and the capital letters to the reservoir mesh.

Remeshing
Each target cell (large reservoir cell), is refined with more cells than the source cells (high-resolution geological cells). The refining in I, J and K direction is done with an even cell number (cf. Figure 1).
For each source cells (high-resolution geological cells), sampling new high-resolution geological cells (large reservoir cell refined) is whose centers are contained in the source cell.
Water 2020, 12, x FOR PEER REVIEW 3 of 11 geometric mean of these two bounds. This geometric mean degenerates well in the case where the arithmetic average or the harmonic average is the correct solution. A ratio of the bounds above this given threshold indicates that the griding must be redesigned.
Combining the half-block approach and the finite volume principles with algebraic methods, the introduced "local" method consists of: (1) Remeshing each reservoir cell with a high-resolution geological mesh finer than the geological mesh; (2) Decomposing each reservoir cell in local sub-domain containing geological cells; (3) Calculate an upper and lower admissible half-transmissibility values on each sub-domain of each reservoir cell, by defining for each shifted block, 1D local problems with linear pressure boundary conditions [11] and apply a finite-volume numerical scheme to discretize these local problems; (4) Calculate upscaled transmissibility values between each reservoir cells using halftransmissibility values of corresponding cells weighted by the surface exchange between these cells.
Let us assume that the geological mesh refers to the ensemble of three-dimensions in space (3D) cells with their associated porosities and permeability tensor used to represent the geological model, that the simulator mesh can be built independently of the geological mesh, without constraint between both, and that the geological and reservoir mesh could be both built using cells represented by eight points in space without restrictions on their alignments. Throughout this paper, the lowercase refers to variables associated with the geological mesh and the capital letters to the reservoir mesh.

Remeshing
Each target cell (large reservoir cell), is refined with more cells than the source cells (highresolution geological cells). The refining in I, J and K direction is done with an even cell number (cf. Figure 1).
For each source cells (high-resolution geological cells), sampling new high-resolution geological cells (large reservoir cell refined) is whose centers are contained in the source cell.

Sub-Domains Decomposition
The half parts of the blocks on each side of the boundaries of the two grid blocks are extracted to compute the upscaled transmissibility values. Six spatial domains are defined for each cell. Right and Left domains allow computing the I-Transmissibility values. North and South allow computing the J-Transmissibility values. Top and bottom domains allow computing the K-Transmissibility values. Figure 2 presents a 2D illustration of the domains. Then following a spatial clustering of permeability, a one-dimensional local flow-based transmissibility calculation is done to define the lower and upper bound transmissibility values. Below, details are given to explain how these bounds are obtained, and a final estimator for the transmissibility values associated with reservoir mesh is proposed.

Sub-Domains Decomposition
The half parts of the blocks on each side of the boundaries of the two grid blocks are extracted to compute the upscaled transmissibility values. Six spatial domains are defined for each cell. Right and Left domains allow computing the I-Transmissibility values. North and South allow computing the J-Transmissibility values. Top and bottom domains allow computing the K-Transmissibility values. Figure 2 presents a 2D illustration of the domains. Then following a spatial clustering of permeability, a one-dimensional local flow-based transmissibility calculation is done to define the lower and upper bound transmissibility values. Below, details are given to explain how these bounds are obtained, and a final estimator for the transmissibility values associated with reservoir mesh is proposed.

Lower and Upper Bounds Transmissibility
The transmissibility between two cells is given as follows: where is the half-transmissibility for the neighboring cells and ( ) , defined as follows: where is the permeability of the cell in the direction , a normal vector to the interface, the surface between the cells A and B, the vector between the two centers: the cell center and the interface center (cf. Figure 3).

Implementation of the Lower Bound Transmissibility
For each local domain, with n high-resolution geological cells at the cell edge, the computation of the lower estimated bound value for the transmissibility is done considering n 1D steady-state incompressible flow.
The finite volume discretization of Darcy's law for these one-dimensional systems can be written as follows: Dirichlet boundary conditions are defined to solve this system with a pressure equal to unity on one side and zero on the other side. After discretization in space, the following matrix system needs to be solved. It is particularly easy as it is a tri-diagonal symmetric matrix: where is the matrix inclosing the transmissibility terms, the unknowns represented by a vector, and includes the boundary conditions. The dimension of the matrix is ( is the number of cells in each 1D system (cf. Figure 4).

Lower and Upper Bounds Transmissibility
The transmissibility between two cells is given as follows: where τ is the half-transmissibility for the neighboring cells i and n(i) , defined as follows: where k α is the permeability of the cell in the direction α , n α a normal vector to the interface, s AB the surface between the cells A and B, r α the vector between the two centers: the cell center and the interface center (cf. Figure 3).

Lower and Upper Bounds Transmissibility
The transmissibility between two cells is given as follows: where is the half-transmissibility for the neighboring cells and ( ) , defined as follows: where is the permeability of the cell in the direction , a normal vector to the interface, the surface between the cells A and B, the vector between the two centers: the cell center and the interface center (cf. Figure 3).

Implementation of the Lower Bound Transmissibility
For each local domain, with n high-resolution geological cells at the cell edge, the computation of the lower estimated bound value for the transmissibility is done considering n 1D steady-state incompressible flow.
The finite volume discretization of Darcy's law for these one-dimensional systems can be written as follows: Dirichlet boundary conditions are defined to solve this system with a pressure equal to unity on one side and zero on the other side. After discretization in space, the following matrix system needs to be solved. It is particularly easy as it is a tri-diagonal symmetric matrix: where is the matrix inclosing the transmissibility terms, the unknowns represented by a vector, and includes the boundary conditions. The dimension of the matrix is ( is the number of cells in each 1D system (cf. Figure 4).

Implementation of the Lower Bound Transmissibility
For each local domain, with n high-resolution geological cells at the cell edge, the computation of the lower estimated bound value for the transmissibility is done considering n 1D steady-state incompressible flow.
The finite volume discretization of Darcy's law for these one-dimensional systems can be written as follows: Dirichlet boundary conditions are defined to solve this system with a pressure equal to unity on one side and zero on the other side. After discretization in space, the following matrix system needs to be solved. It is particularly easy as it is a tri-diagonal symmetric matrix: where A is the matrix inclosing the transmissibility terms, u the unknowns represented by a vector, and b includes the boundary conditions. The dimension of the matrix A is (n ct + 2) × (n ct + 2) , n ct is the number of cells in each 1D system (cf. Figure 4).
The half-transmissibility can be written as follows: where is the index of a reservoir large cell, is the domain: right, left, north, south, top or bottom, is the number of high-resolution geological cells at the interface between the two large reservoir cells, and is the high-resolution scale flux. The lower estimated bound value for the transmissibility between two given large reservoir cells A and B (cf. Figure 5) is defined as follows: where is the cell interface area between the cells A and B, is the total surface area between the cells A and the connected cells on the domain, ̅ is the opposite domain of the domain .  The Darcy's flux at the interface between the two cells gives: The half-transmissibility can be written as follows: where β is the index of a reservoir large cell, d is the domain: right, left, north, south, top or bottom, n is the number of high-resolution geological cells at the interface between the two large reservoir cells, and q n is the high-resolution scale flux. The lower estimated bound value for the transmissibility between two given large reservoir cells A and B (cf. Figure 5) is defined as follows: where s AB is the cell interface area between the cells A and B, s A d is the total surface area between the cells A and the connected cells on the d domain, d is the opposite domain of the domain d .
The half-transmissibility can be written as follows: where is the index of a reservoir large cell, is the domain: right, left, north, south, top or bottom, is the number of high-resolution geological cells at the interface between the two large reservoir cells, and is the high-resolution scale flux. The lower estimated bound value for the transmissibility between two given large reservoir cells A and B (cf. Figure 5) is defined as follows: where is the cell interface area between the cells A and B, is the total surface area between the cells A and the connected cells on the domain, ̅ is the opposite domain of the domain .

Implementation of the Upper Bound Transmissibility
To compute the upper bound transmissibility value, the cells are clustered to generate slices perpendicular to the considered direction for each local domain (cf. Figure 6).
The computation of the half-transmissibility value (Equation (2)) between two slices is done as follow: • k α is the arithmetic mean of the permeability for each slice perpendicular to the considered direction. The upper bound of the admissible permeability value and the exact value for a stratified medium considering a parallel flow is given by arithmetic mean of the permeability values; • n α is a normalized vector through each cell; • s is the sum of the areas of the interfaces of the high-resolution geological cells; • r α is the arithmetic mean of the vector connecting the centers of each high-resolution geological cell with the center of the interface.
Water 2020, 12, x FOR PEER REVIEW 6 of 11

Implementation of the Upper Bound Transmissibility
To compute the upper bound transmissibility value, the cells are clustered to generate slices perpendicular to the considered direction for each local domain (cf. Figure 6).
The computation of the half-transmissibility value (Equation (2)) between two slices is done as follow: • is the arithmetic mean of the permeability for each slice perpendicular to the considered direction. The upper bound of the admissible permeability value and the exact value for a stratified medium considering a parallel flow is given by arithmetic mean of the permeability values; • is a normalized vector through each cell; • is the sum of the areas of the interfaces of the high-resolution geological cells; • is the arithmetic mean of the vector connecting the centers of each high-resolution geological cell with the center of the interface. Equation (1) allows calculating the transmissibility value across each slice. Then, the same strategy as the lower bound is applied, i.e., the computation of the upper estimated bound value for the transmissibility is done considering one steady-state incompressible 1D single-phase flow model without gravity. The pressure equation is solved considering Dirichlet boundary conditions with a pressure equal to unity on one side and zero on the other side.
The Darcy's flux at the interface between the two cells can be written as follows: The half-transmissibility can be written as follows: where is the index of the large reservoir cell, is the domain: right, left, north, south, top or bottom, and is the flux at the interface.
The upper estimated bound value for the transmissibility between two given large reservoir cells A and B (cf. Figure 7) can be defined as follows: where is the cell interface area between cell A and B, is the total surface area between cell A and the connected cells on the domain. Equation (1) allows calculating the transmissibility value across each slice. Then, the same strategy as the lower bound is applied, i.e., the computation of the upper estimated bound value for the transmissibility is done considering one steady-state incompressible 1D single-phase flow model without gravity. The pressure equation is solved considering Dirichlet boundary conditions with a pressure equal to unity on one side and zero on the other side.
The Darcy's flux at the interface between the two cells can be written as follows: The half-transmissibility can be written as follows: where β is the index of the large reservoir cell, d is the domain: right, left, north, south, top or bottom, and q is the flux at the interface. The upper estimated bound value for the transmissibility between two given large reservoir cells A and B (cf. Figure 7) can be defined as follows: where s AB is the cell interface area between cell A and B, s A d is the total surface area between cell A and the connected cells on the d domain.
Water 2020, 12, x FOR PEER REVIEW 7 of 11

Quality Control of the Upscaling and Gridding
The quality of the upscaling and the gridding can be evaluated considering the ratio of the upper and lower bounds. For each large reservoir cell where the ratio of the upper and the lower bounds is lesser than a given value , the upscaled transmissibility value is calculated considering the geometric mean of these bounds. A ratio of these bounds higher than a given value (larger than 10) highlights that the gridding must be revisited.
• For cells with > : reconsidering of gridding.

Application to the SPE-10 Model
The presented method is applied to a two-dimensional synthetic case corresponding to an extracted layer of the "SPE-10 2nd comparative solution" project [26]. This map corresponding to the 51st layer is located in the lower part (Upper-Ness formation) of the model. Figure 8 shows the permeability distribution in I and J-direction and the porosity distribution. The range of the permeability varies by 8 orders of magnitude (1.5 × 10 −3 mD to 2 × 10 5 mD). The range of the porosity varies from 10 −3 to 0.4.
One injector well and one producer well are considered. The location of the wells is given in Table 1. Both wells are vertical with perforation throughout the formation of the reservoir.  Injector  3  3  Producer  58  218 The scenario is an eight-year CO2 injection. The fluid injection rate is equal to 100 m 3 /day constraints by a maximum Bottom hole pressure at 690 bar. The fluid is produced at a constant bottom hole pressure of 275 bar.

Wells I-Coordinate J-Coordinate
The breakthrough curves obtained with the high-resolution geological scale and large-scale reservoir simulation have been compared. Figure 9 shows that the results of geological highresolution scale simulation and large-scale reservoir simulation with the local transmissibility upscaling method are similar. The arriving time is the same. The large-scale reservoir simulation with a classical transmissibility upscaling approach based on flow calculation shows an offset on the breakthrough time.

Quality Control of the Upscaling and Gridding
The quality of the upscaling and the gridding can be evaluated considering the ratio of the upper and lower bounds. For each large reservoir cell where the ratio of the upper and the lower bounds is lesser than a given value ε, the upscaled transmissibility value is calculated considering the geometric mean of these bounds. A ratio of these bounds higher than a given value (larger than 10) highlights that the gridding must be revisited.

•
For cells with •

Application to the SPE-10 Model
The presented method is applied to a two-dimensional synthetic case corresponding to an extracted layer of the "SPE-10 2nd comparative solution" project [26]. This map corresponding to the 51st layer is located in the lower part (Upper-Ness formation) of the model. Figure 8 shows the permeability distribution in I and J-direction and the porosity distribution. The range of the permeability varies by 8 orders of magnitude (1.5 × 10 −3 mD to 2 × 10 5 mD). The range of the porosity varies from 10 −3 to 0.4.
One injector well and one producer well are considered. The location of the wells is given in Table 1. Both wells are vertical with perforation throughout the formation of the reservoir. The scenario is an eight-year CO 2 injection. The fluid injection rate is equal to 100 m 3 /day constraints by a maximum Bottom hole pressure at 690 bar. The fluid is produced at a constant bottom hole pressure of 275 bar.
The breakthrough curves obtained with the high-resolution geological scale and large-scale reservoir simulation have been compared. Figure 9 shows that the results of geological high-resolution scale simulation and large-scale reservoir simulation with the local transmissibility upscaling method are similar. The arriving time is the same. The large-scale reservoir simulation with a classical transmissibility upscaling approach based on flow calculation shows an offset on the breakthrough time.

Conclusions
Reservoir engineering procedures require the upscaling of geological grid properties to perform reservoir simulations. The upscaling process affects these properties, especially in highly heterogeneous reservoirs.
In this paper, a new approach has been developed to upscale transmissibility values. This method uses the half-block approach and combines the finite volume principles with algebraic methods to provide an upper and a lower bound of the upscaled transmissibility values. The upscaled transmissibility calculation is done considering the half parts of the blocks on each side of the boundaries of the two grid blocks to keep the heterogeneity contrast observed at the highresolution geological scale. The combination of 1D numerical solutions and algebraic averaging leads to a rapid approach.
In this paper, we have shown on an extracted map of the SPE10 model that: (1) This approach is more accurate than classical transmissibility upscaling method based on flow simulation;

Conclusions
Reservoir engineering procedures require the upscaling of geological grid properties to perform reservoir simulations. The upscaling process affects these properties, especially in highly heterogeneous reservoirs.
In this paper, a new approach has been developed to upscale transmissibility values. This method uses the half-block approach and combines the finite volume principles with algebraic methods to provide an upper and a lower bound of the upscaled transmissibility values. The upscaled transmissibility calculation is done considering the half parts of the blocks on each side of the boundaries of the two grid blocks to keep the heterogeneity contrast observed at the highresolution geological scale. The combination of 1D numerical solutions and algebraic averaging leads to a rapid approach.
In this paper, we have shown on an extracted map of the SPE10 model that: (1) This approach is more accurate than classical transmissibility upscaling method based on flow simulation; Figure 9. Comparisons of breakthrough curves for the production well.

Conclusions
Reservoir engineering procedures require the upscaling of geological grid properties to perform reservoir simulations. The upscaling process affects these properties, especially in highly heterogeneous reservoirs.
In this paper, a new approach has been developed to upscale transmissibility values. This method uses the half-block approach and combines the finite volume principles with algebraic methods to provide an upper and a lower bound of the upscaled transmissibility values. The upscaled transmissibility calculation is done considering the half parts of the blocks on each side of the boundaries of the two grid blocks to keep the heterogeneity contrast observed at the high-resolution geological scale. The combination of 1D numerical solutions and algebraic averaging leads to a rapid approach.
In this paper, we have shown on an extracted map of the SPE10 model that: (1) This approach is more accurate than classical transmissibility upscaling method based on flow simulation; (2) This approach keeps the contrast of transmissibility values observed at the high-resolution geological scale and improves the accuracy of field-scale flow simulation for highly heterogeneous reservoirs; (3) The upper and lower bounds delivered by the algebraic method allow checking the quality of the upscaling and the gridding.
In this paper, the problematic of the transmissibility upscaling in fractured reservoirs [27][28][29][30] has not been addressed. An extension of the proposed method will be proposed to upscale the "matrix-fracture" transmissibility values for the dual-continuum models.
Author Contributions: D.G. and J.B. proposed the algorithm, worked together on the numerical implementation in Petrel and wrote the manuscript.
Funding: The publication of this article was funded by the Qatar National Library.