# Vadose Zone Journal

- Soil Science Society of America

## Abstract

This work describes an approach to porous flow modeling in which the “micro-length scale to macro-length scale” physical descriptions are addressed as Lagrangian, pore-level flow and transport. The flow features of the physical domain are solved by direct numerical simulation (DNS) with a grid-free, hybrid smoothed particle hydrodynamics (SPH) numerical method (Berry, 2002) based on a local Riemann solution. In addition to being able to handle the large deformation, fluid–fluid and fluid–solid interactions within the contorted geometries of intra- and inter-pore-scale modeling, this Riemann–SPH method should be able to simulate other complexities, such as multiple fluid phases and chemical, particulate, and microbial transport with volumetric and surface reactions. A simple model is presented for the transfer of a contaminant from a carrier fluid to solid surfaces and is demonstrated for flow in a simulated porous media.

Current textbook literature describes the phenomenological parameters and spatial averaging of material properties and velocity used to generate basic “mixture” continuum field theories for describing and simulating flow in a porous media (e.g., Slattery, 1978; Whitaker, 1999). The many successes of this approach are well documented, but relatively little discussion is devoted to those cases for which predictions with this approach fail. A growing body of evidence is now showing substantial deficiencies in the existing models and predictive capabilities for unsaturated porous flow and transport, particularly in fractured rock vadose zones (Pruess et al., 1999; Glass et al., 2002). In attempts to better describe the transport of water and contaminants through the vadose zone, specialized mixture continuum field theories have been recently extended to include preferential flow, “where water and solutes move along certain pathways, while bypassing a fraction of the porous matrix” (Hendrickx and Flury, 2001), and nonequilibrium flow, where “for various reasons, infiltrating water does not have sufficient time to equilibrate with slowly moving resident water in the bulk of the soil matrix” (Jarvis, 1998). Hendrickx and Flury (2001) discussed preferential flow mechanisms and processes at various scales, ranging from pore to areal scales, while Simunek et al. (2003) gave a review of models describing macropore and nonequilibrium flows resulting from processes at the pore and pedon scales.

The physical description of the disparate length and time scales inherent in geological multiphase fluid transport is traditionally accomplished using continuum theory. This theory views the fluids and solids as interpenetrating mixtures, each being governed by conservation laws either postulated or derived from the method of volume averaging. This Eulerian continuum approach results in field equations for the flow properties of all phases of the multiphase mixture system. As currently practiced, there are at least two major drawbacks with this approach. First, the averaging process produces a loss of information. When manipulating the averaged equations, covariance terms are typically neglected because pore- and fracture-scale models are needed for their construction. The resulting equation system is not complete because the pore- and fracture-scale information is filtered in the averaging process. If these covariance terms are to be included directly, suitable constitutive relations must be developed, in their stead, to obtain closure. Second, averaging leads to unknown terms representing the interactions between the phases (including the solid phase). These terms must be modeled to close the description of the system. The detailed nature of the interactions between solid and fluid phases cannot be understood by the application of mixture theory alone. It is the pore- and fracture-scale physics and the phasic interactions that determine the nature of these constitutive relations. This important information is seldom used, however, because of the difficulty in obtaining it.

As a first step toward an increased understanding of the pore- and fracture-scale physical behavior and phasic interactions, we describe a DNS approach with which a Lagrangian description of fluid flows can be simulated through a broad range of representative porous and/or fracture configurations. The DNS approach provides the capability of studying the nonlinear and geometrically complicated phenomena of multiphase inter- and intra-pore flow wherein preferential and nonequilibrium flow and transport are accounted for as an integral part of the physical description. We present a novel approach that uses a grid-free numerical method to simulate a simplified Lagrangian, microscopic flow. Our approach, however, should extend to include multiple fluid phases with surface tension, liquid bridging, capillary attachment–detachment, chemical constituents with volumetric and surface reactions and pore-surface adsorption–desorption and transport, suspended particulate transport and deposition, and microbial transport, deposition, and buildup. Here we modify the SPH grid-free method, which employs no permanent connectivity between computational nodes or particles, to incorporate local Riemann solutions, to produce our Riemann–SPH grid-free method. It should apply to problems where traditional Eulerian and Lagrangian techniques fail, such as those involving large deformation, fluid–fluid or fluid–solid interactions, and contorted geometries, precisely the situation present in detailed, intra- and inter-pore-scale (and fracture-scale) modeling.

