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

Log Cosimulation

This document goes over development of a LMC and cosimulation of the log data, conditioned on the vertical stress exhaustive data. Similar to the previous notebook that addressed stress data.

# geospatial packages
library(sf)
library(gstat)
library(raster)
library(sp)

# statistical normalization and analysis
library(ggpubr)
library(moments)
library(bestNormalize)

# data munging & plotting
library(tidyverse)
library(viridis)

source('../R/geospatial.R')
source('../R/geostats.R')

Load spatial dataframes

See 'geostats_data_prep.Rmd' to generate these files and dump into the cache. Read GeoJSONs from the cache.

We run an aggregation function to take the mean gradient across measurement points because a) we aren't doing 3D kriging, and b) because geostatistics doesn't work with zero distance variance that isn't equal to zero

study_raster = raster('../data/study_raster.tif')
sv_sf = st_read('../data/sv_sf.GeoJSON')
logs_sf = st_read('../data/facies_data.GeoJSON')
shale_sf = logs_sf %>% filter(facies == 'dv_shale')
shale_sp = as_Spatial(shale_sf)
carb_sf = logs_sf %>% filter(facies == 'dv_carb')
carb_sp = as_Spatial(carb_sf)

Variography

This section investigates the normality, anisotropy, and variography of each of the log variables (bulk density, vp, vs). We also repeat the SV variography minus the exploration (see geostats_stress.Rmd)

Since Sequential Gaussian Cosimulation is used, we transform variables to normal distributions using the ordered quantile technique

Descriptive Variography - SV

sv_sp = as_Spatial(sv_sf)
sv_sp = sv_sp[!is.na(sv_sp$svgrad_kpa_m),]
sv_norm = bestNormalize::orderNorm(sv_sp$svgrad_kpa_m)
sv_sp$svgrad_norm = sv_sp$svgrad_norm = sv_norm$x.t
sv_v = variogram(svgrad_norm~1, sv_sp)
sv_vgm = vgm(0.8, "Sph", 40E3, nugget = 0.01)
sv_vgm = vgm(0.2, "Sph", 10E3, anis=c(67, 0.75), add.to=sv_vgm)
plot(sv_v, sv_vgm)

Descriptive Variography - RHOB / Shale

• Density actually has a relatively well behaved variogram
• No anisotropy noted
• Good coverage
• Nothing of note in variogram map

shale_rhob = shale_sp[!is.na(shale_sp$rhob_mean),]
# investigate normality

# make anisotropic variograms
ggdensity(shale_sf, x='rhob_mean') +
  stat_overlay_normal_density(color='red', linetype='dashed')
moments::skewness(shale_sf$rhob_mean)
moments::kurtosis(shale_sf$rhob_mean)

# normal transform
shale_rhob_norm = bestNormalize::orderNorm(shale_rhob$rhob_mean)
shale_rhob$rhob_norm = shale_rhob_norm$x.t
ggdensity(as.data.frame(shale_rhob), x='rhob_norm') +
  stat_overlay_normal_density(color='red', linetype='dashed')
moments::skewness(shale_rhob_norm$x.t)
moments::kurtosis(shale_rhob_norm$x.t)

# what is the global variance (should be 1...)
var(shale_rhob$rhob_norm)

# make an experimental variogram
shale_rhob_v = variogram(rhob_norm~1, shale_rhob)
plot(shale_rhob_v)

# check directional variograms
rhob_v_dir = variogram(rhob_norm~1, shale_rhob, alpha=c(0,22,45,67,90,112,135,157))
plot(rhob_v_dir)

# check for geometric anisotropy
ggvarmap = ggvariogram_map(
  formula = rhob_norm~1, 
  sp_df = shale_rhob,
  cutoff = 100E3, 
  width = 5E3,
  threshold = 1
  )

ggsave('../output/ggvarmap_rhob_shale.pdf', ggvarmap)

Descriptive Variography - RHOB / Carb

• The carbonate density appears to be almost all nugget (60% of variance?)
• Some short range predictability is noted
• Some anisotropy is noted

carb_rhob = carb_sp[!is.na(carb_sp$rhob_mean),]
# investigate normality

# make anisotropic variograms
ggdensity(carb_sf, x='rhob_mean') +
  stat_overlay_normal_density(color='red', linetype='dashed')
moments::skewness(carb_sf$rhob_mean)
moments::kurtosis(carb_sf$rhob_mean)

# normal transform
carb_rhob_norm = bestNormalize::orderNorm(carb_rhob$rhob_mean)
carb_rhob$rhob_norm = carb_rhob_norm$x.t
ggdensity(as.data.frame(carb_rhob), x='rhob_norm') +
  stat_overlay_normal_density(color='red', linetype='dashed')
