A critique of Hubbert’s model for peak oil

In his seminal paper, prepared for a meeting of the American Petroleum Institute in Texas in 1956, Shell Oil geologist Marion King Hubbert attempted to answer the question, “how far along have we come on our way to complete exploration [of world oil fields]?” (Hubbert 1956, p. 2). The paper presented a resource model constructed to answer this question. Hubbert writes that,

as an essential part of our analysis, we can assume with complete assurance that the industrial exploitation of the fossil fuels will consist in the progressive exhaustion of an initially fixed supply to which there will be no significant additions. (Hubbert 1956, p. 4).

Consequently, the two fundamental questions that drive the analysis, both for Hubbert and again here in this critique, are as follows:

As per Morehouse (1997), we will use the term “ultimate resource” (UR) for the amount in question 1, and note that several synonyms are often used: “estimated ultimate recovery” (EUR), “ultimate recoverable resources”, and “ultimate recovery” (Morehouse 1997; Maggio and Cacciola 2012).

In 2005, to mark its 125th anniversary, Science published a special issue introduction titled “What Don’t We Know?” (Kennedy and Norman 2005). Among this list of the most profoundly vexing inquiries currently confounding the human race is an updated variant of the Hubbert questions: “What can replace cheap oil—and when?” (Kerr and Service 2005, p. 101). In this form of the question there appears to be an implicit assumption that progressive exhaustion of world oil reserves results in “expensive” (or even, “no more”) oil. Thus, this question leads back to the twin questions of ultimate reserve amount and prediction of future production levels, now nuanced by additional questions concerning: the cost and difficulty of extraction, the availability of cost-competitive energy alternatives, and society’s willingness and ability to curb its appetite for cheap oil. These latter questions, though increasingly relevant, are not addressed. Rather, we attempt to explore, update, and critique the answers to the two fundamental questions provided by the original Hubbert model from 60 years ago.

A first estimate of the UR of oil is to use the petroleum engineering profession’s estimate of “world oil reserve”, R(t), which is the amount of technically and economically recoverable oil from all known world oil fields at a particular time, t ≥ 0, measured in years. Such estimates are fraught with uncertainty, depend heavily on reliable geologic and engineering data (not always available), require expert and experienced interpretation of those data, and are hopefully untainted by geo-political interference. “Proven” reserves, based on the definition approved by the Society of Petroleum Engineers, are those claimed to have a reasonable certainty (statistically, at least 90% confidence) of being recoverable. Such proven reserves are known in the industry as “1P” (Society of Petroleum Engineers 2005).

In this paper we will use the consistent time series oil data as published annually by British Petroleum (BP) (BP 2017a, 2017b). These BP data, both for “reserve” and for “production”, purport to be inclusive regarding the type of oil referenced. According to BP,

Oil reserves include field condensate and natural gas liquids as well as crude oil. This inclusive approach helps to develop consistency with the oil production numbers published in the Review, which also include these categories of oil. Liquid hydrocarbon fuels from non-hydrocarbon sources, such as ethanol from corn or sugar or synthetic oil derived from natural gas (so-called GTL or gas-to-liquids), are not included in either the reserves or production series. (BP 2017c)

We will use BP data for the world oil production rate at time t, here denoted P(t) gigabarrels (Gbbl)/year. From the same source, R(t) consists of 37 data points, beginning in 1980 (t = 0) to the most recent values (t = 36) from 2016. Data for P(t) go back to 1965, yielding 52 data points. The cumulative amount of oil extracted, Q(t) Gbbl, since the discovery of oil in Pennsylvania in 1859, is then constructed by adding the yearly production amounts to an estimate from Gallagher (2011) for Q(0), herein labelled Q_o = Q(0), setting t = 0 to a convenient year.

