Click
Here

WATER RESOURCES RESEARCH, VOL. 44, W12435, doi:10.1029/2008WR006853, 2008

for

Full
ArticleDiffusion ellipsoids of anisotropic porous rocks calculated by X-ray
computed tomography-based random walk simulations
Yoshito Nakashima,1 Susumu Kamiya,1 and Tsukasa Nakano1
Received 21 January 2008; revised 31 July 2008; accepted 1 October 2008; published 23 December 2008.

[1] Water molecules and contaminants migrate in water-saturated porous strata by

diffusion in systems with small Péclet numbers. Natural porous rocks possess theanisotropy for diffusive transport along the percolated pore space. An X-ray computed
tomography (CT) based approach is presented to quickly characterize anisotropic diffusion
in porous rocks. High-resolution three-dimensional (3-D) pore images were obtained for a
pumice and three sandstones by microfocus X-ray CT and synchrotron microtomographysystems. The cluster-labeling process was applied to each image set to extract the 3-D
image of a single percolated pore cluster through which diffusing species can migrate along distance. The nonsorbing lattice random walk simulation was performed on the
percolated pore cluster to obtain the mean square displacement. The self-diffusion
coefficient along each direction in the 3-D space was calculated by taking the time
derivative of the mean square displacement projected on the corresponding direction. A
diffusion ellipsoid (i.e., polar representation of the direction-dependent normalized
self-diffusivity) with three orthogonal principal axes was obtained for each rock sample.
The 3-D two-point autocorrelation was also calculated for the percolated pore cluster of
each rock sample to estimate the pore diameter anisotropy. The autocorrelation ellipsoids
obtained by the ellipsoid fitting to the high correlation zone were prolate or oblate inshape, presumably depending on the eruption-induced deformation of magma and regional
stress during sandstone diagenesis. The pore network anisotropy was estimated by
calculating the diffusion ellipsoid for uniaxially elongated or compressed rock images.
The degree and direction of the geological deformation of the samples estimated by thepore diameter anisotropy analysis agreed well with those estimated by the pore
network anisotropy analysis. We found that the direction of the geological deformation
coincided with the direction of the major (or minor) principal axis of the prolate (or oblate)diffusion ellipsoid for each sample. Thus, it can be concluded that the deformationinduced pore structure anisotropy is responsible for the anisotropy of the diffusivetransport properties.
Citation: Nakashima, Y., S. Kamiya, and T. Nakano (2008), Diffusion ellipsoids of anisotropic porous rocks calculated by X-ray
computed tomography-based random walk simulations, Water Resour. Res., 44, W12435, doi:10.1029/2008WR006853.

1. Introduction[2] Diffusion or random walk is a dominant mechanismof material transport in water-saturated porous strata systems with small Péclet numbers. Examples include the porescale water/rock interactions [e.g., Sorai et al., 2007]
and meter-scale contaminant migration [e.g., Ohlsson and
Neretnieks, 1995] in deep (1 km) strata with small
hydraulic head difference relevant to the undergroundCO2 storage [e.g., Gherardi et al., 2007] and high-level
radioactive waste disposal [e.g., Madsen, 1998]. Because
contaminants and water molecules diffuse along the microscopic pore space, pore-scale characterization is needed tounderstand the diffusive migration mechanism in porous
rocks. Natural porous rocks often possess anisotropic dif1
National Institute of Advanced Industrial Science and Technology,
Geological Survey of Japan, Tsukuba, Japan.

Copyright 2008 by the American Geophysical Union.
0043-1397/08/2008WR006853$09.00

fusive transport properties [e.g., Farver and Yund, 1999].
Characterization of the pore-scale anisotropy is thus essential to understanding the diffusive transport properties ofreal porous rocks. Because conventional laboratory methods
using diffusion cells or electron microscopy [e.g., Van Loonet al., 2004; Yokoyama and Nakashima, 2005] are timeconsuming and exhausting, a new quick method should be
developed to examine the direction-dependent diffusivity ofanisotropic rocks. In relation to the diffusive anisotropy,
there is a likely hypothesis that the pore structure anisotropy
causes the anisotropic diffusive transport [Yokoyama andNakashima, 2005]. For example, geometrically tortuous
pipes in rocks lead to high diffusion tortuosity (i.e., low
self-diffusion coefficient), as illustrated by Figure 1. Thishypothesis needs to be tested to understand the cause of the
diffusive anisotropy in rocks. The diffusion ellipsoid (i.e.,polar representation of the direction-dependent normalized
diffusivity) is usually convex ellipsoidal [e.g., Westin et al.,2002]. Recently, however, a few studies on the water selfdiffusion for synthetic clay gels using nuclear magnetic

W12435

1 of 19

W12435

NAKASHIMA ET AL.: DIFFUSION ELLIPSOIDS OF ANISOTROPIC POROUS ROCKS

W12435pore structure of the geologically deformed samples wasanalyzed in terms of the pore diameter anisotropy and pore
network anisotropy. The pore diameter anisotropy wasestimated by computing the two-point autocorrelation ellipsoid of the pore cluster [Ikeda et al., 2000]. The pore
network anisotropy was estimated by calculating the diffusion ellipsoid for uniaxially elongated or compressed rock
images. To test the relationship between pore structure
anisotropy and diffusive anisotropy, the deformation direction inferred from the analysis of pore diameter anisotropy
and pore network anisotropy was compared with the direction of the principal axes of the diffusion ellipsoids. As forthe unusual constriction of the diffusion ellipsoid, we firstfound that it occurs for natural porous rocks in the present
study. To discuss the generation mechanism of this constriction, diffusion ellipsoids for highly anisotropic porous
media were calculated using synthetic pore CT images
produced by systematic uniaxial elongation/compression
using affine transformation. The effect of the degree ofthe elongation/compression on the ellipsoid shape was
examined.
Figure 1. Two-dimensional schematic of the pore anisotropy. The vertical pipes are more tortuous, yielding a low
self-diffusion coefficient for pore fluid molecules comparedwith the horizontal, straight pipes. This is the pore network
anisotropy. The system involves another type of anisotropy(i.e., pore diameter anisotropy); while the vertical pipes are
thick, the horizontal ones are thin in diameter. Both types of
the anisotropy can be caused by the uniaxial compression in
the vertical direction. An isolated pore, which cannot
contribute to the long-distance material transport, is also
shown.
resonance spectroscopy [Porion et al., 2001, 2003; deAzevedo et al., 2007] reported the unusual shape, namely
constricted in the direction of minimum self-diffusivity (nolonger convex ellipsoidal). This breakdown of the convex
ellipsoid expression likely occurs also for natural porous
rocks. Because the mechanism of the generation of thispinched-in ellipsoid has not been studied in detail, this is
another point to be addressed in relation to the diffusion
anisotropy study on rocks.
[3] X-ray computed tomography (CT) is a powerful tool
for obtaining accurate three-dimensional (3-D) images of
microscopic pores rapidly [e.g., Auzerais et al., 1996;
Lindquist et al., 2000; Karacan and Mitchell, 2003; Arnset al., 2004, 2005 Ketcham and Iturrino, 2005; Al-Raoushand Willson, 2005; Altman et al., 2005; Liu et al., 2005; Duet al., 2007; Goldstein et al., 2007; Jones et al., 2007], andisotropic diffusion has been simulated by random walk indigitized CT images of porous sandstones [e.g., Schwartz etal., 1994; Nakashima et al., 2004]. In the present study, we
extended this CT-based method to the quick evaluation ofanisotropic diffusion. Digital 3-D pore images wereobtained by X-ray CT for one pumice and three sandstone
samples. Nonsorbing lattice random walk simulations were
performed in a percolated pore cluster of the CT images tocalculate the mean square displacement along an arbitrary
direction in 3-D space. The diffusion ellipsoid was
calculated by taking the time derivative of the directiondependent mean square displacement. The anisotropy of the

2. Calculation of the Diffusion Ellipsoidby Random Walk in CT Images[4] We developed a Mathematica1 version 6 program,
DellipsoidM6.nb, to calculate a diffusion ellipsoid using a
3-D pore image set. The method for calculation of the
diffusion ellipsoid by 3-D lattice random walk in the porespace using the program is described. The random walksimulation should only be performed in a percolated porecluster responsible for the macroscopic or long-distancetransport of materials (not in isolated small pores, see
Figure 1). Thus, the raw 3-D CT images of porous samples
need to be preprocessed by cluster labeling to exclude
isolated pores [e.g., Ikeda et al., 2000; Nakashima andKamiya, 2007]. The DellipsoidM6.nb program is available
at http://staff.aist.go.jp/nakashima.yoshito/progeng.htm tofacilitate the X-ray CT approach to the pore-scale analysisof anisotropic diffusive transport in porous media. For
further information about data preparation including cluster
labeling and program output, readers are requested to referto the ‘‘readme’’ file available at the URL above.
[5] DellipsoidM6.nb is an advanced version of Rwalk.nb
developed by Nakashima and Kamiya [2007]. The random
walk algorithm used is the nonsorbing 3-D lattice-walk on apercolated pore cluster in a simple cubic lattice image
system with a periodic boundary condition using a mirror
operation. An excellent pseudorandom number generator,
Mersenne Twister, developed by Matsumoto and Nishimura[1998] was implemented in DellipsoidM6.nb. The random
walker should be nonsorbing because the purpose is tocalculate the geometrical tortuosity of the pore structure andthe undesirable effects of the sorption of walkers on thesolid surface should be eliminated. The difference between
DellipsoidM6.nb and Rwalk.nb is that while the former cancalculate the self-diffusivity along any direction in the 3-D
space, the latter calculates the self-diffusivity only along thethree orthogonal directions.
[6] The self-diffusivity along any direction in the 3-D
space is calculated as follows. Each 3-D random walktrajectory is projected on the a-axis originating from thecenter of the 3-D image system (Figure 2a). The displace-

