Abstract

We have developed an analytic formulation for quick and more accurate volumetric estimations of subsurface resource potential. Our formulation is conceptually based on a structurally conformable model built deterministically using known and interpreted reservoir properties from wells, such as net-to-gross, porosity, and hydrocarbon saturation, along with oil-water contact or lowest known oil depth and interpreted seismic top-of-pay depth horizon. We have evaluated an important function, the hydrocarbon pore capacity (HPC), which is a product of net-to-gross, porosity, and hydrocarbon saturation. HPC reflects the heterogeneity of key reservoir rock and fluid properties, particularly in the vertical direction. We calculated the total hydrocarbon pore volume directly from HPC and the top-of-pay horizon, without the explicit need for building a geologic model first. Our efficient solution form can preserve the vertical resolution of wireline logging with transparent parameterization and the least amount of averaging and upscaling. We also provided additional formulation for incorporating different HPC regimes for cases of multiple existing wells in the reservoir. We demonstrate the practical application of the formulation with a data example from a Cretaceous carbonate reservoir in the southern Gulf of Mexico offshore and with comparisons to other common approaches. In the application example, we determine that the most common approach of using single-average-values can underestimate the reserve upside by as much as 30%, whereas the stochastic modeling approach provided improved estimates when simulating porosity with a lognormal distribution and preserving the net-to-gross log in its original vertical resolution.

Introduction

Volumetric estimation is a critical step in petroleum exploration and development. It involves multiple disciplines in geoscience and reservoir engineering. In particular, there are three key input parameters closely related to the volumetric estimation of subsurface resources: net-to-gross (NTG), porosity, and hydrocarbon saturation. Historically and quite often in early stages of exploration, single average values of these parameters are used for bracketing possible drilling and development outcomes and forecasting hydrocarbon production for financial planning purposes. With advancements in geoscience technology, more reliable measurements and interpretation of the subsurface have become possible (e.g., Zeng, 2015; Wu and Hale, 2016; Wang et al., 2018). The subsurface can be evaluated more accurately, in part at least around the wellbore, and completely with help from contemporary seismic imaging and interpretation. The common approach of single-average-values can be further refined to quantitatively examine the impact of spatial variations of these input parameters to the subsurface volumetric estimation process and their integral effect on the estimation of reservoir resource potential.

Historically, volumetric estimations are based commonly on estimated averaged values of area, thickness, NTG, porosity, and hydrocarbon saturation (Craft and Hawkins, 1959; Tearpock and Bischke, 2003; Demirmen, 2007). However, numerous studies have reported serious pitfalls with the use of averages in scientific and economic decision making (e.g., Welsh et al., 1988; Claerbout, 1992; Savage et al., 2009). The key conclusion of these studies and many others is that conducting calculations based on averages can create erroneous and misleading outcomes. There are mathematical proofs and real case studies confirming that there can be serious pitfalls in making decisions or drawing conclusions based on averages of individual variables, even though the averages may have been thought of as our friends for “best guesses.”

This study presents an analytic formulation for subsurface volumetric estimation. The calculation of hydrocarbon pore volume (HPV), which is defined as a product of the abovementioned five quantities, is considered. In this manner, a deterministic and quantitative approach for calculating a product function of NTG, porosity, and hydrocarbon saturation, which is defined as the hydrocarbon pore capacity (HPC), is proposed. HPC accounts for the spatial variations of key reservoir rock and fluid properties, particularly in the vertical direction. By integrating HPC in 3D, guided by an interpreted top-of-pay horizon, the heterogeneity of reservoir rock properties, hence the resource potential, can be accounted for more efficiently and accurately. The formulation provides a direct analytic solution for HPV given HPC and top-of-pay horizon.

Method

