4 Stochastic Reconstruction of Microstructures

4.1 Description of Experimental Sample

[218.0.2] The experimental sample, denoted as $\mathbb{S}_{{\sf EX}}$ , is a threedimensional microtomographic image of Fontainebleau sandstone. [218.0.3] This sandstone is a popular reference standard because of its chemical, crystallographic and microstructural simplicity [14, 13]. [218.0.4] Fontainebleau sandstone consists of monocrystalline quartz grains that have been eroded for long periods before being deposited in dunes along the sea shore during the Oligocene, roughly 30 million years ago. [218.0.5] It is well sorted containing grains of around $200\ \mu\text{m}$ in diameter. [218.0.6] The sand was cemented by silica crystallizing around the grains. [218.0.7] Fontainebleau sandstone exhibits intergranular porosity ranging from $0.03$ to roughly $0.3$ [13].

Table 2: Overview of geometric properties of the four microstructures displayed in Figures 1 through 4

Properties	$\mathbb{S}_{{\sf EX}}$	$\mathbb{S}_{{\sf DM}}$	$\mathbb{S}_{{\sf GF}}$	$\mathbb{S}_{{\sf SA}}$
$M_{1}$	300	255	256	256
$M_{2}$	300	255	256	256
$M_{3}$	299	255	256	256
$\phi(\mathbb{P}\cap\mathbb{S})$	0.1355	0.1356	0.1421	0.1354
$\phi _{2}(\mathbb{P}\cap\mathbb{S})$	10.4 mm ${}^{{-1}}$	10.9 mm ${}^{{-1}}$	16.7 mm ${}^{{-1}}$	11.06 mm ${}^{{-1}}$
$L^{*}$	35	25	23	27
$1-\lambda _{c}(0.1355,L^{*})$	0.0045	0.0239	0.3368	0.3527

[218.1.1] The computer assisted microtomography was carried out on a micro-plug drilled from a larger original core. [218.1.2] The original core from which the micro-plug was taken had a measured porosity of $0.1484$ , permability of $1.3D$ and formation factor $22.1$ . [218.1.3] The porosity $\phi(\mathbb{S}_{{\sf EX}})$ of the microtomographic data set is only $0.1355$ (see Table 2). [218.1.4] The difference between the porosity of the original core and that of the final data set is due to the heterogeneity of the sandstone and to the difference in sample size. [218.1.5] The experimental sample is referred to as EX in the following. [218.1.6] The pore space of the experimental sample is visualized in Figure 1.

[page 219, §0]

Figure 1: Sample EX: Threedimensional pore space $\mathbb{P}_{{{\sf EX}}}$ of Fontainebleau sandstone. The resolution of the image is $a=7.5\ \mu\text{m}$ the sample dimensions are $M_{1}=300$ , $M_{2}=300$ , $M_{3}=299$ . The pore space is indicated opaque, the matrix space is transparent. The right image shows the front plane of the sample as a twodimensional thin section (pore space black, matrix grey).

[page 220, §1]

4.2 Sedimentation, Compaction and Diagenesis Model

[220.1.1] Fontainebleau sandstone is the result of complex physical, chemical and geological processes known as sedimentation, compaction and diagenesis. [220.1.2] It is therefore natural to model these processes directly rather than trying to match general geometrical characteristics. [220.1.3] This conclusion was also obtained from local porosity theory for the cementation index in Archie’s law [27]. [220.1.4] The diagenesis model abbreviated as DM in the following, attempts model the main geological sandstone-forming processes [4, 48].

[220.2.1] In a first step porosity, grain size distribution, a visual estimate of the degree of compaction, the amount of quartz cement and clay contents and texture are obtained by image analysis of backscattered electron/cathodo-luminescence images made from thin sections. [220.2.2] The sandstone modeling is then carried out in three main steps: grain sedimentation, compaction and diagenesis described in detail in [4, 48].

[220.3.1] Sedimentation begins by measuring the grain size distribution using an erosion-dilation algorithm. [220.3.2] Then spheres with random diameters are picked randomly according to the grain size distribution. [220.3.3] They are dropped onto the grain bed and relaxed into a local potential energy minimum or, alternatively, into the global minimum.

