1 Introduction

The spatial distribution of petrophysical properties in subsurface fluid reservoirs is controlled by the geologic processes that formed the reservoir and typically has a significant impact on fluid flow (Gómez-Hernández and Wen 1998; Miller et al. 1998; Jackson et al. 2003, 2005; Janković et al. 2006; Deveugle et al. 2011; Graham et al. 2015b) It is, therefore, important to capture geologic heterogeneity accurately and realistically in reservoir models. In most modelling workflows, the model volume is subdivided early in the workflow into a grid or mesh that contains a large number (typically hundreds-of-thousands to millions) of discrete volumes (grid blocks/cells/elements or similar, heretoforward termed cells) and geostatistical methods are used to populate these cells with petrophysical properties (Janssen et al. 2006; Feyen and Caers 2006). The equations governing fluid flow are then discretized and solved numerically (e.g. Helmig and Huber 1998; Jackson et al. 2015).

In some studies, the same grid is used to model the petrophysical properties and solve the flow equations; more typically, the flow equations are solved using a coarser grid and petrophysical properties upscaled from the finer grid. In either case, the grid used to simulate flow contains petrophysical properties that vary on a cell-to-cell basis (e.g. Fig. 1a).

Fig. 1
figure 1

a Reservoir models typically contain cell-to-cell variations in petrophysical properties; b The geological concept captures the architecture of key heterogeneities. Adapted from Christie and Blunt 2001

Full size image

There are several problems with this workflow. First, models are typically computationally expensive to create and maintain, as they contain a large number of cells. Second, upscaling from one grid to another can result in a loss of model fidelity (Renard and De Marsily 1997). Finally, the petrophysical properties vary from cell-to-cell at a length-scale that is controlled by the chosen grid geometry and resolution, rather than the underlying geological heterogeneity. While petrophysical properties of rocks do vary on a point-by-point basis [where a point here corresponds to a Representative Elementary Volume (e.g. Nordahl and Ringrose 2008)], this variability is at a different (and usually much finer) scale than the cell-to-cell variability captured in reservoir models. That petrophysical properties and flow equations are discretized on the same grid appears to be an accidental outcome of the gradual evolution of modelling workflows. Discretization of the petrophysical properties should be designed to preserve those features that are most important to flow; discretization of the flow equations should be designed to yield the most accurate solutions for the lowest computational cost. It is highly improbable that a single grid will meet both of these requirements.

An alternative to the grid-based modelling workflows described above is to use surface-based modelling (e.g. Xie et al. 2001; Pyrcz et al. 2005; Jackson et al. 2005, 2009, 2015; Sech et al. 2009; Zhang et al. 2009; Caumon et al. 2009; Graham et al. 2015a, b; Massart et al. 2016; Jacquemyn et al. 2018). In this approach, all geological heterogeneity that impacts flow is modelled as discrete volumes, or ‘geological domains’, bounded by surfaces. A nested hierarchy of surfaces can represent multi-scale heterogeneity. The surfaces are modelled without reference to an underlying grid; a grid or mesh is created only when required for numerical calculations. The grid cells inherit the petrophysical properties of the parent geological domain, which may be constant, or vary according to locally defined trends consistent with the underlying geologic concept: for example, to reflect upwards or downwards grading in grain size, or proximal to distal trends in reservoir quality. In both cases, the model can be constructed without reference to an underlying grid or mesh: petrophysical properties are assigned to each grid cell or mesh element only when the grid or mesh is created for numerical calculations.

There are several potential advantages to surface-based modelling: surfaces can accurately and realistically capture complex geological features at low computational cost; the model geometry can be created without reference to a grid or mesh defined early in the workflow, so the grid resolution does not place restrictions on model resolution, and the surfaces can be used to constrain the geometry of a grid or mesh that is optimally designed for the numerical solution of interest. However, in surface-based modelling, the short length-scale, cell-to-cell variability in petrophysical properties observed in grid-based models is omitted (e.g. Fig. 1b). It is not yet clear whether this cell-to-cell variability is required to properly capture flow. If cell-to-cell variability can be omitted from reservoir models while still correctly predicting fluid flow, then the route is open to focus modelling effort on capturing geologically defined, surface-bounded domains in models that are ‘fit for purpose’ (Ringrose and Bentley 2015; Bentley and Ringrose 2017).

No previous studies have specifically investigated the importance of cell-to-cell scale variability on predicted fluid flow, although some have tested the relative importance of modelling petrophysical property variations within, and between, geologically defined domains that correspond to lithofacies or lithofacies assemblages. Bianchi and Zheng (2016) found that a key control on contaminant flow was the connectivity of different lithofacies types, which were associated with significant contrasts in petrophysical properties. Similarly, Feyen and Caers (2006) found that lithofacies geometry (and the related permeability variation) was more important to solute migration than the distribution of permeability within lithofacies. However, in both studies, cell-to-cell variability in petrophysical properties within lithofacies was retained, so the hypothesis that cell-to-cell variability is not important to flow was not tested. White and co-workers (White and Barton 1999; Willis and White 2000; White et al. 2004) examined a number of petrophysical modelling approaches including homogenized lithofacies, assigning lithofacies-specific geometric mean values of permeability. Other approaches tested used geostatistical methods following Deutsch and Journel (1998) to populate petrophysical properties within lithofacies, yielding cell-to-cell variability. Comparison of the models using flow-based measures showed that if the length scale of petrophysical variability was small with respect to the length scale of the problem of interest, the impact on flow was small.