## RIEMANN–SPH METHOD

The SPH method is a Lagrangian, particle based approach to hydrodynamics. Unlike finite-difference, finite-element, and finite-volume methods, which utilize grid cells to cover the spatial domain of interest, particle methods use computational points which are moving in space where computational data in the material is sampled. The motion of the particles is usually governed by the equations of the material. The SPH method is Lagrangian in the hydrodynamics sense in that the particles follow the material (fluid and perhaps microbes and solids) flow. There is no need for a computational grid because the particles carry all of the computational information, and the particles themselves form the computational framework on which the material flow equations are solved. Originally introduced in the context of astrophysics (Lucy, 1977; Gingold and Monaghan, 1977), the SPH method has evolved into multiple variations with a plethora of applications ranging from high-strain-rate solid deformation mechanics to ice-jammed river flows. This approach offers a substantial advantage for problems that contain moving discontinuities and for those characterized by large deformations as would be encountered in a DNS simulation of saturated and unsaturated, inter- and intra-pore flows. The Lagrangian feature of SPH may also prove ideal for chemical species and biological organism transport, and microscopic surface transport, as well as chemical reaction and biofilm interaction and growth.

With SPH, a moving continuum field is simulated by a large group of moving computational particles. We then replace the point mechanical and material thermomechanical variables (e.g., fluid velocity or fluid pressure at a specific point), which vary rapidly on a scale comparable with the particle spacing, by smoothed variables obtained by averaging over regions that are large compared with the particle spacing but small compared with the complete system. Intuitively, a continuum is modeled as a system of particles where, in the limit, the particles' sizes vanish and the distances between adjacent particles approach zero. In this limit, the system of particles approaches a deformable continuum, that is, an integral over a continuous body that may be considered as the limiting extension of a summation over a corresponding system of particles. A formal mathematical definition of local mean variables is used to translate the point equations for the material flow into computationally attractive representations of the conservation–balance equations.

The basis for SPH is a regularization or smoothing of any point property of the material, *f*. The local mean value is defined by1

The weighting or smoothing kernel function *w* is usually a centrally peaked, spherically symmetric function, with compact support radius *h* (the smoothing length) representing the effective width of the kernel. The properties of the weighting function are further specified to ensure that the ordinary equations of vector calculus can be performed on the local mean variables, and also that certain other mathematical manipulations, such as differentiating under an integral, can be used. Typical kernels have a Gaussian-like shape, continuous derivatives of specified order, and compact support so that only the volume within a constant multiple of *h* contributes to the integration. The SPH method is usually viewed as an approximation method for information at unordered points, but it can also be loosely interpreted as a Monte Carlo technique for evaluating Eq. [1] [in which *w* is a probability distribution function]. A Monte Carlo summation approximation of Eq. [1] is 2where the factor *m _{j}*/ρ

*is the volume element that SPH associates with particle*

_{j}*j*, and

*n*is the number of particles within the influence range of the kernel. The spatial gradient can be directly calculated without a computational grid by appropriate differentiation of the approximation formula, Eq. [2], producing 3Divergence of a vector function is found similarly.

The basics of the SPH method are well described in Monaghan (1992). For ease of exposition, the Euler equations with a barotropic equation of state (pressure dependence only on density) will be used to illustrate development of the SPH method on realistic, fluid-dynamic field equations and to show the hybrid extension to incorporate a Riemann solution method. For the porous flow demonstration that follows, viscosity terms are added to produce the Navier–Stokes equations governing viscous Newtonian fluid dynamic motions. A variety of SPH approximations of the inviscid, isothermal hydrodynamic equations (Euler equations) have appeared in the literature. For example, 4and, for the two-dimensional analysis to be demonstrated here, the weighting kernel (Wingate and Fisher, 1993) is selected. These equations represent a system of ordinary differential equations for each particle that are time-integrated numerically with the explicit-Euler method. Time step size is adjusted dynamically to satisfy the minimum Courant condition, Δ*t* ≤ 0.4, where *c* is the sound speed.

