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

Stress Geostatistical Model Generation

This notebook covers the first of two geostatistical model generation processes that are linked together by an exhaustive secondary variable (vertical stress). This notebook covers data loads from previous notebooks, exploratory variography, and development of a linear model of coregionalization.

# 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)
library(duvernaygeomechmodel)

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::raster('../data/study_raster.tif')
sv_sf = st_read('../data/sv_sf.GeoJSON') %>%
  drop_na()

pp_sf = st_read('../data/pp_sf.GeoJSON') 
pp_sf = pp_sf %>% 
  aggregate(., by=pp_sf %>% distinct(), mean) %>% 
  drop_na()

shmin_sf = st_read('../data/shmin_sf.GeoJSON') 
shmin_sf = shmin_sf %>% 
  aggregate(., by=shmin_sf %>% distinct(), mean) %>% 
  drop_na()

Variography

This section investigates the normality, anisotropy, and variography of each of the stress variables. 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)
# investigate normality

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

sv_norm = bestNormalize::orderNorm(sv_sf$svgrad_kpa_m)
sv_sf$svgrad_norm = sv_sp$svgrad_norm = sv_norm$x.t
ggdensity(sv_sf, x='svgrad_norm') +
  stat_overlay_normal_density(color='red', linetype='dashed')
moments::skewness(sv_norm$x.t)
moments::kurtosis(sv_norm$x.t)

# what is the global variance (should be 1...)
var(sv_sp$svgrad_norm)

# make an experimental variogram
sv_v = variogram(svgrad_norm~1, sv_sp)
plot(sv_v)

# check directional variograms
sv_v_dir = variogram(svgrad_norm~1, sv_sp, alpha=c(0,22,45,67,90,112,135,157))
plot(sv_v_dir)

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

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

Variogram Modelling - SV

Vertical stress shows geometric anisotropy with a major axis direction of 157 degrees, and a trend in the shorter direction. Because we are trying to fit an LMC, we try to keep variogram structures simple. A anisotropy of 75% reflects the variogram map nicely. Good enough to start.

sv_v = variogram(svgrad_norm~1, sv_sp, alpha=c(67,157)) 
plot(sv_v)

sv_vgm = vgm(0.69, "Sph", 50E3, nugget = 0.01)
sv_vgm = vgm(0.3, "Sph", 15E3, anis=c(67, 0.75), add.to=sv_vgm)
plot(sv_v, sv_vgm)

Descriptive Variography - SV

Looking at pore pressure, there is an outlier in the data because two points are way above the duvernay. A couple observations: - The low pressures towards the east are creating a binomial distribution but I don't think they should be thrown out. - The data is quite skewed (-0.82) with heavy tails (k=3.6) - The variograms are awful. I honestly think a pure nugget effect might be the best way to look at them!? - You can see the data sparsity in the variogram map - patchy with no discernable anisotropy.

pp_sp = as_Spatial(pp_sf)
# investigate normality

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

# normal transform
pp_norm = bestNormalize::orderNorm(pp_sf$ppgrad_kpa_m)
pp_sf$ppgrad_norm = pp_sp$ppgrad_norm = pp_norm$x.t
ggdensity(pp_sf, x='ppgrad_norm') +
  stat_overlay_normal_density(color='red', linetype='dashed')
moments::skewness(pp_norm$x.t)
moments::kurtosis(pp_norm$x.t)

# what is the global variance (should be 1...)
var(pp_sp$ppgrad_norm)

# make an experimental variogram
pp_v = variogram(ppgrad_norm~1, pp_sp)
plot(pp_v)

# check directional variograms
pp_v_dir = variogram(ppgrad_norm~1, pp_sp, alpha=c(0,22,45,67,90,112,135,157))
plot(pp_v_dir)

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

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

Variogram Modelling - PP