moments::skewness(carb_rhob_norm$x.t)
moments::kurtosis(carb_rhob_norm$x.t)

# what is the global variance (should be 1...)
var(carb_rhob$rhob_norm)

# make an experimental variogram
carb_rhob_v = variogram(rhob_norm~1, carb_rhob)
plot(carb_rhob_v)

# check directional variograms
rhob_v_dir = variogram(rhob_norm~1, carb_rhob, alpha=c(0,22,45,67,90,112,135,157))
plot(rhob_v_dir)

# check for geometric anisotropy
ggvariogram_map(
  formula = rhob_norm~1, 
  sp_df = carb_rhob,
  cutoff = 100E3, 
  width = 5E3,
  threshold = 1
  )

ggsave('../output/ggvarmap_rhob_carb.pdf', ggvarmap)

Descriptive Variography - DTC / Shale

• Density actually has a relatively well behaved variogram
• No anisotropy noted
• Good coverage
• Nothing of note in variogram map

shale_dtc = shale_sp[!is.na(shale_sp$dtc_mean),]
# investigate normality

# make anisotropic variograms
ggdensity(shale_sf, x='dtc_mean') +
  stat_overlay_normal_density(color='red', linetype='dashed')
moments::skewness(shale_dtc$dtc_mean)
moments::kurtosis(shale_dtc$dtc_mean)

# normal transform
shale_dtc_norm = bestNormalize::orderNorm(shale_dtc$dtc_mean)
shale_dtc$dtc_norm = shale_dtc_norm$x.t
ggdensity(as.data.frame(shale_dtc), x='dtc_norm') +
  stat_overlay_normal_density(color='red', linetype='dashed')
moments::skewness(shale_dtc_norm$x.t)
moments::kurtosis(shale_dtc_norm$x.t)

# what is the global variance (should be 1...)
var(shale_dtc$dtc_norm)

# make an experimental variogram
shale_dtc_v = variogram(dtc_norm~1, shale_dtc)
plot(shale_dtc_v)

# check directional variograms
shale_dtc_v_dir = variogram(dtc_norm~1, shale_dtc, alpha=c(0,22,45,67,90,112,135,157))
plot(shale_dtc_v_dir)

# check for geometric anisotropy
ggvariogram_map(
  formula = dtc_norm~1, 
  sp_df = shale_dtc,
  cutoff = 100E3, 
  width = 5E3,
  threshold = 1
  )

ggsave('../output/ggvarmap_dtc_shale.pdf', ggvarmap)

Descriptive Variography - DTC / Carb

• The carbonate DTC appears shorter range than shale
• Some anisotropy is noted, but nothing enough to be significant or justify the use of anisotropic variograms
• The range of 16 km from Carb/Rhob matches really nicely with the DTC.

carb_dtc = carb_sp[!is.na(carb_sp$dtc_mean),]
# investigate normality

# make anisotropic variograms
ggdensity(carb_sf, x='dtc_mean') +
  stat_overlay_normal_density(color='red', linetype='dashed')
moments::skewness(carb_dtc$dtc_mean)
moments::kurtosis(carb_dtc$dtc_mean)

# normal transform
carb_dtc_norm = bestNormalize::orderNorm(carb_dtc$dtc_mean)
carb_dtc$dtc_norm = carb_dtc_norm$x.t
ggdensity(as.data.frame(carb_dtc), x='dtc_norm') +
  stat_overlay_normal_density(color='red', linetype='dashed')
moments::skewness(carb_dtc_norm$x.t)
moments::kurtosis(carb_dtc_norm$x.t)

# what is the global variance (should be 1...)
var(carb_dtc$dtc_norm)

# make an experimental variogram
carb_dtc_v = variogram(dtc_norm~1, carb_dtc)
plot(carb_dtc_v)

# check directional variograms
dtc_v_dir = variogram(dtc_norm~1, carb_dtc, alpha=c(0,22,45,67,90,112,135,157))
plot(dtc_v_dir)

# check for geometric anisotropy
ggvariogram_map(
  formula = dtc_norm~1, 
  sp_df = carb_dtc,
  cutoff = 100E3, 
  width = 5E3,
  threshold = 1
  )

ggsave('../output/ggvarmap_dtc_carb.pdf', ggvarmap)

Descriptive Variography - DTS / Shale

• DTS shale is pretty damn normal
• We only have 24 observations here... low numbers, result in some strange variograms
• We use the same scale as used on the DTC shale (40km), with a decent nugget

shale_dts = shale_sp[!is.na(shale_sp$dts_mean),]
# investigate normality

# make anisotropic variograms
ggdensity(shale_sf, x='dts_mean') +
  stat_overlay_normal_density(color='red', linetype='dashed')
