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".

SOURCE CODE

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
   begin
      if I in Arr'Range then
         return Arr (I);
      else
         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;
   begin
      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);
   begin
      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, 
                                                              Avg_Mature)));
   -- Safe function for reading forcing function
   function Discovery_Window (Year : Float; Start : Float) return Float is
      Y : constant Float := Year - Start;
   begin
      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;
   begin
      for J in  S'First + 1 .. S'Last loop
         if S (J).Year > Year then
            K := J - 1;
            exit;
         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;
begin
   loop
      -- 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;
         Text_IO.Put_Line
           (Integer'Image (Integer (Time)) & "," &
            P'Img & "," &
            Data'Img & "," &
            Float'Image (V));
      end if;
      Time := Time + DT;
      exit when Time > Finish;
   end loop;
end ASPO;