[220.4.1] Compaction occurs because the sand becomes buried into the subsurface. [220.4.2] Compaction reduces the bulk volume (and porosity). [220.4.3] It is modelled as a linear process in which the vertical coordinate of every sandgrain is shifted vertically downwards by an amount proportional to the original vertical position. [220.4.4] The proportionality constant is called the compaction factor. [220.4.5] Its value for the Fontainebleau sample is estimated to be $0.1$ from thin section analysis.

[220.5.1] In the diagenesis part only a subset of known diagenetical processes are simulated, namely quartz cement overgrowth and precipitation of authigenic clay on the free surface. [220.5.2] Quartz cement overgrowth is modeled by radially enlarging each grain. [220.5.3] If $R_{0}$ denotes the radius of the originally deposited spherical grain, its new radius along the direction ${\boldsymbol{r}}$ from grain center is taken to be [59, 48]

$R({\boldsymbol{r}})=R_{0}+\min(b\ell({\boldsymbol{r}})^{\gamma},\ell({\boldsymbol{r}}))$

(59)

where $\ell({\boldsymbol{r}})$ is the distance between the surface of the original spherical grain and the surface of its Voronoi polyhedron along the direction ${\boldsymbol{r}}$ . [220.5.4] The constant $b$ controls the amount of cement, and the growth exponent $\gamma$ controls the type of cement overgrowth. [220.5.5] For $\gamma>0$ the cement grows preferentially into the pore bodies, for $\gamma=0$ it grows concentrically, and for $\gamma<0$ quartz cement grows towards the pore throats [48]. [220.5.6] Authigenic clay growth is simulated by precipitating clay voxels on the free mineral surface. The clay texture may be pore-lining or pore-filling or a combination of the two.

[220.6.1] The parameters for modeling the Fontainebleau sandstone were $0.1$ for the compaction factor, and $\gamma=-0.6$ and $b=2.9157$ for the cementation parameters. [220.6.2] The resulting model configuration of the sample DM is displayed in Figure 2.

[page 221, §0]

Figure 2: Sample DM: Threedimensional pore space $\mathbb{P}_{{{\sf DM}}}$ of the sedimentation and diagenesis model described in the text. The resolution is $a=7.5\ \mu\text{m}$ , the sample dimensions are $M_{1}=255$ , $M_{2}=255$ , $M_{3}=255$ . The pore space is indicated opaque, the matrix space is transparent. The right image shows the front plane of the sample as a twodimensional thin section (pore space black, matrix grey)

[page 222, §1]

4.3 Gaussian Field Reconstruction Model

[222.1.1] A stochastic reconstruction model attempts to approximate a given experimental sample by a randomly generated model structure that matches prescribed stochastic properties of the experimental sample. [222.1.2] In this and the next section the stochastic property of interest is the correlation function $G_{{\sf EX}}({\boldsymbol{r}})$ of the Fontainebleau sandstone.

[222.2.1] The Gaussian field (GF) reconstruction model tries to match a reference correlation function by filtering Gaussian random variables [49, 2, 1, 69]. [222.2.2] Given the reference correlation function $G_{{\sf EX}}({\boldsymbol{r}})$ and porosity $\phi(\mathbb{S}_{{\sf EX}})$ of the experimental sample the Gaussian field method proceeds in three main steps:

Initially a Gaussian field $X({\boldsymbol{r}})$ is generated consisting of statistically independent Gaussian random variables $X\in\mathbb{R}$ at each lattice point ${\boldsymbol{r}}$ .
[222.2.3] The field $X({\boldsymbol{r}})$ is first passed through a linear filter which produces a correlated Gaussian field $Y({\boldsymbol{r}})$ with zero mean and unit variance. [222.2.4] The reference correlation function $G_{{\sf EX}}({\boldsymbol{r}})$ and porosity $\phi(\mathbb{S}_{{\sf EX}})$ enter into the mathematical construction of this linear filter.
[222.2.5] The correlated field $Y({\boldsymbol{r}})$ is then passed through a nonlinear discretization filter which produces the reconstructed sample $\mathbb{S}_{{\sf GF}}$ .

