[[ Check out my Wordpress blog Context/Earth for environmental and energy topics tied together in a semantic web framework ]]

Thursday, October 06, 2005

Shock Model with ASPO Discovery Data

Update: This page has had enough references that I should point to an overview.

PeakOil analyst khebab pointed to the ASPO global discovery curve which intrigued me enough to use it as an input forcing function to my oil shock model.

(note that the curve shown above includes a 3-year moving average)

I had to really understand all the latencies involved from the actual point of discovery, i.e. the initial "oil strike", to the point at which each discovery starts producing oil at a mature clip. So based on the elements of the Micro peak oil model, I thought it best to create linear composable latencies that served to shift the point of discoveries in time by the 36 years or so that others empirically observe.
The discovery curve mirrors approximately the production curve with a lag that varies from country to country. The US-48, for example, had a lag of 41 years whilst the UK North Sea production, with its urgency and technological basis had a lag of 25 years. The World's lag is estimated to be 36 years.
The composable stochastic latencies include:
  1. Mean time from discovery to decision to extract -- The "fallow" period
  2. Mean time from decision to extract to completing rig construction -- The "build" period
  3. Mean time from construction complete to maturity of production -- The "maturation" period
If these times remain independent, we can use the convolution technique to generate a mature discovery window. I know that the fallow period can range well over 10 years; for example the industry knew about Alaskan oil well before they made the decision to start extracting. As for the build period, I have a few references that suggest that it takes a minimum of 3 years to construct an oil rig on land and 5 years for an off-shore platform. And the maturity period includes all sorts of extraneous considerations, such as support features (building pipelines, etc) and the possibility of dry wells caused by improper placement. I ended up choosing a Markovian 8-year average latency for each of these phases. This caused the mature discovery curve to match fairly well to the triangular discovery curve that I used in the past for the oil shock model (of course this curve does not show the sharp peaks of the triangular curve).

Finally, I added the previously described extraction phase to the model. The extraction rate basically relates the mean time to deplete a reservoir to 1/e of its original value (or to 36.8% of its initial volume). As suggested before, I use as an implicit assumption that any rate of extraction or flow is proportional to the amount available and nothing more; past and future history do not apply.

And as the final necessary ingredient to match the spiky/notchy behavior of global production, I add the oil shocks which tend to suppress the production during critical geo-political periods. The following fit to the BP data assumes the initial ASPO data as the forcing function, the 3 mean latencies occurring before mature production, an initial extraction rate, and the 3 shocked (or perturbed) extraction rates to match the dips. As the final frill, I then added a reverse shock starting after 2001, on which Bubba previously posted some suspicions :
"The orange squares are the BP cumulative data, right? What is causing the upturn in the orange square trend around 2004?"

Because of my use of external data for discoveries, which actually constitute someone else's estimates and therefore may turn out low, take the curve at face value. I don't know the history behind the ASPO numbers, but just as a cautionary note, I can summarize with this statement : "gulp".


Prints to standard out "Year", "ModelResults", "BPData"

with Text_IO;
with Ada.Numerics.Elementary_Functions;

procedure ASPO is
type Flt is array (Natural range <>) of Float;

-- Safe array retrieval function
function Get (Arr : Flt; I : Integer) return Float is
if I in Arr'Range then
return Arr (I);
return 0.0;
end if;
end Get;

-- Discretized convolution function
function Convolve (A, B : in Flt) return Flt is
Total : constant Natural := A'Length + B'Length;
C : Flt (0 .. Total);
V : Float;
for J in 0 .. Total loop
V := 0.0;
for I in 0 .. J loop
V := V + Get (A, I) * Get (B, J - I);
end loop;
C (J) := V;
end loop;
return C;
end Convolve;

-- Discretized exponential function
function Exponential (L : Natural; Alpha : Float) return Flt is
use Ada.Numerics.Elementary_Functions;
R : Flt (0 .. L - 1);
for I in 0 .. L - 1 loop
R (I) := Alpha * Exp (-Alpha * Float (I));
end loop;
return R;
end Exponential;

-- Data transcribed from ASPO discovery curve
-- http://wolf.readinglitho.co.uk/mainpages/discoveries.html
Strikes : constant Flt := -- GBBls/year starting year 1932
( 9.0, 4.0, 4.5, 3.0, 5.0, 5.5, 42.0, 42.5, --1930's
52.0, 18.0, 16.0, 5.0, 3.0, 7.4, 7.6, 7.0, 49.0, 53.0, --1940's
54.0, 19.5, 16.0, 26.0, 20.0, 28.0, 24.0, 35.0, 40.0, 40.5, --1950's
42.0, 44.0, 50.0, 41.0, 50.0, 48.5, 47.5, 32.0, 30.5, 30.4, --1960's
29.5, 40.5, 37.0, 38.5, 24.0, 26.5, 31.0, 37.0, 38.0, 36.0, --1970's
26.0, 24.0, 20.0, 20.0, 22.0, 21.0, 20.5, 18.0, 16.2, 17.5, --1980's
16.5, 20.0, 17.0, 16.0, 9.0, 8.5, 9.0, 9.0, 9.5, 14.0, --1990's
18.0, 17.5, 13.0, 9.0, 9.0); --2004

