Type Curve Economics

In Part 1 we downloaded production data and manipulated it. In Part 2, we learned how to create type-curves. Now, let’s create a quick single-well economic model, and I suppose play around with creating functions for the first time.

Follow My Blog

Functions

Not as bad as you think. Generally the rule of thumb is that if you write the same bit of code more than 3 times, then you should be using a function. Of course, if you’ve been following along with the tutorials so far then you’ve already been using functions. The pre-built R packages are really just a lot of bundled functions.

But now, we’ll just build one ourself.

Revisiting aRpsDCA

You will recall our single-segment Arps Modified Hyperbolic Decline curve from Part 2.

library(aRpsDCA)

qiGas <- 20000 #mcfd
bGas <- 1 #B-factor
DiGas <- 0.9 #Effective Annual Decline (Initial)
DfGas <- 0.08 #Effective Annual Decline (Terminal)
curtail <- 5 #Months of curtailment
wellLife <- 30 #Years

gasFcst =  curtailed.q(arps.decline(
  qiGas*365.25, as.nominal(DiGas), bGas,
  as.nominal(DfGas)),
  curtail/12.0,seq(1/12, wellLife, by= (1/12)))/12

There are actually three functions here, curtailed.q, arps.decline, and as.nominal. All are defined within the package documentation for aRpsDCA. These allow me to forecast well production. What I want to do is create a single function that will spit out a dataframe of oil, gas, and water production in monthly increments.

Our Function

There are two things I need to do:

Determine the variables within function()
Create the calculation after function()

So I will just use the above code as my guide.

prodFunc <- function(wellLife, 
                     qiGas, DiGas, bGas, DfGas, curtailGas,
                     qiOil, DiOil, bOil, DfOil, curtailOil,
                     qiWater, DiWater, bWater, DfWater, curtailWater) {
  
  df1 <-  data.frame(
    #Months of Prouction
    monthsOn = seq(1, wellLife*12, 1),
    #Define Gas
    gas = curtailed.q(arps.decline(
      qiGas*365.25, as.nominal(DiGas), bGas,
      as.nominal(DfGas)),
      curtailGas/12.0,seq(1/12, wellLife, by= (1/12)))/12,
    #Define Oil
    oil = curtailed.q(arps.decline(
      qiOil*365.25, as.nominal(DiOil), bOil,
      as.nominal(DfOil)),
      curtailOil/12.0,seq(1/12, wellLife, by= (1/12)))/12,
    #Define Water
    water = curtailed.q(arps.decline(
      qiWater*365.25, as.nominal(DiWater), bWater,
      as.nominal(DfWater)),
      curtailWater/12.0,seq(1/12, wellLife, by= (1/12)))/12
  )
  
  return(df1)
}

Calculation

Now that our function is build, let’s test it out. To avoid confusion, I’m going to throw a 1 behind each variable I define.

qiGas1 <- 2000 #mcfd
bGas1 <- 1.2 #B-factor
DiGas1 <- 0.8 #Effective Annual Decline (Initial)
DfGas1 <- 0.08 #Effective Annual Decline (Terminal)
curtailGas1 <- 3 #Months of curtailment
qiOil1 <- 1000 #bopd
bOil1 <- 1 #B-factor
DiOil1 <- 0.9 #Effective Annual Decline (Initial)
DfOil1 <- 0.08 #Effective Annual Decline (Terminal)
curtailOil1 <- 1 #Months of curtailment
qiWater1 <- 1000 #bwpd
bWater1 <- 1 #B-factor
DiWater1 <- 0.9 #Effective Annual Decline (Initial)
DfWater1 <- 0.08 #Effective Annual Decline (Terminal)
curtailWater1 <- 0 #Months of curtailment
wellLife1 <- 30 #Years

#Create a dataframe with our function
prod <- prodFunc(wellLife = wellLife1, 
                     qiGas = qiGas1, DiGas = DiGas1, bGas = bGas1, 
                     DfGas = DfGas1, curtailGas = curtailGas1,
                     qiOil = qiOil1, DiOil = DiOil1, bOil = bOil1, 
                     DfOil = DfOil1, curtailOil = curtailOil1,
                     qiWater = qiWater1, DiWater = DiWater1, bWater = bWater1, 
                     DfWater = DfWater1, curtailWater = curtailWater1)
head(prod)

I get this.

monthsOn	gas	oil	water
1	60875	30438	25537
2	60875	25537	21996
3	60875	21996	19317
4	53756	19317	17220
5	48244	17220	15534
6	43839	15534	14148

Table 1: Production Data Frame (1st 6 columns)

If we sum each column, that should give us our EUR (Estimated Ultimate Recovery).

sum(prod$oil)/1000 #Convert Oil EUR from BBL to MBO