2 of 19

W12435

NAKASHIMA ET AL.: DIFFUSION ELLIPSOIDS OF ANISOTROPIC POROUS ROCKS

ment on the a-axis is used to calculate the walkers’ meansquare displacement ha(t)2i in the a-directionhaðt Þ2 i ¼

n
1Xðai ðt Þ  ai ð0ÞÞ2 ;
n i¼1

ð1Þ

W12435where n is the number of random walkers, t is thedimensionless integer time (t = 0, 1, 2, 3, . . .), ai(t) is thecoordinate of the walker’s position on the a-axis at time t
for the ith walker, and ai(0) is the initial position of the ith
walker. The exact solution for a lattice walk in free space(i.e., porosity = 100 vol %), ha2ifree, is given by ha2ifree =‘2t/3 where ‘ is the lattice constant of the simple cubic
lattice (i.e., the dimension of a cubic CT voxel). The selfdiffusion coefficient is proportional to the time derivative of
the mean square displacement [e.g., Nakashima andKamiya, 2007; Ohkubo, 2008]. ‘‘Normalized diffusivity’’in the present study refers to the self-diffusivity in porous
media divided by that in free space, and the tortuosity is its
reciprocal. Thus, we havedhaðt Þ2 i
1
dt
normalized diffusivity ¼¼
tortuosity‘2
3

as t ! 1:
ð2Þ

The normalized diffusivity is unity for the random walk infree space, and less than unity for that in porous media
owing to the obstructing effect of the solid. The asymptote
value (i.e., t ! 1) is essential because the normalized
diffusivity can be overestimated owing to the short diffusion
distance compared with the pore network size if the time
derivate is taken at a small t value [e.g., Nakashima and
Watanabe, 2002]. For an identical set of random walktrajectories, the a-axis is scanned on the upper-half plane of
the 3-D polar coordinate system (Figure 2b) to express the
directional mean square displacement as a function of theangles q and f. Equation (2) is again applied to the meansquare displacement to obtain the normalized self-diffusion
coefficient depending on q and f. Typically q covers 0, 1, 2,. . ., 90°, and f covers 0, 1, 2, . . ., 359°. Thus, the number ofthe data points we obtain for the normalized diffusivity is aslarge as 90 360 + 1 = 32,401, which is likely sufficient toevaluate the diffusive anisotropy of porous rocks.
[7] The normalized self-diffusion coefficient depending
on q and f can be visualized by a polar representation,
namely, diffusion ellipsoid [e.g., Westin et al., 2002] a 2Aþ 2  
b
c 2
þ
¼ 1;
B
C

ð3Þ

where A, B, and C are the radii on the a-, b-, and c-axes,
respectively (Figure 2c). The order of the three principal
self-diffusion coefficient values is A B C. The data set
Figure 2. Principles of quantification of the anisotropic
self-diffusivity using DellipsoidM6.nb. (a) A random walktrajectory is projected on the a-axis originating from thecenter of the 3-D image system to calculate the walkers’
mean square displacement in the a-direction. (b) The a-axis
is scanned on the upper-half plane of the 3-D polar
coordinate system to express the directional mean squaredisplacement as a function of q and f. (c) The directiondependent mean square displacement is converted into a
diffusion ellipsoid. Generally, the principal axes (i.e., a-, b-,
and c-axes) of the ellipsoid do not coincide with the x-, y-,and z-axes of the original 3-D CT image system.
3 of 19

W12435

NAKASHIMA ET AL.: DIFFUSION ELLIPSOIDS OF ANISOTROPIC POROUS ROCKSof the normalized self-diffusion coefficients calculated byequation (2) was then fitted to equation (3) in terms of thesecond moment of the object as described by Ikeda et al.

W12435[2000] and Ohtani et al. [2001]. First, the search for a
principal axis yielding the minimum value of the second
moment was performed to assign it to the c-axis (i.e., themajor axis). Second, a principal axis yielding the maximum
value of the second moment (i.e., a-axis or minor axis) was
searched in the direction perpendicular to the c-axis. Finally,
the third principal axis perpendicular to both the a- and c-axes
was automatically chosen to be the b-axis (i.e., intermediateaxis). The distances from the origin of Figure 2c to the crossing
points of the ellipsoid (i.e., raw data on the normalized selfdiffusion coefficients calculated by equation (2)) and thethree principle axes were determined to be radii A, B, and C.
[8] With respect to the determination method of theprincipal axes and diffusivity eigenvalues, our method using
the second moment of the ellipsoidal object is basically
identical to the method conventionally used in diffusiontensor magnetic resonance imaging studies [e.g., Westin etal., 2002]. However, there is a difference between the two
methods. The latter method a priori assumes that the
diffusion ellipsoid is convex (i.e., without constriction).
While this assumption is correct, only seven directions of
the magnetic field gradients including a zero gradient areneeded [e.g., Bhagat and Beaulieu, 2004]. However, theassumption breaks down for highly anisotropic porous
rocks. Thus, directions as many as 90 360 + 1 = 32,401
were required for the a-axis scanning to measure the 3-D
shape of the normalized diffusivity accurately in the present
study. The simple shape-fitting method of the normalized
diffusivity to the convex ellipsoidal surface is inappropriate
for diffusion ellipsoids with constriction. It should be notedthat our search method for the principal axes works
properly even for diffusion ellipsoids having a constriction
because the fitting algorithm is based on the minimization/
maximization of the second moment (not ellipsoidal shapefitting).

3. Application to CT Images of Pumice,
Sandstones, and a Bead Pack
[9] DellipsoidM6.nb was applied to 3-D CT image sets offour natural rocks having anisotropic pore structures (onepumice from Ito and three sandstones from Berea, Shirahama, and Tako) and a pack of synthetic monosized beads
to calculate the diffusion ellipsoid. The Ito pumice sample
from Kagoshima, Japan [Yamamoto et al., 1978] has anelongated pore structure formed by deformation during
eruption. A cylindrical sample (diameter 6 mm, length

Figure 3. The Berea sandstone sample used in the present
study. The z-direction is normal to the bedding plane. (a) X-ray
CT image after the cluster-labeling process. Black areas are
the largest or percolated pore cluster, white areas are solidmineral grains, and red areas are isolated pore clusters(7379 isolated pore clusters were identified). The dimensionof the image system is 4203 voxels = 4.13 mm3. (b) Crossplotof the volume and surface area of 7380 pore clusters. The
theoretical upper and lower limits [Nakashima and Kamiya,2007] for sufficiently large clusters are shown by dotted lines.
(c) Wireframe rendering of the autocorrelation ellipsoid
projected on the x-z plane. The principal axes of the ellipsoidfrom Table 1 are shown. Note that the shape of the ellipsoid is
prolate in the z-direction.
4 of 19

W12435

NAKASHIMA ET AL.: DIFFUSION ELLIPSOIDS OF ANISOTROPIC POROUS ROCKS

Figure 4

5 of 19

W12435

W12435

NAKASHIMA ET AL.: DIFFUSION ELLIPSOIDS OF ANISOTROPIC POROUS ROCKS

W12435

Table 1. Comparison of the Diffusion Ellipsoid (Normalized Diffusivity) and Autocorrelation Ellipsoida
SampleBerea Sandstone

Minor axis (A, q, f)Intermediate axis (B, q, f)
Major axis (C, q, f)(0.142, 87°, 189°)(0.144, 71°, 98°)
(0.160, 19°, 288°)

Minor axis (A, q f)
Intermediate axis (B, q, f)
Major axis (C, q, f)(34.8, 82°, 36°)(36.6, 78°, 128°)(40.3, 14°, 272°)
1.13
1.05

f0
f1

ShirahamaSandstone

Ito Pumice

Tako Sandstone

Bead Pack

Diffusion Ellipsoid(0.019, 1°, 167°)(0.165, 90°, 104°)(0.031, 89°, 347°)(0.181, 89°, 14°)(0.032, 90°, 257°)(0.467, 1°, 186°)(0.015, 2°, 86°)
(0.040, 90°, 188°)(0.044, 88°, 278°)(0.577, 70°, 101°)(0.587, 81°, 8°)(0.595, 22°, 255°)

Autocorrelation Ellipsoid(9.86, 10°, 225°)(50, 88°, 283°)(15.1, 90°, 315°)(63, 88°, 193°)(15.2, 80°, 45°)
(130, 3°, 64°)(33.7, 12°, 72°)
(42.9, 80°, 283°)(44.5, 84°, 192°)(560, 82°, 305°)(566, 88°, 35°)
(587, 9°, 141°)

0.77
0.61

1.00
1.00

0.65
0.77

2.32
2.60

a

Refer to equation (3) and Figure 2b for (A, B, C) and (q, f), respectively. The unit of (A, B, C) for diffusion ellipsoids is dimensionless, while that for
autocorrelation ellipsoids is mm. The degree of the deformation for the pore diameter ( f0) and for the pore network ( f1) is also indicated.

10 mm) was prepared so that the longitudinal direction of
the cylinder would be parallel to the pore elongation
direction. Three cylindrical sandstone samples from Berea(USA), Shirahama (Wakayama, Japan), and Tako (Gunma,
Japan) of similar dimensions were cored normal to thebedding plane for CT imaging. Some petrophysical dataon the sandstones (e.g., pore size distribution and mineralogy) have been published elsewhere [Xue et al., 1992; Linand Takahashi, 2000]. A random pack of monosized beads(bead diameter 2.11 mm, porosity 37 vol %) was chosen as
an analog of the unconsolidated well-sorted sandy sediment
before diagenesis, which yields an isotropic diffusion ellipsoid (i.e., sphere).
[10] Cylindrical samples of Ito pumice and Berea sandstone having a diameter of 6 mm were imaged by microfocus X-ray CT (Toscaner-30000, Toshiba IT and Control
System Corporation) at a cubic voxel size of ‘ = 9.8 mm
with 50 kV acceleration voltage and 0.35 mA tube current.
The original reconstructed 3-D images (10241024voxels, 691 slices) were trimmed to extract a cubic image
system of 4203 voxels = 4.13 mm3. The two cylindricalsandstone samples with a diameter of 4 mm from Shirahama and Tako were imaged by the SP-mCT system(SPring-8 synchrotron radiation facility [Uesugi et al.,2001]) at a cubic voxel size of ‘ = 5.8 mm with 25 keV
monochromatic X-ray energy. The reconstructed original
3-D images (1000 1000 voxels, 700 slices) were trimmed
to extract a cubic image system of 4353 voxels = 2.53 mm3.
As for the monosized bead pack sample, a 3-D cubic image
system of 2563 voxels = 13.63 mm3, which was obtained byanother microfocus X-ray CT scanner (Nittetsu Elex Co.,
Ltd.) and reported elsewhere [Nakashima and Watanabe,2002], was reused in the present study.
[11] Cluster labeling was applied to each trimmed cubic
image system before the random walk simulations toeliminate isolated small pores (see, for example, Figure 1)
and to extract a single percolated pore cluster responsible