In this paper, we test the hypothesis that cell-to-cell variability can be omitted in reservoir models. It is not an upscaling study; rather, it is a study designed to test whether it is necessary to model cell-to-cell scale variability in petrophysical properties in the first place. If this variability is omitted, then it does not need to be upscaled. Nor is it a study designed to test whether small length-scale geologic heterogeneity is important to flow. Numerous previous studies have used nested, hierarchical models to test the impact on flow of heterogeneity at small- to large length-scales, and shown that both can be significant (e.g. Jones et al. 1995; Howell et al. 2008; Choi et al. 2011). However, the cell-to-cell variability in grid-based reservoir models does not capture small length-scale heterogeneity: petrophysical properties in real rocks, unlike grid-based models, are not homogeneous over close-to-rectangular volumes measuring of order 10s to 100s m laterally and 10s to 100s cm vertically. The cell-to-cell variability observed in grid-based reservoir models is an artefact of applying pixel-based geostatistical modelling methods on a grid of arbitrary resolution and geometry that is unrelated to the resolution and geometry of geologic heterogeneity present in the reservoir.

To test the hypothesis, we develop a workflow that removes cell-to-cell variability from conventional grid-based models, each of which initially contain several hundred thousand unique values of porosity and permeability, to yield models containing just a few tens of unique porosity and permeability values grouped into a few hundreds to a few thousands of domains. The same initial grid is used to solve the flow equations and the flow behavior of the original model is used as a reference. Note that we do not propose this workflow as a practical method to construct reservoir models; it is implemented here only to test the hypothesis in question.

2 Test Cases

Three reservoir-scale models were used as test cases. The first two models are synthetic and represent unfaulted reservoirs deposited in fluvial and shallow marine environments. The third is a real reservoir comprising fluvial, deltaic and shallow marine deposits divided into several different fault blocks. The synthetic models were taken from the SPE 10 Model 2, a classic benchmark model (see Christie and Blunt 2001), heretoforward referred to as the ‘SPE 10 Model’. This model has been widely used to represent a challenging problem for reservoir simulation with a highly heterogeneous distribution of petrophysical properties. All three models simulated multiphase flow where oil is displaced by water, although we suggest that the key findings of the study generalize to other types of flow relevant to aquifers and CO\(_{2}\) storage. The impact of heterogeneity on flow in the synthetic models was tested using two different injector/producer well configurations and three sets of oil/water fluid properties. The real reservoir case used the properties of the fluids in place and the well configuration applied in the field development.

Fig. 2
figure 2

Petrophysical properties within the Lower SPE 10 Model. Plots a to c show the spatial distribution of horizontal permeability in (a) 3D; b vertical cross section (X–\(X'\)) and c horizontal cross section (Y–\(Y'\)). Dimensions in metres. Plots d to f show the frequency distribution of d horizontal permeability; e vertical permeability and f porosity. Plots d to f also show the frequency distribution after division into three segments; see “Method” section for more details

Full size image
Fig. 3
figure 3

Petrophysical properties within the Upper SPE 10 Model. Plots a to c show the spatial distribution of horizontal permeability in (a) 3D; b vertical cross section (X–\(X'\)) and c horizontal cross section (Y–\(Y'\)). Dimensions in metres. Plots d to f show the frequency distribution of d horizontal permeability; e vertical permeability and f porosity

Full size image
Fig. 4
figure 4

The H-Field Model. Plot a shows the 3D spatial distribution of horizontal permeability; plot b shows the fault blocks, where red indicates the fault plane; plot c shows the location of production wells in green and water injection wells in blue. Dimensions in metres. Plots d to f show the frequency distribution of d horizontal permeability; e vertical permeability and f porosity

Full size image

2.1 Lower SPE 10: Fluvial Model

The first synthetic model comprised the lower 50 layers of the SPE 10 Model, which was originally constructed to represent reservoirs in the Upper Ness formation, part of the Middle Jurassic Upper Brent Group in the Northern North Sea (Mitchener et al. 1992). The Upper Ness formation was largely deposited in a non-marine, fluvial environment and reservoirs are typically characterised by high permeability channelized sandbodies surrounded by low permeability mudstone (Mitchener et al. 1992). The synthetic model with dimensions of \(365 \times 670 \times 30\) m is described on a regular \(60 \times 220 \times 50\) Cartesian grid (Fig. 2a to c). The model contains high permeability channelized sandbodies surrounded by low permeability mudstone and the modelling of these different lithofacies types yields bimodal porosity and permeability distributions (Fig. 2d to f). The model has a variable permeability anisotropy (\(K_{\mathrm{V}}/K_{\mathrm{H}}\)) ratio of 0.3 in the sandstone and \(10^{-3}\) in the mudstone.

2.2 Upper SPE 10: Shallow Marine Model

