Henderson  Petrophysics

Limitations of Statistical Methods For Estimating Hydrocarbon Reserves

Abstract

Statistical (stochastic) methods are widely used to calculate hydrocarbon reserves. With modern software it is very quick and easy to generate probability distributions that describe uncertainty in the reserves calculations.

Without exercising extreme care, the statistical method can lead to serious errors in estimating reserves. Close examination of the data is essential to ensure appropriate use of parameter probability distribution functions (PDFs). Extensive documentation is also required that fully explains the derivation of these PDFs. It is also essential that any parameter dependence is fully documented - parameter independence must be demonstrated.

There should be no expectation that statistical methods provide either improved accuracy or better assessment of uncertainty. The statistical method is unlikely to have a rigorous audit trail and the methodology gives no information about the spatial distribution of the hydrocarbon reserves.

Introduction

In this article we investigate some limitations of statistical methods of estimating hydrocarbon reserves.

Statistical reserves modelling has been promoted as a means to provide a more accurate calculation and to better characterize uncertainty. In the following sections we will show that the probability distribution of calculated hydrocarbon reserves is dramatically affected by the methods and assumptions used to analyse and process the available reservoir data. There should be no expectation that statistical methods provide either improved accuracy or better assessment of uncertainty. The statistical method is unlikely to have a rigorous audit trail and the methodology gives no information about the spatial distribution of the hydrocarbon reserves.

Only at the time of final field abandonment are the hydrocarbon reserves known with certainty. Until that time there is always uncertainty in calculating the volumes of hydrocarbon that exist in reservoirs and the quantities that can be economically extracted.

Uncertainty in reserves estimates results from limitations of the input parameter data where the data is incomplete either because of the high cost of obtaining the data or the inability to access existing proprietary reservoir data. A related source of uncertainty is bias, resulting from inappropriate use of the available data. Bias is a predisposition of the interpreter and is not necessarily conscious or intentional.

Hydrocarbon reserves can be estimated using either deterministic or statistical methods. The deterministic method uses a single value, or a very limited range of values, for each input parameter in the reserves calculation. The values for each parameter are those that are deemed to be the most appropriate. Limitations of the deterministic method are the logistics in revising the estimates and the very limited scope for describing the uncertainty of the resulting reserves estimate. However, the deterministic method is likely to have a robust audit trail and an excellent visualisation of the spatial hydrocarbon distribution.

Statistical methods are very widely used to model the uncertainty is estimating oil and gas reserves, especially in the exploration and early appraisal stages. Each parameter in the reserves calculation is assigned a range of possible values with their associated probability. Monte Carlo simulation is used to repeatedly sample random values from the parameter probability distributions. These values are then used to calculate a reserves volume. The results of the Monte Carlo simulation are then sorted to yield a probability distribution for the reserves. The expectation is that the process of incorporating uncertainty functions will result in a correct analysis of the range of possible reserves.

Data for Stochastic Reserves Estimation

In this article we have used, for illustation and modelling purposes, data from;

Benmore, R., Cooper, M., Wells, B., "Stochastic Modelling for Reducing Risk in Prospect Evaluation", AAPG Annual Convention, Houston, Tx., March 1995.

The Benmore et al. paper describes the statistical analysis of expected reserves in a North Sea satellite development prospect. Since the area is a mature basin, there is extensive reservoir data. Table 1 and figure 1 show the reservoir parameters from the Benmore et al, 1995 paper.

Table 1: Drillable Prospect Reservoir Parameter Table

Parameter	Lowside	Most Likely	High Side
Gross Rock Volume(Acre-feet)	7600	70150	246400
Porosity (%)	14	18	24
Net to Gross(%)	26	40	58
Hydrocarbon Saturation(%)	15	48	78
Formation Volume Factor(RB/STB)	1.38	1.25	1.12
Recover Factor(%)	10	15	30

Figure 1: Reservoir Parameter Triangular Probability Distributions

Porosity Model

Figure 2 shows a histogram of the porosity data from the Benmore et al., 1995 paper on stochastic reserves modelling. The histogram shows that the data is very negatively skewed - it has a long 'tail' towards low porosity values. The average porosity value is 20.69 percent and the modal (most common) value is about 22 percent.

Figure 2: Histogram of porosity data

