Forecasting Fossil Fuel Production Through Curve-Fitting Models: An Evaluation of the Hubbert Model

Several methodologies have been developed to represent and forecast oil production curves for oilfields, regions, and countries. These methodologies are primarily categorized into three distinct approaches: economic, geophysical, and hybrid. The economic approach relies on factors like oil prices, costs, regulations, and other variables to elucidate oil production and supply. Conversely, the geophysical approach typically employs curve-fitting models, among which the Hubbert model of oil depletion is widely recognized [1].

Hubbert postulated that the production rate from a finite resource-producing region ascends to a peak before declining as the resources gradually deplete. His model, which popularized the concept of a symmetrical, 'bell-shaped' production cycle (a 'Hubbert' curve), aligns remarkably well with the historical experience of the USA [2].

Hubbert's model, which uses a logistic growth curve to represent cumulative production and the first derivative of the logistic curve for the production cycle, is most accurate in natural domains undisturbed by political or significant economic interference. It is particularly suitable for regions and areas with a multitude of independent oilfields [3].

The original, symmetric, and monocyclic logistic growth curve model introduced by Hubbert in 1959 has undergone various modifications and enhancements, spawning related curves such as asymmetric and multicyclic Hubbert, and Gaussian [4, 5].

Several elements of Hubbert's model have been altered, including the monocyclic character of the curve. The model hinges on the assumption of a single exploration cycle, with subsequent cycles either not occurring or having a negligible impact on total resource production due to their relative insignificance. However, for many oilfields, regions, or countries, a multi-cycle model provides a better fit for the production data [6, 7].

The symmetric form of Hubbert's model has also been scrutinized and improved. Critics argue that the symmetric form neglects the interplay of geological, technical, and economic factors, and it is impractical to assume that the same pattern will apply across a broad range of oilfield production cycles. Hallock et al. [8] proposed a modified version of the bell-shaped curve, peaking at 60% of ultimate production, as opposed to the conventional 50% peak.

Linear models have also proved to be realistic and reliable in some contexts; for instance, the production in the United States between 1945-2000 was better fitted by a linear production profile than a bell-shaped curve.

An alternative simple model, the exponential model, was utilized by Hubbert in his own paper. Data analysis revealed a 2% exponential growth in global oil production, followed by a decline of 10% per annum [9].

Hubbert's model has been adapted to represent the production of rare co-minerals such as zinc and indium. A related model, termed the Copula-Hubbert, was tested and selected as the best fit. Differences in peak values have been observed between the major minerals and co-minerals when compared to the Hubbert model [10].

Maggio and Cacciola [11, 12] employed the multiple Hubbert model (two cycles) to predict future trends in global oil production. By considering three potential scenarios for oil reserves, they were able to forecast peak oil production of crude oil and NGL (natural gas liquid).

A version of the multicyclic Hubbert model was used to forecast OPEC crude oil production. The multi-cyclic model was applied using a nonlinear least-squares numerical computation technique, with an initial guess for the three unknown parameters for each production cycle: k - the coefficient, Pmax - the maximum production (bbl/year, tons/year), and Tmax - the year of maximum production [13].

Brandt [14] evaluated six curve models utilizing 139 sets of oil production data from various regions and countries. The models consisted of three symmetric three-parameter models (Hubbert, linear, and exponential) and three asymmetric four-parameter models (asymmetric Hubbert, asymmetric linear, and asymmetrical exponential). After examining oil production data from regions in the USA and globally, Brandt concluded that oil production curves were significantly asymmetric in one direction, with a median rate of increase at 7.8%, while the median rate of decline was 2.6%.

Saraiva et al. [15] projected Brazil’s oil production curves per different URR scenarios (P90, P50, P10) by employing a modified multi-Hubbert model. The findings suggested that, excluding recent discoveries in pre-salt layers, Brazil’s peak oil should range from 2.37 Mbbl/day (2015), 3.33 Mbbl/day (2022), to 6.59 Mbbl/day (2035), contingent on URR scenarios. The model's accuracy, in relation to the data from 1954 to 2012, resulted in a relative standard deviation of less than 2.5%.

Regarding studies on Albania, a country with over 100 years of history in oil and gas exploration and production, substantial scientific research has been conducted and numerous articles have been published by both Albanian and international scholars, predominantly by geologists and engineers specializing in oil and gas exploration and production. These experts were primarily affiliated with the Institute of Oil and Gas in Fier, the Institute of Technology in Patos, and later at the Faculty of Geology-Mining, Albanian Geological Services, and other related institutions [16, 17].

