[1]

$$\begin{eqnarray} { \enq{\end{eqnarray}$$

r x M A C D E F G K L P R S T V W Int Bdy

Stereology as inverse problem

James G. Berryman

berryman@sep.stanford.edu

ABSTRACT

Stereology is the part of imaging science in which the three-dimensional structure of a body is determined from two-dimensional views. Although it is relatively easy to determine volume information from 2-D slices, it is nontrivial in general to determine other physical properties such as internal surface areas unless the medium is known to have some simple symmetry such as isotropy. For this reason, stereology can be viewed as a type of inverse problem. In earlier work I showed that an anisotropic spatial correlation function of a random porous medium could be used to compute the specific surface area when it is stationary as well as anisotropic by first performing a three-dimensional radial average and then taking the first derivative with respect to lag at the origin. This result generalized the earlier result for isotropic porous media of Debye et al. (1957). Here I provide more detailed information about the use of spatial correlation functions for anisotropic porous media and in particular I show that, for stationary anisotropic media, the specific surface area can be related to the derivative of the two-dimensional radial average of the correlation function measured from cross sections taken through the anisotropic medium. The main concept is first illustrated using a simple pedagogical example for an anisotropic distribution of spherical voids. Then, a general derivation of formulas relating the derivative of the planar correlation functions to surface integrals is presented. When the surface normal is uniformly distributed (as is the case for any distribution of spherical voids), my formulas can be used to relate specific surface area to easily measureable quantities from any single cross section. When the surface normal is not distributed uniformly (as would be the case for an oriented distribution of ellipsoidal voids), my results show how to obtain valid estimates of specific surface area by averaging measurements on three orthogonal cross sections. One important general observation for porous media is that the surface area from nearly flat cracks may be underestimated from measurements on orthogonal cross sections if any of the cross sections happen to lie in the plane of the cracks. This result is illustrated by taking the very small aspect ratio (penny-shaped crack) limit of an oblate spheroid, but holds for other types of flat surfaces as well.

INTRODUCTION

Seismologists are very familiar with the concepts of tomography, whether X-ray or seismic, in which details of some physical property such as density or velocity are displayed for one plane of a body after the solution of a moderately complex inversion problem. A much older part of imaging science is stereology, in which the three-dimensional structure of a body is determined from two-dimensional views. In the modern laboratory, these two imaging and analysis methods may be (and sometimes are) used in tandem; X-ray tomography may be used to obtain a 2-D image which can then be analyzed using stereological methods. While it is relatively easy to determine volume information from 2-D slices, it is nontrivial in general to determine other physical properties such as internal surface areas unless the medium (which we assume here to be random as will be the case for rocks) is known to have some simple symmetry such as isotropy or transverse isotropy. Specific surface area is a key physical parameter of a rock. This parameter is known to be related to fluid permeability of rocks and is therefore of great practical interest to the oil industry. Stereology for stationary anisotropic random media has not been explored in much detail and is an important new area of analysis for rock physics. Furthermore, like tomography, stereology is a type of inverse problem, using incomplete sets of measurements to deduce internal (statistical) structure of the material being analyzed. The solution of this inverse problem for real materials is a challenge. In this paper, I will focus on the use of spatial correlation functions to address the issue of measurement of specific surface area from cross sections of anisotropic porous rocks. The paper should be viewed as a first step towards solution of the inverse problem, as it shows what the expected range of uncertainty will be for some typical microgeometries and provides a method of determining the specific surface area with reasonable accuracy if enough orthogonal views are obtained.

When images of cross sections of a random porous material are available for analysis, the techniques of stereology (Greenwood, 1970) can be applied to determine various physical properties such as the porosity and the specific surface area. The tools used are typically counting methods based on lengths of intersections of lines with voids for porosity and numbers of intersections of lines with interfaces for specific surface area. With the advent of modern computing capabilities, the possibility of obtaining significantly more information about the statistics and characteristics of the microgeometry of random composite materials has not escaped notice. There have been extensive efforts already published in the literature providing details of various means of making such measurements, particularly through the computation and analysis of spatial correlation functions (Debye et al., 1957; Berryman, 1985; Berryman and Blair, 1986; Berryman, 1987; Berryman and Blair, 1987; Martys and Garboczi, 1992; Blair and Berryman, 1992; Adler, 1992; Coker and Torquato, 1995; Coker et al., 1996; Blair et al., 1996) for isotropic media. Although such applications are more difficult in general, correlation functions may nevertheless also be applied to anisotropic random media (Berryman, 1987; Berge et al., 1997).

In this present paper, I focus on clarifying and generalizing my previous work (Berryman, 1987) on stationary and anisotropic random porous media. An earlier result showed that an anisotropic spatial correlation function of a random porous medium could be used to compute the specific surface area when it is both stationary and anisotropic by first performing a three-dimensional radial average and then taking the first derivative with respect to lag at the origin. This result generalized the earlier result for isotropic porous media of Debye et al. (1957). This result was of limited use because of the idealization that the anisotropic spatial correlation function was known exactly. If the correlation function is measured from images of cross sections of the random medium, then some important deviations from the ideal case can arise in practice, depending on the nature of the anisotropy in the material.

The present paper provides more detailed information about the use of spatial correlation functions for anisotropic porous media and in particular shows that, for stationary anisotropic media, the specific surface area can be related to the derivative of the two-dimensional radial average of the correlation function measured from any cross-section through the anisotropic medium. The fundamental concepts being presented are illustrated first using a simple pedagogical example in the next section. Then the main result is derived. I show how the correlation functions may depend on the distribution of surface normals in cases when it is not uniformly distributed, and how this can cause the results obtained from cross sections to deviate from the ideal. A discussion of implications and limitations of the work is presented in the final section.

PEDAGOGICAL EXAMPLE: SPHERICAL VOIDS

One of the basic results of this paper may not seem intuitively obvious to the reader (the result seemed counter-intuitive to the author at first, which is why an illustrative model was sought). I will first present a simple example to illustrate this fundamental concept. The rigorous proof of the general result is presented in the next section.

Suppose we have an anisotropic porous medium containing N spherical holes, each of radius R. The distribution of hole centers is different in the three directions x,y,z with the average number density of hole centers in each direction given respectively by n_x/L, n_y/L, and n_z/L, so that in a cube of volume L³ there will be a total of N = n_xn_yn_z spherical holes. This example is slightly oversimplified because I am going to specify that the spheres are not overlapping which implies some correlation among sphere centers that I will not take into account in this model. And I also assume that the sphere radius $R \ll L$.So the model is only valid for very dilute concentrations of spheres. For this model, the porosity $\phi$ is given by

_true = n_xn_yn_zL^343R^3,   and the specific surface area s is given by

s_true =n_xn_yn_zL^34R^2.  

Now suppose we take a slice through this cube somewhere in the interior at a distance greater than the radius R from the edges so that edge effects do not play an important role. If I take this slice perpendicular to the z-axis, then the average number of holes that intersect the plane of the slice is determined by n_xn_y (2 R n_z/L). These intersecting spheres all produce intersections that are circles in the plane of the slice, but the radius r(h) of these intersecting circles can lie anywhere in the range $0 \le r(h) \le R$ depending on the distance h of the sphere center from the plane of intersection. The circle of intersection has radius $r(h) = \sqrt{R^2-h^2}$. The values of h are uncorrelated with the location of the plane and therefore may be treated as evenly distributed, so the average circular circumference of intersection is computed by considering

1R_0^R dh 2R^2-h^2 = 2R _0^/2 d ^2 = ^22R.   Combining this result with the probabilities of occurence quoted in the second paragraph of this section and dividing by the total area of the slice L², I find

s_planar = n_xn_yn_zL^3^2 R^2 = 4s_true.   This result differs from the true specific surface area (s) by a factor of $\pi/4$, but this difference is a standard result of stereology [see Greenwood (1970)] and therefore expected. In contrast, the average circular area of intersection is determined by

1R_0^R dh (R^2 - h^2) = 23R^2.   Multiplying by the number of sphere centers within R of the plane determines the apparent porosity $\phi_{planar}$ as

_planar = n_xn_yn_zL^3 43R^3 = _true,   in complete agreement with (porosity) which is also expected from stereology.

Now, if we consider taking the slices in other directions, we quickly discover that the results are the same as (splanar) and (Gvplanar) in each case. For example, if we take a slice perpendicular to the y-axis, then the probablities are $n_xn_z\times 2 R n_y/L$, which is simply a permutation of the first case considered (but nevertheless giving the same result). And similarly, the result does not change for a slice perpendicular to the x-axis.

For slices at arbitrary angles to the x-, y- or z-axes, the only major modification comes from the fact that the total area of the resulting slice will generally be larger than L² and the total number of spherical hole centers within a distance R on either side will also be larger than n_xn_yn_z2R/L by the same factor. Otherwise, the calculations are identical for both $\phi_{planar}$ and s_planar.

Thus, I find that, for this rather simple model, the specific surface area computed by considering the circumferential intersections of the spheres with arbitrary planes through the cube will always give the same result for any slice even though the medium is highly anisotropic. I show how this result is related to the general result in the next section, where I find that the key simplifying feature of the spherical model is its uniform distribution of surface normals.

SURFACE AREA FROM CORRELATION FUNCTIONS ON PLANE SLICES

The characteristic function f of a porous material is defined as $f(\vecx) = 1$ when the position vector lies inside the pore space and $f(\vecx) = 0$ when lies outside the pore space. The two-point spatial correlation function is then defined by

S_2(_1,_2) = <f(+_1)f(+_2)>,   where the brackets denote a volume average over the spatial coordinate . When $\vecr_1 = \vecr_2$, it is easy to see that the value of the two-point function is just equal to the volume fraction of void space in the total space. This volume fraction is called the porosity $\phi$. The definition (S2) is meaningful when the random medium is stationary in space (Beran, 1968), so that the result obtained by averaging is actually independent of the choice of the origin in the three-dimensional space occupied by the material. Then the function cannot depend independently on the vector arguments $\vecr_1$ and $\vecr_2$ but only on their difference $\vecr_2-\vecr_1$, as is easily seen by moving the origin from $\vecx = 0$ to $\vecx = - \vecr_1$.

It is therefore adequate to write

S_2() 1V_V d^3x f()f(+)   for a random, spatially stationary medium. If the medium is also isotropic, then S₂ depends only on the magnitude $r=\vert\vecr\vert$ of its argument. However, I want to study anisotropic stationary media for which (S2vol) is the most useful statement of the definition of S₂. It will also prove helpful to express the argument explicitly as $\vecr = r\hat{r}(\theta,\varphi)$, where $\hat{r}$ is the unit radial vector with $\theta$ and $\varphi$ being the usual polar and azimuthal angles in spherical coordinates.

Planar analysis

Now suppose that we take a slice through the porous medium in a plane perpendicular to the z-axis. Computing S₂ in this plane provides one realization of $S_2(r\hat{r}({{\pi}\over{2}},\varphi))$.If the slice is statistically representative of the whole, then I may consider taking the planar radial average of S₂ and the result should closely approximate

A_xy(r) = 12_0^2 d S_2(rr)                      = 12V_0^2 d_Vd^3x f() f(+ rr(2,)).   Following Debye et al. (1957) and Berryman (1987), I now take the r derivative of this expression and consider the limit as $r \to 0$:

_r0 dA_xy(r)dr = _r0 12V _0^2 d r(2,)_V d^3x f() f(+ rr) = _r0 12V _0^2 d r(2, )_V_p d^3x f(+ rr) = 1Vda_s 12_0^2 d r(2,)n_s f(_s + 0^+r).   The result has been converted to a surface integral on the interface, where $\vecx_s$ is the surface point and da_s is the infinitesimal of local surface area. The unit normal at the surface point $\vecx_s$ is $\hat{n}_s$, with this vector pointing outward from the void space. V_p appearing in the intermediate result is the pore volume.

Two-dimensions: perpendicular to axis of symmetry

For completeness, I now consider the possibility that the anisotropic stationary medium might actually be two-dimensional, i.e., the random variations might occur (for example) only in the xy-plane and have no dependence on z. In this particular case, the evaluation of (limit) is especially simple because the normal vector $\hat{n}_s$ must necessarily lie in the same plane as $\hat{r}$.Then, each point on the interface will contribute exactly the same amount to the integral since it is easy to see that

_0^2 d r(2,)n_s f(_s + 0^+r) = -_0^ d= -2,   so

A_xy^(0) = - s.   Thus, if the cross section is taken perpendicular to the axis of symmetry for a two-dimensional medium, the formula $s = -\pi A_{xy}'(0)$ maybe used to compute the true specific surface area from the planar image.

General analysis for three dimensions

Now, if the medium is truly three-dimensional as well as anisotropic, then a complication arises because the normal vector $\hat{n}_s$ at the interface does not necessarily lie in the plane of $\hat{r}$.Thus, unlike (2Dcase), the contribution to the total integral in (limit) from the interface point $\vecx_s$ will differ from position to position. In general it is clear that the two-dimensional case gives the largest magnitude contribution and therefore that it must always be true that

- A_xy^(0) s,   for three-dimensional anisotropy. If I express the normal vector $\hat{n}_s$ as a radial vector with polar and azimuthal angles $\theta_s$ and $\varphi_s$ at position $\vecx_s$, then

n_s = _s(x_s + y_s) + z_s,   and

r(2,)n_s = _srp,   where $\hat{p} = \hat{x}\cos\varphi_s + \hat{y}\sin\varphi_s$ is a unit vector in the plane of $\hat{r}(\pi/2,\varphi)$.Thus, the analysis reduces at this point to the analysis given for the two-dimensional case except for a factor of $\sin\theta_s$ and the result that follows immediately is

A_xy^(0) = - 1Vda_s_s.   If I assume that these polar angles $\theta_s$ are evenly distributed over the full $4\pi$ range of solid angles, then I have the average contribution

<_s> = 14d_s d_s ^2_s = 4.   Combining (Axyprime3D) and (sinsqave), I have

A_xy^(0) = - s4,   recovering the result (splanar) of the preceding section. Also, note that (Axyprime3D) can be used to recover the result (Axyprime2D) for the two-dimensional case by setting $\theta_s = \pi/2$ and therefore $\sin\theta_s = 1$, as is necessarily true for 2-D.

Next I consider a cross section taken in the zx-plane. The analysis goes through as before, but the pertinent radial vector is given by $\hat{r}(\theta, 0)$ and the part of the normal vector $\hat{n}_s$ that contributes to the result is proportional to the unit vector

p = (x_s_s + z_s)/ ^2_s + ^2_s^2_s   and its magnitude is most conveniently expressed as $\sqrt{1-\sin^2\theta_s\sin^2\varphi_s}$. The result analogous to (Axyprime3D) is therefore

A_zx^(0) = - 1Vda_s 1-^2_s^2_s.   To check the result for evenly distributed normal vectors, I compute

<1-^2_s^2_s> = 14d_s d_s_s 1-^2_s^2_s,   which can be evaluated using tabulated integrals (Gradshteyn and Ryzhik, 1965). Performing the $\theta_s$ integration first gives

_0^d_s_s1-^2_s^2_s = 1 + (1+_s1-_s) ^2_s2_s.   Then, using the fact that

_0^1 dxx1-x^2(1+x1-x) = ^22 - ,   I obtain the final result that

<1-^2_s^2_s> = 4,   as expected since the result for an even distribution of normal vectors should not depend on an arbitrary choice of axis labels. A similar derivation for the yz-plane shows that

A_yz^(0) = - 1Vda_s 1-^2_s^2_s,   and that

<1-^2_s^2_s> = 4,   also as expected.

Two dimensions: parallel to axis of symmetry

To illustrate the use of these formulas for cases in which the anisotropy is so special that the 3-D analysis reduces exactly to 2-D, consider again that for such cases $\theta_s = \pi/2$. Then, $\sqrt{1-\sin^2\theta_s\sin^2\varphi_s} = \vert\cos\varphi_s\vert$ and $\sqrt{1-\sin^2\theta_s\cos^2\varphi_s} = \vert\sin\varphi_s\vert$.It follows from (Azxprime3D) and (Ayzprime3D) that

A_zx^(0) = A_yz^(0) = - 2s^2,   since

12d_s |_s| = 12d_s |_s| = 2.  

It is easy to check that this result is correct for the case of a distribution of nonoverlapping circular cylinders of radius R aligned with the z-axis and having center point densities n_x/L, n_y/L. The true specific surface area is $s = n_xn_y2\pi R/L^2$, while the specific surface area that would be measured in either the zx-plane or the yz-plane is $s_{apparent} = n_xn_y(2R/L)(2/L) = 2s/\pi$. Thus, I confirm that the factor of $2/\pi$ found in (sincosave) is generic for such two-dimensional media.

NONUNIFORM DISTRIBUTION OF SURFACE NORMALS

It is easy to see that the results in my pedagogical example cannot be true for general stationary anisotropic media. Suppose that the porosity is composed entirely of oblate spheroidal holes aligned so that their minor axes are always parallel to the z-axis. If these spheroids have major axes a and minor axes c, then I can consider the situation when the minor axis satisfies $c \ll a$.This limit is often called the limit of ``penny-shaped cracks'' (Walsh, 1969). In this limit, the surface of the spheroid rigorously approaches (see the formulas in the next subsection) $2\pi a^2$, which is just the sum of the area of two circles each of radius a. So in the limit $c/a \to 0$,cross sections parallel to the z-axis will sample this surface area, while any cross sections taken perpendicular to the z-axis will miss this surface area (being a set of measure zero in this direction). Thus, this counterexample shows that it cannot be true that the specific surface area is always obtainable from any arbitrary single cross section when the distribution of surface normals is not uniform, and especially so when the distribution is strongly skewed. I therefore need to analyze this situation to understand how the surface area can be correctly obtained from sets of two or three orthogonal cross sections.

The first subsection that follows describes the idealized results that are expected when the spatial correlation function can be measured without resorting to measurements on cross sections, while the second subsection shows how the results differ when measurements on cross sections are used.

Oriented oblate and prolate spheroids

For an oblate spheroid with major axes a and minor axis c, it is straightforward to show that the volume is given by $4\pi a^2c/3$and the surface area by

A_obl = 2a^2 _0^d 1-k^2^2 = 2a^2 + c^2k1+k1-k,   where the eccentricity k is determined by

k^2 = 1 - c^2a^2.   For a prolate spheroid with major axis a and minor axes c, the volume is similarly given by $4\pi ac^2/3$ and the surface area by

A_prol = 2ac _0^d1-k^2^2 = 2c^2 + 2ack^-1k,   where k is again defined by (eccentricity).

Since the integrals in both (Sobl) and (Sprol) are just exactly $\int da_s$ in the sense of (Axyprime3D), with $\theta= \theta_s$,the pertinent integrals for the idealized derivatives of the correlation functions in the xy-plane are given by

_0^/2 d ^21-k^2^2 = 4F(32, -12; 2; k^2)   for oblate spheroids and

_0^/2 d^21-k^2^2 = _0^/2 d1-k^2^2 - _0^/2 d^21-k^2^2 = E(k) - 4F(32, -12; 2; k^2) =4F(12, -12; 2; k^2)                for prolate spheroids. The special functions appearing in (Soblxy) and (Sprolxy) are Gauss' hypergeometric function $F(\alpha, \beta; \gamma; k^2)$ and the complete elliptic integral of the second kind E(k). The result on the right of (Soblxy) can be obtained by expanding the square root on the left and integrating term by term, or it may be found in standard references (Gradshteyn and Ryzhik, 1965; Whittaker and Watson, 1965; Abramowitz and Stegun, 1970). The elliptic integral E(k) may also be written in terms of the hypergeometric function as

E(k) = 2F(12, -12; 1; k^2).   The hypergeometric function is relatively easy to evaluate numerically as a power series in k². In the numerical examples that follow I computed these integrals using 200 terms in double precision. Note that the final identity on the far right of (Sprolxy) follows by transforming the integral on the far left to $\int d\theta\cos^2\theta\sqrt{1 - k^2\sin^2\theta}$,using the properties of the trigonometric functions.

For slices in the zx-plane, the pertinent integrals for oblate spheroids are

dd1 - ^2^2 1 - k^2^2                      = 8 _0^/2 dE() 1 - k^2^2,   while for prolate spheroids they are

dd1 - ^2^2 1 - k^2^2                      = 8 _0^/2 dE() 1 - k^2^2.   These integrals are substantially more difficult to analyze than any of the other formulas that I have included in the main text so far or that I include in the following discussions, so - in order to avoid a major digression - I present my method of evaluation of (Soblyz) and (Sprolyz) in the appendix.

Figures 1 and 2 compare the results (Sobl) with (Soblxy), and (Sprol) with (Sprolxy), respectively. Because of the standard stereological correction of $\pi/4$, I have divided (Soblxy), (Sprolxy), (Soblyz), and (Sprolyz) by this factor before plotting. Observe that the results for cross sections of oblate spheroids diverge from the exact result (Sobl) much more so than the results for prolate spheroids diverge from (Sprol). If in addition I were to plot the simple averages for the results for the three orthogonal slices for comparison on the plots for oblate and prolate surface areas, I would find that the averages overlay the solid curves and would be indistinguishable from them. This result shows that the earlier result (Berryman, 1987) correctly predicts that for ideal data from perfect spatial correlation functions the full $4\pi$ radian average of the anisotropic S₂ can simply be replaced by the average of results from three mutually orthogonal planes. Furthermore, information from two orthogonal planes may be sufficient if it is known that one of the planes is orthogonal to an axis of symmetry while the other is parallel to it.

 

fig1
Figure 1 Normalized surface area A_obl/a² in (29) for oriented oblate spheroids (solid line), compared to the apparent surface area (that would be obtained from the exact correlation function S₂) in the xy-plane (dot-dash line) computed using (32) and to the apparent surface area in either the zx- or the yz-plane (dashed line) using (44). The eccentricity is given by $k = \sqrt{1-c^2/a^2}$.The apparent surface areas for all slices have been divided by the standard stereological correction factor $\pi/4$.When the average of the results from these three types of mutually orthogonal (idealized) slices are averaged, the result is indistinguishable from the solid line.

 

fig2
Figure 2 Normalized surface area A_prol/a² in (31) for oriented prolate spheroids (solid line), compared to the apparent surface area (that would be obtained from the exact correlation function S₂) in the xy-plane (dot-dash line) computed using (33) and to the apparent surface area in either the zx- or the yz-plane (dashed line) using (46). The apparent surface areas for all slices have been divided by the standard stereological correction factor $\pi/4$.When the average of the results from these three types of mutually orthogonal (idealized) slices are averaged, the result is indistinguishable from the solid line.

Planar cross sections and observability

The preceding calculations summarize the results for the ideal case when S₂ is known precisely. When S₂ is obtained instead from measurements of individual realizations of cross sections, I can expect there to be some deviations from the ideal results. To illustrate this situation, consider that when a cross section is made, an image of the surface shows the changes from solid to void across this surface; but note further that these observed changes actually measure interfaces that intersect the plane of the surface at some finite (nonzero) angle. If there happens to be an interface that lies exactly (i.e., within the tolerances of the sectioning procedure) in the plane of the imaged surface, its presence will not be observed. For example, consider again (as in the opening paragraph of this section) the oblate spheroid in the limit when $c/a \to 0$. In this limit, the oblate spheroid degenerates into a penny-shaped crack with surface area $2\pi a^2$. If the resulting crack lies exactly (practically speaking) in the xy-plane, then cross sections taken perpendicular to the z-axis will never observe this interface area, while cross sections taken parallel to the z-axis will intersect it at a finite angle and thereby allow observations to be made. The mathematical reason for this is that the distribution function has degenerated into a delta function along z in this limit; while the delta function is correctly included in the ideal S₂, it will sometimes be missed when making cross-sectional measurements.

The results that are presented next are those that would be obtained by computing the zero lag derivative of the spatial correlation function of an individual realization of a representative cross section through these anisotropic media.

First consider the same distribution of oblate spheroids, oriented so their minor axis is aligned with the z axis. Then, it is easy to see that any cross section perpendicular to the z-axis that intersects one of these spheroids will produce a circular intersection of radius $r(h) = a\sqrt{1-h^/c^2}$, so the average area is given by

1c_0^c dh a^2(1 - h^2/c^2) = 2a^23,   while the average circumference is

1c_0^c dh 2a1-h^2/c^2 = ^2 a2.   As in the case of the pedagogical example for spherical voids, the resulting measured porosity and apparent specific surface area can be obtained by multiplying these values by n_xn_yn_z(2c/L)/L². Doing so, I find that the porosity is given exactly (as expected) and the surface area differs both from (Sobl) and from (Soblxy).

For cross sections taken parallel to the z-axis, the average area is given by

1a_0^a dh a c (1-h^2/a^2) = 2a c3.   The average circumference is obtained by first noting that the circumference of intersection is determined multipying 4 E(k) by the major ellipsoidal axis in the plane which will be $a\sqrt{1-h^2/a^2}$. Thus, the average circumference is given by

1a_0^a dh 4aE(k)1-h^2/a^2 = a E(k).   The measured porosity and specific surface area can be obtained in this case by multiplying these values by n_xn_yn_z(2a/L)/L². The porosity is again given correctly by the result (oblateyzvol), while the surface area obtained from (oblateyzsurf) differs. These results for specific surface area are shown in comparison to the exact results in Figure 3. I divided (oblatexysurf) and (oblateyzsurf) by $\pi/4$ to account for the standard stereological correction for the surface area of spheres, or for uniform distribution of surface normals. This correction is handled automatically in the correlation function analysis by setting the specific surface area equal to $-4\times S_2^\prime(0)$.Note that this Figure shows that (oblatexysurf) and (oblateyzsurf) taken separately are both rather poor approximations to the correct result (Sobl), even for relatively small deviations from sphericity (e.g., notice that $c/a \le 0.99$ corresponds to $k \ge 0.14$).

Results for a distribution of prolate spheroids, all oriented with major axis along z, can be simply obtained by interchanging the roles of a and c in the three formulas (oblatexyvol)-(oblateyzvol). A little more care is required in translating from (oblateyzsurf) to the result for prolate spheroids sliced in the yz-plane, because the result (oblateyzsurf) becomes

1c_0^c dh 4aE(k)1-h^2/c^2 = a E(k),   which has the same right hand side as (oblateyzsurf). The reason for this is just that the circumference of the ellipse is always defined in terms of E(k) and the major axis which is proportional to a in both cases. This result is then multiplied by n_xn_yn_z(2c/L)/L² and this final result is not what I would have obtained by simply interchanging a and c in the corresponding result for oblate spheroids. These results for specific surface area are shown in comparison to the exact results in Figure 4. Again I divided the formulas for apparent specific surface area of oriented prolate spheroids by $\pi/4$ to account for the standard stereological correction. Note that, in contrast to the results for oblate spheroids, the individual cross-sectional results for prolate spheroids do not deviate drastically from the correct result and could be used to provide estimates accurate within about 5% or better over the entire range of k values.

 

fig3
Figure 3 Comparison of the apparent surface areas that would be obtained on average from individual realizations of S₂ from cross sections of random medium composed of a distribution of oriented oblate spheroids. Solid line is the exact result (29). Dot-dash line is the apparent surface area obtained for an xy-slice according to (38) times the appropriate weighting factor for the oblate case. Dashed line is the apparent surface area for either zx- or yz-slices according to (40). Apparent surface areas have been divided by $\pi/4$.

 

fig4
Figure 4 Comparison of the apparent surface areas that would be obtained on average from individual realizations of S₂ from cross sections of random medium composed of a distribution of oriented prolate spheroids. Solid line is the exact result (31). Dot-dash line is the apparent surface area obtained for an xy-slice according to (38) times the appropriate weighting factor for the prolate case. Dashed line is the apparent surface area for either zx- or yz-slices according to the discussion following (41). Apparent surface areas have been divided by $\pi/4$.

The final question is how to combine the cross-sectional results to obtain proper estimates of the specific surface area for anisotropic media. The obvious combination to try is the simple average of the results obtained in three orthogonal directions. It is reasonable to suppose this might work because of the earlier results on the ideal case [see Berryman (1987)] and because the averaging did produce the correct results in Figures 1 and 2. Since the present results in Figures 3 and 4 show clearly that the measurement bias is distributed on both sides of the correct result (as was the case in Figs.1 and 2), the average might be expected to give good results. I show the computed averages in Figures 5 and 6. We see in Figure 6 that the average obtained from three orthogonal cross sections gives excellent results for the case of prolate spheroids for all values of the eccentricity. The average for oblate spheroids in Figure 5 is not quite as good as that for prolate spheroids for $k \ge 0.6$, but I have already pointed out why there should in fact be deviations for this case as the eccentricity goes to zero. In any case, the deviations even for the oblate spheroids are small and the largest deviations occur for cases where the overall surface area is relatively small. Thus, a large deviation in a small quantity still results in a small deviation in the final result if I assume that real materials will in fact often contain a complex mixture of the results presented here.

 

fig5
Figure 5 Comparison of the exact surface area for oblate spheroids (29) with the average of the results displayed in Fig.3. These results differ as the eccentricity k rises due to the effect discussed in the text that flat cracks will be undersampled in real cross-sectional imaging, leading to an underestimate of the apparent surface area.

 

fig6
Figure 6 Comparison of the exact surface area for prolate spheroids (31) with the average of the results displayed in Fig.4. These results are identical within the precision of our calculations.

Of course, voids in real materials are not restricted to the spheroids considered here, but the oblate and prolate spheroids at least qualitatively span the range of the types of behavior expected. Therefore, I anticipate that real materials will show behavior very similar to that covered in my analysis.

DISCUSSION

The main results of the paper are summarized in the three formulas (Axyprime3D), (Azxprime3D), and (Ayzprime3D) for ideal measurements of the spatial correlation function S₂. They apply to an arbitrary stationary random three-dimensional anisotropic porous medium. Indeed, it has been shown that the formulas apply equally well to three-dimensional cross sections of two-dimensional materials [see (Axyprime2D) and (Azxprime2D)], since such materials are also stationary in three dimensions.

One important general conclusion that is easily derived from my formulas is that the specific surface area can be correctly computed for anisotropic, stationary random media from any single representative cross section as long as the distribution of surface normals is uniform. However, if the distribution of surface normals is not uniform, then care must be taken to compute the apparent s in three orthogonal planes and then take the average.

The examples show that these results change slightly if the correlation functions are obtained from individual realizations of cross sections. Flat cracks will generally be missed in the analysis of such cross sections and therefore the measured specific surface area may be smaller than the true specific surface area. This issue has been illustrated here for oblate spheroids, but the same issue is equally important for other types of flat, or nearly flat (within the precision of the sectioning method), internal surfaces.

Finally, it is also important to recognize that not all anisotropic random systems of practical interest are stationary. One particularly important example that arises commonly when studying earth materials is a layered medium with transverse isotropy. If the z-axis is taken as the vertical (into the earth), then it is very common to assume that layers (planes) perpendicular to the z-axis (of symmetry) are isotropic and stationary in the two orthogonal directions (xy-plane). However, stationarity may not be valid for such earth materials because it is often the case that the behavior along the z-axis has some trend with depth and is therefore nonstationary. Thus, the results of the present paper do not apply to such materials. However, it would be incorrect to conclude that spatial correlation functions cannot be used for such transversely isotropic media. A different set of results applies to such materials and these have been discussed recently by Berge et al. (1997).

ACKNOWLEDGMENTS

I thank P. A. Berge and S. C. Blair for helpful discussions.

REFERENCES

• Abramowitz, M., and Stegun, I. A., (eds.), 1970, Handbook of Mathematical Functions, Dover, New York, pp.556-562, 590-592, 608-609.
• Adler, P. M., 1992, Porous Media: Geometry and Transports, Butterworth-Heinemann, Boston, Massachusetts, pp.46, 503-504.
• Beran, M. J., 1968, Statistical Continuum Theories, Interscience, New York, p.39.
• Berge, P. A., Berryman, J. G., Blair, S. C., and Peña, C., 1997, Scalar properties of transversely isotropic tuff from images of orthogonal cross sections: to appear in Proceedings of the Second International Conference on Imaging Technologies: Techniques and Applications in Civil Engineering, edited by S. McNeil, D. Frost, and L. Vuillet Engineering Foundation, New York, Davos, Switzerland, May 25-30, 1997.
• Berryman, J. G., 1985, Measurement of spatial correlation functions using image processing techniques: J. Appl.Phys., 57, 2374-2384.
• Berryman, J. G., 1987, Relationship between specific surface area and spatial correlation functions for anisotropic porous media: J. Math.Phys., 28, 244-245.
• Berryman, J. G., and Blair, S. C., 1986, Use of digital image analysis to estimate fluid permeability of porous materials: Application of two-point correlation functions: J. Appl.Phys., 60, 1930-1938.
• Berryman, J. G., and Blair, S. C., 1987, Kozeny-Carman relations and image processing methods for estimating Darcy's constant: J. Appl.Phys., 62, 2221-2228.
• Blair, S. C., Berge, P. A., and Berryman, J. G., 1996, Using two-point correlation functions to characterize microgeometry and estimate permeabilities of sandstones and porous glass: J. Geophys.Res., 101, 20359-20375.
• Blair, S. C., and Berryman, J. G., 1992, Permeability and relative permeability in rocks: in Fault Mechanics and Transport Properties of Rocks, edited by B. Evans and T.-f. Wong, Academic Press, San Diego, Chapter 7, pp.169-186.
• Coker, D. A., and Torquato, S., 1995, Extraction of morphological quantities from a digitized medium: J. Appl.Phys., 77, 6087-6099 .
• Coker, D. A., S. Torquato, and J. H. Dunsmuir, 1996, Morphology and physical properties of Fontainebleau sandstone via a tomographic analysis: J. Geophys.Res., 101, 17497-17506.
• Debye, P., Anderson, H. R., Jr., and Brumberger, H., 1957, Scattering by an inhomogeneous solid. II. The correlation function and its application: J. Appl.Phys., 28, 679 .
• Gradshteyn, I. S., and Ryzhik, I. W., 1965, Table of Integrals, Series, and Products, Academic, New York, pp.386, 557, 904-909, 1039-1045.
• Greenwood, E. E., 1970, Quantitative Stereology, Addison-Wesley, Reading, Massachusetts, pp.23-70.
• Martys, N., and Garboczi, E. J., 1992, Length scales relating the fluid permeability and electrical conductivity in random two-dimensional model porous media: Phys.Rev.B, 46, 6080-6090.
• Walsh, J. B., 1969, New analysis of attenuation in partially melted rock: J. Geophys.Res., 74, 4333-4337.
• Whittaker, E. T., and Watson, G. N., 1965, A Course of Modern Analysis, Cambridge University Press, Cambridge, pp.281-301, 514-518.

A

EVALUATION OF INTEGRALS

In this appendix, I will evaluate the integrals required to compute the ideal apparent surface area that would be observed in the zx-plane for both oblate and prolate spheroidal voids. The formulas to be evaluated are (Soblyz) and (Sprolyz).

To proceed in evaluating (Soblyz), I first note the identity (Gradshteyn and Ryzhik, 1965)

_0^/2 d^2n+11 - k^2^2 = 2^n n!(2n+1)!!F(n+1,-12; n + 32; k^2),   which can be derived by expanding the square root as a power series in k² and then evaluating the resulting integrals of powers of $\sin\theta$term by term. The special function appearing in (sinidentity) is the hypergeometric function $F(\alpha, \beta; \gamma; k^2)$. To make use of (sinidentity), I must first expand $E(\sin\theta)$ according to

E() = 2F(12, -12; 1; ^2) = 2{1 - _n=1^ [(2n-1)!!2^nn!]^2^2n2n-1}.   Substituting (Eexpansion) into the integral on the right hand side of (Soblyz), I find that

_0^/2 dE() 1 - k^2^2 = 2 {F(1,-12; 32; k^2)      - _n=1^ [(2n-1)!!2^nn!(4n^2-1)] F(n+1,-12; n + 32; k^2)}.   Equation (Soblintyz) shows that the desired integral is expressible as an infinite series, each of whose terms contains a hypergeometric function as a factor.

To evaluate (Sprolyz), I note the identity (Gradshteyn and Ryzhik, 1965)

_0^/2 d^2n+11 - k^2^2 = _0^/2 d^2n+11 - k^2^2 = 2^n n!(2n+1)!!F(12, -12; n + 32; k^2).   Using (Eexpansion) as before, I find easily that

_0^/2 dE() 1 - k^2^2 = 2 {F(12,-12; 32; k^2)      - _n=1^ [(2n-1)!!2^nn!(4n^2-1)] F(12,-12; n + 32; k^2)}.  

I have found these expressions can be computed accurately by using about 40 terms of the series (Soblintyz) and (Sprolintyz), while taking about 200 terms in the expansion of the hypergeometric functions in double precision.