text:
# 1 Introduction

## 1 Introduction

It has been proposed that there may exist extra dimensions which are almost the size of a millimeter . In this case, the gravitational potential would go as $`1/r^{1+n}`$ for sub-millimeter measurements, where $`n`$ is the number of extra millimeter dimensions. The four-dimensional potential is recovered at larger scales. If the four-dimensional metric is multiplied by a ”warp” factor which is a rapidly changing function of an extra dimension, then four-dimensional gravity is recovered at low-energy scales even for an infinite extra dimension . This is due to a massless graviton state that is bound to a 3-brane embedded in the higher-dimensional space. Large-energy modifications to the four-dimensional potential would result from the Kaluza-Klein modes propagating in the extra dimension.

It is possible that nature has chosen some extra dimensions to be infinite and others to be compact, especially if our spacetime is embedded within M-theory. For example, five-dimensional domain walls which localize gravity may arise from a sphere reduction from ten or eleven dimensions, as the near-horizon of extremal $`p`$-brane configurations 111It has been found that a non-extremal brane origin yields a localized graviton of nonzero mass .. In particular, it was found that such domain walls can be derived from a D3, D4 (or M5) and D5-brane, as well as two M5-branes intersecting over a 3-brane222Other $`p`$-brane origins are possible with the inclusion of a naked singularity ..

In embedding five-dimensional brane world scenarios in higher dimensions, one generally only considers the s-wave Kaluza-Klein graviton mode. That is, as part of the dimensional reduction ansatz, one discards modes moving about the compact directions. It would appear to be a reasonable assumption that, unless the compact dimensions are quite large, the effects of such rotation modes on the corrections to the effective four-dimensional gravitation would be over-shadowed by the s-wave Kaluza-Klein modes propagating along the infinite extra dimension. However, it once appeared to be reasonable to assume that a theory with more than 1+3 large dimensions did not agree with our direct observations.

It was found that the ”warp” factor in certain five-dimensional spaces is such that the effects of gravitational modes propagating in the extra dimension is dampened, yielding four-dimensional gravity at low energies. Embedding the brane world in ten or eleven dimensions brings the possibility of additional ”warp” factors. In particular, the effects of modes rotating about compact dimensions could be amplified if the corresponding metric elements are multiplied by a ”warp” factor which increases rapidly along the infinite extra dimension. Depending on the geometry of the higher-dimensional origin, rotational Kaluza-Klein modes may actually be as significant as the s-wave Kaluza-Klein modes333Throughout this paper, we assume that all compact dimensions are the same size, fixed by the length scale 1/$`k`$ in the ”warp” factors..

This paper is organized as follows. In section $`p1`$ we analyse how the rotational graviton modes modify the four-dimensional gravitational potential for a brane world originating from a D$`p`$-brane, where $`p=3,4,5`$. In these cases, we find that the rotational modes do not contribute significantly to the modification of the gravitational potential relative to the corrections due to the s-wave Kaluza-Klein modes. It is interesting to note that, in some cases, there are rotational modes that are localized on the brane world. In section 5, we are pleasantly surprised to find that the metric for the near-horizon of the M5/M5 system contains a ”warp” factor which amplifies the effect of the rotational graviton modes. In this case, the corrective contributions from the s-wave Kaluza-Klein graviton modes and the rotational Kaluza-Klein graviton modes are of the same magnitude. We offer Conclusions in section 6.

## 2 D3-brane origin

The D3-brane metric is

$$ds_{10}^2=H^{1/2}(dt^2+dx_i^2)+H^{1/2}(dr^2+r^2d\mathrm{\Omega }_5^2),$$
(2.1)
where

$$H=1+\frac{R^4}{r^4},$$
(2.2)
and $`i=1,2,3`$. In the near-horizon limit $`r<<R`$, the above metric is that of $`AdS_5\times S^5`$. This can be expressed as

$$ds_{10}^2=(1+k|z|)^2(dt^2+dx_i^2+dz^2)+\frac{1}{k^2}d\mathrm{\Omega }_5^2,$$
(2.3)
where

$$\frac{r}{R}=(1+k|z|)^1,$$
(2.4)
and $`R=k^1`$. Note that the absolute value in the above coordinate transformations imposes $`Z_2`$ symmetry, which has no known source in ten dimensions444Recently, it has been found that resolved branes on an Eguchi-Hanson instanton dimensionally reduce to five-dimensional domain walls which localize gravity, without the need of an extra source term..

The equation of motion for a graviton fluctuation is that of a minimally-coupled scalar:

$$_M\sqrt{g}g^{MN}_N\mathrm{\Phi }=0.$$
(2.5)
We take $`\mathrm{\Phi }=P_{\mathrm{}}(S^5)(1+k|z|)^{3/2}\psi (z)M(t,x_i)`$, where $`\text{    }\text{ }_{(4)}M=m^2M`$ and $`\text{    }\text{ }_{(4)}`$ is the Laplacian on $`t,x_i`$. $`P_{\mathrm{}}`$ are the Legendre polynomials of $`S^5`$. For graviton modes on the metric (2.3) we obtain

$$_z^2\psi +\left[\frac{(\mathrm{}+3/2)(\mathrm{}+5/2)k^2}{(1+k|z|)^2}3k\delta (z)\right]\psi =m^2\psi ,$$
(2.6)
where $`\mathrm{}`$ is an integer parametrizing the rotational dynamics about $`S^5`$. For the massless modes, the solution is

$$\psi =(1+k|z|)^{\mathrm{}3/2}+\frac{\mathrm{}}{\mathrm{}+4}(1+k|z|)^{\mathrm{}+5/2}.$$
(2.7)
For $`\mathrm{}=0`$, this wave function is normalizable and there is a corresponding localized graviton. However, for $`\mathrm{}>0`$ the wave function diverges for large $`z`$, and thus there is no permissible massless solution. The massive wave functions are given by

$$\psi _m=N_m(1+k|z|)^{1/2}[(Y_{1+\mathrm{}}\left(\frac{m}{k}\right)\mathrm{}\frac{k}{m}Y_{2+\mathrm{}}\left(\frac{m}{k}\right))J_{2+\mathrm{}}\left(\frac{m}{k}(1+k|z|)\right)+$$

$$(\mathrm{}\frac{k}{m}J_{2+\mathrm{}}\left(\frac{m}{k}\right)J_{1+\mathrm{}}\left(\frac{m}{k}\right))Y_{2+\mathrm{}}\left(\frac{m}{k}(1+k|z|)\right)],$$
(2.8)
where

$$N_m=\frac{1}{\sqrt{2k}}\left[\left(Y_{1+\mathrm{}}\left(\frac{m}{k}\right)\mathrm{}\frac{k}{m}Y_{2+\mathrm{}}\left(\frac{m}{k}\right)\right)^2+\left(\mathrm{}\frac{k}{m}J_{2+\mathrm{}}\left(\frac{m}{k}\right)J_{1+\mathrm{}}\left(\frac{m}{k}\right)\right)^2\right]^{1/2}.$$
(2.9)
The Newtonian gravitational potential between masses $`M_1`$ and $`M_2`$ can be estimated by

$$U(r)\underset{\mathrm{}=0}{\overset{\mathrm{}}{}}\left(\frac{G_4M_1M_2}{r}\mathrm{e}^{m_{\mathrm{}}r}+\frac{G_5M_1M_2}{r}_{m_0^2}^{\mathrm{}}d(m^2)|\psi _m(z=0,\mathrm{})|^2\mathrm{e}^{mr}\right),$$
(2.10)
where $`m_{\mathrm{}}`$ is the mass of a bound state and $`m_0`$ is the minimum mass for non-localized states, each for a given $`\mathrm{}`$. In the present case, there are no massive bound states or mass gaps due to $`\mathrm{}>0`$ states, so that $`m_0=m_{\mathrm{}}=0`$. The four and five-dimensional Newton constants are related by $`G_4=kG_5`$. This yields

$$U(r)\frac{G_4M_1M_2}{r}\left(\left[1+\frac{c_0}{(kr)^2}+\mathrm{}\right]+\underset{\mathrm{}=1}{\overset{\mathrm{}}{}}\left[\frac{c_{\mathrm{}}}{(kr)^{6+2\mathrm{}}}+\mathrm{}\right]\right),$$
(2.11)
where $`c_0`$ and $`c_{\mathrm{}}`$ are constants of order one. The $`\mathrm{}=0`$ terms were found in . The terms in the second pair of square brackets are the dominant corrections from the rotational graviton modes ($`\mathrm{}>0`$). As can be seen, the dominant contribution from the $`\mathrm{}>0`$ modes is smaller than that of the s-wave Kaluza-Klein modes by the rather small factor $`1/(kr)^6`$.

## 3 D4-brane origin

The metric for the D4-brane is

$$ds_{10}^2=H^{3/8}(dt^2+dx_i^2)+H^{5/8}(dr^2+r^2d\mathrm{\Omega }_4^2),$$
(3.1)
where

$$H=1+\frac{R^3}{r^3},$$
(3.2)
and $`i=1,..,4`$. In the near-horizon limit $`r<<R`$, the above metric can be expressed as

$$ds_{10}^2=(1+k|z|)^{9/4}(dt^2+dx_i^2+dz^2)+(1+k|z|)^{1/4}\frac{1}{k^2}d\mathrm{\Omega }_5^2,$$
(3.3)
where