When we use R(t) as a first estimate of UR, the simplest estimate of the “time to depletion” (i.e., an answer to the second posed question) is to compute the ratio, R(t)/P(t) years, which measures the number of years remaining given the current reserve level and present rate of extraction (production) until the reserve is completely exhausted. Measurements of this R/P ratio have been employed for decades, popping up in diverse places. For instance, in the environmental science text (Harte 1988), the computation based on 1980 data¹, reportedly extracted from Hubbert (1969), is as follows:

The Institution of Mechanical Engineers in the UK posted a 2014 article to their website, which has since been removed, titled “When will oil run out?”. Employing the R/P ratio, they stated:

There are an estimated 1.3 trillion barrels of proven oil reserve left in the world’s major fields, which at present rates of consumption will be sufficient to last 40 years. By 2040, production levels may be down to 15 million barrels per day—around 20% of what we currently consume. (Institution of Mechanical Engineers 2014)

Other journalists use the same information in a less alarmist fashion. For instance Maclean’s from February 2014:

According to BP’s Statistical Review of World Energy, global oil reserves at the end of 2012 were 1.7 trillion barrels. Given that the world consumes about 86 million barrels of crude oil per day, it would be easy to conclude we’ll run out of oil in 55 years, or sooner if we increase consumption. (Leach 2014)

The author then goes on to explain why, despite generally increasing consumption rates, the oil reserves amount has increased over the past 30 years rather than declined, concluding with a definitive prediction, “So, what’s the bottom line? Are we going to run out of oil? No”. (Leach 2014)

So what do the data show? Figures 1a and 1b show the most recent world oil reserve and production data from BP, with their regression lines in each case. Both reserve and production exhibit strongly increasing trends. The regression line for reserve has a very strong correlation coefficient of 0.9895, and reserves have been increasing at a rate of almost 30 Gbbl/year since 1980 (slope of 29.17). The correlation coefficient for production is also strong at 0.9564, and world oil production has approximately tripled since 1964.

Figure 2 shows the time history of the R/P ratio for the world with the regression line for the computed data.

Clearly the data show an upward trend in the R/P ratio as well; the regression line has a slope of 0.5004. Thus, the number of years of oil remaining has been increasing by half a year per year since 1980. Further, the growth in the R/P ratio has occurred despite steadily increasing production.² The correlation coefficient of 0.8932, corresponding to a coefficient of determination of r² = 0.7979 for the R/P ratio data, indicates that the growth trend is strong. As a consequence we conclude that the R/P ratio does not provide a reasonable predictor for the end of the age of oil.

M. King Hubbert was interested in the answers to the two fundamental questions in the context of oil, coal, and natural gas, and their application to various regions, such as individual oil fields, Texas, the US, and finally, the world. In this paper, we focus on world oil. To address his questions, Hubbert proposed a more sophisticated approach than the simplistic R/P ratio. He realized that production rate curves, as functions of time, for oil fields, whole geographic regions, and other fossil fuels share the feature in which the “curve starts slowly and then rises more steeply until finally an inflection point is reached after which it becomes concave downward” (Hubbert 1956, p. 6). He deduced that “the significance of this is that during the initial stages [of resource exploitation] all of these rates of production tend to increase exponentially with time” (Hubbert 1956, p. 7). However, “no finite resource can sustain for longer than a brief period such a rate of growth of production; therefore, although production rates tend initially to increase exponentially, physical limits prevent their continuing to do so” (Hubbert 1956, p. 8). He proposed two requirements for his models (Hubbert 1956, p. 9):

Consequently, to satisfy requirements 1 and 2, Hubbert settled on a resource production rate model that is a family of curves bearing the characteristic shape of a bell (Hubbert 1956, fig. 11, p. 10). These logistic curves depend on three parameters, Q_∞, r, and b. They are described algebraically by the formula

From eq. (1) Q(t) can be recovered as:where Q_o = Q(0) and $Q_{\infty }={\lim }_{t\to \infty }Q(t)$ . Equation (2) follows from eq. (1) because the production/extraction rate, P(t), is the time derivative of the cumulative extracted amount, Q(t).