Avg_Fallow : constant Flt := Exponential (200, 0.125); -- 8 years
Avg_Build : constant Flt := Exponential (200, 0.125);
Avg_Mature : constant Flt := Exponential (200, 0.125);
Discovery : constant Flt := Convolve (Strikes,
Convolve (Avg_Fallow,
Convolve (Avg_Build,
-- Safe function for reading forcing function
function Discovery_Window (Year : Float; Start : Float) return Float is
Y : constant Float := Year - Start;
if Y < 0.0 then
return 0.0;
end if;
return Discovery (Integer (Y));
end Discovery_Window;

-- British Petroleum data
BP_data : array (1965 .. 2500) of Integer := -- in million Bls/day
(1965 => 31803,
1966 => 34568,
1967 => 37118,
1968 => 40436,
1969 => 43633,
1970 => 48061,
1971 => 50844,
1972 => 53666,
1973 => 58463,
1974 => 58617,
1975 => 55824,
1976 => 60412,
1977 => 62713,
1978 => 63331,
1979 => 66049,
1980 => 62946,
1981 => 59533,
1982 => 57296,
1983 => 56598,
1984 => 57683,
1985 => 57468,
1986 => 60461,
1987 => 60785,
1988 => 63160,
1989 => 64051,
1990 => 65470,
1991 => 65288,
1992 => 65788,
1993 => 66046,
1994 => 67116,
1995 => 68103,
1996 => 69895,
1997 => 72158,
1998 => 73586,
1999 => 72333,
2000 => 74950,
2001 => 74828,
2002 => 74443,
2003 => 77054,
2004 => 80260,
others => 0);

type Shock is record
Year : Float; -- Year of shock
Rate : Float; -- Rate of extraction
end record;
type Shocks is array (Natural range <>) of Shock;

S : constant Shocks :=
((1932.0, 0.070), -- Start of data
(1974.0, 0.070), -- Start of oil embargo
(1974.5, 0.065), -- End of oil embargo
(1980.0, 0.065), -- Start of Iran hostage crisis
(1983.5, 0.044), -- End of recession
(1990.0, 0.044), -- Start of Gulf War
(1992.0, 0.042), -- End of recesssion
(2001.0, 0.042), -- Last good year!
(2003.0, 0.046), -- Running Out??? Why is extraction going up?
(2100.0, 0.046));

-- Make the shock a continuous function
-- by interpolating between points in the list
function Interpolate_Shocks (Year : Float) return Float is
K : Integer := S'Last - 1;
for J in S'First + 1 .. S'Last loop
if S (J).Year > Year then
K := J - 1;
end if;
end loop;
return S (K).Rate +
(S (K + 1).Rate - S (K).Rate) * (Year - S (K).Year) /
(S (K + 1).Year - S (K).Year);
end Interpolate_Shocks;

C : Float := 0.0; -- Reserve
Rate : Float; -- Extraction Rate
P : Float; -- Production Rate
DT : constant Float := 0.01; -- 1/100th of year
Start : constant Float := S (S'First).Year;
Finish : constant Float := S (S'Last).Year;
Time : Float := Start;
V : Float := 0.0; -- Volume of oil extracted
Data : Float := 0.0;
-- Integration of discovery window against extraction rate
Rate := Interpolate_Shocks (Time);
C := C + Discovery_Window (Time, Start) * DT - C * Rate * DT;
P := Rate * C;
V := V + Rate * C * DT;
if Time > Float (Integer (Time)) - DT / 2.0 and
Time < Float (Integer (Time)) + DT / 2.0
then -- Print output only once per year
if Time >= Float (BP_data'First) then
-- BP data only available beyond a certain date
Data := Float (BP_data (Integer (Time))) * 365.0 / 1.0e6;
end if;
(Integer'Image (Integer (Time)) & "," &
P'Img & "," &
Data'Img & "," &
Float'Image (V));
end if;
Time := Time + DT;
exit when Time > Finish;
end loop;
end ASPO;


Professor Blogger Big Gav said...

Errr - "gulp" does kind of sum up that graph.

Could you please fix your code so it shows a "bumpy plateau" for 10 years while we work out WTF to do ?

3:02 AM  
Professor Blogger Mr. Sprang said...

Good match!

A few questions:

1. Is that a "backdated" discovery datebase from ASPO? Cause that's a potential source of error - it would systematically undercount more recent discoveries.

2. The projection for future discoveries looks a little low. The extrapolation line seems to come off the troughs in discovery, instead of a mean value between the troughs and peaks. It looks like the problem is that the extrapolation comes off a three-year moving average, the last value of which happens to be centered in a discovery trough (coincidentally missing the higher discovery rates of 2000-2001).

3. I'm skeptical because as soon as the historical data run out, the model takes the curve in a totally different direction. The fact that your model seems to be swamping the data - and is even in disagreement with the most pessismistic recent estimates from Campbell - is suspicious.

4. The older post mentions four pre-2000 shocks: 1973, 1979, 1984, and 1991. Which did you drop? And while 1973 and 1979 are obvious, the effects of the Gulf war were more transient, and the 1984 "end of recession" seems truly odd. The recession ended in 1982 - it lasted 3 quarters. What about the 1981 decision to end US domestic price controls (which had strangled domestic E&P), or the 1985 Saudi decision to abandon high prices and finally open their taps?

5. An unexplained fifth "reverse shock" has been added in 2001. Why? What does it represent?

10:42 AM  
Professor Anonymous Anonymous said...

Good Job! I like the way you deal with the delay between time of discovery and production start. This delay should also depend on the size of the field and the current balance between supply and demand.

I'm planning to implement your algorithm in the R language. I will probably have questions for you!

Be careful on the ASPO's definition of regular oil which excludes oil-from-coal, heavy, bitumen, shale, deepwater, polar. Their URR is therefore low because of that. The BP production data includes all liquids I believe.

11:28 AM  
Professor Blogger Mr. Sprang said...


If the ASPO data doesn't include non-conventionals and the BP data does, that explains a large part of the sharp crash after the most recent data!

1:29 PM  
Professor Anonymous Anonymous said...

that explains a large part of the sharp crash after the most recent data!

exactly, the low URR is pushing the curve down faster.

3:47 PM  
Professor Blogger SW said...

Nevertheless, I think we will find, in retrospect, that the discontinuity was real.

4:07 PM  
Professor Blogger @whut said...

It could be low for a number of other reasons as well. As I said, these are only estimates of oil discovered. Also take a look at this C. Campbell graph:
This includes the non-conventional oil sources, and seems to match well with the BP data.

By the way, the upward reverse shock indicates a sudden increase in extraction rate. So it's either the start of a bumpy ride or ?

Since the use of non-conventional oils only started increasing in use after 1970, it should not be too hard to refit and just assume a deeper suppressive shock in the last 30 years.

4:53 PM  
Professor Anonymous Anonymous said...

What happened in 2003-2004? In 2004 the Russian oil production increase accouinted for about 40% of of the total production growth. There was also using some spare capacity of Saudi-Arabia and all of OPEC.

Without these effects there probable would not be that kind of "reverse schock" we see in the curve. Now we know that Russia is not going to do the same again - production will be flat this year. I strongly suspect destructive overproduction there lately.

Basically this kind of technique cannot really predict production. It is always possible to make the curve fit the data afterwards but this doesn't tell the future. I remember one try to predict stock prices with a similar technique. First it seemed convincing: the historical data was fitted nicely with a complex equation but then the new real data jumped completely off the predicted course. With this I don't want to dismiss your work - only to be cautious.

12:42 PM  
Professor Blogger @whut said...

You don't seem to understand the way the stochastic differential equations work. Whatever is current is current, the past curve does not change. When new data comes along, you can make a snapshot of the previous data and then just keep on convolving with the new data. If you want to predict with data that does not yet exist, that is a problem with the unpredictability of the future.

I really don't see the difference between what this technique does versus what a reliability engineer does; he looks at historical data of failure rates and predicts the lifetime of a bunch of parts. If a bad batch comes along in the future, it doesn't invalidate the approach; he just has to incorporate, perhaps through Bayesian techniques, the new date into the estimates.

About the stock market ... tee-hee.

3:44 PM  
Professor Anonymous Anonymous said...

This is really great stuff. I've spent a while looking over it, and i'm really impressed. However, I believe there is one problem you'll want to address.

Campbell usually estimates a pre-1930 discovery of 122 Gb of regular oil. This appears to be unaccounted for in your method. It could explain why the curve seems so unnaturally constrained.

Also, isn't it the case that your extrapolation into the future using this method assumes no future discoveries? This may explain the apparent discontinuity at 2004.

Good Luck with your modelling.

8:59 AM  
Professor Anonymous Anonymous said...

This is a very interesting model.

Try to use Jean Laherrère's data for discovery instead, it is similar, but contains also the latest discovery cycle in deepwater.

11:35 AM  
Professor Blogger @whut said...

I have used Laherrere's data, but not for world, just for various regions.

I would use this one:

6:36 PM  

Post a Comment

<< Home

"Like strange bulldogs sniffing each other's butts, you could sense wariness from both sides"