moments::skewness(shale_dts$dts_mean)
moments::kurtosis(shale_dts$dts_mean)

# normal transform
shale_dts_norm = bestNormalize::orderNorm(shale_dts$dts_mean)
shale_dts$dts_norm = shale_dts_norm$x.t
ggdensity(as.data.frame(shale_dts), x='dts_norm') +
  stat_overlay_normal_density(color='red', linetype='dashed')
moments::skewness(shale_dts_norm$x.t)
moments::kurtosis(shale_dts_norm$x.t)

# what is the global variance (should be 1...)
var(shale_dts$dts_norm)

# make an experimental variogram
shale_dts_v = variogram(dts_norm~1, shale_dts)
plot(shale_dts_v)

# check directional variograms
shale_dts_v_dir = variogram(dts_norm~1, shale_dts, alpha=c(0,22,45,67,90,112,135,157))
plot(shale_dts_v_dir)

# check for geometric anisotropy
ggvariogram_map(
  formula = dts_norm~1, 
  sp_df = shale_dts,
  cutoff = 100E3, 
  width = 5E3,
  threshold = 1
  )

ggsave('../output/ggvarmap_dts_shale.pdf', ggvarmap)

Descriptive Variography - dts / Carb

• Carbonate DTS is skewed towards the lower end of the scale
• We only have 17 observations...
• Once again, we see a smaller range (16 km should do fine)
• With very limited information, we start with the same variogram as we use for DTC/Carb.

carb_dts = carb_sp[!is.na(carb_sp$dts_mean),]
# investigate normality

# make anisotropic variograms
ggdensity(carb_sf, x='dts_mean') +
  stat_overlay_normal_density(color='red', linetype='dashed')
moments::skewness(carb_dts$dts_mean)
moments::kurtosis(carb_dts$dts_mean)

# normal transform
carb_dts_norm = bestNormalize::orderNorm(carb_dts$dts_mean)
carb_dts$dts_norm = carb_dts_norm$x.t
ggdensity(as.data.frame(carb_dts), x='dts_norm') +
  stat_overlay_normal_density(color='red', linetype='dashed')
moments::skewness(carb_dts_norm$x.t)
moments::kurtosis(carb_dts_norm$x.t)

# what is the global variance (should be 1...)
var(carb_dts$dts_norm)

# make an experimental variogram
carb_dts_v = variogram(dts_norm~1, carb_dts)
plot(carb_dts_v)

# check directional variograms
dts_v_dir = variogram(dts_norm~1, carb_dts, alpha=c(0,22,45,67,90,112,135,157))
plot(dts_v_dir)

# check for geometric anisotropy
ggvariogram_map(
  formula = dts_norm~1, 
  sp_df = carb_dts,
  cutoff = 100E3, 
  width = 5E3,
  threshold = 1
  )

ggsave('../output/ggvarmap_dts_carb.pdf', ggvarmap)

Cosimulation

We are going to cosimulate stress values using a Linear Model of Coregionalization (LMC) via gstat. We could get fancier, but that's not what is important, what is important is that we have heterotopic data underpinned by an exhaustive secondary variable and need to quantify the uncertainty in our stress measurements, and how that compares to our model uncertainty.

We first try with two variables to develop intuition about cross-variogram fitting

# make a prediction grid
st_grid <- rasterToPoints(study_raster, spatial = TRUE)
gridded(st_grid) <- TRUE
st_grid <- as(st_grid, "SpatialPixels")

For the logs we have SV, RHOB, DTC, DTS. SV is exhaustive and the same for both carbonate and shale (the exhaustive link raster). RHOB, DTC, DTS are simulated separately for carbonate and shale

We use the same structures as employed for the stress model for consistency, despite some improvement being seen from adding a mid-range exponential model and adding more variograms. We believe this increases the parsimony of the approach and reduces bias, given all the other uncertainties we are addressing.

Shale first

# setup model structures
vmodel = vgm(0, "Sph", 40E3)
vmodel = vgm(0, "Sph", 20E3, add.to=vmodel)
vmodel = vgm(0, "Lin", 10E3, add.to=vmodel)
vmodel = vgm(0, "Nug", 0, add.to=vmodel)

# setup problem in gstat
g = gstat(id='sv', formula=svgrad_norm~1, data=sv_sp, beta=0, nmax=50, model=vmodel)
g = gstat(g, id='rhob', formula=rhob_norm~1, data=shale_rhob, beta=0, nmax=50, model=vmodel)
g = gstat(g, id='dtc', formula = dtc_norm~1, data=shale_dtc, beta=0, nmax=50, model=vmodel)
g = gstat(g, id='dts', formula = dts_norm~1, data=shale_dts, beta=0, nmax=50, model=vmodel)
g = gstat(g, id=c('sv','rhob'), model=vmodel, beta=0, nmax=50)
g = gstat(g, id=c('sv','dtc'), model=vmodel, beta=0, nmax=50)
g = gstat(g, id=c('sv','dts'), model=vmodel, beta=0, nmax=50)
g = gstat(g, id=c('rhob','dtc'), model=vmodel, beta=0, nmax=50)
g = gstat(g, id=c('rhob','dts'), model=vmodel, beta=0, nmax=50)
g = gstat(g, id=c('dtc','dts'), model=vmodel, beta=0, nmax=50)

