VZJ sign up for etocs
HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS
QUICK SEARCH:   [advanced]


     

Published online 13 May 2005
Published in Vadose Zone J 4:380-388 (2005)
DOI: 10.2136/vzj2004.0114
© 2005 Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA
This Article
Right arrow Abstract Freely available
Right arrow Figures Only
Right arrow Full Text (PDF) Free
Right arrow Alert me when this article is cited
Right arrow Alert me if a correction is posted
Services
Right arrow Similar articles in this journal
Right arrow Similar articles in ISI Web of Science
Right arrow Alert me to new issues of the journal
Right arrow Download to citation manager
Right arrow reprints & permissions
Citing Articles
Right arrow Citing Articles via ISI Web of Science (8)
Right arrow Citing Articles via Google Scholar
Google Scholar
Right arrow Articles by Vogel, H.-J.
Right arrow Articles by Roth, K.
Right arrow Search for Related Content
PubMed
Right arrow Articles by Vogel, H.-J.
Right arrow Articles by Roth, K.
GeoRef
Right arrow GeoRef Citation
Agricola
Right arrow Articles by Vogel, H.-J.
Right arrow Articles by Roth, K.
Related Collections
Right arrow Soil Physics
Right arrow Tomography
Right arrow Pore-Scale Modeling

ORIGINAL RESEARCH

Comparison of a Lattice-Boltzmann Model, a Full-Morphology Model, and a Pore Network Model for Determining Capillary Pressure–Saturation Relationships

H.-J. Vogela,*, J. Tölkeb, V. P. Schulza, M. Krafczykb and K. Rotha

a Institute of Environmental Physics, Univ. of Heidelberg, INF 229, 69120 Heidelberg, Germany
b Institute for Computer Applications in Civil Engineering, TU Braunschweig, Pockelsstr. 3, 38106 Braunschweig, Germany
* Corresponding author (hjvogel{at}iup.uni-heidelberg.de)

Received 5 August 2004.



   ABSTRACT
TOP
 ABSTRACT
INTRODUCTION
 MATERIALS AND METHODS
 RESULTS AND DISCUSSION
 CONCLUSIONS
 REFERENCES
 
Effective hydraulic properties of porous media such as the capillary pressure–saturation relation and the hydraulic conductivity function are a direct manifestation of the underlying pore geometry. The porous structure of a macroscopically homogeneous porous medium (sintered glass) was measured in detail using X-ray microtomography. We investigated the possibility of deriving the water characteristic on the basis of structure analysis. We compared three approaches differing in the amount of required input data, the effort of data processing, and their predictive potential. With increasing complexity, these were (i) a simple pore network model based on a few structural parameters and a simplified process model, (ii) a direct simulation of the pressure–saturation relation based on the full morphology of the porous structure together with a simplified process model, and (iii) Lattice-Boltzmann (LB) simulations based on the full morphology of the porous structure and modeling the complete multiphase fluid dynamics. We found that the network model, the most simple approach, may be sufficient for estimating the main drainage curve, thus suggesting that the complex pore structure can be reduced to a small set of geometric properties. If dynamic effects are considered, the LB approach provides a reliable representation, but at the cost of substantial computational efforts. Since very thin water films exist in the dry range, the continuity of the water phase cannot be described correctly with the LB approach, which presently uses uniform grids.

Abbreviations: D3Q19, three-dimensional 19 velocity LB model • FE, free energy [model] • FM, full-morphology [model] • LB, Lattice-Boltzmann [model] • NW, network [model] • NWP, nonwetting phase • REV, representative elementary volume • RK, Rothman–Keller • SC, Shan–Chen


   INTRODUCTION
TOP
 ABSTRACT
 INTRODUCTION
MATERIALS AND METHODS
 RESULTS AND DISCUSSION
 CONCLUSIONS
 REFERENCES
 
MODELING FLUID FLOW through porous media at a scale considerably larger than the size of individual pores is typically based on constitutive relations or macroscopic material properties. Such properties represent an averaged "effective" description of the complicated geometry of the pore space. The constitutive relations are the capillary pressure–saturation relation sw,n(pc), also termed the soil water characteristic, and the conductivity for the wetting and the nonwetting fluid as a function of fluid saturation Kw,n(sw,n). The subscripts "w" and "n" denote the wetting and the nonwetting phases (NWP), respectively.

The use of constitutive relations implies that the sample is large enough to allow a meaningful average over the microscopic heterogeneities at the pore scale in the sense of a representative elementary volume (REV). Once the constitutive relations are known, the dynamics of fluids at the larger scale can be described using appropriate partial differential equations. Depending on the specific situation, only the wetting phase may need to be considered, (e.g., the Richards equation), or the NWP is considered additionally.