666.5 MBO for this example.

Putting it All Together

Prices

Now that I can create a production function, I can start tacking on additional functionality. For this example, I am going to use a forward strip price from the Wall Street Journal for my forecast, as I covered in Scraping Strip with R (Oil & Gas Coding Series). But this time I’ll make it a function!

library(rvest)
library(httr)
library(dplyr)
library(zoo)

priceFunc <- function() {
  crude = 'https://quotes.wsj.com/futures/CRUDE%20OIL%20-%20ELECTRONIC/contracts'
  webpage <- read_html(crude)
  tbls <- rvest::html_nodes(webpage, 'table')
  
  tbls_ls <- webpage %>%
    rvest::html_nodes('table') %>%
    .[1] %>%
    rvest::html_table(fill = TRUE)
  
  wti <- tbls_ls[[1]] %>% select(MONTH, SETTLEMENT)
  wti <- wti %>% filter(MONTH != 'Front Month')
  wti$MONTH <- gsub('Crude Oil', '', wti$MONTH)
  
  wti$MONTH <- trimws(wti$MONTH, which = 'both')
  
  wti$Year <- substr(wti$MONTH, 5, 8)
  wti$Month <- substr(wti$MONTH, 1, 3)
  wti$Date <- as.Date(paste0(wti$Month, '/01/', wti$Year), format = '%b/%d/%Y')
  wti <- wti %>% select(DATE = Date, WTI = SETTLEMENT)
  
  gas = 'https://quotes.wsj.com/futures/NATURAL%20GAS/contracts'
  webpage <- read_html(gas)
  tbls <- rvest::html_nodes(webpage, 'table')
  
  tbls_ls <- webpage %>%
    rvest::html_nodes('table') %>%
    .[1] %>%
    rvest::html_table(fill = TRUE)
  
  hh <- tbls_ls[[1]]%>% select(MONTH, SETTLEMENT) %>% filter(MONTH != 'Front Month')
  hh$MONTH <- gsub('Natural Gas', '', hh$MONTH)
  hh$MONTH <- trimws(hh$MONTH, which = 'both')
  hh$Year <- substr(hh$MONTH, 5, 8)
  hh$Month <- substr(hh$MONTH, 1, 3)
  hh$Date <- as.Date(paste0(hh$Month, '/01/', hh$Year), format = '%b/%d/%Y')
  
  hh <- hh %>% select(DATE = Date, HH = SETTLEMENT)
  
  wti <- wti %>% left_join(hh)
  wti <- wti %>% filter(!is.na(WTI)) %>% filter(!is.na(HH))
  
  return(wti)
}

Then, to call my price file, I just run:

prices <- priceFunc()

So, now if I want to rerun my price file, I can just run priceFunc() and then join it to my production dataframe.

You could potentially throw this function into our original prodFunc(), but that is not preferable. You will eventually be running a lot of different variations of type-curves in the future, so it is preferable to only run the priceFunc() at the beginning of your session. Otherwise you will be trying to scrape prices each time you change a type-curve, and that is SLOW.

Below would be the preferred workflow.

prod <- prodFunc()
prices <- priceFunc()

#Add similar monthsOn column to prices
prices <- prices %>% mutate(monthsOn = seq(1, n(), 1))

#Join prices to production file on monthsOn column
prod <- prod %>% left_join(prices)

#Drag missing values downwards
prod$WTI <- zoo::na.locf(prod$WTI)
prod$HH <- zoo::na.locf(prod$HH)

#View first rows
head(prod)

monthsOn	gas	oil	water	DATE	WTI	HH
1	60875	30438	25537	2020-10-01	38.05	2.406
2	60875	25537	21996	2020-11-01	38.41	2.873
3	60875	21996	19317	2020-12-01	38.87	3.241
4	53756	19317	17220	2021-01-01	39.34	3.361
5	48244	17220	15534	2021-02-01	39.78	3.323
6	43839	15534	14148	2021-03-01	40.18	3.201

Table 2: Production Data Frame with Price (1st 6 columns)

Revenue

Now I will just add in a quick revenue model based on absolutely nothing but probably close enough. What will I need for this?

Working Interest/Net Revenue Interest
Shrink (or as I define it Gas Sold as % of Gas Produced)
NGL Yield (BBL/MMCF of wellhead gas)
Diff’s (I’ll just use simple percents here but can be per-unit)

wi <- 1
nri <- 0.75 #To the 100% of WI
shrink <- 0.75
nglYield <- 100
oilDiff <- 0.95
gasDiff <- 0.8
nglDiff <- 0.3 #%WTI

