- Email: [email protected]
Pore-morphology-based simulation of drainage in totally wetting porous media
Advances in Water Resources 24 (2001) 243±255 www.elsevier.com/locate/advwatres
Pore-morphology-based simulation of drainage in totally wetting porous media Markus Hilpert *, Cass T. Miller Center for Advanced Study of the Environment, Department of Environmental Sciences and Engineering, University of North Carolina, CB 7400, 104 Rosenau Hall, Chapel Hill, NC 27599-7400, USA Received 11 December 1999; received in revised form 29 May 2000; accepted 31 August 2000
Abstract We develop and analyze a novel, quasi-static, pore-scale approach for modeling drainage in a porous medium system. The approach uses: (1) a synthetic, non-overlapping packing of a set of spheres, (2) a discrete representation of the sphere packing, and (3) concepts from pore morphology and local pore-scale physics to simulate the drainage process. The grain-size distribution and porosity of two well-characterized porous media were used as input into the drainage simulator, and the simulated results showed good agreement with experimental observations. We further comment on the use of this simulator for determining the size of a representative elementary volume needed to characterize the drainage process. Ó 2001 Elsevier Science Ltd. All rights reserved. Keywords: Capillarity; Morphology; Network model; Percolation
1. Introduction Many engineering applications, such as groundwater remediation and enhanced oil recovery (EOR), involve multiphase ¯uid ¯ow in porous media. Though most of these applications use a continuum approach to describe the ¯ow processes, pore-scale modeling provides an important means to improve our understanding of the underlying physical processes and to determine macroscale constitutive relationships, such as the capillary pressure±saturation relation. Due to the natural pore space's complex morphology, ¯uid ¯ow has been modeled at the pore scale primarily using network models that rely upon idealized representations of the pore morphology. For example, a common approach is to represent the pore morphology by spheres that are connected by cylindrical throats (e.g. [14]). The most dicult task when using network models as a quantitative predictive tool is identifying and specifying coordination numbers and size distributions for pore bodies and pore throats [3]. Some noteworthy contributions, which derive the network structure from a pore-morphological analysis of detailed three-dimen*
Corresponding author. E-mail addresses: [email protected] [email protected] (C.T. Miller).
(M.
Hilpert),
sional pore geometries, can be found in [1,2,11,19,26]. But there are also pore-scale modeling approaches, which work with digital representations of the pore morphology. Lattice-gas and lattice-Boltzmann [22] have been increasingly used over the last two decades. While these approaches better represent the porous medium morphology and ¯ow process than idealized network models, they are computationally much more demanding [17]. Hazlett [9] recently suggested an approach for simulating quasi-static two-phase ¯ow based upon a size and connectivity analysis of the digital pore space. This heuristic approach is computationally very attractive and yields reasonable agreement with experimental data [5] but has not been widely used. The overall goal of this work is to develop an ecient and accurate method of simulating drainage of a wetting phase in a real porous medium system. The speci®c objectives are: (1) to develop an approach to link easily quanti®able macroscopic measures of a porous medium to an accurate mapping of the pore morphology; (2) to develop a methodology for using detailed information of pore morphology as a direct input into a drainage simulator without idealizing the morphology; (3) to compare simulations based upon the developed method with experimental data from well-characterized systems and with Hazlett's method; and (4) to utilize the developed method to investigate certain aspects of drainage, such
0309-1708/01/$ - see front matter Ó 2001 Elsevier Science Ltd. All rights reserved. PII: S 0 3 0 9 - 1 7 0 8 ( 0 0 ) 0 0 0 5 6 - 7
244
M. Hilpert, C.T. Miller / Advances in Water Resources 24 (2001) 243±255
sn snp
Notation NWP REV WP
non-wetting phase representative elementary volume wetting phase
sw X Vt Z U q c h
Variables B structuring element D diameter grain diameter Dg sphere diameter Ds d fractal dimension g gravitational acceleration capillary pressure head hc L domain length M mass of the percolating NWP cluster N number of spheres P pore space capillary pressure pc principal radius of curvature R1 principal radius of curvature R2 S sphere
Mathematical operations C that component of a set that is connected to the NWP reservoir D morphological dilation E morphological erosion O morphological opening Vol volume D standard deviation Minkowski addition Minkowski subtraction
as the representative elementary volume needed to characterize the process.
2. Background We ®rst introduce some concepts from mathematical morphology, which we will use to formulate our simulation approach. The textbooks of Matheron [15] and Serra [24] provide introductions to mathematical morphology that may be of interest to the reader; we present only a few fundamental notions. The morphological erosion E of a set X by a structuring element B is the locus of the centers ~ r of the B~r , which are included in X: n o r : B~r X : 1 EB X ~ E, B, and X are subsets of the underlying space (e.g., R2 or R3 ). The subscript ~ r denotes the translate of B by the vector ~ r. Another way of writing Eq. (1) is EB X X B;
NWP saturation NWP saturation at percolation point WP saturation set in R2 or R3 toroidal volume coordination number of the solid phase porosity density interfacial tension contact angle
2
where stands for the Minkowski subtraction and B for the re¯ected set of B with respect to the origin, B fÿ~ r :~ r 2 Bg. A distinction between B and B is only necessary for non-symmetric B. The actual choice of B depends on the speci®c application, but spheres are a common choice. Note that, unlike real spheres, their digital representations are not necessarily symmetric [8]. Fig. 1(a) shows a two-dimensional set X and a sym-
metric two-dimensional sphere S (circle). Fig. 1(b) shows the erosion of X by S S. The morphological dilation D of X by B is the locus of the centers of the B~r , which hit X: n o r : B~r \ X 6 ; ; 3 DB X ~ which can also be written as DB X X B;
4
where stands for the Minkowski addition. Fig. 1(c) shows the dilation of X, which was shown in Fig. 1(a), by S S. The mathematical opening O of X by B is the domain swept out by all the translates of B that are included in X: OB X [ B~r : B~r X :
5
The opening can be re-expressed by B: OB X X B
6
Fig. 1(d) shows the opening of X, which was shown in Fig. 1(a), by S S. The opening O is closely related to the morphological grain-size distribution if X represents the solid phase of a porous medium. This leads to the mathematical concept of a granulometry, which seeks to capture the result of a sieve analysis for granular media. A granulometry is essentially an opening with variable size structuring elements, kB, where k is a positive real number, i.e., the structuring elements are self-similar. The opening of X by kS is assumed to be that part of X that would be
M. Hilpert, C.T. Miller / Advances in Water Resources 24 (2001) 243±255
245
Fig. 1. Morphological operations in two dimensions: (a) set X and sphere S; (b) erosion E of X by S; (c) dilation D of X by S; (d) opening O of X by S.
withheld by a sieve with a mesh whose size equals the size of k > 0. Thus, one expects that a sieve with a small mesh size retains more of a given medium than a sieve with a larger mesh size, i.e., k0 6 k ) OkB Ok0 B :
7
This requirement only holds true if B is convex [15]. Then, the cumulative morphological grain-size distribution f k
VolOkB X ; VolX
8
is a monotonically decreasing function of k, which quanti®es the grain size. Vol stands for the volume of its
argument, being a subset of the underlying space. If X is a pore space, then f k represents the morphological pore-size distribution. For a digital representation, where the pore space is represented by voxels on a cubic lattice, the structuring elements B are naturally given by digital representations as well. There are many ways of de®ning digital spheres with integer diameter D as structuring elements. One possibility is 2
S D f~ r 2 N 3 : ~ r ÿ~ c 6 D2 =4g;
9
where the center point of the sphere is ci D=2 1=2 for i 1; 2; 3. Fig. 2 shows these symmetric digital spheres for D 1 to 4. The shape of S is not self-similar
246
M. Hilpert, C.T. Miller / Advances in Water Resources 24 (2001) 243±255
Fig. 2. Digital spheres S D for D 1±4.
Fig. 3. Modi®ed digital spheres S 0 D for D 1±4 and S 0 D S 0 D 1.
and not convex. Hence, D0 6 D generally does not imply OS D OS D0 . For example, it would not hold true for a cubic cavity, the walls of which contain cross-shaped holes with extension 3 voxels. Then, structuring spheres of diameter 3 would penetrate into the holes but not structuring spheres of diameter 2. Other examples can be easily constructed. In order to alleviate the violation of Eq. (7), Glantz [8], who investigated the percolation of solid particles in digital pore spaces, suggested the use of structuring elements S 0 with the property S 0 D S 0 D 1. The modi®ed structuring elements S 0 , being de®ned by S 0 1 S 1 and S 0 n S n [ S n ÿ 1 possess this property. The use of the S 0 weakens but does not prevent the violation of axiom (7) [8]. Fig. 3 shows these digital spheres for D 1 to 4. 3. Approach The physical model is a digital porous medium with the bottom connected to a non-wetting phase (NWP) reservoir and the top connected to a wetting phase (WP) reservoir. This is a model of a tempe cell, which is commonly used to determine capillary pressure±saturation relations. Our algorithm for modeling quasi-static primary drainage works as follows: 1. At the beginning, the porous medium is saturated with WP. NWP exists only in its reservoir. The capillary pressure pc is zero. 2. Then, pc is increased incrementally. The diameter Ds of a sphere is calculated from Laplace's equation for a spherical interface, pc 4c=Ds ;
10
where c is the interfacial tension. 3. The sphere serves as a probe. It makes new pore space accessible to the NWP, if it can be moved without intersecting with the solid phase from a location that is totally NWP-saturated to a neighboring location that might be partially WP-saturated. All probe
locations that are totally included in the NWP-®lled portion of the pore space are tested. This step is repeated until an equilibrium state is reached. 4. The NWP saturation sn is computed. 5. The algorithm proceeds with step (2). The methods imply a vanishing contact angle, h 0° and thus a totally wetting system. The approach also assumes that there is no trapped or irreducible WP. This requires either (1) that WP be connected through the edges of the pore space; (2) WP ®lms; or (3) a groove system on the solid surface. Our algorithm does not resolve WP ®lm ¯ow, and the groove system can only be modeled for an unrealistically ®ne resolution of the pore space. There is another slight inconsistency in our approach: a throat that is invaded from two sides by NWP (connected to its reservoir) does not get ®lled totally with NWP if the two menisci just touch each other as would happen in reality [20]. But we expect the error in sn introduced by this shortcoming to be only of minor importance because of the small volume of the throats. We use a pore-morphological framework in order to express sn formally during primary drainage. First, we i.e., the erosion of the determine ES D P P S D, pore space by a sphere of diameter D, which corresponds to the capillary pressure given by Eq. (10). We call CX the part of X that is connected to the NWP reservoir. We then identify CES D P S D with the NWP-®lled portion of P. The NWP saturation becomes h h i i Vol C P S D S D : 11 sn D VolP Similar to Eq. (7) for the morphological grain-size distribution, we want sn D 6 sn D0 if D0 6 D. For reasons already discussed, this requirement is often but not always ful®lled if digital spheres are used as structuring elements. But the use of the modi®ed spheres S 0 instead of the S should alleviate violations of it. For illustrative purposes, we applied the pore-morphology-based drainage algorithm to a two-dimensional pore space P. Fig. 4(a) shows ES P . Fig. 4(b) then
M. Hilpert, C.T. Miller / Advances in Water Resources 24 (2001) 243±255
247
Fig. 4. Two-dimensional illustration of drainage simulation in a two-dimensional pore space P. The solid phase is black. In red are shown (a) ES P , the erosion E of P by the circle S; (b) CES P , which is the component of ES P that is connected to the NWP reservoir; and (c) CES P S, which represents the NWP-®lled portion of P.
shows CES P , which is the component of ES P that is connected to the NWP reservoir. Finally Fig. 4(c) shows CES P S, which represents the NWP-®lled portion of P. The capillary pressure pc is inversely proportional to the diameter of the circle S that was used as the structuring element for the morphological operation. We implemented a Fortran90 code for pore-morphology-based drainage in three-dimensional, digital porous media. The centers of the solid phase voxels are located at subsets of Z 3 f 1=2; 1=2; 1=2g. Then, possible centers of the structuring elements are Z 3 f 1=2; 1=2; 1=2g and Z 3 for odd and even diameters, respectively. Connected components of digital sets were determined based upon the six nearest neighbors. Other centers for the structuring elements were not considered; this ensured that the structuring elements extended to the restricting solid phase. Fig. 5 illustrates the simulations for a three-dimensional digital porous medium with 353 voxels. Fig. 5(a) shows the porous medium. Figs. 5(b) and (c) show the NWP for sphere diameters D 6 and D 4 voxels, respectively. One can see from Eq. (11) or from Fig. 4 that one point on the drainage curve is obtained in a three-step approach, which involves an erosion, a connected component analysis, and a dilation. Both the erosion and the dilation involve boolean operations, and the values of these operations at any point in space can be calculated independently. The connected component analysis can either be performed recursively with boolean operations or in a ®xed number of steps with operations on integer variables [13,21]. Thus, our approach is, like lattice-Boltzmann, highly suitable for parallel computer platforms. But our algorithm requires only a small and ®xed number of operations on boolean or integer variables in order to determine a point on the primary drainage curve, whereas lattice-Boltzmann
methods need many time steps, involving real number operations, in order to achieve convergence. The computations of this work, for example, were performed on a SGI Origin 2000 using only one processor. The simulations on the largest domain considered (8003 voxels) required 2 GB of system memory.
4. Porous medium systems We simulated random sphere packings, which follow the grain-size statistics of experimental porous media, using the sphere-packing computer code developed by Yang et al. [28]. For the porous media, we used a uniform glass bead packing, labeled as GB1b, and a less uniform sand, labeled as C-109. See Table 1 for the porous medium properties. For both porous media, measured primary drainage curves were available, which were inferred from equilibrium tetrachloroethylene± water distributions in long vertical columns [10]. The surface tension was 36:23 0:21 dyn/cm. Although not being an input parameter for the simulator, we also report the densities of the ¯uids, which were 1:613 0:002 g=cm3 for tetrachloroethylene and 0:9978 0:002 g=cm3 for water. The contact angle was assumed to be zero, h 0°. The grain-size statistics were obtained using image analysis [4]. We assumed the distribution functions to be lognormal. Table 2 shows the properties of the simulated packings that represent the GB1b and C-109 porous media. Both packings contain approximately 10,000 spheres. The output of the sphere packing code (centers and diameters of the spheres, coordination number of the grains) was then used to generate digital representations of the porous media. We generated discretizations with 2003 , 4003 and 8003 voxels.
248
M. Hilpert, C.T. Miller / Advances in Water Resources 24 (2001) 243±255
Fig. 5. Three-dimensional illustration: (a) solid phase; (b) NWP for D 6 voxels; (c) NWP for D 4 voxels.
Table 1 Properties of the two experimental multiphase systems Arithmetic mean grain diameter Dg (mm) Geometric mean grain diameter Dg (mm) Arithmetic standard deviation of Dg (mm) Porosity U a
GB1b
C-109
0.1156a 0.115 0.0121 0:356 0:002
0.24 0.2182a 0.11 0:346 0:002
Based on the assumption of a lognormal grain-size distribution but not part of the experimental data set.
5. Simulation results 5.1. Discretization eects Fig. 6 shows the primary drainage curves for discretization containing 2003 , 4003 and 8003 voxels for both the simulated GB1b and C-109 systems described in Table 2. In all plots presented, the capillary pressures pc are expressed as heights of a corresponding water column, hc pc = qw g, where hc is the capillary pressure head, qw the density of water, and g 9:81 m=s2 is the gravitational acceleration. We used the modi®ed spheres S 0 as structuring elements. The dierences in the WP saturation sw for the same hc values due to dierent spatial discretization can be quite high in the ¯at parts of the primary drainage curves. The dierences may be on the order of 0.3 for the C-109 system, for example. Discretization eects occur because both the digital structuring elements and the digital pore space are not
self-similar as the resolution changes. The overall shape and entry pressure of the primary drainage curve, however, are not aected that much by the discretization. The data density increases with the resolution. The mean throat radius in a random packing of uniform spheres is given by hRt i 0:21Dg [18]. From that, one can estimate the entry pressure head hc 2c= 0:21Dg qw g. We tested the validity of this equation for non-uniform sphere packings by assuming that Dg stands for the arithmetic mean grain diameter hDg i. For the GB1b and C-109 systems, we estimated NWP entry at hc 31 cm and hc 15 cm, respectively, which is in acceptable agreement with both experimental and simulated data. In order to investigate the in¯uence of the structuring element on the primary drainage curve, we also used the symmetric structuring elements S D for the simulations in the GB1b medium. Fig. 7 shows hc versus the dierence in sw between the simulation with the modi®ed
M. Hilpert, C.T. Miller / Advances in Water Resources 24 (2001) 243±255
249
Table 2 Sphere-packing realizations for the GB1b and C-109 porous media
120 100 80 60 40 20
(a)
0 0
0.2
0.4
0.6
w
0.8
Saturation s
0.1149 0.1155 0.0116 2.35 0.356 9532 5.95
0.2397 0.2198 0.1024 5.78 0.345 10,666 5.99
60
experiment 2003 voxel 4003 voxel 3 800 voxel
2
experiment 3 200 voxel 4003 voxel 3 800 voxel
C-109
Capillary pressure head (cm H O)
140
2
Capillary pressure head (cm H O)
Arithmetic mean grain diameter Dg (mm) Geometric mean grain diameter Dg (mm) Arithmetic standard deviation of Dg (mm) Domain length L (mm) Porosity U Number of spheres N Coordination number Z
GB1b
1
50 40 30 20 10
(b) 0 0
0.2
0.4
0.6
Saturation sw
0.8
1
Capillary pressure head (cm H2O)
Fig. 6. In¯uence of spatial resolution on primary drainage curves: (a) GB1b system; (b) C-109 system.
140
from below. Interestingly, the coarsest discretization with 2003 voxels yields the smallest sw dierence, because at NWP breakthrough (hc 63 cm), the diameter of the structuring element is Ds 2 voxels, for which S and S 0 are identical. Overall, the structuring element does not have too much impact on the simulated primary drainage curve. The ®nal choice of the structuring element is somewhat arbitrary and not crucial for our application. Eventually, we followed Glantz [8] and used the modi®ed structuring elements S 0 for the remaining simulations.
2003 voxel 4003 voxel 3 800 voxel
120 100 80 60 40 20 0 0
0.025
w
0.05
s [S’] s
w
0.075
0.1
[S]
Fig. 7. hc versus sw S 0 ÿ sw S, the dierence in sw between the simulation with the modi®ed spheres S 0 D and the symmetric spheres S D for the GB1b system.
spheres S 0 D and the S D, sw S 0 ÿ sw S. This dierence is always positive, although S and S 0 are not convex. For our digital porous media, the condition S D S 0 D is sucient to ensure that the smaller S makes more space accessible to NWP than the larger S 0 . Large values of sw S 0 ÿ sw S only occur after NWP breakthrough; for small sn values, the porous medium does not experience the dierence between S and S 0 , because they look identical from the top and enter the porous medium
5.2. WP at very low saturations The open circles in Fig. 8 show again the simulation results for the discretization with 8003 voxels. The largest pc corresponds to a diameter Ds 1 voxel of the sphere S 0 . The simulation always predicts the minimum WP saturation to be zero, because it does not account for WP trapping, except if residual pore space exists between the grains, which is not the case for our digital porous media. Assuming that there is no WP trapping, the simulation overpredicts sw for low sw values, where all WP is pendular. This is because Eq. (10) for spherical interfaces does not hold true for pendular WP. The general form of the Young±Laplace equation
400
estimate for saddleshaped interface estimate for spherical interface experiment simulation corrected simulation
300
200
100
0 0
0.2
0.4
0.6
w
0.8
250
2
500
Capillary pressure head (cm H O)
M. Hilpert, C.T. Miller / Advances in Water Resources 24 (2001) 243±255
Capillary pressure head (cm H2O)
250
1
Saturation s (a)
200
estimate for saddleshaped interface estimate for spherical interface experiment simulation corrected simulation
150
100
50
0 0
0.2
0.4
0.6
Saturation sw
0.8
1
(b)
Fig. 8. Solid line ± pendular WP saturation based on the coordination number of the simulated packing and the Laplace equation for saddle-shaped interfaces. Dashed line ± pendular WP saturation based on the coordination number of the simulated packing and the Laplace equation for spherical interfaces. Closed circles ± experiment; open circles ± simulation; (a) GB1b system; (b) C-109 system.
1 1 p c ; R1 R 2 c
12
must be used, where R1 and R2 are the principal radii of curvature. For pendular WP, one of the principal radii is negative, the other positive. Pendular WP saturation predicted by Eq. (10) is larger than the one predicted by the correct Eq. (12). Signi®cant error is introduced by assuming implicitly a spherical interface when calculating pc . To show this, we estimated the primary drainage curve in the range of small WP saturations, where all WP is pendular, only from the sphere pack parameters, namely the mean grain size hDg i and the coordination number Z. Assuming that a WP ring between two solid spheres with diameter Dg has the shape of a torus, its volume is given by Vt 2p fj2 ÿ f 2 j cot j ÿ R21 jj fR21 ÿ Dg f 2 =2 where sin j Dg = Dg 2R1 , j Dg =2 R1 cos j, and f R1 sin j [16]. See Fig. 9 for the geometry. The negative principal radius of curvature, R2 , is given by q R2 R1 ÿ R21 Dg R1 : 13
Fig. 9. Geometry of pendular WP.
Assuming that Z of the simulated sphere packing corresponds to that of the experimental system and that all grains have the same diameter hDg i, the WP saturation is then given by sw NZVt = 2L3 U [16], which varies only with respect to R1 for a given porous media system. The capillary pressure can then be calculated from Eq. (12) and (13). Note that the torus approximation neglects the variation in R2 when calculating pc . The solid lines in Fig. 8 show the resulting primary drainage curves. If we falsely assume that pc 2c=R1 , we obtain the dashed lines, which overpredict sw signi®cantly. But then the match between simulated and estimated data is good for low sw values, because both methods rely on the use of Laplace's equation for a spherical interface. As a comparison between the solid lines and the ®lled circles in Fig. 8 shows, the experimental sw values for both porous media systems are larger in the range of high pc than the ones predicted above, where we assumed all WP to be pendular and the absence of WP trapping. These observations suggest the existence of WP other than pendular ± for example WP ®lms ± or WP trapping. We believe that a WP ®lm was likely to exist because the experiments started out with the WPsaturated, water-wet porous medium [23]. This fact and the surface roughness, which was likely greater for the sand grains of the C-109 medium than for the smoother glass beads of the GB1b system, exclude the occurrence of trapped WP [6]. The failure of the simulation in the pc calculation for pendular WP can be compensated by rescaling the pc axis. If one assumes again that all WP is pendular and that all grains have the same diameter hDg i, one can estimate the negative principal radius of curvature R2 , which is given by Eq. (13), and substitute this expression into Eq. (12). The positive radius of curvature is given
M. Hilpert, C.T. Miller / Advances in Water Resources 24 (2001) 243±255
The triangles in Fig. 8 show the modi®ed simulated primary drainage curves. As expected there is a good agreement between the corrected simulation and the estimate based upon both the sphere packing parameters and the general form of Laplace's equation for low sw values, as a comparison between the triangles and the solid lines shows. We do not suggest the rescaling of the pc axis as a means to obtain better simulation results, since there is no sharp and identi®able transition from pendular WP saturation to the one with spherical menisci. 5.3. Domain size eects
Capillary pressure head (cm H2O)
We used the discretizations with 4003 voxels and cut them into non-overlapping, rectangular, equal-sized subdomains by generating four subdomains with 200 200 400 voxels, 32 with 100 100 200 voxels, 256 with 50 50 100 voxels, and 2048 with 25 25 50 voxels. The subdomains were rectangular for reasons concerning percolation theory. We averaged the simulation results of the subdomains with equal size. Fig. 10 shows the primary drainage curve for the various subdomain sizes for both the GB1b and C-109 systems. For increasing domain size, the shoulder of the primary drainage curve at the entry point becomes sharper, consistent with the pore-network model simulations and experiments of Larson and Morrow [12]. The volume of NWP, which tries to invade the domain per area of the NWP reservoir, is independent of the domain size, resulting in smaller sn values for larger domains. Next, we sought to determine the fractal dimension of the percolating NWP cluster, which represents the NWP
experiment 25 x 25 x 50 voxels 50 x 50 x 100 voxels 100 x 100 x 200 voxels 200 x 200 x 400 voxels
120 100 80 60 40 20
(a) 0 0
0.2
0.4
0.6
Saturation sw
0.8
1
Capillary pressure head (cm H2O)
that ®rst extends to the boundary opposed to the NWP reservoir when increasing pc incrementally. An object possesses the spatial dimension d if its mass M scales as M Ld , where L is the domain length. We used the number of NWP-occupied voxels as a measure for M. The object is said to be fractal if d is a non-integer. Following Wilkinson and Willemsen [27], who investigated the fractal properties of the percolating cluster formed by an invasion percolation process, we used the rectangular domains of size L L 2L and determined the percolating cluster over the central region with size L L L. This geometry suppresses boundary eects resulting from the NWP reservoir and the sample outlet. The outlet causes a boundary eect because the NWP can advance opposite to the external capillary pressure. Contrary to the investigations of [27], we did not have periodic boundary conditions on the sides because the sphere packing was not periodic. Fig. 11 shows a double-logarithmic plots of M versus L including linear regression curves. The slope of these curves equals the fractal dimension of the percolating cluster. We obtained d 2:84 0:09 for the GB1a system, and d 2:89 0:12 for the C-109 system. We also determined the local fractal dimension, d d log M= d log L [25], and evaluated this expression only based on the results of the two largest subdomains. We obtained d 3:05 for both the GB1b and C-109 systems. The local fractal dimension is larger than the d obtained by use of all four data points, likely because of boundary eects: for example, the length of the smallest investigated subdomain with 25 25 50 voxels amounts to only very few sphere diameters. Thus, boundary eects dominate and reduce d [7]. Hence, the local fractal dimension provides a more faithful estimate for d, but error bounds are not provided by our analysis. As d 3:05 is not smaller than the spatial dimension 3, there is no evidence that the percolating NWP cluster is a fractal object. A local fractal dimension d 3:05 for both porous media systems suggests that the fractal dimension of the
by the radius of the structuring element, R1 Ds =2. The corrected capillary pressure then becomes # " 2c 1 0 p : 14 pc 1 Ds 1 ÿ 1 2hDg i=Ds
140
251
60
experiment 25 x 25 x 50 voxels 50 x 50 x 100 voxels 100 x 100 x 200 voxels 200 x 200 x 400 voxels
50 40 30 20 10
(b) 0 0
0.2
0.4
0.6
Saturation sw
0.8
1
Fig. 10. In¯uence of domain size on primary drainage curve keeping the spatial resolution constant: (a) GB1b system; (b) C-109 system.
252
M. Hilpert, C.T. Miller / Advances in Water Resources 24 (2001) 243±255 7
6.5
6
6
5.5
5.5
log10 (M)
log10 (M)
6.5
7
(a)
5 4.5
5 4.5
4
4
3.5
3.5
3 1
1.5
log
10
(L)
2
2.5
(b)
3 1
1.5
log10 (L)
2
2.5
Fig. 11. Mass M of the percolating cluster versus the domain length L. The fractal dimension d is given by the slope of the curve. (a) GB1b system. (b) C-109 system.
percolating NWP cluster predicted by the morphologybased approach is larger than d 2:52, the value for an invasion percolation process in three-dimensional networks [27]. This is not too surprising because the two approaches dier in some important respects. In a bond invasion percolation model, the NWP invasion starts from one side of the lattice. Then, the largest radius throat on the NWP-WP interface is searched and invaded by NWP. This whole process is repeated until no change in the ¯uid distribution occurs. For the case of the no WP trapping rule, the entire network is eventually NWP saturated. Many pore-network model formulations (e.g. [29]) are based upon the invasion percolation concept. In our algorithm, we increase the capillary pressure pc (decrease the size of the structuring element) in steps. We then check in the entire pore space as to whether the structuring element can access new pore space. Then, the entire accessible pore space is invaded by NWP simultaneously, not just the largest radius throat, as for the case of invasion percolation. Many quasi-static pore-network models are based on the assumption that the entire accessible pore space is invaded simultaneously (e.g. [14]). The morphological approach may yield a larger fractal dimension than invasion percolation theory, because the invasion of all accessible throats at once has more the character of a macroscopic piston ¯ow. The fact that the throat sizes are discretized increases this surge when increasing pc . Invasion percolation generates more fractal patterns, because ®ngers are more likely to develop when NWP invades the throats one by one. Invasion percolation mimics the dynamics of the displacement process. It assumes that throats with the least resistance (largest radius) tend to be invaded ®rst. This assumption is reasonable, but examples can be easily constructed where a simultaneous invasion into all accessible throats yields better results. Both invasion percolation and the morphological approach must
be seen as approximations of the actual dynamic processes. Whereas we observed that sn increases with L for large pc , Zhou and Stenby [29] observed the opposite using an invasion-percolation approach: their pore-network model included a trapping mechanism for WP, which we believe was not important for our porous medium systems. 5.4. Randomness of pore space Fig. 12 shows the standard deviation of the WP saturation, Dsw , versus hc for the four sizes of the simulated subdomain. Clearly, smaller subdomain sizes generate a larger scatter in sw . The C-109 system has a larger Dsw value than the GB1b system because of the larger variability of the grain sizes. Fig. 13 shows the maximum standard deviation of the WP saturation, max fpc : Dsw g, and the standard deviation of the porosity, DU, versus the height of the subdomain. For both the GB1b and the C-109 system, max fpc : Dsw g is larger than DU. This suggests that a REV with respect to the porosity is smaller than that with respect to the primary drainage curve. The randomness of the pore space has a negligible impact on the primary drainage curve if one considers domains with 200 200 400 voxels, which corresponds approximately to 2500 spheres. 5.5. Comparison to Hazlett's method Hazlett [9] suggested a model for quasi-static, twophase ¯ow in porous media, the drainage part of which is similar to our approach. In the language of poremorphology, Hazlett assumes that h hh i ii Vol C P S D S D : 15 sn D VolP
Capillary pressure head (cm H2O)
Capillary pressure head (cm H O)
M. Hilpert, C.T. Miller / Advances in Water Resources 24 (2001) 243±255 120
250
2
25 x 25 x 50 voxels 50 x 50 x 100 voxels 100 x 100 x 200 voxels 200 x 200 x 400 voxels
200
253
25 x 25 x 50 voxels 50 x 50 x 100 voxels 100 x 100 x 200 voxels 200 x 200 x 400 voxels
100 80
150
60
100
40
50
20
(a) 0 0
0.1
0.2
∆ sw
0.3
0.4
(b) 0 0
0.1
0.2
0.3
∆ sw
0.4
Fig. 12. hc versus standard deviation of the WP saturation, Dsw : (a) GB1b system; (b) C-109 system.
0.4
0.4 w
w
max { ∆ s } ∆Φ
0.3
max { ∆ s }
0.2
0.2
0.1
0.1
0 0 (a)
100
200
300
∆Φ
0.3
400
Domain height (voxel)
0 0 (b)
100
200
300
400
Domain height (voxel)
Fig. 13. Standard deviation of the porosity, DU, and maximum standard deviation of WP saturation, maxfpc : Dsw g, versus the height 2L of the rectangular subdomain: (a) GB1b system; (b) C-109 system.
acceptable agreement between experimental and simulated data. Their simulation overpredicted pc unlike ours. Uncertainties in interfacial tension c, contact angle h, and the CMT data seem to outweigh the negative bias introduced by the approximations used to represent drainage events.
NWP Reservoir
The deviation from our Eq. (11) seems minor but turns out to be crucial. We illustrate the dierence in a two-dimensional example. Fig. 14 shows two grains. The NWP reservoir is on the left side. We assume that pc and thus the circle diameter Ds is chosen such that the menisci formed by the opening O just touch the vertical symmetry axis. Hazlett's approach would assume that the entire opening O ES D P S D represents NWP. Clearly, this is not true since the meniscus on the left side still has to pass the throat. The meniscus does not know about the hypothetical NWP on the right side, to which it can connect. We implemented Hazlett's approach and simulated the drainage curve for the GB1b medium. The simulated pc in the horizontal portion of the primary drainage curve underpredicted the experimental pc on the order of 20%. This phenomenon is not a discretization eect. Coles et al. [5] used Hazlett's model to simulate capillary pressure-saturation curves in sandstone. The three-dimensional pore space was obtained by computed microtomography (CMT). Coles et al. obtained only an
Fig. 14. Drainage simulation in a two-dimensional pore space. The solid phase is black. Hazlett's method assumes that the opening O (shown in gray) represents NWP, because O is completely connected to the NWP reservoir. This is wrong, because the left meniscus cannot pass the throat.
254
M. Hilpert, C.T. Miller / Advances in Water Resources 24 (2001) 243±255
6. Discussion and summary We suggested a pore-morphology-based approach for modeling quasi-static drainage. We compared our simulations to experimental data and found very good agreement between simulated and experimental results in the horizontal part of the primary drainage curves. We attribute deviations at the entry point to size eects and mismatches in the minimum WP saturation to WP ®lms and grooves on the solid surface, neither of accounted for. The inappropriate use of Laplace's equation for spherical interfaces overpredicts sw in the range of high pc . Lattice-Boltzmann methods do not have this shortcoming, but they are computationally much more expensive. We obtained the digital representation of experimental porous media by using a sphere-packing algorithm [28], which only needs the grain-size distribution and the porosity as input parameters. The very good agreement between simulated and experimental data testi®es to the ability of the sphere-packing algorithm [28] used to model natural unconsolidated porous media accurately, for which a grain-size distribution can be readily obtained. Further, the numerical simulation of sphere packings is much less expensive than computed microtomography, and thus much larger porous medium domains can be investigated. The suggested simulator yields good predictions, if the menisci are spherical, which is the case for sphere packings, unconsolidated porous media, and sand stones in the range of high WP saturations. Other media, such as fractured rocks and clays, will likely yield worse predictions. The morphology-based approach is currently restricted to quasi-static drainage. But in conjunction with the ability to simulate natural porous media, it bears an important application: compared to lattice-Boltzmann simulations, the model is computationally very inexpensive, because it requires less memory and much less CPU time and thus allows the simulation of much larger domains. Thus, many random realizations of statistically identical pore spaces can be investigated and the random error of primary drainage curves estimated. This information can then be used to perform latticeBoltzmann simulations wisely, which also can be used to model imbibition and dynamic displacement processes. Acknowledgements This work was supported by the the National Science Foundation (NSF) grant EAR-9901660, the National Institute of Environmental Health Sciences (NIEHS) grant 5 P42 ES05948-02, and the Department of Energy (DOE) grant DE-FG07-96ER14703. The authors acknowledge the insightful comments and suggestions of
William G. Gray. The authors are also grateful to HansJ org Vogel, Martin Blunt, and one anonymous reviewer for their valuable comments. References [1] Bakke S, Oren P-E. 3-D pore-scale modelling of sandstones and ¯ow simulations in the pore networks. SPE J 1997;2:136±49. [2] Bryant SL, King PR, Mellor DW. Network model evaluation of permeability and spatial correlation in a real random sphere packing. Transport in Porous Media 1993;11(1):53±70. [3] Celia MA, Reeves PC, Ferrand LC. Recent advances in pore scale models for multiphase ¯ow in porous media. Rev Geophysics 1995;33(Suppl):1049±57. [4] Chen C-H, McBride JF. Measurement of mean particle size distribution of three samples using Optimas5 Image Analysis Software. Technical report. Chapel Hill: University of North Carolina; 1996. [5] Coles ME, Hazlett RD, Muegge EL, Jones KW, Andrews B, Dowd B, Siddons P, Peskin A, Spanne P, Soll W. Developments in synchrotron X-ray microtomography with applications to ¯ow in porous media. SPE Reservoir Eval Eng 1998;1(4): 288±96. [6] Dullien FAL, Zarcone C, MacDonald IF, Collins A, Bochard RDE. The eects of surface roughness on the capillary pressure curves and the heights of capillary rise in glass bead packs. J Colloid Interface Sci 1989;127(2):362±72. [7] Gabrielli A, Ca®ero R, Caldarelli G. Theory of boundary eects in invasion percolation. J Physics A 1998;31(37):7429±46. [8] Glantz R. Porennetzwerke von Erdsto-Filtern: Mathematischmorphologische Beschreibung kernspintomographischer Aufnahmen. PhD thesis. Karlsruhe: Universit at Karlsruhe; 1997. [9] Hazlett RD. Simulation of capillary-dominated displacements in microtomographic images of reservoir rocks. Transport in Porous Media 1995;20(1±2):21±35. [10] Hilpert M, McBride JF, Miller CT. Investigation of the residualfunicular nonwetting-phase-saturation relation. Adv Water Resour 2000 (in press). [11] Kwiecien MJ, Macdonald IF, Dullien FAL. Three-dimensional reconstruction of porous media from serial section data. J Microsc 1990;159(3):343±59. [12] Larson RG, Morrow NR. Eects of sample size on capillary pressures in porous media. Powder Technol 1981;30:123±38. [13] Lohmann G. Volumetric image analysis. Chichester: Wiley; 1998. [14] Lowry MI, Miller CT. Pore-scale modeling of non-wetting-phase residual in porous media. Water Resour Res 1995;31(3):455±73. [15] Matheron G. Random sets and integral geometry. New York: Wiley; 1975. [16] Mayer RP, Stowe RA. Mercury porosimetry: ®lling of toroidal void volume following breakthrough between packed spheres. J Phys Chem 1966;70(12):3867±73. [17] Miller CT, Christakos G, Imho PT, McBride JF, Pedit JA, Trangenstein JA. Multiphase ¯ow and transport modeling in heterogeneous porous media: challenges and approaches. Adv Water Resour 1998;21(2):77±120. [18] Ng KM, Davis HT, Scriven LE. Visualization of blob mechanics in ¯ow through porous media. Chem Eng Sci 1978;33:1009±17. [19] Oren P-E, Bakke S, Arntzen OJ. Extending predictive capabilities to network models. SPE Journal 1998;3:324±36. [20] Reeves PC, Celia MA. A functional relationship between capillary pressure, saturation, and interfacial area as revealed by a pore-scale network model. Water Resour Res 1996;32(8):2345± 58. [21] Ronse C, Devijver PA. Connected components in binary images: the detection problem. Research Studies Press, 1984.
M. Hilpert, C.T. Miller / Advances in Water Resources 24 (2001) 243±255 [22] Rothman DH, Zaleski S. Lattice-gas cellular automata: simple models of complex hydrodynamics. Cambridge: Cambridge University Press; 1997. [23] Schwille F. Migration of organic ¯uids immiscible with water in the unsaturated zone. In: Yaron B, Dagan G, Goldshmid J, editors. Pollutants in porous media. Berlin: Springer Verlag; 1984. [24] Serra J. Image analysis and mathematical morphology. London: Academic Press; 1982. [25] Stauer D. Introduction to percolation theory. Washington, DC: Taylor & Francis; 1992.
255
[26] Vogel HJ, Roth K. A new approach for determining eective soil hydraulic functions. Eur J Soil Sci 1998;49(4):547±56. [27] Wilkinson D, Willemsen JF. Invasion percolation: a new form of percolation theory. J Phys A 1983;16:3365±76. [28] Yang A, Miller CT, Turcoliver LD. Simulation of correlated and uncorrelated packing of random size spheres. Phys Rev E 1996;53(2):1516±24. [29] Zhou DG, Stenby EH. Interpretation of capillary pressure curves using invasion percolation theory. Transport in Porous Media 1993;11(1):17±31.
Pore-morphology-based simulation of drainage in totally wetting porous media Markus Hilpert *, Cass T. Miller Center for Advanced Study of the Environment, Department of Environmental Sciences and Engineering, University of North Carolina, CB 7400, 104 Rosenau Hall, Chapel Hill, NC 27599-7400, USA Received 11 December 1999; received in revised form 29 May 2000; accepted 31 August 2000
Abstract We develop and analyze a novel, quasi-static, pore-scale approach for modeling drainage in a porous medium system. The approach uses: (1) a synthetic, non-overlapping packing of a set of spheres, (2) a discrete representation of the sphere packing, and (3) concepts from pore morphology and local pore-scale physics to simulate the drainage process. The grain-size distribution and porosity of two well-characterized porous media were used as input into the drainage simulator, and the simulated results showed good agreement with experimental observations. We further comment on the use of this simulator for determining the size of a representative elementary volume needed to characterize the drainage process. Ó 2001 Elsevier Science Ltd. All rights reserved. Keywords: Capillarity; Morphology; Network model; Percolation
1. Introduction Many engineering applications, such as groundwater remediation and enhanced oil recovery (EOR), involve multiphase ¯uid ¯ow in porous media. Though most of these applications use a continuum approach to describe the ¯ow processes, pore-scale modeling provides an important means to improve our understanding of the underlying physical processes and to determine macroscale constitutive relationships, such as the capillary pressure±saturation relation. Due to the natural pore space's complex morphology, ¯uid ¯ow has been modeled at the pore scale primarily using network models that rely upon idealized representations of the pore morphology. For example, a common approach is to represent the pore morphology by spheres that are connected by cylindrical throats (e.g. [14]). The most dicult task when using network models as a quantitative predictive tool is identifying and specifying coordination numbers and size distributions for pore bodies and pore throats [3]. Some noteworthy contributions, which derive the network structure from a pore-morphological analysis of detailed three-dimen*
Corresponding author. E-mail addresses: [email protected] [email protected] (C.T. Miller).
(M.
Hilpert),
sional pore geometries, can be found in [1,2,11,19,26]. But there are also pore-scale modeling approaches, which work with digital representations of the pore morphology. Lattice-gas and lattice-Boltzmann [22] have been increasingly used over the last two decades. While these approaches better represent the porous medium morphology and ¯ow process than idealized network models, they are computationally much more demanding [17]. Hazlett [9] recently suggested an approach for simulating quasi-static two-phase ¯ow based upon a size and connectivity analysis of the digital pore space. This heuristic approach is computationally very attractive and yields reasonable agreement with experimental data [5] but has not been widely used. The overall goal of this work is to develop an ecient and accurate method of simulating drainage of a wetting phase in a real porous medium system. The speci®c objectives are: (1) to develop an approach to link easily quanti®able macroscopic measures of a porous medium to an accurate mapping of the pore morphology; (2) to develop a methodology for using detailed information of pore morphology as a direct input into a drainage simulator without idealizing the morphology; (3) to compare simulations based upon the developed method with experimental data from well-characterized systems and with Hazlett's method; and (4) to utilize the developed method to investigate certain aspects of drainage, such
0309-1708/01/$ - see front matter Ó 2001 Elsevier Science Ltd. All rights reserved. PII: S 0 3 0 9 - 1 7 0 8 ( 0 0 ) 0 0 0 5 6 - 7
244
M. Hilpert, C.T. Miller / Advances in Water Resources 24 (2001) 243±255
sn snp
Notation NWP REV WP
non-wetting phase representative elementary volume wetting phase
sw X Vt Z U q c h
Variables B structuring element D diameter grain diameter Dg sphere diameter Ds d fractal dimension g gravitational acceleration capillary pressure head hc L domain length M mass of the percolating NWP cluster N number of spheres P pore space capillary pressure pc principal radius of curvature R1 principal radius of curvature R2 S sphere
Mathematical operations C that component of a set that is connected to the NWP reservoir D morphological dilation E morphological erosion O morphological opening Vol volume D standard deviation Minkowski addition Minkowski subtraction
as the representative elementary volume needed to characterize the process.
2. Background We ®rst introduce some concepts from mathematical morphology, which we will use to formulate our simulation approach. The textbooks of Matheron [15] and Serra [24] provide introductions to mathematical morphology that may be of interest to the reader; we present only a few fundamental notions. The morphological erosion E of a set X by a structuring element B is the locus of the centers ~ r of the B~r , which are included in X: n o r : B~r X : 1 EB X ~ E, B, and X are subsets of the underlying space (e.g., R2 or R3 ). The subscript ~ r denotes the translate of B by the vector ~ r. Another way of writing Eq. (1) is EB X X B;
NWP saturation NWP saturation at percolation point WP saturation set in R2 or R3 toroidal volume coordination number of the solid phase porosity density interfacial tension contact angle
2
where stands for the Minkowski subtraction and B for the re¯ected set of B with respect to the origin, B fÿ~ r :~ r 2 Bg. A distinction between B and B is only necessary for non-symmetric B. The actual choice of B depends on the speci®c application, but spheres are a common choice. Note that, unlike real spheres, their digital representations are not necessarily symmetric [8]. Fig. 1(a) shows a two-dimensional set X and a sym-
metric two-dimensional sphere S (circle). Fig. 1(b) shows the erosion of X by S S. The morphological dilation D of X by B is the locus of the centers of the B~r , which hit X: n o r : B~r \ X 6 ; ; 3 DB X ~ which can also be written as DB X X B;
4
where stands for the Minkowski addition. Fig. 1(c) shows the dilation of X, which was shown in Fig. 1(a), by S S. The mathematical opening O of X by B is the domain swept out by all the translates of B that are included in X: OB X [ B~r : B~r X :
5
The opening can be re-expressed by B: OB X X B
6
Fig. 1(d) shows the opening of X, which was shown in Fig. 1(a), by S S. The opening O is closely related to the morphological grain-size distribution if X represents the solid phase of a porous medium. This leads to the mathematical concept of a granulometry, which seeks to capture the result of a sieve analysis for granular media. A granulometry is essentially an opening with variable size structuring elements, kB, where k is a positive real number, i.e., the structuring elements are self-similar. The opening of X by kS is assumed to be that part of X that would be
M. Hilpert, C.T. Miller / Advances in Water Resources 24 (2001) 243±255
245
Fig. 1. Morphological operations in two dimensions: (a) set X and sphere S; (b) erosion E of X by S; (c) dilation D of X by S; (d) opening O of X by S.
withheld by a sieve with a mesh whose size equals the size of k > 0. Thus, one expects that a sieve with a small mesh size retains more of a given medium than a sieve with a larger mesh size, i.e., k0 6 k ) OkB Ok0 B :
7
This requirement only holds true if B is convex [15]. Then, the cumulative morphological grain-size distribution f k
VolOkB X ; VolX
8
is a monotonically decreasing function of k, which quanti®es the grain size. Vol stands for the volume of its
argument, being a subset of the underlying space. If X is a pore space, then f k represents the morphological pore-size distribution. For a digital representation, where the pore space is represented by voxels on a cubic lattice, the structuring elements B are naturally given by digital representations as well. There are many ways of de®ning digital spheres with integer diameter D as structuring elements. One possibility is 2
S D f~ r 2 N 3 : ~ r ÿ~ c 6 D2 =4g;
9
where the center point of the sphere is ci D=2 1=2 for i 1; 2; 3. Fig. 2 shows these symmetric digital spheres for D 1 to 4. The shape of S is not self-similar
246
M. Hilpert, C.T. Miller / Advances in Water Resources 24 (2001) 243±255
Fig. 2. Digital spheres S D for D 1±4.
Fig. 3. Modi®ed digital spheres S 0 D for D 1±4 and S 0 D S 0 D 1.
and not convex. Hence, D0 6 D generally does not imply OS D OS D0 . For example, it would not hold true for a cubic cavity, the walls of which contain cross-shaped holes with extension 3 voxels. Then, structuring spheres of diameter 3 would penetrate into the holes but not structuring spheres of diameter 2. Other examples can be easily constructed. In order to alleviate the violation of Eq. (7), Glantz [8], who investigated the percolation of solid particles in digital pore spaces, suggested the use of structuring elements S 0 with the property S 0 D S 0 D 1. The modi®ed structuring elements S 0 , being de®ned by S 0 1 S 1 and S 0 n S n [ S n ÿ 1 possess this property. The use of the S 0 weakens but does not prevent the violation of axiom (7) [8]. Fig. 3 shows these digital spheres for D 1 to 4. 3. Approach The physical model is a digital porous medium with the bottom connected to a non-wetting phase (NWP) reservoir and the top connected to a wetting phase (WP) reservoir. This is a model of a tempe cell, which is commonly used to determine capillary pressure±saturation relations. Our algorithm for modeling quasi-static primary drainage works as follows: 1. At the beginning, the porous medium is saturated with WP. NWP exists only in its reservoir. The capillary pressure pc is zero. 2. Then, pc is increased incrementally. The diameter Ds of a sphere is calculated from Laplace's equation for a spherical interface, pc 4c=Ds ;
10
where c is the interfacial tension. 3. The sphere serves as a probe. It makes new pore space accessible to the NWP, if it can be moved without intersecting with the solid phase from a location that is totally NWP-saturated to a neighboring location that might be partially WP-saturated. All probe
locations that are totally included in the NWP-®lled portion of the pore space are tested. This step is repeated until an equilibrium state is reached. 4. The NWP saturation sn is computed. 5. The algorithm proceeds with step (2). The methods imply a vanishing contact angle, h 0° and thus a totally wetting system. The approach also assumes that there is no trapped or irreducible WP. This requires either (1) that WP be connected through the edges of the pore space; (2) WP ®lms; or (3) a groove system on the solid surface. Our algorithm does not resolve WP ®lm ¯ow, and the groove system can only be modeled for an unrealistically ®ne resolution of the pore space. There is another slight inconsistency in our approach: a throat that is invaded from two sides by NWP (connected to its reservoir) does not get ®lled totally with NWP if the two menisci just touch each other as would happen in reality [20]. But we expect the error in sn introduced by this shortcoming to be only of minor importance because of the small volume of the throats. We use a pore-morphological framework in order to express sn formally during primary drainage. First, we i.e., the erosion of the determine ES D P P S D, pore space by a sphere of diameter D, which corresponds to the capillary pressure given by Eq. (10). We call CX the part of X that is connected to the NWP reservoir. We then identify CES D P S D with the NWP-®lled portion of P. The NWP saturation becomes h h i i Vol C P S D S D : 11 sn D VolP Similar to Eq. (7) for the morphological grain-size distribution, we want sn D 6 sn D0 if D0 6 D. For reasons already discussed, this requirement is often but not always ful®lled if digital spheres are used as structuring elements. But the use of the modi®ed spheres S 0 instead of the S should alleviate violations of it. For illustrative purposes, we applied the pore-morphology-based drainage algorithm to a two-dimensional pore space P. Fig. 4(a) shows ES P . Fig. 4(b) then
M. Hilpert, C.T. Miller / Advances in Water Resources 24 (2001) 243±255
247
Fig. 4. Two-dimensional illustration of drainage simulation in a two-dimensional pore space P. The solid phase is black. In red are shown (a) ES P , the erosion E of P by the circle S; (b) CES P , which is the component of ES P that is connected to the NWP reservoir; and (c) CES P S, which represents the NWP-®lled portion of P.
shows CES P , which is the component of ES P that is connected to the NWP reservoir. Finally Fig. 4(c) shows CES P S, which represents the NWP-®lled portion of P. The capillary pressure pc is inversely proportional to the diameter of the circle S that was used as the structuring element for the morphological operation. We implemented a Fortran90 code for pore-morphology-based drainage in three-dimensional, digital porous media. The centers of the solid phase voxels are located at subsets of Z 3 f 1=2; 1=2; 1=2g. Then, possible centers of the structuring elements are Z 3 f 1=2; 1=2; 1=2g and Z 3 for odd and even diameters, respectively. Connected components of digital sets were determined based upon the six nearest neighbors. Other centers for the structuring elements were not considered; this ensured that the structuring elements extended to the restricting solid phase. Fig. 5 illustrates the simulations for a three-dimensional digital porous medium with 353 voxels. Fig. 5(a) shows the porous medium. Figs. 5(b) and (c) show the NWP for sphere diameters D 6 and D 4 voxels, respectively. One can see from Eq. (11) or from Fig. 4 that one point on the drainage curve is obtained in a three-step approach, which involves an erosion, a connected component analysis, and a dilation. Both the erosion and the dilation involve boolean operations, and the values of these operations at any point in space can be calculated independently. The connected component analysis can either be performed recursively with boolean operations or in a ®xed number of steps with operations on integer variables [13,21]. Thus, our approach is, like lattice-Boltzmann, highly suitable for parallel computer platforms. But our algorithm requires only a small and ®xed number of operations on boolean or integer variables in order to determine a point on the primary drainage curve, whereas lattice-Boltzmann
methods need many time steps, involving real number operations, in order to achieve convergence. The computations of this work, for example, were performed on a SGI Origin 2000 using only one processor. The simulations on the largest domain considered (8003 voxels) required 2 GB of system memory.
4. Porous medium systems We simulated random sphere packings, which follow the grain-size statistics of experimental porous media, using the sphere-packing computer code developed by Yang et al. [28]. For the porous media, we used a uniform glass bead packing, labeled as GB1b, and a less uniform sand, labeled as C-109. See Table 1 for the porous medium properties. For both porous media, measured primary drainage curves were available, which were inferred from equilibrium tetrachloroethylene± water distributions in long vertical columns [10]. The surface tension was 36:23 0:21 dyn/cm. Although not being an input parameter for the simulator, we also report the densities of the ¯uids, which were 1:613 0:002 g=cm3 for tetrachloroethylene and 0:9978 0:002 g=cm3 for water. The contact angle was assumed to be zero, h 0°. The grain-size statistics were obtained using image analysis [4]. We assumed the distribution functions to be lognormal. Table 2 shows the properties of the simulated packings that represent the GB1b and C-109 porous media. Both packings contain approximately 10,000 spheres. The output of the sphere packing code (centers and diameters of the spheres, coordination number of the grains) was then used to generate digital representations of the porous media. We generated discretizations with 2003 , 4003 and 8003 voxels.
248
M. Hilpert, C.T. Miller / Advances in Water Resources 24 (2001) 243±255
Fig. 5. Three-dimensional illustration: (a) solid phase; (b) NWP for D 6 voxels; (c) NWP for D 4 voxels.
Table 1 Properties of the two experimental multiphase systems Arithmetic mean grain diameter Dg (mm) Geometric mean grain diameter Dg (mm) Arithmetic standard deviation of Dg (mm) Porosity U a
GB1b
C-109
0.1156a 0.115 0.0121 0:356 0:002
0.24 0.2182a 0.11 0:346 0:002
Based on the assumption of a lognormal grain-size distribution but not part of the experimental data set.
5. Simulation results 5.1. Discretization eects Fig. 6 shows the primary drainage curves for discretization containing 2003 , 4003 and 8003 voxels for both the simulated GB1b and C-109 systems described in Table 2. In all plots presented, the capillary pressures pc are expressed as heights of a corresponding water column, hc pc = qw g, where hc is the capillary pressure head, qw the density of water, and g 9:81 m=s2 is the gravitational acceleration. We used the modi®ed spheres S 0 as structuring elements. The dierences in the WP saturation sw for the same hc values due to dierent spatial discretization can be quite high in the ¯at parts of the primary drainage curves. The dierences may be on the order of 0.3 for the C-109 system, for example. Discretization eects occur because both the digital structuring elements and the digital pore space are not
self-similar as the resolution changes. The overall shape and entry pressure of the primary drainage curve, however, are not aected that much by the discretization. The data density increases with the resolution. The mean throat radius in a random packing of uniform spheres is given by hRt i 0:21Dg [18]. From that, one can estimate the entry pressure head hc 2c= 0:21Dg qw g. We tested the validity of this equation for non-uniform sphere packings by assuming that Dg stands for the arithmetic mean grain diameter hDg i. For the GB1b and C-109 systems, we estimated NWP entry at hc 31 cm and hc 15 cm, respectively, which is in acceptable agreement with both experimental and simulated data. In order to investigate the in¯uence of the structuring element on the primary drainage curve, we also used the symmetric structuring elements S D for the simulations in the GB1b medium. Fig. 7 shows hc versus the dierence in sw between the simulation with the modi®ed
M. Hilpert, C.T. Miller / Advances in Water Resources 24 (2001) 243±255
249
Table 2 Sphere-packing realizations for the GB1b and C-109 porous media
120 100 80 60 40 20
(a)
0 0
0.2
0.4
0.6
w
0.8
Saturation s
0.1149 0.1155 0.0116 2.35 0.356 9532 5.95
0.2397 0.2198 0.1024 5.78 0.345 10,666 5.99
60
experiment 2003 voxel 4003 voxel 3 800 voxel
2
experiment 3 200 voxel 4003 voxel 3 800 voxel
C-109
Capillary pressure head (cm H O)
140
2
Capillary pressure head (cm H O)
Arithmetic mean grain diameter Dg (mm) Geometric mean grain diameter Dg (mm) Arithmetic standard deviation of Dg (mm) Domain length L (mm) Porosity U Number of spheres N Coordination number Z
GB1b
1
50 40 30 20 10
(b) 0 0
0.2
0.4
0.6
Saturation sw
0.8
1
Capillary pressure head (cm H2O)
Fig. 6. In¯uence of spatial resolution on primary drainage curves: (a) GB1b system; (b) C-109 system.
140
from below. Interestingly, the coarsest discretization with 2003 voxels yields the smallest sw dierence, because at NWP breakthrough (hc 63 cm), the diameter of the structuring element is Ds 2 voxels, for which S and S 0 are identical. Overall, the structuring element does not have too much impact on the simulated primary drainage curve. The ®nal choice of the structuring element is somewhat arbitrary and not crucial for our application. Eventually, we followed Glantz [8] and used the modi®ed structuring elements S 0 for the remaining simulations.
2003 voxel 4003 voxel 3 800 voxel
120 100 80 60 40 20 0 0
0.025
w
0.05
s [S’] s
w
0.075
0.1
[S]
Fig. 7. hc versus sw S 0 ÿ sw S, the dierence in sw between the simulation with the modi®ed spheres S 0 D and the symmetric spheres S D for the GB1b system.
spheres S 0 D and the S D, sw S 0 ÿ sw S. This dierence is always positive, although S and S 0 are not convex. For our digital porous media, the condition S D S 0 D is sucient to ensure that the smaller S makes more space accessible to NWP than the larger S 0 . Large values of sw S 0 ÿ sw S only occur after NWP breakthrough; for small sn values, the porous medium does not experience the dierence between S and S 0 , because they look identical from the top and enter the porous medium
5.2. WP at very low saturations The open circles in Fig. 8 show again the simulation results for the discretization with 8003 voxels. The largest pc corresponds to a diameter Ds 1 voxel of the sphere S 0 . The simulation always predicts the minimum WP saturation to be zero, because it does not account for WP trapping, except if residual pore space exists between the grains, which is not the case for our digital porous media. Assuming that there is no WP trapping, the simulation overpredicts sw for low sw values, where all WP is pendular. This is because Eq. (10) for spherical interfaces does not hold true for pendular WP. The general form of the Young±Laplace equation
400
estimate for saddleshaped interface estimate for spherical interface experiment simulation corrected simulation
300
200
100
0 0
0.2
0.4
0.6
w
0.8
250
2
500
Capillary pressure head (cm H O)
M. Hilpert, C.T. Miller / Advances in Water Resources 24 (2001) 243±255
Capillary pressure head (cm H2O)
250
1
Saturation s (a)
200
estimate for saddleshaped interface estimate for spherical interface experiment simulation corrected simulation
150
100
50
0 0
0.2
0.4
0.6
Saturation sw
0.8
1
(b)
Fig. 8. Solid line ± pendular WP saturation based on the coordination number of the simulated packing and the Laplace equation for saddle-shaped interfaces. Dashed line ± pendular WP saturation based on the coordination number of the simulated packing and the Laplace equation for spherical interfaces. Closed circles ± experiment; open circles ± simulation; (a) GB1b system; (b) C-109 system.
1 1 p c ; R1 R 2 c
12
must be used, where R1 and R2 are the principal radii of curvature. For pendular WP, one of the principal radii is negative, the other positive. Pendular WP saturation predicted by Eq. (10) is larger than the one predicted by the correct Eq. (12). Signi®cant error is introduced by assuming implicitly a spherical interface when calculating pc . To show this, we estimated the primary drainage curve in the range of small WP saturations, where all WP is pendular, only from the sphere pack parameters, namely the mean grain size hDg i and the coordination number Z. Assuming that a WP ring between two solid spheres with diameter Dg has the shape of a torus, its volume is given by Vt 2p fj2 ÿ f 2 j cot j ÿ R21 jj fR21 ÿ Dg f 2 =2 where sin j Dg = Dg 2R1 , j Dg =2 R1 cos j, and f R1 sin j [16]. See Fig. 9 for the geometry. The negative principal radius of curvature, R2 , is given by q R2 R1 ÿ R21 Dg R1 : 13
Fig. 9. Geometry of pendular WP.
Assuming that Z of the simulated sphere packing corresponds to that of the experimental system and that all grains have the same diameter hDg i, the WP saturation is then given by sw NZVt = 2L3 U [16], which varies only with respect to R1 for a given porous media system. The capillary pressure can then be calculated from Eq. (12) and (13). Note that the torus approximation neglects the variation in R2 when calculating pc . The solid lines in Fig. 8 show the resulting primary drainage curves. If we falsely assume that pc 2c=R1 , we obtain the dashed lines, which overpredict sw signi®cantly. But then the match between simulated and estimated data is good for low sw values, because both methods rely on the use of Laplace's equation for a spherical interface. As a comparison between the solid lines and the ®lled circles in Fig. 8 shows, the experimental sw values for both porous media systems are larger in the range of high pc than the ones predicted above, where we assumed all WP to be pendular and the absence of WP trapping. These observations suggest the existence of WP other than pendular ± for example WP ®lms ± or WP trapping. We believe that a WP ®lm was likely to exist because the experiments started out with the WPsaturated, water-wet porous medium [23]. This fact and the surface roughness, which was likely greater for the sand grains of the C-109 medium than for the smoother glass beads of the GB1b system, exclude the occurrence of trapped WP [6]. The failure of the simulation in the pc calculation for pendular WP can be compensated by rescaling the pc axis. If one assumes again that all WP is pendular and that all grains have the same diameter hDg i, one can estimate the negative principal radius of curvature R2 , which is given by Eq. (13), and substitute this expression into Eq. (12). The positive radius of curvature is given
M. Hilpert, C.T. Miller / Advances in Water Resources 24 (2001) 243±255
The triangles in Fig. 8 show the modi®ed simulated primary drainage curves. As expected there is a good agreement between the corrected simulation and the estimate based upon both the sphere packing parameters and the general form of Laplace's equation for low sw values, as a comparison between the triangles and the solid lines shows. We do not suggest the rescaling of the pc axis as a means to obtain better simulation results, since there is no sharp and identi®able transition from pendular WP saturation to the one with spherical menisci. 5.3. Domain size eects
Capillary pressure head (cm H2O)
We used the discretizations with 4003 voxels and cut them into non-overlapping, rectangular, equal-sized subdomains by generating four subdomains with 200 200 400 voxels, 32 with 100 100 200 voxels, 256 with 50 50 100 voxels, and 2048 with 25 25 50 voxels. The subdomains were rectangular for reasons concerning percolation theory. We averaged the simulation results of the subdomains with equal size. Fig. 10 shows the primary drainage curve for the various subdomain sizes for both the GB1b and C-109 systems. For increasing domain size, the shoulder of the primary drainage curve at the entry point becomes sharper, consistent with the pore-network model simulations and experiments of Larson and Morrow [12]. The volume of NWP, which tries to invade the domain per area of the NWP reservoir, is independent of the domain size, resulting in smaller sn values for larger domains. Next, we sought to determine the fractal dimension of the percolating NWP cluster, which represents the NWP
experiment 25 x 25 x 50 voxels 50 x 50 x 100 voxels 100 x 100 x 200 voxels 200 x 200 x 400 voxels
120 100 80 60 40 20
(a) 0 0
0.2
0.4
0.6
Saturation sw
0.8
1
Capillary pressure head (cm H2O)
that ®rst extends to the boundary opposed to the NWP reservoir when increasing pc incrementally. An object possesses the spatial dimension d if its mass M scales as M Ld , where L is the domain length. We used the number of NWP-occupied voxels as a measure for M. The object is said to be fractal if d is a non-integer. Following Wilkinson and Willemsen [27], who investigated the fractal properties of the percolating cluster formed by an invasion percolation process, we used the rectangular domains of size L L 2L and determined the percolating cluster over the central region with size L L L. This geometry suppresses boundary eects resulting from the NWP reservoir and the sample outlet. The outlet causes a boundary eect because the NWP can advance opposite to the external capillary pressure. Contrary to the investigations of [27], we did not have periodic boundary conditions on the sides because the sphere packing was not periodic. Fig. 11 shows a double-logarithmic plots of M versus L including linear regression curves. The slope of these curves equals the fractal dimension of the percolating cluster. We obtained d 2:84 0:09 for the GB1a system, and d 2:89 0:12 for the C-109 system. We also determined the local fractal dimension, d d log M= d log L [25], and evaluated this expression only based on the results of the two largest subdomains. We obtained d 3:05 for both the GB1b and C-109 systems. The local fractal dimension is larger than the d obtained by use of all four data points, likely because of boundary eects: for example, the length of the smallest investigated subdomain with 25 25 50 voxels amounts to only very few sphere diameters. Thus, boundary eects dominate and reduce d [7]. Hence, the local fractal dimension provides a more faithful estimate for d, but error bounds are not provided by our analysis. As d 3:05 is not smaller than the spatial dimension 3, there is no evidence that the percolating NWP cluster is a fractal object. A local fractal dimension d 3:05 for both porous media systems suggests that the fractal dimension of the
by the radius of the structuring element, R1 Ds =2. The corrected capillary pressure then becomes # " 2c 1 0 p : 14 pc 1 Ds 1 ÿ 1 2hDg i=Ds
140
251
60
experiment 25 x 25 x 50 voxels 50 x 50 x 100 voxels 100 x 100 x 200 voxels 200 x 200 x 400 voxels
50 40 30 20 10
(b) 0 0
0.2
0.4
0.6
Saturation sw
0.8
1
Fig. 10. In¯uence of domain size on primary drainage curve keeping the spatial resolution constant: (a) GB1b system; (b) C-109 system.
252
M. Hilpert, C.T. Miller / Advances in Water Resources 24 (2001) 243±255 7
6.5
6
6
5.5
5.5
log10 (M)
log10 (M)
6.5
7
(a)
5 4.5
5 4.5
4
4
3.5
3.5
3 1
1.5
log
10
(L)
2
2.5
(b)
3 1
1.5
log10 (L)
2
2.5
Fig. 11. Mass M of the percolating cluster versus the domain length L. The fractal dimension d is given by the slope of the curve. (a) GB1b system. (b) C-109 system.
percolating NWP cluster predicted by the morphologybased approach is larger than d 2:52, the value for an invasion percolation process in three-dimensional networks [27]. This is not too surprising because the two approaches dier in some important respects. In a bond invasion percolation model, the NWP invasion starts from one side of the lattice. Then, the largest radius throat on the NWP-WP interface is searched and invaded by NWP. This whole process is repeated until no change in the ¯uid distribution occurs. For the case of the no WP trapping rule, the entire network is eventually NWP saturated. Many pore-network model formulations (e.g. [29]) are based upon the invasion percolation concept. In our algorithm, we increase the capillary pressure pc (decrease the size of the structuring element) in steps. We then check in the entire pore space as to whether the structuring element can access new pore space. Then, the entire accessible pore space is invaded by NWP simultaneously, not just the largest radius throat, as for the case of invasion percolation. Many quasi-static pore-network models are based on the assumption that the entire accessible pore space is invaded simultaneously (e.g. [14]). The morphological approach may yield a larger fractal dimension than invasion percolation theory, because the invasion of all accessible throats at once has more the character of a macroscopic piston ¯ow. The fact that the throat sizes are discretized increases this surge when increasing pc . Invasion percolation generates more fractal patterns, because ®ngers are more likely to develop when NWP invades the throats one by one. Invasion percolation mimics the dynamics of the displacement process. It assumes that throats with the least resistance (largest radius) tend to be invaded ®rst. This assumption is reasonable, but examples can be easily constructed where a simultaneous invasion into all accessible throats yields better results. Both invasion percolation and the morphological approach must
be seen as approximations of the actual dynamic processes. Whereas we observed that sn increases with L for large pc , Zhou and Stenby [29] observed the opposite using an invasion-percolation approach: their pore-network model included a trapping mechanism for WP, which we believe was not important for our porous medium systems. 5.4. Randomness of pore space Fig. 12 shows the standard deviation of the WP saturation, Dsw , versus hc for the four sizes of the simulated subdomain. Clearly, smaller subdomain sizes generate a larger scatter in sw . The C-109 system has a larger Dsw value than the GB1b system because of the larger variability of the grain sizes. Fig. 13 shows the maximum standard deviation of the WP saturation, max fpc : Dsw g, and the standard deviation of the porosity, DU, versus the height of the subdomain. For both the GB1b and the C-109 system, max fpc : Dsw g is larger than DU. This suggests that a REV with respect to the porosity is smaller than that with respect to the primary drainage curve. The randomness of the pore space has a negligible impact on the primary drainage curve if one considers domains with 200 200 400 voxels, which corresponds approximately to 2500 spheres. 5.5. Comparison to Hazlett's method Hazlett [9] suggested a model for quasi-static, twophase ¯ow in porous media, the drainage part of which is similar to our approach. In the language of poremorphology, Hazlett assumes that h hh i ii Vol C P S D S D : 15 sn D VolP
Capillary pressure head (cm H2O)
Capillary pressure head (cm H O)
M. Hilpert, C.T. Miller / Advances in Water Resources 24 (2001) 243±255 120
250
2
25 x 25 x 50 voxels 50 x 50 x 100 voxels 100 x 100 x 200 voxels 200 x 200 x 400 voxels
200
253
25 x 25 x 50 voxels 50 x 50 x 100 voxels 100 x 100 x 200 voxels 200 x 200 x 400 voxels
100 80
150
60
100
40
50
20
(a) 0 0
0.1
0.2
∆ sw
0.3
0.4
(b) 0 0
0.1
0.2
0.3
∆ sw
0.4
Fig. 12. hc versus standard deviation of the WP saturation, Dsw : (a) GB1b system; (b) C-109 system.
0.4
0.4 w
w
max { ∆ s } ∆Φ
0.3
max { ∆ s }
0.2
0.2
0.1
0.1
0 0 (a)
100
200
300
∆Φ
0.3
400
Domain height (voxel)
0 0 (b)
100
200
300
400
Domain height (voxel)
Fig. 13. Standard deviation of the porosity, DU, and maximum standard deviation of WP saturation, maxfpc : Dsw g, versus the height 2L of the rectangular subdomain: (a) GB1b system; (b) C-109 system.
acceptable agreement between experimental and simulated data. Their simulation overpredicted pc unlike ours. Uncertainties in interfacial tension c, contact angle h, and the CMT data seem to outweigh the negative bias introduced by the approximations used to represent drainage events.
NWP Reservoir
The deviation from our Eq. (11) seems minor but turns out to be crucial. We illustrate the dierence in a two-dimensional example. Fig. 14 shows two grains. The NWP reservoir is on the left side. We assume that pc and thus the circle diameter Ds is chosen such that the menisci formed by the opening O just touch the vertical symmetry axis. Hazlett's approach would assume that the entire opening O ES D P S D represents NWP. Clearly, this is not true since the meniscus on the left side still has to pass the throat. The meniscus does not know about the hypothetical NWP on the right side, to which it can connect. We implemented Hazlett's approach and simulated the drainage curve for the GB1b medium. The simulated pc in the horizontal portion of the primary drainage curve underpredicted the experimental pc on the order of 20%. This phenomenon is not a discretization eect. Coles et al. [5] used Hazlett's model to simulate capillary pressure-saturation curves in sandstone. The three-dimensional pore space was obtained by computed microtomography (CMT). Coles et al. obtained only an
Fig. 14. Drainage simulation in a two-dimensional pore space. The solid phase is black. Hazlett's method assumes that the opening O (shown in gray) represents NWP, because O is completely connected to the NWP reservoir. This is wrong, because the left meniscus cannot pass the throat.
254
M. Hilpert, C.T. Miller / Advances in Water Resources 24 (2001) 243±255
6. Discussion and summary We suggested a pore-morphology-based approach for modeling quasi-static drainage. We compared our simulations to experimental data and found very good agreement between simulated and experimental results in the horizontal part of the primary drainage curves. We attribute deviations at the entry point to size eects and mismatches in the minimum WP saturation to WP ®lms and grooves on the solid surface, neither of accounted for. The inappropriate use of Laplace's equation for spherical interfaces overpredicts sw in the range of high pc . Lattice-Boltzmann methods do not have this shortcoming, but they are computationally much more expensive. We obtained the digital representation of experimental porous media by using a sphere-packing algorithm [28], which only needs the grain-size distribution and the porosity as input parameters. The very good agreement between simulated and experimental data testi®es to the ability of the sphere-packing algorithm [28] used to model natural unconsolidated porous media accurately, for which a grain-size distribution can be readily obtained. Further, the numerical simulation of sphere packings is much less expensive than computed microtomography, and thus much larger porous medium domains can be investigated. The suggested simulator yields good predictions, if the menisci are spherical, which is the case for sphere packings, unconsolidated porous media, and sand stones in the range of high WP saturations. Other media, such as fractured rocks and clays, will likely yield worse predictions. The morphology-based approach is currently restricted to quasi-static drainage. But in conjunction with the ability to simulate natural porous media, it bears an important application: compared to lattice-Boltzmann simulations, the model is computationally very inexpensive, because it requires less memory and much less CPU time and thus allows the simulation of much larger domains. Thus, many random realizations of statistically identical pore spaces can be investigated and the random error of primary drainage curves estimated. This information can then be used to perform latticeBoltzmann simulations wisely, which also can be used to model imbibition and dynamic displacement processes. Acknowledgements This work was supported by the the National Science Foundation (NSF) grant EAR-9901660, the National Institute of Environmental Health Sciences (NIEHS) grant 5 P42 ES05948-02, and the Department of Energy (DOE) grant DE-FG07-96ER14703. The authors acknowledge the insightful comments and suggestions of
William G. Gray. The authors are also grateful to HansJ org Vogel, Martin Blunt, and one anonymous reviewer for their valuable comments. References [1] Bakke S, Oren P-E. 3-D pore-scale modelling of sandstones and ¯ow simulations in the pore networks. SPE J 1997;2:136±49. [2] Bryant SL, King PR, Mellor DW. Network model evaluation of permeability and spatial correlation in a real random sphere packing. Transport in Porous Media 1993;11(1):53±70. [3] Celia MA, Reeves PC, Ferrand LC. Recent advances in pore scale models for multiphase ¯ow in porous media. Rev Geophysics 1995;33(Suppl):1049±57. [4] Chen C-H, McBride JF. Measurement of mean particle size distribution of three samples using Optimas5 Image Analysis Software. Technical report. Chapel Hill: University of North Carolina; 1996. [5] Coles ME, Hazlett RD, Muegge EL, Jones KW, Andrews B, Dowd B, Siddons P, Peskin A, Spanne P, Soll W. Developments in synchrotron X-ray microtomography with applications to ¯ow in porous media. SPE Reservoir Eval Eng 1998;1(4): 288±96. [6] Dullien FAL, Zarcone C, MacDonald IF, Collins A, Bochard RDE. The eects of surface roughness on the capillary pressure curves and the heights of capillary rise in glass bead packs. J Colloid Interface Sci 1989;127(2):362±72. [7] Gabrielli A, Ca®ero R, Caldarelli G. Theory of boundary eects in invasion percolation. J Physics A 1998;31(37):7429±46. [8] Glantz R. Porennetzwerke von Erdsto-Filtern: Mathematischmorphologische Beschreibung kernspintomographischer Aufnahmen. PhD thesis. Karlsruhe: Universit at Karlsruhe; 1997. [9] Hazlett RD. Simulation of capillary-dominated displacements in microtomographic images of reservoir rocks. Transport in Porous Media 1995;20(1±2):21±35. [10] Hilpert M, McBride JF, Miller CT. Investigation of the residualfunicular nonwetting-phase-saturation relation. Adv Water Resour 2000 (in press). [11] Kwiecien MJ, Macdonald IF, Dullien FAL. Three-dimensional reconstruction of porous media from serial section data. J Microsc 1990;159(3):343±59. [12] Larson RG, Morrow NR. Eects of sample size on capillary pressures in porous media. Powder Technol 1981;30:123±38. [13] Lohmann G. Volumetric image analysis. Chichester: Wiley; 1998. [14] Lowry MI, Miller CT. Pore-scale modeling of non-wetting-phase residual in porous media. Water Resour Res 1995;31(3):455±73. [15] Matheron G. Random sets and integral geometry. New York: Wiley; 1975. [16] Mayer RP, Stowe RA. Mercury porosimetry: ®lling of toroidal void volume following breakthrough between packed spheres. J Phys Chem 1966;70(12):3867±73. [17] Miller CT, Christakos G, Imho PT, McBride JF, Pedit JA, Trangenstein JA. Multiphase ¯ow and transport modeling in heterogeneous porous media: challenges and approaches. Adv Water Resour 1998;21(2):77±120. [18] Ng KM, Davis HT, Scriven LE. Visualization of blob mechanics in ¯ow through porous media. Chem Eng Sci 1978;33:1009±17. [19] Oren P-E, Bakke S, Arntzen OJ. Extending predictive capabilities to network models. SPE Journal 1998;3:324±36. [20] Reeves PC, Celia MA. A functional relationship between capillary pressure, saturation, and interfacial area as revealed by a pore-scale network model. Water Resour Res 1996;32(8):2345± 58. [21] Ronse C, Devijver PA. Connected components in binary images: the detection problem. Research Studies Press, 1984.
M. Hilpert, C.T. Miller / Advances in Water Resources 24 (2001) 243±255 [22] Rothman DH, Zaleski S. Lattice-gas cellular automata: simple models of complex hydrodynamics. Cambridge: Cambridge University Press; 1997. [23] Schwille F. Migration of organic ¯uids immiscible with water in the unsaturated zone. In: Yaron B, Dagan G, Goldshmid J, editors. Pollutants in porous media. Berlin: Springer Verlag; 1984. [24] Serra J. Image analysis and mathematical morphology. London: Academic Press; 1982. [25] Stauer D. Introduction to percolation theory. Washington, DC: Taylor & Francis; 1992.
255
[26] Vogel HJ, Roth K. A new approach for determining eective soil hydraulic functions. Eur J Soil Sci 1998;49(4):547±56. [27] Wilkinson D, Willemsen JF. Invasion percolation: a new form of percolation theory. J Phys A 1983;16:3365±76. [28] Yang A, Miller CT, Turcoliver LD. Simulation of correlated and uncorrelated packing of random size spheres. Phys Rev E 1996;53(2):1516±24. [29] Zhou DG, Stenby EH. Interpretation of capillary pressure curves using invasion percolation theory. Transport in Porous Media 1993;11(1):17±31.