Selecting an Appropriate Upscaled Reservoir Model Based on Connectivity Analysis
Sélection de méthodes d’upscaling par analyse de connectivité
IFP Energies nouvelles, 14 avenue de BoisPréau, 92852
RueilMalmaison Cedex  France
email: christophe.preux@ifpen.fr  mickaele.leravalec@ifpen.fr  guillaume.enchery@ifpen.fr
^{*} Corresponding author
Received:
2
December
2014
Accepted:
21
July
2016
Reservoir engineers aim to build reservoir models to investigate fluid flows within hydrocarbon reservoirs. These models consist of threedimensional grids populated by petrophysical properties. In this paper, we focus on permeability that is known to significantly influence fluid flow. Reservoir models usually encompass a very large number of fine grid blocks to better represent heterogeneities. However, performing fluid flow simulations for such fine models is extensively CPUtime consuming. A common practice consists in converting the fine models into coarse models with less grid blocks: this is the upscaling process. Many upscaling methods have been proposed in the literature that all lead to distinct coarse models. The problem is how to choose the appropriate upscaling method. Various criteria have been established to evaluate the information loss due to upscaling, but none of them investigate connectivity. In this paper, we propose to first perform a connectivity analysis for the fine and candidate coarse models. This makes it possible to identify shortest paths connecting wells. Then, we introduce two indicators to quantify the length and trajectory mismatch between the paths for the fine and the coarse models. The upscaling technique to be recommended is the one that provides the coarse model for which the shortest paths are the closest to the shortest paths determined for the fine model, both in terms of length and trajectory. Last, the potential of this methodology is investigated from two test cases. We show that the two indicators help select suitable upscaling techniques as long as gravity is not a prominent factor that drives fluid flows.
Résumé
Les ingénieurs de réservoir construisent des modèles de réservoir pour comprendre les écoulements de fluides dans les réservoirs d’hydrocarbures. Ces modèles se composent de grilles en trois dimensions peuplées par des propriétés pétrophysiques. Dans cet article, nous nous concentrons sur la perméabilité qui est connue pour influencer de manière significative l’écoulement du fluide. Les modèles de réservoir ont généralement un très grand nombre de mailles fines afin de mieux représenter les hétérogénéités. Cependant, les simulations d’écoulement sur ces modèles fins sont très coûteuses en temps CPU. Une pratique courante consiste à convertir les modèles fins en modèles grossiers avec moins de mailles : c’est le processus de mise à l’échelle. De nombreuses méthodes de changement d’échelle ont été proposées dans la littérature qui mènent à des modèles grossiers différents. Le problème est de savoir comment choisir la méthode de mise à l’échelle appropriée. Différents critères ont été établis pour évaluer la perte d’information due à l’upscaling, mais aucun d’entre eux ne s’intéresse à la connectivité. Dans cet article, nous nous proposons d’abord d’effectuer une analyse de connectivité pour les modèles fins et grossiers. Cela nous permet d’identifier les chemins les plus courts reliant les puits. Puis, nous présentons deux indicateurs pour quantifier la différence de longueur et de trajectoire entre les chemins du modèle fin et des modèles grossiers. La méthode de mise à l’échelle qui sera recommandée est alors celle qui fournit le modèle grossier pour lequel le plus court chemin est le plus proche de celui déterminé pour le modèle fin, à la fois en termes de longueur et de trajectoire. Enfin, le potentiel de cette méthodologie est étudié à partir de deux cas de test. Nous montrons que les deux indicateurs permettent de sélectionner les techniques de changement d’échelle appropriées tant que la gravité n’est pas un facteur important de l’écoulement.
© C. Preux et al., published by IFP Energies nouvelles, 2016
This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Introduction
Reservoir models are commonly built to represent hydrocarbon reservoirs. They consist of grids with petrophysical properties assigned to every cells. They can be given as inputs to flow simulators in order to predict fluid displacements and hydrocarbon production. A necessary condition to ensure the reliability of the simulated flow responses is that the reservoir models closely replicate the true spatial distribution of petrophysical properties. These models are usually finely gridded in order to properly capture heterogeneity occurring also at the fine scale. Therefore, they encompass a very significant number of grid blocks, which contribute to strongly increase the computational overburden. A possibility to circumvent this drawback is upscaling: the fine geological models are converted into coarser reservoir models. This results in a decrease in the number of grid blocks and makes flow simulation feasible in a reasonable amount of time. Clearly, the coarse reservoir models do not reproduce exactly the dynamic behavior of the fine geological models as upscaling induces a loss of information. Many upscaling techniques have been described in the literature and some of them may be more appropriate depending on the case studied. The interested reader can refer to [1] for a thorough review. However, selecting the most suitable technique is far from straightforward. To date, distinct criteria [26] have been proposed to quantify the information loss due to upscaling so that the quality of an upscaled model can be estimated without performing any fluid flow simulation. For instance, it was recommended to define the information loss from the ratio of the number of fine cells to the number of coarse cells or from the comparison of the histograms of the studied property at the fine and coarse scales. However, a few methods only have focused on connectivity, which is known to have a strong impact on fluid flow. A recent comprehensive review of the influence of connectivity on fluid flow is available in [7]. There are several works in the literature, which use streamlines for upscaling and its evaluation [810]. Even though faster than classical volume finite simulation, streamline simulation still calls for the solution of a transport equation. In this paper, we aim to develop an alternative method to streamline simulation based on reservoir connectivity analysis, which does not solve any transport equation. This method can also be used with any transmissibility considered. It is beyond the scope of this paper to perform any statistical analysis of the upscaled properties. However the reader interested in this topic can refer to [11]. We consider hereafter that the permeability field is known and far above the percolation threshold [12]. We introduce two new quality indicators based upon reservoir connectivity. First, we determine the shortest paths between injectors and producers. These shortest paths are assumed to approximate the main flow paths. Then, we focus on the lengths and trajectories of these paths. The idea behind this connectivity analysis is to identify the main flow paths without performing any full fluid flow simulation. Such tools can help identify the most appropriate upscaling technique because we aim to build coarse models that preserve the same level of connectivity between wells as the fine models. The potential of such criteria to evaluate the suitability of upscaled reservoir models is investigated from two numerical experiments inspired by the SPE10 case. We apply different upscaling methods to build coarse reservoir models from a given fine reservoir model. As the upscaling techniques are not equivalent, the resulting coarse models are different. Then, we compute the connectivitybased indicators for each of them. This helps select the most suitable upscaled models for the field case studied. Last, we verify the soundness of this approach by inputting the coarse reservoir models into a fluid flow simulator.
1 Definition of a ConnectivityBased Criterion
1.1 Connectivity Analysis
A reservoir model is defined as a grid with porosity and permeability values assigned to grid blocks. Following graph theory, the centers of the grid blocks are considered as nodes connected by weighted links. The value of the weight characterizing the link between two nodes provides the distance between these two nodes or the cost to move from one node to the other one. Reservoir models are usually based on simple Cartesian grids or Corner Point Grids (CPG). In this case, the maximal number of nodes connected to one node is known and small: for a twodimensional grid, it is at most 4 (Fig. 1). Within the context of reservoir modeling, the distance between two nodes i and j can be considered as the ratio of the pore volume to transmissibility T _{ i↔j } [13]. The suitability of such a definition is questionable as it involves only static properties. To incorporate a flavor of dynamic information without performing the costly full two or three phase flow simulation, we suggest a model using a simplified pressure field. This is obtained by solving the steady state single phase conservation equation. Therefore we just solve this equation and use the resulting pressures. The distance we use comes from the calculation of discrete times of flight. Indeed, one possible way to compute these times consists in solving numerically the following system of equations:(1)where v is Darcy velocity, Φ the porosity, λ _{t} the total mobility, k the permeability and P the pressure. The total mobility λ _{t} is equal to where Kr stands for the relative permeability and μ for the viscosity. As boundary conditions, we assume that the time of flight is equal to 0 in all perforated cells crossed by an injector well. System (1) can be solved by means of a classical finite volume scheme [14]. For each grid cell i, the first equation of (1) is discretized as follows:(2)where stands for the set of all neighbouring cells of i and for the set of the perforations of cell i, if any. Let us consider more precisely the fluxes between two cells. In Equation (2), τ _{ i↔j } represents a discrete upwind time of flight for the face i ↔ j which is given by:
Figure 1. Representation of a reservoir model as a weighted graph. 
In the right hand side, corresponds to the pore volume of the grid block i. For Cartesian grids, the flux v _{ i→j } oriented from cell i to cell j can be computed by using a two point flux approximation:where T _{ i↔j } is the transmissivity and the term is upwinded as previously according to the sign of (P _{ i } − P _{ j }). The discrete system (2) corresponds to a linear system that should be inverted in order to obtain the discrete times of flight τ _{ i } in all cells i. Instead of solving this linear system, we chose in this work to resort to a graph search algorithm solving the singlesource shortest path problem for a graph with nonnegative face path costs. The costs that should be defined on each face i ↔ j have been defined based on the following remarks. By dividing Equation (2) by we can observe that the term in factor to the time of flight τ _{ i↔j } is equal to:which corresponds to the inverse of the discrete time needed to move from cell i to cell j. Conversely,represents the inverse of the discrete time needed to move from cell j to cell i. Since two different pore volumes appear in both previous expressions, we can take a geometric mean of the two and leave out the mobility term which can only be known if at least a twophase model has been previously solved. This led us to define the following surface rate on each interface i ↔ j:(3) d _{ i→j } is a rate considered hereafter as a distance or a surface between two connected nodes i and j and moving from node i to node j (in m^{2}/s). Fluid flow being driven by pressure, d _{ i→j } is positive when P _{ i } > P _{ j }; otherwise, the link is removed (the fluid can flow only from high to low pressures). Once we know the distances between the connected nodes, we can determine the shortest paths between injectors and producers using Dijkstra’s algorithm [15]. Its complexity is in , but a variant of it can be used in the case studied. The algorithm developed by Fredman and Tarjan [16] is appropriate when the number of links per node is small. Its complexity reduces to . In such conditions, identifying the shortest paths is very fast. An example stressing why it is important to account for pressure, even approximately, is depicted in Figure 2. We consider a permeability field with 220 grid blocks along the X axis and 60 along the Y axis. For simplicity, we assume that the reservoir is initially saturated with oil and that production is based upon water injection. There are two wells: the injector in grid block (1, 1) and the producer in the opposite corner. Figure 2a, shows the shortest path derived from static information only [13] while Figure 2b shows the shortest path derived from Equation (3) with pressure information. For comparison purposes only, we performed the full fluid flow simulation. The resulting water saturation is reported at the bottom for the whole field at breakthrough time. The two shortest paths are visibly distinct. We verify that the shortest path that accounts for pressure is much more consistent with the main flow path derived from the full fluid flow simulation.
Figure 2. a) Permeability field and associated shortest path derived from connectivity analysis without accounting for pressure. b) Permeability field and associated shortest path derived from connectivity analysis with pressure accounted for. c) Water saturation at breakthrough time. 
1.2 ConnectivityBased Indicators
We now propose to refer to the shortest paths defined as explained above to assess the ability of a coarse reservoir model to reproduce the dynamic behavior of a fine reservoir model. In a nutshell, we consider that the connectivity properties of a “good” coarse model must be quite similar to the ones of the fine model. We focus hereafter on the determination of the shortest paths connecting injectors to producers. Keeping this goal in mind, we introduce two connectivitybased indicators to measure the closeness of the shortest paths determined for both the coarse and fine models. If we consider E the set of cells included in a shortest path connecting an injector to a producer in a grid G, the length of the shortest path is computed as follows:(4)
1.2.1 First Indicator: Length of the Shortest Path
We consider a path connecting an injector to a producer. The first indicator compares the length of the shortest path identified for the coarse model to the length of the shortest path for the fine model:(5)
This indicator tends to 0 when the lengths of the two shortest paths tend to the same value. As a result, an appropriate coarse model must be characterized by an I _{ L } indicator as small as possible. Our purpose is not to exactly reproduce the length of the shortest path. Differences between the shortest path for the fine and the coarse models are expected. Instead, we aim at ranking the upscaled models by considering the paths. The idea is to identify the upscaled model for which the shortest path is the closest to the one obtained for the fine model. Let us consider a simple two dimensional reservoir model with 2 layers only (Fig. 3a). This model is then upscaled so as to get a single layer (Fig. 3b). Theoretically, in such conditions, the arithmetic average is the most appropriate upscaling method. The volume and porosity are assumed to be identical for all grid cells. Thus, we obtain:(6)
Figure 3. Simple example a) fine grid with 2 layers, b) coarse grid with one layer, (c) pressure gradient (same for fine and coarse grids). 
Each well works at constant bottom hole pressure (Fig. 3c) and we have a linear pressure gradient. So, . To compare the paths, we can also consider:(7)
However . As Δx _{ F } = Δx _{ G } and , we have:(8)
For this simple case, we want to compare:(9)
The results for the shortest path length indicator are given in Table 1. They make it possible to check that the best upscaled model is the one derived from the arithmetic average, as expected (and the 2D isotropic Bound Combination (BC), which comes down to the arithmetic method in this simple test case).
Shortest path length indicator for the simple case
1.2.2 Second Indicator: Trajectory of the Shortest Path
The comparison of the lengths of the shortest paths is not enough as two paths may be different while having the same length. We therefore introduce a second indicator to compare also their trajectories. We note the coordinates of a cell k that belongs to the path identified at the fine scale. Likewise, is the coordinates of a cell K that belongs to the corresponding path at the coarse scale. Superscript i indicates that the path under consideration is the one connecting the injector i to producer p. The measure of difference between the fine and the coarse paths is given as (Fig. 4):(10)
Figure 4. Sketch illustrating the computation of the distance between the shortest path determined for the fine model (red/yellow cells) and the shortest path defined for the coarse model (blue cells). 
The measure M is the sum of the minimum distances calculated between each coarse cell and the cells on the fine grid path. That minimum distance is measured using the Euclidean distance between the centers of coarse and fine grid cells. Examples are represented in green in Figure 4. It has to be noticed that the measure between a coarse grid block and a fine one is set to zero when the fine grid block is included into the coarse one: in this case we get the resolution limit. These fine grid blocks are in yellow in Figure 4.
2 Upscaling Techniques
As already mentionned, we focus on the upscaling of the absolute permeability. In the more general case of multiphase flow, the upscaling of the relative permeabilities has also to be considered. However, in many cases, it is possible to build reasonably accurate coarse models with only the permeability (and porosity) upscaled. Renormalization is the most effective upscaling technique as far as connectivity through scales is considered as a concern [17, 18]. However, this technique is not common in commercial upscaling software. This is the reason why it is not studied here. In addition, in [18], the authors explain that this method gives results very close to the geometric mean which is one of the methods considered hereafter. This section recaps the basics of the upscaling techniques we applied to build coarse reservoir models. The idea is to refer to local or intrinsic techniques that do not require full or even simplified fluid flow simulations over the whole reservoir model. The “local” and “intrinsic” terms mean that the equivalent property attributed to a given coarse cell depends only on the properties on the fine cells included in the coarse one [19, 20]. We note V and v the volumes of the coarse and fine cells, K and k the permeabilities of the coarse and fine cells and J the set of fine cells included in the coarse cell:

Arithmetic mean. When referring to the arithmetic mean, the permeability of coarse cell i is derived from:

Harmonic mean. The harmonic mean is given by:
As is wellknown, the arithmetic and harmonic means provide upper and lower bounds for the coarse permeability:

Geometric mean. The geometric mean is given by:
μ stands for mean while subscripts a, h and g indicate whether the mean is arithmetic, harmonic or geometric.

2D isotropic BC method. Referring to the same idea according to which the coarse permeability can be derived from the averaging of some bounds, one possibility consists in defining bounds closer than the arithmetic and harmonic means and in averaging them. As the interval between the bounds gets narrower, the error in the coarse permeability is expected to decrease. Considering a uniform flow in a given direction, Cardwell and Parsons [21] showed that an upper bound K ^{1} for the coarse permeability is derived from the harmonic mean of the arithmetic means of the fine permeabilities, calculated over each slice of cells perpendicular to the flow direction:
A lower bound K ^{2} is estimated the same way considering slices parallel to the flow direction. K ^{2} is then given by the arithmetic mean of the harmonic means of the fine permeabilities:(15)
Then, Le Loch’ [22] approximated the coarse permeability using the geometric average of these two bounds:(16)

3D isotropic BC method. The above method was extended to the case of three dimensions by Lemouzy [23]. This entails the introduction of K ^{3} and K ^{4}:
In this case, the coarse permeability is given by:(19)