The second synthetic model comprised the upper 35 layers of the SPE 10 Model with dimensions of \(365 \times 670 \times 21\) m, described on a regular \(60 \times 220 \times 35\) Cartesian grid. This model was originally constructed to represent reservoirs in the Tarbert formation, which directly overlies the Upper Ness formation. The Tarbert formation was largely deposited in a shallow marine environment and reservoirs are typically characterized by an overall upwards coarsening and general distal to proximal coarsening succession (Mitchener et al. 1992). However, there is no strongly correlated heterogeneity present in the model and no clear evidence of upwards coarsening or distal to proximal coarsening trends (Fig. 3a to c). The porosity and horizontal permeability show uni-modal distributions; however, the vertical permeability exhibits a bimodal distribution. Although not clear from Christie and Blunt (2001), it appears that two different \(K_{\mathrm{V}}/K_{\mathrm{H}}\) ratios of 0.3 and \(10^{-3}\) were applied in the upper SPE 10 model (Fig. 3d to f).

2.3 H-Field: North Sea Reservoir

The real field example is located in the East Shetland Basin. The field comprises a number of N-S trending, westerly dipping, Mesozoic fault blocks. The reservoir is provided by fluvial, deltaic and shallow marine sandstones of the Brent Group, Lower Jurassic Statfjord Formation and Upper Triassic Upper Lunde Formation (Ritchie 2003). Bitumen horizons observed at specific levels suggest a staged filling of the reservoir, protecting up-dip regions from continued diagenesis (Baillie et al. 1996). The reservoir therefore exhibits significant permeability contrasts between layers due to this complex diagenetic phenomenon (Fig. 4a to c), with marked improvement in reservoir quality above the paleo-oil–water contact (Baillie et al. 1996). However, the petrophysical properties exhibit a uni-modal frequency distribution (Fig. 4d to f). The structure of the reservoir model is complex, comprising a number of fault-bounded blocks, described on a \(58 \times 145 \times 40\) corner point grid with stair-stepped faults (Fig. 4a to c).

2.4 Rock and Fluid Properties

The rock properties used in the SPE 10 models were the same as Christie and Blunt (2001). The relative permeability model is based on the Brooks–Corey correlation, defined by

$$\begin{aligned} K_{\mathrm{rw}}(S_{\mathrm{w}})&= \bigg (\frac{S_{\mathrm{w}}-S_{\mathrm{wirr}}}{1-S_{\mathrm{wirr}}-S_{\mathrm{or}}} \bigg )^{n_{\mathrm{w}}} \end{aligned}$$
(1a)
$$\begin{aligned} K_{{\mathrm{ro}}}(S_{\mathrm{o}})&= \bigg (\frac{S_{{\mathrm{o}}}-S_{\mathrm{nwr}}}{1-S_{\mathrm{wirr}}-S_{\mathrm{or}}} \bigg )^{n_{\mathrm{o}}} , \end{aligned}$$
(1b)

where \(K_{\mathrm{rw}}(S_{\mathrm{w}})\) and \(K_{\mathrm{ro}}(S_{\mathrm{o}})\) are, respectively, the water and oil relative permeabilities, \(S_{\mathrm{wirr}}\) and \(S_{\mathrm{or}}\) are the irreducible water and oil phase saturations, and \(n_{\mathrm{w}}\) and \(n_{\mathrm{o}}\) are the exponents for the water and oil phases. The irreducible saturations used are \(S_{\mathrm{wirr}} = S_{\mathrm{or}} = 0.2\) and the exponents are \(n_{\mathrm{w}} = n_{\mathrm{o}} = 2\). Gravity and capillary pressure were neglected for the SPE 10 models, following Christie and Blunt (2001). The modelled water viscosity was \(1 \times 10^{-3}\) Pa s and displacements were simulated using three different modelled oil viscosities of \(1 \times 10^{-3}\), \(3 \times 10^{-3}\) and \(10 \times 10^{-3}\) Pa s respectively. This allowed the effect of mobility ratio to be tested.

The oil composition in the H-Field varies with depth, with volatile oil at the base of the hydrocarbon column evolving to gas condensate at the top of the column with no sharp change in composition. Consequently, three phase relative permeability and capillary pressure curves were required in the model. The model consists of two regions in which the oil–water capillary pressure curves are markedly different. The water–oil and gas–oil relative permeability curves are uniform in the model across both regions

The Leverett J-Function was used to scale the water–oil and gas–oil capillary functions to the rock porosity and permeability. The modelled water viscosity is \(0.35 \times 10^{-3}\) Pa s while the hydrocarbon viscosity is modelled using the compositional viscosity correlation of Lohrenz et al. (1964) ranging from 0.07 to \(1.91 \times 10^{-3}\) Pa s. The remaining model parameters are summarized in Table 1. The water/oil displacement process in all models was simulated using commercial software.

Table 1 Model data
Full size table

2.5 Well Placements and Controls

Two different injector/producer well configurations were used in the synthetic models to investigate whether the hypothesis tested here holds regardless of the dominant flow direction relative to the modelled heterogeneity. The first induced predominantly radial flow, using an inverted five spot pattern, with a central injector and four producers at the corners of the model. The second induced predominantly linear flow, using a line drive pattern, with a row of injectors on one face and a row of producers on the other face of the model. The production wells flowed at a fixed bottom hole pressure and the injection well(s) flowed at a fixed injection rate. For the H-Field case, the flow directions were dictated by the well placements chosen for the field development. The production and injection wells were set to flow at the historic (measured) rates.

3 Methodology

