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

This notebook loads the raw stress data from Shen et al. (2018)'s publication, produces some figures for the paper, and outputs GeoJSON files for the geostatistics workflow.

library(sf)
library(tidyverse)
library(janitor)
library(plot3D)
library(viridis)
library(raster)
library(duvernaygeomechmodel)

# From geospatial_layers_prep.Rmd
setwd('notebooks')
study_bounds = st_read('../data/study_bounds.GeoJSON')
study_proj = st_crs(study_bounds)$proj4string
study_raster = raster::raster('../data/study_raster.tif')
study_dv_bounds = st_read('../data/study_dv_bounds.GeoJSON')
ground_surface = raster::raster('../data/ground_elev_masl.tif')

Load vertical stress data. Because the datum is provided in masl, we need depth below ground surface to calculate a proper gradient. We use the raster package to extract this from AGS's 3D geomodel using

sv_sf <- read_tsv('../data-raw/shen2018/DIG_2018_0013_vertical_stress.txt') %>%
  st_as_sf(.,coords = c('Long_NAD83','Lat_NAD83'), crs = 4269) %>%
  st_transform(crs = 26911) %>%
  st_crop(study_bounds) %>%
  janitor::clean_names() %>%
  rename(datum_elev_masl = datum_el) %>%
  dplyr::mutate(ground_elev_masl = raster::extract(ground_surface, .,)) %>%
  dplyr::mutate(datum_mbgs = ground_elev_masl - datum_elev_masl) %>%
  mutate(sv_mpa = vertical_stress/1000) %>%
  mutate(svgrad_kpa_m = vertical_stress/datum_mbgs) %>%
  dplyr::select(datum_elev_masl, datum_mbgs, sv_mpa, svgrad_kpa_m)

st_write(sv_sf, "../data/sv_sf.GeoJSON", delete_dsn=TRUE)

Show 2D and 3D gradients using viridis scale. The 3D plot demonstrates vertical variability across study area and shows why stress gradient is so important.

*note - needed to manually save 3D plots

ggplot(sv_sf) +
  geom_sf(data = study_dv_bounds) +
  geom_sf(data = study_bounds, fill = NA) +
  geom_sf(aes(colour=svgrad_kpa_m)) +
  theme_minimal() +
  scale_colour_viridis(option='viridis') +
  coord_sf(datum=st_crs(26911)) +
  ggsave('../output/svgrad_2d.pdf')

make_3d_plot(
  sv_sf, 
  z_col='datum_elev_masl', 
  col_col = 'svgrad_kpa_m',
  col_func = viridis(50),
  size=0.35
  )

Load Shen et al. (2018)'s shmin data. We assume a kelly bushing height of 7 m, which should be pretty fair and at least consistent across the study.

shmin_sf <- read_tsv('../data-raw/shen2018/DIG_2018_0013_DFITSh.txt') %>%
  st_as_sf(.,coords = c('Long_NAD83','Lat_NAD83'), crs = 4269) %>%
  st_transform(crs= 26911) %>%
  dplyr::filter(st_contains(study_bounds, ., sparse = FALSE)) %>%
  janitor::clean_names() %>%
  dplyr::mutate(shmin = ifelse(
    is.na(sh_reported), 
    (p_vs_log_t_high + p_vs_log_t_low)/2,
    sh_reported)) %>%
  dplyr::mutate(datum_elev_masl = kb_el_masl - datum_tvd) %>%
  dplyr::mutate(datum_mbgs = datum_tvd - 7) %>%
  dplyr::mutate(uwi = str_remove_all(well_id,"-")) %>%
  dplyr::mutate(uwi = str_remove_all(uwi,"/")) %>%
  dplyr::mutate(shmin_mpa = shmin/1000) %>%
  dplyr::mutate(shmingrad_kpa_m = shmin / (datum_tvd - 7)) %>%
  dplyr::select(datum_elev_masl, datum_mbgs, shmin_mpa, shmingrad_kpa_m)

st_write(shmin_sf, "../data/shmin_sf.GeoJSON", delete_dsn = TRUE)

Load Shen's 2018 pore pressure data. Following their approach, we take static gradients when available, followed by the mean of the buildup analysis methods, followed by the initial reported pore pressure.

pp_sf <- read_tsv('../data-raw/shen2018/DIG_2018_0013_pore_pressure.txt') %>%
  st_as_sf(.,coords = c('Long_NAD83','Lat_NAD83'), crs = 4269) %>%
  st_transform(crs= 26911) %>%
  dplyr::filter(st_contains(study_bounds, ., sparse = FALSE)) %>%
  janitor::clean_names() %>%
  dplyr::mutate(pp = case_when(
    !is.na(pore_static_gradient) ~ pore_static_gradient,
    test_type != 'static gradient' ~ mean(
      c(buildup_linear,buildup_radial,buildup_bi_linear),na.rm=TRUE),
    TRUE ~ reported_pore_initial,
  )) %>%
  dplyr::mutate(datum_elev_masl = kb_el_masl - datum_tvd) %>%
  dplyr::mutate(datum_mbgs = datum_tvd - 7) %>%
  dplyr::mutate(uwi = str_remove_all(well_id,"-")) %>%
  dplyr::mutate(uwi = str_remove_all(uwi,"/")) %>%
  dplyr::mutate(pp_mpa = pp/1000) %>%
  dplyr::mutate(ppgrad_kpa_m = pp/(datum_tvd-7)) %>%
  dplyr::select(datum_elev_masl, datum_mbgs, pp_mpa, ppgrad_kpa_m)

st_write(pp_sf, "../data/pp_sf.GeoJSON", delete_dsn=TRUE)

Show 2D and 3D gradients using magma scale. Really shows the sparsity of the stress data. We overlay the pp/shmin on the other plot for reference to show the observation heterotopy.

ggplot(shmin_sf) +
  geom_sf(data = study_dv_bounds) +
  geom_sf(data = study_bounds, fill = NA) +
  geom_sf(data=pp_sf, colour='black', shape=21, size=3) +
  geom_sf(aes(colour=shmingrad_kpa_m),size=1.5) +
  theme_minimal() +
  scale_colour_viridis(option='magma') +
  coord_sf(datum=st_crs(26911)) +
  ggsave('../output/shmingrad_2d.pdf')

make_3d_plot(
  shmin_sf, 
  z_col='datum_elev_masl', 
  col_col = 'shmingrad_kpa_m',
  col_func = magma(50),
  size=1
  )

Show 2D and 3D gradients using magma scale. Really shows the sparsity of the pore pressure data.

ggplot(pp_sf) +
  geom_sf(data = study_dv_bounds) +
  geom_sf(data = study_bounds, fill = NA) +
  geom_sf(data=shmin_sf, colour='black', shape=21, size=3) +
  geom_sf(aes(colour=ppgrad_kpa_m), size=1.5) +
  theme_minimal() +
  scale_colour_viridis(option='plasma') +
  coord_sf(datum=st_crs(26911)) +
  ggsave('../output/ppgrad_2d.pdf')

make_3d_plot(
  pp_sf, 
  z_col='datum_elev_masl', 
  col_col = 'ppgrad_kpa_m',
  col_func = plasma(50),
  size=1
  )

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