for the long-distance material transport. A threshold forvoxel intensity is needed to distinguish the pore and solid incluster labeling. We chose threshold values that yield anexperimentally determined effective porosity (volume fraction of a percolated pore cluster) values measured by theweight difference between water-saturated and dried rocksamples (i.e., 20, 12, 21, and 54 vol % for the Berea,
Shirahama, Tako, and Ito samples, respectively).
[12] An example of cluster labeling is shown in Figure 3
for the Berea sandstone sample. The percolated pore cluster
determined extends through the whole image system(Figure 3a), and serves as the 3-D network for diffusivetransport. All pore clusters including the single percolated
pore cluster and 7379 isolated small pore clusters are plottedin Figure 3b. The surface area in Figure 3b refers to the area
sum of numerous squares (each dimension, ‘‘) on thesurface of each blocky pore cluster. The surface area value
was calculated without considering the contribution of thepore faces exactly on the surface of the 4203 image system(i.e., ‘‘definition A’’ of Nakashima and Kamiya [2007])
because such pore faces disappear in random walk simulations with a mirror boundary condition employed. Cluster
labeling revealed that the total porosity (i.e., sum of thevolume fraction of a percolated pore and of numerous
isolated small pores) is nearly equal to the effective porosity(21 versus 20 vol %), implying that almost all of the poresare connected to form a single large pore cluster in the Berea
sandstone. Results of cluster labeling for other rock samples
are shown in Figure 4. The Shirahama sandstone sample
shows a finer pore structure than the Berea sample, aspreviously revealed by mercury intrusion porosimetry [Linand Takahashi, 2000]. There are many prolate pores visibly
elongated in the z-direction of the image of the Ito pumicesample. The Tako sample contains large pores partly filledwith precipitated minerals [Nakashima and Kikuchi, 2007]
and narrow throats connecting those large pores. The
number of isolated pore clusters was as large as 245,985

Figure 4. Porous rock samples used in the present study. A pair of X-ray CT images after (a,c, and e) cluster labeling and(b, d, and f) wireframe rendering of the autocorrelation ellipsoid (Table 1) projected on the x-z plane is shown. The
dimensions of the labeled 3-D image systems are 4353 voxels = 2.53 mm3 for the Shirahama sandstone, 4203 voxels = 4.13
mm3 for the Ito pumice, and 4353 voxels = 2.53 mm3 for the Tako sandstone. The coloring rule of the 3-D image is identicalto that of Figure 3a. As for the two sandstone samples, the z-direction is normal to the bedding plane. Note that the
autocorrelation ellipsoids are oblate (Shirahama sandstone and Tako sandstone) or prolate (Ito pumice) in the z-direction.
6 of 19

W12435

NAKASHIMA ET AL.: DIFFUSION ELLIPSOIDS OF ANISOTROPIC POROUS ROCKS

Figure 5. Anisotropic pore structure of a bead packelongated or compressed in the z-direction. The coloring
rule of the 3-D image is identical to that of Figure 3a. The
elongation/compression factor f is indicated. The image
systems for f = 1.9 and 0.54 are composed of 2562 486
138 cubic voxels, respectively. The originaland 2562
image (i.e., f = 1) of 2563 cubic voxels = 13.63 mm3 having
an isotropic pore structure is also shown.
for Shirahama sandstone, 14,988 for Ito pumice, and172,454 for Tako sandstone. However, the total porosity
was nearly equal to the effective porosity (15 versus 12 vol% for Shirahama sandstone, 54 versus 54 vol % for Itopumice, and 23 versus 21 vol % for Tako sandstone),
implying again that almost all of the pores were connected
to form a single large pore cluster.

W12435[13] In accordance with Ikeda et al. [2000], the two-point
autocorrelation was calculated for the percolated pore cluster images before random walk simulations. The voxel valueof the largest (i.e., percolated) pore cluster was set to unity;that of the solid voxels and isolated pore voxels was zero.
The two-point autocorrelation was calculated as (average ofthe product of the voxel value and the spatially shifted voxel
value)/(effective porosity calculated by the volume fractionof the percolated pore cluster). The autocorrelation
obtained, a function of the space shift, is a 3-D image,and its voxel value decreases from unity to the effective
porosity value with increasing spatial shift. As used by
Ikeda et al. [2000], an autocorrelation value of 0.5 waschosen as a threshold to extract a high correlation zone forthe three sandstone samples and bead pack; the threshold
was taken to be 0.7 for Ito pumice because the two-point
autocorrelation value never decreased to 0.5 owing to theextraordinarily high porosity of pumice. We have confirmedthat the results (i.e., the fitted q and f values and relative
magnitude of radii of autocorrelation ellipsoids in Table 1)
do not show a strong dependence on the chosen thresholdvalue. The ellipsoid fitting algorithm using the second
moment described above was applied to an object with ahigh correlation value of 0.5 (or 0.7) to output theellipsoid defined by equation (3). The fitting results arelisted in Table 1, and shown in Figures 3 and 4. Although
the cluster-labeled pore CT images of Figures 3 and 4 show
no obvious pore anisotropy except for Ito pumice, theautocorrelation analysis (Table 1) clearly reveals that thefour rock samples are all anisotropic in the z-direction(prolate for Berea and Ito, oblate for Shirahama and Tako).
[14] The comparison of the pore diameter in the labeledimages (Figures 3– 5) with A, B, and C values in Table 1 forthe two-point autocorrelation ellipsoid indicates that they
are almost of the same order. Thus, the anisotropy inferredfrom the autocorrelation ellipsoid is the anisotropy of thepore diameter averaged over numerous individual pores.
The degree of the deformation of the pore diameter can be
estimated quantitatively using the autocorrelation data. It isreasonable to assume that initially the pore structure was as
isotropic as that in the randomly packed bead pack ( f = 1.0
in Figure 5), and was significantly elongated or compressedin the z-direction during the geological deformation processes. Then the deformation degree of the pore diameter,f0, can be calculated by f0 = C(AB)1/2 for prolate autocorrelation ellipsoids (Berea and Ito samples) and by f0 =A(BC)1/2 for oblate autocorrelation ellipsoids (Shirahama
and Tako samples). The geometric mean, (AB)1/2 or (BC)1/2,
was taken to conserve the volume of the ellipsoid. In
general, the pore diameter anisotropy is not exactly identicalto the pore network anisotropy (Figure 1) responsible forthe diffusion anisotropy. It is thus difficult to connect theanisotropy of the autocorrelation ellipsoid directly to the
diffusion anisotropy. However, there should be positive
correlation between the pore diameter anisotropy and pore
network anisotropy because they are degrees of geologicaldeformation recorded in the same rock mass. The obtained
f0 values listed in Table 1 were used to discuss the diffusive
anisotropy via pore network anisotropy.
[15] We examined the diffusion anisotropy of the porous
samples as follows. The program DellipsoidM6.nb was run
by importing each labeled CT image in Figures 3– 4 to

