Determination of Spatial Heterogeneity Index for Naturally Fractured Reservoirs

Abstract

Modeling variations in reservoir properties as a function of space is an important task in evaluating reservoir heterogeneity. While it provides considerable value to capture the relative change in properties in the entire reservoir, it is an overwhelming task to perform in 3D. Various scientific approaches exist that facilitate measuring static reservoir heterogeneity along the wellbore. However, they are very limited and confined to small static datasets that are very often not representative. This chapter introduces a new approach to quantify a metric of 3D static reservoir heterogeneity. It describes the methodology to calculate a spatial heterogeneity index that encapsulates complex features to benchmark formations, zones, reservoirs, and fields. The introduced workflow expands on Lorenz’s scientific approach as the basis for this methodology to develop the heterogeneity index. The systematic approach excels in rendering a universal systematic 3D spatial heterogeneity index where typically data are not available, representative of overall prolificacy for a particular region of interest, zone, or the entire reservoir. The chapter also envisions possible applications, where the spatial heterogeneity index is utilized as a 3D conditioning attribute in dynamic model calibration and optimization.

1. Introduction

Spatial heterogeneity is defined as the variation in rock properties in a subsurface reservoir. In subsurface modeling, it is extremely important to capture the dispersion of rock properties to obtain a representative 3D static model for dynamic calibration. Static models are constructed in 1D using geophysical, petrophysical, and geological data which are then distributed using geostatistical principles based on interpolation (e.g., kriging) or simulation (e.g., sequential Gaussian simulation, SGS) to 3D models. So, all modeling techniques intend to capture intra-well uncertainties.

Dykstra and Parsons [1] were among the first to capture reservoir variations (V) to describe the degree of heterogeneity for a reservoir. They proposed a normalized scale (Figure 1) that ranged between zero and one; with zero being a completely homogeneous system, and one for a completely heterogeneous system as per Eq. (1):

Figure 1.

The Dykstra-Parsons coefficient calculation method [2].

$$\mathit{\text{Permeability Variation}}\left(V\right)=\frac{\mathit{\text{stdev}}\left(logk\right)}{\mathit{avg}\left(logk\right)}=\frac{log{\left(k\right)}_{P50}-log{\left(k\right)}_{P84.1}}{log{\left(k\right)}_{P50}}$$E1

Many enhanced oil recovery studies used Dykstra-Parsons’s coefficient as an indicator for the level of heterogeneity [2].

In 1950, Schmalz and Rahme [3] proposed a single value to represent a reservoir’s permeability distribution defined as the Lorenz coefficient. The Lorenz coefficient is constructed from permeability, thickness, and porosity and is quantified as the area between the unity line and the actual reservoir data on a plot between the fraction of total storativity (cum fraction ØH) and the fraction of total flow capacity (cum fraction KH). Here, Ø and K correspond to porosity and permeability, respectively and H represents the thickness. A Lorenz coefficient of 0 describes a perfectly uniform homogenous reservoir (i.e., y = x line), while a value of 1 means a completely heterogeneous reservoir with a big contrast in property variations. Figure 2 shows synthetic examples of two datasets, the black dots represent a reservoir with higher reservoir property variations, deviating from the unit line, while the blue dots represent a reservoir with homogenous characteristics.

Figure 2.

Cum ØH fraction vs. cum KH fraction (Lorenz plot).

However, the Lorenz coefficient method is performed on available datasets, which are typically limited to certain parts of the reservoir that are targeted for development. In addition, applications of these approaches are limited in their extension to full 3D static and dynamic models.

While Liu et al. [4] introduce dispersion variance as a spatial-engineered feature that accounts for heterogeneity within the spatial context, to improve data-driven predictive ML-based models for example, pre-drill prediction and uncertainty quantification, the presented 3D heterogeneity approach uses a modified version of the Lorenz’s coefficient as basis for the vertical and lateral heterogeneity and connectivity mapping in full physics subsurface simulations. We first introduce the heterogeneity index workflow and then present applications to baseline and advanced versions of reservoir simulation models. We conclude with observations and forward-looking statements.