Clearly, the constitutive relations are a direct manifestation of the complicated geometry of the underlying pore space. However, sw,n(pc) and Kw,n(sw,n) are typically not derived from detailed representations of the pore space, but are determined experimentally (Klute, 1986).

The experimental approach is subject to various difficulties. The direct measurement of the conductivity function is expensive, and the results are notoriously inaccurate. As a remedy, the hydraulic conductivity function is typically derived from sw,n(pc) following the approach of Mualem (1976) based on a statistical distribution of pores. The same concept can be used to estimate the conductivity of the NWP (Fischer et al., 1997). Given a parameterization of the constitutive relations, they can be obtained from multistep outflow experiments together with inverse parameter optimization (Van Dam et al., 1994). This approach may be affected by dynamic effects, as pointed out by Hassainzadeh et al. (2002). Another drawback is that ad hoc parameterizations of sw,n(pc) and Kw,n(sw,n) are then required. An additional difficulty arises when multiphase phenomena lead to hysteretic material properties, which are difficult to quantify using classical experiments. Hysteretic properties are typically approximated by means of simplified models (e.g., Kool and Parker, 1987).

An alternative way to determine the constitutive relations is to start from the microscopic structure of the pore space and to follow a process-oriented approach. Nondestructive methods are currently available which allow for direct analysis of the porous structure. This includes X-ray tomography using medical devices with a typical resolution of a few millimeters (Hopmans et al., 1994) but also fine-focus sources with resolutions down to a few microns (Tidwell et al., 2000). A recent review on applications of X-ray tomography in hydrology was given by Wildenschild et al. (2002).

Since the constitutive relations including all multiphase phenomena are a direct consequence of the microscopic porous structure and the physicochemical properties of surfaces, it is tempting to derive those properties from a direct analysis of the structure of the porous medium.

We focus on the shape of the pore space and presume simple and uniform surface properties. Following such a path to determine hydraulic material properties is attractive, since structure is directly observable. Concerning the material properties, we focus on the soil water characteristic and investigate three different approaches to derive this relationship from the underlying porous structure. The three approaches differ in the amount of required input data, the effort of data processing, and in the predictive potential. The most sophisticated and demanding approach is to measure the porous structure in all detail with high spatial resolution and simulate the fluid dynamics of the multiphase system. This is done using a LB simulation (Succi, 2001). A simpler approach is to neglect the dynamic process and to consider only static situations with constant capillary pressures. This is done using the known relation between pc and the curvature radius of the water–gas interface and by direct measurement of the three-dimensional pore geometry, taking into account the connectivity of the wetting and nonwetting fluid. This approach will be referred to as the full-morphology (FM) approach. Finally, the most simple but attractive approach is to reduce the description of the pore morphology to a small set of meaningful parameters from which the constitutive relationships can be derived. This approach is implemented by means of a network (NW) model, which is adjusted to a statistical description of the pore-size distribution and the pore connectivity and permits the prediction of hydraulic properties.

The LB and FM approaches are based on a complete description of the pore structure and hence can only be applied to porous media for which the REV is small enough to allow an appropriate resolution of the smallest single pores (e.g., by high-resolution X-ray tomography). These two requirements, REV and sufficient resolution, are satisfied for macroscopically homogeneous and isotropic materials with a structure having short characteristic lengths as homogeneous sands or sandstones. We chose rigid sintered glass for this study so that the different approaches can be compared without possible restrictions due to nonstationarity of the structure in space and time. The different approaches were evaluated with respect to their specific potential, their limitations, and the related costs.


   MATERIALS AND METHODS
TOP
 ABSTRACT
 INTRODUCTION
 MATERIALS AND METHODS
RESULTS AND DISCUSSION
 CONCLUSIONS
 REFERENCES
 
All simulations were based on data from a sintered borosilicate glass, which contains approximately 81% SiO2 and 13% BO and has a density of about 2.2 x 103 kg m–3.

The choice of the material was motivated by its rigidity, which leads to a well-defined pore space, and its isotropy, which greatly simplifies the analysis. Furthermore, the material was sufficiently uniform so that a geometric REV had a characteristic length of just a few grain diameters. Consequently, the complete structure of the pore space could be measured by high-resolution X-ray tomography with high contrast due to the large difference between the absorption coefficients of air and SiO2.