7 of 19

W12435

NAKASHIMA ET AL.: DIFFUSION ELLIPSOIDS OF ANISOTROPIC POROUS ROCKS

Figure 6
8 of 19

W12435

W12435

NAKASHIMA ET AL.: DIFFUSION ELLIPSOIDS OF ANISOTROPIC POROUS ROCKS

obtain the mean square displacement for nonsorbing
walkers diffusing in the percolated pore cluster. The numberof the random walkers n in equation (1) was as large as
400,000 to reduce the undesirable stochastic fluctuation of
the mean square displacement [Ohkubo, 2008]. The starting
position of the random walk was chosen randomly among
the voxels of the percolated pore cluster. A personal
computer (PC) with an Intel Core21 Duo T7600 CPU(2.33 GHz) and 2 GB RAM running Windows XP1 wasused for the random walk simulations. In order to avoid an
undesirable overestimate of the normalized diffusivity due
to the short diffusion distance compared with the pore
network size [Nakashima and Kamiya, 2007], sufficiently
long values were employed for the random walk duration
to allow the walkers to travel long distances to obtain thelinear (not convex) mean square displacement. The directional mean square displacement was obtained for a number of discrete integer (q, f) pairs (i.e., q = 0, 1, 2, . . . 90°
and f = 0, 1, 2, . . . 359°) in Figure 2b. Then, the time
derivative of the linear directional mean square displacement was calculated according to equation (2) to obtainthe normalized diffusivity as a function of q and f. The 3-D
shell-like object obtained was the raw data on the diffusion
ellipsoid, and was used for the ellipsoid fitting process asfollows. A (q, f) pair that minimized the second moment of
this shell-like object was chosen as the c-axis; another (q, f)pair that was normal to the c-axis and yielded the maximum
second moment was taken to be the a-axis. The intermediateaxis (b-axis) was automatically determined to be perpendicular to both the a- and c-axes. The distances from theorigin of Figure 2c to the crossing points of the 3-D shelllike object and the a-, b-, and c-axes were taken to be theradii A, B, and C in equation (3), respectively.
[16] One of the purposes of the present study is to discuss
the effect of pore structure anisotropy on the anisotropy ofthe diffusion through the percolated pore cluster. To this end,
we calculated the diffusion ellipsoids for 3-D CT images
systematically deformed by an affine transformation. The
bead pack image [Nakashima and Watanabe, 2002] has no
significant pore structure anisotropy in terms of the autocorrelation ellipsoid (Table 1), and yields a diffusion ellipsoid ofan almost perfect sphere. We generated a series of synthetic
bead pack images that were uniaxially elongated or compressed in the z-direction. The case of f > 1 implies
elongation and f < 1 implies compression, where f is theelongation/compression factor. Examples for some f values
are shown in Figure 5. The program DellipsoidM6.nb wasapplied to the image sets to obtain the diffusion ellipsoids
averaged over 200,000 walkers for various f values. We
analyzed the results and discussed the effects of f on theellipsoid shape, focusing on the emergence of the constric-

W12435tion that has been reported for water self-diffusion in clay
gels having anisotropic (nematic) pore structure.
[17] We also applied an affine transformation in thez-direction with a factor of f to the natural rock CT images
in Figures 3 and 4. Unfortunately, according to the autocorrelation ellipsoid analysis shown in Table 1, all rocksamples possessed intrinsic pore structure anisotropy withrespect to the z-axis, presumably caused by geologicalprocesses during the diagenesis in sedimentary basins anderuption of magma. Thus, corrections must be made toestimate the true strain from the state having an isotropic
pore structure. The isotropic pore state yields an isotropic
diffusion ellipsoid (i.e., a sphere). Thus, if a diffusion
ellipsoid becomes spherical (i.e., A = B = C) for a 3-D
image that is elongated or compressed with a specificfactor, 1/f1, in the z-direction, it is reasonable to assume thatthe 3-D pore image is in the no-strain, or isotropic, state.
After some trial and error, we found the intrinsic elongation/
compression factor, f1, for each sample, as listed in Table 1.
The product, f1 f, was used to express the true degree ofpore elongation/compression from the isotropic pore structure state.
[18] The quantity, f1, is an important measure to characterize the pore network anisotropy, which is responsible forthe diffusive anisotropy. Because the isotropic pore image
set was successfully obtained by the affine transformation
applied in the z-direction, the direction of the rock deformation is the z-direction for all samples in terms of the pore
network anisotropy. The estimated magnitude (i.e., f1 values
listed in Table 1) and the direction of the rock deformation
were used to discuss the correlation with the diffusive
anisotropy and pore diameter anisotropy.

4. Results[19] Results for the Berea sandstone sample are shown in
Figure 6. The fitted parameters of the diffusion ellipsoid aretabulated in Table 1. The CPU time of DellipsoidM6.nb
required for the Berea sandstone sample was 4.4 h for the
PC used. The mean square displacements along the x-, y-,and z-axes in Figure 3a, as calculated by equation (1), areshown in Figure 6a. Equation (1) becomes simple if the
cubic voxel dimension, ‘, (9.8 mm for the Berea sandstone)
is taken to be unity. Thus, ‘ = 1 was assumed in Figure 6a,yielding a dimensionless mean square displacement. Thereis a slightly convex portion of the mean square displacement
curve in the early stage (i.e., t < 10,000). This time interval
is derived from the finite size of the pore network [e.g.,
Nakashima and Kamiya, 2007]. The pore network of mostsedimentary rocks consists of large subspherical pools andnarrow throats bridging the pools [Lindquist et al., 2000].
The percolated pore cluster of Figure 3a suggests that the

Figure 6. Results of the random walk simulation for the percolated pore cluster in Figure 3a (Berea sandstone).
(a) Dimensionless mean square displacement along the directions of x-, y-, and z-axes averaged over 400,000 walkers as afunction of t. (b) Snapshot of the 400,000 walkers at t = 30,000. The unit of the axes is the cubic voxel dimension (i.e., ‘ =9.8 mm for the Berea sandstone). The origin of the 3-D image is the starting position of the random walk of each walker.
(c) Unwrapped expression of the Northern Hemisphere (i.e., q 90°) of the normalized diffusivity calculated by equation(2) as a function of q and f (in units of degrees). This is the raw data (i.e., before ellipsoid fitting) and normalized
diffusivity values for 32,401 integer (q, f) pairs (i.e., q = 0, 1, 2, . . . 90° and f = 0, 1, 2, . . . 359°) are plotted. (d) Bird’seye-view of the wireframe rendering of the wrapped normalized diffusivity. The three principal axes (a, b, and c) of thefitted diffusion ellipsoid are superimposed. (e) Wrapped diffusion ellipsoid viewed along the principal axes.
9 of 19

W12435

NAKASHIMA ET AL.: DIFFUSION ELLIPSOIDS OF ANISOTROPIC POROUS ROCKS

pool diameter is typically 10 to 30 voxels and the typicaldistance between the pools is 20 to 40 voxels. Random
walkers should travel more than 40 voxels to reach theasymptote in terms of equation (2). On the other hand, the3-D mean square displacement of the walkers, hx2i + hy2i +hz2i, is 1800 at t = 10,000 (Figure 6a). Thus, the characteristic travel distance of the walkers is 18001/2  42 voxels
at t = 10,000. Because 42 voxels are larger than the typicalvalues of the pool diameter and the distance between
pools, it is reasonable to conclude that the mean squaredisplacement reaches the asymptote for t > 10,000. The
undesirable early stage (i.e., t < 10,000) having a convex
mean square displacement should be eliminated from thecalculation of normalized diffusivity because the time
derivative (equation (2)) has not yet reached the asymptote
value. Thus, to obtain a reliable diffusion ellipsoid, thetime derivative was calculated using the very last portionof the mean square displacement time series data for all
simulations of the present study, namely, the slope value ofha(t)2i between the maximum time step and 75% of it was
employed. For example, the normalized diffusivity wascalculated using the difference of ha(t)2i at t = 30,000
and 30,000
0.75 = 22,500 for the Berea sandstone.
[20] Figure 6a suggests that the diffusion ellipsoid is
prolate in the z-direction because the mean square displacement in the z-direction is significantly larger than those in
the x- and y-directions. The final positions of all the walkers
are shown in Figure 6b. The overall shape of the diffusion
front can be seen. Although Table 1 indicates that the
diffusion ellipsoid is prolate for the Berea sample, it isdifficult to see the prolate nature in Figure 6b because thedegree of the diffusive anisotropy is small (i.e., C  A). The
unwrapped and wrapped diffusion ellipsoids are shown in
Figures 6c and 6d, respectively. The prolate shape of the
diffusion ellipsoid in Figure 6a is confirmed in Figures 6c
and 6d. The prolate shape of the ellipsoid can be seen moreclearly when it is viewed from the principal axes (Figure 6e).
[21] We carried out a random walk simulation in the freelattice space (i.e., porosity 100 vol %), which theoreticallyyields a complete sphere with a radius of unity for the
diffusion ellipsoid. The results of DellipsoidM6.nb showedthat the fitted major and minor axis values of the normalized
diffusivity ellipsoid averaged over 400,000 walkers wereC = 1.01 and A = 0.99, respectively. These are very close
to the theoretical value of unity, demonstrating that ourrandom walk simulations were performed reliably.
[22] Results for the Shirahama sandstone sample areshown by Figure 7 and Table 1. The CPU time of DellipsoidM6.nb required for this sample was 14 h for the PC
used. The normalized diffusivity was calculated using thedifference between ha(t)2i at t = 300,000 and 300,000
0.75 = 225,000. The distinguishing feature of Figure 7 ascompared with Figure 6 is the oblate shape of the diffusion
ellipsoid compressed in the z-direction. The oblate shape isevident from the significantly small mean square displacement in the z-direction as compared with those in the x- andy-directions (Figure 7a), from the oblate shape of the
diffusion front (Figure 7b), from the small normalized
diffusivity value at the high-latitude (i.e., low-q) zone
(Figure 7c), and from the shape of the diffusion ellipsoid(Figures 7d and 7e).