A uninformative variogram here might be best solution. We use a spherical variogram with a short range (10 km) to fit our first observation. We need a model, so this is a pretty broad assumption. The other alternative is full nugget, but that simply imparts random noise and isn't all that useful for cosimulation.

pp_v = variogram(ppgrad_norm~1, pp_sp, alpha=c(67,157))

pp_vgm = vgm(0.1, "Sph", 50E3, nugget = 0.1)
pp_vgm = vgm(0.85, "Sph", 10E3, anis=c(67, 0.75), add.to=pp_vgm)
plot(pp_v, pp_vgm)

pp_vgm

Descriptive Variography - Minimum Horizontal Stress

• Actually, shmin isn't horribly non-normal (surprisingly), with nearly no skewness (-0.10) and light tails (k=2.07). But once again (as with pp), the variogram looks like garbage and the data is very sparse.

sh_sp = as_Spatial(shmin_sf)
# investigate normality

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

# normal transform
sh_norm = bestNormalize::orderNorm(shmin_sf$shmingrad_kpa_m)
shmin_sf$shgrad_norm = sh_sp$shgrad_norm = sh_norm$x.t
ggdensity(shmin_sf, x='shgrad_norm') +
  stat_overlay_normal_density(color='red', linetype='dashed')
moments::skewness(sh_norm$x.t)
moments::kurtosis(sh_norm$x.t)

# what is the global variance (should be 1...)
var(sh_sp$shgrad_norm)

# make an experimental variogram
sh_v = variogram(shgrad_norm~1, sh_sp)
plot(sh_v)

# check directional variograms
sh_v_dir = variogram(shgrad_norm~1, sh_sp, alpha=c(0,22,45,67,90,112,135,157))
plot(sh_v_dir)

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

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

Variogram Modelling - SV

We apply the same variogram as pore pressure

sh_v = variogram(shgrad_norm~1, sh_sp, alpha=c(67,157))

sh_vgm = vgm(0.5, "Sph", 50E3, nugget = 0.1)
sh_vgm = vgm(0.5, "Sph", 15E3, anis=c(157, 0.5), add.to=sh_vgm)
plot(sh_v, sh_vgm)
sh_vgm

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 use a linear short range model with nugget and two spherical variograms. This seems to provide a relatively okay quantification of our sparse data.

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

Setup the analysis in gstat for the LMC fitting.

# setup model structures
vmodel = vgm(0, "Sph", 20E3)
vmodel = vgm(0, "Sph", 40E3, 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='shmin', formula=shgrad_norm~1, data=sh_sp, beta=0, nmax=50, model=vmodel)
g = gstat(g, id='pp', formula = ppgrad_norm~1, data=pp_sp, beta=0, nmax=50, model=vmodel)
g = gstat(g, id=c('sv','shmin'), beta=0, nmax=50, model=vmodel)
g = gstat(g, id=c('sv','pp'), beta=0, nmax=50, model=vmodel)
g = gstat(g, id=c('shmin','pp'), beta=0, nmax=50, model=vmodel)

# make variogram
v = variogram(g, cressie=TRUE)

Fit the linear model of coregionalization. Since the random initiation point of the model governs the outcome, a seed is set to fix the random number generator and allow reproducible results. We iterate from 1 to 50 to find a model to visually minimizes the loss along the SV variogram while providing reasonable fits to PP and SH. A single model is carried forward into cosimulation.

Seeds with reasonable fits (+ indicates preferred): 1, 2, 7, 8, 15, 16, 22, 23+, 27, 30, 40, 46

# setup lmc problem
set.seed(23)
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, lmc_in$vartype,
  weighting_method=6
  )

g = assign_lmc_sills(g, res$sills)

plot(v,g)

Save variogram + gstat object for simulation

norm_objs = list(sv_norm, pp_norm, sh_norm)
names(norm_objs) = c('sv','pp','sh')

data_objs = list(sv_sp, pp_sp, sh_sp)
names(data_objs) = c('sv','pp','sh')

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/stress_vgmobj.rds')

```

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