2. Methods and models

This work builds on the same foundations set by Schmalz and Rahme [3] in 1950 for the classical Lorenz approach that provides a systematic methodology to evaluate the change in reservoir properties. The Lorenz approach has become the standard especially for upscaling and retaining log heterogeneity; however, it is limited to the wellbore and does not have any spatial representation, nor does it really represent a reservoir’s heterogeneity globally due to lack of data. This often requires the model to be defined as a “layer cake” model, with fairly uniform layer properties because of not having a representative dataset (i.e., vertical penetrations when the majority of wells are drilled horizontally) to properly capture the spatial heterogeneity. In contrast, a modified version of the Lorenz plot derived herewith incorporates depth as well as placement in space. The implementation of the modified approach allows for more flexibility to characterize heterogeneity vertically and horizontally.

The reference model used to showcase the approach is the SPE10 model [5] that is contained within a regular Cartesian grid with 60 × 220 × 85 (1.1 million cells). The base SPE10 model consists of two formations that manifest large permeability contrast (approx. 10 orders of magnitude) and are qualitatively different. The reservoir model resembles water flooding from a vertical well in the center and four vertical producers at each corner as depicted in Figure 3.

Figure 3.

SPE10 porosity, permeability, and well locations.

The process of generating the 3D spatial heterogeneity index starts with the degree of heterogeneity. For a highly heterogeneous reservoir, data sampling frequency should be very high, and for a homogenous reservoir, a representative 3D spatial heterogeneity index can be evaluated on a much lower sampling frequency. The sampling frequency refers to the availability of data, gathered either from existing vertical wells or extracted from a 3D static model. If data is scarce due to a limited number of vertical wells, the 3D static model can be utilized to extract the data systematically for the evaluation.

Utilizing the 3D static model, data are extracted vertically across the entire reservoir, including the aquifer (as applicable). These data are used to construct the typical Lorenz plot (porosity, permeability, thickness, and zonation). Additional data on naturally occurring features or phenomena that dominantly impact fluid flow such as faults, fractures, high permeability streaks, baffles, should be extracted as well. The classical Lorenz coefficient is calculated at every data point in the reservoir which could be controlled at certain zones, formations, or expanded to cover multiple reservoirs within the same field. After calculating the vertical heterogeneity indices at every data point, a 2D map is created by correlating all the vertical values together. The magnitude of the vertical heterogeneity index represents how properties along that vertical reservoir column vary.

The workflow allows for simultaneous heterogeneity evaluation for any pre-defined zone, combined neighboring zones, combined scattered zones, reservoir, or field. Applications of this methodology can aid in well placement, field development, and advanced completion optimization. The workflow is schematically depicted in Figure 4.

Figure 4.

Heterogeneity index workflow.

3. Lorenz plot

As illustrated in Figure 4, the Lorenz coefficient represents the heterogeneity along the wellbore per degree of vertical resolution and is used to quantify how much the corresponding properties change compared to the entire wellbore. An example of synthetic porosity and permeability logs, calculated along the wellbore, is shown in Figure 5. The following equations are used to calculate the parameters needed to create the Lorenz plot.

Figure 5.

Synthetic porosity & permeability logs.

$$\mathit{\text{Storage Capacity}}\left(ØH\right)=Ø\ast \mathit{\text{Thickness}}$$E2

$$\mathit{\text{Flow Capacity}}\left(\mathit{KH}\right)=K\ast \mathit{\text{Thickness}}$$E3

$$\mathit{\text{Storage Capacity Fraction}}\left(Ø\mathit{HF}\right)=\frac{Ø\ast \mathit{\text{Thickness}}}{\mathit{\text{Total}}ØH}$$E4