W12435[23] Results for the Ito pumice sample are shown in
Figure 8 and Table 1. Because the pore size of the vesicularsample is much larger than that of the Berea sample, as
many as 400,000 time steps were required to make the meansquare displacement reach the asymptote line (Figure 8a).
Figures 8a– 8c show that the diffusion ellipsoid is prolate in
the z-direction, as can be easily predicted by the clearly
elongated shape of the individual pores in Figure 4c. The
large elongation factor (f0 = 2.32 and f1 = 2.60, see Table 1)of the pores yields the following important point regarding
the shape of the prolate diffusion ellipsoid. The constrictionof the ellipsoid occurs in the direction of the minor andintermediate principal axes (i.e., a- and b-axes), as depicted
in Figures 8d and 8e. Thus, the diffusion ellipsoid is nolonger convex, like what has been reported for highly
anisotropic porous clay gels [Porion et al., 2001, 2003; deAzevedo et al., 2007].
[24] Results for the Tako sandstone sample are shown in
Figure 9 and Table 1. Figures 9a– 9c show that the diffusion
ellipsoid is oblate, compressed in the z-direction. The large
degree of the compression (f0 = 0.77 and f1 = 0.61, see
Table 1) of the sample reveals a remarkable feature of the
shape of the oblate diffusion ellipsoid. The constriction of
the ellipsoid occurs in the direction of the minor principal
axis (i.e., a-axis), as depicted in Figures 9d and 9e. The
diffusion ellipsoid is no longer convex, as in the case of theIto pumice sample.
[25] Results for the bead pack sample without the affine
transformation (i.e., f = 1) are shown in Figure 10 and
Table 1. The isotropic diffusion front in Figure 10b and thespherical diffusion ellipsoid in Figures 10d and 10e indicatethat the pore structure of the bead pack is isotropic.
Although Figures 10a and 10c show slight anisotropy, it
is reasonable to neglect the anisotropy because A = 0.577 isvery nearly equal to C = 0.595 (Table 1).
[26] A relatively simple thresholding in which all voxels
below (or above) a chosen X-ray linear absorption coefficient value are assigned to pore space (or solid) was
employed to extract a binary image set in the present
study. Although there are more advanced methods for
extracting a pore space [e.g., Oh and Lindquist, 1999; Arnset al., 2001], we chosen this simple method because it
yields acceptable pore characteristics as follows. Following
Nakashima [2000], Polak et al. [2003], and Zhu et al. [2007],
a medical (i.e., low-resolution) X-ray CT scanner was used
to measure the tortuosity of each rock sample preparedfrom the identical block from which the small cylindrical
CT sample was cored (Y. Nakashima, unpublished data,2008). Invasion of nonsorbing iodine ions from a KI-bearing
filter paper on one face of the water-saturated sample into
the sample was monitored by the medical CT. For example,
the Gaussian distribution of iodine concentration within thesample was fitted to obtain a tortuosity of 6 for the Berea
sample. No significant anisotropy was detected, presumably
owing to the inevitable low precision of this method.
According to Table 1, the random walk simulations suggestthat the tortuosity of the Berea sample falls within 1/C 
6.25 and 1/A  7.04. Thus, results of the random walksimulations agree well with those of the laboratory diffusion
test, indicating the spatial resolution and thresholding methodof the present study are acceptable. As for the bead pack
image set of ‘ = 53 mm, we have shown that our thresholding

10 of 19

W12435

NAKASHIMA ET AL.: DIFFUSION ELLIPSOIDS OF ANISOTROPIC POROUS ROCKS

Figure 7. Same as Figure 6 but for the Shirahama sandstone and with snapshot of the 400,000 walkers
at t = 300,000.
11 of 19

W12435

W12435

NAKASHIMA ET AL.: DIFFUSION ELLIPSOIDS OF ANISOTROPIC POROUS ROCKS

Figure 8. Same as Figure 6 but for the Ito pumice and with snapshot of the 400,000 walkers at t =400,000.

12 of 19

W12435

W12435

NAKASHIMA ET AL.: DIFFUSION ELLIPSOIDS OF ANISOTROPIC POROUS ROCKS

Figure 9. Same as Figure 6 but for the Tako sandstone and with snapshot of the 400,000 walkers at t =600,000.

13 of 19

W12435

W12435

NAKASHIMA ET AL.: DIFFUSION ELLIPSOIDS OF ANISOTROPIC POROUS ROCKS

Figure 10. Same as Figure 6 but for a bead pack and with snapshot of the 400,000 walkers at t =50,000.
14 of 19

W12435

W12435

NAKASHIMA ET AL.: DIFFUSION ELLIPSOIDS OF ANISOTROPIC POROUS ROCKS

method and spatial resolution are sufficient to reproduce theexperimentally measured tortuosity and surface-to-volume
ratio values by the random walk simulations [Nakashima

W12435and Watanabe, 2002]. The laboratory KI diffusion tests for
the Ito pumice sample showed that the tortuosity parallel tothe pore elongation direction and that normal to the poreelongation direction were 4 and 13, respectively.
According to the random walk simulations for the pumicesample (Table 1), the tortuosity values for the corresponding
directions are 1/C  2.1 and 1/(AB)1/2  5.8, respectively,
which are less than those obtained by the laboratory tests.
The discrepancy is probably derived from the simple
thresholding method and relatively low spatial resolutionof the CT system employed. Similar discrepancy of the
tortuosity was also observed for Shirahama and Takosamples, suggesting that submicrometer pore structure
was not accurately imaged. However, the maximum-tominimum diffusivity ratio by the KI diffusion experimentagrees well with the (BC)1/2A1 or C(AB)1/2 values of theoblate/prolate diffusion ellipsoids in Table 1 (13/4  3
versus 5.8/2.1  2.8 for the Ito sample, 2 versus 1.7 for
the Shirahama sample, and 3 versus 2.8 for the Takosample). Thus, it is reasonable to conclude that the grossanisotropy of the pore structure derived from the detectablelarge pores was successfully obtained by the employed CT
apparatus with the simple thresholding method, which is
sufficient for the present study because the purpose is topresent the principle of the new method to calculate thediffusion ellipsoids quickly, to evaluate the general trends inthe diffusive anisotropy in relation to the degree of the rockdeformation, and to discuss the cause of the constriction ofdiffusion ellipsoids rather than to reproduce the individual
diffusivity values measured by the laboratory experiments.
[27] When the bead pack image was elongated or compressed by a factor f in the z-direction, as shown in Figure 5,
the spherical diffusion ellipsoid of Figure 10d broke down.
The dependence of the shape of the diffusion ellipsoid on thef value is shown in Figure 11. The direction of the affine
transformation is along the z-direction. Therefore, the resultant diffusion ellipsoid for the transformed CT image
should be symmetric with respect to the z-axis, andindependent of f. The q dependence of the ellipsoid
averaged over f values is shown in Figures 11a and 11b
for elongation and compression, respectively. The polar
representations of the diffusion ellipsoids, Figures 11a and11b, clearly show that the sphere becomes distorted (prolate
or oblate) if the f value differs from unity. It should be
noted that when the f value exceeds a critical value, theconstriction of the ellipsoid appears in the direction ofminimum diffusivity (i.e., q = 90° for Figure 11a and 0° for
Figure 11b). We express the deviation of the shape of the
Figure 11. Diffusion ellipsoid for random walk in the
percolated pore cluster of the bead pack for various f factors.
The ellipsoid was averaged over 200,000 walkers. (a) Polar
representation of the averaged diffusion ellipsoid for theelongated (i.e., f = 1.0, 1.3, 1.9, 2.9, and 5.0) bead pack
images. (b) Polar representation of the averaged diffusion
ellipsoid for the compressed (i.e., f = 1.0, 0.7, 0.54, 0.4, 0.3,
and 0.2) bead pack images. (c) Degree of anisotropy of selfdiffusivity Dz/Dxy plotted for various values of elongation/
compression f for the bead pack CT image. The case for f = 1(i.e., Figure 10) is indicated by an arrow. The solid line is a
guide for the eye. The four ellipsoid regimes with threeboundaries (horizontal broken lines) are illustrated.