Because the basic SPH method exhibits what has become known in the literature as “tensile instability” (can also occur when pressure is positive) (Swegle et al., 1995), and because these compressible formulations also require an artificial viscosity for stabilization, a hybrid SPH method was developed which incorporates a Riemann solution technique. Following the work of Parshikov et al. (2000), for high-velocity impact of solids, we insert (into SPH) approximations for the velocities and pressures obtained from Riemann solutions to break up discontinuities between interacting particles (Berry, 2002). No artificial viscosity is explicitly added, and tensile instability appears to be dramatically reduced or eliminated. In this hybrid approach, the interaction of two particles *i* and *j* is taken to be equivalent to that of the acoustic Riemann solution that results from assuming an initial discontinuity between *i* and *j*. The mass continuity and momentum balance Eq. [4] become 5and *u*^{*R}_{ij} and *p*^{*}_{ij} are normal velocities and pressures resulting from the Riemann solution between particles *i* and *j*. Furthermore, the solution between particles *i* and *j* is allowed to mimic fluid cavitation or fracture behavior and thereby produce dual values of velocity corresponding to a specified cavitation or fracture pressure, or tension cut-off value (space limitations preclude further discussion here) (Stanyukovich, 1960).

The SPH method described is a compressible method. An approach to represent nearly incompressible flow is given by Rogers et al. (1987). We pick, for example, a stiff, barotropic equation of state such as6where ρ_{0} is the reference density and *p*_{0} is the magnitude of the pressure. If these parameters are properly picked to yield (in the simulations) density variations of 0.25 to 1% with corresponding Mach numbers of around 0.1 and a Courant number of about 0.1, then the flow can be regarded as approximately incompressible. This however, may not always be an adequate way to produce the incompressible constraint. For such cases, the method should be modified to update pressure semi-implicitly in a manner consistent with other incompressible fluid dynamic algorithms.

## CONTAMINANT MODEL

A contaminant transport and interaction model is developed here to show the rich dynamical behavior that can result from even the simplest of physical assumptions. It is assumed that a contaminant may be present in a carrier liquid (water) as it flows, and that this action occurs solely by advection (i.e., without diffusion). All solid surfaces are assumed to have a much higher affinity for the contaminant than does the carrier liquid. Furthermore, they have the ability to pull the contaminant out of the carrier liquid if the liquid is within a prescribed distance, from solid surfaces (thus making contaminant pick-up a localized phenomena). This model is simplified by assuming that, when contaminated fluid is within this prescribed distance of a solid surface, all of that carrier liquid's contaminant is transferred to the solid on a time scale that is very small compared with the transport time scale, and that the carrier liquid is then contaminant-free or clean. For the demonstration example that follows, no limits are placed on how much contaminant may be accumulated by a solid surface and it is assumed that the solid's affinity for the contaminant is so high that it cannot transfer back to the carrier liquid (one-way transfer). It is assumed that contaminated carrier fluid has a contaminant concentration value of 1.0 and that the cleaned carrier fluid has a contaminant concentration of 0.0. Thus, solid surfaces continually accumulate contaminant in unit increments as carrier liquid particles approach solid particles (surfaces) within the prescribed distance for contaminant transference.

## DEMONSTRATION EXAMPLE

We illustrate this Riemann–SPH method with an example from Berry (2002). Utilizing a particle neighbor-tracking method, edge detection algorithm (Berry, 2001), and basic particle code structures (again space limitations preclude further details here), it is relatively easy to add increasingly more complex physical behaviors. We have constructed a basic model for unsaturated liquid flow through a high-porosity media with a small particulate contaminant transport and deposition model for mock subsurface flows. Though easily extendable to three dimensions, as a first attempt at DNS of porous flows with the Riemann–SPH method, this initial effort solves the Navier–Stokes equations in two dimensions with only one fluid phase and no surface tension. However, it clearly shows the merit of this approach. Because we are ultimately interested in an increased understanding of the pore- and fracture-scale physical behaviors we pick a representative-scale problem to demonstrate.