$$\mathit{\text{Flow Capacity Fraction}}\left(\mathit{KHF}\right)=\frac{K\ast \mathit{\text{Thickness}}}{\mathit{\text{Total}}\mathit{KH}}$$E5

$$\mathit{\text{Cumulative Storage Capacity Fraction}}\left(Ø\mathit{HF}\right)=\sum _1^nØ\mathit{HF}$$E6

$$\mathit{\text{Cumulative Flow Capacity Fraction}}\left(\mathit{KHF}\right)=\sum _1^n\mathit{KHF}$$E7

For that, storage capacity and flow capacity are used to calculate the fractions for each interval n as well as the cumulative flow fraction and cumulative storage fraction (Figure 6), summed over all intervals n. The heterogeneity index is then the area under the curve and y = x line. The closer the data to the unity line, the more homogeneous the properties along the wellbore and vice versa.

Figure 6.

Cumulative ØH fraction vs. cumulative KH fraction for synthetic wells.

Utilizing the Lorenz plot, specific zones, formations, or multiple reservoirs can be calculated. For example, Figure 7 shows different scenarios that can be assessed using this approach. The vertical solid line extracts data for the entire reservoir, while the dashed arrows show scenarios where only specific parts of the well penetrating that zone/reservoir are considered. Moreover, a combined effect of different indices can be acquired using a modified workflow where different separated zones within the same reservoirs can be assessed and measured.

Figure 7.

The heterogeneity index approach could be applied to the full reservoir section, specific individual zone(s), or multiple separated zones.

Figure 8 shows a lower frequency data extraction compared to a higher frequency which corresponds to the resolution and accuracy for the 3D model. It is essential to extract all available data as well as structure information to control the level of details for various reservoirs, zones, or the entire model.

Figure 8.

Comparison between low-frequency data sampling (left), and high-frequency data sampling (right) which is the first step in the heterogeneity index approach.

By utilizing the 3D static model, data is extracted at different parts of the reservoir and the Lorenz assessment of the heterogeneity is done at every part of the reservoir. To better evaluate heterogeneity, it is more suitable to assess each zone individually, in comparison with the entire reservoir. For demonstration, the SPE10 model was assigned different zones, three in the top part, and four in the bottom part, as shown in Figure 9. The top three zones are part of the “Shallow-marine Tarbert formation” and the bottom four zones are part of the “Fluvial Upper Ness formation,” with subzone classification presented in Table 1.

Figure 9.

SPE10 3D porosity array (left) with corresponding formation assignment (right).

Layer	Zone	Formation
1–10	Zone 1	Shallow-marine Tarbert
11–21	Zone 2
22–35	Zone 3
36–47	Zone 4	Fluvial Upper Ness
48–59	Zone 5
60–70	Zone 6
71–85	Zone 7

Table 1.

SPE10 formation assignment per layer and per zone.

The Lorenz plot can be constructed for each individual zone, multiple zones, formation, or the entire reservoir. The more data that are available in the spatial domain, the more detailed the spatial heterogeneity index. Hence, the detailed assessment of the change in properties can be for each zone, rather than as an overall index for the entire reservoir.

4. Results and discussion

4.1 Scenario 1: SPE10, original properties

To demonstrate the methodology, the SPE10 model is utilized with various scenarios. The first scenario utilizes the original model properties without any modification. Figures 10 and 11 show the 3D property histograms for porosity and permeability for top and bottom formations separately, while Table 2 shows the 3D property statistics for porosity and permeability arrays in the Shallow-marine Tarbert and Fluvial Upper Ness formations.

Figure 10.

Permeability and porosity histograms for the shallow-marine Tarbert formation.

Figure 11.

Permeability and porosity histograms for the Fluvial Upper Ness formation.

Formation	Property	Min	Max	Mean	Std
Shallow-marine Tarbert	Permeability	0.0007	20,000	356.5	1596.5
	Porosity	0	0.5	0.1911	0.0803