15 of 19

W12435

NAKASHIMA ET AL.: DIFFUSION ELLIPSOIDS OF ANISOTROPIC POROUS ROCKS

W12435ples using the same quantity, Dz/Dxy, as in Figure 11. The
f curves,pore structure, and hence the Dz/Dxy versus f1
differ from sample to sample. However, all curves include( f1 f, Dz/Dxy) = (1, 1) owing to the properly chosen f1 values
listed in Table 1. All rock samples showed a common
feature: the same four regimes as in Figure 11 appeared
depending on the f1 f values. It is noteworthy that the regime
transition occurs at a common value with respect to Dz/Dxy.
The transition from the prolate ellipsoid without constriction
to that with constriction occurs at Dz/Dxy = 1.49, 1.50, 1.50,
1.54, and 1.37 for the bead pack, Berea, Ito, Shirahama, andTako samples, respectively (average 1.5). The transition
from the oblate ellipsoid without constriction to that withconstriction occurs at Dz/Dxy = 0.64, 0.66, 0.68, 0.60, and0.67 for the bead pack, Berea, Ito, Shirahama, and Takosamples, respectively (average 0.65). The boundaries at
Dz/Dxy = 1.5 and 0.65 are indicated in Figure 12. Because
the reciprocal of 0.65 is nearly equal to 1.5, we canconclude that the constriction occurs when the maximum
to minimum diffusivity ratio of the diffusion ellipsoid
exceeds 1.5.

5. Discussion

Figure 12. Degree of anisotropy of self-diffusivity Dz/Dxy
averaged over 200,000 walkers as a function of the degreeof elongation/compression f1 f of CT images for the five
porous samples. The case for f = 1 (i.e., Figures 6 –10) areindicated by an arrow for each sample. The solid lines are
guides for the eye. Four ellipsoid regimes with threeboundaries (horizontal broken lines), namely Dz/Dxy = 0.65,
1.0, and 1.5, are illustrated.

diffusion ellipsoid from the isotropic sphere by a single
quantity, Dz/Dxy, where Dz is the normalized diffusivity inthe z-direction (i.e., normalized diffusivity for q = 0°) andDxy is the normalized diffusivity averaged on the x-y plane(i.e., normalized diffusivity averaged over f = 0 to 359° at
q = 90°). The quantity, Dz/Dxy , is greater than unity forimage elongation, and less than unity for image compression. The Dz/Dxy versus f relationship is shown in Figure 11c
for a wide range of f values. Four regimes were identified interms of the shape of the diffusion ellipsoid. The ellipsoidsfor f = 1.9 and 0.54 in Figure 5 are located at theboundary of the regimes for with and without constriction.
[28] Results for rock samples having nonunity f1 values
are shown in Figure 12. Both the direction of the intrinsic
pore diameter anisotropy inferred from the autocorrelation
ellipsoid (Table 1) and the direction of the affine transformation lie along the z-direction. Thus, it is reasonable tocharacterize the diffusion ellipsoid shape of the rock sam-[29] One of the purposes of the present study is toconfirm the hypothesis [Yokoyama and Nakashima, 2005]
that the pore structure anisotropy is responsible for the
anisotropic diffusivity. It should be noted that there are
two independent measures for the pore structure anisotropyof porous media: anisotropy of the pore diameter and that ofthe pore network. For example, in Figure 1, the former isthe anisotropy of the pipe diameter, and the latter is theanisotropy of the geometrical tortuosity of the pipes. The
two measures can be expressed as f0 and f1, respectively.
The pore network anisotropy (not pore diameter anisotropy)
is responsible for the long-distance diffusive anisotropy
calculated by the asymptote of equation (2). It is reasonable
to assume that initially the pore diameter and pore network
were as isotropic as those in the randomly packed beadpack, and were significantly elongated or compressed in thez-direction during the geological deformation processes.
This assumption is reasonably supported by that (1) the f0
value is almost identical to f1 value for each sample (Table 1)and (2) the direction of the elongation/compression of the
autocorrelation ellipsoid is the z-direction identical to thedirection of the pore network deformation. Thus, it can beconcluded that the pumice and sandstone samples areelongated or compressed in the z-direction, and the degreeof the geological deformation is given by the almost
identical values f0 and f1.
[30] The hypothesis that the pore structure anisotropy is
responsible for the anisotropic diffusivity predicts thatthe direction of the geological deformation agrees with thedirection of the major (or minor) principal axis of the prolate(or oblate) diffusion ellipsoid. As described, the direction of
the geological deformation is the z-direction for all rocksamples. The parameters of the diffusion ellipsoids listed inTable 1 clearly show that the direction of the major or minor
principal axis of the diffusion ellipsoid is almost parallel tothe z-direction (i.e., q  0°). For example, q = 19° for the
major axis of the Berea sample, q = 1° for the minor axis of
the Shirahama sample, q = 1° for the major axis of the Itosample, and q = 2° for the minor axis of the Tako sample.

16 of 19

W12435

NAKASHIMA ET AL.: DIFFUSION ELLIPSOIDS OF ANISOTROPIC POROUS ROCKS

Figure 13. Schematic diffusion path of pore fluid
molecules in porous media containing highly anisotropic
solid (shaded) obstacles. P1 is a straight diffusion path along
the elongated solid. P2 is a tortuous path avoiding the solid
obstacles.
Thus, the results of all the rock samples show that the
deformation-induced pore structure anisotropy is responsible
for the anisotropy of the diffusion of the nonsorbing pore
fluid molecules along a percolated pore.
[31] The anisotropy of the sandstone samples is not very
obvious by visual inspection in the 3-D CT images in
Figures 3 and 4. However, we clearly detected the anisotropyof the pore structure, which is seen to be oblate or prolatewith respect to the bedding plane in Figures 6, 7, and 9 andTable 1. For example, a prolateness in diffusion of as smallas (C  A)/C = 11% was successfully detected for the Berea
sandstone sample. The smooth Dz/Dxy versus f1f curves of
Figure 12 free of undesirable stochastic fluctuation showthat a much lower anisotropy can be detected if a sufficiently large number of random walkers are employed. Thishigh sensitivity to pore anisotropy suggests that the present
CT-based method would be particularly effective in evaluating rocks with relatively low anisotropy, for which theconventional technique [e.g., Ohlsson and Neretnieks,
1995] fails to detect the anisotropy because of the low
accuracy.
[32] A diffusion ellipsoid with constriction was observed
for water self-diffusion in synthetic clay gels using nuclear
magnetic resonance spectroscopy [Porion et al., 2001,
2003; de Azevedo et al., 2007]. However, the mechanismof the generation of this pinched-in ellipsoid has not beenstudied in detail. In the present study, we report for the first

W12435time that the constriction occurs for natural porous rocks as
well. We also discuss the generation mechanism in detailby analyzing highly anisotropic pore CT images producedby systematic uniaxial elongation/compression operations
using affine transformation. According to Table 1, thediffusion ellipsoid constriction occurs for samples with avery large f1 value (Ito pumice) or very small f1 value (Takosandstone). This suggests that the constriction appears when
the pore is substantially elongated or compressed. This isstrongly supported by the systematic calculations of diffusion
ellipsoids for uniaxially deformed pore images (Figures 11
and 12). Figures 11 and 12 demonstrate that (1) the constriction occurs in the direction of the minor principal axis evenfor the Berea, Shirahama, and bead pack samples when they
are elongated or compressed significantly and that (2) the
constriction occurs when the maximum to minimum selfdiffusivity ratio of the diffusion ellipsoid exceeds a criticalvalue of 1.5. Thus, the constriction is not particular to the
Ito and Tako samples, but can be found in any porous
geomaterial showing strong diffusive anisotropy.
[33] The constriction of the diffusion ellipsoid for porous
rocks is a remarkable phenomenon because there are few
experimental/theoretical studies on the cause for the onset ofthe constriction. The porous media are, by definition,
heterogeneous because they are composed of solids (impermeable barriers) and pores. Random walkers diffuse inporous media through the pore space by avoiding the solid
obstacles. While a straight diffusion path of tortuosity 1 isensured for diffusion parallel to the elongation, random
walkers diffusing normal to the elongation are forced to goaround solid obstacles with a high aspect ratio, yielding a
tortuous diffusion path (Figure 13). With increasing poreanisotropy, the obstructing effect of the solid selectively
reduces the self-diffusivity exactly normal to the solid
elongation. When the degree of the pore anisotropy exceeds
a critical value, the ellipsoid constriction appears. This is thequalitative explanation for the appearance of the ellipsoid
constriction observed in the direction of minimum diffusivity. A further theoretical study is needed to quantitatively
explain the constriction of the diffusion ellipsoid.
[34] On the basis of the conventional magnetic resonanceimaging studies, a microscopic approach (i.e., pore-scalerandom walk simulating self-diffusion) was employed in thepresent study for the evaluation of the diffusion anisotropy.
There is another possible approach for the diffusion anisotropy study, namely, macroscopic or continuum method
using a parabolic partial differential equation with a diffusion tensor of rank 2 [e.g., Crank, 1975]. In the latterapproach, the directional mean square displacement is not
calculated, but a macroscopic concentration gradient isintroduced as a boundary/initial condition to drive thediffusion. The introduced concentration gradient would bias
the movement of diffusing species to follow the gradient,and thus care must be taken when one compare the random
walk method with the continuum method. For example, it ispossible to obtain an anisotropic diffusion front without
constriction such as Figure 8b by numerically solving theparabolic partial differential equation. However, it is difficult to calculate the 3-D shell of the direction-dependent diffusivity with constriction (e.g., Figure 8d) by the continuum
method because the macroscopic concentration gradient distorts the direction of the diffusion flux vector.