Picking the parameters to match physical data and (or) estimates is an exercise in curve fitting. The first parameter, Q_∞, represents the UR. Hubbert used his expertise as an oil geologist to carefully justify a global estimate for Q_∞, the “initial amount present” for the whole world, namely Q_∞ =1250 Gbbl (Hubbert 1956, p. 22). This value becomes the area under the world oil production curve, in agreement with eq. (2). Next, the value b can be computed from the initial amount via b = ln(Q_∞/Q_o−1) (c.f. eq. (4), evaluated at t = 0). Finally, r is calculated from Hubbert’s estimate of the peak (maximum) production rate, “two and one-half times the present [1956] rate” (Hubbert 1956, p. 22). Thus, P_max = 12.5 Gbbl/year, and r comes from the relation P_max = rQ_∞/4. The result of this model by Hubbert (1956, fig. 20) “… places the [predicted] date of the peak at about the year 2000” (Hubbert 1956, p. 22).

In the mathematical modelling of physical phenomena, simple curve fitting has inherent weaknesses. A wide variety of curve families with simple algebraic structure can be chosen to fit a given set of data over time, but the usefulness of such models as predictors of the future is suspect, unless the choice is based on more fundamental dynamics or first principles. In this case, the family of curves Hubbert employed arises from the solution of a differential equation called the logistic equation. In this context, the logistic equation relates the production rate P(t) and the cumulative amount Q(t) dynamically in time, allowing for a more fundamental justification of the model. In the authors’ opinion, it is this more fundamental justification that gives the Hubbert peak oil model its credibility and enduring impact. We therefore attempt the following detailed presentation.

Noting that P is the time rate of change of Q (i.e., P(t) = Q′(t)), the logistic differential equation for Q(t) is:

This model is to be satisfied for all t > 0 and for the given (small, positive) initial value, Q_o = Q(0). The solution of the differential equation is the unknown Q(t), which satisfies the initial value at t = 0 for given (positive, constant) parameters, r and Q_∞. As the differential equation depends only on Q and the model parameters, and not explicitly on t, it is termed time-autonomous. The model from eq. (3) requires that P = Q′ be jointly proportional to a “growth term”, Q, and to a “limiting term”, (1 − Q(t)/Q_∞). We examine each of these model features in order.

Suppose the limiting term had little effect and could be ignored. This situation would occur for the time period near the outset at t = 0, when Q ≈ 0, which gives the limiting term a value approximately equal to one. In this initial period the model essentially behaves like the easier relation, Q′(t)/Q(t) = r (take the limiting term to be 1 in eq. (3) and divide through by the remaining operative term, Q). This latter equation says that the relative growth rate, Q′/Q, is a constant, r, allowing us a first principles understanding of the r parameter as the relative growth rate. The solution would now be exponential, Q(t) = Q_oe^rt and P(t) = Q′(t) = rQ_oe^rt. As desired, for a brief time near the outset, the production rate starts close to zero and increases exponentially at relative growth rate r. However such a model cannot be valid for an extended time period, as P(t) has no maxima and never again declines (much less to zero). Rather, P(t) → ∞ as t → ∞, violating the physical limits for a finite resource. To get the production rate to reach a maximum, and decline again toward zero (producing the characteristic bell shape), the model must limit the growth term. Thus, the “limiting term” begins to take effect as Q increases, thereby reducing the effective relative growth rate.