Fluvial Upper Ness	Permeability	0.0007	20,000	351.2	1306.2
	Porosity	0	0.4	0.1587	0.0975

Table 2.

3D property statistics for porosity and permeability for the SPE10 model.

4.2 Shallow-marine Tarbert (top formation)

4.3 Fluvial upper ness (bottom formation)

5. Construction of the spatial heterogeneity index

For the first scenario, three horizontal wells are designed and completed in zones 40–86 as well as two vertical wells completed across the entire reservoir, from top to bottom, Figure 12. This is done to establish a benchmark case to illustrate how heterogeneity is typically measured using existing wells.

Figure 12.

Generated vertical and horizontal well locations.

Conventional approaches to determine the heterogeneity of a reservoir typically consider vertical wells. The corresponding heterogeneity index might be suboptimal in case there are not enough wells to provide vertical coverage. Figure 13 shows synthetic logs generated from the intersection between wells and the penetration with the 3D grid. Wells 1, 2, and 3 represent slanted wells that are completed in the bottom zones, while wells V1 and V2 represent vertical wells going through the entire reservoir.

Figure 13.

Synthetically generated logs for the five wells (two vertical, three horizontals).

In this approach, the entire wellbore section is considered even if the well is only completed in the bottom zones given that properties are available through logs or a 3D model. Since the entire well is considered, the well heel typically starts at 0 on the cumulative ØH Fraction vs. Cumulative KH Fraction, and the toe would be at the far end of the scale at 1. In Figures 13 and 14, all wells were considered, regardless of whether they were vertical, horizontal, or slanted.

Figure 14.

Lorenz plots calculated for the five wells.

Typically, in conventional methods, particularly for modeling spatial property distribution, only vertical wells are considered. The horizontal wells are not considered, as they would skew the probability distribution since they are by design placed in regions with prolific reservoir quality.

If we analyze the Lorenz plot for well WELL0002 (Figure 15), we can make a few observations such as:

• On the cumulative ØH fraction vs. cumulative KH fraction plot, the shallower zones are displayed closer to the (0.0) point.

• Increase in KH from the heel of the well to the deeper part, which is denoted by the jumps, or kicks in the cumulative ØH Fraction vs. cumulative KH fraction plot over a negligible increase in ØH.

• Toward the toe of the well, no observed increase in KH or ØH. So, in this part of the reservoir, the top zones are more prolific than the deeper zones.

Figure 15.

Lorenz plot for WELL0002 (left), and mapping using WELL0002 data only (right).

For the SPE10 case, with two vertical and three horizontal wells, for any heterogeneity analysis to be done, only available data across that zone/reservoir would be considered. For example, the heterogeneity evaluation for the top formation, Shallow-marine Tarbert, considers only two wells with available data across this formation, while the bottom formation, Fluvial Upper Ness, will consider the available data for the five wells. However, this typically provides the wrong notion, that when there are not enough data, then the available data provide good representations of the formation, or the reservoir. Figure 16 shows heterogeneity mapping for the Shallow-marine Tarbert formation considering two wells and the Fluvial Upper Ness formation that was built using the available five wells.

Figure 16.

Heterogeneity index for shallow-marine Tarbert formation (left) using available data from two wells, while heterogeneity index for Fluvial Upper Ness formation from available five wells data (right).

The question remains, will five data points be enough to construct a representation of the heterogeneity changes for any reservoir? The answer to that depends on many parameters such as the rate of change of porosity and permeability, fracture density, tar that need to be factored in the heterogeneity assessment. The more complexity, the more data points are needed to capture their impact.