For their stochastic prospect evaluation, Benmore et al. modelled porosity using a nearly symmetrical triangular probability distribution function (Low = 14%, ML = 18%, High = 24%). The green dashed line in figure 3 shows that the cumulative probability distribution function (CDF) for this triangular distribution is a very poor match for the actual data (blue line in figure 3).

The orange line in figure 3, derived from a highly asymmetric triangular probability distribution (Low = 14%, ML = 22.2%, High = 24%) , matches the available data quite well.

Porosity data is very frequently misused in statistical reserves calculations. We are interested in modelling the uncertainty of the average value in the Gross Rock Volume. It is absolutely incorrect to use the range of values seen in a data set. The porosity average uncertainty function is expected to have a normal probability density function (bell-shaped curve) centered on the population mean and having a standard deviation (Standard Error) that is far less than the sample standard deviation. The red curve in figure 3 is a CDF showing the uncertainty of the average porosity value. This CDF has a mean of 20.69 percent (average of all of the available porosity data) and a standard deviation of 0.576 percent (the sample Standard Error).

Standard Error is calculated by dividing the sample standard deviation by the square root of the number of samples. As the number of samples increases the Standard Error (uncertainty of the average value) will decrease.

Figure 3: Porosity Cumulative Distribution Functions (CDF)

Hydrocarbon Saturation Model

The saturation-height functions in the upper portion of figure 4 show that there is a strong correlation between porosity and hydrocarbon saturation. The lower part of figure 4 shows two functions relating hydrocarbon pore-volume to porosity. The blue curve is derived from the saturation-height curves. The purple curve is derived using the porosity and hydrocarbon saturation triangular distribution data in table 1 and assuming that there is a correlation between porosity and hydrocarbon saturation. There is a very significant difference between these two functions.

Figure 4 : Saturation-height Functions Showing a Strong Porosity-Sh Dependence

Figure 5 shows four possible hydrocarbon saturation CDF curves. The black curve shows the hydrocarbon saturation CDF calculated using the triangular probability density function data in table 1. This curve appears to be a very poor match with the saturation-height derived Sh data.

The blue, red and green curves incorporate the Sh-porosity dependence derived from the saturation-height curves in figure 4.

The blue curve uses the porosity CDF calculated from the individual porosity measurements (data is sorted by increasing value and plotted against the corresponding cumulative probability).

The red curve is a reasonable match for the sample data. It was derived using the porosity triangular probability density data in table 1. The apparent good match with the data is misleading, since it was shown previously that the porosity triangular probability distribution data in table 1 is a very poor match for the available porosity data (blue curve in figure 3).

The green curve in figure 5 is an estimate of the probability distribution of the Sh average value. Porosity data was first used to calculate hydrocarbon saturation. The average and Standard Error of these estimates was used to generate the green curve. This should be the best model to describe uncertainty in the average Sh value in the Gross Rock Volume. However, it will be shown below that the porosity-Sh dependence must be specifically modelled in the Monte Carlo reserves estimate.

Figure 5 : Hydrocarbon Saturation Functions

Figure 6 shows hydrocarbon pore-volume (porosity times hydrocarbon saturation) plotted against porosity for three elevations above the free water level (FWL). These show that the average saturation decreases towards the FWL and that low porosity sands will contain little or no hydrocarbons.

The hydrocarbon reserves impact of the transition zone are enormous, in part because of the increasingly large GRV increments toward the FWL. Stochastic processes are probably not appropriate for modelling potential reserves in the transition zone.

Figure 6 : Decreasing Hydrocarbon Saturation towards the free water level.

Net To Gross Model

Figure 7 shows a crossplot of porosity and net to gross data. For this data the correlation between porosity and net to gross is very weak and is if no practical use for modelling hydrocarbon reserves.

Note: In some data sets there is a strong porosity-NTG correlation that should be modelled for Monte Carlo simulation.

The Net to Gross data shown in figure 7 is clearly derived from Petrophysical zone summaries (coupled porosity and NTG values). This Net to Gross data is of only very limited use for statistical prospect reserves calculations;
   - We have no way of knowing the location of the net reservoir within the stratigraphic interval of interest.
   - We do not know if the data is for contiguous reservoir units. Petrophysics summary tables rarely list the properties of non-net intervals.
   - We do not know the relationship between unit thickness and net to gross. It is quite possible that the high NTG values correspond to relatively thin stratigraphic intervals.

