Frac Ban in the Permian (Maps in R)

With the potential for regime change in the White House, shale operators and investors alike are rightly concerned about the potential of a frac ban on federal lands. When quantifying the impact, the first area we need to look at is the Permian’s Delaware Basin, and in particular New Mexico. Why?

Most productive part of the Wolfcamp play is in Lea County, New Mexico
Large percentage of Bureau of Land Management (BLM) acreage in New Mexico

In this insight, we are going to show how to do some high-level analysis using publicly available data and the R programming language to help quantify the impact of a frac ban. Of course, if you want to skip me being a nerd, just skip the coding sections and head on down to the Analysis portion.

Data Sources

Generally, I need two pieces of data to do this analysis.

Well Data
Federal Acreage

Luckily for me, New Mexico provides shapefiles for both.

Ideally, I would have individual lease ownership by operator, but that is difficult to get publicly. I would typically grab the acreage image from an investor presentation and throw it into QGIS, but that can be time-consuming. And I am trying to show how to do things quickly and efficiently.

So it won’t be perfect, and won’t accurately estimate non-drilled acreage, but it should be indicative of total acreage splits by operator.

Packages

We will need some R packages to do this, which should be installed already if you have been following my tutorials.

sf – Working with spatial data in R
dplyr – Maniuplating data in R
downloader – Data downloading from R
tigris – Census data in R

If none of these are installed, please go follow my tutorials! But if you’re being lazy.

install.packages(c('sf', 'dplyr', 'downloader', 'tigris'))

I will install more when I get to some more advanced mapping.

Well Data

The New Mexico Oil Conservation Division keeps a pretty up to date well shapefile. Let’s download it into R and put it into an sf dataframe.

library(sf)
library(dplyr)
library(downloader)
library(tigris)

#Download zip file of wells
downloader::download('ftp://164.64.106.6/Public/OCD/OCD%20GIS%20Data/Shapefiles/WellGIS.zip', destfile = 'WellGIS.zip', mode = 'wb')
#Unzip it to your directory
unzip('WellGIS.zip')
#Delete zip file
unlink('WellGIS.zip')

#Read into a sf dataframe
sf2 <- sf::read_sf('WellGIS.shp')

Upon inspection, it’s not perfect but it is workable for our purposes.

If anyone is reading from the OCD please include bottom-hole locations. And maybe fracture stimulation data. Also make production data easily downloadable. Thank you for attending my TED talk.

What data does it have?

Operator Info
Well Identifiers
County Info
Surface Location (Lat/Long and Township-Range-Section)
Well Status/Type
Link to main well file page
Spud Date
Measured Depth/True Vertical Depth
Reservoir

It is not laid out perfectly, so we have to do some data cleanup.

The steps I will take are:

Filter to Lea and Eddy counties, which is where the Delaware Wolfcamp/Bone Spring is located.
Identify the reservoir by using the pool_id_li column.
Clean up the Section-Township-Range column, which is identified as ulstr.
Filter to only include wells that are not P&A’d.
Rename operators to make the data more readable.

#Filter counties
sf2 <- sf2 %>% filter(county == 'Eddy'|county == 'Lea')

#Create new reservoir name column; I only really care about Wolfcamp/Bone Spring
sf2$reservoir[grepl('WOLFC', sf2$pool_id_li)] <- 'Wolfcamp'
sf2$reservoir[grepl('BONE S', sf2$pool_id_li)] <- 'Bone Spring'

#Clean up Section/Twn/Rng column
sf2$ulstr[nchar(sf2$ulstr) == 12] <- substr(sf2$ulstr[nchar(sf2$ulstr)==12], 3, 12)

#Filter to non P&A'd wells
sf2 <- sf2 %>% filter(status == 'Active'|status == 'New'|grepl('emporary', status))

#Create new operator column
sf2$operator <- sf2$ogrid_name
sf2$operator[grepl('CHEVRON', sf2$operator)] <- 'CHEVRON'
sf2$operator[grepl('Chevron', sf2$operator)] <- 'CHEVRON'
sf2$operator[grepl('DEVON', sf2$operator)] <- 'DEVON'
sf2$operator[grepl('OXY', sf2$operator)] <- 'OXY'
sf2$operator[grepl('OCCIDENTAL', sf2$operator)] <- 'OXY'
sf2$operator[grepl('XTO', sf2$operator)] <- 'XTO'
sf2$operator[grepl('CIMAREX', sf2$operator)] <- 'CIMAREX'
sf2$operator[grepl('COG O', sf2$operator)] <- 'CONCHO'
sf2$operator[grepl('COG P', sf2$operator)] <- 'CONCHO'
sf2$operator[grepl('EOG', sf2$operator)] <- 'EOG'
sf2$operator[grepl('MATADOR', sf2$operator)] <- 'MATADOR'
sf2$operator[grepl('Spur', sf2$operator)] <- 'SPUR'
sf2$operator[grepl('TAP ROCK', sf2$operator)] <- 'TAP ROCK'
sf2$operator[grepl('CONOCO', sf2$operator)] <- 'CONOCO'
sf2$operator[grepl('MARATHON', sf2$operator)] <- 'MARATHON'