The area under the curve is calculated to establish the Lorenz coefficient value at a particular vertical data column. This step is then performed across all available data columns, and the outcomes denote a vertical heterogeneity index at the specified frequency. The horizontal heterogeneity map is created using interpolation between all the anchor points for a particular subzone, zone, formation, reservoir, multiple reservoirs, or field. Figure 15 also shows the anchor vertical heterogeneity index point used to map the heterogeneity over the entire reservoir. However, one point or multiple points are not sufficient to properly map a representative index map. The same approach is repeated with a higher data frequency, and the resultant map is displayed in Figure 17. Figure 17 has three maps, the left one shows the heterogeneity map across the entire reservoir, while the middle one and right ones show the heterogeneity map across the top and bottom formations, respectively.

Figure 17.

Heterogeneity Index for the entire reservoir (left), and split heterogeneity index map top formation (middle) and bottom formation (right). While the heterogeneity indices have relatively similar values, the mapping provides a completely different understanding.

It is important to note that the spatial heterogeneity index described in scenario 1 only incorporates porosity and permeability values without any consideration to various reservoir complexities, such as faults, fractures, baffles.

Looking at the SPE10 model closely, the Shallow-marine Tarbert formation has a relatively smooth permeability, while the fluvial Upper-Ness formation has a more structured distribution. So, it is expected to have a higher heterogeneity index in the bottom formation to reflect the irregular permeability distribution of that formation. However, while inspecting the values of the heterogeneity indices for top and bottom formations, one notices quite the opposite. Bottom formation exhibits a lower heterogeneity index value when compared with the top formation. There is a good explanation for that, which also suggests that the index should not be used as is in certain situations. Subzonation within the same formation could mask the heterogeneity for the entire formation. Another reason is the averaging. For example, in the bottom zone, there are a lot of variations denoted by the contrast in colors. This contrast also masks the heterogeneity assessment, and in these situations, the evaluation should be done on a zone-by-zone basis.

Figure 18 shows the zone-by-zone heterogeneity index mapping for individual zones separately. While there are not a lot of variations within each zone, the major heterogeneity variations would be attributed to the change from one zone to the other reflecting the sedimentological and depositional environments.

Figure 18.

Zone by zone heterogeneity indices and maps starting from the Shallow-marine Tarbert formation and its associated individual three zones (left), to the Fluvial Upper Ness formation and its associated individual four zones (right).

Applications of the approach can be realized during the model construction, which accounts for spatial subsurface uncertainty and as such considers many probable realizations. Figure 19 contains three different porosity realizations that were created from varying variogram parameters, such as seed, minor, major, and vertical variogram directions. The formations (Shallow-marine Tarbert and Fluvial Upper Ness) were not retained during the construction of the cases below, which are constructed for illustration purposes. In an actual scenario, as compared to the synthetic SPE10 case, hundreds of possible scenarios could be created with equal probability. While porosity realizations in Figure 19 are different, permeability realizations are identical due to the independence of volume calculations from permeability.

Figure 19.

Layer 1 view of sample porosity realizations that were generated by varying seed, minor, major, and vertical variogram directions while retaining the same permeability realization.

To focus on the impact of different realizations, heterogeneity indices were calculated for the three scenarios from Figure 19, and the result can be seen in Figure 20 for the Shallow-marine Tarbert formation. While the results have major similarities, they are not totally unexpected. Since only porosity realization is slightly different, and the main synthetic input data is the same, the results show a close resemblance. However, even if the heterogeneity indices have had similar values due to averaging, the mapping in return provides considerable value in evaluating the reservoir heterogeneity.

Figure 20.

Shallow-marine Tarbert heterogeneity evaluation for a sample of porosity realizations. Permeability realization is identical in the three cases. The three cases exhibit very similar characteristics because the impact of the change in porosity alone is not as significant while the permeability array is the same.

6. Scenario 2: advanced applications of the Lorenz plot

Up until this section, all heterogeneity representations are attributed to more continuous subsurface properties, such as porosity, permeability, zonation. Herewith we introduce the advanced application of the Lorenz approach and the spatial heterogeneity index method.