3D anisotropic BC method. The BC method was refined later on in order to account for anisotropy [24, 25], thus resulting in:
The θ exponents are used to describe anisotropy.

Numerical local method. The numerical local method is a flowbased method: the fine pressure equation is solved for every coarse grid blocks [2628]. It accounts for the spatial distribution of permeabilities at the fine scale and is usually considered as more accurate than the simple algebraic means presented above. However, a significant issue is the choice of the boundary conditions to be imposed [29]. Noflow boundary conditions are considered hereafter.
Other very interesting methods based on transmissibility upscaling exist, but they are not used in this paper [30].
3 Numerical Experiment
Two types of numerical experiments are performed to verify whether the proposed connectivitybased indicators can help identify the best coarse reservoir models.
3.1 First SPE10 Case
The case studied in this first numerical experiment is the one developed within the framework of the Tenth SPE Comparative Solution [31]. It is a twophase (oil and gas) reservoir with a simple 2D vertical crosssectional geometry and no dipping or faults (Fig. 4). The reservoir is 762 m long, 7.62 m wide and 15.24 m thick. The fine grid used for the fine geological model encompasses 100 × 1 × 20 grid blocks, all with identical dimensions. The top of the model is at 0.0 m and the initial pressure at this level is 6.89 bar (100 psia). Initially, the field is fully saturated with oil, meaning that there is no mobile water. For simplicity, porosity is assumed to be 20% over the entire reservoir. On the other hand, the spatial distribution of permeability is characterized by an anisotropic variogram with the main range along the horizontal axis. In addition, fluids are assumed to be incompressible and immiscible. Their viscosities and densities are reported in Table 2. The relative permeabilities are given in [31]. Last, we consider that capillary pressures are negligible. The reservoir is produced from two wells, a vertical injector on the left and a vertical producer on the right (Fig. 5). Both of them have an internal diameter of 0.3048 m (1.0 ft) and are perforated over the entire reservoir. The injection rate is set to 6.97 m^{3}/day with a bottom pressure limit of 47.55 bar (690 psia). The producer is submitted to a constant bottom pressure limit of 6.55 bar (95 psia). Two injection schemes are investigated in the following subsections: gas injection and water injection. The coarse reservoir grid consists of 25 × 1 × 5 grid blocks. In other words, a coarse grid block includes 4 × 4 fine grid blocks. The coarse grid is then populated with coarse permeability values derived from the 5 following upscaling techniques: arithmetic mean, harmonic mean, geometric mean, 2D isotropic BC method and numerical local method. In this specific case, the 2D and 3D isotropic BC methods are equivalent as the permeability field is a twodimensional section. The following step consists in identifying the shortest path between the injector and the producer for the fine and the 5 coarse models. It is worth mentioning that the distance between two grid blocks (Eq. 3) belonging to the same well is set to 0 when performing the connectivity analysis. This implies that there is a zero distance between two points of the same well. This yields a correct path between two perforated cells. At this stage, there is no need to perform any full fluid flow simulation, but a single pressure solve.
Figure 5. First SPE test case: permeability field. 
Fluid properties
3.1.1 First Indicator: Length of the Shortest Path
The shortest path determined for the coarse model built from the arithmetic mean is the one whose length is closer to the length of the shortest path at the fine scale. The values calculated for the length indicator are listed in Table 3.
Shortest path length indicator for the first Tenth SPE case
3.1.2 Second Indicator: Trajectory of the Shortest Path
The shortest paths determined for the 5 upscaling methods are displayed in Figure 6. They are compared to the shortest path defined for the fine model. The corresponding D distances (Eq. 10) are provided in Table 4. As can be seen in Figure 6, the coarse models yielding the shortest path the closest to the one obtained at the fine scale is derived from the 2D isotropic BC methods. In this case, the superimposition is perfect. Otherwise, there is more or less the same discrepancy when using the arithmetic mean or the numerical local method. This discrepancy gets stronger with the geometric and harmonic means. According to the first length indicator, the upscaling technique to be used in the case studied is the arithmetic mean. On the other hand, the trajectory indicator pleads in favor of the isotropic BC approximation. These results show that the two indicators must be jointly considered to find a good upscaling technique. In view of these indicators both together, the recommended upscaling techniques could be either the bound combination, the numerical local methods or the arithmetic mean. For sure, the harmonic and geometric means have to be rejected.
Figure 6. Comparison of the shortest path determined for the fine model (red cells) to the shortest paths determined for the 5 upscaled models (blue cells). The 5 upscaled models are built from the arithmetic mean, the geometric mean, the numerical local method, the harmonic mean, and, the isotropic BC method. 
Shortest path trajectory indicator for the first case of the Tenth SPE Comparative Solution
3.1.3 Validation
We now run the full fluid flow simulations for the fine and the 5 coarse models, the objective being to compare the simulated production responses. We focus first on gas injection, second on water injection.
Gas Injection
The simulated bottomhole pressure and gas injection rate at surface condition are displayed in Figure 7 for the injector. Clearly, the bottom hole pressure is drastically overestimated for the coarse model derived from the harmonic mean. In addition, the injection rate for this coarse model does not respect the required constraint of 6.97 m^{3}/day. A first conclusion, given the fine and coarse simulation results, is that the harmonic average is not appropriate in this first example for building the coarse reservoir model. This is actually consistent with the analysis based upon the shortest path length and trajectory indicators (Tab. 3 and 4). On the other hand, the bottomhole pressure computed for the coarse model derived from the geometric mean is also overvalued. This is again in agreement with the two length and trajectory indicators. Tables 3 and 4 visibly indicate that this upscaling technique is not suitable for the case studied. Therefore, we focus hereafter on the other upscaling techniques. The oil rate simulated for the producer is now shown in Figure 8. We observe that the responses simulated for the coarse models constructed from the bound combination and the numerical local methods are pretty close, but less good than the ones obtained for the arithmeticbased coarse model. Thus, considering the results reported in Figures 7 and 8, we conclude that the best upscaling method to be used here is the arithmetic mean. This yields a coarse model reproducing the dynamic behavior simulated for the fine model. Besides, the analysis of the shortest path length indicator (Tab. 3) points out that the most suitable upscaling technique is the arithmetic mean, closely followed by the numerical local and BC methods. Therefore, the flow simulation results seem to validate the information content provided by this length indicator. However, the analysis of Table 4 leads to different conclusions. The shortest path trajectory indicator recommends the use of the BC method that provides better results than the numerical local method and arithmetic mean. These apparently puzzling results can be explained by the nature of the fluid injected. Gas being lighter than oil, it tends to move upward in the reservoir. However, we do not account for gravity when computing the distance (Eq. 3) between two adjacent grid blocks. The shortest path identified for the fine model goes down: it is not capable of capturing the gravity effect. Therefore, it departs from the main flow path (Fig. 9). This is the same for the shortest paths computed for the coarse models. The reason why the arithmetic mean is more suitable in this specific case is just pure coincidence.
Figure 7. Gas rate a) and bottom hole pressure b) at the injector simulated for the fine first SPE tenth case model and the corresponding coarse models in the case of gas injection. 
Figure 8. Oil rate at surface conditions simulated for the fine first SPE tenth case model and the corresponding coarse models in the case of gas injection. 
Figure 9. Shortest path (gray) superposed to the gas saturation map at the fine scale. 
Water Injection
In order to verify our previous conclusions and get rid of the prominent effect of gravity, we now consider the same case, but with water injection (Tab. 2) instead of gas injection. The bottomhole pressure and water rate simulated for the injector are displayed in Figure 10. Again, the harmonic and geometric means provide coarse models that overestimate bottomhole pressure. It is better to avoid these upscaling techniques as suggested by the length and trajectory indicator values (Tab. 3 and 4). When focusing on the oil rate at the producer (Fig. 11), we verify that the studied upscaling methods, the harmonic and geometric means being disregarded, finally yield coarse models that behave more or less the same way. This is now consistent with the conclusions derived from the analysis of the length and trajectory indicators. Therefore, provided gravity is not prominent, the proposed length and trajectory indicators can be used to make recommendations about the upscaling techniques to apply or not.
Figure 10. Water rate a) and bottom hole pressure b) at the injector simulated for the fine first SPE tenth case model and the corresponding coarse models in the case of water injection. 
Figure 11. Oil rate at surface conditions simulated for the fine first SPE tenth case model and the corresponding coarse models in the case of water injection. 
3.2 Second SPE10 Case
A second numerical experiment is considered to investigate the interest of the proposed connectivity indicators. The case studied in this section is the second one developed within the framework of the Tenth SPE Comparative Solution [31]. It is characterized by a very simple geometry, with no top structure or faults. The fine geological model (Fig. 12) is constructed on a regular Cartesian grid with dimensions 1200 × 2200 × 170 ft^{3}. It describes part of a Brent sequence with the Tarbert formation on the top and the Upper Ness one at the bottom. The Tarbert formation represents a prograding nearshore environment: it is 70 ft thick and includes 35 layers. The bottom 100 ft, that is the Upper Ness formation, is fluvial and encompasses 50 layers. The dimension of a grid block is 20 × 10 × 2 ft^{3} at this fine scale. The fine grid consists of 60 × 220 × 85 cells, that is 1.122 million cells.
Figure 12. Second tenth SPE case: permeability field. 
3.2.1 Flow Model
The flow model is briefly described below. The interested reader can refer to [31] for more details. The formation volume factor, compressibility, viscosity and density for water are 1.01 bbl/bbl, 3 × 10^{−6} psi^{−1}, 0.3 cp and 1025 kg/m^{3} (64 lb/ft^{3}), respectively. The oil properties are characterized by the PVT table reported in Table 5. Oil density is assumed to be 849 kg/m^{3} (53 lb/ft^{3}). The coarse grid built from the fine grid contains 30 × 110 × 17 cells (56 100 cells). In such conditions, the modeling of the top formation involves 7 layers and the bottom one 10 layers. The upscaling methods used to populate the coarse grid with coarse permeability values are the arithmetic, harmonic and geometric means, the 2D and 3D isotropic bound combination, the 3D anisotropic bound combination and numerical local methods. Porosity being additive, coarse values are computed using the arithmetic mean. The reservoir is produced from 4 producers located in the 4 corners of the model (Tab. 6), each producing at 275.7 bar (4000 psi) bottomhole pressure. There is also a central water injection with a injection rate of 795.9 m^{3}/day (5000 bbl/day) and a maximum pressure of 689 bar (10 000 psi). All wells are vertical and completed throughout the entire formation.
PVT table for the oil phase
Well locations
3.2.2 Length and Trajectory Indicators
The singlesource shortest path problem was separately solved for the two units of the second tenth SPE case. This choice was preferred due to the very distinct geological features of these two units. Thus, we identified the 4 shortest paths connecting the injector to the producers for the fine model and the 7 coarse ones, first for the top unit and second for the bottom unit. We then calculated the shortest paths length and trajectory indicators for all coarse models. The length indicator is reported in Tables 7 and 8 for the top and bottom units, respectively. The trajectory indicator is also provided for the two units in Tables 9 and 10. Given the length indicator, we conclude that the coarse models derived from the harmonic and geometric means are not able to duplicate the lengths of the shortest paths determined for the fine model. As such, they can be considered as unsuitable. Otherwise, the arithmetic mean, the 2D and 3D isotropic bound combination, 3D anisotropic bound combination and numerical local method yield coarse models for which the lengths of the shortest paths are close to the ones obtained for the fine model. This conclusion holds whatever the unit considered. At this stage, we focus on the trajectory indicator. The harmonic and geometric means are no longer considered as disqualified by the length indicator. In the case of the top unit (Tab. 9), all results seem to be more or less equivalent. This can be explained by the fact that the coarse models in this case are all very smooth. They do not exhibit any special kind of contrast or connectivity. None of the candidate upscaling methods turns out to be better than the others. The analysis for the bottom unit is different (Tab. 10). The high values of the trajectory indicator calculated for the coarse arithmetic meanbased model suggest that the arithmetic mean has to be avoided. Finally, in this case, the recommendation would be to apply either the 2D or 3D isotropic BC method, the 3D anisotropic BC method or the numerical method.
Shortest path length indicator for the TOP unit
Shortest path length indicator for the BOTTOM unit
Shortest path trajectory indicator for the TOP unit
Shortest path trajectory indicator for the BOTTOM unit
3.2.3 Validation
Again, to validate the above analysis, we performed the full flow simulation for the fine model and the 7 corresponding coarse models. The bottomhole pressures and water injection rates simulated for the injector at surface condition are reported in Figure 13. As indicated by the shortest path study previously performed, the harmonic mean turns out to be undesired. It leads to an overestimated pressure and an injection rate significantly below the required constraint of 794.93 m^{3}/day (5000 bbl/day). The same criticism can be formulated when considering the oil rates: they are underestimated for the four producers (Fig. 14: the oscillations shown here result from the difficulty to converge. It could be removed by adjusting the time step. However, for comparison purposes, we decide to use the same time step for all simulations). This definitely excludes the harmonic mean. Another obvious result is that the bottomhole pressure simulated for the injector when considering the coarse model derived from the geometric mean is very different from the expected one (Fig. 13). Therefore, we also eliminate this upscaling technique. Let us now focus on the arithmetic mean. Figure 13 shows that the pressure is also underestimated in this case. In addition, it does not properly reproduce the oil rates computed for the fine model. This again suggests that the arithmetic mean has to be rejected. Finally, just as suggested by the analysis of the two length and trajectory indicators, the recommended upscaling techniques in this studied case are the 2D or 3D isotropic BC method, the 3D anisotropic BC method or the numerical method.
Figure 13. Bottomhole pressure b) and water rate a) simulated at the injector for the fine second SPE tenth case model and corresponding coarse models. 
Figure 14. Oil rates at surface conditions simulated at producers for the fine second SPE tenth case model and the corresponding coarse models. 
Conclusion
In this paper, we presented a novel methodology based upon connectivity analysis that can be seen as a firstpass to evaluate how appropriate an upscaling technique is for a given study case. It actually relies on the definition of two indicators that are easy and fast to estimate as they do not require any full fluid flow simulation. The leading idea is to determine the shortest paths between injectors and producers, these paths being considered as approximations of the main flow paths. The two indicators make it possible to compare the lengths and trajectories of the shortest paths defined for the fine model and the coarse model under consideration. Two numerical experiments were presented to investigate the potential of the proposed methodology. The first one was the first tenth SPE case. In the case of gas injection, we noted that the conclusions derived from the indicator study were not consistent with the results of the flow simulations. This was related to the fact that the definition of the shortest paths does not account for gravity. Therefore, the indicators can not properly capture the flow behavior. The same case were run with water injection instead of gas injection. We then verified the consistency between the flow simulation results and the study of the shortest path length and trajectory indicators. The second example considered was the second tenth SPE case. Again, it stressed the validity of the proposed method. The analysis of the two indicators made it possible to disqualify the arithmetic, harmonic and geometric means while it showed the suitability of the 2D and 3D isotropic BC method, the 3D anisotropic method and the numerical local method. An important future step to prove the soundness of this methodology will consist in studying more complex field cases with more wells. Moreover, the shortest paths as defined in this paper do not account for gravity. This is actually the subject of ongoing research, the idea being to be able to extend the definition of our connectivity criteria to the case of gas injection. It is also worthwhile emphasizing the interest in entending the above analysis to a stochastic context.
References
 Renard Ph., De Marsily G. (1997) Calculating equivalent permeability: A review, Adv. Water Resour. 20, 56, 253–278. [CrossRef] [Google Scholar]
 Sablok C., Aziz K. (2005) Upscaling and discretization errors in reservoir simulation, SPE Simulation Symposium, Houston, Texas, SPE 93372, 31 January2 February. [Google Scholar]
 Preux C. (2011) Study of evaluation criteria for reservoir upscaling, 73rd EAGE Conference & Exhibition Incorporating SPE EUROPEC, Vienna, Austria, 2326 May. [Google Scholar]
 Preux C. (2012) Coarsening criteria for reservoir upscaling, 74th EAGE Conference & Exhibition Incorporating SPE EUROPEC 2012, Copenhaguen, Denmark, 47 June. [Google Scholar]
 Preux C. (2014) About the use of quality indicators to reduce information loss when performing upscaling, Oil Gas Sci Technol. 71, 7. [CrossRef] [EDP Sciences] [Google Scholar]
 Qi D., Hesketh T. (2004) Quantitative evaluation of information loss in reservoir upscaling, Pet. Sci. Technol. 22, 11, 1625–1640. [CrossRef] [Google Scholar]
 Renard Ph., Allard D. (2013) Connectivity metrics for subsurface flow and transport, Adv. Water Resour. 51, 168–196. [CrossRef] [Google Scholar]
 Ates H., Bahar A., ElAbd S., Charfeddine M., Kelkar M., DattaGupta A. (2003) Ranking and upscaling of geostatistical reservoir models using streamline simulation: A field case study, Society of Petroleum Engineers Middle East Oil Show, Bahrain, 912 June, SPE81497MS. [Google Scholar]
 Maschio C., Schiozer D.J. (2003) A new upscaling technique based on dykstra parsons coefficient: Evaluation with streamline reservoir simulation, Journal of Petroleum Science and Engineering 40, 12, 27–36. [CrossRef] [Google Scholar]
 Samier P., Quettier L., Thiele M. (2002) Applications of streamline simulations to reservoir studies, SPE Reserv. Eval. Eng. 5, 4. [CrossRef] [Google Scholar]
 Boschan A., Noetinger B. (2012) Scale dependence of effective hydraulic conductivity distributions in 3D heterogeneous media: A numerical study, Transp. Porous Media 94, 1, 101–121. [CrossRef] [Google Scholar]
 Hunt A.G., Ewing R. (2009) Percolation Theory for Flow in Porous Media, Lecture Notes in Physics, 2nd edn., Springer. [Google Scholar]
 Li D., Wu X.H., Sun T., Goulding F.J., Stuart R.M., Chartrand T.A., Ramage C.J. (2010) Method for quantifying reservoir connectivity using travel times, US Patent, US 2010/0057418 A1. [Google Scholar]
 Eymard R., Gallouët T., Herbin R. (2000) Finite volume methods, in P.G. Ciarlet, J.L. Lions (eds), Handbook of numerical analysis, Vol. VII, NorthHolland, Amsterdam, pp. 7131020. [Google Scholar]
 Dijkstra E.W. (1959) A note on two problems in connection with graphs, Numer. Math. 1, 269271. [CrossRef] [MathSciNet] [Google Scholar]
 Fredman M.L., Tarjan R.E. (1984) Fibonacci heaps and their uses in improved network optimization algorithms, 25th Annual Symposium on Foundations of Computer Science (IEEE), pp. 338–346. [Google Scholar]
 King P.R. (1989) The use of renormalization for calculating effective perme ability, Transp. Porous Media 4, 1, 37–58. [Google Scholar]
 Gautier Y., Noetinger B. (1997) Preferential flowpaths detection for heterogeneous reservoirs using a new renormalization technique, Transp. Porous Media 26, 1–23. [CrossRef] [Google Scholar]
 SanchezVila X., Guadagnini A., Carrera J. (2006) Representative hydraulic conductivities in saturated groundwater flow, Rev. Geophys. 44, 3. [CrossRef] [Google Scholar]
 Wen X.H., GomezHernandez J.J. (1996) Upscaling hydraulic conductivities in heterogeneous media: An overview, J. Hydrol. 183, 9–32. [Google Scholar]
 Cardwell W.T., Parsons R.L. (1944) Average permeabilities of heterogeneous oil sands, Trans. Am. Inst. Mining Metal. Eng. 160, 34–42. [Google Scholar]
 Le Loch G. (1987) Étude de la composition des perméabilités par des méthodes variationnelles, PhD Thesis, École Nationale Supérieure des Mines de Paris. [Google Scholar]
 Lemouzy P. (1991) Calcul de la perméabilité absolue effective, Technical report, Institut Français du Pétrole, Paris. [Google Scholar]
 Duquerroix J.P.L., Lemouzy P., Noetinger B., Romeu R.K. (1993) Influence of the permeability anisotropy ratio on largescale properties of heterogeneous reservoirs, 68th Annual Technical Conference and Exhibition of the Society of Petroleum Engineers, Houston, Texas, SPE 26612, October. [Google Scholar]
 Romeu R. (1994) Écoulement en milieu poreux hétérogènes: prise de moyenne de perméabilités en régime permanent et transitoire, PhD Thesis, University of Paris VI. [Google Scholar]
 Christie M.A. (1996) Upscaling for reservoir simulation, J. Pet. Technol. 48, 11, 1004–1010. [CrossRef] [Google Scholar]
 Pickup G.E., Hern C.Y. (2002) The development of appropriate upscaling procedures, Transp. Porous Media 46, 119–138. [CrossRef] [Google Scholar]
 Pickup G.E., Ringrose P.S., Jensen J.L., Sorbie K.S. (1994) Permeability tensors for sedimentary structures, Math. Geol. 26, 2, 227–250. [CrossRef] [Google Scholar]
 Durlofsky L. (2005) Upscaling and gridding of fine scale geological models for flow simulation, 8th International Forum on Reservoir Simulation, Iles Borromees, Stresa, Italy, 2024 June. [Google Scholar]
 Guérillot D., Bruyelle J. (2014) A fast and accurate upscaling of transmissivities for field scale reservoir simulation, ECMOR XIV  14th European Conference on the Mathematics of Oil Recovery, 8 September. [Google Scholar]
 Christie M.A., Blunt M.J. (2001) Tenth SPE comparative solution project: A comparison of upscaling techniques, SPE Reservoir Simulation Symposium, Houston, Texas, SPE 66599, February. [Google Scholar]