All of the above factors are absolutely critical for a robust modelling of drillable prospect reserves uncertainty. It is highly likely that non-reservoir (zero net-to-gross) intervals exist in the stratigraphic interval of interest. These intervals are rarely included in Petrophysics summary tables. It is also essential to know the vertical distribution of the reservoirs. Low (or high) NTG zones near the base of the interval of interest will have a very dramatic effect on the expected hydrocarbon reserves (net sands have a large area). Zone thickness information corresponding to the porosity and NTG data is also vital.

Figure 7 : Crossplot of Net to Gross and Porosity Data

Figure 8 shows three possible NTG Cumulative Probability Distribution Functions. The blue CDF is derived by sorting the available data. The green dashed CDF shows that the available NTG data is quite well described by the triangular distribution data in table 1. However, the probability distribution of the actual NTG values is of very limited use in reserves calculation. In fact it is both incorrect and misleading. We are most concerned with our ability to describe the uncertainty of the average value in the entire gross rock volume. The red curve in figure 8 describes the CDF for the sample average NTG. It is the cumulative probability distribution of a normal PDF with an average value of 42.5 percent and a standard error of 1.88 percent.

Figure 8 : Net to Gross Cumulative Distribution Functions

Recovery Factor

Rarely would it be appropriate to use a recovery factor PDF to estimate potential reserves in a multi-reservoir exploration prospect. Recovery factors are calculated using specific information on several reservoir factors, for example; in-place reserves contained in a modelled drainage area, expected reservoir continuity, drive mechanism, proximity of water or gas contacts, depletion strategy and individual well economics. In a following section we have used a very simple example to show the effects of modelling recovery factor as a dependent function of gross rock volume. However, the recovery factor PDF must also be linked to an actual or conceptual field development plan.

The Impact of Parameter Dependence on Stochastically Generated Reserves

The green curve in figure 9 shows the probability distribution of reserves calculated using a restricted Monte-Carlo simulation. The Gross Rock Volume, Net to Gross, Formation Volume Factor and Recovery Factor parameters were held constant at the most likely values in table 1. The Monte-Carlo simulation assumed that porosity and hydrocarbon saturation are independent variables.

The porosity and Sh CDF functions were sampled by finding the value corresponding to random numbers 0 <= u <= 1 with each parameter being sampled with a different random number. For each iteration the reserves are calculated using;

Reserves = Constant[GRV_ML x NTG_ML x 1/FVF_ML x RF_ML] x Porosity x Hydrocarbon Saturation

The Monte-Carlo simulation results are sorted and plotted by cumulative probability.

Figure 9 : Stochastic Simulation of Reserves Showing Impact of Porosity-Sh Dependence

The purple curve in figure 9 shows the results of modelling the Sh dependence on porosity. As with the previous Monte-Carlo simulation, the Gross Rock Volume, Net to Gross, Formation Volume Factor and Recovery Factor parameters were held constant at the most likely values in table 1.

The porosity CDF function was sampled by finding the value corresponding to random numbers 0 <= u <= 1. Hydrocarbon pore volume (PhiSh) was calculated as a function of porosity using the relationship derived from the saturation-height functions (blue curve in the lower portion of figure 4). For each iteration the reserves are calculated using;

Reserves = Constant[GRV_ML x NTG_ML x 1/FVF_ML x RF_ML] x Phish=f(porosity)

In the above simple example, there is about 10 percent difference in calculated reserves at the P50 level by incorporating a porosity-Sh dependence into the Monte-Carlo simulation.

Monte Carlo Hydrocarbon Reserves Modelling Results

Case 1: Triangular Probability Distributions with No Parameter Dependence

The black line in figure 10 shows the probability distribution of reserves calculated using the triangular probability distribution parameters in table 1 and figure 1. All of the parameters are assumed to be independent variables.

Case 2: Reserves Probability Distribution With Modelled Porosity-Sh Dependence

The green line in figure 10 shows the probability distribution of reserves calculated with hydrocarbon saturation being a dependent function of porosity. The porosity CDF is derived directly from the available data. The Sh-porosity dependence is derived from the saturation-height curves in figure 4. For each iteration of the Monte Carlo simulation the porosity CDF was randomly sampled. Hydrocarbon saturation (Sh) was then calculated using;

HPV = -0.4001 x Porosity ² + 1.5275 x Porosity - 1.27
Sh = HPV/Porosity

All of the other hydrocarbon reserves parameters in table 1 were treated as independent variables.