A well penetrating a reservoir pay interval was assumed, with htop as the top-of-pay and hOWC as the oil-water contact (OWC) in vertical depth below the datum. The OWC depth can be replaced with that of the base-of-pay or the lowest known oil (LKO) hLKO, if the OWC is not encountered at the well. This study further assumed that the hydrocarbon type is crude oil. The well has available petrophysically interpreted reservoir property logs, among others, with r as NTG, φ as porosity, and (1Sw) as hydrocarbon saturation, where Sw refers to the interpreted water saturation. These three quantities are typically derived from integrated well-log and core analysis. For the purpose of this study, these functions were assumed to have zero values above the top of pay and below the OWC. A scalar function q is defined to represent the HPC at the well with a subscript “0” as follows:  
q0(z)={r(z)φ(z)(1Sw(z))if  z[htop,hOWC];0if  z[htop,hOWC].
(1)
Conceptually, this function is best described as, on a per unit volume basis, the ultimate in situ amount of hydrocarbons at the wireline logging scale. It is noteworthy that the unit of HPC, the q-function, is m3/m3, that is, the HPV per unit volume of rock. This is related to the fact that it is the product of NTG, porosity, and hydrocarbon saturation, which are all volume fractional quantities having values between 0 and 1. Hence, the values of HPC are also between 0 and 1. In the imperial unit system and practically, HPC has the unit of reservoir barrels per acre-foot (R.B./Ac. Ft.), which is calculated by multiplying an approximate conversion constant of “7758.” Note that R.B. is different from the more often used S.T.B., that is, stock-tank barrels, by a factor called the formation volume factor B0 in the petroleum engineering community (Craft and Hawkins, 1959).

An interpreted top-of-pay depth structure horizon h(x,y), which is typically derived from surface seismic, was also assumed. Here, h refers to the depth along the z-axis (assumed pointing downward vertically) and x and y are the horizontal Cartesian coordinates, typically lined up with the seismic acquisition grid’s inlines and crosslines. Note that h(x,y) is a single-z valued surface. The well location is denoted as (x0, y0), assuming a vertically straight hole, for convenience of the discussion.

With the assumption of structure conformance, the HPC log of the well can then be “walked” up and down the structure along the surface h(x,y), above the OWC (Tearpock and Bischke, 2003). Pictorially, this is similar to the process of “hanging” the HPC function of the well q0 from the horizon surface h(x,y). The q-function can now be populated laterally throughout the seismically interpreted reservoir extent, defined as follows:  
q(h(x,y),z)=q0(zh(x,y)+htop)H(hOWCz),
(2)
where H(hOWCz) is the Heaviside step function or unit step function (Abramowitz and Stegun, 1972), which is used to ensure zero values below the OWC.
Another important quantity is the HPV per unit area of rock (m3/m2, or R.B./Ac. in the imperial unit system), which is hereby defined as the hydrocarbon pore efficiency (HPE). Arps et al. (1967) report that this quantity tends to be proportional to the reservoir recovery efficiency. HPE is calculated at every (x,y) location by vertically integrating the q-function between the top-of-pay horizon and the OWC:  
Q(h(x,y))=h(x,y)hOWCq(h(x,y),z)dz.
(3)
With the assumptions of structural conformance and the single-z horizon surface, HPE is uniquely defined by the depth of the horizon and should have the same value along each horizon contour line. It should be noted that the Q-function has a maximum value of  
Qmax=htophOWCq0(z)dz,
(4)
which is an integration over the entire pay-column encountered at the well. The Q-function is a monotonically decreasing function of h(x,y) with a minimum value of zero at the OWC depth. This is because the q-function is always positive or zero. For the portion of the h(x,y) surface shallower than the htop, the HPE is a constant equal to Qmax. The total HPV can then be expressed as  
HPV=AQ(h(x,y))dxdy,
(5)
where area A refers to parts of the horizon h(x,y) above the OWC. The unit for HPV is m3 or R.B. in the imperial unit system. Equations 15 serve as steps for the calculation and quality control of HPC, HPE, and HPV.
Equation 5 can also be written in a more compact form by combining it with equations 14:  
HPV=A0hOWCh(x,y)q0(z+htop)dzdxdy.
(6)
Equation 6 represents a single-well direct solution, given the HPC and top-of-pay horizon.
For a reservoir with multiple-well penetrations, equation 6 can be applied to an individual well’s HPC regime. An HPC regime is defined as subarea(s), including the well’s take point, where its HPC is interpreted as applicable. The following is a generalization of equation 6 to N HPC regimes:  
HPV=i=1NA(i)0hOWCh(x,y)qi(z+htop(i))dzdxdy,
(7)
where A(i) is the subarea of the ith HPC regime, qi is the ith well’s HPC, and htop(i) is the corresponding well’s top-of-pay depth. Typically, the HPC regime can be interpreted based on geology, geophysics, and geometry. Local geology such as faults and stratigraphy can be treated as guides for inferring boundaries separating the wells. Geophysically, seismic waveform analysis (e.g., Chopra et al., 2019) can provide further support in defining HPC regimes. A workflow of applying an unsupervised seismic waveform classification was successfully tested to help in defining HPC regimes in a two-well situation. Mathematically, the wells can also be partitioned according to geometry, using approaches such as the Voronoi diagram (Aurenhammer and Klein, 2000).

Equations 6 and 7 are key results of this report; equations 15 represent their expanded version for implementing a practical workflow, with more detailed quality analysis and control.

The implementation of the formulation can be summarized as follows:

  1. 1)

    From the well, measure, derive, and interpret the three key input logs, NTG, porosity and hydrocarbon saturation, and OWC (LKO) depth.

  2. 2)

    Calculate the HPC log at the well, using equation 1.

  3. 3)

    From the processed surface seismic, interpret the top-of-pay depth horizon.

  4. 4)

    Compute the HPE horizon from steps 2 and 3, using equation 3.

  5. 5)

    Estimate the total resource potential by summing the live samples of HPE scaled by the seismic bin size, using equation 5.

  6. 6)

    Optional step: Derive the HPC attribute volume for further QC and analysis, using equation 2.

  7. 7)

    Optional step after step 3: Compute HPV directly using equations 6 and 7.

