Abstract
We demonstrate how modelling decisions for a giant carbonate reservoir with a thick transition zone in the Middle East, most notably the approach to reservoir rock typing and modelling the initial fluid saturations, impact the hydrocarbon distributions and oil-in-place estimates in the reservoir. Rather than anchoring our model around a single base case with an upside and downside, we apply a comprehensive 3D multiple deterministic scenario workflow to compare-and-contrast how modelling decisions and geological uncertainties influence the volumetric estimates. We carry out a detailed analysis which shows that the variations in STOIIP estimates can be as high as 28% depending on the preferred modelling decision, which could potentially mask the impact of other geological uncertainties. These models were validated through repeated and randomised blind tests. We hence present a quantitative approach that helps us to assess if the static models are consistent in terms of the integration of geological and petrophysical data. Ultimately, the decision which of the different modelling options should be applied does not only influence STOIIP estimates, but also subsequent history matching & forecasts.
Carbonate reservoirs are estimated to hold 60% of the world's oil reserves and 40% of the world's natural gas reserves (Montaron 2008). However, they can be geologically complex in terms of sedimentology, diagenesis, wettability, and structural deformation and hence can exhibit substantial variations in reservoir properties within a scale of a few metres (e.g. Choquett and Pary 1970; Hollis et al. 2010; Burchette 2012; Agar and Hampson 2014; Chandra et al. 2015). A significant portion of these oil reserves are contained in the giant carbonate reservoirs of the Middle East, many of which have extensive transition zones with complex oil and water distributions that contribute to reserve estimates (Masalmeh 2000; Harrison and Jing 2001). The transition zone generally refers to the interval between the free-water level (where the capillary pressure is zero) or oil–water contact (at 100% water saturation) to the height where the water saturation reaches its irreducible level (Masaimeh et al. 2007; Spearing et al. 2014; Bera and Belhaj 2016). Understanding and modelling the distribution and associated uncertainties of the petrophysical properties in these reservoirs, including the transition zones, is therefore crucial for reliable STOIIP estimation, well planning, and reservoir performance forecasting.
In reservoir modelling, a single base case (also referred to as mid-case) model has been traditionally used. This base case is commonly considered as the anchor point for defining the high- and low-case end members or generating further stochastic realizations of the base case to address geological uncertainties. However, the base case can be biased by several factors; bias can be due to over/under-sampling of data (e.g. drilling wells in high-permeability zones or high-permeability zones being more fragile and hence under-represented in the core data), interpretation bias such as anchoring or prior experience, or the choice of specific modelling decisions (Baddeley et al. 2004; Bond et al. 2007; Bentley and Smith 2008; Hollis et al. 2010; Chellingsworth et al. 2011; Refsgaard et al. 2012; Arnold et al. 2013; Deveugle et al. 2014;Chandra et al. 2015; Wellmann and Caumon 2018; Costa Gomes et al. 2019). Moreover, it is also not uncommon to build reservoir models to support development decisions that have already been taken, instead of challenging the understanding of the reservoir by assessing different potential scenarios, which has been referred to as ‘modelling for comfort’ (Bentley 2016).
An alternative workflow that aims to limit the impact of bias in reservoir modelling and explore a more realistic range of geological uncertainties is to consider multiple deterministic scenarios. The fundamental idea of a multiple deterministic scenario workflow is to challenge our understanding of the reservoir by generating an ensemble of geological models that assesses the impact of uncertainty in the imperfect input data (Bentley and Smith 2008; Ringrose and Bentley 2015). Hence different geological scenarios are explored which consider, for example, different interpretations of the same data, different modelling approaches such as the decision which reservoir rock typing method to apply or which probabilistic algorithm to use when modelling the spatial distribution of petrophysical properties. New reservoir modelling approaches are being developed to facilitate the rapid design and testing of such scenarios (Jackson et al. 2015; Jacquemyn et al. 2021).
Figure 1 illustrates the conceptual differences between a workflow that aims to perturb a single base-case model and a workflow that considers multiple deterministic scenarios. While the latter workflow can yield a broader range of reservoir performance forecasts that need to be constrained as production data becomes available using the appropriate history matching techniques such as Bayesian methods (Christie et al. 2006), it is also more likely to include the actual performance of the reservoir behaviour.
Illustration of the conceptual differences between reservoir modelling workflows that consider perturbations around a single base case and workflows that consider multiple deterministic scenarios. Modified after Bentley (2016).
In the context of using multiple deterministic scenarios to quantify the impact of geological uncertainties, reservoir rock typing (RRT) schemes are a key concept that needs to be included. RRT comprises an essential element for any reservoir modelling workflow because it is used to propagate variability of petrophysical properties beyond wellbores and impacts volumetric estimates and the predicted dynamic behaviour of the reservoir (Lucia 1983; Amaefule et al. 1993). RRT can be defined as a way of categorizing and grouping rocks with similar geological or/and petrophysical or/and dynamic characteristics. These categories can be based on small-scale (e.g. pore-scale) data or large-scale (e.g. log data or even seismic data) (Ghanbarian et al. 2019). In most carbonate reservoirs, small-scale variations in petrophysical properties usually exhibit weak spatial correlations, while large-scale heterogeneities tend to be spatially controlled by stratigraphy and diagenesis (Jennings et al. 2002).
RRT for carbonate reservoirs is particularly challenging because petrophysical properties are not only controlled by the sedimentary fabric but also diagenetic and structural overprints that cross facies boundaries (Gomes et al. 2008; Rebelle and Lalanne 2014; Skalinski and Kenter 2015). Traditional RRT methods for carbonate reservoirs focus on depositional features (e.g. depositional environment, facies, or diagenesis). One of the first examples of such a classification was presented by Dunham (1962), who proposed a simple method used to describe rock textures based on depositional settings while ignoring diagenesis. Another RRT approach is based on Lucia's method (Lucia 1983, 1995, 1999), which focuses on the concept of pore-size distributions in non-vuggy carbonates. Lucia's method identifies three rock quality classes (grainstone, grain-dominated packstone and mud-dominated) based on petrophysical properties (e.g. porosity, permeability and particle size groups). While this method is a convenient approach, it usually breaks down for more complex porous structures such as carbonates containing abundant micro-porosity (Rebelle and Lalanne 2014). Flow zone indicators (FZIs) (Amaefule et al. 1993) or global hydraulic elements (GHEs) (Corbett and Potter 2004) define petrophysical groups using the porosity and the permeability data of the reservoir. Although both methods are easy to apply, they do not link back to rock textures, facies, or diagenesis as they are purely based on petrophysical data. As a result, different depositional environments could end up being grouped into one rock type and, hence, could introduce additional uncertainties in the reservoir model, for example when the depositional environments that are grouped into a single rock type have different relative permeability and/or capillary pressure curves (Gomes et al. 2008).
Additional RRT approaches for carbonate reservoirs can be based on the pore type characteristics obtained from capillary pressure data, pore throat distributions, or nuclear magnetic resonance data. Winland introduced an empirical correlation for defining RRT groups; the so-called Winland R35 method uses the mercury injection capillary pressure curves and the pore throat radius corresponding to the 35th percentile of the mercury saturation to classify the RRT groups (Kolodzie 1980; Pittman (1992) suggested altering the empirical equation of Winland to improve the correlation by presenting regression equations for pore throat radii corresponding to mercury saturations ranging from 10% to 75%. Hollis et al. (2010) incorporated depositional lithofacies, pore network properties, and relative permeability measurements to define RRT groups. Recently, Tawfik et al. (2019) proposed the ternary rock typing plot, which uses the interrelationship between porosity, permeability, and irreducible water saturation for defining RRT groups.
Closely linked to RRT is the modelling of initial fluid saturation distribution because RRT allows the assignment of saturation height function for each rock type which influences the 3D saturation model. The way initial fluid saturations are modelled ultimately impacts oil-in-place estimates and the dynamic behaviour of the reservoir model, i.e. the quality of the history match and reliability of the reservoir performance forecast (e.g. AlAmeri et al. 2020). The initial hydrocarbon saturation can be modelled using a variety of techniques. For example, fluid saturations can be distributed statistically between wells using the water saturations from the well-log as anchor points (Helle and Bhatt 2002; Kumar et al. 2002; Cuddy et al. (1993) suggest calculating the water saturation as the product of porosity and log-derived water saturation based on the height above the free water level. While this approach is simple, fast, and unaffected by rock type, its application to reservoirs with a thick transition zone has disadvantages because the approach does account for rock quality changes beyond porosity or the calibrating wellbore interval or dynamic capillarity effects, and is biased towards the curve fitting of the log-derived water saturation (Harrison and Jing 2001).
Drainage capillary pressure curves from special core analysis (SCAL) data can also be used for defining saturation height functions. Leverett's J-function approach attempts to normalize all capillary pressure data into a universal curve as a function of water saturation (Leverett 1941). This approach is widely used and simple; however, the averaging inherent to Leverett J-functions produces poorer results when the permeability range is large. Skelt and Harrison (1995) suggest tuning an empirical correlation to fit capillary pressure curves for each rock type, which can be difficult and time-consuming. More recently, dedicated software packages have been developed to model the saturation distributions using a combination of the Leverett J-functions obtained from SCAL data and the Leverett J-function obtained from well log saturation (O'Meara 2019). This method adjusts saturations at the wellbore and in the model, and both are returning an intermediate solution with some changes to both sets of initial assumptions.
Considering the uncertainties and inherent advantages and disadvantages of the different RRT and saturation modelling approaches discussed above, the aim of this study is to evaluate and quantify the impact of different RRT and saturation modelling decisions in a multi-deterministic scenario workflow when estimating hydrocarbons in place and reservoir connectivity. We hypothesize that the choice of RRT and saturation modelling approach is a major uncertainty when calculating fluid volumes and ultimately influences preferential flow paths in giant carbonate reservoirs. More specifically, we analyse if the choice of RRT and saturation modelling method could mask other geological uncertainties in carbonate reservoirs with extensive transition zones. We therefore apply a reservoir modelling workflow that combines multiple deterministic scenarios and multiple stochastic realizations in a tight and giant carbonate reservoir from the Middle East with the specific objectives to:
Evaluate the impact on STOIIP estimates when using multiple deterministic scenarios compared to multiple stochastic realizations.
Explore how RRT and petrophysical modelling approaches alter the permeability model and reservoir connectivity.
Investigate how different saturation modelling methods impact the initial hydrocarbon distribution, particularly in the transition zone.
This paper is structured as follows. In the next section, we provide a brief overview of the field reservoir and available reservoir data used in this study. We then describe the general workflow methodology, starting with structural modelling, followed by property modelling and saturation modelling. This section is followed by model comparisons and validations. Finally, we evaluate and compare the results for each modelling step and discuss their impact on hydrocarbon estimates.
Field description and available data
Field X is a giant heterogeneous carbonate reservoir in the Middle East with a low-relief anticlinal structure. It is a late Cretaceous limestone reservoir deposited in a carbonate ramp environment. The reservoir is characterized by poor reservoir quality with an average porosity of 19% at the crest and 12% at the flank, and permeability ranging from 3 to 9 millidarcies at the crest and one millidarcy at the flank. The reservoir's thickness is around 45 metres, with a significant transition zone reaching up to 25 metres in thickness. The proprietary data for this study were collected from more than 100 wells and include conventional well logs, routine core analysis (RCA), special core analysis (SCAL), and image logs. In addition, data from more than ten well tests, including fall-off and build-up tests, and a 3D pre-stack depth migration (PSDM) cube were available. Field X has a good quality and reliable dataset, except for the seismic cube. The quality of the seismic data is poor and has low vertical resolution caused by the complex surface topography as the field is located in a coastal area covering land Sabkha and shallow and deep marine environments.
Methodology
Using the available data for Field X, we applied a deterministic multi-scenario reservoir modelling workflow that considers a hierarchy of uncertainties (Fig. 2). Throughout the paper the term ‘scenario’ refers to the part of the reservoir model design where we generate geologically consistent models using different interpretations of the available input data and making different choices as to which modelling approaches to apply (i.e. different RRT and saturation height function). We use the term ‘realization’ to refer to equiprobable versions of the same scenario that were generated by changing the seed number in the stochastic part of the modelling process.
Overview of the deterministic multi-scenario reservoir modelling workflow for Field X.
We started the process by analysing the structural uncertainty associated with the ambiguity inherent in the seismic data acquisition and the construction of velocity models. Hence, the uncertainties in the 3D seismic cube arise when identifying the seismic horizon corresponding to the top of the reservoir. This horizon was used to build the first scenario for the structural model. The second structural scenario was then generated deterministically using the same seismic horizon but with a revised velocity model. The third and fourth structural scenarios were generated using an error map that consisted of an ensemble of 1200 stochastically generated structural models. The error map was calculated based on the difference between the actual and the predicted depth of the drilled wells. Four different structural model scenarios were generated this way. Porosity model realizations were generated by Sequential Gaussian Simulation for each structural model scenario.
Next, three different RRT approaches, namely, porosity–permeability cloud transforms, global hydraulic elements, and Lucia pore fabric classes were applied to each structural model scenario to create the RRT models and corresponding permeability distributions. This modelling step resulted in 12 different permeability models. For each permeability model, the water saturation was modelled using three different approaches, namely Leverett J-functions from core data (Leverett 1941), a combination of Leverett J-functions from the core and log data (O'Meara 2019), and the Skelt and Harrison method, which involves fine-tuning an empirical correlation to fit capillary pressure curves for each rock type (Skelt and Harrison 1995). The combination of permeability and saturation models resulted in 36 reservoir models. Finally, the robustness of all 36 reservoir models was analysed using cross-validation, blind testing, and through comparison with well-test data. Each step of the workflow will be described in more detail below.
Structural modelling
Four structural models of Field X were created to represent different geological scenarios. These models were generated based on the uncertainty in the depth analysis of the top horizon in the 3D seismic cube, arising, as noted above, from the uncertainties in the velocity model itself, as well as the poor quality of the seismic data caused by the complex topography and the combination of different types of equipment during the data collection process. The first two structural model scenarios were generated deterministically based on two different velocity models (old and revised velocity models); they are called structural Scenario A and B, respectively.
Additionally, we created an uncertainty envelope for the top horizon of structural Scenario A in order to capture uncertainties in the depth of individual reservoir tops away from the seismic well ties. This envelope was estimated stochastically by calculating the standard deviation between the predicted and actual well tops (well markers) for the reservoir top of structural Scenario A. In total, the envelope contains 1200 stochastically generated top horizons. Two of these 1200 models have been selected to quantify the effect on STOIIP estimates; these are referred to as structural Scenarios A1 and A2, and correspond to the P10 and P90 end members, respectively (Fig. 3).
The cross-section for structural Scenario A displaying 1200 stochastically generated reservoir tops, two of which were selected to create structural Scenarios A1 and A2. The thin vertical and horizontal lines are the well trajectories.
We performed the model zonation after reviewing all available data based on the vertical heterogeneity and lateral continuity of petrophysical properties. To begin, we analysed the core data heterogeneity using the Lorenz plot and Lorenz coefficients. Then, over 100 wells were correlated across the field and the results strongly suggest that the reservoir is layer-cake; neutron logs, density logs, and gamma-ray logs were used for this correlation. Following these analyses, the reservoir was subdivided into sixteen stratigraphic subzones. Then, the layering scheme was designed based on the heterogeneity and pore volume in the subzones, as well as the amount of oil present within each subzone. The thickness of the layers ranges from 1–6 ft, and generally a coarser resolution was used in low pore volume, more homogeneous, and high-water saturation intervals. The lower limit of the resolution was selected based on the vertical resolution of the logs, which was c. 1 ft during the logging-while-drilling data acquisition. Table 1 summarizes the static model of Field.
Summary of the static model of Field X
Property modelling
Multiple stochastic realizations were used to develop porosity models for each scenario, i.e. equiprobable models were generated by Sequential Gaussian Simulation based on a variogram analysis for each zone. The variogram has been calculated using the upscaled porosity data to describe the spatial continuity of the porosity data. The Sequential Gaussian Simulation is a commonly used algorithm for creating property models of continuous variables such as porosity. By changing the seed number and variograms parameters, 4200 different porosity model realizations were generated. Furthermore, we randomly changed the input well data, holding back 37% of the wells containing porosity data for later cross-validation and blind testing of each model. The resulting 3D porosity models were designed to honour the porosity logs at the wells. It is worth noting that the upscaled porosity log in each realization was calculated from the density-neutron logs, which were calibrated and validated with core data.
Permeability models were then generated by correlating porosity to permeability using different RRT approaches, namely
Porosity–permeability cloud transforms (cloud);
Global hydraulic elements (GHE);
Lucia's pore fabric classes.
The porosity–permeability cloud transforms were generated based on a porosity–permeability relationship that captures the spread of the data points for each individual subzone (Fig. 4). Porosity–permeability cloud transforms have the advantage of preserving the uncertainty in the relationship between porosity–permeability datasets compared to traditional linear porosity–permeability correlations (Waite et al. 2004). In the cloud transform, a specific porosity value may correspond to a range of permeability values that can reach up to two orders of magnitude (Waite et al. 2004). These porosity–permeability clouds use a bivariate distribution method that was based on the porosity model as a secondary variable. The porosity model was used as a secondary variable to control the spatial permeability distribution because the porosity model has a lower uncertainty compared to other petrophysical properties considering that data from 100 well logs were available to generate the porosity model.
Example porosity–permeability cloud transforms for two different subzones in Field X. The black dots correspond to the measured porosity and permeability from the core data, and the coloured dots represent the recalculated porosity and permeability data for each subzone.
Ten GHEs were generating using FZIs that are defined as (Amaefule et al. 1993)
Porosity and permeability values from core data are plotted on the ten predefined Global Hydraulic Elements (GHE) for Field X.
Permeability was modelled separately for each subzone based on the assigned GHE using the following relation which ensures that the spatial correlation of the porosity field is inherited by the permeability field (Corbett and Potter 2004)
Porosity–permeability values from core data for one of the subzones in Field X plotted against Lucia's rock fabric classes (left), and all porosity–permeability data for Field X plotted against Lucia's rock fabric classes (right).
In addition, we also used the Winland R35 method (Kolodzie 1980) to support the saturation modelling discussed below but did not use it for permeability modelling because there was considerable uncertainty in assigning Winland R35 groups to uncored wells and intervals. Using a subset of the data from 20 cored wells and clustering nearly 300 corresponding mercury injection capillary pressure data, we identified five rock types based on Winland R35 (Fig. 7 and Table 2). The water saturation for each rock type was then calculated using Skelt and Harrison's (1995) method. The rock types based in Winland R35 are as follows:
RRT1 is considered to have the best reservoir rock quality because it is a coarse-grained peloidal skeletal grainstone that suggests a high energy depositional environment (inner ramp).
RRT2 and RRT3 consist of fine- to medium-grained peloidal skeletal wackestone/packstone with a coarser fraction of coarse-grained to granule-size skeletal fragments (e.g. orbitolinids), deposited in a moderate to high-energy, shallow-subtidal environment.
RRT4 is a heterogeneous packstone texture, and the grain composition suggests a moderate to high energy environment (shallow subtidal); it represents moderate to poor reservoir quality.
RRT5 is the lowest quality reservoir rock with moderate to low porosity and low permeability below 1 md.
Five reservoir rock types for Field X based on the Winland R35 method. The porosity–permeability plot (left) with capillary pressure curves (middle) shows the RRT scheme's clustering effect for each rock type.
Map view showing the distribution of the five Winland R35 rock types in one of the reservoir layers.
Boundaries for each reservoir rock type based on the Winland R35 classification
The Winland 35 boundaries for each rock type are given by
Saturation modelling
Three different approaches were employed for modelling the initial fluid saturations, i.e. Leverett J-functions from core data (Leverett 1941), a combination of Leverett J-functions from the core and log data (O'Meara 2019), and the Skelt and Harrison (1995) method.
The Leverett J-function equation is given by (Leverett 1941)
Single J-function curve fitted to 18 capillary pressure measurements obtained using the porous plate method.
In the second approach, we used a commercial software package (Geo2Flow) to generate a model of the initial water saturation by combining the water saturation from log data and capillary pressure curves from the core data (Fig. 10). The specific steps to combine these data sources are as follows:
Quality control and then import all the capillary pressure curves.
Convert all capillary pressure curves to dimensionless Leverett J-functions.
Specify the Leverett J-function's average, maximum and minimum curves.
Convert the well-log water saturation from the depth domain to the Leverett J-function's domain.
Create synthetic permeability logs for each well using the 3D permeability model.
Define the error ranges and the cumulative distribution functions of these errors for the following parameters: irreducible water saturation, permeability, and Archie exponent data.
Calculate synthetic saturation logs which honour the log-derived and core-data derived Leverett J-function saturation by modifying permeability values, irreducible water saturation and Archie exponents within specified error ranges.
Produce the 3D saturation model using geostatistical interpolations (kriging) of the synthetic saturation log.
Water saturation modelling in Field X using a combination of the Leverett J-functions derived from core data and log data.
In the last approach, saturation height functions were generated for the five Winland R35 reservoir rock types using Skelt and Harrison (1995) (Fig. 11). The Skelt and Harrison (1995) correlation is given by
Saturation height functions for the five Winland 35 reservoir rock types based on the Skelt and Harrison (1995) method.
For simplicity, we ignore three other uncertainties in the reservoir data. First, we assume a constant formation volume factor of 1.39 RB/STB, but note that the formation volume factor varies between 1.3 and 1.75 due to variations in fluid composition across the reservoir. Second, we assume a constant free water level, which varies between ±10 ft across the reservoir. Third, we assume a constant gas–oil contact which varies between ±5 ft across the reservoir. These simplifications allow us better to analyse the impact of modelling decisions on STOIIP estimates.
Model comparison and validation
The porosity models were compared and validated using three different approaches. First, we quantified the impact of the stochastic nature of the Sequential Gaussian Simulation on the model's pore volume by running the algorithm multiple times for different seed numbers; the resulting porosity models are referred to as PHIE_SEED_1, PHIE_SEED_2, and PHIE_SEED_3. Secondly, we carried out repeated blind tests by excluding 37% of the wells containing saturation data from the input data and using the remainder of the data to model the porosity distribution. Porosity logs predicted by the model for the excluded well data were then compared to these available data for validation. This process was repeated randomly several times in order to remove any bias during the selection of the wells. The resulting porosity models are referred to as PHIE_Blind_1, PHIE_Blind_2, and PHIE_Blind_3. Lastly, we re-ran the Sequential Gaussian Simulation multiple times for model PHIE_Blind_1 to quantify the impact of the stochastic nature of this algorithm again; the resulting porosity models are referred to as PHIE_Blind_1_SEED_1, PHIE_Blind_1_SEED_2, and PHIE_Blind_1_SEED_3.
The permeability models were validated in two ways. First, we compared permeability model predictions to the core permeability data for each well where such data exist. Secondly, we calculated the geometric average for the permeability model over the drainage radius for each well where a well-test was carried out and compared the resulting value to the permeability obtained from the well-test. We note that such comparisons are only semi-quantitative because permeability data from different scales are compared. However, this exercise still provides an indication if there are significant inconsistencies between the permeability model and the available static and dynamic reservoir data.
The saturation models were compared to the water saturations obtained from the well logs. The water saturation log was calculated using Archie's equation and validated with core samples using the Dean-Stark method (Dean and Stark 1920); therefore, water saturations from well logs are assumed to be correct and used as a reference to estimate the error. The saturation models were also validated through repeated and randomised blind tests, where 30% of the wells containing saturation data were held back for validation purposes, and the other 70% of the wells provided the input data. Saturation logs predicted by the model for the excluded wells were then compared to these available data for validation. As with the porosity models, this process was repeated in order to remove any bias during the selection of the wells.
Two error norms were selected to calculate the model accuracy and the average error value of the porosity and saturation models for all the scenarios. We used the mean absolute percentage error ε given by
Results and discussion
This section presents the results obtained from the structural modelling, porosity modelling, RRT and permeability modelling, and from modelling the saturation distributions. The impact on the estimation of the reservoir and fluid volumes will be discussed.
Structural model
The variation in pore volumes as a result of the geological uncertainties encountered during the structural modelling is summarized in Figure 12. The calculation of the pore volume was limited to the reservoir volume between the gas-oil contact and the free water level. Recall that the fluid contacts are assumed to be fixed to eliminate any impact of the uncertainty in the precise location of the fluid contact. There is only a ±5% variation in the pore volume calculated across all 1200 different structural models where the reservoir top was varied stochastically between the wells (Fig. 3). This relatively minor variation is due to the fact that all these stochastically generated reservoir tops were anchored to a single base case model (Scenario A). In contrast, creating a different deterministic structural scenario while honouring the same well data (Scenario B) yields an increase of 12% in pore volume. This result may be explained by the fact that while much of the reservoir top is constrained by the over 100 wells, the area of high pore volume mainly lies in the undeveloped part of the field where the seismic quality is poor and less well data is available, which leaves more room for different interpretations and therefore results in larger changes in pore volume.
Variability in pore volume estimates resulting from the different structural models. Note that Scenario B represents a single different deterministic scenario while the distribution reflects the uncertainty arising from modelling the location of the reservoir top in Scenario A away from the well.
Porosity model
The different stochastic realizations for the porosity models show that the seed number has a minimal impact on pore volume calculation; the total pore volume variability is less than 0.5% when changing the seed number (Fig. 13). We point out that although changing the seed number hardly changes the total pore volume, it changes the spatial distribution of the pore volumes, which then influences the permeability distribution, and hence impacts the dynamic simulation and forecasting. During the blind tests where only 63% the wells containing saturation data were used as input, the variability in total pore volumes was less than 1%; this variability was caused by using random combinations of different wells as the input data. When the seed number was varied together with the blind test, the variability in the total pore volumes increased to c. 2%. Figure 13 summarizes the variability in pore volumes for the different porosity modelling approaches. Table 3 summarizes the error metrics for the porosity model. The absolute error ε can be as high as 8% but the RMSE was always below c. 0.02 (i.e. two porosity units). This relatively small variability in total pore volume was expected because of the high well density. The wells are also relatively evenly distributed across the developed part of the field, which provides more reliable input data for porosity modelling and hence reduces uncertainty. Additionally, the variance of the porosity data is relatively consistent across scales, ranging from 0.003 for the core data to 0.002 for the well logs used across the entire model. Based on the observation that the total pore volume hardly changes, a single porosity model for each structural scenario was taken forward for the subsequent permeability and saturation modelling steps. Figure 14 gives an example of a 3D porosity model, and Figure 15 shows a comparison of the seven porosity models across Well F by changing the seed numbers (i.e. for different stochastic realizations) and selecting different input wells during blind testing.
Comparison of the relative changes in pore volumes for the porosity model of Scenario A using all the wells and the porosity models created from different stochastic realizations by changing the seed number (PHIE_SEED), the porosity models created for the blind testing (PHIE_Blind), and the porosity models created for a combination of blind testing and different seed numbers (PHIE_Blind_SEED).
Example porosity model for structural Scenario A.
Representative comparison of the porosity model predictions for a single well (Well F) in structural Scenario A using different stochastic realizations (PHIE_SEED 1, 2 and 3) and during the blind testing (PHIE_Blind 1, 2 and 3).
Summary of the error metrics for the porosity model in structural Scenario A
RRT and permeability model
Figure 16 shows the resulting permeability models that were created from the porosity model using the three different RRT approaches, i.e. the cloud transfer, the GHE method, and Lucia's method. Recall that more than 20 wells were available with relevant porosity and permeability measurements and well-log data for this modelling step. It is apparent and not surprising that the different RRT approaches lead to different permeability distributions, which ultimately will impact reservoir performance predictions during dynamic modelling. Table 4 and Figure 17 show the variation and averages of the 3D permeability data. The permeability distribution using the porosity–permeability cloud transforms exhibits significant vertical heterogeneity. Recall that the cloud transforms resulted in 16 unique porosity–permeability relationships compared to the five petrotype groups in GHE and the three Lucia classes.
Three illustrative permeability models were obtained from three different RRT approaches.
Permeability histograms for each of the three RRT approaches for each structural model scenario.
Arithmetic permeability average k¯ for each of the three RRT approaches for each structural model scenario
Figure 18 compares the geometrically averaged model permeability to the permeability obtained from the well tests that were carried out for 12 wells. These well tests data were used for qualitative comparisons only as both, the permeability averaging and the well-test data are prone to uncertainties. We note that more consistent modelling approaches exist that aim to use well-test data to calibrate reservoir models (Hamdi et al. 2014; Egya et al. 2021). The arithmetic permeability average was calculated based on the estimated reservoir thickness and the drainage radius, which contribute to the flow during the well test for each well. The drainage radius has been defined as a circular shape (vertical well) or ellipsoid shape (horizontal well) around the wellbore and varies between 250 meters to 1000 meters. The correlation between the well test permeability and the averaged model permeability in the vicinity of the wells is variable, and there are some noticeable outliers where the averaged permeability overpredicts the well-test permeability or vice versa. These outliers occur at wells where the uncertainty in the input parameters for the well test (e.g. estimating the actual reservoir thickness contributing to the flow, not allowing enough testing time to reach the linear flow region in a horizontal well, or poor calibration for the pressure sensors). Furthermore, some wells intersect regions where the reservoir contains disconnected fractures, which can lead to higher well-test permeabilities, and the matrix permeability must be adjusted to account for the presence of such fractures (AlAmeri et al. 2020).
Comparison between the averaged modelled permeability from the three different RRT methods to the observed well test permeability obtained at 12 wells. Note that the largest outliers are caused by the presence of fractures or uncertainties in the well-test analysis itself.
Saturation modelling
Thirty-six different saturation models were generated for every possible combination of the structural models, RRT, and saturation modelling approaches. The resulting nine different saturation models for structural Scenario A are shown in Figure 19. Figure 20 shows a comparison of these nine saturations models at three wells to illustrate the variability in the saturation distribution. The saturation model generated using the core J-functions is dependent on porosity, permeability, and the height above free water level (HAFWL). In contrast, the log- and core-derived Leverett J-functions that are used to populate the saturation model are intended to honour the saturation log data at the well location and depend on porosity, permeability, and HAFWL observed at the wells. Finally, the Skelt-Harrison saturation model is dependent on the HAFWL and Winland R35 RRT groups, which define the porosity and permeability of the RRT. Figure 21 depicts the results of an example blind test in which 30% of the wells containing saturation data were not used when constructing the reservoir model but instead used for analysing how accurately the different permeability and saturation models predict the properties at these wells. Note that the saturation logs that were generated using the combination of the log- and core-data derived Leverett J-functions yield a better prediction during the blind test than the saturation logs generated using the Skelt-Harrison and core-data derived Leverett J-functions.
Nine different saturation models for structural Scenario A arising from the different possible configurations of saturation models and RRT methods.
Comparison of the nine saturation models at three illustrative wells (A, B and C). The background colour with the black line represents the saturation measured in the logs. The first track depicts the saturation model derived from three different RRT approaches employing the core- and log-derived J-functions which always honour the saturation logs and hence there is no difference between the saturation model and well log. The second track depicts the saturation model derived from three different RRT approaches employing the Skelt-Harrison method. The third track depicts the saturation model derived from three different RRT approaches employing the core-data derived J-functions.
Results of an example blind test in which 30% of the wells containing saturation data were excluded from the saturation modelling process for validation purposes. The resulting tracks compare the saturation models from the log- and core-derived Leverett J-functions and the Skelt-Harrison approach to the available log data at three illustrative wells that were held back from the saturation modelling process.
The saturation models that combine the log- and core-derived Leverett J-functions show an excellent agreement between the modelled saturations at the wells and the saturation data from the logs (Fig. 22). This approach also shows a more consistent saturation profile in the transition zone compared to the saturation data at the wells that were modelled using core-derived Leverett J-functions or the Skelt-Harrison method (Fig. 23).
Comparison of water saturation data observed at the well logs and the predicted water saturation models for structural Scenario A.
Comparison of water saturation data observed at the well logs and the predicted water saturation for the blind-tested wells using a saturation model based on core- and log-derived Leverett J-functions. The underlying permeability model is based on the cloud transform.
Saturation models using a combination of log- and core-data derived Leverett functions have the lowest errors, with RMSE and ε both below 0.04 and 2%, respectively. The reason for this excellent agreement is that this particular method honours the saturation logs by altering the permeability, irreducible water saturation, and Archie exponent based on the predefined error range from the available data. For permeability, these input ranges vary by approximately one order of magnitude while the irreducible water saturation ranges vary between 0.05 to 0.15, and the Archie exponents ranges vary between 1.8–2.2. In contrast, RMSE and ε are above 0.1 and 25%, respectively, for the Skelt-Harrison method and the core-data derived Leverett-J-functions which are based on independent data sources from the well logs with their own underlying assumptions. Therefore, a combination of log- and core-data derived Leverett functions is considered the most robust saturation modelling approach. Further blind tests were carried out using the porosity–permeability cloud transform with a combination of log- and core-data derived Leverett saturation methods; again, only 70% of the wells containing saturation data were used as input for the saturation model and data for the remaining 30% of the wells was used for validation. In this blind test, the RMSE remains below 0.07, but ε increases to 7%, which is still significantly lower than the RMSE and ε values for the Skelt-Harrison method and core-data derived J-functions. Detailed error metrics for the different scenarios and models are provided in the appendix.
In all other structural scenarios, the saturation models that use a core-data derived Leverett J-functions show the highest error because the broad range of capillary pressures curves inherent to this reservoir is averaged into a single J-function for each rock type. This observation is also apparent in the relative lack of sensitivity of the Skelt-Harrison method to the choice of the RRT approach in comparison to the core-data derived Leverett J-function.
STOIIP estimates
Figure 19 shows the different saturation models which will, of course, impact STOIIP estimates. Figure 24 summarizes the relative changes in STOIIP estimates across all 36 different saturation models. STOIIP estimates increased by up to 27% (structural Scenario B) or decreased by up to 7% (structural Scenario A2) compared to the operator's current STOIIP estimate. The absolute difference between the maximum and minimum STOIIP estimates range between 25–28% in each structural scenario.
Variations in STOIIP estimates across the 36 different saturation models, which consider the uncertainty in the structural model, the RRT approach, and the saturation modelling method. These changes in STOIIP estimate are with reference to the current operator estimate.
Structural Scenario A contains one model where the STOIIP estimates decrease by 7% and another model where the STOIIP estimate increases by 18% compared to the operator estimate. The differences in these two models, that is a change in 25% in STOIIP estimate, are only caused by the different choices made during the modelling workflow and the inherent assumptions of a given workflow, not in the underlying reservoir data. Similarly, structural Scenarios A1 and A2 predict a minimum STOIIP that is c. 7% below the operator estimate and a maximum STOIIP that is 19% and 21%, respectively, above the operator estimate. In contrast, for structural Scenario B, the minimum STOIIP is only 1% below the operator estimate while the maximum STOIIP is 28% above the operator estimate. In all cases, this significant range in STOIIP estimates is only caused by the different choices made during the reservoir modelling workflow.
The impact of choosing a certain RRT approach on the variability in the STOIIP estimates is comparatively small when using GHE or the Lucia classes. Depending on the chosen saturation modelling approach, the absolute difference between the maximum and minimum STOIIP estimate is between 24–14% for GHE, and between 25–9% for the Lucia classes for the different structural scenarios. This difference is in stark contrast to the RRT approach where porosity–permeability cloud transforms are used; here the absolute difference between the maximum and minimum STOIIP estimate is 28% depending on the choice of saturation model.
The variability in STOIIP estimates is small for the saturation modelling approach that combines log- and core-data derived Leverett J-functions; the absolute difference between the maximum and minimum STOIIP estimates is between 3–4% depending on the choice of RRT method. This difference increases to c. 7% for saturation models based on the Skelt-Harrison method and 23% for saturation models based on core-data derived Leverett J-functions.
The saturation model using a porosity–permeability cloud transform and the core-data derived Leverett J-functions yields the lowest STOIIP estimates in all four structural scenarios. This is because the porosity–permeability cloud transforms have the lowest permeability ranges in each grid cell compared to the other permeability modelling methods and use a single average J-function curve, which shifts the water saturation to higher values and decreases STOIIP estimates. On the other hand, the same porosity–permeability cloud transforms yield the highest STOIIP estimates when the saturation model is based on a combination of log- and core-data derived Leverett J-functions, that could be explained as log- and core-data derived Leverett J-functions method used a range of Leverett J-function instead of the single curve and alter the permeability value based on a predefined error ranges to match well log saturation as discussed earlier (Fig. 10).
This analysis confirms the hypothesis which is central to this study, namely that the choice of RRT and saturation modelling approach is a first-order uncertainty when calculating fluid volumes in giant carbonate reservoirs and could potentially mask other geological uncertainties. Importantly, the choice of RRT and saturation modelling approach also directly influence preferential flow paths and the quality of the history match, as shown by AlAmeri et al. (2020) for the same reservoir.
Conclusions
We showed that using a multi-deterministic scenario workflow allowed us to explore a broader range of the uncertainties related to the geological assumptions and modelling approaches (i.e. commonly used RRT methods and saturation height function models) in a giant carbonate reservoir. The uncertainties inherent to the decisions made during the reservoir modelling workflow and its inherent assumptions have a significant impact on the STOIIP estimates, with the absolute difference between the lowest and highest STOIIP estimates being as high as 28%. Such uncertainties could potentially mask other geological uncertainties inherent to the reservoir description, such as defining alternative correlations, considering the presence of fractures, etc. The uncertainties related to reservoir modelling decisions are also expected to have a significant influence on the resulting reservoir dynamics as indicated by AlAmeri et al. (2020).
For this particular reservoir, we found that the combination of core- and log-derived Leverett J-functions cause the lowest error in the saturation model and generally yield the highest STOIIP estimates, independent of the chosen RRT approach. If the saturation model is based on core-data derived Leverett J-functions, the saturation model always has the highest error and yields the lowest STOIIP estimates when combined with porosity–permeability cloud transforms. Permeability models generated using the GHE methods show the smaller variations in STOIIP estimates compared to permeability models that are based on Lucia's classes or porosity–permeability cloud transforms.
While it is not clear if our observations are universally applicable to other carbonate reservoirs and the root causes for the variability in STOIIP estimates need to be analysed case by case, our work highlights the need to design an appropriate a multi-deterministic modelling scenario that carefully considers both, uncertainties inherent to the geological data and reservoir description, and the uncertainties inherent to the choice of reservoir modelling workflows. Finally, we note that the uncertainty analysis discussed here can be readily linked to other uncertainty quantification workflows such as the Monte-Carlo based workflows of Athens and Caers (2019) and Yin et al. (2020), which provide a probabilistic framework to quantify uncertainties in individual model scenarios but have previously mainly been applied to a single base case, not a range of different model scenarios developed here.
Acknowledgements
The authors would like to thank Abu Dhabi National Oil Company (ADNOC) for providing access to data and the scholarship for Mohamed. Sebastian acknowledges support from Energi Simulation for his Chair in Carbonate Reservoir Simulation. We are also grateful to Schlumberger for access to Petrel and to Vicki and Dan O'Meara for access to Geo2Flow.
Author contributions
MA: conceptualization (lead), formal analysis (lead), investigation (lead), methodology (lead), validation (lead), writing – original draft (lead), writing – review & editing (equal); SG: conceptualization (equal), methodology (supporting), supervision (supporting), writing – review & editing (equal); PC: writing – review & editing (supporting)
Funding
This work was funded by the Abu Dhabi National Oil Company.
Data availability
The datasets generated during and/or analysed during the current study are not publicly available due to the confidentiality agreement with the funder.
Appendix
Summary of the error metrics for the saturation models of structural Scenario A
Summary of the error metrics for the saturation models for structural Scenario B
Summary of the error metrics for the saturation models for structural Scenario A1
Summary of the error metrics for the saturation models for structural Scenario A2
- © 2022 The Author(s). Published by The Geological Society of London for GSL and EAGE
This is an Open Access article distributed under the terms of the Creative Commons Attribution 4.0 License (http://creativecommons.org/licenses/by/4.0/)
References
Impact of modelling decisions and rock typing schemes on oil in place estimates in a giant carbonate reservoir in the Middle East