The methodology used in this paper removed the cell-to-cell variability in petrophysical properties while keeping the grid resolution fixed. The flow behavior was then measured using a ‘Goodness of Fit’ between the well responses of the refined models and the corresponding base model. The following well responses were chosen for the synthetic models: (1) oil rate, (2) water cut (the fraction of the produced fluid that is water) and (3) the bottom hole pressure of the injection well(s). For the H-Field model, the well responses used for the Goodness of Fit were (1) gas rate; (2) water cut and (3) the bottom hole pressure of the injection wells. In addition to these ‘global’ parameters, the first arrival time of water (i.e. the time of water breakthrough) at the production wells was also tested in all models. Water breakthrough time is sensitive to the ‘local’ flow around the production wells.

3.1 Goodness of Fit

The normalized root-mean square deviation was used to determine the Goodness of Fit of an individual well response for each well

$$\begin{aligned} E_{\mathrm{r,w}} = 1 - \sqrt{\frac{\sum \nolimits ^T_{t=1}(y^{*}_{\mathrm{r,w},t}-y_{\mathrm{r,w},t})^2 }{ \sum \nolimits ^T_{t=1} (y^{*}_{\mathrm{r,w},t})^2 } }, \end{aligned}$$
(2)

where \(E_{\mathrm{r,w}}\) is the Goodness of Fit for a given response (r) in a given well (w), \(y^{*}_{\mathrm{r,w},t}\) is the well response on the base model and \(y_{\mathrm{r,w},t}\) is the well response on the refined model. The Goodness of Fit ranges between -inf (bad fit) and 1 (perfect fit). The averaged Goodness of Fit, for a given response over all the wells, is then given by \(E_{\mathrm{r}}\)

$$\begin{aligned} {E}_{\mathrm{r}} = \frac{\sum \nolimits ^W_{w=1}E_{\mathrm{r,w}}}{W}. \end{aligned}$$
(3)

Finally, \(\overline{E} \) is the Goodness of Fit for the whole reservoir, averaged over all of the wells and all of the well responses:

$$\begin{aligned} \overline{E} = \frac{\sum \nolimits ^N_{r=1}E_{\mathrm{r}}}{N}. \end{aligned}$$
(4)

3.2 Removal of Cell-to-Cell Variability

The workflow removed cell-to-cell variability in petrophysical properties in a two-step process. First, the base model was segmented into a number of discrete domains with internally homogeneous petrophysical properties. Uniform petrophysical values were assigned to each segment, calculated using the arithmetic average of the cell-based values within that segment. A series of simulations were then run on each segmented model. Flow behavior was compared using the Goodness of Fit described in the previous section (Eq. 4). The optimal number of segments was defined as the number of segments above which the Goodness of Fit between segmented and base models did not significantly increase. Second, the number of discrete domains in the optimum segmented model was reduced by using a clustering algorithm based on neighborhood filtering. The Goodness of Fit of the subsequent clustered model is then determined. The optimal window size (in terms of Goodness of Fit) depends on the connectivity of the domains and was determined by trial and error. The neighborhood algorithm is most suited to a Cartesian grid as it assumes simple neighborhood connections. Consequently, the clustering of the H-Field model (built on a corner point grid) was not optimal, but, as will be shown later, proved sufficient.

The workflow implemented the following steps:

Step 1.1 The log-space distribution of the horizontal permeability was used to divide the data into Z discrete segments. This was done by dividing the range [\(\min (\log _{10}(K_{\mathrm{H}})\)), \(\max (\log _{10}(K_{\mathrm{H}})\))], where \(K_{\mathrm{H}}\) is the horizontal permeability field, into Z non-overlapping, equal (in range) segments.

Step 1.2 An indicator function I(x) was introduced, such that \(I(x)=1\) if the horizontal permeability at grid-block x is a member of Segment 1, \(I(x)=2\) if the horizontal permeability is a member of Segment 2, and so on, i.e.

$$\begin{aligned}&I: X \rightarrow \{ 1,2,\ldots Z\} \nonumber \\&I(x) =\left\{ \begin{array}{ll} 1, &{} \quad \text {if } K_{\mathrm{H}}(x)\in \text {Segment 1 } \\ 2, &{} \quad \text {if } K_{\mathrm{H}}(x)\in \text {Segment 2 } \\ \vdots \\ \text {Z}, &{}\quad \text {if } K_{\mathrm{h}}(x)\in \text {Segment}\;Z \\ \end{array}\right\} , \end{aligned}$$
(5)

where x is the grid cell number and Z is the total number of segments. The same indicator function was used to segment the vertical permeability and porosity.

Step 1.3 The property values within each segment were homogenized to yield a single value for each. Porosity was averaged using the volume-weighted arithmetic mean. Horizontal permeability was averaged using an arithmetic average and vertical permeability using a geometric average.

Step 1.4 The indicator function was used to distribute the averaged properties on the original grid. The resulting model is now called a ‘segmented’ model. Figure 5a shows an example segmentation of the horizontal permeability in the Lower SPE 10 Model using \(Z=3\), i.e. three segments.

Step 1.5 The Goodness of Fit (\(\bar{E}\)) of the flow performance between the segmented and base models was then calculated based on the results of numerical simulations of flow in each model. Steps 1–4 were repeated, increasing the number of segments (Z) after each iteration, until no significant improvement in \(\bar{E}\) was seen. At this stage, this ‘optimal’ segmented model will contain fewer and larger domains than were present in the base model. However, various isolated domains comprising just one or a few cells remained in the model.