17 of 19

W12435

NAKASHIMA ET AL.: DIFFUSION ELLIPSOIDS OF ANISOTROPIC POROUS ROCKS[35] Thermal diffusivity or thermal conductivity also is an
important transport property of rocks [e.g., Guéguen and
Palciauskas, 1994]. The thermal energy diffuses throughrock-forming minerals, and the thermal diffusivity differs
among minerals and pores. The thermal conduction phenomenon obeys a parabolic partial differential equation as
well as the material diffusion phenomenon. Thus, thepresent study predicts that ellipsoid constriction of thermal
diffusivity also would occur if the material possesses highanisotropy with respect to thermal conductivity. Examples
include pumice made of volcanic glass (conductor) and airfilled pore (insulator) and graphite schist made of graphite(conductor) and silicate (insulator). Assuming zero conductivity for the air-filled pore, a preliminary random walksimulation on the conductive solid framework was performed for the Ito pumice image set of Figure 4c using
the same number of walkers and same time steps as Figure 8.
The obtained ellipsoid of the thermal diffusivity clearly
showed the constriction in the direction normal to the z-axis,
supporting the prediction.
[36] It should be noted that the CT-based method foranisotropy evaluation requires much less time than theconventional laboratory method [e.g., Van Loon et al.,
2004; Yokoyama and Nakashima, 2005] using a diffusion
cell and cut samples. For example, as many as 90 360 + 1 =32,401 diffusivity data points were acquired in 4.4 h for
Figure 6. It is practically impossible to obtain a data set ofthe same magnitude by conventional laboratory diffusion
tests using a diffusion cell and cut samples. Unfortunately
although synchrotron-based submicrometer imaging systems are being developed [Promentilla et al., 2008], thecurrent state of the art of X-ray CT is not always satisfactoryin terms of the spatial resolution. Thus, the conventional
diffusion cell method would not be immediately replacedwith the CT-based method described in this study. However,the current CT technology enables us to accurately image
rock pores larger than several or several tens of micrometers. Thus, the CT-based method can contribute to thequick analysis of the anisotropy of the diffusion throughrelatively large pores.

6. Conclusions[37] High-resolution 3-D pore images were obtained by
microfocus and synchrotron X-ray CT for one pumice andthree sandstone samples. Nonsorbing lattice random walksimulations were performed in a percolated pore cluster of
the CT images to calculate the mean square displacement
along an arbitrary direction in 3-D space. Directiondependent self-diffusion coefficients were calculated bytaking the time derivative of the mean square displacement
projected on any given direction. The second moment wascomputed for the direction-dependent self-diffusivity data toobtain the principal axes of the diffusion ellipsoid. The twopoint autocorrelation ellipsoid was also calculated for the
percolated pore cluster of each rock sample to estimate thepore diameter anisotropy. The pore network anisotropy wasestimated by calculating the diffusion ellipsoid for uniaxiallyelongated or compressed rock images. The degree anddirection of the geological deformation of the samples
estimated by the pore diameter anisotropy analysis agreed
well with those estimated by the pore network anisotropy
analysis. We found that the direction of the geological

W12435deformation coincided with the direction of the major (orminor) principal axis of the prolate (or oblate) diffusion
ellipsoid for each sample. Thus, it can be concluded that the
deformation-induced pore structure anisotropy is responsible
for the anisotropy of the diffusive transport properties. For
the highly anisotropic Ito pumice and Tako sandstone samples, the shape of the diffusion ellipsoid was significantly
distorted (constricted in the direction of minimum diffusivity), and no longer convex ellipsoidal. Systematic calculations of the diffusion ellipsoids for synthetic 3-D rock
images uniaxially elongated or compressed using the affine
transformation were performed to explore the generation
mechanism of this remarkable constriction of the diffusion
ellipsoid. The results revealed that the constriction is notparticular to the Ito and Tako samples, but can occur for any
rocks when the maximum to minimum self-diffusivity ratioof the diffusion ellipsoids exceeds a critical value of 1.5.
We offer the DellipsoidM6.nb program used to calculate thediffusion ellipsoids on the Internet to facilitate the X-ray CT
approach to groundwater hydrology.
[38] Acknowledgments. Comments by O. Nishizawa and three anonymous reviewers were very helpful. A part of this study was supported by
METI of the Japanese government. CT images of Berea sandstone and Itopumice were obtained by S. Kawachi, S. Kasama, and M. Nakano (Nippon
Steel Corporation). The SP-mCT experiments were performed with theassistance of K. Nakamura, K. Uesugi, A. Tsuchiyama, and S. Ikeda at theJapan Synchrotron Radiation Research Institute (proposal 2001B0501NOD-np).

References
Al-Raoush, R. I., and C. S. Willson (2005), A pore-scale investigation of amultiphase porous media system, J. Contam. Hydrol., 77, 67 – 89,
doi:10.1016/j.jconhyd.2004.12.001.
Altman, S. J., W. J. Peplinski, and M. L. Rivers (2005), Evaluation ofsynchrotron X-ray computerized microtomography for the visualizationof transport processes in low-porosity materials, J. Contam. Hydrol., 78,
167 – 183, doi:10.1016/j.jconhyd.2005.05.004.
Arns, C. H., M. A. Knackstedt, W. V. Pinczewski, and W. B. Lindquist(2001), Accurate estimation of transport properties from microtomographic images, Geophys. Res. Lett., 28, 3361 – 3364, doi:10.1029/
2001GL012987.
Arns, C. H., M. A. Knackstedt, W. V. Pinczewski, and N. S. Martys (2004),
Virtual permeametry on microtomographic images, J. Petrol. Sci. Eng.,45, 41 – 46, doi:10.1016/j.petrol.2004.05.001.
Arns, C. H., F. Bauget, A. Ghous, A. Sakellariou, T. J. Senden, A. P.
Sheppard, R. M. Sok, W. V. Pinczewski, J. C. Kelly, and M. A. Knackstedt(2005), Digital core laboratory: Petrophysical analysis from 3D imaging ofreservoir core fragments, Petrophysics, 46, 260 – 277.
Auzerais, F. M., J. Dunsmuir, B. B. Ferreol, N. Martys, J. Olson, T. S.
Ramakrishnan, D. H. Rothman, and L. M. Schwartz (1996), Transport insandstone: A study based on three dimensional microtomography, Geophys. Res. Lett., 23, 705 – 708, doi:10.1029/96GL00776.
Bhagat, Y. A., and C. Beaulieu (2004), Diffusion anisotropy in subcorticalwhite matter and cortical gray matter: Changes with aging and the role of
CSF-suppression, J. Magn. Reson. Imaging, 20, 216 – 227, doi:10.1002/
jmri.20102.
Crank, J. (1975), The Mathematics of Diffusion, 2nd ed., Oxford Univ.
Press, Oxford, U. K.
de Azevedo, E. N., M. Engelsberg, J. O. Fossum, and R. E. de Souza(2007), Anisotropic water diffusion in nematic self-assemblies of clay
nanoplatelets suspended in water, Langmuir, 23, 5100 – 5105,
doi:10.1021/la0632629.
Du, D. X., D. X. Zitha, and M. G. H. Uijttenhout (2007), Carbon dioxide
foam rheology in porous media: A CT scan study, SPE J., 12, 245 – 252,
doi:10.2118/97552-PA.
Farver, J. R., and R. A. Yund (1999), Oxygen bulk diffusion measurements
and TEM characterization of a natural ultramylonite: Implications for
fluid transport in mica-bearing rocks, J. Metamorph. Geol., 17,669 – 683, doi:10.1046/j.1525-1314.1999.00223.x.

18 of 19

W12435

NAKASHIMA ET AL.: DIFFUSION ELLIPSOIDS OF ANISOTROPIC POROUS ROCKS