The workflow is demonstrated with a real data example in the next section.

Application example

The formulation and the proposed workflow were applied to data from two wells located in the Bay of Campeche, offshore southern Gulf of Mexico, in water depths of less than 100 m. The stratigraphic column consists of Mesozoic and Cenozoic strata (e.g., Bartolini and Ramos, 2009). Wells drilled within the area encountered approximately 4600 m of Tertiary, 300 m of Cretaceous, 230 m of Jurassic Tithonian, and 300 m of Jurassic Kimmeridgian sections, respectively. The Tertiary section can be best characterized as a clastic environment dominated by shale and sandstones. The Cretaceous section is dominated by a carbonate platform composed of clean limestones, mudstones, dolomitized limestones, and breccias. The Tithonian section is the primary source rock in the basin, with organic-rich carbonates, fine-grained clastics, shales, and mudstones as seen in the wells. The Kimmeridgian section is primarily dolomites consisting of oolites, grainstones, wackestones, and mudstones. Data shown below are from logged intervals in the Cretaceous section and consist primarily of limestones, with some dolomites and mudstones. Figure 1 presents summary plots of these two wells. The insets in these figures show the input logs, NTG, porosity and water saturation, and the derived HPC log (the q-function), varying with respect to vertical depth. The bold curves in both figures represent HPE (Q-function) derived through vertical integration of HPC (q-function).

Well 1 logged oil pay from approximately 4723 to 4955 m, whereas well 2 logged pay in the interval of approximately 4682–4930 m. The porosity and permeability of well 2 are primarily attributable to a fracture network with very little secondary porosity. It has an average porosity of approximately 4%. Well 1 logged a higher average porosity of approximately 8%, which can be best characterized as a highly fractured, brecciated, vuggy limestone. The dual-porosity nature of this rock creates enhanced permeability due to fractures in the reservoir and the vug-to-vug connectivity. The rock properties differences were also observed by a surface seismic waveform analysis. Gross oil pay thicknesses of wells 1 and 2 are 232 m (approximately 758 ft) and 248 m (approximately 814 ft), respectively. Figure 1 also shows that the HPE derived from well 1 is almost twice that of well 2 at almost all structural positions.

To further illustrate the workflow, a synthetic top-of-pay structure map with two fault blocks was generated (Figure 2). Well 1 is located on a down-thrown fault block to the northwest, whereas well 2 is located on an up-thrown fault block to the southeast. To aid the visualization of the reservoir model derived according to equation 2, an HPC attribute volume (i.e., the q-functions) was populated. Figure 3 shows cross-sectional views of the HPC volume in perpendicular directions, indicated as lines 1 and 2 in Figures 2 and 3. The final seismic depth-imaging grid spacing here is 12.5 m in the inline and crossline directions.