[222.2.6] Step 4.3 is costly because it requires the solution of a very large set of non-linear equations. [222.2.7] A computationally more efficient method uses Fourier Transformation [1]. [222.2.8] The linear filter in step 4.3 is defined in Fourier space through

$Y(\boldsymbol{k})=\alpha(G_{Y}(\boldsymbol{k}))^{{\frac{1}{2}}}X(\boldsymbol{k}),$

(60)

where $M=M_{1}=M_{2}=M_{3}$ is the sidelength of a cubic sample, $\alpha=M^{{\frac{d}{2}}}$ is a normalisation factor, and

$X(\boldsymbol{k})=\frac{1}{M^{d}}\sum _{{{\boldsymbol{r}}}}X({\boldsymbol{r}})e^{{2\pi i\boldsymbol{k}\cdot{\boldsymbol{r}}}}$

(61)

denotes the Fourier transform of $X({\boldsymbol{r}})$ . [222.2.9] Similarly $Y(\boldsymbol{k})$ is the Fourier transform of $Y({\boldsymbol{r}})$ , and $G_{Y}(\boldsymbol{k})$ is the Fourier transform of the correlation function $G_{Y}({\boldsymbol{r}})$ . [222.2.10] $G_{Y}({\boldsymbol{r}})$ has to be computed by an inverse process from the correlation function $G_{{\sf EX}}({\boldsymbol{r}})$ and porosity of the experimental reference (details in [1]).

[222.3.1] The Gaussian field reconstruction requires a large separation $\xi _{{\sf EX}}\ll N^{{1/d}}$ where $\xi _{{\sf EX}}$ is the correlation length of the experimental reference, and $N=M_{1}M_{2}M_{3}$ is the number of sites. [222.3.2] $\xi _{{\sf EX}}$ is defined as the length such that $G_{{\sf EX}}(r)\approx 0$ for $r>\xi _{{\sf EX}}$ . [222.3.3] If the condition $\xi _{{\sf EX}}\ll N^{{1/d}}$ is violated then step 4.3 of the reconstruction fails in the sense that the correlated Gaussian field $Y({\boldsymbol{r}})$ does not have zero mean and unit variance. [222.3.4] In such a situation the filter $G_{Y}(\boldsymbol{k})$ will differ from the Fourier transform of the correlation function of the $Y({\boldsymbol{r}})$ . [222.3.5] It is also difficult to calculate $G_{Y}(r)$ accurately near $r=0$ [1]. [222.3.6] This leads to a discrepancy at small $r$ between $G_{{\sf GF}}(r)$ and $G_{{\sf EX}}(r)$ . [222.3.7] The problem can be overcome by choosing large $M$ . [222.3.8] However, in $d=3$ very large $M$ also demands prohibitively large memory. [222.3.9] In

[page 223, §0]

Figure 3: Sample GF: Threedimensional pore space $\mathbb{P}_{{{\sf GF}}}$ with $G_{{{\sf GF}}}(r)\approx G_{{{\sf EX}}}(r)$ constructed by filtering Gaussian random fields. The resolution is $a=7.5\ \mu\text{m}$ , the sample dimensions are $M_{1}=256$ , $M_{2}=256$ , $M_{3}=256$ . The pore space is indicated opaque, the matrix space is transparent. The right image shows the front plane of the sample as a twodimensional thin section (pore space black, matrix grey)

[page 224, §0]

earlier work [2, 1] the correlation function $G_{{\sf EX}}({\boldsymbol{r}})$ was sampled down to a lower resolution, and the reconstruction algorithm then proceeded with such a rescaled correlation function. [224.0.1] This leads to a reconstructed sample $\mathbb{S}_{{\sf GF}}$ which also has a lower resolution. [224.0.2] Such reconstructions have lower average connectivity compared to the original model [9] For a quantitative comparison with the microstructure of $\mathbb{S}_{{\sf EX}}$ it is necessary to retain the original level of resolution, and to use the original correlation function $G_{{\sf EX}}({\boldsymbol{r}})$ without subsampling. [224.0.3] Because $G_{{\sf EX}}(r)$ is nearly $0$ for $r>30a$ $G_{{\sf EX}}(r)$ was truncated at $r=30a$ to save computer time. [224.0.4] The resulting configuration $\mathbb{S}_{{\sf GF}}$ with $M=256$ is displayed in Figure 3.