Step 2.1 A sliding neighborhood algorithm (Thompson and Shure 1995) was applied on the segmented model (Fig. 5). The algorithm scans the model, cell-by-cell, creating a ‘window’ containing a chosen number of cells around each input cell (in this example, the window size is \(3 \times 3 \times 3\)). The input cell is then replaced by the modal value of all cell values in the cube. This was repeated for each cell of the model. Figure 5c shows an example of a clustered model. The homogeneous regions are now referred to as domains. This step reduces the total number of domains by removing isolated groups of grid cells that can be a few cells wide.

Step 2.2 The flow performance of the clustered model was compared to that of the parent segmented model. The optimal window size was defined as the one which gave the closest fit (based on Goodness of Fit) to the optimal segmented model.

Step 2.3 The final Goodness of Fit (\(\bar{E}\)) of the flow performance between the optimal clustered and base models was then calculated based on the results of numerical simulations of flow in each model.

Fig. 5
figure 5

Segmentation and clustering of the Lower SPE 10 Model. a Segmented model with three segments. b \(3 \times 3 \times 3\) cube around an input cell during clustering. The central cell’s value is then replaced with the modal value of the neighboring cells; here a green permeability value is replaced by a blue one. c Clustered model. The cell-to-cell variability is removed, creating more geologically meaningful domains

Full size image

4 Results

Results from the Lower SPE 10 Model at each step of the workflow are first reported to illustrate application of the workflow and show how the final Goodness of Fit results were obtained. For the other test cases, only the final Goodness of Fit results are reported, as they show similar behavior to the Lower SPE 10 Model and the individual workflow steps were the same.

4.1 Lower SPE 10 Model: Segmentation and Clustering

We first consider the Goodness of Fit for the segmented (but not clustered) Lower SPE 10 Model for the various fluid properties tested and the inverted five spot well pattern (inducing primarily radial flow) (Fig. 6a). The solid lines denote the Goodness of Fit including only ‘global’ measures of well response i.e. production rates and pressures. The dashed lines denote the Goodness of Fit with both ‘global’ and ‘local’ well responses (water breakthrough time) included.

The results show that, even with only two segments, the ‘global’ Goodness of Fit for the segmented model (i.e. calculated without including water breakthrough) is greater than 88%, except for the case with the highest oil/water viscosity ratio (Fig. 6a). A two segment model captures the contrast between high permeability, channelized sandstones and low permeability mudstones. However, increasing the number of segments beyond two further improves the fit and, at 10 segments or more, the ‘global’ Goodness of Fit for the segmented model is greater than 93%; moreover, it is insensitive to the chosen oil/water viscosity ratio. Further increasing the number of segments beyond 10 yields a very slow rate of increase in the fit, tending towards one as the number of segments approaches the number of cells in the model. Note that 10 segments corresponds to a very small number of unique values of porosity and permeability compared to the 660,000 cells in the base model.

The inclusion of water breakthrough to produce a ‘global \(+\) local’ Goodness of Fit measure causes a large increase in the mismatch with a small number of segments, but the mismatch is again small at 10 segments or more (compare solid and dashed lines in Fig. 6a).

Fig. 6
figure 6

Goodness of Fit (GOF) as a function of segment number. a Lower SPE 10 Model with the inverted five spot well configuration; b Lower SPE 10 model with the line drive well configuration; c Upper SPE 10 Model with the inverted five spot well configuration; d Upper SPE 10 Model with the line drive well configuration; e H-Field model. Points with lines denote the segmented models without clustering. The different line markers denote the different fluid viscosity ratios. Dashed and solid lines denote the ‘global’ GOF (excluding water breakthrough) and ‘global \(+\) local’ GOF (including water breakthrough), respectively. Unfilled and filled triangles denote the ’global’ and ‘global \(+\) local’ GOF for the clustered models, with and without water breakthrough respectively

Full size image
Fig. 7
figure 7

Well and reservoir responses of the base, segmented and clustered models of the Lower SPE 10 Model and the inverted five spot well configuration with uniform oil and water viscosities. Plots a to d show oil production rate as a function of time for each producer well. Plots e to h show water production rate as a function of time for each producer well. Plot i shows the bottom hole pressure for the water injection well. Plots j to l show, respectively, the total reservoir oil production rate, water production rate and cumulative oil production as a function of time

Full size image

The ‘global’ Goodness of Fit for the same test case with the line drive well pattern (inducing primarily linear flow) shows similar behavior, except that the fit for the two segment model is lower, especially with the higher oil/water viscosity ratio (Fig. 6b). Including water breakthrough has a more limited effect on the ‘global \(+\) local’ Goodness of Fit in this case with linear flow than it did with primarily radial flow. However, at 10 segments or more, the ‘global \(+\) local’ Goodness of Fit for the segmented model is again high and identical irrespective of the oil/ water viscosity ratio. Further increasing the number of segments beyond 10 yields a very slow rate of increase in the fit. The 10-segment model was therefore considered optimal for the Lower SPE 10 Model and was taken into the clustering step.