Figure 4 shows a map view of HPE, derived using equation 3 based on the HPC of the wells (Figure 1) and the horizon surface (Figure 2). The calculated HPVs for the up-thrown and down-thrown fault blocks are 4.84 and 14.25 MM m3, respectively (“MM” stands for millions). In the imperial unit, these values translate to approximately 30.4 and 89.6 MM R.B., respectively. For reference, the total enclosed areas above the LKO are also displayed in Figure 4 for the down-thrown and up-thrown fault blocks.

The total elapsed computation time in this application example, including HPCs, HPEs, and HPV, is less than 3 s on an HP Z820 workstation equipped with dual Intel Xeon CPUs of 3.10 GHz and 64 GB of RAM. This efficiency of near real time can provide a flexible interactive platform for HPV scenario testing.

Discussion

The proposed method, abbreviated as AVE, was compared with commonly used and established approaches, including (1) the use of single-average-values (SAV) for all three input logs, NTG, porosity, and hydrocarbon saturation, (2) stochastic modeling (SCM; see e.g., Peebles, 1987; Limpert et al., 2001) of all three inputs, and (3) the retention of the NTG log of the well, whereas stochastically modeling porosity and hydrocarbon saturation (hereby annotated as NTG_SCM). As shown below, geologic models built with these three approaches represent improved estimates, in increasing order, that fit the subsurface data within some prior statistical measures.

Figure 5 shows comparisons of AVE with SAV for wells 1 and 2, respectively. When averaged values of r, φ, and (1Sw) are used as input, the resulting HPC for well 1 is approximately 0.03. This constant HPC ignored the fact that the top half of well 1 has a much better capacity profile than that of the bottom; a similar observation is made for well 2. The resulting HPE function (the dashed line) significantly underestimates reservoir efficiency. The color fill between the two HPE functions reflects the portion of the reserve potential missed by the SAV approach.

Figure 6 shows further comparisons with SCM. Stochastic models can be realized based on probability distributions for r, φ, and (1Sw). Here normal, lognormal, and uniform distributions are used for these three input quantities, with key distribution parameters closely matching their respective actual data distributions from the two wells. Because of the improved estimation of rock properties, the resulting HPE function behaves better than that derived from SAV, but it still underestimates reservoir efficiency compared with AVE, illustrated by the gray fill. Nevertheless, the stochastic approach provides accurate estimates at the full pay columns of the wells.

The improvement of SCM over SAV can be mainly attributed to the fact that porosity exhibits a lognormal distribution. It should be noted that if a normal distribution is used in the calculation, SCM shows minimal improvement over SAV.

The estimation can be further improved by providing a good description of the vertical heterogeneity of reservoir properties. For example, the top half of well 1 has a better capacity profile than the bottom half. For the experiment, the actual input NTG function r was retained, and the remaining two parameters φ and (1Sw) were stochastically modeled, with lognormal and uniform distributions, respectively. Figure 7 shows comparisons of this approach with AVE. The HPC functions resembled the actual input, with better top rock properties. In addition, the corresponding HPEs were similar to those of AVE.

To summarize the impact of various approaches on HPV, Table 1 shows the estimated resource potentials using these approaches, and their percentage differences compared with that of AVE. For this particular example, SAV and SCM approaches underestimated resource potentials by approximately 30% and 20%, respectively, whereas NTG_SCM underestimated by approximately 7%.

In geologic model building typically facilitated by large commercial software, a process called upscaling is required to represent rock properties on length scales of the order of the seismic wavelength or a fraction of (e.g., Jones et al., 2009; Li and Hu, 2015; Lie, 2016). Upscaling is in essence a process of averaging. For the Cretaceous carbonate reservoir in the above example, a quarter of the seismic wavelength is estimated to be on the order of 100 m. Calculating HPV using a geologic model built based on upscaling the wireline logs to this scale or even just 5%–30% of it effectively reduces the total reserve upside, similar to the case of SAV shown in Figure 5, in which the reduction is approximately 30% for an equivalent of approximately 200 m upscaling. An upscaling from a fraction of a meter to several meters or several tens of meters would place the reduction between 30% and 0%.

