In ScottHMcKean/duvernay_geomech_model: A Geostatistical Model of Geomechanical Properties in the Duvernay Play

This notebook loads the .las files from the raw data folder, intersects them with the top and base of the Duvernay, and saves the logs back to data-raw folder. Upscaled point values will be provided for reproducibility once the Duvernay B Carbonate is identified.

library(tidyverse)
library(janitor)
library(data.table)
library(raster)
library(sf)
library(plotly)
library(duvernaygeomechmodel)

# From geospatial_layers_prep.Rmd
study_bounds = st_read('../data/study_bounds.GeoJSON')
study_proj = st_crs(study_bounds)$proj4string
study_raster = raster('../data/study_raster.tif')
study_dv_bounds = st_read('../data/study_dv_bounds.GeoJSON')
surface_locs = st_read('../data/study_surf_locs.GeoJSON')
dv_top = raster('../data/dv_top_masl.tif')
dv_base = raster('../data/dv_base_masl.tif')
dv_iso = raster('../data/dv_iso_m.tif')

Pull up las files

Due to copyright reasons, we cannot provide the las files. These are stored in a data-raw folder which is omitted from the repository (sorry)

las_file_paths <- list.files('../data-raw/las_files/', full.names = TRUE)
las_file_names <- list.files('../data-raw/las_files/', full.names = FALSE)

las_file_uwis <- las_file_names %>%
    str_replace_all(., '.las', '') %>%
    str_replace_all(., 'w', 'W')

Get unique study locations

We take mean surface location for point locations, since there are multiple coordinates for many uwis. Select pertinent columns and make sure the projection is right.

study_locs = surface_locs %>% 
  filter((surface_locs$uwi %in% las_file_uwis)) %>%
  dplyr::select(uwi, latitude, longitude, kbe, ground_elev, geometry) %>% 
  group_by(uwi) %>% 
  summarize_all(mean) %>%
  st_simplify()

# verify projections
ggplot() +
  geom_sf(data=study_dv_bounds) +
  geom_sf(data=study_locs)

Extract the Duvernay Top and Base

We lose 4 wells on the outside of the study area because the AGS layer doesn't have values for them, but this shouldn't affect the results too much, and ordinary kriging created way too much variance (~500 m) when trying to predict the tops, so it wasn't feasible to interpolate the top and base.

study_points <- study_locs %>%
  dplyr::select(uwi,geometry) %>%
  as(., 'Spatial')

tops <- study_locs %>%
  mutate(dv_top = raster::extract(dv_top, study_points, method='bilinear')) %>%
  mutate(dv_base = raster::extract(dv_base, study_points, method='bilinear')) %>%
  drop_na()

Get Logs

Okay, so we have 104 logs in our study area with tops. Not bad. Now we take those logs, subset using the AER 3D model (duvernay top and base).

study_las_file_paths = las_file_paths[las_file_uwis %in% tops$uwi]

This function retrieves all the log data and take quite a while. Don't run it.. We instead load the las_df from file (but can't provide this data sorry)

las_dfs <- map(.x = study_las_file_paths, get_log_data, study_locs, dv_base, dv_top)

log_names <- map(.x = las_dfs, .f = colnames) %>%
  unlist() %>%
  unique()

comb_las_df <- bind_rows(las_dfs) %>%
  group_by(uwi) %>%
  mutate(rel_depth_m = depth_m - min(depth_m)) %>%
  dplyr::ungroup()

colnames(comb_las_df) <- colnames(comb_las_df) %>% str_replace(., '_qcgeo', '')

comb_las_df %>% fwrite(., '../data-raw/logs_df.csv')

Do some basic filtering of logs here to: 1. Replace NAs with zero for proper filtering 2. remove badhole observations 3. filter very low DTC values (< 100 us/m) 4. filter very high DTC values (> 350 us/m)

The DTC values were causing some major non-physical results downstream, so we filter them using the IQR outlier method. This has the added effect that when we remove these observations, outlying DTS / RHOB values are also removed.

comb_las_df <- fread('../data-raw/logs_df.csv') %>%
  mutate_at(c("badhole","dtc_us_m","dtc_us_m","gr_gapi"), ~replace(., is.na(.), 0)) %>%
  filter(badhole == 0) %>%
  filter(dtc_us_m > 100 & dtc_us_m < 350)

log_summary = comb_las_df %>%
  group_by(uwi) %>%
  dplyr::summarise(rows = n(), 
            rhob = sum(na.omit(rhob_k_m) < 0),
            badhole = sum(na.omit(badhole) > 0),
            dtc = sum(na.omit(dtc_us_m) < 0),
            dts = sum(na.omit(dts_us_m) < 0),
            gr = sum(na.omit(gr_gapi) < 0)
            )

Manually ID the Duvernay Carbonate so we can get the Upper Duvernay and Duvernay Carbonate for geostatistical simulation. Log these picks in data/dv_b_tops.csv

plot_ly(
  comb_las_df %>% filter(uwi=="100020906418W500"), x = ~elev_masl, y = ~dtc_us_m, 
  name = 'DTC', type = 'scatter', mode = 'lines', orientation='h') %>% 
  add_trace(y = ~gr_gapi, name = 'GR', mode = 'lines+markers')

Read in the top picks and intersect the logs with each top and bottom for further isolation and flattening across the study area. Now aggregation of the logs into 2D summaries is straightforward.

dv_carb_tops = read_csv('../data/dv_carb_picks.csv')

facies = map_df(
    unique(dv_carb_tops$uwi), 
    classify_log_intervals, comb_las_df, dv_carb_tops) %>%
  filter(facies != 'na') %>%
  group_by(uwi, facies) %>%
  summarise(
    mean_elev_masl = mean(elev_masl),
    rows = n(), 
    thickness = max(rel_depth_m) - min(rel_depth_m),
    rhob_mean = mean(rhob_k_m, na.rm=TRUE),
    rhob_sd = sd(rhob_k_m, na.rm=TRUE),
    dtc_mean= mean(dtc_us_m, na.rm=TRUE),
    dtc_sd = sd(dtc_us_m, na.rm=TRUE),
    dts_mean = mean(dts_us_m, na.rm=TRUE),
    dts_sd = sd(dts_us_m, na.rm=TRUE),
    gr_mean = mean(gr_gapi, na.rm=TRUE),
    gr_sd = sd(gr_gapi, na.rm=TRUE))

facies %>% write_csv('../data/facies_values.csv')

In addition to saving the raw log values, we also want a geometry-aware file so we can load the data plus the locations. So we take the simplified study locations and merge with the facies summary, writing as as GeoJSON.

facies %>% merge(study_locs, by ='uwi') %>% 
  st_write('../data/facies_data.GeoJSON', delete_dsn=TRUE)

ScottHMcKean/duvernay_geomech_model documentation built on Sept. 12, 2022, 10:08 a.m.