$$\frac{r}{R}=(1+k|z|)^2,$$
(3.4)
The analysis is similar to the case of the D3-brane. An additional element is that the worldvolume direction $`x_4`$ is taken to be compact, in order for the worldvolume of the braneworld to have 1+3 dimensions. Rotational dynamics around $`x_4`$ are parametrized by the integer $`n`$. For simplicity, we will assume that all compact dimensions are of the same size. From (2.10) we find that

$$U(r)\frac{G_4M_1M_2}{r}\underset{n=0}{\overset{\mathrm{}}{}}\mathrm{e}^{nkr}\left(\left[1+\frac{b_0}{(kr)^4}+\mathrm{}\right]+\underset{\mathrm{}=1}{\overset{\mathrm{}}{}}\left[\frac{b_{\mathrm{}}}{(kr)^{2\sqrt{\mathrm{}^2+3\mathrm{}+9}+2}}+\mathrm{}\right]\right),$$
(3.5)
where $`b_0`$ and $`b_{\mathrm{}}`$ are constants of order one. The $`\mathrm{}=0,M=0`$ terms were found in . The terms in the second pair of square brackets are the dominant corrections for the rotational graviton modes ($`\mathrm{}>0`$). As can be seen, the dominant contribution from the $`\mathrm{}>0`$ modes is smaller than that of the s-wave Kaluza-Klein modes by a rather small factor of approximately $`1/(kr)^5`$. An additional element here, as opposed to the D3-brane origin, is the presence of massive bound states and additional Kaluza-Klein states, both due to the $`n>0`$ modes. These contribute Yukawa-like terms to the gravitational potential. In the case of an $`M5`$-brane, the result is the same as above, except that now there is the possibility of modes rotating around two compact worldvolume dimensions, so that the $`n`$ in (3.5) is replaced by $`\sqrt{n_1^2+n_2^2}`$.

## 4 D5-brane origin

The metric for the D5-brane is

$$ds_{10}^2=H^{1/4}(dt^2+dx_i^2)+H^{3/4}(dr^2+r^2d\mathrm{\Omega }_3^2),$$
(4.1)
where

$$H=1+\frac{R^2}{r^2},$$
(4.2)
and $`i=1,..,5`$. In the near-horizon limit $`r<<R`$, the above metric can be expressed as

$$ds_{10}^2=\mathrm{e}^{k|z|/4}(dt^2+dx_i^2+dz^2+\frac{1}{k^2}d\mathrm{\Omega }_3^2),$$
(4.3)
where

$$\frac{r}{R}=\mathrm{e}^{k|z|/2},$$
(4.4)
In order to yield a braneworld with 1+3 worldvolume dimensions, we take $`x_4`$ and $`x_5`$ to be compact. We take $`\mathrm{\Phi }=\mathrm{e}^{i(n_1x_4+n_2x_5)}P_{\mathrm{}}(S^3)\mathrm{e}^{k|z|/2}\psi (z)M(t,x_i)`$. For graviton modes on the metric (4.3) we obtain

$$_z^2\psi +\left[\left(\frac{1}{4}+\alpha ^2\right)k^2k\delta (z)\right]\psi =m^2\psi ,$$
(4.5)
where $`\alpha ^2\mathrm{}(\mathrm{}+2)+n_1^2+n_2^2`$. $`\mathrm{}`$, $`n_1`$ and $`n_2`$ are integers parametrizing the rotational dynamics around $`S^3`$, $`x_4`$ and $`x_5`$ respectively. For the massless modes, the solution is

$$\psi =\mathrm{e}^{\sqrt{\alpha ^2+1/4}k|z|}+A\mathrm{e}^{\sqrt{\alpha ^2+1/4}k|z|},$$
(4.6)
where

$$A\frac{2\sqrt{\alpha ^2+1/4}1}{2\sqrt{\alpha ^2+1/4}+1}.$$
(4.7)
For $`\mathrm{},n_1,n_2=0`$, this wave function is normalizable and there is a corresponding localized graviton. However, for all higher angular momentum modes, the wave function diverges for large $`z`$ and thus there is no permissible massless solution. Nonzero $`\mathrm{},n_1,n_2`$ effectively decrease the mass term in the wave equation (4.5). Thus, a higher angular momentum mode has a localized state on the brane world of mass $`m^2=(\mathrm{}(\mathrm{}+2)+n_1^2+n_2^2)k^2`$. Also, each mode has an equal mass gap between the bound state and a continuum of Kaluza-Klein modes.

The massive wave functions are given by

$$\psi _m=N_m(k\mathrm{sin}q|z|2q\mathrm{cos}q|z|),$$

$$N_m=\frac{1}{2\sqrt{\pi q(m^2\alpha ^2k^2)}},$$
(4.8)
$$q=\sqrt{m^2k^2(1/4+\alpha ^2)}.$$

From (2.10) we find that

$$U(r)\frac{G_4M_1M_2}{r}\underset{\mathrm{},n_1,n_2=0}{\overset{\mathrm{}}{}}\left[\mathrm{e}^{\alpha kr}+\frac{2}{\sqrt{\pi }}\left(1+4\alpha ^2\right)^{3/4}\mathrm{e}^{\sqrt{\alpha ^2+1/4}kr}\left(\frac{1}{(kr)^{3/2}}+\mathrm{}\right)\right],$$
(4.9)
The $`\mathrm{},n_1,n_2=0`$ terms were found in . As can be seen, the rotational graviton modes yield Kaluza-Klein potential terms whose dominant correction is smaller than that of the s-wave Kaluza-Klein terms by a relative Yukawa-like factor of approximately $`\mathrm{e}^{kr}`$. Notice that there are also contributions from rotational bound states, whose dominant term is smaller than that of the s-wave Kaluza-Klein modes by a relative Yukawa-like factor of $`\mathrm{e}^{kr/2}`$.

## 5 M5/M5 origin

The metric for two M5-branes intersecting over 1+3 dimensions is

$$ds_{11}^2=H_1^{1/3}H_2^{1/3}\left(dt^2+dx_i^2+H_1(dy_1^2+dy_2^2)+H_2(dy_3^2+dy_4^2)+H_1H_2(dr^2+r^2d\mathrm{\Omega }_2^2)\right),$$
(5.1)
where

$$H_{1(2)}=1+\frac{R_{1(2)}}{r},$$
(5.2)
and $`i=1,2,3`$. For $`R_1=R_2R`$ and in the near-horizon limit $`r<<R`$, the above metric can be expressed as

$$ds_{11}^2=\mathrm{e}^{\frac{2}{3}k|z|}(dt^2+dx_i^2+dz^2+\frac{1}{k^2}d\mathrm{\Omega }_2^2)+\mathrm{e}^{\frac{1}{3}k|z|}dy_j^2,$$
(5.3)
where

$$\frac{r}{R}=\mathrm{e}^{k|z|},$$
(5.4)
and $`j=1,..,4`$. We take $`y_j`$ to be compact. Notice that the metric elements corresponding to the compact M5-brane worldvolume directions $`y_j`$ are multiplied by a ”warp” factor, which increases rapidly with $`|z|`$. It will be shown how this ”warp” factor amplifies the effect that the graviton modes rotating about $`y_j`$ have on the corrections to the four-dimensional gravitational potential.

We take $`\mathrm{\Phi }=\mathrm{e}^{i(n_jy_j)}P_{\mathrm{}}(S^2)\mathrm{e}^{k|z|/2}\psi (z)M(t,x_i)`$, where we shall implicitly sum over $`j`$ hereafter. For graviton modes on the metric given by (5.3), we obtain

$$_z^2\psi +\left[(\mathrm{}+1/2)^2k^2+n_j^2k^2\mathrm{e}^{k|z|}k\delta (z)\right]\psi =m^2\psi ,$$
(5.5)
The solution is

$$\psi =N_m[((i\gamma +1)Y_{i\gamma }(2in)2inY_{i\gamma 1}(2in))J_{i\gamma }\left(2in\mathrm{e}^{k|z|/2}\right)+$$

$$(2inJ_{i\gamma 1}(2in)(i\gamma +1)J_{i\gamma }(2in))Y_{i\gamma }\left(2in\mathrm{e}^{k|z|/2}\right)],$$
(5.6)
where $`\gamma \sqrt{4m^2/k^2(2\mathrm{}+1)^2}`$ and $`n^2n_j^2`$. The higher $`\mathrm{}`$ modes contribute massive bound states at $`m^2=\mathrm{}(\mathrm{}+1)k^2`$ and equal mass gaps between the bound states and continua of Kaluza-Klein modes. The $`n>0`$ modes only contribute the latter. For $`m>(\mathrm{}+2)k`$, the Kaluza-Klein modes can be plane-wave normalized with

$$N_m=\frac{1}{2\sqrt{k}}\left[\mathrm{sinh}(\pi \gamma )\right|(i\gamma +1)Y_{i\gamma }(2in)2inY_{i\gamma 1}(2in)|^2+$$

$$\frac{1}{\mathrm{sinh}(\pi \gamma )}|2inJ_{i\gamma 1}(2in)(i\gamma +1)J_{i\gamma }(2in)|^2]^{1/2}.$$
(5.7)
From (2.10) we find that