Gherardi, F., T. F. Xu, and K. Pruess (2007), Numerical modeling of selflimiting and self-enhancing caprock alteration induced by CO2 storage ina depleted gas reservoir, Chem. Geol., 244, 103 – 129, doi:10.1016/
j.chemgeo.2007.06.009.
Goldstein, L., S. O. Prasher, and S. Ghoshal (2007), Three-dimensional
visualization and quantification of non-aqueous phase liquid volumes
in natural porous media using a medical X-ray computed tomographyscanner, J. Contam. Hydrol., 93, 96 – 110, doi:10.1016/j.jconhyd.
2007.01.013.
Guéguen, Y., and V. Palciauskas (1994), Introduction to the Physics ofRocks, Princeton Univ. Press, Princeton, N. J.
Ikeda, S., T. Nakano, and Y. Nakashima (2000), Three-dimensional study
on the interconnection and shape of crystals in a graphic granite by X-ray
CT and image analysis, Min. Mag., 64, 945 – 959, doi:10.1180/
002646100549760.
Jones, K. W., H. Feng, S. Tomov, W. J. Winters, M. Prodanovic, and
D. Mahajan (2007), Characterization of methane hydrate host sediments
using synchrotron-computed microtomography (CMT), J. Petrol. Sci.
Eng., 56, 136 – 145, doi:10.1016/j.petrol.2006.03.029.
Karacan, C. O., and G. D. Mitchell (2003), Behavior and effect of different
coal microlithotypes during gas transport for carbon dioxide sequestration into coal seams, Int. J. Coal Geol., 53, 201 – 217, doi:10.1016/
S0166-5162(03)00030-2.
Ketcham, R. A., and G. J. Iturrino (2005), Nondestructive high-resolution
visualization and measurement of anisotropic effective porosity in complex lithologies using high-resolution X-ray computed tomography,
J. Hydrol. Amsterdam, 302, 92 – 106, doi:10.1016/j.jhydrol.2004.06.037.
Lin, W., and M. Takahashi (2000), Mercury intrusion porosimetry and its
application for determination of pore volume distribution of rock (inJapanese), Chishitsu Nyusu, 549, 61 – 68.
Lindquist, W. B., A. Venkatarangan, J. Dunsmuir, and T.-F. Wong (2000),
Pore and throat size distributions measured from synchrotron X-ray tomographic images of Fontainbleau sandstone, J. Geophys. Res., 105,
21,509 – 21,527, doi:10.1029/2000JB900208.
Liu, J. S., A. Polak, D. Elsworth, and A. Grader (2005), Dissolution-induced preferential flow in a limestone fracture, J. Contam. Hydrol., 78,
53 – 70, doi:10.1016/j.jconhyd.2005.03.001.
Madsen, F. T. (1998), Clay mineralogical investigations related to nuclear
waste disposal, Clay Miner., 33, 109 – 129, doi:10.1180/000985598545318.
Matsumoto, M., and T. Nishimura (1998), Mersenne Twister: A 623dimensionally equidistributed uniform pseudorandom number generator,
ACM Trans. Model. Comput. Simul., 8, 3 – 30, doi:10.1145/272991.
272995.
Nakashima, Y. (2000), The use of X-ray CT to measure diffusion coefficients of heavy ions in water-saturated porous media, Eng. Geol., 56,
11 – 17, doi:10.1016/S0013-7952(99)00130-1.
Nakashima, Y., and S. Kamiya (2007), Mathematica programs for the analysis of three-dimensional pore connectivity and anisotropic tortuosity ofporous rocks using X-ray computed tomography image data, J. Nucl. Sci.
Technol., 44, 1233 – 1247, doi:10.3327/jnst.44.1233.
Nakashima, Y., and T. Kikuchi (2007), Estimation of the apertures of watersaturated fractures by nuclear magnetic resonance well logging, Geophys.
Prospect., 55, 235 – 254, doi:10.1111/j.1365-2478.2006.00596.x.
Nakashima, Y., and Y. Watanabe (2002), Estimate of transport properties ofporous media by micro-focus X-ray computed tomography and random
walk simulation, Water Resour. Res., 38(12), 1272, doi:10.1029/
2001WR000937.
Nakashima, Y., T. Nakano, K. Nakamura, K. Uesugi, A. Tsuchiyama, andS. Ikeda (2004), Three-dimensional diffusion of non-sorbing species inporous sandstone: Computer simulation based on X-ray microtomography using synchrotron radiation, J. Contam. Hydrol., 74, 253 – 264,
doi:10.1016/j.jconhyd.2004.03.002.
Oh, W., and W. B. Lindquist (1999), Image thresholding by indicator kriging, IEEE Trans. Pattern Anal. Mach. Intell., 21, 590 – 602,
doi:10.1109/34.777370.

W12435

Ohkubo, T. (2008), Tortuosity based on anisotropic diffusion process in
structured plate-like obstacles by Monte Carlo simulation, Transp. PorousMedia, 72, 339 – 350, doi:10.1007/s11242-007-9153-0.
Ohlsson, Y., and I. Neretnieks (1995), Literature survey of matrix diffusion
theory and of experiments and data including natural analogues, SKB
Tech. Rep. TR 95 – 12, Swed. Nucl. Fuel and Waste Manage. Co., 1-89.
Ohtani, T., T. Nakano, Y. Nakashima, and H. Muraoka (2001), Threedimensional shape analysis of miarolitic cavities and enclaves in the
Kakkonda granite by X-ray computed tomography, J. Struct. Geol., 23,
1741 – 1751, doi:10.1016/S0191-8141(01)00033-5.
Polak, A., A. S. Grader, R. Wallach, and R. Nativ (2003), Chemical diffusion between a fracture and the surrounding matrix: Measurement bycomputed tomography and modeling, Water Resour. Res., 39(4), 1106,
doi:10.1029/2001WR000813.
Porion, P., S. Rodts, M. Al-Mukhtar, A. M. Faugére, and A. Delville(2001), Anisotropy of the solvent self-diffusion tensor as a probe ofnematic ordering within dispersions of nanocomposites, Phys. Rev. Lett.,
87, 208302, doi:10.1103/PhysRevLett.87.208302.
Porion, P., M. Al-Mukhtar, A. M. Faugére, R. J. M. Pellenq, S. Meyer, and
A. Delville (2003), Water self-diffusion within nematic dispersions of
nanocomposites: A multiscale analysis of 1H pulsed gradient spin-echo
NMR measurements, J. Phys. Chem. B, 107, 4012 – 4023, doi:10.1021/
jp022161q.
Promentilla, M. A. B., T. Sugiyama, T. Hitomi, and N. Takeda (2008),
Characterizing the 3D pore structure of hardened cement paste withsynchrotron microtomography, J. Adv. Concr. Technol., 6, 273 – 286,
doi:10.3151/jact.6.273.
Schwartz, L. M., F. Auzerais, J. Dunsmuir, N. Martys, D. P. Bentz, and S.
Torquato (1994), Transport and diffusion in three-dimensional composite
media, Physica A, 207, 28 – 36, doi:10.1016/0378-4371(94)90351-4.
Sorai, M., T. Ohsumi, M. Ishikawa, and K. Tsukamoto (2007), Feldspar
dissolution rates measured using phase-shift interferometry: Implications
to CO2 underground sequestration, Appl. Geochem., 22, 2795 – 2809,
doi:10.1016/j.apgeochem.2007.06.019.
Uesugi, K., Y. Suzuki, N. Yagi, A. Tsuchiyama, and T. Nakano (2001),
Development of high spatial resolution X-ray CT system at BL47XU in
SPring-8, Nucl. Instrum. Methods Phys. Res. Sect. A, 467 – 468, 853 –
856, doi:10.1016/S0168-9002(01)00491-0.
Van Loon, L. R., J. M. Soler, W. Müller, and M. H. Bradbury (2004),
Anisotropic diffusion in layered argillaceous rocks: A case study with
Opalinus Clay, Environ. Sci. Technol., 38, 5721 – 5728, doi:10.1021/
es049937g.
Westin, C. F., S. E. Maier, H. Mamata, A. Nabavi, F. A. Jolesz, andR. Kikinis (2002), Processing and visualization for diffusion tensor MRI,
Med. Image Anal., 6, 93 – 108, doi:10.1016/S1361-8415(02)00053-1.
Xue, Z., Y. Ishijima, and M. Takahashi (1992), Permeability and microgeometry of sandstones under hydrostatic pressure (in Japanese withEnglish abstract), Shigen to Sozai, 108, 769 – 775.
Yamamoto, M., K. Oki, and S. Hayasaka (1978), On the Ito and the
Yoshino pyroclastic flows in Kagoshima prefecture, Kagoshima Daigaku
Rigakubu Kiyo, 11, 121 – 127.
Yokoyama, T., and S. Nakashima (2005), Diffusivity anisotropy in a rhyolite and its relation to pore structure, Eng. Geol., 80, 328 – 335,
doi:10.1016/j.enggeo.2005.07.002.
Zhu, W. C., J. Liu, D. Elsworth, A. Polak, A. Grader, J. C. Sheng, and J. X.
Liu (2007), Tracer transport in a fractured chalk: X-ray CT characterization and digital-image-based (DIB) simulation, Transp. Porous Media,
70, 25 – 42, doi:10.1007/s11242-006-9080-5.



S. Kamiya, T. Nakano, and Y. Nakashima, National Institute of Advanced
Industrial Science and Technology, Geological Survey of Japan, Central 7,
Higashi 1-1-1, Tsukuba 305-8567, Japan. (nakashima.yoshito@aist.go.jp)

19 of 19