Cite this article as: C. Preux, M. Le Ravalec and G. Enchéry (2016). Selecting an Appropriate Upscaled Reservoir Model Based on Connectivity Analysis, Oil Gas Sci. Technol 71, 60.
All Tables
Shortest path trajectory indicator for the first case of the Tenth SPE Comparative Solution
All Figures
Figure 1. Representation of a reservoir model as a weighted graph. 

In the text 
Figure 2. a) Permeability field and associated shortest path derived from connectivity analysis without accounting for pressure. b) Permeability field and associated shortest path derived from connectivity analysis with pressure accounted for. c) Water saturation at breakthrough time. 

In the text 
Figure 3. Simple example a) fine grid with 2 layers, b) coarse grid with one layer, (c) pressure gradient (same for fine and coarse grids). 

In the text 
Figure 4. Sketch illustrating the computation of the distance between the shortest path determined for the fine model (red/yellow cells) and the shortest path defined for the coarse model (blue cells). 

In the text 
Figure 5. First SPE test case: permeability field. 

In the text 
Figure 6. Comparison of the shortest path determined for the fine model (red cells) to the shortest paths determined for the 5 upscaled models (blue cells). The 5 upscaled models are built from the arithmetic mean, the geometric mean, the numerical local method, the harmonic mean, and, the isotropic BC method. 

In the text 
Figure 7. Gas rate a) and bottom hole pressure b) at the injector simulated for the fine first SPE tenth case model and the corresponding coarse models in the case of gas injection. 

In the text 
Figure 8. Oil rate at surface conditions simulated for the fine first SPE tenth case model and the corresponding coarse models in the case of gas injection. 

In the text 
Figure 9. Shortest path (gray) superposed to the gas saturation map at the fine scale. 

In the text 
Figure 10. Water rate a) and bottom hole pressure b) at the injector simulated for the fine first SPE tenth case model and the corresponding coarse models in the case of water injection. 

In the text 
Figure 11. Oil rate at surface conditions simulated for the fine first SPE tenth case model and the corresponding coarse models in the case of water injection. 

In the text 
Figure 12. Second tenth SPE case: permeability field. 

In the text 
Figure 13. Bottomhole pressure b) and water rate a) simulated at the injector for the fine second SPE tenth case model and corresponding coarse models. 

In the text 
Figure 14. Oil rates at surface conditions simulated at producers for the fine second SPE tenth case model and the corresponding coarse models. 

In the text 