Most research and studies are primarily concentrated in the realm of geological and geophysical research, where significant successes have been achieved [18]. In contrast, there are fewer studies in the field of oil reserves assessment, historical performance analysis, and future production forecasting based on historical data analysis. Mathematical models and geostatistical methods, or daily production optimization methods, are largely absent, indicating limited collaboration between specialists from different disciplines within the Albanian oil industry. A persistent issue has been, and continues to be, the lack of publicly available data regarding the annual production of individual oilfields, regions, or national production [19].

The online data from various institutions are fragmented, making it challenging to get a comprehensive view of the history of each oilfield, particularly the most significant ones like Marinza, Patos, Kucova, Cakran, and more.

The completion of data for the entire period of oil production, 1929-2019, was achieved by sourcing from various national and local sites as well as data from international institutions and private oil companies operating in several locations in Albania.

As a result, a comprehensive picture of historical oil production in Albania since 1929 has been compiled. This data is invaluable for engineers and oil specialists for analysing production performance over the years, evaluating and explaining changes, and identifying the best model to represent the production curve.

In this paper, we conduct an extensive analysis of Albanian oil production data from the initial years of exploration in 1929 up to recent years. The data is utilized to identify the most suitable models to represent and forecast future production. These models are reliable and can be used for other oilfields and regions, both old and new.

Table 2. Oil production in Albania, 1929-2019

3.1 The Hubbert model

Hubbert forecasted the future US oil production by fitting a curve to historical data on annual US production and projecting this forward in time under the assumptions that:

• production must eventually decline,
• the area under the curve must equal the URR of the USA—or the amount of oil that is both technically possible and economically feasible to extract over the full production cycle.

In his first presentation, in 1949, Hubbert did not produce a formula. In 1956, he produced the logistic formula for his bell curve model as shown in Figure 3.

3.png

Figure 3. Hubbert’s bell curve forecasting a peak in US oil production between 1965 and 1970

The key features of Hubbert’s logistic model are:

Cumulative production is modelled by a logistic function. Yearly production is modelled as the first derivative of the logistic function.

The production profile is symmetric (i.e., maximum production occurs when the resource is half-depleted and its functional form is symmetric on both sides of the curve).

Production increases and decreases in a single cycle without multiple peaks.

For any production curve of a finite resource, the production rate must begin at zero, and then after passing through one or several maxima, it must decline again to zero.

The area under the production rate curve (above the time axis), between the start time and the end time, is the URR amount (equivalent to EUR-Estimated Ultimate Recovery) [28].

Hubbert’s model of oil production is defined mathematically as follows:

$Q(t)=\frac{Q_{\infty}}{1+e^{-a(t-t m)}}\qquad, \quad q(t)=\frac{d Q(t)}{d t}=\frac{a Q_{\infty} \cdot e^{-a(t-t m)}}{\left(1+e^{-a(t-t m)}\quad\,\right)^2}$              (1)

where,

Q(t)-Cumulative production (tons, barrels- bbl),

q(t)-Production rate (tons/year, barrels- bbl/year),

$Q_{\infty}$-URR (Ultimate Recovery Reserves) (tons, barrels-bbl),

a- a coefficient that defines the “steepness” of the cumulative production curve,

$t m$-The time when cumulative production reaches one-half of the $Q_{\infty}$=URR, (year).

The multi-cycle model, introduced by Laherrere has been proven to fit the oil production data of several regions. In practice, the additional cycles can also result from efficient improvements in exploration and production methods, which can increase oil production significantly [29, 30].

The multi-cycle model is:

$Q(t)=\sum_{i=1}^k Q(t)_i=\sum_{i=1}^k\left\{\frac{2 Q_{\max }}{1+\cosh (b(t-t m))}\quad\right\}i$              (2)

where,

k-Total number of logistic curves that fit the production data;

Q_max- peak production rate (tons, barrels) and;

t_m- Peak production time, respectively for each cycle (year).

3.2 The Gompertz model

The Gompertz function is another useful model that can successfully fit the world’s oil production from the year 2009 onward, resulting in an asymmetrical production profile [31]. It was developed by Gompertz in 1825, as a modification of the symmetric logistic function to an asymmetrical growth curve.

Gompertz and Logistic models generate very similar curves. But when Y is low, the Gompertz model grows more quickly than the logistic model. Conversely, when Y is large, the Gompertz model grows more slowly than the Logistic model. The cumulative and production rate functions are:

$Q(t)=P_{\max } \cdot e^{-e^{-k(t-T p)}}$              (3)

$q(t)=\frac{d Q}{d t}=P_{\max } \cdot k \cdot e^{-k\left(t-T_p\right)} \cdot e^{-e^{-k(t-T p)}}$              (4)

where,

Q(t)- cumulative oil production, (tons, barrels- bbl);

q(t)- production rate, (tons/year, barrels- bbl/year);

P_max- maximum production rate, (tons/year, barrels- bbl/year);

T_p-year of peak production rate (year).

3.3 The Gaussian model

Another mathematical model, which may fit best in several cases, is the normal distribution, or the Gaussian model, because of its similar curve shape compared to Hubbert’s logistic model [32].

The Gaussian model works best with numerous independent wells and few regulatory constraints. The central limit theorem (CLT) is the justification for applying the Gaussian curve in statistical applications, as in this case. Based on oil production data in the US, it is concluded that over 20,000 oil producers in the USA are acting at random, "leading to a Gaussian curve", and the normal distribution model is the best-fit model [33].

The probability density function for a normal distribution is given by the formula:

$p(x)=\frac{1}{\sigma \sqrt{2 \pi}} \cdot e^{-\frac{(x-\mu)^2}{2 \sigma^2}},-\infty<x<\infty$              (5)

where,

μ- the mean, (the peak oil of the bell curve, time- year);

σ- the standard deviation (a parameter of the width of the bell curve);

If the three oil production rate parameters $Q_{\infty}, t_m, s$, are included in the original mathematical formula (5) then an adaptive Gaussian symmetrical model is produced:

$q(t)=\frac{Q_{\infty}}{s \sqrt{2 \pi}} \cdot e^{-\frac{\left(t-t_m\right)\,^2}{2 s^2}}$              (6)

where,

$Q_{\infty}$- the Ultimate Oil Recovery, (tons, barrels- bbl);

t_m- time, year, at the maximum production rate, (time- year);

s- the width parameter of the bell curve.

The asymmetric model based on a Gaussian curve is given by the formula:

$q(t)=s_{d e c}-\frac{s_{d e c}-s_{i n c}}{1+e^{k\left(t-t_m\right)}}$              (7)

$Q(t)=q_{\max } \cdot e^{-\frac{\left(t-t_m\right)^2}{2(q(t))^2}}$              (8)

where,

Q(t)- the cumulative production rate in the year t, (tons/year, barrels- bbl/year);

$q_{\max }$- maximum production rate, (tons, barrels- bbl);

t_m- time of peak oil production rate, - (year);

q(t)- the sigmoid function that changes the standard deviation near t=t_m;

s_inc- the standard deviation of the pre-peak production curve;

s_dec- the standard deviation of the post-peak production curve.

3.4 DCA (Decline Curve Analyses) models

DCA is another method of predicting the future of oil and gas production rates. It is used to describe and analyse declining production rates and predict the future performance of oilfield production [34].

The production rate of oil and gas decreases over time, for several reasons, mainly due to the loss of reservoir pressure as well as the decrease in the amount of oil in the reservoir.

The DCA method is based on fitting a line to performance history and assuming that the same trend will continue going forward.

3.4.1 Exponential decline

The exponential decline represents the case when the production rate declines by the same percentage each period.

The equation is:

$q=q_0(1-d)^t$              (9)

where,

q=production rate at the time t;

q_o=initial producing rate at the time to;

d=decline rate per period, (year, month);

t=time at which the calculation of q is required;

The cumulative production to a future time t is:

$Q(t)=\frac{q o\left(1-(1-d)^t\right)}{-\ln (1-d)}$              (10)

3.4.2 Hyperbolic decline

Production rates and cumulative functions are:

$q(t)=\frac{q_0}{(1+b d t)^{\frac{1}{b}}}$              (11)

$Q(t)=\frac{q_o{ }^b}{d \cdot(1-b)}\left[q_0{ }^{1-b}-q^{1-b}\right]$              (12)

where,

q=production rate at time t;

q_o=initial production rate;

b=factor;

d=initial hyperbolic decline rate;

t=time;

3.4.3 Harmonic decline

When b=1, the equations will show the harmonic decline function. Production rate and cumulative functions are:

$q(t)=\frac{q_0}{1+d . t}$              (13)

$Q(t)=\frac{q_o}{d} * \ln (1+d t)$              (14)