The reserves calculated using a porosity-dependent hydrocarbon saturation (green line) are far larger than the reserves calculated using independent variables (black line). For this data set, the difference results primarily because the porosity and saturation triangular distribution data in table 1 is a very poor match for the available data.

Figure 10 : Reserves Probability Distribution

Case 3: Porosity Probability Distributions with Modelled Porosity-Sh Dependence

The green curve in figure 11 shows the results of using the porosity CDF derived directly from the available data (blue curve in figure 3)

The purple curve in figure 11 shows the results of modelling the uncertainty of the average porosity in the gross rock volume. The available porosity data has an average of 20.69 percent and a standard deviation of 2.57 percent. The distribution of averages is therefore modelled as having an average of 20.69 percent (assume average = population mean) and a Standard Error of 0.575 percent.

As with Case 2, hydrocarbon saturation is modelled as a dependent function of porosity.

As expected, the resulting reserves probability distribution has lower variance than the results from Case 2.

Figure 11 : Monte Carlo Reserves Simulation

Case 4: Modelling the Results of Recover Factor Dependent On GRV

For the data described by Benmore et al. 1995 it is very likely that recovery factor is inversely related to gross rock volume. The GRV vs depth plot from Benmore et al. 1995 is shown in the top left panel of figure 12. The top right panel shows Recovery Factor plotted against GRV. The recovery factor data is from the table 1 triangular probability distributions. The bottom panel in figure 12 shows the impact of modelling RF as a GRV-dependent variable. As expected, there is a significant reduction in the reserves because of the low recovery factor interpreted in the large area towards the lowest spill point.

Figure 12 : Modelling the effects of a relationship between GRV and Recovery Factor

Figure 13 and table 2 show the results of the above Monte Carlo simulation models. It is appealing to conclude that inclusion of parameter dependence into the Monte Carlo simulation resulted in a robust evaluation of the expected reserves probability distribution for this satellite development prospect. This is incorrect because we have shown that the available data is both incomplete and provides no information on the spatial distribution of the reservoir parameters.

Figure 13 : Summary of Monte Carlo Reserves Simulation

Table 2: Summary of Stochastic Modelling

Model	P90 Lowside	P50 Most_Likely	P10 Upside	Spread (P90-P50)
Independent Triangular Distributions	1.40	3.84	8.32	6.92
Porosity CDF from dataset Sh as a porosity-dependent function using Saturation-height functions	2.90	7.54	15.91	13.01
Model as above but with Recovery Factor modelled as a function of GRV	2.77	5.61	9.53	6.76
Model as above but using the porosity average uncertainty function	3.32	5.94	9.72	5.95

Summary and Conclusions

With modern software it is both quick and easy to implement statistical (stochastic) methods for estimating hydrocarbon reserves. However, the method is not a robust and reliable method for modelling the uncertainty in the hydrocarbon reserves calculations.

Humans have a natural inclination to;
- Underestimate the complexity of systems (in this case the geological parameters required to calculate reserves),
- Be over confident in their ability to characterize uncertainty, and,
- Propose that extreme values have a possibility of being representative of the entire system (porosity, for example).

The vertical distribution of hydrocarbons in a reservoir is a complex function of porosity, grain size, permeability, fluid density/viscosity and distance from fluid contacts. The saturation-height functions shown in figure 4 are gross simplifications. Statistical (stochastic) methods are typically inadequate for modelling these complexities. For example, it would be extremely rare to find an example where a stochastic model adequately characterized the hydrocarbon-water transition zone. This is part of the hydrocarbon system with the largest gross rock volume. In low porosity/permeability reservoirs, where the transition zone can be extremely thick, failure to correctly model hydrocarbon saturation can lead to gross overestimation of reserves.

Without exercising extreme care the statistical method can lead to serious errors in estimating reserves. Close examination of the data is essential for a robust reserves assessment. Extensive documentation is required that fully explains the derivation of the input parameter probability distribution functions.

It is imperative that any parameter dependence is adequately modelled. It also must be proven that parameters assumed to be independent are not correlated with any other input parameter.

The Statistical method should address only the uncertainty of the expected (mean) value of the input parameter in the entire gross rock volume of interest. The CDF must define the uncertainty in the estimate of the average. This is not the same as the actual range of values observed in the reservoir.

It is crucial that the parameter probability distributions are appropriate for the rock volumes being evaluated. For example, it is completely incorrect to use the highest measured porosity value as having a realistic possibility in the gross rock volume being evaluated. An example would be the decision to use the highest fracture porosity measured in a hand or core sample as having a possibility of being realised in a fracture reservoir play.