Microtomography and Three-Dimensional Representation of the Pore Space
A cylindrical sample of 7-mm diameter and 8.5-mm height was scanned with X-ray tomography at a resolution of 15 µm. We used an industrial micro-CT scanner (Frauenhofer-Institute for Nondestructive Testing, IZFP, Saarbrücken) with a polychromatic cone beam having a focal spot size of some 4 µm. The three-dimensional structure of X-ray attenuation was obtained as a gray scale image after tomographic reconstruction (Fig. 1 , left). A 3 x 3 x 3 median filter was then applied to the gray-scale image of absorption coefficients to eliminate artifacts from the reconstruction. Next, the image was segmented into pore space and solid by a bilevel thresholding method (Vogel and Kretzschmar, 1996). The two thresholds define a narrow region of the gray scale where it is not obvious if a pixel belongs to the pore space or not. This decision depends on the nature of the adjacent pixels. Hence, the bilevel thresholding takes into account the connectivity of the resulting image and thus further reduces artifacts by removing small isolated pores. Such pores are implausible in view of the way the sintered glass was generated. The thresholds were chosen such that the porosity of the binarized image was the same as the true value, 0.43, which was calculated from the sample volume and the known density of SiO2.



View larger version (75K):
[in this window]
[in a new window]
  		Fig. 1. (left) Original gray-scale image of reconstructed X-ray absorption coefficients for a horizontal section through the sintered borosilicate glass (natural width 3.2 mm) and (right) three-dimensional representation of the subsample used for further analysis after segmentation (solid phase is white, natural width 1.47 mm). The size of the subsample is also indicated in the horizontal section.

 
All simulations were based on the binary structure of the cubic subsample shown in Fig. 1 (right). The size of the subsample was 983 voxels, or 3.18 mm3. This sample size was considered to cover a REV of the structure. This is required to identify the soil water retention characteristic as a macroscopic material property. The existence of a REV could be confirmed not only for geometrical properties, such as the porosity and the connectivity of the pore space, but also for the saturated hydraulic conductivity, which was simulated using a single-phase LB approach. With increasing sample size these parameters all converged toward a single value (Fig. 2) . Hence, a cube with a side length of about 100 voxels (i.e., 1.5 mm) was large enough to represent a REV. It was presumed that this implies the existence of a REV with respect to the soil water retention characteristic, which is a multiphase property. This hypothesis was later corroborated by our results. Similar analyses for the REV were made by Clausnitzer and Hopmans (1999) for a random packing of glass beads. They found the REV to be about five times the characteristic lengths of the structure, which agrees with our sample.



View larger version (19K):
[in this window]
[in a new window]
  		Fig. 2. Porosity (solid line), Euler number (dashed line), and hydraulic conductivity (bars) vs. sample size. The values are normalized to the result for 1503 voxels. The hydraulic conductivity was modeled for pressure gradients along the three main axes. Bars indicate the range of the results.

 
The isotropy of the porous structure within the representative sample was examined through simulation of the hydraulic conductivity in the three principal directions using the single-phase LB approach. We found deviations of at most 5.5%, which clearly demonstrates the almost perfect isotropy of the sample.

Lattice-Boltzmann Approach
The LB method is an efficient simulation tool for multiphase and multicomponent flows through porous media. There are three standard models for simulating two-phase systems: the Rothman–Keller (RK) model (Rothmann and Keller, 1988; Gunstensen and Rothman, 1991; Grunau et al., 1993), the Shan–Chen (SC) model (Shan and Chen, 1993), and the free energy (FE) approach (Swift et al., 1996). Advantages of the RK model are the local mass conservation of each phase and the possibility to adjust the surface tension and the ratio of densities and viscosities independently. An advantage of the SC model is the possible use of high density ratios, although the surface tension and the ratio of densities and viscosities cannot be adjusted independently. Also, some parameters can only be determined by numerical experiments. A major drawback of the FE approach is the unphysical non-Galilean invariance for the viscous terms in the macroscopic Navier–Stokes equation (Luo, 1998). Because for the simulation of flows in porous media the inertia forces can be neglected, we found the RK model with some modifications according to Tölke et al. (2002) the most appropriate.

D3Q19 LBE Model for Immiscible Binary Fluids
The first lattice gas model for immiscible binary fluids was proposed by Rothmann and Keller (1988), and an equivalent lattice BGK model was developed by Gunstensen and Rothman (1991). Grunau et al. (1993) modified the model for binary fluids with different density ratio and viscosities on a triangular lattice in two dimensions. We use here a three-dimensional 19 velocity LB model (D3Q19) for immiscible binary fluids developed by Tölke et al. (2002), which allows independent adjustment of the surface tension and the ratio of densities and viscosities. A method for adjusting the contact angle between the wetting and NWP was proposed by Tölke (2001). While the RK model cannot be rigorously based on thermodynamics (e.g., Pan et al., 2004), this should not be a problem as long as no major thermodynamical effects are present (e.g., phase changes). Here we are interested in an efficient numerical scheme to solve the two-phase Navier–Stokes equation. A detailed analysis of the RK model can be found in Kehrwald (2002).