The clustering algorithm was applied using a window size of \(3 \times 3 \times 1\) twice sequentially, yielding approximately 2,000 individual domains. The Goodness of Fit of the clustered models is shown by the triangles in Fig. 6. Clustering causes the Goodness of Fit to decrease by a few percent (up to 4%), as confirmed by examination of the individual well and reservoir responses (Fig. 7). For example, the model with just 10 segments (and, in Fig. 7, equal oil and water viscosities and the inverted five spot well pattern) provides a close match to the oil and water rates in each production well (Fig. 7a to h). Not only are the production rates closely matched, but the water breakthrough time is also captured in each well. The bottom hole pressure in the injection well (Fig. 7i) shows very good agreement with the base model. Most of the discrepancy in a given response occurs at relatively early time (first 400 days, or 20% of simulated production). The 10-segment model also provides a good match to the reservoir responses of oil and water production rate, and cumulative oil produced (Fig. 7j to l).

The key geological heterogeneity that the Lower SPE 10 Model needed to capture was the connectivity of the high permeability channelized sandbodies (Fig. 8; see also Fig. 1). The two segment model did reasonably well in this respect, capturing the high permeability channelized sandstones and low permeability mudstones. However, the geostatistical method used in the original model to populate permeability within lithofacies created some connections between channelized sandstones via one or a few cells in the mudstone that were assigned intermediate permeability values. It is not clear that these connections are realistic and consistent with the underlying geologic concept, or are merely an artefact of the geostatistical modelling algorithm. Nonetheless, increasing the number of segments to 10 preserved these intermediate permeability connections, yielding very similar connectivity and, therefore, sweep patterns in the base and 10-segment models (Fig. 8). The clustering step did not significantly modify the permeability distribution. Thus a model with just 10 unique values of porosity and permeability, grouped into a few thousand domains, was able to replicate the two-phase flow behavior of a model with 660,000 cells, each of which contains a different value of porosity and permeability to its neighbor, irrespective of the well configuration or fluid properties tested.

Fig. 8
figure 8

Horizontal and vertical sections through the Lower SPE 10 Model showing horizontal permeability and water saturation at a given timestep during production via the inverted five spot well configuration with uniform oil and water viscosities. Plots a and b show, respectively, the permeability at a horizontal section through the base and 10-segment and clustered models. Plots c and d show the water saturation at a given timestep on the same horizontal section. Plots e and f show the vertical section through the base and 10-segment and clustered models, respectively. Plots g and h show the water saturation at the same vertical section

Full size image

4.2 Upper SPE 10 Model: Summary of Goodness of Fit Results

The Goodness of Fit results for the same study on the Upper SPE 10 Model show similar behavior to those of the Lower SPE 10 Model (Fig. 6c to d). The ‘global’ Goodness of Fit for two segments is lower than for the Lower SPE 10 Model, but increases as the number of segments increases. At 10 segments or more, the ‘global’ Goodness of Fit for the segmented model is high (\(> 90\)) irrespective of the oil/water viscosity ratio or well configuration. Further increasing the number of segments beyond 10 yields a very slow rate of increase in the fit. Including water breakthrough in the ‘global \(+\) local’ Goodness of Fit in this case had a larger effect for the highest viscosity ratio and the inverted 5 spot well pattern. The ‘global \(+\) local’ Goodness of Fit at 20 or more segments was lower. This was found to be due to a few cells very close to the wells, and the effect for the line drive well configuration was negligible. Overall, for both well patterns, the 10-segment model was considered optimal for the Upper SPE 10 Model and was taken into the clustering step.

The clustering algorithm was run with a window size of \(6 \times 6 \times 1\) yielding approximately 850 individual domains. The base and final clustered models are shown in Fig. 9a, b. The Goodness of Fit of the clustered models is shown by the triangles in Fig. 6. Clustering reduced the Goodness of Fit by a few percent but, similar to the Lower SPE 10 Model, the decrease was small.

Fig. 9
figure 9

Horizontal sections through the Upper SPE 10 Model and the H-Field Model showing horizontal permeability. Plot a shows the base Upper SPE 10 model; plot b shows the Upper SPE 10 model with 10 segments after clustering. Plot c shows the base H-Field model; plot d shows the H-Field model with 20 segments after clustering

Full size image

Unlike the Lower SPE Model, the Upper SPE 10 shows no strongly correlated permeability distribution and there are no clear geological heterogeneities that should be captured in the model. Thus the two segment model yields a significantly lower Goodness of Fit for the Upper SPE 10 Model as compared to the Lower SPE 10 Model. However, at 10 segments, the segmented and clustered models again yield a close match to the base model. A model with just 10 unique values of porosity and permeability, grouped into hundreds of domains, was able to replicate the two-phase flow behavior of a model with 462,000 cells, each of which contains a different value of porosity and permeability to its neighbor, irrespective of the well configuration or fluid properties tested.

4.3 H-Field Model: Summary of Goodness of Fit Results

We finish by reporting the Goodness of Fit results for the same study on the H-Field Model (Fig. 6e). Despite the significantly higher complexity of the H-Field, in terms of geological heterogeneity, number of wells, and multiphase flow, the results are similar to the synthetic SPE 10 test cases. For both the ‘global’ and ‘global \(+\) local’ Goodness of Fit measures, the match for two segments is low, but increases as the number of segments increases. At 20 segments or more, both the ‘global’ and ‘global \(+\) local’ Goodness of Fit for the segmented model are high. Further increasing the number of segments beyond 20 yields a very slow rate of increase in the fit. The 20-segment model was therefore considered optimal for the H-Field Model and was taken into the clustering step. The model was clustered using three iterations of the algorithm, with a window size of \(5 \times 5 \times 5\), then two with \(5 \times 5 \times 1\), yielding approximately 850 individual domains. The H-Field model was more sensitive to clustering than the synthetic SPE 10 test cases; however, the mismatch was still small.