$$U(r)\frac{G_4M_1M_2}{r}\underset{\mathrm{}=0}{\overset{\mathrm{}}{}}[\mathrm{e}^{\sqrt{\mathrm{}(\mathrm{}+1)}kr}+(\frac{2}{\sqrt{\pi }}(2\mathrm{}+1)^{3/2}+$$

$$\sqrt{2}\pi ^{3/2}(2\mathrm{}+1)\underset{j=1}{\overset{4}{}}\underset{n_j=1}{\overset{\mathrm{}}{}}n^2\left|\frac{J_1(2in)Y_0(2in)J_0(2in)Y_1(2in)}{2inY_1(2in)J_0(2in)}|^2\right)\mathrm{e}^{(\mathrm{}+1/2)kr}(\frac{1}{(kr)^{3/2}}+\mathrm{})],$$
(5.8)
For the D5 and M5/M5 origins, the contributions to the gravitational potential are the same for the $`\mathrm{},n=0`$ modes  but different for the rotational modes. As can be seen, the $`\mathrm{}>0`$ modes yield potential terms that are sub-leading to the s-wave terms by a relative Yukawa-like factor, as in the case of the D5-brane origin. Notice there are contributions from massive bound states as well as additional Kaluza-Klein modes.

The $`n>0`$ Kaluza-Klein modes contribute to the Newtonian potential on the same order in $`kr`$ as the s-wave Kaluza-Klein modes! In particular, if we consider the $`\mathrm{}=0`$ modes, we have

$$U(r)\frac{G_4M_1M_2}{r}\left[1+\frac{2}{\sqrt{\pi }}\beta \mathrm{e}^{kr/2}(\frac{1}{(kr)^{3/2}}+\mathrm{})\right],$$
(5.9)
where the factor $`\beta =1.2`$ is due to the presence of the $`n=1`$ modes. Of course this factor is larger if the length scale of the $`y_j`$ directions is larger than $`1/k`$. Modes with $`n>1`$ do not make such a significant contribution.

## 6 Conclusions

We have found that, for brane worlds originating from the near-horizon of a D$`p`$-brane where $`p=3,4,5`$, the dominant correction of the four-dimensional gravitational potential arises from s-wave Kaluza-Klein modes propagating in the infinite extra dimension. The effect of Kaluza-Klein modes rotating about compact dimensions is sub-leading to that of the s-wave Kaluza-Klein modes, regardless of whether the former states are bound to the brane world (with Yukawa-like factors) or are propagating in the infinite extra dimension. This result is perhaps to be expected.

On the other hand, for a brane world scenario embedded in the near-horizon region of an M5/M5 intersection, we have found that there is a ”warp” factor which amplifies the effect of graviton modes rotating about compact directions that lie in the worldvolumes of the M5-branes. In this case, in fact, the modification of the four-dimensional gravitational potential due to these rotational modes is on the same order as that due to the s-wave Kaluza-Klein states. This can be seen explicitly by the presence of $`\beta `$ in (5.9). It is rather surprising that, the effects of modes around compact directions can be considerably amplified by ”warp” factors which depend on an extra large or infinite dimension. More complicated geometries and compactification schemes could yield additional amplifications in the gravitational sector.

## Acknowledgments

We are grateful to Malcolm Fairbairn and Fernando Quevedo for useful discussions.


format: text

text:
# Lattice-Boltzmann and finite-difference simulations for the permeability for three-dimensional porous media

## I Introduction

Almost all investigations of porous media focus on the prediction of effective material properties such as fluid permeability, electric or thermal conductivity, or elastic constants sah95 ; hil95d . The knowledge or at least a reliable prediction of these properties is of great interest in a wide field of technical applications ranging from petroleum engineering KBDHLPS99 ; lak89 , to paper manufacturing Kop98 , and contaminant transport hel97 . The predictions are obtained either from approximate theories which link the physical properties to geometrical observables, from geometrical models for which the physical problem can be solved more easily, or from various cross property relations, which relate the parameter in question to other physical transport parameters.

In this context, the exact numerical calculation of transport properties serves three purposes: (i) testing and validation of theories and theoretical predictions, (ii) comparison of geometrical models, and (iii) testing of faithfulness of computerized tomographic imaging by comparing numerically calculated transport parameters with their experimental values.

Of particular interest for porous media is the permeability, and more precisely its fluctuations. These fluctuations are important because they dominate the large scale permeability. For this reason it is important to collect as many micropermeametry measurements as possible. Experimental micropermeametry is costly and inaccurate. Hence exact numerical calculations are becoming an interesting alternative for studying fluctuations in permeability.

In this article we compare the permeabilities of a three-dimensional computerized tomographic image of Fontainebleau sandstone and its three physical and stochastic reconstruction models. We find that especially the stochastic models fail in reconstructing the fluid permeability of the original sandstone. This finding is in good agreement with conclusions drawn previously from a purely geometrical characterization of the same microstructures hil99c .

The exact numerical calculation of transport parameters for digitized three-dimensional samples remains a computationally demanding task, and only few studies exist which test the accuracy of such calculations. Here, we compare the results obtained by means of a finite difference (FD) method and a lattice-Boltzmann (LB) algorithm. We begin our study by calibrating both simulation methods against exact solutions of the Stokes equation for straight tubes and cubic arrays of spheres. These calculations also serve to compare the speed of both methods. We then proceed to apply both algorithms to the experimental sample and its three models.

## II Definition of the problem

The problem to be solved is that of slow laminar flow through a three-dimensional porous medium on a microscopic level. The three-dimensional microstructure of a two-phase porous medium $`𝕊`$ consisting of a pore phase $``$ and a matrix or rock phase $`𝕄`$ with $`𝕊=𝕄`$ is described in detail by the characteristic function $`\chi \text{}_𝔾`$ of a single phase $`𝔾\{𝕄,\}`$ with