That should be good for now on the well data. I will manipulate it more when we get to the analysis portion.

Federal Acreage

This time, we will be downloading a shapefile from the BLM that shows surface ownership in New Mexico.

We don’t need to do too much work, but I will append county data to the shapefile using the tigris package so I can filter just to Eddy/Lea.

#Download data
downloader::download('https://www.nm.blm.gov/shapeFiles/state_wide/SURFACE_OWN.zip', destfile = 'SURFACE_OWN.zip', mode = 'wb')
unzip('SURFACE_OWN.zip')
unlink('SURFACE_OWN.zip')

#Read in sf dataframe
surface <- read_sf('Surface_Ownership.shp')

#Grab New Mexico county data
counties <- tigris::counties(state = 'NM')
counties <- st_as_sf(counties)

#Transform both to WGS84
surface <- surface %>% st_transform(crs = 4326)
counties <- counties %>% st_transform(crs = 4326)

#Filter to Eddy/Lea
counties <- counties %>% filter(NAME == 'Lea'|NAME == 'Eddy') %>% select(county = NAME)

#Join county file and remove any na's
surface <- surface %>% st_join(counties)
surface <- surface %>% filter(!is.na(county))
surface <- surface %>% filter(!duplicated(globalid))

Analysis

Now we can get to the good stuff. We will need to install more packages, largely for cool maps (I definitely haven’t covered a lot of these before so fresh installs needed).

leaflet
tidyr
data.table
sp
rgdal
KernSmooth
raster

Federal Exposure

How much federal acreage is there? Let’s look at it with leaflet.

library(leaflet)

leaflet(options = leafletOptions(zoomControl = FALSE)) %>% 
  addProviderTiles(providers$OpenStreetMap) %>%
  addPolygons(data = counties, fillColor = 'transparent', color = 'black', weight = 2) %>%
  addPolygons(data = surface %>% filter(own == 'BLM'), fillColor = 'red', fillOpacity = 0.25, color = 'black', weight = 1)%>%
  addMiniMap(position = 'topleft')

Figure 1: BLM acreage in Eddy/Lea

In the words of my 2-year old, “Ouchie Boo Boo”. That is pretty hefty exposure, and it even infiltrates a lot of the best acreage (Southwestern Lea).

Reservoir Exposure

I can also look at both Bone Spring and Wolfcamp historical activity to help quantify impact.

wcList <- sf2 %>% filter(reservoir == 'Wolfcamp') %>% 
  st_transform(4326) %>% st_join(surface) %>% 
  filter(!duplicated(ulstr)) %>%
  mutate(own = replace(own, own != 'BLM', 'Other')) %>%
  data.frame() %>% group_by(own) %>% summarise(count = n()*640) %>%
  filter(!is.na(own))

boneSp <- sf2 %>% filter(reservoir == 'Bone Spring') %>% 
  st_transform(4326) %>% st_join(surface) %>% 
  filter(!duplicated(ulstr)) %>%
  mutate(own = replace(own, own != 'BLM', 'Other')) %>%
  data.frame() %>% group_by(own) %>% summarise(count = n()*640) %>%
  filter(!is.na(own))

What I am doing is:

Filter my original well dataframe to be either Wolfcamp or Bone Spring.
Transform coordinates to be the same projection as the Federal acreage dataframe.
Remove duplicated Twn/Rng/Sections. I do this to remove bias. If operators have concentrated a lot of their drilling in one area (which they have), it will not accurately gauge the acreage impact unless I condense it down to sections. Bias is also a huge problem with type curves, but I will leave that post for another day.
Rename non BLM acreage sections to Other
Group by acreage type and count the sections, and then multiply by 640 to get an estimate of gross acreage footprint.

Figure 2a: BLM acreage in existing Wolfcamp

Figure 2b: BLM acreage in existing Bone Spring

Looks like a roughly even split for the Bone Spring so far, while Wolfcamp is closer to 40%. 40-50% of acreage looks to be impacted in this cursory look, which is definitely concerning.

Let’s condense the data a little bit and look at individual operators by reservoir.

Wolfcamp

To guesstimate operator footprint, I need to get a little creative. I will just make an assumption that if an operator has a producing well in the section regardless of reservoir, then they own depth rights to all of it. While there is probably some depth severance, it’s likely not THAT much.

I also want to just look at areas where the Wolfcamp horizontal play is really active. There is some shallow carbonate stuff in the north that may be impacted, but it is not as concerning to me as the traditional Wolfcamp unconventional development. For this, I am going to look at some density maps.

library("data.table")
library("sp")
library("rgdal")
library("KernSmooth")
library("raster")

#Convert sf dataframe to data table after filtering to just Wolfcamp wells
dat<- data.table(sf2 %>% filter(reservoir == 'Wolfcamp'))

