SimulatingMicrostructures.Rmd

title: "Simulating microstructures with gmGeostats"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{my-vignette}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}

knitr::opts_chunk$set(
  collapse = TRUE,
  echo = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 6
)

# ToDo: anistropoy: what means the ratio exactly?

library(gmGeostats)
library(ggplot2)
library(data.table) # only needed, because Solveig is faster with melting in data.table than with dplyr and Co. can be replaced by dplyr, if someone changes the code from data.table to dplyr later.

Part 1: Granite/Gneiss

Granites or Gneiss consist mainly of Quartz, Feldspars and micas, e.g. Biotite and/or Muscovite. Due to the brittle characteristics these types of rocks have often fractures or joint systems. With gmGeostats mineral composition, anisotropy of certain mineral phases and features like joints and fractures can be simulated, e.g. if simulations of crystalline rocks are required for testing code or work flows.

General idea - Turning bands algorithm

here comes some nice background from Raimon

Example for simulation

The simulation consists of five steps:

1. decide for the number of phases to be simulated. Fracture and joint systems count as a "phase", e.g. if there is one joint system, this would be one phase, if there are two joint systems this would be simulated by two additional phases.
2. define the characteristics of each phase, i.e. the crystal elongation (or void shape) and orientation/foliation
3. generate the variograms by using a function to define for each phase the model of the phase-variogram, the nugget, the sill, the range and the anistropy with respect to the range.
4. use the algorithm of Turning Bands to calculate the potential for each phase to occur.
5. Visualize the result.

In the following example, we will show the implementation for 2D. Because the code expects 3D the third dimension is always "muted" by setting it as a constant.

Number of phases

For the example let's assume five phases, Quartz (qu), K-Feldspar (fsp), Biotite (bi), a fracture system (frac) and a joints system (joints).

Generate a template matrix with only zeros, which has the dimensions n x n:

n_phases = 5

template = matrix(0, nrow = n_phases, ncol = n_phases)
phases = c("qu", "fsp", "bi", "frac", "joints")
colnames(template) <- rownames(template) <- phases

In a later step when the five variograms are stacked, these five matrices define the sill for each phase. Therefore, on the diagonal we replace the 0 with a 1 for the respective phase:

s1 = template
s1[1,1] <-1
s2 = template
s2[2,2] <-1
s3 = template
s3[3,3] <-1
s4 = template
s4[4,4] <-1
s5 = template
s5[5,5] <-1

The example for Quartz looks like this:

s1

Characteristics of phases

The anisotropy of the elongation for each phase can be obtained by the function anis2D_._par2A. The function parameter ratio provides the ratio of shortest to longest range, and should therefore be smaller than 1. E.g. ratio = 1 is a sphere (isotrop), while ratio = 0.5 would describe a mineral which grains are typicall double as long as thick. The function parameter angle gives the orientation of this anisotropy, measured in degrees counterclockwise from the east (E). E.g. angle = 45 with ratio = 0.5 would result into NE-SW elongated crystals.

For Quartz, let's assume we have mostly roundish/quadratic grains and only a minor preferred orientation of the crystals, so that they are slightly oriented along NE-SW:

p1_qu = anis_GSLIBpar2A(ratios = .9, angles = 45)

The Feldspars have a more pronounced elongated shape, also with preferred orientation along the same direction as the Quartz grains:

p2_fsp = anis_GSLIBpar2A(ratios = .6, angles = 45)

Function anis_GSLIBpar2A returns the inverse Cholesky decomposition of the anisotropy matrix. So, if the anistropy of the mineral phases had been generated by this function, then we need to re-calculate the anistropy range matrices A:

A1_qu = as.AnisotropyRangeMatrix(p1_qu)   
A2_fsp = as.AnisotropyRangeMatrix(p2_fsp)

This results into

A1_qu

for Quartz and into

A2_fsp

for Feldspar.

There is also the option to directly generate the anisotropy matrix. In this case, the matrix can be written manually: The Biotite occurs as thin plates and has a slightly different preferred orientation then the two main minerals:

A3_bi = AnisotropyRangeMatrix(
    matrix(data = c(1, .6, 0, .6, 1, 0, 0, 0, 1), ncol = 3)
  )
A3_bi