$$\chi \text{}_𝔾(\stackrel{}{x})=\{\begin{array}{cc}1\hfill & \text{for}\stackrel{}{x}𝔾\hfill \\ 0\hfill & \text{for}\stackrel{}{x}𝔾\hfill \end{array}.$$
(1)
In the following $`\stackrel{}{x}`$ is the position vector of a cubic lattice, $`\stackrel{}{x}=ax_1\stackrel{}{e}_1+ax_2\stackrel{}{e}_2+ax_3\stackrel{}{e}_3`$, with $`x_i=0,1,\mathrm{},M_i1`$, the unit vectors $`\stackrel{}{e}_i`$ of the Cartesian coordinate system, and the grid spacing $`a`$. The total number of lattice points is given by $`N=M_1M_2M_3`$.

The Reynolds numbers of interest in geophysical and petrophysical applications are usually much smaller than unity hil94c  and hence it suffices to solve the Stokes equation. In the pore space geometry described by the characteristic function $`\chi \text{}_{}(\stackrel{}{x})`$, the steady-state Stokes equation and the constraint of incompressibility read

$`\eta \mathrm{\Delta }\stackrel{}{v}(\stackrel{}{x})\stackrel{}{}p(\stackrel{}{x})`$ $`=`$ $`0;\stackrel{}{x},`$ (2)
$`\stackrel{}{}\stackrel{}{v}(\stackrel{}{x})`$ $`=`$ $`0;\stackrel{}{x}.`$ (3)
On the pore-matrix interface $``$ we apply no-flow boundary conditions,

$$\stackrel{}{v}(\stackrel{}{x})=0;\stackrel{}{x}.$$
(4)
Darcy’s law permits us to compute the $`j`$-th column of the macroscopic permeability tensor from the microscopic solution of the hydrodynamic problem Eqs. (2) - (4) for a pressure gradient applied along the $`\stackrel{}{e}_j`$ direction according to

$$k_{ij}=v_i(\stackrel{}{x})_{\stackrel{}{x}𝕊}\frac{\eta (M_j+1)a}{p_{\mathrm{in}}p_{\mathrm{out}}}.$$
(5)
Here the pressures at the inlet and the outlet surface are given by $`p_{\mathrm{in}}`$ and $`p_{\mathrm{out}}`$ as defined below in Eqs. (16) and (19) 111Due to the spatial discretization $`p_{\mathrm{in}}`$ and $`p_{\mathrm{out}}`$ are the pressure values at a distance $`a/2`$ from the sample surface. Hence, we have to divide by $`(M_j+1)a`$ to obtain the mean pressure gradient. , respectively, $`\mathrm{}_{\stackrel{}{x}𝕊}`$ denotes an average over all lattice points, and $`v_i=\stackrel{}{v}\stackrel{}{e}_i`$.

## III Samples

We investigate four different samples $`𝕊_{\text{EX}}`$, $`𝕊_{\text{DM}}`$, $`𝕊_{\text{GF}}`$, and $`𝕊_{\text{SA}}`$. The first sample $`𝕊_{\text{EX}}`$, or in abbreviated form EX, was obtained experimentally by means of computerized tomography from a core of Fontainebleau sandstone. This sandstone is a popular reference standard because of its exceptional chemical, crystallographic and microstructural simplicity BZ85 ; BCZ87 . Fontainebleau sandstone consists of crystalline quartz grains that have been eroded for long periods before being deposited in dunes along shore lines during the Oligocene, i.e. roughly 30 million years ago. It is well sorted containing grains of around $`200`$ $`\mu `$m in diameter. During its geological evolution, which is still not fully understood, the sand was cemented by silica crystallizing around the grains. Fontainebleau sandstone exhibits intergranular porosity ranging from $`0.03`$ to roughly $`0.3`$ BCZ87 . The computer-assisted microtomography was carried out on a micro-plug drilled from a larger original core. The original core from which the micro-plug was taken had porosity $`\varphi ^{}=0.1484`$, permeability $`k^{}=1.3`$ D, ($`1`$ D=0.987 $`\mu `$m2) and formation factor $`22.1`$ (dimensionless electrical resistivity dul92 ). The microtomographic dataset has dimension $`M_1\times M_2\times M_3=299\times 300\times 300`$ with a resolution of $`a=7.5\mu \mathrm{m}`$, and porosity $`\varphi =0.1355`$. The pore space $`_{\text{EX}}`$ is visualized in Fig. 1 of Ref. hil99c .

The three remaining samples are physical and stochastic reconstruction models for the Fontainebleau sample EX. All have the same lattice resolution, $`a=7.5\mu `$m, and approximately the same porosity. The “diagenesis model” DM tries to mimic the geological formation process of the natural sandstone in three steps: the sedimentation of spherical grains, the compaction of the sediment, and the simulation of quartz cement overgrowth. The sample dimensions are $`M_1\times M_2\times M_3=255\times 255\times 255`$.

The SA and GF samples are stochastic models with dimensions $`M_1\times M_2\times M_3=256\times 256\times 256`$. Both models reconstruct the porosity $`\varphi `$ and the two-point correlation function of the original sandstone EX. This implies the reconstruction of the specific surface $`S_V`$. However, due to problems in the reconstruction procedure of the Gaussian field method, the porosity and the specific surface of the GF model do not match exactly those of the original sandstone. We find $`\varphi =0.1354`$ for SA, and $`\varphi =0.1421`$ for GF.

For a more detailed description of the modeling procedures and a visualization of the microstructures, the reader is referred to Ref. hil99c  and the references therein. In the same article results of an extensive geometrical investigation of the four samples are presented, which use both classical geometric quantities but also concepts introduced in local porosity theory hil91d ; hil96 . The main findings in Ref. hil99c  were as follows:

1. None of the three models can reproduce the visual appearance of the original sandstone. EX shows a granular structure of the matrix phase where single sand grains can easily be identified. The matrix phase of DM is also clearly granular but with artificial, spherical grains. In both models the matrix and the pore phase are very well connected, and the pore-matrix interface is smooth. In contrast, the pore-matrix interface of the stochastic models SA and GF is very rough. Here, both phases are strongly scattered and exhibit isolated clusters.
2. The two-point correlation functions, $`\text{S}_2(r)=\chi \text{}(\stackrel{}{x}_1)\chi \text{}(\stackrel{}{x}_2)`$ with $`r=|\stackrel{}{x}_1\stackrel{}{x}_2|`$, of SA and GF show good agreement with the original sandstone except for minor deviations at small $`r`$ in the case of GF. The correlation function of DM clearly deviates from that of EX. Moreover, it shows strong anisotropy with respect to directions $`\stackrel{}{e}_1`$, $`\stackrel{}{e}_2`$ and $`\stackrel{}{e}_3`$. A modified correlation function giving the conditional probability to find two points in pore space that are also connected by a path inside the pore space was measured and discussed in Ref.hil00f . Small but significant differences exist between the samples. The experimental sample EX is more stable under the morphological operations of erosion and dilation hil00f ; hil00 .
3. The differences between the samples are most pronounced when comparing the geometrical connectivity of the pore space. As a measure for the geometrical connectivity, we use the total fraction of percolating cells $`p_3(L)`$ at scale $`L`$, introduced in local porosity theory hil00 , which is defined as the probability for a cubic subblock of size $`L`$ of the sample to percolate in all three coordinate directions $`\stackrel{}{e}_i`$. Here, percolation in direction $`\stackrel{}{e}_i`$ means that there exists a path lying entirely in the pore space, which spans from one face of the cubic subsample perpendicular to $`\stackrel{}{e}_i`$ to the opposite face. For EX and DM the curves of $`p_3(L)`$ nearly coincide while for the stochastic models SA and GF the curves of $`p_3(L)`$ fall well below that of the original sandstone (see Fig. 13 in Ref. hil99c ). Table 1 gives the values of $`p_3(L)`$ for $`L=60a`$. Again we find anisotropy in the DM model when we measure the probability $`p_{\stackrel{}{e}_i}(L)`$ for a subblock to percolate in directions $`\stackrel{}{e}_i`$.

Geometrical connectivity is an indispensable precondition for dynamical connectivity and physical transport. Hence, we expect to find a strong correlation between the total fraction of percolating cells $`p_3(L)`$ and the macroscopic permeability whose calculation is discussed next.

## IV The numerical methods

### IV.1 The finite difference method

#### IV.1.1 Algorithm

Numerically we obtain the solution of Eqs. (2) - (4) from the infinite time limit of the time dependent Stokes problem using an iterative pressure-correction algorithm B:pat80 ; B:hir88 . Discretization in time of the time-dependent Stokes equation yields

$`{\displaystyle \frac{\stackrel{}{v}^{n+1}(\stackrel{}{x})\stackrel{}{v}^n(\stackrel{}{x})}{\mathrm{\Delta }t}}`$ $`=`$ $`\eta \mathrm{\Delta }\stackrel{}{v}^n(\stackrel{}{x})\stackrel{}{}p^{n+1}(\stackrel{}{x}),`$ (6)
$`\stackrel{}{}\stackrel{}{v}^{n+1}(\stackrel{}{x})`$ $`=`$ $`0,`$ (7)
where the superscript $`n`$ denotes the iteration step. In our case the discretized time derivative on the left hand side of Eq. (6) has no physical meaning. In the long-time limit $`\stackrel{}{v}^{n+1}(\stackrel{}{x})=\stackrel{}{v}^n(\stackrel{}{x})`$ holds, and we recover Eq. (2).

Given the solutions $`\stackrel{}{v}^n`$ and $`p^n`$ at iteration step $`n`$, an approximate solution $`\stackrel{}{v}^{}`$ for the velocity field is obtained from

$$\frac{\stackrel{}{v}^{}(\stackrel{}{x})\stackrel{}{v}^n(\stackrel{}{x})}{\mathrm{\Delta }t}=\eta \mathrm{\Delta }\stackrel{}{v}^n(\stackrel{}{x})\stackrel{}{}p^n(\stackrel{}{x}).$$
(8)
Subtracting Eq. (8) from Eq. (6), we find that

$$\frac{\stackrel{}{v}^{n+1}(\stackrel{}{x})\stackrel{}{v}^{}(\stackrel{}{x})}{\mathrm{\Delta }t}=\stackrel{}{}\left(p^{n+1}(\stackrel{}{x})p^n(\stackrel{}{x})\right).$$
(9)
On the right hand side of this equation appears the pressure correction $`p^{}(\stackrel{}{x})p^{n+1}(\stackrel{}{x})p^n(\stackrel{}{x})`$. Applying the $`\stackrel{}{}`$ operator to Eq. (9), and using the incompressibility constraint Eq. (7) we obtain Poisson’s equation for $`p^{}`$,

$$\mathrm{\Delta }p^{}(\stackrel{}{x})=\frac{1}{\mathrm{\Delta }t}\stackrel{}{}\stackrel{}{v}^{}(\stackrel{}{x}).$$
(10)
Thus, we arrive at the following algorithm:

1. Let $`\stackrel{}{v}^n`$ and $`p^n`$ be the solution of the velocity and the pressure field, respectively, at iteration step $`n`$ with the maximum absolute error $`ϵ^n=\mathrm{max}_\stackrel{}{x}|\eta \mathrm{\Delta }\stackrel{}{v}^n(\stackrel{}{x})p(\stackrel{}{x})|`$. From $`\stackrel{}{v}^n`$ and $`p^n`$ an approximate solution $`\stackrel{}{v}^{}`$ of the velocity field is calculated using Eq. (8).
2. Using the definition of the pressure correction, a new pressure $`p^{n+1}=p^n+p^{}`$ is obtained from a solution of Eq. (10). This part of the algorithm consumes most of the computation time, because Eq. (10) has to be solved for each iteration step. However, we found that it suffices to solve Eq. (10) only up to an error

$$\underset{\stackrel{}{x}}{\mathrm{max}}|\mathrm{\Delta }p^{}(\stackrel{}{x})\frac{1}{\mathrm{\Delta }t}\stackrel{}{}\stackrel{}{v}^{}(\stackrel{}{x})|\gamma ϵ^n$$

where $`\gamma `$ is an empirical factor. For the calculations presented here, a value of $`\gamma `$ in the range $`0.01<\gamma <0.1`$ seems to be appropriate.

We use a successive over-relaxation method to solve Eq. (10). Of course it would be desirable to use more sophisticated methods such as e.g. a multigrid method, but we could not find a general procedure to restrict the microstructure of the porous medium to a coarser grid without changing the topology of the pore space.
3. From $`p^{}`$ and $`\stackrel{}{v}^{}`$ a new velocity $`\stackrel{}{v}^{n+1}`$ is calculated using Eq. (9). The algorithm terminates when $`ϵ^{n+1}`$ is smaller than some given value.

The equations are spatially discretized using a marker-and-cell (MAC) grid HW65 . The pressure values are placed at the centers of the grid cells. They coincide with the lattice points of the discretized characteristic function $`\chi \text{}_{}`$. On each face of a grid cell the velocity component perpendicular to this face is located. The pore-matrix interface $``$ follows the surface of the cubic grid cells. For velocity components perpendicular to the interface the boundary condition Eq. (4),

$$v_{}(\stackrel{}{x})=0,$$
(11)
is implemented exactly. For parallel velocity components at distance $`a/2`$ from the interface, Eq. (4) is fulfilled to second order accuracy,

$$v_{}(\stackrel{}{x})=v_{}(\stackrel{}{x}+a\stackrel{}{e}_{})+O(a^2),$$
(12)
in which the interface is located as shown in Fig. 1. Inserting Eq. (4) into the Stokes equation Eq. (2) one finds as the boundary conditions for the pressure

$$_{}p(\stackrel{}{x})|_{}=0.$$
(13)
Using Eqs. (11) - (13), the $`\mathrm{\Delta }`$ operators in Eqs. (8) and (10) take the form

$$\begin{array}{cc}\hfill \mathrm{\Delta }v_i(\stackrel{}{x})& =v_i(\stackrel{}{x}a\stackrel{}{e}_i)+v_i(\stackrel{}{x}+a\stackrel{}{e}_i)2v_i(\stackrel{}{x})\hfill \\ & +\underset{ji}{}[\chi \text{}_{}(\stackrel{}{x}+a\stackrel{}{e}_j)\chi \text{}_{}(\stackrel{}{x}+a\stackrel{}{e}_j+a\stackrel{}{e}_i)[v_i(\stackrel{}{x}+a\stackrel{}{e}_j)+v_i(\stackrel{}{x})]\hfill \\ & +\chi \text{}_{}(\stackrel{}{x}a\stackrel{}{e}_j)\chi \text{}_{}(\stackrel{}{x}a\stackrel{}{e}_j+a\stackrel{}{e}_i)[v_i(\stackrel{}{x}a\stackrel{}{e}_j)+v_i(\stackrel{}{x})]4v_i(\stackrel{}{x})];\stackrel{}{x},\hfill \end{array}$$
(14)
and

$$\begin{array}{cc}\hfill \mathrm{\Delta }p^{}(\stackrel{}{x})=\underset{i}{}[& \chi \text{}_{}(\stackrel{}{x}+a\stackrel{}{e}_i)\left(p^{}(\stackrel{}{x}+a\stackrel{}{e}_i)p^{}(\stackrel{}{x})\right)\hfill \\ & \chi \text{}_{}(\stackrel{}{x}a\stackrel{}{e}_i)(p^{}(\stackrel{}{x}a\stackrel{}{e}_i)p^{}(\stackrel{}{x}))];\stackrel{}{x},\hfill \end{array}$$
(15)
where we used second order accurate central differences to discretize the spatial derivatives.

#### IV.1.2 Boundary conditions

On the sample surface an additional outer layer of grid cells, the so-called shadow row, is added. It provides the neighboring pressure and velocity values needed for the evaluation of Eqs. (14) and (15) for the grid cells on the sample surface. The pressure and velocity values of the shadow row are set according to the macroscopic boundary conditions. Let $`\stackrel{}{e}_j`$ denote the direction of the applied pressure gradient. For grid cells of the shadow row on the outflow boundary of the sample, i.e. $`\stackrel{}{x}\{\stackrel{}{x}:x_j=M_j\}`$, we choose

$`p(\stackrel{}{x})`$ $`=`$ $`p_{\mathrm{out}},`$ (16)
$`v_i(\stackrel{}{x})`$ $`=`$ $`0\text{for }ij,`$ (17)
$`v_j(\stackrel{}{x})`$ $`=`$ $`v_j(\stackrel{}{x}a\stackrel{}{e}_{}),`$ (18)
as the boundary conditions. For $`\stackrel{}{x}\{\stackrel{}{x}:x_j=1\}`$, i.e. on the inflow boundary, we set

$`p(\stackrel{}{x})`$ $`=`$ $`p_{\mathrm{in}},`$ (19)
$`v_i(\stackrel{}{x})`$ $`=`$ $`0\text{for }ij,`$ (20)
$`v_j(\stackrel{}{x})`$ $`=`$ $`v_j(\stackrel{}{x}a\stackrel{}{e}_{}),`$ (21)
as the boundary conditions. The last condition has to be introduced directly into Eq. (14) because the position $`\stackrel{}{x}a\stackrel{}{e}_j`$ lies outside the simulation lattice. In the simulations reported below we used always $`p_{\mathrm{in}}=1`$, $`p_{\mathrm{out}}=1`$ unless indicated otherwise.

On the remaining sample surfaces the grid cells of the shadow row are assigned to the matrix phase. The velocity and pressure values are set to zero.

### IV.2  The lattice-Boltzmann method

In this section we introduce the lattice-Boltzmann method used here, in particular the LBGK (lattice-Bhatnagar-Gross-Krook) model. We discuss then the basic hydrodynamics of the model and the relevant boundary conditions. The numerical accuracy of the lattice-Boltzmann results for permeability will be considered in terms of finite-size effects.

#### IV.2.1 Lattice-Boltzmann hydrodynamics

The lattice-Boltzmann method Qia92 ; Ben92 ; Rot97 ; CD98  is a mesoscopic approach for computational fluid dynamics in which the basic idea is to solve a discretized Boltzmann equation. The macroscopic dynamics of the system can be shown to obey the Navier-Stokes equation. One of the most successful applications of the method has been to flow in porous media Fer95 ; Mar96 ; Kop98 .

In this method the fluid is modeled by particle distributions that move on a regular lattice. In our implementation each lattice point is connected to its nearest and next nearest neighbors. Together with a rest particle, each lattice point is then occupied by 19 different particles (the D3Q19 model). At each time step particles propagate to their adjacent lattice points, and re-distribute their momenta in the subsequent collisions. The dynamics of the LBGK model is given by the equation Qia92 ; Ben92

$$f_i(𝐫+𝐜_i,t+1)=f_i(𝐫,t)+\frac{1}{\tau }[f_i^{\mathrm{eq}}(𝐫,t)f_i(𝐫,t)],$$
(22)
where $`𝐜_i`$ is a vector pointing to an adjacent lattice site, $`f_i(𝐫,t)`$ is the density of the particles moving in the $`𝐜_i`$-direction, $`\tau `$ is the BGK relaxation parameter, and $`f_i^{\mathrm{eq}}(𝐫,t)`$ is the equilibrium distribution towards which the particle populations are relaxed. Hydrodynamic quantities like density $`\rho `$ and velocity $`𝐮`$ are obtained from the velocity moments of the distribution $`f_i`$ in analogy with the kinetic theory of gases. The equilibrium distribution can be chosen in many ways. A common choice is

$$f_i^{\mathrm{eq}}=t_i(1+\frac{1}{c_s^2}(𝐜_i𝐮)+\frac{1}{2c_s^4}(𝐜_i𝐮)^2\frac{1}{2c_s^2}u^2),$$
(23)
in which $`t_i`$ is a weight factor that depends on the length of the link vector $`𝐜_i`$, and $`c_s`$ is the speed of sound in the fluid. The weights $`t_i`$ we choose here in accordance with the 19-link LBGK model, and they are $`\frac{1}{3},\frac{1}{18}`$ and $`\frac{1}{36}`$ for the rest particle and the particles moving to the nearest and next-nearest neighbor sites, respectively. The speed of sound is $`c_s=\frac{1}{\sqrt{3}}`$ for this model, and the kinematic viscosity of the simulated fluid is $`\eta =\frac{2\tau 1}{6}`$. (Here and in the following, lattice units are always used if the units are not specified.) The fluid pressure is given by

$$p(𝐫,t)=c_s^2(\rho (𝐫,t)\overline{\rho })c_s^2\mathrm{\Delta }\rho (𝐫,t),$$
(24)
where $`\overline{\rho }`$ is the mean density of the fluid.

The Stokes equation, Eq. (2), is produced directly by the linearized lattice-Boltzmann method, in which the quadratic velocity terms in the equilibrium distribution function Eq. (23) are neglected. To be consistent with the finite difference method, we use in what follows the linearized lattice-Boltzmann method if not stated otherwise.

#### IV.2.2 Boundary conditions

The physical boundary condition at solid-fluid interfaces is the no-flow condition Eq. (4), which in lattice-Boltzmann simulations is usually realized by the so-called bounce-back rule Gal97 ; He97 . In this approach the momenta of the particles that meet a solid wall are simply reversed.

In simple shear flows the bounce-back condition assumes that the location of the wall is exactly halfway between the last fluid point and the first wall point. In more complicated cases the no-flow boundary lies somewhere in between these two points, the exact place depending on the relaxation parameter and the geometry of the system He97 ; Kan99 . In Poiseuille flow, e.g., the bounce-back rule gives velocity fields that deviate from the exact solution, for no-flow boundaries at exactly halfway between the last fluid point and the first solid point, by He97

$$\mathrm{\Delta }u=u_{\mathrm{sim}}u_{\mathrm{exact}}=u_{\mathrm{max}}\frac{48\eta ^24\eta 1}{L^2},$$
(25)
where $`u_{\mathrm{sim}}`$ and $`u_{\mathrm{exact}}`$ are the simulated and the exact velocities, respectively, $`u_{\mathrm{max}}`$ is the velocity at the center of the channel, and $`L`$ is channel width. This implies that the simulated permeability will depend somewhat on viscosity especially at low discretization levels. This viscosity dependence can practically be eliminated by using the so-called second order boundaries, in which case the desired location of the no-flow boundary is determined by extrapolating the distribution function from the last fluid points. Some of these more sophisticated solid-fluid boundaries are restricted to regular geometries Sko93 ; Nob95 , but there are also general boundary-fitted models Fil98 ; Che98  available. For practical simulations the bounce-back boundary is however very attractive, because it is a simple and computationally efficient method for imposing no-flow conditions on irregularly shaped walls. Also, the error created by the bounce-back boundary does not destroy the spatial second-order convergence of the method Kan99 ; Mai96 .

In simulating fluid flow it is important that the velocity and pressure boundary conditions of the system have been imposed in a consistent way. However, general velocity and pressure boundaries are still under development for the lattice-Boltzmann method Mai96 ; Che96 ; Gin96 . So far in most of the practical simulations a body force has been implemented Fer95 ; Kan99 ; Bui2000  instead of pressure or velocity boundaries.

When the body force is used, the pressure gradient acting on the fluid is replaced with a uniform external force. The use of a body force is based on the assumption that, on the average, the effect of an external pressure gradient is constant throughout the system, and that it can thus be replaced with a constant force that adds at every time step a fixed amount of momentum on the fluid points. Conditions that are close to pressure boundaries can be obtained by averaging the velocity and pressure fields over the planes of the inlet and outlet of the simulated system Kan99 .

During one iteration step, the fluid momentum oscillates in the stationary state by an amount given to each fluid point by the body force. For this reason the fluid velocity is now defined as the average of the pre-collision and post-collision values Shan93 ; Shan95 .

Pressure fields generated by the body force are obtained from the effective pressure $`p_{\mathrm{eff}}`$,

$$p_{\mathrm{eff}}(𝐫,t)=c_s^2\mathrm{\Delta }\rho (𝐫,t)\overline{\rho }gx,$$
(26)
where $`x`$ is the distance from the inlet of the system measured in the flow direction, and $`g`$ is the acceleration the body force gives to the fluid.

It is a well-known fact that, due to staggered invariants, the fluid momentum may oscillate in a time scale of a few time steps CD98 , even in the stationary state. In open areas this effect is usually unimportant, but in closed pores this effect may become visible as the fluid momentum may oscillate around zero, with a magnitude determined by the body force. This effect may lead to some corruption of the fluid-velocity distributions as can be seen in Fig.6. Notice that staggered momenta can be eliminated by averaging the quantities over a few time steps.

Notice finally that the diagonal links allow the fluid to leak to neighboring lattice points which have only a single edge in common. For this reason the ’standard’ lattice-Boltzmann model is not expected to be accurate very close to the percolation threshold. For the 3D checkerboard structure, e.g., we found that the permeability of the system was about 0.036 lattice units for all the six lattice resolutions that were used, although the structure is not percolating. If better accuracy is needed, diagonal leaks can be eliminated by applying the bounce-back rule on such diagonal links which actually cross a solid boundary (like diagonals in the checkerboard structure).

#### IV.2.3 Finite-size effects and the saturation time

The accuracy of lattice-Boltzmann simulations depends on the ratio of the mean free path $`\lambda _{\mathrm{mfp}}`$ of the fluid particles to the representative size $`\lambda _0`$ of the obstacles and pores Rot97 ; Fer95 ; Gin96 . The simulated flow field does not describe the true hydrodynamic behavior unless this ratio is small. For increasing $`\lambda _{\mathrm{mfp}}/\lambda _0`$ ratio Knudsen-flow behavior is found, which is also true in real fluids Sch57 . These effects must always be considered when lattice-Boltzmann simulations are performed. In this way the maximum size of the lattice spacing can be estimated together with the accuracy of the simulations.

Finite-size effects restrict to some extent the use of LB-methods based on regular lattices. In porous media close to the percolation threshold, e.g., many pores are very small, and very big lattices may be needed for realistic simulations. It is still an open question whether the finite-size effects are always dominated by the minimum pore size or the average pore size. The effect can be estimated by simulating the system with several different lattice spacings, but occasionally it is difficult to distinguish the finite-size effects from other sources of numerical error.

Practice has shown Fer95 ; Kop98  that smaller values of the relaxation parameter $`\tau `$ tend to decrease the finite-size effects (see also our simulations below). Equation (25) can be used to explain this: due to Knudsen-flow effects, low-discretization simulations regularly give too high permeabilities, whereas decrease of $`\tau `$ has the opposite effect down to $`\tau =0.625`$, at least for tube flows. On the other hand, the lattice-Boltzmann algorithm may become CD98  unstable with values of $`\tau `$ close to $`0.5`$. In practical permeability simulations the relaxation parameter has usually been chosen to be bigger than 0.6. The effect of $`\tau `$ on the behavior of the lattice-Boltzmann model is thus quite complicated and not yet fully understood.

In permeability simulations a simple dimensional analysis shows that, with a constant body force, the saturation time $`t_{\mathrm{sat}}`$ needed to reach the steady state is of the form

$$t_{\mathrm{sat}}R_{\mathrm{pore}}^2/\eta ,$$
(27)
where $`R_{\mathrm{pore}}`$ is the characteristic length of the void pores in the system. For systems with high porosity, the saturation times can therefore be very long. In some cases, tens of thousands of time steps may be needed. It is thus evident that a constant body force may be computationally inefficient, especially when one is only interested in the steady-state solution. The saturation time can be reduced by using e.g. the iterative momentum-relaxation (IMR) method, where the applied body force is adjusted during the iteration in a definite relation to the change of the fluid momentum during iteration steps Kan99 . For other ways to reduce the saturation time see Ver99 .

## V Results

### V.1 Tube with quadratic cross section

One of the few cases for which the analytical solution of the hydrodynamic problem Eqs. (2) - (4) is known is Poiseuille flow, the flow through a linear tube with constant cross section. We will consider a tube with quadratic cross section, because here the geometry can be discretized on a cubic lattice without discretization error.

We consider a tube directed along the $`\stackrel{}{e}_1`$ direction with quadratic cross section of side length $`B=32a`$. We compare the velocity component $`v_1(x_2,x_3)`$ with its exactly known reference value $`v_1^{\mathrm{ref}}(x_2,x_3)`$ given in Ref. B:wie74 . Figure 2a shows the relative error $`(v_1v_1^{\mathrm{ref}})/v_1^{\mathrm{ref}}`$ for the LB solution with relaxation parameter $`\tau =0.688`$. In Fig. 2b we show the same relative error for the LB solution with $`\tau =1.0`$ (lower surface) and the FD solution (upper surface). The pressure gradient in the FD simulation was $`2/33`$ while in the LB simulations it was around $`10^4`$.

Around the center of the tube the analytical flow profile is very well recovered. Near the boundaries we find deviations which are biggest in the corners. The LB solution with $`\tau =0.688`$ underestimates the reference value $`v^{\mathrm{ref}}`$ while the solution with $`\tau =1.0`$ overestimates the true value. Hence, the relaxation parameter $`\tau `$ or equivalently the viscosity $`\eta `$ could be adjusted to find better agreement with the analytical velocity field. From Fig. 2 one expects to find a value $`0.688<\tau <1.0`$ for which the numerical solution closely matches the analytical velocity profile.

The computation time needed by the FD method which terminated when $`\mathrm{max}_\stackrel{}{x}|\mathrm{\Delta }\stackrel{}{v}(\stackrel{}{x})p(\stackrel{}{x})|<10^8`$ was 951 s on a DEC Alpha workstation. In LB simulations, the relative error of permeability was below $`10^5`$ in 754 s on a Cray T3E for $`\tau =1.0`$, but the simulations were continued for over 5000 s to make sure the saturation of the velocity fields.

### V.2 Cubic array of spheres

To test the accuracy and efficiency of our two algorithms in a more complicated geometry with narrow constrictions we computed flow past a cubic array of spheres. This problem has become a reference system for checking hydrodynamic algorithms because accurate reference values for the permeability, and the drag coefficient, are available over a wide range of porosities LH89 ; lad88 .

The solution of this problem proceeds by solving the problem in a single unit cell of the cubic lattice. We generated six different unit cells of size $`L\{20a,36a,56a,63a,71a,89a\}`$. A sphere is placed at the center of each cell whose radius is chosen such that the porosity matches as close as possible to $`\varphi =0.15`$. Thus, the porosity is close to the porosity of the sandstones investigated later.

The FD solution of the flow field was computed using periodic boundary conditions on those faces of the unit cell that are parallel to the macroscopic flow direction. On the faces of the unit cell perpendicular to the macroscopic flow we applied the conditions Eqs. (16) - (21) with the standard pressure gradient $`p_{\mathrm{in}}=1`$, $`p_{\mathrm{out}}=1`$. In the LB solution the flow field was computed using periodic boundary conditions in all directions. This difference of the boundary conditions arises from the fact that in the LB simulations the density fluctuations around an average density are calculated while in the FD simulations the pressure field enters directly.

Once the velocity field was known we calculated the permeability from Eq. (5). Following LH89  we then utilize the expression

$$\frac{k}{R^2}=\frac{1}{6\pi C_D}\left(\frac{L}{R}\right)^3,$$
(28)
to obtain the reference value $`k_{\mathrm{ref}}`$ of the permeability from the drag coefficient $`C_D`$ given in Ref. LH89 . The radius $`R(\varphi ,L)`$ of the spheres depends on $`\varphi `$ and $`L`$ and is given implicitly by the expression

$$\varphi =\frac{8\pi }{3}\left(\frac{R}{L}\right)^33\pi \left(\frac{R}{L}\right)^2+\frac{\pi }{4}+1.$$
(29)
To calculate the reference value $`k_{\mathrm{ref}}`$ we solve this equation and find $`(L/R)1.6011`$ for $`\varphi =0.15`$. Using the drag coefficient $`C_D=1.020\times 10^3`$ LH89 , we find $`k_{\mathrm{ref}}/R^20.0002135`$.

In Fig. 3 the relative error $`(kk_{\mathrm{ref}})/k_{\mathrm{ref}}`$ of the permeability is plotted as a function of the linear dimensionless system size $`L/a`$. With increasing resolution the results of both methods converge to each other and the error predominantly decreases. In the lattice-Boltzmann simulations the relaxation parameter $`\tau =1.0`$ is seen to give regularly better results than $`\tau =0.688`$. It thus appears that for $`\tau =1.0`$ the effective location of the no-flow boundary is more satisfactory than for $`\tau =0.688`$, and, consequently, the relative error is smaller for $`\tau =1.0`$ even though the finite-size effects are similar in both cases. The deviation for $`L=89a`$ is negative and varies between $`3`$% and $`6`$% depending on the method.This discrepancy might possibly result mainly from discretization errors as there is a similar oscillatory trend in the LB as well as the FD results. Further work on larger systems is however needed to answer the question whether the discrepancy might also result from other sources.

The curves in Fig. 3 are all nonmonotonous. This results most likely from the discretization of the cross-sectional area of the pore throats between the spheres. The curve of the discretized cross-sectional area as a function of $`L/a`$ shows a similar nonmonotonic behavior.

Besides the accuracy, the demand of computation time is the second important characteristic of a numerical method. Comparison of the computation time of the LB and the FD algorithms is difficult. The FD method iterates the physical velocity field $`\stackrel{}{v}`$ and the physical pressure field $`p`$. The iteration is terminated when $`\stackrel{}{v}`$ and $`p`$ fulfill the Stokes equation Eq. (2) with a predefined accuracy. In the LB method on the other hand, the particle distributions $`f_i`$ are iterated. The velocity and pressure fields are calculated from the final density distribution. In our implementation there is no deterministic stopping criterion, although such criterion could be included. In practice the calculations were terminated after a fixed number of iteration steps. This number was determined for each system by comparing runs of different length.

All calculations with the FD code were performed on a Cray T3E-900/512 at the HLRS computing center of the University of Stuttgart with a peak performance of 461 Gflops. The LB code was run on a Cray T3E-750/512 at Center for Scientific Computing (CSC) in Espoo, Finland, with a peak performance of 384 Gflops. In order to compare the run times for the two codes we took the actual time required to execute the program corrected by the ratio of peak performances. The run time required for the LB code was calculated from the number of iterations multiplied by a conversion factor. The conversion factor was 0.039448 for $`L=20a`$, 0.208096 for $`L=36a`$, 0.628828 for $`L=56a`$, 0.885556 for $`L=63a`$, 1.275176 for $`L=71a`$, and 2.21636 for $`L=89a`$. It was determined by the wall time spent for one iteration step computed from averages over several 100-step iterations. The memory requirements of the two algorithms are different. The FD algorithm requires to store 8 numbers per lattice node in the version used here. The D3Q19 model used for the LB-algorithm requires to store 19 numbers per lattice node

In Fig. 4 we compare the time evolution of the numerical value of the permeability $`k`$ for different system sizes. Plotted on the $`x`$ axis is the total time in seconds needed on two Cray-processors. For $`L=20a`$ and $`L=36a`$, both methods reach the final value of $`k`$ in approximately the same time. For large $`L`$ the LB method seems to be faster. Notice that the results of the LB simulations shown in Fig. 4 were performed for $`\tau =0.688`$. For $`\tau =1.0`$ the simulations were about $`45`$ % faster.

The convergence of $`k(t)`$ towards its asymptotic value is monotonic for the FD method, while in the LB case $`k(t)`$ shows strong oscillations for all $`L`$. The reason for these oscillations is probably the slight compressibility error inherent in the model Rot97 ; CD98 . The horizontal line in Fig. 4 gives the reference value $`k_{\mathrm{ref}}`$ towards which the asymptotic values of both algorithms converge with $`L`$.

We also compared the permeabilities given by the Navier-Stokes and Stokes (linearized) versions of the LB method for the cubic arrays of spheres. Very much as expected, the two results were the same within the first seven digits (body force about $`10^4`$).

### V.3 Three-dimensional sandstones

We now apply both algorithms to the solution of the hydrodynamic problem Eqs. (2) - (4) within the irregular pore space geometries of the whole experimental sample EX and the model samples DM, GF, and SA.

The FD algorithm used the boundary conditions as described above. The iteration scheme was terminated when the condition $`\mathrm{max}_\stackrel{}{x}|\mathrm{\Delta }\stackrel{}{v}(\stackrel{}{x})p(\stackrel{}{x})|<10^6`$ for the dimensionless left hand side of Eq. (2) was fulfilled for the first time. Thus, the relative error $`ϵ(k_{ii})/k_{ii}`$ of the diagonal elements of the permeability tensor is estimated to be smaller than $`0.012`$ in the case of EX, and $`ϵ(k_{ii})/k_{ii}<0.36`$ in the case of SA. The relative error for the samples DM and GF lies in between these two extreme values.

In the LB simulations no flow boundary conditions were applied on the sample surfaces parallel to the main flow direction. At the inlet and outlet (i.e. the sample surfaces perpendicular to the main flow direction) an additional fluid layer with a thickness of $`1921`$ lattice spacings was added and then periodic boundary conditions were applied. The body force did not act in the additional fluid layer. These additional fluid layers increased the total number of lattice points by about $`8`$ % in comparison with the FD method. The relaxation parameter was $`\tau =0.688`$. The simulation stopped after a predefined number of iterations which was estimated to suffice for the permeability to converge.

In Table 2 we give the components of the permeability tensors for all four samples and for both algorithms.

The permeability results confirm the predictions from a previous purely geometrical analysis based on local porosity theory hil99c . The analysis in Ref. hil99c  emphasized the importance of local connectivity. The permeabilities are strongly correlated with the geometrical connectivity of the pore space, measured by means of the total fraction of percolating cells $`p_3(L)`$. In accordance with our discussion of $`p_3(60a)`$ given in Table 1 of Section III, we find that the permeability of the original sandstone EX and that of the process model DM are in good agreement, while the permeabilities of the stochastic models GF and SA are an order of magnitude smaller. It seems as if the stochastic reconstruction models cannot reproduce the high degree of geometrical connectivity present in the original sandstone. The reconstructed two-point correlation function lacks information about the geometrical connectivity of the pore space. A correct description of the geometrical connectivity is however an indispensable precondition for the correct dynamical connectivity that determines the transport properties.

We now proceed to compare the numerically obtained value of the permeability $`k_{\text{EX}}`$ of the Fontainebleau sandstone EX with the experimental value $`k^{}=1.3`$ D. Such a comparison requires a correction due to the difference between the porosity of the EX sample and the porosity $`\varphi ^{}=0.1484`$ of the original core sample on which the experiment was performed. There exists a well-known experimental correlation between porosity and permeability of Fontainebleau sandstone BCZ87 . This correlation is usually approximated in the form

$$k=A\varphi ^b,$$
(30)
in which $`A`$ and $`b`$ are constants. In the porosity range of interest $`\varphi 0.13\mathrm{}0.15`$, this correlation has $`b4`$, with however a large uncertainty due to the scatter in the measured results. Hence, we can extrapolate the numerically determined permeabilities $`\overline{k}=(k_{11}+k_{22}+k_{33})/3`$ into the prediction

$$\overline{k}^{}=\overline{k}\left(\frac{\varphi ^{}}{\varphi }\right)^b,$$
(31)
where $`\varphi ^{}`$ is the previously defined porosity of the core sample and the constant $`A`$ has dropped out. From Eq. (31) we obtain $`\overline{k}^{}=1150\mathrm{mD}`$ for the FD method and $`\overline{k}^{}=1015\mathrm{mD}`$ for the LB algorithm. These values are surprisingly close to the experimental value $`k^{}=1300\mathrm{mD}`$. Such an excellent agreement is not common. This will in fact be seen in the following when we determine the permeability of a subsample.

We also checked by the LB methods the difference between the Navier-Stokes and Stokes permeability of sample EX. The relative difference was found to be $`0.00036`$ for the parameters specified above, with an about $`17\%`$ longer simulation time in the full Navier-Stokes case. The smallness of this difference only demonstrates that, for the small pressure differences considered here, we indeed are in the Stokes regime.

### V.4 Fine Graining

We investigate the permeability of subsamples of the original samples for two reasons. Firstly this gives us an estimate of the magnitude of permeability fluctuations, and secondly it allows us to estimate the dependence of numerically obtained permeability on the lattice resolution $`a`$. The subsamples have dimensions of $`M_1\times M_2\times M_3=100\times 100\times 100`$. To test the dependence of $`k`$ on $`a`$, we apply a fine graining scheme man01 . The fine graining algorithm replaces each lattice point by $`n\times n\times n`$ lattice points of the same phase with $`n\{2,3,4\}`$. Thus, we get systems of dimensions $`M_1\times M_2\times M_3=n100\times n100\times n100`$ with lattice spacing $`a_n=\frac{7.5}{n}\mu \mathrm{m}`$. For each system the hydrodynamic problem Eqs. (2) - (4) is solved, and the permeability $`k^{(n)}`$ is determined.

Figure 5 shows $`k^{(n)}`$ for a cubic subsample of EX with an applied pressure gradient in direction $`\stackrel{}{e}_1`$. The permeabilities obtained from the LB simulations are significantly higher for all $`\tau `$. This is due to the extra fluid layer with a thickness of $`10`$ % of $`M_1`$ outside the sample and the periodic boundary conditions used in the LB simulations. For this smaller piece of the sample, the effect of the boundary conditions is more significant than for the whole sample. We also performed with LB a couple of simulations with conditions similar to those used in FD. The pressure boundary condition was imitated using the body force combined with density and momentum averaging at the inlet and the outlet in the adjacent free fluid layer with thickness of one lattice spacing. No-flow boundary conditions were applied on the other cube sides. We found that the results of the LB and the FD method were again very close. This is indicated in Fig. 5 by two isolated points (filled square and star) at $`a_1=7.5\mu \mathrm{m}`$ which were obtained from LB simulations with these boundary conditions for $`\tau =0.688`$ and $`\tau =1.0`$, respectively. However, a complete recomputation of the data of Fig. 5 would have exhausted the available computation time. Moreover, the results for the free fluid layer and for periodic boundary conditions allow us to estimate the influence of the boundary conditions on the result.

For the original resolution the permeabilities obtained from the LB algorithm differ with varying $`\tau `$ by nearly a factor $`2`$. This again shows that the accuracy of the LB results for low-porosity (and low-discretization) systems strongly depends on the relaxation parameter. While for $`\tau =0.6`$ the LB results are nearly independent of the lattice spacing, the changes in the permeability are drastically increased with increasing $`\tau `$. As we have already discussed, in the permeability simulations, decrease in the relaxation parameter $`\tau `$ (or viscosity) can be used to compensate an inadequate grid resolution (finite-size effects) within certain limits. By testing the structures of the subsamples for several grid resolutions and values of $`\tau `$, we chose $`\tau =0.688`$ for the main simulations of the permeabilities. This is not necessarily the optimal choice and it differs from the results for pipe flow, but a more accurate calibration of $`\tau `$ would require further simulations or other independent results.

The permeability $`k^{(\mathrm{})}`$ for an infinite resolution of the lattice is obtained from a linear extrapolation of the data displayed in Fig. 5. Extrapolation of the FD results yields $`k^{(\mathrm{})}=224\mathrm{mD}`$. Thus, we obtain an estimate for the relative error of our permeabilities in the case of EX, $`(k^{(1)}k^{(\mathrm{})})/k^{(1)}0.33`$. For LB the extrapolated permeability is $`k^{(\mathrm{})}=550\mathrm{mD}`$. Analogously, for $`\tau =0.688`$ we obtain $`(k^{(1)}k^{(\mathrm{})})/k^{(1)}0.15`$. For DM this error of LB is about the same size as for EX, and approximately $`0.30.4`$ for SA and GF.

The simulations of the subsamples show also that the magnitude of fluctuations in permeability, within one and the same sample, can be 100 % even when the sample is extremely homogeneous. An effect of similar size may be induced by inaccurate boundary conditions. This is important when comparing two micropermeameter experiments in which different boundary conditions may have applied.

When analyzing a subsample of sample SA, we encountered another difference between the FD and LB simulation. Geometrical analysis of the considered subsample of SA reveals that this subsample is not percolating in direction $`\stackrel{}{e}_1`$. We found nevertheless that its LB permeability is $`k_{11}^{(1)}=50\mathrm{mD}`$ for $`\tau =0.688`$, and for the original resolution $`a_1=7.5\mu \mathrm{m}`$. We attribute this result to diagonal leaks that are present in the LB model used.

### V.5 Velocity distributions

Next we consider the velocity fields in more detail, and analyze the histograms of velocity magnitudes. Figure 6 shows the scaled distributions of the magnitude of the velocity $`|\stackrel{}{v}(\stackrel{}{x})|`$ with $`\stackrel{}{x}`$ for samples EX and SA. The distributions were sampled using the solutions of the flow fields for an applied pressure gradient in direction $`\stackrel{}{e}_1`$. The higher permeability of the original sandstone EX is reflected by a higher probability density of regions with average flow velocity. The distribution of sample SA on the other hand exhibits a higher peak at $`|\stackrel{}{v}|=0`$, and it extends to higher velocities. The former observation indicates large stagnant areas where no transport is taking place. The increased probability at high velocities may be related to a large number of narrow pore throats through which the fluid has to move. The velocity results for sample GF resemble closely to those of sample SA while those for sample DM resemble those for EX.

Unlike the results for the permeabilities, which are close to each other, the velocity distributions reveal more significant differences between FD and LB calculations (see Fig. 6). The distributions obtained from the LB solution exhibit a maximum at small velocities which is not present in the distributions obtained from the FD solution. The differences could perhaps be attributed to the slip velocities at the boundary due, e.g., to the bounce-back boundary condition and the diagonal leak flows in the LB algorithm, which could lead to a systematic deficit of zero velocities near the boundary.

For SA the LB simulations show an additional spurious (double) peak at $`|\stackrel{}{v}|/|\stackrel{}{v}|1.5`$. Similar, although smaller peaks were also found for GF and EX. The locations of these peaks were found to be $`\mathrm{\Delta }P/2\rho |\stackrel{}{v}|`$, where $`\mathrm{\Delta }P`$ is the body force. We attribute these peaks to the staggered momenta found in small closed pores, the number of which is very high for SA and GF. We have checked that these peaks disappear when time-averaged velocities are used, and a corresponding increase appears at zero velocity, which is the expected velocity for small closed pores.

## VI conclusion

We have performed numerical micropermeametry on three dimensional porous microstructures with a linear size of approximately 3 mm and a resolution of 7.5 $`\mu `$m. One of the samples has been a microtomographic image of Fontainebleau sandstone. Two of the samples were stochastic reconstructions with the same porosity, specific surface area, and two-point correlation function as the Fontainebleau sample. The fourth sample was a physical model, which mimics the processes of sedimentation, compaction and diagenesis of Fontainebleau sandstone. The permeabilities of these samples were determined by numerically solving at low Reynolds numbers the appropriate Stokes equations in the pore spaces of the samples, using standard finite differences methods and the lattice-Boltzmann method. Our work shows that both methods, the LB method as well as standard FD methods, are applicable to the solution of the steady-state Stokes equation within the microstructure of a three-dimensional porous medium. We investigated systems with sizes of up to $`400^3`$ lattice points. The solution of even larger systems seems possible. Hence, numerical micropermeametry is becoming a feasible technique for studying permeability fluctuations.

An accurate, quantitative comparison of the two numerical methods is difficult due to the different approaches underlying these methods. The memory requirements of the algorithms used in this study differ by roughly a factor of 2.5. The FD algorithm requires to store $`8`$ real numbers per lattice node, while the LB algorithm needs $`19`$ real numbers per node. However, in this LB model 15 real numbers could also be used Qia92 . Considering the time consumption, both methods are quite similar. Our comparison has shown that there are some features in the standard application of the LB method, which need special attention. These include the $`\tau `$-dependence of the no-flow boundary, the compressibility of the fluid, staggered invariants, and diagonal leak flows. One should notice that these occasionally inconvenient features can usually be eliminated if so needed. The compressibility of the fluid can be eliminated for stationary flows Lin96 , and the effects of staggered invariants can be eliminated with proper averaging. The diagonal leak flow, which becomes noticeable for rough surfaces and near the percolation threshold and limits there the accuracy of the method, can also be eliminated by fairly straightforward means. The $`\tau `$-dependence of the no-flow boundaries can as well be eliminated, entirely by introducing a modified LB model Gin94 , or almost entirely by using second-order boundaries as discussed above. We could of course have implemented here all these corrections in the LB code, and achieved thereby a more favourable comparison with the FD code and experiment, but we wanted to show the points of concern in a “standard” implementation of the LB method.

Our results provide a quantitative comparison of various models for porous rocks. We show for the first time that stochastic reconstruction models for Fontainebleau are less accurate than originally believed adl92 ; YT98a . In addition, our results for the permeabilities of the Fontainebleau sandstone and its models confirmed previous predictions hil99c  of a purely geometrical investigation of the same samples based on local porosity theory. The numerical value of $`k`$ for the EX sample was found to be in good agreement with the experimental value. Nevertheless, one has to keep in mind that “numerically exact” results for the fluid permeability must be handled with care. We have shown here that they may depend on the numerical method, the boundary conditions, the size of the sample, and the resolution of the microstructure. Errors of as much as 100 % cannot typically been ruled out. In summary, numerically exact determination of permeability is difficult to achieve. On the other hand, the same is true for the precision measurements in an experiment.

## Acknowledgments

U.A., A.K. and J.T. thank the Center for Scientific Computing in Finland for providing computational resources, and the Academy of Finland for financial support under the MaDaMe Programme and the Finnish Center of Excellence Programme 2000-2005 (Project No. 44875). C.M. and R.H. thank the HLRS at Universität Stuttgart for providing computing time and the Deutsche Forschungsgemeinschaft for financial support.