Capillary Pressure–Saturation Relation
Computations of the sw,n(pc) relationship based on LB simulations can be found in Adler and Thovert (1998), Bekri and Adler (2002), Tölke et al. (2002), and Pan et al. (2004). Pan et al. (2004) calculated the capillary pressure–saturation relations using the SC model and compared them with experimental data, obtaining good results. One drawback of the SC method is that the primary physical parameters, such as the fluid–fluid and fluid–solid interaction coefficients, must be determined by model calibration using numerical experiments.

In our approach we do not use model calibration and compute the capillary pressure–saturation relationship using a system without gravity. Initially the entire pore space is filled with the wetting phase (i.e., water). To simulate the primary drainage process, we started from zero capillary pressure by imposing a constant pressure P0 at the top and bottom of the domain. We then increased the capillary pressure by decreasing the pressure at the bottom in several pressure steps Pa = P0 – {Delta}Pa. The simulation for one constant pressure step was run until the system reached a static equilibrium, thus obtaining one point on the sw,n(pc) curve. As criterion for determining that the system has reached static equilibrium we used the cumulative mass outflow during the last 100 time-steps. If the mass outflow was 1010 times smaller than the total mass of water in the sample, the criterion was assumed satisfied.

The most important parameters for two-phase systems in porous media without gravity are:

where U is the Darcy velocity of the flow, µw is the dynamic viscosity of the wetting phase, and {sigma} is the surface tension. The capillary number represents the ratio of viscous and surface tension forces. For flows in the porous medium considered here, Ca is in the range 10–7 to 10–8, indicating that the process is dominated by surface tension forces.

If {gamma} >> 1, as is the case of the water–air system considered here, the wetting phase behaves like a (moving) no-slip boundary condition for the NWP, whereas for {gamma} << 1, the wetting phase behaves like a (moving) slip boundary condition, and internal circulation is generated in the wetting phase.

The contact angle changes the geometric shape of the interface and can affect the development of dynamic process. A dynamic contact angle changes the curvature of the interface and alters the capillary forces. For our simulations we adjusted the contact angle to {alpha} = 0°, thus assuming homogeneous physicochemical properties of the solid surface.

Static capillary pressure–saturation values were determined assuming negligible influence of the capillary number and the ratio of viscosities for our type of porous medium. Since we have a very regular and connected structure, the history of the flow should not affect static equilibrium.

Hence, we adjusted the parameters to obtain maximum numerical efficiency. An example is shown in Fig. 3 , where the simulated distribution of the wetting and the NWP corresponds to {Delta}Pa = –60 hPa.



View larger version (159K):
[in this window]
[in a new window]
  		Fig. 3. Distribution of residual water (blue) and air (red) in the porous medium (gray) at the end of the drainage process simulated with the Lattice-Boltzmann approach.

 
Full-Morphology Approach
Given the complete three-dimensional structure of a porous medium, the stationary distribution of water and gas for an arbitrary capillary pressure pc can be calculated. This is true if pc is reached in a drainage process starting from complete water saturation. The FM approach was initially used by Hazlett (1995) and later also by Hilpert and Miller (2001).

At a given pc the pore space accessible for the NWP is primarily determined by the size of the pores. Within the complex structure of porous media, the hydraulic size of the pore space at a given location can be defined by the mean curvature radius of the interface between water and air that can be established at this location. The capillary pressure is identified by the principal radii of curvature of this interface, r1 and r2:

[1]
for given values of the surface tension {sigma} between the wetting and the NWP and the contact angle {alpha} between the wetting and the solid phase.

To determine the part of the pore space that will be accessible for the NWP at a given pressure during drainage, we first decompose the image with the pore radius as ordering parameter. This was done using a morphological opening, which is defined as

[2]
where X represents the pore space and B is a so-called "structuring element." In other words, OB(X) determines that part of the pore space where the structuring element fits in, and X – OB(X) denotes that part of the pore space which is smaller than B. Hence, OB(X) can be associated with the pore volume occupied by the NWP, and X – OB(X) with the pore volume occupied by the wetting phase where the size of B is associated with the capillary pressure. The shape of the structuring element B implicitly prescribes the shape of the interface. We used spherical structuring elements for the morphological opening, assuming spherical interfaces. Thus, Eq. [1] is simplified to