For the fracture and joint systems the anisotropy is very pronounced, because they are very thin but passing through longer distances through the rock. Accordingly, the anisotropy matrix should be chosen similar to this:

A4_frac = AnisotropyRangeMatrix(
  matrix(data = c(1, -.999, 0, -.999, 1, 0, 0, 0, 1), ncol = 3)
    )
A4_frac

And let's assume, that the joint system is orthogonal to the fracture system:

A5_joints = AnisotropyRangeMatrix(
    matrix(data = c(1, .999, 0, .999, 1, 0, 0, 0, 1), ncol = 3)
  )
A5_joints

Generate variograms

Fore each phase, variogram functions will construct the potentials that control which phase will actually occur. The function setCgram generates these variograms. It has four parameters:

• The type provides the variogram model, e.g. Gaussian or Exponential. Gaussian models provide smooth transitions to other minerals, and may be used if bands of minerals in equilibrium should be modeled. Exponential models provide rather abrupt changes between mineral phases and should be preferred model for fractal contacts, such as non-equilibrium contacts, dissolution/porosity, alteration, etc. In this example, Quarts and Feldspar are modeled by Gaussian types, and Biotite and the fractures by the Exponential type.
• The nugget gives the nugget effect, but in the case of microstructures it should be omitted.
• For the sill, the partial sill of the correlation function, the microstructure require simply the matrices with one 1 in the diagonal. We produced these matrices in the first part, where the number of phases was defined.
• The anisRages use the information from the anisotropy matrices generated in the previous part. This anisotropy range has to be downscaled, therefore a factor is applied prior to the the anisotropy matrices:

v1_qu = setCgram(type = vg.Gau, sill = s1, anisRanges = 0.0005*A1_qu)
v2_fsp = setCgram(type = vg.Gau, sill = s2, anisRanges = 0.0005*A2_fsp)
v3_bi = setCgram(type = vg.Exp, sill = s3, anisRanges = 0.005*A3_bi)
v4_frac = setCgram(type = vg.Exp, sill = s4, anisRanges = 0.05*A4_frac)
v5_joints = setCgram(type = vg.Exp, sill = s5, anisRanges = 0.05*A5_joints)

This factor controls the typical size of each phase, but in a complex, indirect way. Plotting one of these variograms shows this:

plot(v3_bi)

The variograms need now to be nested:

vm = v5_joints + v4_frac + v3_bi  + v2_fsp + v1_qu

Please note, that the order of the stacking is not important.

The visualisation of the stacked variograms looks as follows:

plot(vm)

Apply Turning Bands algorithm

In a first step, a grid of nodes has to be generated, to determine the number and extend of the pixels of the finalized "image". This can be done by expand.grid:

# make grid
resol=200
Xgrid = as.matrix(expand.grid(x = seq(0,1, length.out = resol),
                              y = seq(0,1, length.out = resol),
                              z = 0))

Pleas note: Xgrid needs to have the class matrix, and the third dimension is set to constant by using only zeros for the z-direction. However, the whole algorithm here presented is directly suitable for 3D simulations: the complexity rather lies in the representation of the outcomes.

The data are generated by:

set.seed(122487)
rock1 = gsi.TurningBands(X = Xgrid, 
                         vgram = vm, 
                         nBands = 1000)
colnames(rock1)[-(1:3)] = phases # apply the names

This gives back a large matrix which looks like this:

rock = data.frame(rock1)
head(rock)

For each 'row' in the resulting image, i.e. keep a y-value constant, and look along the x-axis, the five phases have a certain potential to appear. For example for y = 0.10025 :

rock_melt = melt(data.table(rock), id.vars = c("x", "y", "z"), 
                            variable.name = "phase", 
                            value.name = "relation")
rock_melt[, phase := factor(phase, levels = phases)]
ggplot(rock_melt[y > 0.1 & y < 0.101], aes(x = x, y = relation, color = phase)) + 
  theme_bw() +
  #geom_hline(yintercept = 0, color = "red") +
  geom_line() +
  scale_color_brewer(type = "qual", palette = "Set1")

For each pixel, i.e. for each x/y combination respectively each row in rock, the phase with maximum potential is chosen:

rocksims = apply(rock1[, -(1:3)],1, function(x) which.max(x))
dim(rocksims) = c(resol, resol)