One major modification to the approach is to account for reservoir abnormalities, such as high permeability streaks, faults, fractures, baffles, tar, and other features impacting fluid flow of which spatial appearance in the subsurface is more discrete by nature. Every feature is brought to a common yardstick, which is the impact on fluid flow. For example, high permeability streaks, fractures, and discrete fracture networks (DFNs) would represent a preferential flow path and as such favor fluid flow, while tar would act as a barrier and impact fluid flow negatively. These features are considered and introduced based on their impact on porosity and permeability.

To demonstrate the application, the base SPE10 model is used with the exact same parameters as presented earlier (Scenario 1). A fracture realization is created stochastically and embedded in the model (Figure 21).

Figure 21.

Stochastically generated fracture realization for the SPE10 model.

In this realization, the vertical wells in the base SPE 10 model do not intersect any fractures. So, even if conventional heterogeneity approaches account for fractures, they would not impact the heterogeneity evaluation due to intersection, or the lack of it. Hence, the sampling methodology comes into effect (Figure 8). To demonstrate the difference between wells that only account for porosity and permeability, and those accounting for complex reservoir features. Figure 22 shows a well that is completed along a fracture plane and the corresponding heterogeneity index for both scenarios, original properties, and fracture scenario.

Figure 22.

Vertical well completed across a major fracture across the entire reservoir (left), and the corresponding Lorenz plot (right).

Figure 22 shows that accounting for the fractures would completely change the heterogeneity index for that well. The flow is dominated by the intersecting fracture, and the change in matrix properties is insignificant. Therefore, it is very important to consider all reservoir complexities to properly capture the spatial variability. The shown example assumes that the fracture plane crosses the entire reservoir, from top to bottom, but the same reasoning can be applied to layer-bound fractures or DFNs that may occur in certain zones. However, for DFNs, the flow would be dominated by the part of the reservoir section crossing the DFN. Figure 23 shows a comparison between the original case and a case that has a DFN in the top formation (Shallow-marine Tarbert).

Figure 23.

The Lorenz plot showing a comparison between the original (SPE10) case, and a case that has a DFN embedded in the top formation (shallow-marine Tarbert).

Once all heterogeneity indices are evaluated at every part of the reservoir, mapping in 2D is done to have a proper index that considers every feature or parameter. Figure 24 shows a comparison between a high-resolution heterogeneity index of the original SPE10 case with the original properties with another case where fractures are incorporated. While the SPE10 model is a simplified case, notable differences can be observed, and clear differences would be observed if a DFN is present, as opposed to a herewith-used fracture realization that cuts through the entire reservoir.

Figure 24.

Matrix-only heterogeneity mapping for entire reservoir (left) compared to matrix-and-fracture heterogeneity mapping for entire reservoir (right).

7. Observations and conclusions

The methodology described in this chapter provides an in-depth look at reservoir heterogeneity and presents a better way to account for parameters that impact fluid flow. The elaboration and quantitative analysis of the spatial heterogeneity index enhances the understanding of static reservoir properties and consequently the impact on dynamic field performance. The 2D mapping provides major insights regarding reservoir prolificity and can be utilized to derive field development plans with improved sweep. For example, the method does generate maps zone-by-zone or formation-by-formation and could even benchmark zones that do not appear consecutively. It offers the flexibility to generate a single heterogeneity map for the entire model or multiple heterogeneity maps, representing a specific zone or formation. In such instances, the mapping becomes representative of the vertical (z) direction and positioning, even conditioned structurally. The advanced applications of the Lorenz methodology, therefore, provide a common yardstick to standardize any property depending on its impact on fluid flow.

An immediate application of the method is the dynamic calibration (or history matching) of reservoir simulation models, where the automated workflow for the calculation of spatial heterogeneity index can be leveraged to render quantitative representation of reservoir connectivity. Such attributes can be, for example, utilized as an additional spatial conditioning parameter in optimization algorithms that minimize the misfit between the observed data and simulated response.