[3]
For a three-dimensional digital representation of the structure, a morphological opening using spherical structuring elements Br with different radii r [0, rmax] can be used to obtain a size distribution since it fulfills the relation

[4]
This approach corresponds to the concept of granulometry (Soille, 1999) and has been successfully applied to porous media (Vogel and Roth, 1998; Hilpert and Miller, 2001). For r = 0, OBr measures the total pore space. The upper limit, rmax, corresponds to the radius of the maximum sphere that fits in the pores space at some point. In between these extremities, r measures the constant curvature radius of the interface; as such it defines the capillary pressure according to Eq. [3]. This capillary pressure is associated with OBr, the fraction NWP.

As pointed out by Hilpert and Miller (2001), the choice of the structuring elements is critical for digitized data since a digital "sphere" is approximated by rectangular voxels. To respect Eq. [4], we use the straightforward definition of digital spheres given by

[5]
where c is the center point and r is the radius of the sphere. The coordinates in this approach are always located at the center point of the voxel. This method, which is in contrast to the approach of Hilpert and Miller (2001), ensures that Eq. [4] is fulfilled. To obtain a better resolution of the size distribution, we generated structuring elements with increasing r at increments of 0.5 {lambda}, where {lambda} is the length of the cubic voxel.

With decreasing r, the deviation between the digital sphere and its continuous counterpart increases substantially. Hence, for small radii r, it is not obvious how the capillary pressure is to be derived from the size of the digital sphere. We calculated the radius of the spherical structuring elements Br as the radius of the ideal sphere having the same volume as the digital representations. The latter is known for each Br from the number of voxels involved.

The steps of our algorithm to determine the quasistatic primary drainage curve sw(pc) are as follows:

  1. At the beginning, the complete pore space is water filled (i.e., pc = 0). At one side, the medium is connected to a water phase reservoir while at the opposite side it is connected to the air phase. The other four sides of the cubic domain are considered to be impermeable for water and gas.
  2. The pore space X is eroded by spheres with increasing radius r starting with the smallest (i.e., r = 1). The erosion is defined as

    [6]

  3. where Bx,r is the structuring element centered at point x.
  4. At a given capillary pressure, pores can be filled with air only if the erosion of the pore space has a continuous connection to the air reservoir. Hence, all pores are removed from the eroded pore space, which are separated from this reservoir.
  5. To determine the hydraulic size distribution of pores, the eroded set is dilated using the same structuring element as for erosion:

    [7]

  6. which finally leads to the opening of the pore space with respect to Br as defined in Eq. [2]. We identify the diameter r with the corresponding pressure according to the Young–Laplace Eq. [3].
  7. To determine the phase saturations sn and sw related to the given capillary pressure pc, only continuous pores are considered in the dilation procedure. In that case, sn is simply the volume fraction of OBr (X) related to the total pore volume, and sw = 1 – sn.
  8. We proceed with Step 2 for the next pressure step using the next larger structuring element.

As an illustration of the model results, Fig. 4 shows the calculated distribution of interfaces at various capillary pressures for the example of a single constricted pore. This approach neglects the phenomenon of entrapped water due to a discontinuous water phase. Implicitly we assume that water can always move via thin films on the solid surface. Moreover, we assume a vanishing contact angle, {alpha} = 0°.



View larger version (61K):
[in this window]
[in a new window]
  		Fig. 4. Visualization of the wetting and the nonwetting phase distribution at different capillary pressures for a single constricted pore. The connectivity is evaluated for the (left) "eroded" image, the pore-size distribution is obtained by the (middle) "opened" image, and the phase distribution is also obtained from the "opened" image while taking the (right, the air reservoir is on top) connectivity into account. Within the pore the various steps of erosions and openings are illustrated by colors from green to red with increasing radius of the structuring element.

 
It should be noted that the model produces a systematic error if the cross section of the pores is not isotropic and therefore the interface between air and water is not spherical. In this case relation Eq. [3] is not valid anymore. For an elliptical interface described by the minimum and maximum radius of curvature, r1 and r2, the critical capillary pressure for drainage is given by Eq. [1]. In the FM approach such an elliptical pore drains in several steps beginning at pressure pc = 2{sigma}/r1, which demonstrates the systematic underestimation of the pore size in this case. If the anisotropy of the interface is defined as ß = 1 – r1/r2, the underestimation of the pore size increases linearly within ß [0,1] up to a factor of 2 for ß = 1, which would correspond to the pore shape of cracks. Hilpert and Miller (2001) used the same argument to explain significant errors of the model for small water saturations where the form of the spherical interface is transformed into pendular rings. For our material, however, this error is expected to be small because of the isotropy of the sintered glass. Figure 5 is a visualization of the distribution of the wetting and the NWP in the investigated medium at three different capillary pressures.