The spatial distribution of reservoir parameters must be modelled and documented. Using tabulated reservoir summaries of porosity, net to gross and hydrocarbon saturation is absolutely useless without knowing the stratigraphic (vertical and areal) distribution of these properties.

Recovery factor probability distributions must be used with extreme caution in Monte Carlo reserves simulations. Simulator derived recovery factors that ignore economic flow rates or capital investment returns can be very misleading.

---

Appendix 1: Calculation of a CDF From Triangular Probability Distribution Data

Figure 14 shows how to calculate the cumulative probability distribution function (CDF) from the triangular probability density function (PDF).

• a is the minimum value,
• b is the maximum value,
• c is the most likely value, and,
• x is any value between a and b.
• CDF is the cumulative probability distribution function (CDF)

Figure 14 : Triangular Probability Density Function

---

Appendix 2: Porosity Population Modelling

Reservoir porosity data sets are often composed of sub-populations that each typically have normal probability density functions (the PDF's have a characteristic bell shape). The histogram in figure 2 indicates that the porosity data is multi-modal (the data consists of a composite of sub populations). Examination of porosity data and the derived histogram indicates that there are at least 3 overlapping sub-populations with the following properties;

Sub-Population	Fraction of Total Population (%)	Average (%)	Standard Deviation (%)
1	15	15.9	1.9
2	20	19.0	1.2
3	65	22.1	0.94

Figure 15 shows the normal probability distributions for these 3 sub-populations in the porosity data.

Figure 15 : Stochastic Simulation of Porosity CDF

The porosity population CDF can be can be stochastically simulated if the following parameters are known, or can be estimated;

• Proportion of each sub-population and,
• The porosity average and standard deviation for each of the sub-populations,

The stochastic simulation program steps are;

1. Generate a random number from 0 <= u <= 1,
2. If u <= the fraction of sub-population 1, then,
• Generate another random number 0 <= u <= 1
• Calculate the porosity value corresponding to the cumulative probability u of the normal curve with the sub-population 1 average and standard deviation
3. If u <= the sum of sub-populations 1 and 2, then,
• Generate another random number 0 <= u <= 1
• Calculate the porosity value corresponding to the cumulative probability u of the normal curve with the sub-population 2 average and standard deviation
4. If u > the sum of sub-populations 1 and 2, then,
• Generate another random number 0 <= u <= 1
• Calculate the porosity value corresponding to the cumulative probability u of the normal curve with the sub-population 3 average and standard deviation
5. Repeat steps 1 - 4 many times
6. Sort the resulting porosity values and calculate the cumulative probability with each iteration having the probability of 1 divided by the number of iterations.

The following Visual Basic code fragment shows how the porosity CDF can be simulated from the sub-population data.

i = 0
While (i < 5000)
  u = Rnd()
'Sub-population 1
   If (u <= .15) Then
      u = Rnd()
      Porosity = WorksheetFunction.NormInv(u, 0.159, 0.019)
'Sub-population 2
   ElseIf u <= .35 Then
      u = Rnd()
      Porosity = WorksheetFunction.NormInv(u, 0.19, 0.01163)
'Sub-population 3 (Highest porosity and most common) 
  Else
      u = Rnd()
      Porosity = WorksheetFunction.NormInv(u, 0.221, 0.0094)
   End If
i = i + 1
wend

The stochastically generated porosity CDF (red line in figure 15) matches very well with the available porosity data CDF (light blue line with data points in figure 15).

---

Appendix 3: Porosity - Probability Distribution of Sample Averages

Figure 16 shows that the histogram of average values for samples randomly selected from a highly skewed population. The expected sample mean is symmetrical around the population mean.

Figure 16 : Porosity Sample Probability Density Functions

---

 

How to Contact Us

Henderson Petrophysics offers specialised formation evaluation services to the upstream oil and gas business.

In additional to informative articles and useful downloads, our Website contains a lot of information on the services and expertise offered by Henderson Petrophysics.

Henderson Petrophysics Home Page

For more information on our services or our fee schedule you can send an e-mail (click on the icon below) or contact us using the following information;

Phone:    +61 7 3300 3980  (61 is the country code for Australia and 7 is the area code for Brisbane)
Mobile:    +61 407 653342
E-mail:    don@hendersonpetrophysics.com