So when Q(t) is not small, the limiting term, (1 − Q(t)/Q_∞) is not one, has effect, and cannot be ignored. Again dividing by Q, eq. (3) can be written Q′(t)/Q(t) = r(1 − Q(t)/Q_∞). Now it is apparent that the relative growth rate is not the constant r, but rather a linearly decreasing function of Q, whose value near the outset is close to r, and decreases to zero as Q(t) grows to Q_∞. Thus, as Q(t) (and therefore its derivative, P(t)) initially grows exponentially, the limiting term (1 − Q(t)/Q_∞) decreases. This decrease in the limiting term causes the production rate Q′ = P to peak and then to decline. Because the production rate declines to zero (as desired) as the cumulative amount grows toward Q_∞, it does so increasingly slowly and never exceeds Q_∞. Thus, Q_∞ represents the UR, the amount Hubbert carefully estimated to answer question 1, before using it in his logistic model to answer question 2. With this interpretation of Q_∞ as UR, the limiting term represents the proportion of world oil remaining in the ground at t. The authors view the magnitude of Q(t) as a proxy measure of the current size of the aggregate global oil economy, including the amounts invested in equipment and oil infrastructure, the number of uses of oil, and the sheer number of end users. In this view, the rate of change, P(t) = Q′(t), would indeed be proportional to Q(t), as per the logistic model in eq. (3).

Thus, the logistic differential equation, eq. (3), has the following fundamental dynamic interpretation in the context of resource extraction: it models a production rate that changes (grows or declines) simultaneously proportional to

The parameter r is the relative growth rate of cumulative extracted resource, in the absence of any resource limit. The parameter, Q_∞, represents the physical resource limit, the UR.

The solution of the logistic equation, eq. (3), is achieved by the standard technique of separating variables and integrating, which yields:from which P(t) = Q′(t) is given by eq. (1) (see Figs. 3a and 3b). The characteristic shape of Q(t) is sigmoidal, rising asymptotically from 0 as t increases from −∞, passing through an inflection point, and asymptotically growing to Q_∞ as t → ∞. It is symmetric about its inflection point, t^⋆ = b/r, where P′(t^⋆) = Q″(t^⋆) = 0. Therefore, at the inflection Q(t^⋆) = Q_∞/2. Consequently P(t) is symmetric about its peak, P(t^⋆) = rQ_∞/4 = P_max. Thus, according to Hubbert’s model, exactly half the resource remains in the ground at the time when the production level peaks. Subsequent declines in the production level precisely mirror the swiftness of the earlier production level rise.

So what do the data show? We make three comparisons: Hubbert’s 1956 model versus subsequent data for the US, Hubbert’s 1956 model versus subsequent world data, and finally world data to 2016 used to “best fit” an updated Hubbert logistic curve.

Given that we have demonstrated the inaccuracy of the Hubbert logistic model, eq. (3), to adequately capture the dynamics driving the data sets for reserve, production, and cumulative extraction of world oil, answers to our questions 1 and 2 seem elusive. The particular problem seems to be the ongoing growth of reserve, making earlier estimates of UR either wrong or, at best, poorly justified. As Morehouse (1997, p. vii) stated, “reserves growth” is an “intricate puzzle”, and he names the phenomenon “ultimate recovery appreciation”, “the generally observed increase of EUR over time”. Reserve growth is driven by two related factors: the discovery of new oil fields and technological innovation (both new and improved technologies).

As an example, we point to the case of the Hibernia oil project off the coast of Newfoundland, Canada. When the original development project was submitted for government approval before production began in 1997, the application included a reserve estimate of “recoverable oil of 563 million barrels” and a lifespan of about 20 years. However, CBC News reported in December 2016 (Roberts 2016), that the billionth barrel was pumped on the 21st of that month, and that the current expectations are for production to continue until about 2040, recovering 1.6 billion barrels cumulative. The reserve growth is due in part to new abilities to extend drilling horizontally.

The report is striking; after exceeding the oil-field production target by 77% and approaching the end of the projected lifespan of the field, the reserve remaining actually exceeds the original projections of the resource. One wonders whether this same report will be repeated in a further 20 years, with an announcement that the remaining oil still exceeds a half-billion barrels. Yet such a possibility cannot be discounted, even while reason suggests the process cannot continue indefinitely.