View larger version (47K):
[in this window]
[in a new window]
  		Fig. 5. Distribution of the wetting and the nonwetting phases during the simulated drainage process using the full-morphology model. From the top to bottom the wetting phase in blue is being expelled by the brown nonwetting phase.

 
Network Model Approach
Network models are idealized representations of the complex pore geometry allowing for an efficient calculation of hydraulic properties. The crucial questions are: Which are the governing morphological characteristics that control hydraulic properties? and how do we realize these properties in a network model? As already discussed in the above section, the most relevant morphological properties are assumed to be the pore-size distribution and the three-dimensional connectivity of the pores of different size. In the following, the geometric characterization of the porous medium is restricted to these two properties.

We use the same approach as for the FM model to quantify the pore-size distribution. The volume V(r) of pores with a radius larger than r is given by the opened set OBr (X) according to Eq. [2] with respect to a spherical structuring element Br. Hence, V(r) represents the cumulative pore-size distribution. An illustration is given in Fig. 4 (middle). The smallest pore radius is determined by the image resolution, which defines the smallest structuring element according to Eq. [5].

The connectivity of the pore space is quantified by the three-dimensional Euler number {chi} as a function of the smallest pore radius {chi}(r) considered. Based on the set of pores larger than a given radius r, the Euler number in three dimensions is defined as

[8]
where N(r) is the total number of isolated objects or clusters of pores, C(r) is the total number of redundant or multiple connections (in graph theory also referred to as "genus"), and H(r) is the number of completely enclosed cavities. Considering the pore space, the latter feature H(r) would be an aggregate of solid material completely surrounded by pores. In our sample, this feature, which could be interpreted only as an artifact of image processing, does not appear. Hence, the Euler number provides a measure of connectivity having positive values for weakly connected pore structures, and decreases to negative values with increasing connectivity. Typically, the Euler number of porous media is positive when only the largest pores are considered. This means that large isolated pores exist which are connected only through narrow necks. The Euler number drops below zero around the point where the smaller pore radii of these necks are considered. The Euler number is calculated using an efficient algorithm presented by Nagel et al. (2000). The complete morphological analysis is described in detail by Vogel and Roth (1998). The results for the sintered glass material are shown in Fig. 6 (right). Isolated large pores apparently exist at r {approx} 50 to 60 µm and become interconnected at a critical pore radius of r {approx} 40 µm, where the sign of {chi} changes.



View larger version (12K):
[in this window]
[in a new window]
  		Fig. 6. Morphological material functions: (left) pore-size distribution V(r) and (right) connectivity function {chi}(r).

 
Note that measurements of V(r) and {chi}(r) can be obtained from small cutouts of the structure since volume and Euler number are additive quantities. They provide a statistical description of the pore structure and hence there is no need to analyze the complete, continuous, three-dimensional structure as for the LB and FM approaches.

The idea of our network model is to represent the morphological material functions V(r) and {chi}(r) by means of a three-dimensional network of cylindrical tubes. We use a simple model introduced by Vogel and Roth (1998), which can be adapted to the morphological material functions. The coordination number of the resulting network in this model is not a constant but is unequivocally defined by V(r) and {chi}(r). The sintered glass material used in our study could be represented by a network having a coordination number of 2.53.

To calculate the main drainage curve sw(pc) we follow exactly the same approach as for the FM model, but now applied to the reconstructed network and not to the original pore geometry. We start from water saturation while the upper surface of the network is connected to the NWP. Next, the water content for decreasing capillary pressure pc is obtained using the Young–Laplace Eq. [3]. For a given pc all pores larger than the corresponding pore radius are drained, provided that the NWP is continuous. For each pc we obtain a stationary phase distribution within the network. This leads to different points on the main drainage curve where the number of points corresponds to the number of pore-size classes present in the network.

A small cutout of 83 nodes of such a network is shown in Fig. 7 . The calculations of sw(pc) are based on network models having 323 nodes, which correspond to a volume of 15.3 mm3. This number turned out to be sufficient to obtain stable results. Actually, the resulting volume is about five times larger than the original sample. The example shown in Fig. 7 illustrates just one realization. The detailed structure of the network was subjected to a random process. To evaluate to what extent the sw(pc) is determined by the V(r) and {chi}(r) we analyzed 20 different realizations.