v = variogram(g, cressie=TRUE)

The ad-hoc fitting process here focuses on the following parameters in order: 1. SV (exhaustive) 2. RHOB 3. VP (DTC) 4. VS (DTS) 5. Cross-variograms

Seeds with good fits (+ indicates preferred): 3, 4, 6+, 9, 15, 19, 22, 25, 32, 33, 48

# setup lmc problem
set.seed(6)
lmc_in = make_lmc_input(g, vmodel)

# run the lmc algorithm
res <- fit_linear_model_of_coregionalization(
  lmc_in$azimuth, lmc_in$dip, lmc_in$nlag, lmc_in$tail, 
  lmc_in$head, lmc_in$gam, lmc_in$model, vartype=lmc_in$vartype,
  weighting_method=6
  )

g = assign_lmc_sills(g, res$sills)

plot(v,g)

Save geostatistical models

norm_objs = list(sv_norm, shale_rhob_norm, shale_dtc_norm, shale_dts_norm)
names(norm_objs) = c('sv','rhob','dtc', 'dts')

data_objs = list(sv_sp, shale_rhob, shale_dtc, shale_dts)
names(data_objs) = c('sv','rhob','dtc','dts')

vgm_obj = list(v,g,lmc_in, norm_objs, data_objs)
names(vgm_obj) = c('v','g','lmc_in','norm_objs', 'data_objs')
saveRDS(vgm_obj, '../output/shale_logs_vgmobj.rds')

Carbonate second

# setup model structures
vmodel = vgm(0, "Sph", 40E3)
vmodel = vgm(0, "Sph", 20E3, add.to=vmodel)
vmodel = vgm(0, "Lin", 10E3, add.to=vmodel)
vmodel = vgm(0, "Nug", 0, add.to=vmodel)

# setup problem in gstat
g = gstat(id='sv', formula=svgrad_norm~1, data=sv_sp, beta=0, nmax=50, model=vmodel)
g = gstat(g, id='rhob', formula=rhob_norm~1, data=carb_rhob, beta=0, nmax=50, model=vmodel)
g = gstat(g, id='dtc', formula = dtc_norm~1, data=carb_dtc, beta=0, nmax=50, model=vmodel)
g = gstat(g, id='dts', formula = dts_norm~1, data=carb_dts, beta=0, nmax=50, model=vmodel)
g = gstat(g, id=c('sv','rhob'), model=vmodel, beta=0, nmax=50)
g = gstat(g, id=c('sv','dtc'), model=vmodel, beta=0, nmax=50)
g = gstat(g, id=c('sv','dts'), model=vmodel, beta=0, nmax=50)
g = gstat(g, id=c('rhob','dtc'), model=vmodel, beta=0, nmax=50)
g = gstat(g, id=c('rhob','dts'), model=vmodel, beta=0, nmax=50)
g = gstat(g, id=c('dtc','dts'), model=vmodel, beta=0, nmax=50)

v = variogram(g, cressie=TRUE)

Seeds with reasonable fits (+ indicates preferred): 3, 4, 9, 19, 25+, 28, 33, 36+, 40

# setup lmc problem
set.seed(25)
lmc_in = make_lmc_input(g, vmodel)

# run the lmc algorithm
res <- fit_linear_model_of_coregionalization(
  lmc_in$azimuth, lmc_in$dip, lmc_in$nlag, lmc_in$tail, 
  lmc_in$head, lmc_in$gam, lmc_in$model, vartype=lmc_in$vartype,
  weighting_method=6
  )

g = assign_lmc_sills(g, res$sills)

plot(v,g)

Save geostatistical models

norm_objs = list(sv_norm, carb_rhob_norm, carb_dtc_norm, carb_dts_norm)
names(norm_objs) = c('sv','rhob','dtc', 'dts')

data_objs = list(sv_sp, carb_rhob, carb_dtc, carb_dts)
names(data_objs) = c('sv','rhob','dtc','dts')

vgm_obj = list(v,g,lmc_in, norm_objs, data_objs)
names(vgm_obj) = c('v','g','lmc_in','norm_objs', 'data_objs')
saveRDS(vgm_obj, '../output/geostatistical_models/carb_logs_vgmobj.rds')

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