Visualization

The resulting image is:

# somehow ggplot doesnt work...
# phase_info = data.table(n_phase = c(1:5), phase = factor(phases, levels = phases))
# rocksims_melt = melt(data.table(data.frame(rocksims), keep.rownames = T), id.vars = "rn", value.name = "n_phase")
# rocksims_melt = rocksims_melt[phase_info, on = "n_phase"]
# # plot
# ggplot(rocksims_melt, aes(y = rn, x = variable, fill = phase)) + 
#   geom_tile() + 
#   coord_equal() + 
#   scale_fill_brewer(type = "qual", palette = "Set1") + 
#   theme(axis.text = element_blank(), axis.title = element_blank(), axis.ticks = element_blank())

image(rocksims, breaks=0.5:6, col=RColorBrewer::brewer.pal(5,"Set1"), asp = 1)
legend("topright", fill=RColorBrewer::brewer.pal(5,"Set1"), legend = phases)

If one of the phases should be less pronounced, e.g. here in the example the joints are very frequent, then the potential can be regulated by weights. The weights simply increase or diminish the potential of each phase. For example:

ph_weights = c(1,1,0,0,-1) # inurease qu and fsp, diminish joints
rocksims = apply(rock1[, -(1:3)],1, function(x) which.max(x + ph_weights))
dim(rocksims) = c(resol, resol)

Then the the line is pushed down:

rock_melt = melt(data.table(rock)[, joints := joints + ph_weights[5]][, qu := qz + ph_weights[1]][, fsp := fsp + ph_weights[2]], id.vars = c("x", "y", "z"), 
                            variable.name = "phase", 
                            value.name = "relation")
rock_melt[, phase := factor(phase, levels = phases)]
ggplot(rock_melt[y > 0.1 & y < 0.101], aes(x = x, y = relation, color = phase)) + 
  theme_bw() +
  #geom_hline(yintercept = 0, color = "red") +
  geom_line() +
  scale_color_brewer(type = "qual", palette = "Set1")

and the resulting pattern looks like:

image(rocksims, breaks=0.5:6, col=RColorBrewer::brewer.pal(5,"Set1"), asp = 1)
legend("topright", fill=RColorBrewer::brewer.pal(5,"Set1"), legend = phases)

Example with more realistic fractures

If single the phases should be changed then only certain steps have to be retraced. For example, in the simulation above, the geologist would expect, that the fractures rather follow the biotite layers, instead of being orthogonal to the joints.

Changing the orientation of the fractures:

p4_frac = anis_GSLIBpar2A(ratio = .01, angle = 95)
A4_frac = tcrossprod(p4_frac)
A4_frac

Then recalculate the variogram for the fractures and stack all variograms again:

v4_frac = setCgram(type = vg.Exp, sill = s4, anisRanges = 0.05*A4_frac)
vm = v5_joints + v4_frac + v3_bi  + v2_fsp + v1_qu

Plot:

plot(vm)

Applying the algorithm of the turning bands:

set.seed(384952) 
rock2 = gsi.TurningBands(X = Xgrid, vgram = vm, nBands = 1000)
colnames(rock2)[-(1:3)] = phases

Get the phases per pixel, including the weights for the joints:

ph_weights = c(0,0,0,0,-1) # diminish only the joints
rocksims = apply(rock2[, -(1:3)],1, function(x) which.max(x + ph_weights))
dim(rocksims) = c(resol, resol)

and then visualize:

image(rocksims, breaks=0.5:6, col=RColorBrewer::brewer.pal(5,"Set1"), asp = 1)
legend("topright", fill=RColorBrewer::brewer.pal(5,"Set1"), legend = phases)

In this plot, the fracture seem to be too big. Again, this can be easily resolved by downweighing the fracturesa (and joints) and increase a bit the Quartz:

ph_weights = c(.5,0,0,-1.2,-1) # diminish joints and fracture, increase quartz
rocksims = apply(rock2[, -(1:3)],1, function(x) which.max(x + ph_weights))
dim(rocksims) = c(resol, resol)

image(rocksims, breaks=0.5:6, col=RColorBrewer::brewer.pal(5,"Set1"), asp = 1)
legend("topright", fill=RColorBrewer::brewer.pal(5,"Set1"), legend = phases)