This simulation models a column of water, initially 1.2 mm deep with hydrostatic pressure head and contaminated with normalized concentration level of 1.0, placed above a column of representative porous material, also 1.2 mm high by 1.0 mm wide but uncontaminated with a normalized value of 0.0. The water is allowed to flow through the porous column under the influence of gravity, and out the bottom, where it falls into a catch-pan and accumulates. The water is represented with 12000 particles, while the solids of the column walls, catch-pan, and the porous region are represented with 5089 particles. The sticking distance for the contaminant to transfer to the solid surfaces is taken to be 0.015 mm. The degree of flow saturation is not only time dependent, but is spatially varying as well. Figure 1 shows the spatial variation of velocity magnitude within the porous column at approximate times of 0.117, 0.230, and 0.459 s. Figure 2 shows the contaminant distribution at the same respective times as those of Fig. 1. Mechanical dispersion effects can be readily observed in these Figures. From Fig. 2c, which shows the distribution of contaminant when the porous column has nearly drained, it is apparent that a significant portion of the contaminant has been sorbed onto the porous solid particles, thereby reducing its presence in the pan's water. It is noted for Fig. 2 that the color scale of the contaminant level indicates normalized values of 0.0 (uncontaminated, shown as blue), 1.0 (contaminated, shown black), and 3.0 (contaminated with normalized values of 3.0 or higher, shown as magenta), while the contaminant level on some of the porous solid particles approaches normalized values of 150. (Normalized values of 3 or 150 indicate contaminant concentrations 3 times or 150 times, respectively, higher than the initial contaminations concentration level in the water column.) This demonstration problem was executed for approximately 5.5 h on a PC Athlon 1.3-GHz computer, completing roughly 40000 time steps. Though the time step sizes were dynamically adjusted, typical time steps were on the order of 12 μs.

A simulation of this type would be very difficult or impossible using traditional numerical methods. To apply an Eulerian approach, a fixed, boundary-fitted grid would first have to be set up, sophisticated boundary conditions applied, and a very good interface-tracking algorithm incorporated. In addition to setting up the initial, boundary-fitted grid, because of the contorted flow paths, a Lagrangian grid-based method would require extensive remeshing, which is computationally expensive and introduces significant error. The simple, history-based, contaminant transport–deposition model exhibited here would be problematic for an Eulerian, grid-based solution method. Whenever fluid or contaminant behavior depends on the history of what that fluid or contaminant has seen in its dynamical evolution, Eulerian methods are at a significant disadvantage.

## CONCLUSIONS

The simple demonstration example presented here shows the benefits of pursuing the Lagrangian, direct numerical simulation approach with the Riemann–SPH method. We have presented an approach for simulating the physics and chemistry of microscopic flow in porous and/or fractured media. Furthermore, we advocate this approach as a means to building the macroscopic averaged flow models and the appropriate constitutive relations necessary to provide the closure for these macroscopic models. A simple, dynamical contaminant model was shown to produce a result rich in structure that would be difficult to obtain with conventional methods. To extend its applicability to the complexities of multiphase flow, we recommend modification of the method to treat pressure implicitly (in time), addition of another set of particles to represent the gaseous phase, and the incorporation surface tension. We further recommend improvement of the models to incorporate diffusion and chemical reactions, and potentially solid matrix compressibility, deformation, and erosion.

## Acknowledgments

The authors would like to acknowledge the U.S. Department of Energy, Office of Environmental Management, for funding this research as part of the Environmental Systems Research Program during Fiscal-Year 2001. This research was conducted under DOE Idaho Operations Office Contract DE-AC07-99ID13727.

- Received May 14, 2003.