4.4 Simulated Annealing Reconstruction Model

[224.0.5] The simulated annealing (SA) reconstruction model is a second method to generate a threedimensional random microstructure with prescribed porosity and correlation function. [224.0.6] The method generates a configuration $\mathbb{S}_{{\sf SA}}$ by minimizing the deviations between $G_{{\sf SA}}({\boldsymbol{r}})$ and a predefined reference function $G_{0}({\boldsymbol{r}})$ . [224.0.7] Note that the generated configuration $\mathbb{S}_{{\sf SA}}$ is not unique and hence other modeling aspects come into play [42]. [224.0.8] Below, $G_{0}({\boldsymbol{r}})=G_{{\sf EX}}({\boldsymbol{r}})$ is again the correlation function of the Fontainebleau sandstone.

[224.1.1] An advantage of the simulated annealing method over the Gaussian field method is that it can also be used to match other quantities besides the correlation function. [224.1.2] Examples would be the linear or spherical contact distributions [42]. [224.1.3] On the other hand the method is computationally very demanding, and cannot be implemented fully at present. [224.1.4] A simplified implementation was discussed in [70], and is used below.

[224.2.1] The reconstruction is performed on a cubic lattice with side length $M=M_{1}=M_{2}=M_{3}$ and lattice spacing $a$ . [224.2.2] The lattice is initialized randomly with $0$ ’s and $1$ ’s such that the volume fraction of $0$ ’s equals $\phi(\mathbb{S}_{{\sf EX}})$ . [224.2.3] This porosity is preserved throughout the simulation. [224.2.4] For the sake of numerical efficiency the autocorrelation function is evaluated in a simplified form using [70]

		$\displaystyle{\widetilde{G}}_{{\sf SA}}(r)\left({\widetilde{G}}_{{\sf SA}}(0)-{\widetilde{G}}_{{\sf SA}}(0)^{2}\right)+{\widetilde{G}}_{{\sf SA}}(0)^{2}=$
		$\displaystyle=\frac{1}{3M^{3}}\sum _{{\mathbf{r}}}\chi\rule[-4.3pt]{0.0pt}{8.6pt}_{{\mathbb{M}}}({\mathbf{r}})\left(\chi\rule[-4.3pt]{0.0pt}{8.6pt}_{{\mathbb{M}}}({\mathbf{r}}+r{\mathbf{e}}_{1})+\chi\rule[-4.3pt]{0.0pt}{8.6pt}_{{\mathbb{M}}}({\mathbf{r}}+r{\mathbf{e}}_{2})+\chi\rule[-4.3pt]{0.0pt}{8.6pt}_{{\mathbb{M}}}({\mathbf{r}}+r{\mathbf{e}}_{3})\right)$		(62)

where ${\mathbf{e}}_{i}$ are the unit vectors in direction of the coordinate axes, $r=0,\dots,\frac{M}{2}-1$ , and where a tilde $\widetilde{~}$ is used to indicate the directional restriction. [224.2.5] The sum $\sum _{{\mathbf{r}}}$ runs over all $M^{3}$ lattice sites ${\boldsymbol{r}}$ with periodic boundary conditions, i.e. $r_{i}+r$ is evaluated modulo $M$ .

[224.3.1] A simulated annealing algorithm is used to minimize the "energy" function

$E=\sum _{r}\left({\widetilde{G}}_{{\sf SA}}(r)-G_{{\sf EX}}(r)\right)^{2},$

(63)

defined as the sum of the squared deviations of ${\widetilde{G}}_{{\sf SA}}$ from the experimental correlation function $G_{{\sf EX}}$ . [224.3.2] Each update starts with the exchange of two pixels, one

[page 225, §0]