A successful oil production model must account for the phenomenon of reserve growth. A wide variety of Hubbert model variants have been proposed in the literature. These include multi-Hubbert, modified Hubbert, and higher-order Hubbert models, such as those proposed by Campbell and Laherrère (1998), Berg and Korte (2008), Sorrell and Speirs (2010), Maggio and Cacciola (2012). Various other curve-fitting processes have also been employed (see Lynch (2003) and Brandt (2010)). Here we propose a simple, new, “Hubbert variant” model that allows for dynamic reserve growth with a reasonable justification. The model is a system of two time-autonomous differential equations, one for cumulative oil, Q(t), and a second for reserve, R(t), including three positive-value parameters:

The first, eq. (8), reproduces the Hubbert logistic model with one modification. Here, the growth-limiting term, the remaining unextracted resource, Q_∞ − Q(t), has been replaced by the current reserve, R. Because reserve is only the best 1P estimate of the remaining resource, as R(t) ≤ Q_∞− Q(t), the impact of this new growth-limiting term is felt earlier in the time history.⁴ The second, eq. (9), says that the rate of change of the reserve has two independent (additive) factors. The first of these is obvious: reserve declines at an annual rate equal to the amount extracted each year by the current production rate, P. The second term also allows reserve to grow logistically, with relative growth rate α and growth-limiting constant β. In this view, reserve, like a water-loaded sponge, can always have some additional percentage α of the current loading R added to R, because we can squeeze the sponge. Here the ability to “squeeze” the sponge mimics the effect of technology-driven resource expansion. However, this process cannot continue indefinitely and growth is therefore limited by the factor (1 − βR).

Note that in the limiting case where reserve is not allowed to grow (i.e., take α = 0), the model reverts to the Hubbert model. Setting α = 0 in eq. (9) makes R′ = −P = −Q′ (with the latter equality by the definition of production rate). Integrating now implies that R(t) = −Q(t) + Q_∞, where the constant of integration is chosen to match the total reserve when cumulative is zero, i.e., when reserve was UR. Using this expression to replace R in eq. (8) and recalling that ρ = r/Q_∞, we recover the Hubbert model. Therefore, when we express the solution of the Hubbert model in (Q, R) coordinates (with Q on the horizontal axis), the curve is the straight line, R + Q = Q_∞, with time increasing along the line to the right—the direction of increasing cumulative. For different assumed values of Q_∞ the Hubbert lines are all parallel, with intercept values Q_∞ on both axis. In Fig. 7 two possible Hubbert lines are plotted, the first using the current, 2016, value of R + Q = 3130 [Gbbl] as the Q_∞.

The Hubbert variant model can be solved analytically: dividing eq. (9) by eq. (8) and using the chain rule yields a linear differential equation for R as a function of the independent variable, Q:which is solved in the standard way. The solution is:where c is the arbitrary constant of integration and a represents the collection of parameters a = αβ/ρ. The model, eq. (11), has three, positive parameters to be determined, a, c, and 1/β. We used the “NonlinearFit” routine of the Maple^® 15 mathematics package to find optimal values (in the least squares sense) for these parameters to best fit this model equation to our 37-point data sets for R and Q. The result is plotted in Fig. 7, in the plane of cumulative versus reserve.

The “best fit” Hubbert variant model curve in Fig. 7 does fit the BP data reasonably well. From the graph, the horizontal intercept where R = 0 is interpreted as the UR. The model predicts a UR of 5484.2 Gbbl, answering question 1. Using this value as Q_∞ produces the second Hubbert line on the plot. Given that all Hubbert model curves (for any chosen value of Q_∞) are lines parallel to the two in Fig. 7, it should be clear that none bear any resemblance to the actual data. This observation again confirms that the Hubbert model for world oil is highly inaccurate and is therefore refuted. To some, the above predicted UR (5.5 terabarrels (Tbbl)) may seem unjustifiably high.