#Create Kernel Density Output
kde <- bkde2D(dat[ , list(longitude, latitude)],
              bandwidth=c(.0045, .0068), gridsize = c(500,500))

#Convert kde to a Raster for plot
KernelDensityRaster <- raster(list(x=kde$x1 ,y=kde$x2 ,z = kde$fhat))
KernelDensityRaster@data@values[which(KernelDensityRaster@data@values < 1)] <- NA

#Color scale for the raster object (11 bins)
palRaster <- colorBin(rev(brewer.pal(n=11, name='Spectral')),bins = 10, 
                      domain = KernelDensityRaster@data@values, na.color = "transparent")

#Create the leaflet Plot

leaflet(options = leafletOptions(zoomControl = FALSE)) %>% addTiles() %>% 
  addRasterImage(KernelDensityRaster, 
                 colors = palRaster, 
                 opacity = .8) %>%
  addLegend(pal = palRaster, 
            values = KernelDensityRaster@data@values, 
            title = "Well Density<br>Wolfcamp")%>%
  addPolygons(data = counties, fillColor = 'transparent', color = 'black', weight = 2) %>%
  addLabelOnlyMarkers(data = counties$pt1, label = as.character(counties$county),
                      labelOptions = labelOptions(noHide = T, direction = 'top', textOnly = T)) %>%
  addMiniMap(position = 'topleft')

Figure 3: Wolfcamp Well Density in Eddy/Lea

Largely, it appears to me that the bulk of activity is south of Carlsbad, so I will filter my original dataframe sf2 to look at just acreage south of there, and try to estimate the Fed/Fee split impact by operator.

library(tidyr)

#Filter to south of Carlsbad
wc1 <- sf2 %>% filter(latitude < 32.42)

#Find sections and operators in the data
operatorList <- wc1 %>% data.frame() %>%
   group_by(ulstr, operator) %>% summarise(wells = n()) %>%
   ungroup() %>% arrange(operator)

#Join the geometry object to the operatorList dataframe
locs <- sf2 %>%  filter(!duplicated(ulstr))
locs <- locs[,c('ulstr', 'geometry')]
operatorList <- operatorList %>% left_join(locs)

#Transform to WGS84
operatorList <- st_as_sf(operatorList) %>% st_transform(crs = 4326)
#Join BLM grid to dataframe
operatorList <- st_as_sf(operatorList) %>% st_join(surface)
#Remove duplicates
operatorList <- operatorList %>% filter(!duplicated(paste0(operator, ulstr)))

#Find operators with at least 40000 total acres and get split of BLM/Other
operatorList <- operatorList %>%
  mutate(own = replace(own, own != 'BLM', 'Other')) %>%
  data.frame() %>% group_by(operator, own) %>% summarise(count = n()*640) %>%
  filter(!is.na(own)) %>% group_by(operator) %>% filter(sum(count) >= 40000) %>% ungroup() %>%
  spread(own, count) %>%
  mutate(blmFrac = BLM/(BLM+Other), total = BLM+Other)

Figure 4: Wolfcamp BLM Exposure by Operator – Eddy/Lea New Mexico

It appears to me that Exxon (XTO), Devon, and Cimarex have some pretty outsized exposure in the play, so I am definitely concerned if I am them. Can this be wrong? Absolutely, but this is what the data seems to indicate. Of course, as I’ve mentioned, this is not ALL of their acreage, it’s just everything that’s got a well sitting on it.

I’ll take a look real quick at the original federal acreage map with Exxon‘s sections superimposed just to make sure I’m not misinterpreting.

Figure 5: Exxon Existing Developed Sections Exposure to BLM – Eddy/Lea New Mexico

Seems like it’s not far off.

Bone Spring

I will do something similar for the Bone Spring, though I will save you the coding and just get to the meat of the analysis.

Figure 6: Bone Spring Well Density in Eddy/Lea

Just going by view, it looks like there is a lot more active development areas for the Bone Spring, so I will expand the coverage area for this one, to anything south of Artesia.

Figure 7: Bone Spring BLM Exposure by Operator – Eddy/Lea New Mexico

Quite a few more operators pass the acreage threshold set earlier, but we do still see the top 3 Wolfcamp companies at the top of the graph (largely because the acreage is shared).

The Bone Spring view will give us a generally better idea of overall Federal exposure in the entire play, while the Wolfcamp view will be constrained more to the better economic areas of the play.

Summary

That is it for this session. While this is a quick and easy way to estimate BLM exposure, it would be a more improved analysis if you use actual acreage for each operator. However, the methodologies discussed here would still be applicable in doing this valuation.

Obviously, a frac ban will be a very big negative for the Permian, and more specifically operators in the New Mexico portion of the Delaware Basin. We’ll see what happens in November!

See my other post on Estimating The Negative Impacts of a Catchment Area Increase in Weld County, Colorado, for what can also be done using publicly available data and the mapping functionalities in R.