# On the Reconstruction of Random Media using Monte Carlo Methods

## I Introduction

A better understanding of the transport properties of random media, such as fluid flow in sandstones or electrical conductivity of composites requires the micro-structure as input . Digitized, three-dimensional micro-structures of natural sandstones are difficult and expensive to obtain . Thus there is a need for simulation algorithms that are able to provide representative micro-structures from statistical probability functions.

Recently various algorithms have been proposed for the reconstruction of random micro-structures . In this paper we investigate a simulated annealing method that enforces agreement between the correlation functions of the original and the reconstructed micro-structure . To save computation time the authors of  have evaluated the correlation functions only in certain directions assuming isotropy of the medium. The objective of this paper is to study the effect of this simplification on the final, reconstructed configurations. We test the effects on two examples: (i) the correlation function of a Fontainebleau sandstone and (ii) an artificial correlation function with damped oscillations. The results show, that at least for the oscillating correlation function this yields significantly different configurations. More importantly, the reconstructions are anisotropic as a result of the simplified evaluation of the correlation functions.

## II The reconstruction method

We follow the reconstruction algorithm proposed in . The reconstruction is performed on a $`d`$-dimensional hyper-cubic lattice $`^d`$ ($`d=2`$ for the results presented below). Whether a lattice point lies within pore space or matrix space is indicated by the characteristic function