As for the variant model, we add the following observations. Time is an (unknown) parameter along the curve. Differentiating eq. (8) with respect to time (having first replaced R with R(Q(t)) from eq. (11)) and numerically finding the zero, allows us to say, definitively, that the best fitting model of this type predicts peak world oil production to occur when the cumulative amount of oil is Q = 2998 Gbbl, for which the corresponding reserve has declined to R = 1199.7 Gbbl. To determine the time at which this event occurs requires us to estimate the parameter ρ in the differential equation, eq. (8). Unfortunately, there is no analytic solution available to allow the parameter to be estimated from the data to minimize the error. However, we can make the following rough estimate of ρ; using the most recent (2016) data value for Q, P, and R, we estimate ρ = P(2016)/[Q(2016)R(2016)]. Now eq. (8) can be solved numerically, using the initial condition of Q(2016), thereby creating a time series for Q. Then peak world oil occurs when Q = 3000, for which t = 73.8 [years since 1979], i.e., 2052. The production level at that time is approximately 50 Gbbl/year, which is almost 50% higher than current production. This estimate for the timing of peak oil is an upper bound on values produced by various estimates of the parameter, ρ. Also, world cumulative exceeds Q = 4 Tbbl about 2075 and production effectively ceases by the time Q = 5.4 Tbbl at 2160, giving a rough answer to question 2. We reiterate that these are just rough estimates: the computed time value and production value for the “peak” are sensitive to the estimate of ρ in eq. (8).

Hubbert need not feel alone in making unrealized predictions concerning world oil production. In a comprehensive review article Brandt (2010, p. 3958) classified 45 “mathematical models of oil depletion of the last century”, all of which attempt to forecast future oil production paths (most with global scope). Concerning the predictive capacity of all of these models, he concluded:

Existing models have fared poorly in predicting global oil production. Even for models that are commonly thought to be successful, after-the-fact interpretation of the success or failure of a predictive effort is not easy (recall the discussion above of Hubbert’s successful prediction). … This author’s judgement with respect to the predictive value of models is as follows: … There is no scientific justification for making specific predictions (e.g., the year of peak production) with any of the surveyed mathematical models: the uncertainties involved make such predictions of little use. (Brandt 2010, p. 3972)

Despite this advice, “Peak Oil” papers and predictions continue to be published and then proven false. For example, Mohr and Evans (2008) made a 2007 prediction for (UR) 2.2 trillion barrels with a world oil production peak in 2013, with the caveat that the peak will be somewhat later if the UR reaches 3 trillion barrels. However, the data (Figs. 1a, 1b, and 5) show that the 2016 value of R + Q = 3130 Gbbl is already >3 trillion barrels, with no production peak in sight, therefore completely refuting the predictions. As a final example, Gallagher (2011) used an “idealized Hubbert curve” to make a 2010 prediction that peak oil had just occurred in 2009 at 30.4 Gbbl/year with UR of 2.24 trillion barrels. Again, these values have already been surpassed, indicating that the prediction has failed.

Brandt ended the earlier quote of his conclusions concerning oil production predictions with the following exhortation (Brandt 2010, p. 3972): “Efforts should move away from making these kinds of predictions, and toward understanding the impacts of the inevitable transition to oil substitutes.”

And so we are back to the beginning: we have not provided credible answers to either of the two fundamental questions. We acknowledge that the answers are elusive because prediction (beyond predicting a regression line trend for the very near term) is extremely difficult. Clearly the underlying dynamics of reserve growth are poorly understood, and though many have made oil production peak and depletion time predictions, none can be trusted. Finally Brandt (2010) points us towards contemplation of the broader aspects of the investigation, along the lines of the question from Science (Kerr and Service 2005): “What can replace cheap oil—and when?”

The authors acknowledge helpful discussions with, and computational aid from, Bishop’s University summer undergraduate research assistants, Jacob Courtemanche (2015), and Courtney Drew (2017).