View larger version (141K):
[in this window]
[in a new window]
  		Fig. 7. Typical realization of the pore network model showing the distribution of air (red) and water (blue) close to the air entry point. The shown volume represents 2.8 mm3, which is approximately the same as the measured sample of sintered glass (3.18 mm3).

 

   RESULTS AND DISCUSSION
TOP
 ABSTRACT
 INTRODUCTION
 MATERIALS AND METHODS
 RESULTS AND DISCUSSION
CONCLUSIONS
 REFERENCES
 
The simulated sw(pc) relationships for primary drainage as obtained with the three different models are shown in Fig. 8 . Overall, the shapes of the different curves are very similar, showing a marked air entry pressure and a rapid drop of water saturation beyond that point. This behavior can be expected for the sintered glass, which resembles a homogeneous granular medium. The most striking difference between the models is the lower air entry pressure and the higher amount of residual water obtained by the LB approach in contrast to the two other approaches, which yield almost identical results.



View larger version (22K):
[in this window]
[in a new window]
  		Fig. 8. Capillary pressure–saturation relationship as determined with the different models: Lattice-Boltzmann model, thick dashed line; full-morphology model, two solid lines representing extreme values for simulation of drainage in six principal directions; network model, gray shaded area enclosing 20 independent realizations.

 
The difference in air entry pressure can be explained by the shape of the water–air interfaces, which are assumed to be spherical for the FM and NM models. In contrast, the shape of the interfaces generated with the LB model depends on the detailed shape of the pore structure and hence can be assumed to be more realistic. The accuracy is related to the spatial discretization of the LB model. As shown above, any deviation from the idealization of the spherical interface leads to a lower air entry pressure, which explains the observed results. The systematic error produced by the FM and NM model would be even more pronounced for planar pores, cracks, for example. To overcome this deficiency a correction factor could be introduced, which could be estimated from the detailed structure of the porous medium.

The considerable amount of residual water obtained with the LB approach is due to the limited spatial resolution of the lattice. The separation of the lattice nodes is large compared with the thickness of water layers connecting the water phase at low saturations. Hence, this connection is lost too early during the simulation of drainage, leading to trapped, residual water. For the FM and NM approaches it is implicitly assumed that the water phase is connected throughout the simulated drainage process.

The almost identical results of the FM and NM models can be explained by the fact that both methods use the same criteria for the hydraulic size of the pores (i.e., the morphological opening with respect to spherical structuring elements). However, the connectivity of the pores, which is another important aspect, is treated differently. In the FM model the spatial arrangement of pores is considered explicitly at the cost of a complete three-dimensional analysis. In contrast, the NM model evaluates the connectivity function using the Euler number and hence is based on a statistical description of the pore structure that can be obtained from local measurements. The simulated water saturation at a given capillary pressure is consistently higher for the NM than for the FM model. This small deviation can be explained by the topology of the pore structure. The network model always contains pores that are connected to other pores only by narrower necks, irrespective of the absolute size of the pores. This obviously is not the case within a real sample, where the smaller pores are always connected to a continuous larger pore, which can be understood as a result of how the sintered glass was fabricated.

The topology of the pore structure is a crucial property with respect to the effective hydraulic behavior of porous media. Using NW models, this has been demonstrated for the shape of the pressure–saturation relation (Jerauld and Salter, 1990), the relative permeability (Vogel, 2000), and residual saturations (Sok et al., 2002). Hence, a correct representation of the topology is essential when using structure models. For our sample of sintered glass, the Euler number as a function of pore size seems to provide the required information.

The three models differ significantly in their complexity, which is reflected by the required computing power. In Table 1 we compare the amount of computer memory (RAM) and the computing time required to calculate the sw(pc) curve. In contrast to the very moderate requirements of the FM and the NW, the demands of the LB approach with respect to calculation time is four to five orders higher and with respect to computer memory two to three orders of magnitude higher. This is the price for an explicit description of the multiphase dynamics. The algorithmic complexity with respect to mesh refinement is O(N3) for the NW and FM models and O(N4) for the LB method, where N is the number of voxels along one axis. Considering the results for sw(pc) this does not appear to be justifiable. However, the potential of the LB approach goes far beyond this specific application. The LB approach has the potential to provide a detailed analysis of multiphase dynamics, including imbibition, relative conductivity, hysteresis, and transient conditions. Moreover, the microscopic velocity field is a direct result of the LB approach which, therefore, can provide valuable insight for solute transport. This is in contrast to the FM, which can only treat stationary phase distributions at a given pc during primary drainage.