$$\chi \text{}(\stackrel{}{x})=\chi \text{}(x_1,x_2,\mathrm{},x_d)=\{\begin{array}{c}0\text{for}\stackrel{}{x}\hfill \\ 1\text{for}\stackrel{}{x}𝕄\hfill \end{array}$$
(1)   
with $`x_i=0,1,\mathrm{},M_i1`$ where $``$ denotes the pore space and $`𝕄`$ the matrix space of a two phase porous medium. The $`x_i`$ are in units of the lattice spacing $`a`$. The porosity $`\varphi `$ is given as $`\varphi =\frac{1}{N}_{j=1}^N\left(1\chi \text{}(\stackrel{}{x}_j)\right)`$ where $`N=_{i=1}^dM_i`$ is the total number of lattice sites.

Simulated annealing is an iterative technique for combinatorial optimization problems. The iteration steps are denoted by a subscript $`t`$. The optimum is found by lowering a fictitious temperature $`T_t`$ that controls the acceptance or rejection of configurations with ”energy” (or cost function) $`E_t`$. The energy function used in our simulations is defined as

$$E_t=\underset{k=1}{\overset{J}{}}w_k\underset{\stackrel{}{r}𝔻_k}{}\left(g_t^k(\stackrel{}{r})g_{\mathrm{ref}}^k(\stackrel{}{r})\right)^2$$
(2)   
where $`𝔻_k^d`$ is a subset of lattice points and $`g_t^k`$ is the $`k`$th function of a set of $`J`$ statistical probability functions calculated for the configuration of step $`t`$. For example $`g_t^k`$ may be a $`k`$-point correlation function. The real valued factor $`w_k`$ is a weight for the $`k`$th function. Hence, the energy can be understood as a measure for the deviations of the probability functions $`g_t^k`$ from predefined reference functions $`g_{\mathrm{ref}}^k`$.

The simulated annealing algorithm consists of the following steps.

1. Initialization: The 0’s and 1’s are randomly distributed with given porosity $`\varphi `$.
2. Two lattice points of different phase are chosen at random and exchanged. In this way the porosity $`\varphi `$ is conserved.
3. The ”energy” $`E_t`$ of the current configuration is calculated according to Equation (2).
4. The ”temperature” $`T_t`$ is adjusted according to a fixed cooling schedule.
5. The new configuration created by the exchange of the two points is accepted with probability

$$p=\mathrm{min}(1,\mathrm{exp}\left(\frac{E_tE_{t1}}{T_t}\right)).$$
(3)   
In case of rejection, the two points are restored and the old configuration is left unchanged.
6. Return to step 2.

As can be seen from Equation (3), configurations with lower energy are immediately accepted while the acceptance of a configuration with higher energy is controlled by the temperature $`T`$. In order to obtain a configuration with minimal energy $`E`$ or, in other words, a configuration with minimal deviations of the probability functions $`g^k`$ from their reference functions $`g_{\mathrm{ref}}^k`$, the temperature $`T`$ has to be lowered in a suitable way.

The algorithm is applicable to various functions $`g^k`$. Most authors propose the two point correlation function . Assuming homogeneity and ergodicity, the two-point correlation function can be defined for our case as

$$g(\stackrel{}{r})=\frac{\chi \text{}(\stackrel{}{x})\chi \text{}(\stackrel{}{x}+\stackrel{}{r})(1\varphi )^2}{\varphi \varphi ^2}.$$
(4)   
where the average $`\mathrm{}`$ indicates a spatial average over all lattice sites $`\stackrel{}{x}`$. If $`g_{t1}(\stackrel{}{r})`$ is known a single update step of the annealing process requires the recalculation of the correlations of the two pixels that are exchanged in step 2 with all other pixels. The numerical effort to obtain $`g_t(\stackrel{}{r})`$ is therefore proportional to $`N`$ where $`N`$ is the total number of lattice points. In the reconstruction of three-dimensional porous media with $`N10^6\mathrm{}10^7`$ this leads to unacceptably long calculation times and therefore one has to find ways to speed up this calculation.

One possibility for reducing the numerical effort is to truncate $`g`$ at a value $`r_c`$ for which $`g(\stackrel{}{r})0`$ for $`r=|\stackrel{}{r}|r_c`$ holds. Below we set $`g(\stackrel{}{r})=0`$ for $`|\stackrel{}{r}|r_c`$ where $`r_c`$ is a parameter in the reconstruction. For isotropic media, where $`g(\stackrel{}{r})=g(r)`$ with $`|\stackrel{}{r}|=r`$, it was suggested in  to calculate $`g(r)`$ only in certain directions by setting $`\stackrel{}{r}=r\stackrel{}{e}_k`$ where $`\stackrel{}{e}_k`$ is an arbitrary unit vector. In the two-dimensional reconstructions presented below $`\stackrel{}{e}_k`$ will be set to the radial unit vector in a polar coordinate system $`\stackrel{}{e}_k=\stackrel{}{e}_{\phi _k}=\stackrel{}{e}_x\mathrm{cos}\phi _k+\stackrel{}{e}_y\mathrm{sin}\phi _k`$ where $`\stackrel{}{e}_x`$ and $`\stackrel{}{e}_y`$ are the unit vectors of the Cartesian coordinate system and $`\phi _k`$ is the angle between $`\stackrel{}{e}_x`$ and $`\stackrel{}{e}_{\phi _k}`$. Hence, instead of Equation (4) we use

$$g^k(r)=\frac{\chi \text{}(\stackrel{}{x})\chi \text{}(\stackrel{}{x}+r\stackrel{}{e}_k)(1\varphi )^2}{\varphi \varphi ^2}$$
(5)   
with $`r=0,1,\mathrm{},r_c`$ in the simplified reconstruction scheme. Since (5) is a set of $`J`$ one-dimensional correlation functions the numerical effort is reduced by a factor of roughly $`Jr_c/N`$ as compared to (4).

The above algorithm is now used to reconstruct two-dimensional, isotropic media with correlation function

$$g_{\mathrm{ref}}(r)=\mathrm{exp}\left(\frac{r}{8}\right)\mathrm{cos}(\omega r)$$
(6)   
with $`r`$ in units of the lattice spacing and $`\omega =1`$. The same function was used by the authors of  to exemplify their algorithm. In the evaluation of the correlation functions periodic boundary conditions are assumed. We use an exponential decrease of the temperature

$$T_t=\mathrm{exp}\left(\frac{t}{1610^5}\right).$$
(7)   
The remaining parameters are the lattice size $`M_1=M_2=400`$ and the porosity $`\varphi =0.5`$. They are the same as in . The following reconstructions have all been initialized with the same random seed. The algorithm terminated when the configuration did not change for 20000 subsequent update steps.

## III Results

In the first reconstruction, the correlation function was calculated only in the horizontal and vertical direction, i.e. Equations (2) and (5) were used with $`J=2`$, $`\stackrel{}{r}_1=r\stackrel{}{e}_0`$ and $`\stackrel{}{r}_2=r\stackrel{}{e}_{\pi /2}`$. Both correlation functions $`g^1,g^2`$ had the same weight $`w_1=w_2=0.5`$ and the same reference function $`g_{\mathrm{ref}}^1=g_{\mathrm{ref}}^2`$ given in Equation (6) was used. The correlation function was truncated at $`r_c=100`$. Figure 1 shows the final configuration of the reconstruction. A similar pattern was found in . The pattern consists of stripes in direction of $`\stackrel{}{e}_{\pi /4}`$ and in direction $`\stackrel{}{e}_{\pi /4}`$. The distance between the stripes is determined by the cosine term. On a larger length scale, the pattern organizes into several regions in which all lines are parallel. The typical size of these regions is given by roughly 20 in view of Equation (6) and Figure 2. Clearly, the pattern is not isotropic as it should be. The stripes are preferentially directed along the directions $`\stackrel{}{e}_{\pi /4}`$ and $`\stackrel{}{e}_{\pi /4}`$. Also, one expects that the oscillation frequency $`\omega `$ of the correlation function in direction $`\stackrel{}{e}_{\pm \pi /4}`$ is not $`\omega =1`$ but $`\omega =\sqrt{2}`$.

Figure 2 shows the correlation functions for the configuration of Figure 1 in the directions of the unit vectors $`\stackrel{}{e}_0,\stackrel{}{e}_{\pi /2}`$ and $`\stackrel{}{e}_{\pi /4}`$. The first and second have been used for the reconstruction and hence show good agreement with the reference function while the latter deviates drastically. The correlation in direction $`\stackrel{}{e}_{\pi /4}`$ is better described by the function

$$f(r)=\frac{1}{2}\left(\mathrm{exp}\left(\frac{r}{8}\right)\mathrm{cos}\left(\sqrt{2}r\right)+\mathrm{exp}\left(\frac{r}{8}\right)\right).$$
(8)   
shown as the dotted line. This may be interpreted as the arithmetic mean of the correlation function for regions (described above) with stripes perpendicular to the direction $`\stackrel{}{e}_{\pi /4}`$ given by the first term in Equation (8) and the correlation function for regions with stripes parallel to the direction $`\stackrel{}{e}_{\pi /4}`$ given by the second term. Of course, the same correlation function is found in direction $`\stackrel{}{e}_{\pi /4}`$.

One step towards an isotropic reconstruction may be to use $`J=4`$, and to force agreement of the correlations not only in two but in four directions $`\stackrel{}{e}_0`$, $`\stackrel{}{e}_{\pi /2}`$, $`\stackrel{}{e}_{\pi /4}`$ and $`\stackrel{}{e}_{\pi /4}`$. The resulting configuration is shown in Figure 3. The pattern is significantly different from the pattern of Figure 1. There are not only stripes parallel to the diagonal directions but more rounded formations and concentric circles at some points. The correlation functions for Figure 3 are plotted in Figure 4. Surprisingly we find that it is not possible to obtain agreement with the reference correlation function. The resulting correlation functions in the horizontal, vertical and diagonal direction all deviate strongly from the reference function especially at the first minimum. Additional simulations with different cooling schedules, damping factors and frequency of the correlation function did not give better agreement. We have also varied the system sizes from $`200\times 200`$ up to $`1000\times 1000`$. The results were identical.

Figure 5 shows the two-dimensional reconstructions of a Fontainebleau sandstone with porosity $`\varphi =0.135`$. The correlation function was truncated at $`r_c=50`$. The reconstruction of Figure 5a used the correlation function only in vertical $`\stackrel{}{e}_0`$ and horizontal $`\stackrel{}{e}_{\pi /2}`$ direction as suggested in . The reconstruction shown in Figure 5b, however, calculates the full two-dimensional correlation function according to (4) where the two-dimensional correlation function was radially binned in the calculations to obtain a one-dimensional function which can be compared to the one-dimensional, isotropic correlation function of the original sandstone. We emphasize that in this calculation the two-dimensional correlation function was evaluated without restrictions or simplifications. Therefore, the calculation needed about 50 times longer than the simplified reconstruction. Again, there are differences visible although they are smaller than those in the reconstructions using (6). The shape of the pores in the full reconstruction (Figure 5b) appear to be smoother and there is not as much ”dust” visible as in the simplified (restricted) reconstruction shown in Figure 5a. The correlation functions plotted in Figure 6 reveal also that the reconstructed micro-structure in Figure 5a is strongly anisotropic while the micro-structure in Figure 5b is isotropic as it should be.

In summary we note that the statistical reconstruction of porous media with predefined two-point correlation function often requires some reduction of numerical effort. Especially the reconstruction of three-dimensional porous media in reasonable time does not seem possible without such simplifications in the calculation of the correlation function. Of course this problem is exacerbated if one wishes to include three-point or even higher order correlation functions. We applied a simplification proposed and used in  which samples the correlation function only in certain directions. The effect of this simplification on the final configurations may in some cases be negligible but in general configurations strong anisotropy and patterns which are significantly different from those of a proper isotropic reconstruction may appear as a result.

## Acknowledgments

We thank B. Biswal for comments and helpful discussions, K. Höfler, S. Schwarzer and M. Müller for technical advise and significant parts of the C++ code, the Deutsche Forschungsgemeinschaft for financial support through the GKKS Stuttgart. 

## Figure captions

* Simplified reconstruction of the damped oscillating correlation function given in Equation (6) by restricting the correlation function evaluation to the horizontal and vertical directions.
* Correlation functions for the reconstruction of Figure 1. The solid line is the reference function of the reconstruction given in Equation (6). $`+`$ and $`\times `$ are the values of the correlation function of Figure 1 along the horizontal and vertical direction, respectively. The values of the correlation function in direction $`\stackrel{}{e}_{\pi /4}`$ are plotted with $``$. The dotted line is the estimate given in Equation (8) for the correlation function along the $`\stackrel{}{e}_{\pi /4}`$-direction.
* Simplified reconstruction of the damped oscillating correlation function of Equation (6) by restricting the correlation function evaluation to the four $`\stackrel{}{e}_0`$, $`\stackrel{}{e}_{\pi /2}`$, $`\stackrel{}{e}_{\pi /4}`$ and $`\stackrel{}{e}_{\pi /4}`$ directions.
* Correlation functions for the reconstruction of Figure 3. The reference function is plotted as solid line. Note the mismatch in the first minimum and maximum.
* (a) Two-dimensional simplified reconstruction of the correlation function of a Fontainebleau sandstone by restricting the correlation function evaluation to the horizontal and vertical directions. (b) Reconstruction with the same reference correlation function as in (a) but complete evaluation of the correlation function, i.e. without directional restrictions.
* Correlation function for the reconstructions of Figure 5. The reference function is plotted as solid line. The correlation functions calculated in horizontal, vertical and diagonal direction refer to the configuration of Figure 5a. The complete correlation function is calculated from the configuration of Figure 5b.

C.Manwart and R.Hilfer Figure 1   

C.Manwart and R.Hilfer Figure 2   

C.Manwart and R.Hilfer Figure 3   

C.Manwart and R.Hilfer Figure 4   

C.Manwart and R.Hilfer Figure 5   
(a)

(b)

C.Manwart and R.Hilfer Figure 6