prod$ngl <- prod$gas*nglYield/1000
prod$oilRevenue <- prod$oil*wi*nri*prod$WTI*oilDiff
prod$gasRevenue <- prod$gas*wi*nri*prod$HH*gasDiff*shrink
prod$nglRevenue <- prod$ngl*wi*nri*prod$WTI*nglDiff
prod$revenue <- prod$oilRevenue+prod$gasRevenue+prod$nglRevenue

Expenses

Expenses are usually based on the Working Interest ownership, so that’s what I will use here. There is a lot you can do, but for my analysis I am completely making it up (or maybe it’s influenced by actuals). For this model I will assume:

Fixed Cost at $10,000 $/month until total fluid rate goes below 100 barrels per day, or 3,000 per month. $5,000 thereafter.
Variable Expenses of $0.5/bbl oil, $0.2/mcf wellhead gas, $1/bbl water.

You could potentially run assumed Workover Expenses as well, but I usually would lump them into fixed cost assumptions for ease.

fixed1 <- 10000
switch <- 100*30.45
fixed2 <- 5000
varOil <- 0.5
varGas <- 0.2
varWtr <- 1

prod$expense <- fixed1 + prod$oil*wi*varOil+prod$gas*wi*varGas + prod$water*wi*varWtr
prod$expense[prod$oil + prod$water < switch] <- prod$expense[prod$oil + prod$water < switch] - fixed2

Taxes

I am going to use Texas severance taxes, assuming an oil well, and then assume 2.5% of revenue for Ad Valorem.

stxOil <- 0.046 #Texas Oil Severance
stxGas <- 0.075 #Texas Gas/NGL Severance
atx <- 0.025 #Ad Valorem

prod$tax <- prod$oilRevenue*stxOil + (prod$gasRevenue+prod$nglRevenue)*stxGas + prod$revenue*atx

Cash Flow

Simply enough, I will just calculate total cash flow. However, I am going to assume that the well is shut in on the last positive cash flow month.

prod$cf <- prod$revenue - prod$expense - prod$tax
last1 <- which(prod$cf > 0)
last1 <- last1[length(last1)]

prod <- prod[1:last1,]

Capex/PNA

A somewhat lazy assumption I will make is that P&A costs are incurred at the last month of the data frame. It is not that big of an impact from a PV-perspective that far out into the future, though you COULD defer it if it so pleases you. I also need to define capex and then how long we go from SPUD to first production date.

drillCapex <- 4000000 #Drill Cost
complCapex <- 6000000 #Capex + Facilities
pna <- 50000 #P&A
spudToProd <- 3 #Production occurs on month 4.

prod <- prod %>% subset(select = -c(DATE)) #Remove useless Date column
prod$capex <- 0  #Create capex column
prod$capex[nrow(prod$capex)] <- pna #Add P&A to last month

capex1 <- prod[1:spudToProd,]
capex1[1:nrow(capex1),1:length(capex1)] <- 0
capex1$capex[1] <- drillCapex
capex1$capex[spudToProd] <- complCapex

prod <- rbind(capex1, prod)

prod$monthsOn <- seq(0, nrow(prod)-1, 1) #New monthly column

Cash Flow/NPV-10

Now that we have the cash flows built out, we can have a quick peek at single-well economics. I am going to show discount rates assuming annual compounding (ie PV = CF / (1 + r/n) ^t*n). Obviously, if you change the discount rate, that changes the PV.

prod$fcf <- prod$cf - prod$capex
r <- .1
n <- 1

prod$pv <- prod$fcf/((1+r/n)^(prod$monthsOn/12*n))

library(scales)
scales::dollar(sum(prod$pv))

Based on my price assumption (strip as of 9/9/2020), then the NPV-10 is $3.7 Million US.

Figure 1: Cumulative Net Present Value (10%) – Example Well

IRR

Instead of an IRR package in R, I’m going to build out an IRR function that loops quickly until we find a value, or returns 0 if it’s negative.

IRRcalc <- function(cf, months){ #Function inputs are our Free Cash Flow profile and our months
  if(sum(cf) > 0){
    IRR1 <- 3 #Start at 300%, because, come on, like it'll ever even be near this.  Just starting high.
    loop <- 1 #Initiate loop at 1
    while(sum(cf/((1+IRR1)^(months/12))) #While with our NPV Formula
          < 0){
      IRR1 <- IRR1 - 0.01 #If not true, knock off 1%
      loop = loop + 1 #Next Loop
    }
    
  }else {
    IRR1 <- 0 #If total cf is less than 0, return 0
  }
  return(IRR1)
}

print(scales::percent(IRRcalc(prod$fcf, prod$monthsOn)))

In this scenario, I get a 24% IRR.

Summary

That’s probably enough brain-hurt for this session. If you want to test yourself, try and combine everything into a single function and see how it goes. I will be doing a pseudo-version when I publish Part 2 of the Shiny Version later this week. Thanks for tuning in!