nmr petrophysical predictions on digitized core images - ANU

1 downloads 0 Views 1MB Size Report
Jun 29, 2005 - 2School of Petroleum Engineering, University of New South Wales, Sydney, Australia ..... tion (Lawson and Hansen, 1974). For the inversion of.
SPWLA 46th Annual Logging Symposium, June 26-29, 2005

NMR PETROPHYSICAL PREDICTIONS ON DIGITIZED CORE IMAGES C. H. Arns1,* , A.P. Sheppard1 , R. M. Sok1 , and M.A. Knackstedt1,2 1 Department of Applied Mathematics, Research School of Physical Sciences and Engineering, Australian National University, Canberra, Australia 2 School of Petroleum Engineering, University of New South Wales, Sydney, Australia * Corresponding Author: [email protected] In a preliminary study we compare predictions of petrophysical properties from the interpretation of the NMR response to direct calculations on the images. The foundational assumption of pore isolation is directly tested by partitioning of the pore space. This allows one to calculate the coupling constants and magnetisation exchange between pores, or between macro- and micro-porous regions. Further, the sensitivity of the responses to variations in relaxivities or the presence of magnetic impurities is studied. Fluid typing is performed on a shaly sandstone sample.

Copyright 2005, held jointly by the Society of Petrophysicists and Well Log Analysts (SPWLA) and the submitting authors. This paper was prepared for presentation at the SPWLA 46th Annual Logging Symposium held in New Orleans, Louisiana, United States, June 26-29, 2005.

ABSTRACT NMR is a popular logging technique used to estimate pore size information, formation permeability, wettability and irreducible water saturation. Quantitative interpretation of NMR data is based on a set of fundamental assumptions (e.g., pore isolation and fast diffusion). These assumptions establish the quantitative link between NMR response and petrophysical predictions. While there is a need to test these assumptions directly, to date no quantitative study on reservoir core material has been undertaken. The ability to digitally image reservoir rock in 3D, calculate petrophysical properties directly from the images coupled with a comprehensive simulation tool to numerically generate a range of NMR response data may help to address this need.

MMM

INTRODUCTION NMR techniques are usually employed in the petroleum industry to either predict permeability or for fluid typing. The former application uses the surface relaxation mechanism or internal fields to derive a length scale (Brownstein and Tarr, 1979; Kenyon et al., 1986; Kenyon et al., 1988; Kenyon, 1992; Song et al., 2000), which can then be used in permeability correlations (Banavar and Schwartz, 1987; Sen et al., 1990). The method relies on the application of the fast diffusion limit, in which the magnetisation of an isolated pore decays over time as a single exponential: (Wayne and Cotts, 1966; Brownstein and Tarr, 1979)   t , (1) M (t) = M0 (t) exp − T2

In this paper we image a large set of reservoir cores including sandstones and carbonates at the pore scale using high resolution micro-CT. A set of petrophysical properties are measured directly on the cores including surfaceto-volume, permeability and pore size distribution. The permeabilities of the cores range from 10 mD to several Darcies. Realistic multiphase fluid distributions are derived by simulation of drainage. We then simulate the NMR responses on the same core images using a comprehensive NMR simulator. The internal magnetic field is derived numerically from applied magnetic fields and susceptibility distributions and the phase evolution of the magnetic spins calculated with a random walk method. NMR responses currently include inversion recovery (T1) and CPMG (T2), and the longitudinal and transversal signals are monitored simultaneously. The interpretation of the signals acquired is done by standard 1D Laplace inversion to calculate the pore size distribution from T2 responses, and with a 2D inverse Laplace transform for fluid typing.

where M0 is the initial magnetization and the transverse relaxation time T2 is given by 1 1 Sp = +ρ , T2 T2b Vp

(2)

with Sp /Vp the surface-to-pore-volume ratio of the pore space, T2b the bulk relaxation time of the fluid that fills the pore space and ρ the surface relaxation strength. For small pores or large ρ the bulk relaxation contribution is considered negligible and 1 1 Sp = =ρ . T2 T2s Vp

1

(3)

SPWLA 46th Annual Logging Symposium, June 26-29, 2005

Real rocks contain a network of pores of different sizes connected by pore throats which restrict diffusion. If interpore diffusion is considered negligible, each pore can be considered distinct and the magnetization within independent pores assumed to decay independently. The decay can then be described as N X

 t , M (t) = M0 (t) ai exp − T2i i=1 

2004; Arns et al., 2005). It has advantages in being able to study petrophysical parameters under very controlled conditions in realistic heterogeneous structures. In this paper we numerically predict NMR responses and a range of petrophysical properties on a suite of 3D images of reservoir rock derived via micro-CT (Arns et al., 2005). This allows us to directly address assumptions used in NMR interpretation and to test common correlations between response and permeability and fluid-typing. The paper is organised as follows: firstly, the experimental setup and datasets underlying this study are briefly introduced. Secondly, the petrophysical calculations are detailed. The result section then presents a preliminary analysis of NMR data on realistic microstructures including studies of;

(4)

where ai is the volume fraction of a pore i decaying with relaxation time T2i . The multi-exponential distribution corresponds to a partition of the pore space into N groups based on the Sp /Vp values of the pores. Analog expressions exist for deriving length scales using higher diffusion eigenmodes (Song et al., 2000; Song, 2003). Given a pore, grain and fluid partitioning, the fast diffusion (Eqn. 3-4) limit can be adapted to varying surface relaxivity and different fluids by accounting for the number of active surfaces N in each pore N 1 X 1 ρ i Si . = T2 N Vp i=1

• length scales used in permeability - relaxation time correlations, showing that prefactors in these correlations can vary by an order of magnitude depending on the sample structure, • the sensitivity of the length scale prediction to heterogeneity in surface relaxivities, which for the sample discussed results in a factor of two difference in predicted permeability,

(5)

In the case of partially saturated pores one would have to treat each fluid partition of a pore separately. However, the information needed to calculate this fast diffusion limit is prohibitive and generally not available.

• the sensitivity of the length scale prediction to diffusion coupling with clay regions, affecting permeability by up to 45% for the considered case,

It is widely recognised that the simple picture of magnetisation decay is not always applicable to reservoir rocks. The length scale estimate can be affected by surface relaxivity heterogeneities(Øren et al., 2002), internal fields (H¨urlimann, 1998; Sen and Axelrod, 1999; Appel et al., 2001; Dunn, 2001; Dunn, 2002), by being outside the fast diffusion limit (Zhang and Hirasaki, 2003), or by being outside the weak coupling regime (Cohen and Mendelson, 1982; McCall et al., 1991; McCall et al., 1993; Zielinski et al., 2002; Toumelin et al., 2003). The correlation to permeability also implies a correlation between NMR response and pore size; this correlation has recently been questioned (Lonnes et al., 2003). Most studies are either experimental or based on simulations in relatively simple geometries.

• partial saturation effects, showing a non-linear relationship of water saturation to log mean relaxation time, • a detailed analysis of pore-pore coupling using a pore partitioning and the concept of topological distance, showing that all rocks considered here are at least in the weak coupling regime. The overall effect on predicted permeability is small for clastics, but more problematic for carbonates. • fluid typing on a shaly sandstone sample at different partial saturations (water/oil). MATERIALS AND METHODS This section describes the samples used in this study, the acquisition of images and the computational techniques used to calculate their physical properties.

The second major area of application of NMR in reservoir rocks is fluid typing (see e.g. (Looyestijn, 1996; H¨urlimann et al., 2003)). This application has been made much more reliable by using 2D NMR techniques. However limitations exist in the Laplace inversion technique required (Sun and Dunn, 2004).