The H-Field model shown in Fig. 4 is only one geologic realization. Applying the same statistical description of the petrophysical properties, many more realizations are possible. Fig. 10 shows the ‘global \(+\) local’ Goodness of Fit for 10 realizations of the H-Field model. While there is significant sensitivity with \(<20\) segments, the increase in Goodness of Fit is very slow with 20 or more segments. Thus the observed behavior is not specific to a single model realisation.

The H-Field model captures many geological heterogeneities, including complex structure with several fault blocks, and diagenetically modified permeability values controlled by a multi-stage filling history during initial charging with oil. Perhaps because of this, more segments and domains were required to reproduce the base model behavior. Nevertheless, a model with just 20 unique values of porosity and permeability, grouped into a few hundred domains, was able to replicate the multiphase, multicomponent flow behavior of a model with 146,000 cells, each of which contains a different value of porosity and permeability, to its neighbor.

Fig. 10
figure 10

Goodness of Fit of segmented models without clustering, as a function of segment number for 10 realizations of the H-Field. a ‘Global’ Goodness of Fit excluding water breakthrough; b ‘Global \(+\) local’ Goodness of Fit including water breakthrough

Full size image

5 Discussion

5.1 Is Cell-to-Cell Scale Variability Necessary in Reservoir Models?

The results presented here indicate that grid-based reservoir models containing many unique values of petrophysical properties (of order hundreds-of-thousands in the examples tested) varying on a cell-to-cell basis can be collapsed into a much smaller number of larger but more geometrically complex domains which are internally homogeneous, irrespective of the reservoir geology, fluid properties or well configuration. Thus, while it is always possible that testing of further cases may contradict this, the evidence presented here suggests that cell-to-cell scale variability is not necessary in reservoir models. We note that the models tested here were of clastic reservoirs. It is well known that carbonate reservoirs often differ from clastic reservoirs in the type, scale and geometry of key geologic heterogeneities and how these impact on flow (e.g. Borgomano et al. 2002; Pranter et al. 2006; Fitch et al. 2014). Nonetheless, as discussed in the Introduction, the cell-to-cell variability observed in reservoir models is an artefact of applying pixel-based modelling methods on a grid of arbitrary resolution and geometry, rather than a consequence of the underlying geologic heterogeneity. We argue that, irrespective of depositional environment, modelling workflows should focus on capturing the geologic domains that represent spatially correlated variability in petrophysical properties (particularly permeability) in the most realistic and computationally efficient manner. The correlated variability reflects the underlying geological heterogeneity that is important to flow.

In the approach adopted here, ‘domain-based’ models were created by segmenting and clustering existing grid-based models. However, we adopted this approach only to test the effect of removing cell-to-cell scale variability from grid-based models, and we do not recommend this workflow for model construction. The important finding here is that capturing domains that represent spatially correlated variability in petrophysical properties is key. This finding offers the prospect of alternative modelling workflows to the grid-based methods used ubiquitously to date.

5.2 Implications for Surface-Based Geologic Modelling

One approach to capture key geologic domains is to model the surfaces that bound them. Surface-based modelling (or boundary representation) has been suggested as an alternative to grid-based modelling to capture a wide variety of different geological heterogeneity types, including faults (Julio et al. 2015), stratigraphy (Sylvester et al. 2015), facies and lithological boundaries (Hassanpour et al. 2013; Ruiu et al. 2016), and fractures (Belayneh et al. 2006; Paluszny et al. 2007). Indeed, all geological heterogeneity that impacts the spatial distribution of petrophysical properties can be modelled using surfaces (Jackson et al. 2013). Moreover, the surfaces can be modelled using a parametrized description without reference to a pre-defined grid or mesh (Graham et al. 2015b; Jacquemyn et al. 2018), significantly increasing computational efficiency. A mesh is created only when required for numerical calculations and preserves the surface-based representation of geological heterogeneity (Jackson et al. 2015). Petrophysical properties in each element of the mesh are inherited from the parent domain i.e. no further property modelling is required. If significant heterogeneity is expected within a domain, additional surfaces can be introduced. Therefore such models omit cell-to-cell scale variability. The models are not less complex or less realistic than grid-based models; as we show here, they simply omit variability that is not important to flow. Indeed, surface-based models can much better capture correlated heterogeneities with complex geometries that are important to flow (Jackson et al. 2009; Deveugle et al. 2014; Graham et al. 2015b).

5.3 Surface- Versus Grid-Based Models: Capturing Key Geologic Heterogeneities