Figure 4: Sample SA: Threedimensional pore space $\mathbb{P}_{{{\sf SA}}}$ with $G_{{{\sf SA}}}(r)=G_{{{\sf EX}}}(r)$ constructed using a simulated annealing algorithm. The resolution is $a=7.5\ \mu\text{m}$ , the sample dimensions are $M_{1}=256$ , $M_{2}=256$ , $M_{3}=256$ . The pore space is indicated opaque, the matrix space is transparent. The right image shows the front plane of the sample as a twodimensional thin section (pore space black, matrix grey)

[page 226, §0] from pore space, one from matrix space. [226.0.1] Let $n$ denote the number of the proposed update step. [226.0.2] Introducing an acceptance parameter $T_{n}$ , which may be interpreted as an $n$ -dependent temperature, the proposed configuration is accepted with probability

$p=\min\left(1,\exp\left(-\frac{E_{n}-E_{{n-1}}}{T_{n}E_{{n-1}}}\right)\right).$

(64)

[226.0.3] Here the energy and the correlation function of the configuration is denoted as $E_{n}$ and ${\widetilde{G}}_{{{\sf SA},n}}$ , respectively. [226.0.4] If the proposed move is rejected, then the old configuration is restored.

[226.1.1] A configuration with correlation $G_{{\sf EX}}$ is found by lowering $T$ . [226.1.2] At low $T$ the system approaches a configuration that minimizes the energy function. [226.1.3] In the simulations $T_{n}$ was lowered with $n$ as

$T_{n}=\exp\left(-\frac{n}{100000}\right).$

(65)

[226.1.4] The simulation was stopped when 20000 consecutive updates were rejected. [226.1.5] This happened after $2.5\times 10^{8}$ updates ( $\approx 15$ steps per site). [226.1.6] The resulting configuration $\mathbb{S}_{{\sf SA}}$ for the simulated annealing reconstruction is displayed in Figure 4.

[226.2.1] A complete evaluation of the correlation function as defined in (29) for a threedimensional system requires so much computer time, that it cannot be carried out at present. [226.2.2] Therefore the algorithm was simplified to increase its speed [70]. [226.2.3] In the simplified algorithm the correlation function is only evaluated along the directions of the coordinate axes as indicated in (62). [226.2.4] The original motivation was that for isotropic systems all directions should be equivalent [70]. [226.2.5] However, it was found in [41] that as a result of this simplification the reconstructed sample may become anisotropic. [226.2.6] In the simplified algorithm the correlation function of the reconstruction deviates from the reference correlation function in all directions other than those of the axes [41]. [226.2.7] The problem is illustrated in Figures 5(a) and 5(b) in two dimensions for a reference correlation function given as

$G_{0}(r)=e^{{-r/8}}\cos r.$

(66)

[226.2.8] In Figure 5(a) the correlation function was matched only in the direction of the $x$ - and $y$ -axis. [226.2.9] In Figure 5(b) the correlation function was matched also along the diagonal directions obtained by rotating the axes 45 degrees. [226.2.10] The differences in isotropy of the two reconstructions are clearly visible. [226.2.11] In the special case of the correlation function of the Fontainebleau sandstone, however, this effect seems to be smaller. [226.2.12] The Fontainebleau correlation function is given in Figure 7 below. [226.2.13] Figure 6(a) and 6(b) show the result of twodimensional reconstructions along the axes only and along axes plus diagonal directions. [226.2.14] The differences in isotropy seem to be less pronounced. [226.2.15] Perhaps this is due to the fact that the Fontainebleau correlation function has no maxima and minima.

[page 227, §0]

Figure 5: Twodimensional stochastic reconstruction for the correlation function of $G_{0}(r)=e^{{-r/8}}\cos r$ (a) for the direction of the $x$ - and $y$ -coordinate axes only, and (b) for the directions of the coordinate axes plus diagonal directions.

Figure 6: A Twodimensional stochastic reconstruction for the correlation function $G_{0}(r)=G_{{{\sf EX}}}(r)$ displayed as the solid line in Figure 7 (a) along the direction of the $x$ - and $y$ -coordinate axes only, and (b) along the directions of the coordinate axes plus diagonal directions.

Author:	R. Hilfer
originally published in:	Räumliche Statistik und Statistische Physik, Lecture Notes in Physics Vol. 554, Springer, Berlin, page 203-241 (2000), ISBN - -
originally published:	2000