View this table:
[in this window]
[in a new window]
  		Table 1. Comparison of the three models investigated in terms of computational time and their potential for predicting hydraulic functions. Checks in parentheses denote "in principle, yes, but additional simplifying assumptions are required."

 
The NW model can be extended to predict absolute and relative permeabilities because the idealization of the pore structure through cylindrical tubes allows for application of a simplified model for water flow based on Hagen Poiseuille's Law (Wise, 1992; Rajaram et al., 1997; Vogel and Roth, 1998). A direct comparison with results obtained with the LB approach will be the subject of future research. The simulation of imbibition, including hysteresis and the entrapment of NWP, requires a correct description of the dynamics of individual interfaces (Blunt and Scher, 1995). For the simpler FM and NW models such a description has to include snap-off effects, which can only be introduced through simplifying rules (Ioannidis and Chatzis, 1993) (indicated by parentheses in Table 1). In contrast, the detailed dynamics of all interfaces is an immediate result of the LB approach. Other simplifying assumptions are required for all three different approaches: (i) homogeneity in the chemical properties of the solid surface (i.e., the contact angle) and (ii) the effect of structure at the subscale (i.e., surface roughness). The latter affects the microscopic movement of the wetting phase with direct consequences for the local entrapment of the NWP. Hence, also the macroscopic phase distribution is expected to be sensitive to microscopic details that are typically unknown. Also notice that imbibition is much more sensitive to the flow dynamics than drainage since the formation of entrapped NWP occurs dominantly in large pores whereas the residual wetting phase is governed by small pores and surface films. Despite these limitations, pore-scale models will continue to provide valuable insight into multiphase phenomena that are also relevant at the larger scale.

Of course it would be desirable to directly compare our simulations with corresponding experiments. Such a comparison, unfortunately, is confronted with two major hurdles. First, in contrast to simulations, the experimental measurements of sw(pc) require much larger samples, which are difficult to produce with the required homogeneity. Second, a direct comparison between experiments and our simulations requires quantification of the physicochemical properties of the solid surfaces (i.e., the water–solid contact angle). In our different simulations we assumed complete wettability. In reality, however, the contact angle is expected to be larger than zero because of construction procedures. To quantify the contact angle one possibility would be to relate the measured pore size to the actual measured air entry point, but this would correspond to fitting our simulations to the experiment. Hence, we restricted this study to comparing the different modeling approaches, which were found to produce consistent results.


   CONCLUSIONS
TOP
 ABSTRACT
 INTRODUCTION
 MATERIALS AND METHODS
 RESULTS AND DISCUSSION
 CONCLUSIONS
REFERENCES
 
We evaluated the ability of three different models to predict the capillary pressure–saturation relationships on the basis of detailed representations of the complex porous structure as measured by X-ray tomography. The models differ considerably in complexity, computational costs, and potential applications. However, the results obtained for the primary drainage curve are consistent, while differences among the models could be explained qualitatively. Our work permits the following conclusions:


   ACKNOWLEDGMENTS
 
This work was supported by the German Research Foundation (DFG).


   REFERENCES
TOP
 ABSTRACT
 INTRODUCTION
 MATERIALS AND METHODS
 RESULTS AND DISCUSSION
 CONCLUSIONS
 REFERENCES
 





This Article
Right arrow Abstract Freely available
Right arrow Figures Only
Right arrow Full Text (PDF) Free
Right arrow Alert me when this article is cited
Right arrow Alert me if a correction is posted
Services
Right arrow Similar articles in this journal
Right arrow Similar articles in ISI Web of Science
Right arrow Alert me to new issues of the journal
Right arrow Download to citation manager
Right arrow reprints & permissions
Citing Articles
Right arrow Citing Articles via ISI Web of Science (8)
Right arrow Citing Articles via Google Scholar
Google Scholar
Right arrow Articles by Vogel, H.-J.
Right arrow Articles by Roth, K.
Right arrow Search for Related Content
PubMed
Right arrow Articles by Vogel, H.-J.
Right arrow Articles by Roth, K.
GeoRef
Right arrow GeoRef Citation
Agricola
Right arrow Articles by Vogel, H.-J.
Right arrow Articles by Roth, K.
Related Collections
Right arrow Soil Physics
Right arrow Tomography
Right arrow Pore-Scale Modeling

HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS
The SCI Journals Agronomy Journal Crop Science
Journal of Natural Resources
and Life Sciences Education		 Soil Science Society of America Journal
Journal of Plant Registrations Journal of
Environmental Quality		 The Plant Genome