Image acquisition A set of 20 sandstones, including five cores from an offshore gas well, four cores from an offshore oil well, four poorly consolidated cores from a prospective oil reservoir, a thinly bedded sandstone, a shaly sandstone, four Fontainebleau sandstone samples, and one Berea sandstone sample were considered as well as two carbonates,

In recent years a new methodology has evolved combining X-ray µCT imaging techniques with numerical calculations of petrophysical properties (Auzerais et al., 1996; Knackstedt et al., 2004; Arns et al., 2004a; Arns,

2

SPWLA 46th Annual Logging Symposium, June 26-29, 2005

a vuggy carbonate and an outcrop limestone. All have been imaged via high-resolution X-ray µCT (Sakellariou et al., 2004b; Sakellariou et al., 2004a). Details of the experimental methodology and most of the numerical procedures has been described in previous papers (Arns et al., 2001; Arns et al., 2002; Arns et al., 2004b; Knackstedt et al., 2004; Arns et al., 2005) as has the permeability analysis for the Fontainebleau (Arns et al., 2004b) and the oil and gas reservoir samples (Arns et al., 2004a). Most images are recorded at a resolution ∼ 5 µm with a field of view of ∼ 8 mm. Here we give information for those samples, which are considered for most extensive analysis in the following sections. Table 1 gives the dimensions of these samples and Fig. 1 slices through the raw tomograms. Sample A is a clean sandstone (Fontainebleau), sample B a shaly sandstone, sample C a vuggy multi-scale carbonate, and sample D another highly porous carbonate (limestone) from an outcrop in South Australia showing large well resolved features.

[a]

[b]

[c]

[d]

Image segmentation and partitioning

Figure 1: Slices through the raw tomograms (2048 × 2048 voxel): [a] Fontainebleau sandstone (segmented, 4803 ), [b] shaly sandstone, [c] vuggy carbonate, [d] feature rich limestone.

Two-phase segmentation: To carry out quantitative analysis on tomographic images one has to phase separate the pore space, and possibly various mineral phases. This is not trivial as the X-ray density does not show a clear distinction between phases and simple thresholding will fail. Here, phase separation into two phases was achieved using a multi-stage procedure (Sheppard et al., 2004).

Table 1: Basic description of the tomographic images in terms of sample type, resolution [µm], segmented and analysed image subsection, and pore volume fraction. Sandstone B is a shaly sandstone with an unresolved clay/silt fraction of 4.83%. Sample Sandstone A Sandstone B Carbonate C Limestone D

Three-phase segmentation: Partitioning of a dataset into three phases (void, clay and solid) is performed on the Xray density data using a modified version of the converging active contours method outlined in (Sheppard et al., 2004). First we identify voxels that can easily be classified into one of the three phases. Any voxel whose gradient value is above a threshold tg is considered to be near an interface and therefore is classified as “undecided” by the initial thresholding. The other voxels are classified as follows: if the voxel intensity is less than t0 then it is considered void, if it is between t1 and t2 , then it is classified clay, while voxels with intensity greater than t3 are set to solid. The remaining voxels are undecided. The voxels classified as void, clay or solid are used as seed regions for the converging active contours method. An example of a three-phase partitioning is given in Fig. 2. Pore partitioning: For the purposes of calculating the coupling constants measured using NMR, we need to partition the pore space into simple geometric cells (pore bodies), separated by narrow constrictions (throats). It is possible to do this by performing a watershed transformation on the Euclidean distance data, but this approach leads to a partitioning, whit topological properties that cannot be controlled.

[a]

Res. 5.68 5.60 1.08 3.02

Size [voxel] 480 × 480 × 480 944 × 1208 × 1480 960 × 960 × 1920 900 × 900 × 1800

φ 17.7 8.93 15.6 50.7

[b]

Figure 2: Illustration of the three-phase segmentation on a slice of the shaly sandstone sample (A). [a] raw tomographic data, [b] phase separated into 3 phases; black (grain), white (pore) and grey (clay).

3

MMM

SPWLA 46th Annual Logging Symposium, June 26-29, 2005

Figure 4: Illustration of the invasion of a non-wetting fluid into the water saturated shaly sandstone sample A. Colours: the solid grain phase is black, the clay region dark grey, the defending fluid light grey, and the invading fluid white. ing a Voronoi tesselation starting from the two medial axis pieces to partition the other voxels in the link. An illustration of the pore partitioning resulting is given in Fig. 3.

Figure 3: Illustration of the pore partitioning on a 2003 subsection of the pore network of sample B.

Drainage simulation Multi-phase fluid distributions are generated by drainage simulations directly on the voxelated image. A fixed capillary pressure is associated with a pore entry radius. Invading the defending fluid by a sphere of that radius from all boundaries mirrors standard mercury intrusion boundary conditions. The center of the invading sphere is allowed to move such that the sphere does not overlap the solid(Hilpert and Miller, 2001). An example of a resulting fluid distribution is given in Fig. 4.

Our approach is significantly more involved. We commence with Euclidean distance data, and perform distanceordered homotopic thinning to determine the voxels that comprise the skeleton, or medial axis, of the pore space. The medial axis defines a network, in which nodes, clusters of highly connected voxels, are connected through links composed of 2-connected voxels. We analyse each link in the network, calculating its length to width ratio, and determining the number and magnitude of the constrictions along its length. If the link is short and represents a minor constriction then the nodes at each end of the link may be merged. On the other hand, if the link consists of a number of distinct constrictions, additional nodes may be inserted along its length. Merging can only occur if it doesn’t alter the topology of the network; if two nodes are already connected to each other by a single link, then merging them will collapse a ring of the network. Node merging and link division is of critical importance in identifying representative pore bodies.

Formation factor The conductivity calculation is based on a solution of the Laplace equation with charge conservation boundary conditions and has been detailed before (Arns et al., 2001). We assign to the matrix phase of the sandstone a conductivity σm = 0 and to the (fluid-filled) pore phase a normalized conductivity σfl = 1. A potential gradient is applied in each coordinate direction, and the system relaxed using a conjugate gradient technique to evaluate the field. The formation factor given by F = σfl /σeff , is used.

Having determined the topology of the network of pore bodies, we now partition the pore space, i.e. determine to which pore body each voxel belongs. This is achieved by first performing a Voronoi tessellation of the pore space, using the link voxel strings as seeds. This associates each voxel in the pore space with the link that it’s closest to, and is done using Dijkstra-Sethian fast marching (Sethian, 1999). Each link region is then divided into two pieces, with the minimum constriction in the link separating the pieces. We do this by breaking the medial axis voxels in the link at the minimum constriction, then us-

Permeability calculation Permeability is calculated using the mesoscopic latticeBoltzmann method (LB) (Martys and Chen, 1996; Qian and Zhou, 1998). It can be shown, that the macroscopic dynamics of the solution of a discretized Boltzmann equation match the Navier-Stokes equation. Due to its simplicity in form and adaptability to complex flow geometries one of the most successful applications of the LB

4

SPWLA 46th Annual Logging Symposium, June 26-29, 2005

method has been to flow in porous media (Chen et al., 1992; Frisch et al., 1986; Rothman, 1988; Ferreol and Rothman, 1995; Martys and Chen, 1996). In this study we applied a pressure gradient by a body force (Ferreol and Rothman, 1995), used closed boundary conditions perpendicular to the flow and mirrored boundaries parallel to the pressure gradient, resulting in a system size of L × L × 2L. Permeability was measured over the L × L × L original image of the simulated system.

0

M(t,g)/M0

first echo intensity [a.u.]

0.10

[a]

−2

10 0

0.0005

0.001

0.0015

tE [s] 0

[b]

10

−1

10

0

10

−1

10

0

M(t,g) / M0

M(t,g)/M0

−1

0.2 T/m 0.5 T/m 1.0 T/m 2.0 T/m 5 T/m 10 T/m Neuman

0.5 T/m 1.0 T/m 2.0 T/m 5.0 T/m 10 T/m Neuman

−1

10

c1=108/175 c2=747/20

−2

[c]

−2

10

10

10

10

−3

10

tE [s]

10

Surface relaxation: The spin relaxation of a saturated porous system is simulated by using a lattice random walk method (Mendelson, 1990; Bergman et al., 1995; Valckenborg et al., 2002). Initially the walkers are placed randomly in the 3D pore space. At each time step the walkers are moved from their initial position to a neighboring site and the clock of the walker advanced by ∆t = 2 /(6D0 ), where  is the lattice spacing and D0 the bulk diffusion constant of the fluid, reflecting Brownian dynamics. The lattice is made periodic by mirroring the structure in all directions. An attempt to go to a site of another phase will kill the walker with probability ν/6, 0 ≤ ν ≤ 1 (Mendelson, 1990). The killing probability ν is related to the surface relaxivity ρ via    ρ 2 ρ +O , (6) Aν = D0 D

0.5 T/m 1.0 T/m 2.0 T/m 5.0 T/m 10 T/m Robertson

−1

10

Simulation Free diffusion eqn. 0.01

NMR Response Simulation

−2

−3

10

−2

−1

10

10

tE [s]

10

0

10

[d]

−3

10

−2

10

0

tE [s]

Figure 5: Comparison of analytical results and simulations for the decay of the spin-echo intensity under [a] free diffusion and [b-d] restricted diffusion inside a ndimensional sphere for various field gradients [b] line, [c] circle, and [d] sphere.

1966; Neuman, 1974; Sukstanskii and Yablonskiy, 2002) 2 −2 4 2 rdg )] , (9) rdg (1 − c2 rsg M (2τ, g) = M0 exp[−c1 rsg

where rsg is the ratio lS /lg , rdg is the ratio lD /lg , and c1 and c2 are geometrical constants. For spheres c1 = 108/175, c2 = 747/20, for circles c1 = 7/36 and c2 = 297/14, and for the planar case c1 = 1/720 and c2 = 51/28. The numerical code is tested by comparison to these analytical results in Fig. 5, using the diffusion constant of water (D0 = 2500 µm2 /s), and d = 2 µm, where d is the diameter of the n-dimensional sphere. The agreement is excellent (10000 random walkers were used for the simulation).

where A is a correction factor of order 1 (here, we take A = 3/2) accounting for the details of the random walk implementation (Bergman et al., 1995). The walkers initially start with zero phase and accumulate phase according to their path. Since in a real experiment one would demodulate the signal with the Lamor frequency of the static applied field, the phase accumulation per step is ∆φ = γ∆t(B(~r) − B0 ),

10

1.00

(7)

Analysis of diffusion coupling: The delineation of the individual pores in some rocks via pore partitioning allows tracking of the interpore movements of the random walkers within the pore space and enables us to calculate the number of crossings between adjacent pores at any Dephasing and code validation: The decay of the spinone time step. To minimise the sensitivity to resolution echo amplitude in porous media is governed by a set of and step size of the random walk, we record a number length scales, in particular the structural length lS = of parameters. First we log the number and direction a V /S, the dephasing length lg = (D0 /(γg))1/3 (with g walker crossed a pore-pore “interface”. We also have a being √ the field gradient strength), and the diffusion length second counter, in which we only record interface crosslD = D0 τ . At short times in the free diffusion regime, ings, which do not reverse the last interface crossing. A where lD is the shortest length and using the pulse sethird counter records the starting pore and end pore of a quence π2 − te − π − te − signal, the spin-echo decay is random walk. For this application we start a fixed numgiven by ber of random walkers at every pore voxel, reflecting the 1 hydrogen index of the fluid present. A fourth counter M (2τ, g) = M0 exp[ γ 2 g 2 D0 (2τ )3 ]. (8) 12 records the mean first passage time across the pore-pore At longer diffusion times, in the motional averaging regime, interfaces. where lS is the shortest length scale, one has (Robertson,

where γ is the gyromagnetic ratio. It is assumed, that any spin-lattice relaxation also destroys phase coherency, thus T2 ≤ T1 .

5

MMM

SPWLA 46th Annual Logging Symposium, June 26-29, 2005

Inverse Laplace transform: The relaxation time distribution is obtained by fitting a multi-exponential decay to M (t) or M (t1 , t2 ). Both the one- and the two-dimensional inversion routines use a bounded least squares solver (Stark and Parker, 1995) combined with Tikhonov regularisation (Lawson and Hansen, 1974). For the inversion of one-dimensional decay curves the L-curve method is used (Hansen, 1992) to find the optimal regularisation parameter. The 2D inversion routine follows closely the implementation of (Venkataramanan et al., 2002) with kernel compression to reduce memory requirements and computing time. Permeability correlation Permeability correlations are usually based on the logarithmic mean T2lm of the relaxation time   P i log(T2i ) i aP , (10) T2lm = exp − ai

which is assumed to be related to an average Vp /Sp or pore size. Commonly used NMR response/permeability correlations include the porosity φ as in (Banavar and Schwartz, 1987; Kenyon et al., 1988) c k = aφb T2lm ,

(11)

with classical factors a = 1, b = 4, c = 2, or the Formation factor F as in c k = aF b T2lm ,

(12)

with standard factors b = −1, c = 2.

The use of c = 2 in Eqns. 11-12 implies a unit of a as surface relaxivity squared. In our fits of the permeability correlations we scale the value of a by 1/(6ρ)2 to make the prefactor dimensionless; this implies that any difference in the prefactor arises only from structural influences. The length scale associated with T2 (the pore size derived from the NMR signal is given by dT 2lm = 6T2lm ρ.

RESULTS NMR Permeability Correlations Structural influences: In these simulations simple surface relaxation was assumed with constant ρ = 16 µm/s for all sandstone samples and no bulk relaxation. For the carbonate samples C and D we used ρ = 5 µm/s. We illustrate data obtained for two sets of samples in Fig. 6 – the gas reservoir samples and the poorly consolidated core. The gas reservoir cores were part of the study from (Arns et al., 2004a) which exhibited a broad range of permeability (varying over three orders of magnitude) across a small range of φ (15% < φ < 25%). While scatter in the data is observed the correlation is quite good. The poorly consolidated sands all have Darcy range permeability and the match is excellent across all values. In Table 2 we summarise the NMR relaxation response correlations with permeability, based on the log-mean of the relaxation time distribution. We immediately note that the use of (Eqn. 12) dramatically increases the correlation coefficient compared to the use of Eqn. 11. This shows the important role a measure of F can make to more accurate permeability estimation. The prefactors in Eqns. 11 and 12 are quite consistent for samples from the same reservoir (particularly the poorly consolidated cores B2-B5). The variation in the prefactor is quite small for all sandstone samples for both correlations. The prefactor of Eqn. 12 is similar for the two carbonates. In contrast the best prefactor for Eqn. 11 is almost an order of magnitude smaller than observed for the sandstone samples. Variability of ρ: We choose sample B to test the sensitivity of the permeability correlations to varying ρ by making use of the three-phase partitioning of the sample into grain phase, clay regions, and pore space. Surface relaxivities are set separately for the grain and clay region surfaces (resolution is 5.6 µm) with ratios 1:1, 1:4, and 1:10 in such a way, that the mean surface relaxivity stays constant. Three sets of simulations are run with hρi ∈ {4, 7.9, 16} µm/s. No bulk relaxation term is included in this simulation, so even random walkers starting far away from a surface will relax eventually on the surface. Accordingly, T2lm is caused by surface relaxation only and can be comparably long. The results are summarised in Table 3 and Fig. 7.

We have previously shown that one can obtain useful estimates of petrophysical properties from simulations at the scale of a few mm3 (Arns et al., 2004a; Arns et al., 2005). Each sample is divided into subregions and for each subregion of all samples the permeability and NMR surface relaxation response is calculated. This means that one obtains 100 individual samples per image and therefore a relationship between k, φ, T2lm for each rock. In all fits the mean residual error P 2 (log10 (kcalc ) − log10 (kemp )) . (13) s2 = n−2 is minimised and the correlation coefficient P (kemp − kemp )(kcalc − kcalc ) R = P 1/2 (14) (kemp − kemp )2 (kcalc − kcalc )2

As can be seen in Table 3, the derived length scale for permeability correlations increases both with higher mean surface relaxivity and with higher heterogeneity of the surface relaxivity at constant mean surface relaxivity. A reason for this could be due to the heterogeneous clay distribution (Fig. 2). The clay is compartmentalized and the regions of different surface relaxivity do not communicate with each other. This may prevent averaging as

calculated.

6

SPWLA 46th Annual Logging Symposium, June 26-29, 2005

3

10

kcalc [mD]

Table 2: NMR response correlations to permeability. The sandstone samples gas 1-4, oil 1-5 and B2-B5 come from three different fields. Added are the standard benchmark rocks Berea and Fontainebleau (Fb, 4 samples), as well as the four samples featured in Fig. 1 and Table 1. Samples C and D are carbonates. The prefactor a, mean residual error s2 and correlation coefficient R have indices. “1” refers to Eqn. 11 and “2” refers to Eqn. 12.

2

10

1

10

A B C D oil 1 oil 2 oil 3 oil 4 gas 1 gas 2 gas 3 gas 4 gas 5 B2 B3 B4 B5 Berea Fb

a1 36ρ2

.088 .114 .029 .016 .252 .386 .301 .045 .096 .036 .177 .134 .143 .088 .091 .078 .115 .152 .145

s21

R1

.273 .191 .357 .066 .181 .146 .097 .833 .191 .019 .054 .034 .010 .0091 .0019 .0044 .0047 .0141 .0818

.772 .875 .056 .510 .622 .791 .797 .727 .521 .315 .728 .701 .725 .524 .681 .444 .458 .804 .846

100a2 36ρ2

.467 .312 .335 .657 .238 .242 .286 .048 .171 .059 .468 .360 .511 .551 .528 .556 .515 1.06 .302

s22

R2

.172 .148 .272 .0293 .221 .0767 .0814 .701 .162 .018 .029 .015 .0024 .0010 .0004 .0005 .0009 .0056 .0897

.929 .962 .714 .770 .866 .960 .952 .948 .945 .600 .956 .973 .891 .961 .908 .959 .972 .919 .951

[a]

ρg [µm/s] 4.000 2.542 1.470 7.869 5.000 2.892 16.00 10.17 5.879

ρc [µm/s] 4.000 10.17 14.70 7.869 20.00 28.92 16.00 40.67 58.79

T2lm [s] 5.30 5.66 6.32 2.77 3.00 3.43 1.42 1.56 1.85

2

10

3

10

4

10

3

10

2

10

[b]

2

10

3

10

4

10

kemp [mD]

Figure 6: Illustration of the deviation of the individual permeability predictions for two sample subsets. The fits according to Eqn. 12 with parameter a2 (see Table 2) for [a] gas reservoir samples and [b] B2-B5. Note the large scatter around the fit for the gas samples (high s2 ) compared with the tight fit for samples B2-B5 (low s2 ).

required by the fast diffusion limit. The range of dT 2lm corresponds to a factor of about two in the permeability prediction.

Table 3: Sensitivity analysis of T2lm upon varying surface relaxivity ρ of the clay (ρc ) and grain surfaces (ρg ) with constant mean hρi for a shaly sandstone sample. The macroporous (void) fraction is .0894 and the fraction of the clay region .0483. The specific void:grain surface area is Sg = 5.64 mm−1 and the specific void:clayregion surface area Sc = 1.33 mm−1 (at image resolution). hρi [µm/s] 4.000 4.000 4.000 7.869 7.869 7.869 16.00 16.00 16.00

1

10

kemp [mD]

kcalc [mD]

Sample

Many clay regions are microporous. In that case spins can diffuse into the clay region and potentially relax much faster. We make a preliminary attempt to estimate this effect by adding a 10% background porosity and an effective background diffusion coefficient in and into the clay region of D0 /10. The results are given in Fig. 8 and Table 4. The effect of this background diffusion is generally a drop in the mean relaxation times, in particular for high contrast and for lower hρi, since the diffusion into the clay region adds surface area. This scenario makes clear, that in the presence of clay a good understanding of diffusion coupling with the clay region is needed to interpret the NMR response in terms of a length scale relevant for permeability correlations. In this particular case, the change in the predicted length scale dT 2lm is about 5-20%, implying a change of 10-45% in permeability. In current work we are extending this analysis to a wider range of clastics with larger clay fractions.

dT 2lm [µm] 127 136 152 130 142 162 136 150 178

7

MMM

SPWLA 46th Annual Logging Symposium, June 26-29, 2005

7

8

4

5

T2lm [s]

T2lm [s]

6

6

ρ = 4, 4/4 ρ = 4, 2.5/10 ρ = 4, 1.5/15 ρ = 8, 8/8 ρ = 8, 5/20 ρ = 8, 3/30 ρ = 16, 16/16 ρ = 16, 10/41 ρ = 16, 5.9/59

4 3

ρ = 4, 4/4 ρ = 4, 2.5/10 ρ = 4, 1.5/15 ρ = 8, 8/8 ρ = 8, 5/20 ρ = 8, 3/30 ρ = 16, 16/16 ρ = 16, 10/41 ρ = 16, 5.9/59

2 2 1

0 0.05

0.06

0.07

0.08

[a]

φ

0.09

0.1

0.11

0 0.05

0.12

160

0.08

φ

0.09

0.1

0.11

0.12

0.1

0.11

0.12

250 ρ = 4, 4/4 ρ = 4, 2.5/10 ρ = 4, 1.5/15 ρ = 8, 8/8 ρ = 8, 5/20 ρ = 8, 3/30 ρ = 16, 16/16 ρ = 16, 10/41 ρ = 16, 5.9/59

200

dT2lm [µm]

dT2lm [µm]

180

0.07

[a]

220 200

0.06

150

ρ = 4, 4/4 ρ = 4, 2.5/10 ρ = 4, 1.5/15 ρ = 8, 8/8 ρ = 8, 5/20 ρ = 8, 3/30 ρ = 16, 16/16 ρ = 16, 10/41 ρ = 16, 5.9/59

140

100 120 100 0.05

0.06

0.07

0.08

[b]

φ

0.09

0.1

0.11

50 0.05

0.12

[b]

Table 4: Sensitivity analysis of T2lm upon varying surface relaxivity ρ of the clay (ρc ) and grain surfaces (ρg ) with constant mean hρi. Same as Table 3, but with 10% background diffusion into the clay region (see text). ρg [µm/s] 4.000 1.470 7.869 2.892 16.00 5.879

ρc [µm/s] 4.000 14.70 7.869 28.92 16.00 58.79

T2lm [s] 4.64 5.26 2.50 3.03 1.44 1.71

0.07

0.08

φ

0.09

Figure 8: Sensitivity of [a] T2lm and [b] dT 2lm upon varying surface relaxivity ρ of the clay-region (ρc ) and grain surfaces (ρg ) compared to constant mean hρi for a shaly sandstone sample assuming a 10% background porosity and an effective background diffusion coefficient in and into the clay region of D0 /10.

Figure 7: Sensitivity of [a] T2lm and [b] dT 2lm upon varying surface relaxivity ρ of the clay-region (ρc ) and grain surfaces (ρg ) compared to constant mean hρi for a shaly sandstone sample. The legends note mean surface relaxivity hρi and the values for the individual surfaces (see also Table 3).

hρi [µm/s] 4.000 4.000 7.869 7.869 16.00 16.00

0.06

Partial saturation effects: We describe briefly a method to study the effect of partial saturation. We use the same sample as in the previous section. The analysis is limited to two global saturations (SW = 35, 60%) resulting from drainage simulations on the full image (sample B, 944×1208×1480 voxel). From that subsection we select a 2 × 3 × 3 matrix of 4003 subsamples. A bulk T2b of 3s was used for both oil and water and the diffusivities of the fluids were identical. Where oil was in contact with the water or clay phase, the surface relaxivity is set to zero. In Fig. 9 we plot the relaxation time distribution of the water phase for various saturations. The relationship is not linear, since increased restriction of the water phase

dT 2lm [µm] 111 126 120 145 138 164

8

SPWLA 46th Annual Logging Symposium, June 26-29, 2005

3

Pore-pore coupling To analyse the NMR response in terms of pore network characteristics, we first derive some intrinsic characteristics of the network. Namely, the number of pores Np , their average coordination number hCi, and for each pair of pores the minimal number of links (throats) needed to traverse from one pore to the other. This latter quantity is called the topological distance dt . The average hdt i of this Np × Np matrix of topological distances is an intrinsic characteristic of the network and dependent on the size of the lattice, increasing with network size.

T2lm (water) [s]

2.5 2 1.5 1:10, SW=35% 1:1, SW=35% 1:10, SW=60% 1:1, SW=60% 1:10, SW=100% 1:1, SW=100%

1 0.5 0

0

0.02

0.04

0.06

φW [%]

0.08

0.1

We choose the four rocks shown in Fig. 1 to analyse the effect of diffusion coupling. Applying the pore partitioning, we track in the NMR simulations for each walker the pore in which it started, the first arrival times to other (ρ) pores, and the travelled topological distance dtt , which we define as the numbers of pores crossed (per walker) excluding crossings which go directly back to the pore visited before. We further count the total number of these crossings at a given topological distance and normalise by the total number of crossings (crossing intensity). This counting avoids noise from multiple crossings at porepore interfaces. While hdt i would increase indefinitely with increasing lattice size, the same is not true for hdtt i, since relaxation events or finite acquisition time limit the achievable topological distance. In Table 6 we report the (ρ) mean of the travelled topological distance hdtt i.. We

0.12

Figure 9: Sensitivity of T2lm upon partial saturation effects for sample B. A mean surface relaxivity of ρ = h16i µm/s and ratios of 1:1 and 1:10 for the relaxivity of the grain surfaces and clay-regions were used (see also Table 3). φW is the water volume fraction.

by the oil phase ( implying “smaller” pore size) competes with blockage of active surfaces by the oil. Table 5 lists the resulting log mean relaxation time of the water phase, and in a separate column the measured mean log relaxation time, taking the oil phase into account as well. For

Table 6: Characteristics of the pore partitioning used to analyse the pore-pore coupling strength influencing the NMR relaxation response. Reduced sample sizes were used, in particular a 4803 section of samples A,C,D and a 6003 section of sample B. Np notes the number of pores and hCi the average coordination number, and hdt i and hdtt i are topological length scales (unitless). The index to dtt indicates the surface relaxivity.

Table 5: Sensitivity analysis of T2lm upon varying surface relaxivity ρ of the clay (ρc ) and grain surfaces (ρg ) with constant mean hρi. Same as Table 3, but varying global water saturation SW (see text). T2w notes the log mean T2 of the relaxing water phase. hρi [µm/s] 16.00 16.00 16.00 16.00 16.00 16.00

ρg [µm/s] 16.00 5.879 16.00 5.879 16.00 5.879

ρc [µm/s] 16.00 58.79 16.00 58.79 16.00 58.79

SW [%] 100 100 59.9 59.9 35.1 35.1

T2w [s] .964 1.14 1.32 1.00 1.23 0.90

T2lm [s] .964 1.14 1.73 1.44 2.31 1.88

Sample A B C D

Np 4009 1210 2349 5349

hCi 4.71 3.51 4.46 5.75

hdt i 12.0 18.2 19.4 9.59

(3)

hdtt i 12.1 12.0 110 71

(16)

hdtt i 3.66 4.74 30.5 31.0

consider two different surface relaxivities, ρ = 3 µm/s, and ρ = 16 µm/s and no bulk relaxation. All isolated pores are removed and the coupling in the resulting network is analysed. Reduced sample sizes were used; a 4803 section of samples A,C,D and a 6003 section of sample B.

constant surface relaxivity the log mean relaxation time of the water phase rises above the 100% water saturated case, while for a 1:10 contrast between grain and clay region surface relaxivities the relationship is inverted. The log mean relaxation time of all diffusing spins, and including bulk relaxation, trends monotonically with water saturation. Heterogeneity in surface relaxivities at partial saturation has a more pronounced effect than in the fully water saturated case.

As expected the surface relaxivity has a large impact on the results. For ρ = 3 µm/s, a large number of pores are traversed by the diffusing spins. For the sands > 10 and

9

MMM

SPWLA 46th Annual Logging Symposium, June 26-29, 2005

for the carbonates ≈ 100. For ρ = 16 µm/s, the average travelled pore distance drops to less than four pores for sands, but remains ≈ 30 for the carbonates. This result leads to a questioning of the foundational assumption of pore isolation — all rocks considered here are at least in the weak coupling regime.

0.25

porosity [a.u.]

0.2

To investigate the effect of pore coupling on NMR responses, we ran two distinct sets of simulations. Coupling “on” is the normal case and walkers can travel freely in the pore space. In the case of coupling “off” we use the pore network to restrict diffusing spins to their respective pores; we isolate all pores based on the pore partitioning algorithm. Fig. 10 shows the resulting T2s distributions for the four rocks for two surface relaxivities. For all samples, the change in ρ causes a large shift in the relaxation time distribution. The log mean relaxation time (see Table 7) stays nearly constant with either coupling on or off. Sharpening of the spectrum occurs for the sandstones only at low surface relaxivity, while for the carbonates (C,D) sharpening of the spectrum is apparent also for high surface relaxivities (in particular Fig. 10d). There is a small but noticeable change in T2lm for sample D ( 9%). From this result we can conclude,

ρ = 3, coupling on ρ = 3, coupling off ρ = 16, coupling on ρ = 16, coupling off

0.15

0.1

0.05

0 2 10

3

10

10

4

T2s [ms]

[a] 0.2

porosity [a.u.]

0.15

ρ = 3, coupling on ρ = 3, coupling off ρ = 16, coupling on ρ = 16, coupling off

0.1

0.05

0 2 10

10

3

4

10

T2s [ms]

[b] 0.2

Table 7: Log mean relaxation time for samples A-D with and without coupling and for two different surface relaxivities (ρ = 3 µm/s, ρ = 16 µm/s). The indices c,n indicate coupling or no coupling, and the indices 3,16 the surface relaxivity. T2lm [s] 5.35 5.00 .981 5.09 (c3)

T2lm [s] 5.41 5.04 .998 5.65 (n3)

T2lm [s] 1.08 1.37 .189 1.77 (c16)

0.15

porosity [a.u.]

Sample A B C D

ρ = 3, coupling on ρ = 3, coupling off ρ = 16, coupling on ρ = 16, coupling off

T2lm [s] 1.08 1.38 .191 1.93 (n16)

0.1

0.05

0

10

100

1000

10000

T2s [ms]

[c] 0.6

that pore-pore coupling does not change the length scale used for permeability correlations if considering constant surface relaxivity, and therefore the permeability prediction is robust under these conditions. Testing this under various other conditions will be undertaken in future studies.

porosity [a.u.]

0.5

Since permeability is controlled by restrictions to flow rather than pore size, it is interesting to ask whether diffusing spins diffuse far enough during the measurement to experience these restrictions. To gain insight into this aspect, we plot the crossing intensity of diffusing spins over the topological distance in Fig. 11, and the first passage time (Kim and Torquato, 1990) over the topological distance in Fig. 12. For sandstone a maximum intensity of crossings is recorded for topological distances of only one pore, and for the number of pores crossed (travelled

ρ = 3, coupling on ρ = 3, coupling off ρ = 16, coupling on ρ = 16, coupling off

0.4 0.3 0.2 0.1 0 2 10

[d]

3

10

10

4

5

10

T2s [ms]

Figure 10: Relaxation time distributions for samples AD ([a]-[d]). Pore-pore coupling effects can be seen for ρ = 3 µm/s, and for the Carbonate samples also for ρ = 16 µm/s.

10

SPWLA 46th Annual Logging Symposium, June 26-29, 2005

5 ρ = 3, topological distance ρ = 3, travelled distance ρ = 16, topological distance ρ = 16, travelled distance

0.6

First passage time [s]

crossing intensity [a.u.]

0.8

0.4

0.2

0

0

5

10

15

[a]

5

10

15

20

0.2 0.1

20

40

60

80

0

5

ρ = 3, topological distance ρ = 3, travelled distance ρ = 16, topological distance ρ = 16, travelled distance

0.4 0.3 0.2 0.1

0

20

40

60

80

10

15

20

25

ρ = 3, topological distance ρ = 3, travelled distance ρ = 16, topological distance ρ = 16, travelled distance

1.5

1

0.5

0

20

Topological distance

[d]

60

80

100

ρ = 3, topological distance ρ = 3, travelled distance ρ = 16, topological distance ρ = 16, travelled distance

15

10

5

0

100

40

Topological distance

First passage time [s]

crossing intensity [a.u.]

[d]

2

[c]

0.5

0

4

0

100

Topological distance

[c]

20

Topological distance

First passage time [s]

crossing intensity [a.u.]

0.3

0

15

2

ρ = 3, topological distance ρ = 3, travelled distance ρ = 16, topological distance ρ = 16, travelled distance

0.4

10

ρ = 3, topological distance ρ = 3, travelled distance ρ = 16, topological distance ρ = 16, travelled distance

[b]

0.5

0

5

6

0

Topological distance

[b]

0

Topological distance

First passage time [s]

crossing intensity [a.u.]

0.2

0

1

8

ρ = 3, topological distance ρ = 3, travelled distance ρ = 16, topological distance ρ = 16, travelled distance

0.4

0

2

[a]

0.8

0.6

3

0

20

Topological distance

ρ = 3, topological distance ρ = 3, travelled distance ρ = 16, topological distance ρ = 16, travelled distance

4

0

20

40

60

80

100

120

Topological distance

Figure 12: Illustration of diffusion coupling for samples A-D ([a]-[d]). The effects are more pronounced for the carbonate samples, in particular D.

Figure 11: Illustration of diffusion coupling for samples A-D ([a]-[d]). The effects are more pronounced for the carbonate samples, in particular D.

11

MMM

SPWLA 46th Annual Logging Symposium, June 26-29, 2005

distance) one has about 5 pores (for ρ = 3 µm/s). For the carbonates this maximum is closer to 10 (30) pores for ρ = 3(16) µm/s. If we look at the corresponding first passage time, we note that the two sandstones again are relatively similar, both showing a plateau in the first passage time over topological distance. In contrast, sample D shows almost a monotonic increase in the first passage time, which would suggest that one is near the tortuosity limit (Sen, 2004); it could also be influenced by the high coordination number of the network. The multiscale carbonate (C) shows an irregular behaviour. The simulations were run with surface relaxation only. To estimate the effect bulk relaxation may have on the crossing intensity over topological distance, we plot Fig. 11 and Fig. 12 next to each other. From the first passage time distribution one can estimate the crossing intensity at a given topological distance by multiplying with an exponential (exp(−t/T2b )).

[a]

[a’]

[b]

[b’]

[c]

[c’]

[d]

[d’]

Fluid typing In a previous section of this paper we discussed partial saturation effects on permeability predictions for sample B and assigning surface relaxivities only. In this section, we extend our analysis to consideration of 2D NMR correlations. Here we report T1 -T2 data. We simulate an inversion recovery sequence to acquire the data varying the diffusion time td and echo spacing te in the sequence (e.g. Fig. 1c of (H¨urlimann and Venkataramanan, 2002)). For both times we consider 150 logarithmically spaced values between 1 ms and 10 s. We use susceptibilities χ of -1 for water, -10 for oil, 50 for the grain phase, and 200 for the shaly region in units of 10−6 . Bulk relaxation times are 3 s for water and 1 s for oil. Further, we apply two different sets of surface relaxivities. The first (1) (1) set is ρwc = ρwg = 16µ m/s, where the subscript “c” is clay, “g” is grain, “w” is water and the superscript indi(1) cates T1 relaxivities. The second set is ρwc = 59µ m/s, (1) ρwg = 5.9, a 1:10 contrast in relaxivities. T2 surface re(2) (1) (2) (1) laxivities are set by ρwc = 4ρwc and ρwg = 2ρwg . The oil phase is assumed to be not in contact with the clay or grain phase (ρoc = ρog = 0). Further, it is assumed that the oil and water phases do not exchange any magnetisation. The diffusion constant of water D0 = 2500 µm/s2 is used for both the water and the oil phase. Fig 13 shows simulation data with no gradient applied and all susceptibility effects suppressed. Thus surface and bulk relaxation are independent of the applied fields. The T1 -T2 correlation spectrum of the 14% water saturation cases (Fig 13a-b) show one main peak caused by the oil phase, which experiences no surface relaxation, and some water in larger pores, while the remaining water signal of smaller pores is shifted towards smaller T1 , T2 , since surface relaxation is active in both for T1 and T2 relaxation. For the 47% water saturation case (Fig 13c-d), the same

Figure 13: T1 -T2 correlation on sample B. All susceptibility effects are suppressed. [a] SW =14%, set one, [b] SW =14%, set two, [c] SW =47%, set one, [d] SW =47%, set two.

can be observed with the difference, that the water peak, shifted by surface relaxation, is now much stronger. This we already noted in the previous section, i.e. surfaces active for relaxation become inaccessible due to the distribution of the oil phase. In the case of high water saturation the shape of the peak region is different between the two different surface relaxivity contrasts. Using exactly the same sample subset as in the previous case, and keeping the two different sets of surface relaxivities and water saturations, we now apply a static magnetic field B0 of 550 Gauss and susceptibilities [SI] of χ = 50×10−6 for the grain phase, χ = 100×10−6 for the shale region, and χ = −10×10−6 for the fluid phases (H¨urlimann, 1998). The internal fields were calculated according to a dipol approximation (Sen and Axelrod,

12

SPWLA 46th Annual Logging Symposium, June 26-29, 2005

partitions, e.g. grain partitions, will be considered in the future to allow greater control of varying surface relaxivity. Particular outcomes are listed in the following:

[a]

[a’]

[b]

[b’]

[c]

[c’]

[d]

1. Two common permeability correlations to NMR responses are tested. The correlation which includes the formation factor, and therefore information about the tortuosity of the sample, works significantly better. The range of prefactors found indicate an order of magnitude variation of permeability due to structure across 20 sandstones and two carbonates. The NMR response formation factor relationships have very high correlation coefficients. 2. Heterogeneity in surface relaxivity is shown to have an effect on permeability predictions of about a factor two for a particular sample. Further, a moderate diffusive coupling between the macro-porous sections of the image (resolved) and the clay regions (unresolved) leads to another 10-45% change in predicted permeability. Calibration of the background porosity might be facilitated by use of Heand Hg-porosimetry equipment. 3. Partial saturation effects can be masked by heterogeneity in surface relaxivity. For the sample studied in this paper they affect the log mean relaxation time in opposite ways. This effect could be severe, if the bulk relaxation time of the oil phase falls within the surface relaxation time distribution of the water phase.

[d’]

4. An analysis of pore-pore coupling using a pore partitioning shows that, for all rocks studied in detail, diffusion between pores occurs to a significant degree. For the sandstones considered here, typically five pores are traversed in 2 s with significant remaining signal, while for the carbonates about 30 pores are traversed in less than 1 s.

Figure 14: T1 -T2 correlation on sample B including susceptibility effects (B0 = 550 G). [a] SW =14%, set one, [b] SW =14%, set two, [c] SW =47%, set one, [d] SW =47%, set two.

1999; Song, 2003). The results are given in Fig. 14. Immediately one observes the large shift of the spectrum, particularly of the region with larger T1 , to smaller T2 values. This can be explained, since the short T1 part already relaxed, before the T2 part of the acquisition sequence starts. For longer T1 , now both fluids have a mechanism to relax close to the surface. This mechanism is stronger for the oil phase, since ∆χ is larger between the oil-grain or oil-clay surfaces.

ACKNOWLEDGEMENTS The authors acknowledge the Australian Government for their support through the ARC grant scheme and the Australian Partnership for Advanced Computing (APAC) for their support through the expertise program and APAC and the ANU Supercomputing Facility for very generous allocations of computer time. We also thank BHPBilliton and Woodside Energy who have provided financial support for the facility.

CONCLUSIONS In this paper we presented a comprehensive NMR response simulation tool for Xray-CT images including simulation of T1 , T2 relaxation and dephasing with a particle tracking method making use of a pore partitioning. Other

13

REFERENCES CITED Appel, M., Freeman, J. J., Gardner, J. S., Hirasaki, G. H., Zhang, Q. G., and Shafer, J. L., 2001, Interpretation of restricted diffusion in sandstones with inter-

MMM

SPWLA 46th Annual Logging Symposium, June 26-29, 2005

nal field gradients: Magnetic Resonance Imaging, 19, 535–537.

Cohen, M. H., and Mendelson, K. S., 1982, Nuclear magnetic relaxation and the internal geometry of sedimentary rocks: J. Appl. Phys., 53, no. 2, 1127–1135.

Arns, C. H., Knackstedt, M. A., Pinczewski, W. V., and Lindquist, W. B., 2001, Accurate estimation of transport properties from microtomographic images: Geophys. Res. Lett., 28, no. 17, 3361–3364.

Dunn, K.-J., 2001, Magnetic susceptibility contrast induced field gradients in porous media: Magnetic Resonance Imaging, 19, 439–442. Dunn, K.-J., 2002, Enhanced transverse relaxation in porous media due to internal field gradients: J. Magn. Res., 156, 171–180.

Arns, C. H., Knackstedt, M. A., Pinczewski, W. V., and Garboczi, E. G., 2002, Computation of linear elastic properties from microtomographic images: Methodology and agreement between theory and experiment: Geophysics, 67, no. 5, 1396–1405.

Ferreol, B., and Rothman, D. H., 1995, Latticeboltzmann simulations of flow through Fontainebleau sandstone: Transport in Porous Media, 20, 3–20.

Arns, C. H., Averdunk, H., Bauget, F., Sakellariou, A., Senden, T. J., Sheppard, A. P., Sok, R. M., Pinczewski, W. V., and Knackstedt, M. A., June 2004, Digital core laboratory: Analysis of reservoir core fragments from 3D images:, paper EEE.

Frisch, U., Hasslacher, B., and Pomeau, Y., 1986, Lattice-gas automata for the Navier-Stokes equation: Phys. Rev. Lett., 56, no. 14, 1505–1508. Hansen, P. C., 1992, Analysis of discrete ill-posed problems by means of the L-curve: SIAM review, 34, 561–580.

Arns, C. H., Knackstedt, M. A., Pinczewski, W. V., and Martys, N., Virtual permeametry on microtomographic images: J. Petr. Sc. Eng., 45, no. 1-2, 41–46.

Hilpert, M., and Miller, C. T., 2001, Pore-morphology based simulation of drainage in totally wetting porous media: Advances in Water Resources, 24, 243–255.

Arns, C. H., Sakellariou, A., Senden, T. J., Sheppard, A. P., Sok, R., Pinczewski, W. V., and Knackstedt, M. A., 2005, Digital core laboratory: Reservoir core analysis from 3D images: Petrophysics, 46, no. 3, to appear.

H¨urlimann, M. D., and Venkataramanan, L., 2002, Quantitative measurement of two-dimensional distribution functions of diffusion and relaxation in grossly inhomogeneous fields: J. Magn. Res., 157, 31–42.

Arns, C. H., 2004, A comparison of pore size distributions derived by NMR and Xray-CT techniques: Physica A, 339, no. 1-2, 159–165.

H¨urlimann, M. D., Flaum, M., Venkataramanan, L., Flaum, C., Freedman, R., and Hirasaki, G. J., 2003, Diffusion-relaxation distribution functions of sedimentary rocks in different saturation states: Magnetic Resonance Imaging, 21, 305–310.

Auzerais, F. M., Dunsmuir, J., Ferr`eol, B. B., Martys, N., Olson, J., Ramakrishnan, T. S., Rothman, D. H., and Schwartz, L. M., April 1996, Transport in sandstone: A study based on three dimensional microtromography: Geophysical Research Letters, 23, no. 7, 705–708.

H¨urlimann, M. D., 1998, Effective gradients in porous media due to susceptibility differences: J. Magn. Res., 131, 232–240.

Banavar, J. R., and Schwartz, L. M., 1987, Magnetic resonance as a probe of permeability in porous media: Phys. Rev. Lett., 58, 1411–1414.

Kenyon, W. E., Day, P., Straley, C., and Willemsen, J., 1986, Compact and consistent representation of rock NMR data from permeability estimation:, Proc. 61st Annual Technical Conference and Exhibition.

Bergman, D. J., Dunn, K.-J., Schwartz, L. M., and Mitra, P. P., 1995, Self-diffusion in a periodic porous medium: A comparison of different approaches: Phys. Rev. E, 51, 3393–3400.

Kenyon, W. E., Day, P., Straley, C., and Willemsen, J., 1988, A three part study of NMR longitudinal relaxation properties of water saturated sandstones: SPE formation evaluation, 3, no. 3, 626–636, SPE 15643.

Brownstein, K. R., and Tarr, C., 1979, Importance of classical diffusion in NMR studies of water in biological cells: Phys. Rev. A, 19, 2446–2453.

Kenyon, W. E., 1992, Nuclear magnetic resonance as a petrophysical measurement: Nucl. Geophys., 6, 153–171.

Chen, H., Chen, S., and Matthaeus, W. H., 1992, Recovery of the Navier-Stokes equations using a latticegas Boltzmann method: Phys. Rev. A, 45, no. 8, R5339–R5342.

Kim, I. C., and Torquato, S., 1990, Effective conductivity of suspensions of hard spheres by brownian motion simulation: J. Appl. Phys., 69, no. 4, 2280–2289.

14

SPWLA 46th Annual Logging Symposium, June 26-29, 2005

Knackstedt, M. A., Arns, C. H., Limaye, A., Sakellariou, A., Senden, T. J., Sheppard, A. P., Sok, R. M., Pinczewski, W. V., and Bunn, G. F., March 2004, Digital core laboratory: Properties of reservoir core derived from 3D images:, Presented at the 2004 SPE Asia Pacific Conference on Integrated Modelling for Asset Management, 1–14.

Sakellariou, A., Senden, T. J., Sawkins, T. J., Knackstedt, M. A., Turner, M. L., Jones, A. C., Saadatfar, M., Roberts, R. J., Limaye, A., Arns, C. H., Sheppard, A. P., and Sok, R. M., August 2004, An x-ray tomography facility for quantitative prediction of mechanical and transport properties in geological, biological and synthetic systems, Developments in X-Ray Tomography IV, 473–484.

Lawson, C. L., and Hansen, R. J., 1974, Solving least squares problems: Prentice-Hall.

Sakellariou, A., Sawkins, T. J., Senden, T. J., and Limaye, A., X-ray tomography for mesoscale physics applications: Physica A, 339, no. 1-2, 152–158.

Lonnes, S., Guzman-Garcia, A., and Holland, R., June 22-25 2003, NMR petrophysical predictions on cores: NMR petrophysical predictions on cores:, 44th Annual Logging Symposium.

Sen, P. N., and Axelrod, S., 1999, Inhomogeneity in local magnetic field due to susceptibility contrast: J. Appl. Phys., 89, no. 8, 4548–4554.

Looyestijn, W. J., June 1996, Determination of oil saturation from diffusion NMR logs: Determination of oil saturation from diffusion NMR logs:, SPWLA 37th Annual Logging Symposium, SS1–SS11.

Sen, P. N., Straley, C., Kenyon, W. E., and Whittingham, M. S., 1990, Surface-to-volume ratio, charge density, nuclear magnetic relaxation, and permeability in clay-bearing sandstones: Geophysics, 55, no. 1, 61–69.

Martys, N. S., and Chen, H., 1996, Simulation of multicomponent fluids in complex three-dimensional geometries by the lattice Boltzmann method: Phys. Rev. E, 53, no. 1, 743–750.

Sen, P. N., 2004, Time-dependent diffusion coefficient as a probe of geometry: Concepts in Magnetic Resonance Part A, 23A, no. 1, 1–21.

McCall, K. R., Johnson, D. L., and Guyer, R. A., 1991, Magnetization evolution in connected pore systems: Phys. Rev. B, 44, 7344–7355.

Sethian, J. A., 1999, Level set methods and fast marching methods: Cambridge University Press.

McCall, K. R., Guyer, R. A., and Johnson, D. L., 1993, Magnetization evolution in connected pore systems. II. pulsed-field-gradient NMR and pore-space geometry: Phys. Rev. B, 48, no. 9, 5997–6006.

Sheppard, A. P., Sok, R. M., and Averdunk, H., 2004, Techniques for image enhancement and segmentation of tomographic images of porous materials: Physica A, 339, no. 1-2, 145–151.

Mendelson, K. S., 1990, Percolation model of nuclear magnetic relaxation in porous media: Phys. Rev. B, 41, 562–567.

Song, Y.-Q., Ryu, S., and Sen, P. N., 2000, Determining multiple length scales in rocks: Nature, 406, no. 13, 178–181.

Neuman, C., 1974, Spin echo of spins diffusing in a bounded medium: J. Chem. Phys., 60, no. 11, 4508– 4511.

Song, Y.-Q., 2003, Using internal magnetic fields to obtain pore size distributions of porous media: Concepts in Magnetic Resonance Part A, 18A, no. 2, 97– 110.

Øren, P. E., Antonsen, F., Ruesl˚atten, H. G., and Bakke, S., September 2002, Numerical simulations of NMR responses for improved interpretations of NMR measurements in reservoir rocks:, Presented at the SPE 77th Annual Technical Conference and Exhibition.

Stark, P., and Parker, R., 1995, Bounded-variable least-squares: an algorithm and applications: Computational Statistics, 10, no. 2, 129–141. Sukstanskii, A. L., and Yablonskiy, D. A., 2002, Effects of restricted diffusion on MR signal formation: J. Magn. Res., 157, 92–105.

Qian, Y.-H., and Zhou, Y., 1998, Complete Galileaninvariant lattice BGK models for the Navier-Stokes equation: Europhys. Lett., 42, no. 4, 359–364.

Sun, B., and Dunn, K.-J., 2004, Methods and limitations of NMR data inversion for fluid typing: J. Magn. Res., 169, 118–128.

Robertson, B., 1966, Spin-echo decay of spins diffusing in a bounded region: Phys. Rev., 151, no. 1, 273– 277.

Toumelin, E., Torres-Verd´ın, C., Chen, S., and Fischer, D. M., 2003, Reconciling NMR measurements and numerical simulations: Assessment of temperature and diffusive coupling effects on two-phase carbonate samples: Petrophysics, 44, no. 2, 91–107.

Rothman, D. H., 1988, Cellular-automaton fluids: A model for flow in porous media: Geophysics, 53, no. 4, 509–518.

15

MMM

SPWLA 46th Annual Logging Symposium, June 26-29, 2005

Valckenborg, R. M. E., Huinink, H. P., v. d. Sande, J. J., and Kopinga, K., 2002, Random-walk simulations of NMR dephasing effects due to uniform magnetic-field gradients in a pore: Phys. Rev. E, 65, 21306.

M.A. Knackstedt: Mark Knackstedt was awarded a BSc in 1985 from Columbia University and a PhD in Chemical Engineering from Rice University in 1990. He is an Associate Professor at the Department of Applied Mathematics at the Australian National University and a visiting Fellow at the School of Petroleum Engineering at the University of NSW. His work has focussed on the characterisation and realistic modelling of disordered materials. His primary interests lie in modelling transport, elastic and multi-phase flow properties and development of 3D tomographic image analysis for complex materials.

Venkataramanan, L., Song, Y.-Q., and H¨urliman, M. D., 2002, Solving Fredholm integrals of the first kind with tensor product structure in 2 and 2.5 dimensions: IEEE Trans. on Signal Process., 50, no. 5, 1017–1026. Wayne, R. C., and Cotts, R. M., 1966, Nuclearmagnetic-resonance study of self-diffusion in a bounded medium: Phys. Rev., 151, no. 1, 264–272. Zhang, G. Q., and Hirasaki, G. J., 2003, CPMG relaxation by diffusion with constant magnetic field gradient in a restricted geometry: numerical simulation and application: J. Magn. Res., 163, 81–91. Zielinski, L. J., Song, Y.-Q., Ryu, S., and Sen, P. N., 2002, Characterization of coupled pore systems from the diffusion eigenspectrum: J. Chem. Phys., 117, no. 11, 5361–5365. ABOUT THE AUTHORS C. H. Arns: Christoph Arns was awarded a Diploma in Physics (1996) from the University of Technology Aachen and a PhD in Petroleum Engineering from the University of New South Wales in 2002. He is a Research Fellow at the Department of Applied Mathematics at the Australian National University. His research interests include the morphological analysis of porous complex media from 3D images and numerical calculation of transport and linear elastic properties with a current focus on NMR responses and dispersive flow. Member: AMPERE, ANZMAG, ARMA, DGG. A.P. Sheppard: Adrian Sheppard received his B.Sc. from the University of Adelaide in 1992 and his PhD in 1996 from the Australian National University and is currently a Research Fellow in the Department of Applied Mathematics at the Australian National University. His research interests are network modelling of multiphase fluid flow in porous material, topological analysis of complex structures, and tomographic image processing. R. M. Sok: Rob Sok studied chemistry and received his PhD (1994) at the University of Groningen in the Netherlands and is currently a Research Fellow in the Department of Applied Mathematics at the Australian National University. His main areas of interest are computational chemistry and structural analysis of porous materials.

16