The formulations in this study provide direct link of HPV to HPC and top-of-pay horizon. It can calculate the reserve upside without the explicit need of building a geologic model first. It can retain the original resolution of wireline logs. Contemporary geologic model building commonly requires upscaling of rock properties, which means a typical one to two orders of magnitude averaging on the input rock properties. Additionally, it demands more intense manpower and expertise in the model-building process, which many independent operators may not possess. The formulation was extended to cover cases with multiple wells penetrating the reservoir. HPC regimes can be defined with additional knowledge and information based on geology, geophysics, and geometry. The preliminary experiments suggest that seismic waveform classification can be successfully applied to support the definition of HPC regimes.

Conceptually, the method can be applied to pays in carbonate and clastic reservoirs, as long as the assumption of structure conformance is within practical tolerance. In more complex geologic settings, one can bring in seismic application elements, such as waveform classification and prestack inversion, as quantitative guides to the HPC regimes (equation 7). Additionally, the HPC volume can be validated, at least partially, through an iterative process of inverse modeling with HPC-reflectivity volume benchmarked with the surface seismic images, automatically and through seismic interpretation.

Although the formulas are for crude oil as the hydrocarbon type, in an output form of HPV under in situ conditions, it can be equally applied to the scenario in which the hydrocarbon type is natural gas, in which case the OWC and LKO annotations become GWC and LKG, respectively, and with an output unit of m3, which is identical to that of the oil case. In the imperial unit system, the unit for in situ natural gas HPV is ft3, with a conversion constant of approximately 43,560 (e.g., Tearpock and Bischke, 2003).

The high efficiency of the formulation allows for quick evaluation of key reservoir rock properties and reservoir geometric parameters of resource potentials. The formulation can fit into an interactive workflow in an integrated geoscience interpretation and reservoir engineering environment. It can be deployed rapidly to iterate multiple reservoir models with different reservoir properties, such as LKO contacts or complex stratigraphic distributions of rock properties, without building separate geomodels for each instance. It can also be ultimately applied to the population of a geomodel into the 3D work environment, as used by more conventional methods. The time-saving and flexible features of this formulation can help in propagating a petrophysical model across a depth structure, with original resolution and parameter transparencies, hence providing an excellent means for further multidisciplinary integration.

Conclusion

This paper presents an analytic formulation for subsurface resource volumetric estimations by integrating data and interpretation from wells and surface seismic. The key parameter log from the well is HPC, which accounts for effects from the NTG, porosity, and hydrocarbon saturation logs. Vertical integration of the HPC log yields HPE that directly impacts the total volumetric estimation. With a seismically mapped top-of-pay horizon and an estimated depth of the OWC, the formulation can efficiently compute the HPV based on an implicit deterministic reservoir model with structural conformance. It also allows for outputting HPC as an attribute volume for additional visualization and quality control analysis.

The formulation was demonstrated with an application example from a carbonate reservoir, offshore the southern Gulf of Mexico. HPCs, HPEs, and HPVs can be efficiently calculated and analyzed in near real time. The ability of direct calculation of HPV from HPC and top-of-pay horizon provides a practical platform for interpretive procedures such as resource scenario testing. The proposed formulation was also compared with some other common approaches. The most commonly used method, single-average-values approach, is shown to underestimate the reserve upside, by as much as 30% in the data example. With progressively less averaging, or upscaling, this difference would decrease. The results show that the stochastic modeling approaches can provide better estimates when simulating the porosity behavior with a lognormal distribution, for this particular Cretaceous carbonate reservoir, yielding an underestimation of approximately 20%. The estimates for the stochastic modeling further improved when the wells’ NTG logs were preserved in their original vertical resolution, in which case the underestimation reached 7%.

Although the assumption of structural conformance can be challenged in real applications, the formulation’s ability to preserve wireline logging resolution can provide deterministic and transparent computation with the least amount of averaging. This in turn can provide more accurate subsurface volumetric estimation and potential reserve upside.

Acknowledgments

We deeply appreciate Fieldwood Energy for permission to publish these results. We thank Carlos Buenrostro for providing petrophysical interpretation of the two wells, especially the three key input logs, net-to-gross, porosity, and water saturation. We thank Eugene Wissinger for his constructive comments on our manuscript. We also thank Gary Janik and Scott Schmidt for their encouragements on our pursuit of technical excellence. We wish to thank all five reviewers and the editorial board for their comments and suggestions.

Data and materials availability

Data associated with this research are confidential and cannot be released.

Biographies and photographs of the authors are not available.

Freely available online through the SEG open-access option.