In the approach adopted here, we assumed that each conventional, grid-based, base model yielded the ‘correct’ flow behavior for the corresponding reservoir, and the removal of cell-to-cell scale variability was tested by comparison against the flow behavior of the base model. Yet the base models may not have captured important geological heterogeneities present in the reservoir. For example, the Upper SPE 10 Model was originally constructed to represent reservoirs deposited in a shallow-marine environment. Numerous studies have shown that facies architecture is a key control on flow in such reservoirs, with facies types corresponding to order-of-magnitude contrasts in horizontal and vertical permeability (Ringrose et al. 2008; Deveugle et al. 2014). The thickness and relative proportion of each facies type strongly controls flow, along with interfingering of facies along inclined surfaces termed clinoforms. Calcite cements are also a common feature of such reservoirs and are often associated with clinoform surfaces (Sech et al. 2009; Jackson et al. 2009). Yet, despite the importance of clinoform surfaces in controlling the spatial arrangement of key geological heterogeneities, they are typically neglected in reservoir models because their inclined geometry relative to reservoir bounding surfaces means they are difficult to represent using grid-based approaches. There is no evidence that the key heterogeneities described above are captured in the Upper SPE 10 Model (Fig. 3) and, despite its ubiquitous use in later studies, it likely provides a poor representation of flow in shallow-marine reservoirs (White and Barton 1999; Sech et al. 2009; Jackson et al. 2009; Enge and Howell 2010; Deveugle et al. 2014; Graham et al. 2015b amongst many others).

Figure 11a, b shows a surface-based model of the SPE 10 geology (Jacquemyn et al. 2018). The upper part of the model represents a reservoir deposited in a shallow-marine environment and corresponds to the Upper SPE 10 Model; the lower part of the model represents a fluvial reservoir and corresponds to the Lower SPE 10 Model. The surface-based model was constructed entirely using parametrized surfaces (in this case, Non Uniform Rational Basis Splines or NURBS) without reference to a pre-defined grid or mesh. The surfaces have been assembled to create watertight volumes that correspond to geologic domains. Domains of different facies/lithologies are shown by the different colors. In the upper part of the model, the domains correspond to different shallow-marine sandstone facies types; in the lower part of the model, domains corresponding to channelized sandstones. Domains representing non-reservoir mudstones have not been visualized, but are present in the model.

Fig. 11
figure 11

Surface-based model inspired by SPE 10 model 2. a Full model; b Lower, fluvial section; from below (Jacquemyn et al. 2018)

Full size image

The surface-based representation of the Upper SPE 10 geology includes clinoforms and the shallow marine facies interfinger along these surfaces; the model also includes patches of calcite cement along clinoform surfaces. Thus, although the model does not include cell-to-cell variability in petrophysical properties, it includes the key correlated heterogeneities observed to impact flow (Jackson et al. 2009; Graham et al. 2015b). Similarly, the surface-based representation of the Lower SPE 10 geology includes channelized sandbodies of varying width and sinuosity and, although the model does not include cell-to-cell scale variability, it captures the spatial distribution and connectivity of the sandbodies, which is the key heterogeneity observed to impact flow. Such a model can be used as the initial input to new generation simulators that used adaptive unstructured meshes to optimize computational effort when solving the governing flow equations (Mostaghimi et al. 2015; Jackson et al. 2015; Gomes et al. 2017b; Salinas et al. 2018).

5.4 Surface-Based Modelling for Complex Structure and Flow Simulation

Surface-based models are also able to capture structural heterogeneity such as recumbent folds, intersecting normal faults or listric growth faults (see Jacquemyn et al. 2018). These are often very difficult to represent using conventional pillar grids. Figure 12 shows a surface-based model of conjugate normal faults offsetting stratigraphic layers of contrasting permeability. If required, the model can be discretised using tetrahedral unstructured meshes. Figure 12d shows an example of fluid flow simulated using the Imperial College Finite Element Reservoir Simulator (ICFERST; Jackson et al. 2015; Salinas et al. 2017. In this example, the model is initially saturated with oil and irreducible water, and injection of water is simulated through the left boundary, displacing oil toward the right boundary with all other boundaries closed. Dynamic adaptive mesh optimisation (Pain et al. 2001) is used to optimize the geometry of the mesh to a solution field of interest; here, the mesh is optimized to the water saturation field. High mesh resolution is therefore focused at the saturation front. Ahead and behind the front, the mesh is kept coarse.

Fig. 12
figure 12

Surface-based model of a layered reservoir offset by conjugate normal faults. a Side view; b Top view. c Perspective view. The model is discretised using a coarse unstructured tetrahedral mesh, with under 20,000 elements. Initially saturated with oil and irreducible water. d Water is injected cross the left boundary with production from the right boundary. All other boundaries are closed. As the displacing water front propagates through the model, dynamic adaptive mesh optimisation is used to change the geometry of the mesh. Computational resources are therefore allocated more efficiently around areas of interest (Jacquemyn et al. 2018)

Full size image

There are a few outstanding issues for surface-based modelling methods. While it is possible to accurately represent the complex geometries of geology, conditioning these models to control data, such as well or seismic data, remains a challenge (Bertoncello et al. 2013). Another challenge is to update surface-based models to better match observed behavior (history match). These challenges are by no means insurmountable (e.g. Trani and Graham 2016); however, further work is needed to develop workflows to condition and history match surface-based models.

6 Conclusions

We show that grid-based reservoir models containing a large number of petrophysical property values that vary on a cell-to-cell basis can be collapsed into a much smaller number of larger and internally homogeneous but more geometrically complex geologic domains, irrespective of the reservoir geology, fluid properties or well configurations tested. We conclude from this that cell-to-cell scale variability is not necessary in reservoir models to capture multiphase flow. Although the geometrically complex domains can be represented in standard grid-based models, surface-based reservoir modelling methods are more flexible and better able to capture the complex geometries of geologic domains in a realistic and computationally efficient manner. Numerical solutions to the governing flow equations can be obtained in such models using adaptive unstructured meshes.