A practical model for mobile, residual, and entrapped NAPL in water-wet porous media

Flow of nonvolatile nonaqueous phase liquid (NAPL) and aqueous phases that account for mobile, entrapped, and residual NAPL in variably saturated water-wet porous media is modeled and compared against results from detailed laboratory experiments. Residual saturation formation in the vadose zone is a process that is often ignored in multifluid flow simulators, which might cause an overestimation of the volume of NAPL that reaches the ground water. Mobile NAPL is defined as being continuous in the pore space and flows under a pressure gradient or gravitational body force. Entrapped NAPL is defined as being occluded by the aqueous phase, occurring as immobile ganglia surrounded by aqueous phase in the pore space and formed when NAPL is replaced by the aqueous phase. Residual NAPL is defined as immobile, nonwater entrapped NAPL that does not drain from the pore spaces and is conceptualized as being either continuous or discontinuous. Free NAPL comprises mobile and residual NAPL. The numerical model is formulated on mass conservation equations for oil and water, transported via NAPL and aqueous phases through variably saturated porous media. To account for phase transitions, a primary variable switching scheme is implemented for the oil-mass conservation equation over three phase conditions: (1) aqueous or aqueous-gas with dissolved oil, (2) aqueous or aqueous-gas with entrapped NAPL, and (3) aqueous or aqueous gas with free NAPL. Two laboratory-scale column experiments are modeled to verify the numerical model. Comparisons between the numerical simulations and experiments demonstrate the necessity to include the residual NAPL formation process in multifluid flow simulators.     

Mixed finite element methods and higher order temporal approximations for variably saturated groundwater flow

Richards’ equation (RE) is commonly used to model flow in variably saturated porous media. However, its solution continues to be difficult for many conditions of practical interest. Among the various time discretizations applied to RE, the method of lines (MOL) has been used successfully to introduce robust, accurate, and efficient temporal approximations. At the same time, a mixed-hybrid finite element method combined with an adaptive, higher order time discretization has shown benefits over traditional, lower order temporal approximations for modeling single-phase groundwater flow in heterogeneous porous media. Here, we extend earlier work for single-phase flow and consider two mixed finite element methods that have been used previously to solve RE using lower order time discretizations with either fixed time steps or empirically based adaption. We formulate the two spatial discretizations within a MOL context for the pressure head form of RE as well as a fully mass-conservative version. We conduct several numerical experiments for both spatial discretizations with each formulation, and we compare the higher order, adaptive time discretization to a first-order approximation with formal error control and adaptive time step selection. Based on the numerical results, we evaluate the performance of the methods for robustness and efficiency.     

Modeling Variably Saturated Multispecies Reactive Groundwater Solute Transport with MODFLOW‐UZF and RT3D

A numerical model was developed that is capable of simulating multispecies reactive solute transport in variably saturated porous media. This model consists of a modified version of the reactive transport model RT3D (Reactive Transport in 3 Dimensions) that is linked to the Unsaturated‐Zone Flow (UZF1) package and MODFLOW. Referred to as UZF‐RT3D, the model is tested against published analytical benchmarks as well as other published contaminant transport models, including HYDRUS‐1D, VS2DT, and SUTRA, and the coupled flow and transport modeling system of CATHY and TRAN3D. Comparisons in one‐dimensional, two‐dimensional, and three‐dimensional variably saturated systems are explored. While several test cases are included to verify the correct implementation of variably saturated transport in UZF‐RT3D, other cases are included to demonstrate the usefulness of the code in terms of model run‐time and handling the reaction kinetics of multiple interacting species in variably saturated subsurface systems. As UZF1 relies on a kinematic‐wave approximation for unsaturated flow that neglects the diffusive terms in Richards equation, UZF‐RT3D can be used for large‐scale aquifer systems for which the UZF1 formulation is reasonable, that is, capillary‐pressure gradients can be neglected and soil parameters can be treated as homogeneous. Decreased model run‐time and the ability to include site‐specific chemical species and chemical reactions make UZF‐RT3D an attractive model for efficient simulation of multispecies reactive transport in variably saturated large‐scale subsurface systems.     

Numerical modelling of multiphase flow in porous media

The simultaneous flow of immiscible fluids in porous media occurs in a wide variety of applications. The equations governing these flows are inherently nonlinear, and the geometries and material properties characterizing many problems in petroleum and groundwater engineering can be quite irregular. As a result, numerical simulation often offers the only viable approach to the mathematical modelling of multiphase flows. This paper provides an overview of the types of models that are used in this field and highlights some of the numerical techniques that have appeared recently. The exposition includes discussions of multiphase, multispecies flows in which chemical transport and interphase mass transfers play important roles. The paper also examines some of the outstanding physical and mathematical problems in multiphase flow simulation. The scope of the paper is limited to isothermal flows in natural porous media; however, many of the special techniques and difficulties discussed also arise in artificial porous media and multiphase flows with thermal effects.     

Visualization and modeling of the colonization dynamics of a bioluminescent bacterium in variably saturated,translucent quartz sand

An experimental and numerical investigation was conducted to study the colonization dynamics of a bioluminescent bacterium, Pseudomonas fluorescens HK44, during growth in a porous medium under steady, variably saturated flow conditions. Experiments were conducted in a thin-slab light transmission chamber filled with uniform, translucent quartz sand. Steady, variably saturated flow conditions were established using drip emitters mounted on the top of the chamber, with glucose applied through a central dripper located directly above an inoculated region of the chamber. Periodic pulses of salicylate and a dye tracer were applied to induce bioluminescence of the bacterium to monitor colony expansion and to track changes in the hydraulic and transport properties of the sand. Changes in the apparent water saturation of the sand were quantified by monitoring light transmission through the chamber with a CCD camera. The colonized region expanded laterally by about 15 cm, and upward against the flow by 7–8 cm during the 6-day experiment while apparent saturations in the colonized region decreased by 7–9% and the capillary fringe dropped by ∼5 cm. The observed data were reproduced approximately using a numerical model that accounted for the processes of water flow, solute and bacterial transport, cell growth and accumulation, glucose and oxygen consumption, and gas diffusion and exchange. The results of this study illustrate some of the complexities associated with coupled flow, reactive transport, and biological processes in variably saturated porous media, such as localized desaturation, capillary fringe lowering effects, and upstream movement of bacterial colonization, that may not readily observable using other experimental techniques.     

FracKfinder: A MATLAB Toolbox for Computing Three-Dimensional Hydraulic Conductivity Tensors for Fractured Porous Media

Fractures in porous media have been documented extensively. However, they are often omitted from groundwater flow and mass transport models due to a lack of data on fracture hydraulic properties and the computational burden of simulating fractures explicitly in large model domains. We present a MATLAB toolbox, FracKfinder, that automates HydroGeoSphere (HGS), a variably saturated, control volume finite-element model, to simulate an ensemble of discrete fracture network (DFN) flow experiments on a single cubic model mesh containing a stochastically generated fracture network. Because DFN simulations in HGS can simulate flow in both a porous media and a fracture domain, this toolbox computes tensors for both the matrix and fractures of a porous medium. Each model in the ensemble represents a different orientation of the hydraulic gradient, thus minimizing the likelihood that a single hydraulic gradient orientation will dominate the tensor computation. Linear regression on matrices containing the computed three-dimensional hydraulic conductivity (K) values from each rotation of the hydraulic gradient is used to compute the K tensors. This approach shows that the hydraulic behavior of fracture networks can be simulated where fracture hydraulic data are limited. Simulation of a bromide tracer experiment using K tensors computed with FracKfinder in HGS demonstrates good agreement with a previous large-column, laboratory study. The toolbox provides a potential pathway to upscale groundwater flow and mass transport processes in fractured media to larger scales.     

Transport and fate of nonaqueous phase liquid (NAPL) in variably saturated porous media with evolving scales of heterogeneity

 A stochastic simulation is performed to study multiphase flow and contaminant transport in fractal porous media with evolving scales of heterogeneity. Numerical simulations of residual NAPL mass transfer and subsequent transport of dissolved and/or volatilized NAPL mass in variably saturated media are carried out in conjunction with Monte Carlo techniques. The impact of fractal dimension, plume scale and anisotropy (stratification) of fractal media on relative dispersivities is investigated and discussed. The results indicate the significance of evolving scale of porous media heterogeneity to the NAPL transport in the subsurface. In general, the fractal porous media enhance the dispersivities of NAPL mass plume transport in both the water phase and the gas phase while the influence on the water phase is more significant. The porous media with larger fractal dimension have larger relative dispersivities. The aqueous horizontal dispersivity exhibits a most significant increase against the plume scale.     

Necessary conditions for inverse modeling of flow through variably saturated porous media

Non-unique solutions of inverse problems arise from a lack of information that satisfies necessary conditions for the problem to be well defined. This paper investigates these conditions for inverse modeling of water flow through multi-dimensional variably saturated porous media. It shows that in order to obtain a unique estimate of hydraulic parameters, along each streamline of the flow field (1) spatial and temporal head observations must be given; (2) the number of spatial and temporal head observations required should be greater or equal to the number of unknown parameters; (3) the flux boundary condition or the pumping rate of a well must be specified for the homogeneous case and both boundary flux and pumping rate are a must for the heterogeneous case; (4) head observations must encompass both saturated and unsaturated conditions, and the functional relationships for unsaturated hydraulic conductivity/pressure head and for the moisture retention should be given, and (5) the residual water content value also need to be specified a priori or water content measurements are needed for the estimation of the saturated water content.For field problems, these necessary conditions can be collected or estimated but likely involve uncertainty. While the problems become well defined and have unique solutions, the solutions likely will be uncertain. Because of this uncertainty, stochastic approaches are deemed to be appropriate for inverse problems as they are for forward problems to address uncertainty. Nevertheless, knowledge of these necessary conditions is critical to reduce uncertainty in both characterization of the vadose zone and the aquifer, and prediction of water flow and solute migration in the subsurface.     

Simulating a lake as a high-conductivity variably saturated porous medium

One approach for simulating ground water–lake interactions is to incorporate the lake into the ground water solution domain as a high-conductivity region. Previous studies have developed this approach using fully saturated models. This study extends this approach to variably saturated models, so that ground water–lake interactions may be more easily simulated with commonly used or public domain variably saturated codes that do not explicitly support coupled lake–water balance modeling. General guidelines are developed for the choices of saturated hydraulic conductivity and moisture retention and relative permeability curves for the lake region. When applied to an example ground water–lake system, model results are very similar to those from a model in which the lake is represented as a specified head boundary continuously updated by a lake mass balance. The high-conductivity region approach is most suitable for relatively simple geometries and lakes with slower and smaller fluctuations when the overall flow pattern and system fluxes, rather than the detailed flow pattern around the intersection of the lake and land surfaces, are of interest.     

A two‐dimensional model of hydraulic performance of stormwater infiltration systems

Stormwater infiltration systems are a popular method for urban stormwater control. They are often designed using an assumption of one‐dimensional saturated outflow, although this is not very accurate for many typical designs where two‐dimensional (2D) flows into unsaturated soils occur. Available 2D variably saturated flow models are not commonly used for design because of their complexity and difficulties with the required boundary conditions. A purpose‐built stormwater infiltration system model was thus developed for the simulation of 2D flow from a porous storage. The model combines a soil moisture–based model for unsaturated soils with a ponded storage model and uses a wetting front‐tracking approach for saturated flows. The model represents the main physical processes while minimizing input data requirements. The model was calibrated and validated using data from laboratory 2D stormwater infiltration trench experiments. Calibrations were undertaken using five different combinations of calibration data to examine calibration data requirements. It was found that storage water levels could be satisfactorily predicted using parameters calibrated with either data from laboratory soils tests or observed water level data, whereas the prediction of soil moistures was improved through the addition of observed soil moisture data to the calibration data set. Copyright © 2012 John Wiley & Sons, Ltd.     

Laboratory Study of Air Sparging: Air Flow Visualization

Laboratory flow visualization experiments, using glass beads as the porous medium, were conducted to study air sparging, an innovative technology for subsurface contaminant remediation. The purpose of these experiments was to observe how air flows through saturated porous media and to obtain a basic understanding of air plume formation and medium heterogeneity effects. The experiments indicate that air flow occurring in discrete, stable channels is the most probable flow behavior in medium to fine grained water saturated porous media and that medium heterogeneity plays an important role in the development of air channels. Several simulated scales of heterogeneities, from pore to field, have been studied. The results suggest that air channel formation is sensitive to the various scales of heterogeneities. Site-specific hydrogeologic settings have to be carefully reviewed before air sparging is applied to remediate sites contaminated by volatile organic compounds.     

New preconditioning strategy for Jacobian-free solvers for variably saturated flows with Richards’ equation

We develop a new approach for solving the nonlinear Richards’ equation arising in variably saturated flow modeling. The growing complexity of geometric models for simulation of subsurface flows leads to the necessity of using unstructured meshes and advanced discretization methods. Typically, a numerical solution is obtained by first discretizing PDEs and then solving the resulting system of nonlinear discrete equations with a Newton-Raphson-type method. Efficiency and robustness of the existing solvers rely on many factors, including an empiric quality control of intermediate iterates, complexity of the employed discretization method and a customized preconditioner. We propose and analyze a new preconditioning strategy that is based on a stable discretization of the continuum Jacobian. We will show with numerical experiments for challenging problems in subsurface hydrology that this new preconditioner improves convergence of the existing Jacobian-free solvers 3-20 times. We also show that the Picard method with this preconditioner becomes a more efficient nonlinear solver than a few widely used Jacobian-free solvers.     

Flow in unsaturated random porous media,nonlinear numerical analysis and comparison to analytical stochastic models

This work presents a rigorous numerical validation of analytical stochastic models of steady state unsaturated flow in heterogeneous porous media. It also provides a crucial link between stochastic theory based on simplifying assumptions and empirical field and simulation evidence of variably saturated flow in actual or realistic hypothetical heterogeneous porous media. Statistical properties of unsaturated hydraulic conductivity, soil water tension, and soil water flux in heterogeneous soils are investigated through high resolution Monte Carlo simulations of a wide range of steady state flow problems in a quasi-unbounded domain. In agreement with assumptions in analytical stochastic models of unsaturated flow, hydraulic conductivity and soil water tension are found to be lognormally and normally distributed, respectively. In contrast, simulations indicate that in moderate to strong variable conductivity fields, longitudinal flux is highly skewed. Transverse flux distributions are leptokurtic. the moments of the probability distributions obtained from Monte Carlo simulations are compared to modified first-order analytical models. Under moderate to strong heterogeneous soil flux conditions (σ²_y≥1), analytical solutions overestimate variability in soil water tension by up to 40% as soil heterogeneity increases, and underestimate variability of both flux components by up to a factor 5. Theoretically predicted model (cross-)covariance agree well with the numerical sample (cross-)covarianaces. Statistical moments are shown to be consistent with observed physical characteristics of unsaturated flow in heterogeneous soils.©1998 Elsevier Science Limited. All rights reserved     

A new saturated/unsaturated model for stormwater infiltration systems

Infiltration systems are widely used as an effective urban stormwater control measure. Most design methods and models roughly approximate the complex physical flow processes in these systems using empirical equations and fixed infiltration rates to calculate emptying times from full. Sophisticated variably saturated flow models are available, but rarely applied owing to their complexity. This paper describes the development and testing of an integrated one‐dimensional model of flow through the porous storage of a typical infiltration system and surrounding soils. The model accounts for the depth in the storage, surrounding soil moisture conditions and the interaction between the storage and surrounding soil. It is a front‐tracking model that innovatively combines a soil‐moisture‐based solution of Richard's equation for unsaturated flow with piston flow through a saturated zone as well as a reservoir equation for flow through a porous storage. This allows the use of a simple non‐iterative numerical solution that can handle ponded infiltration into dry soils. The model is more rigorous than approximate stormwater infiltration system models and could therefore be valuable in everyday practice. A range of test cases commonly used to test soil water flow models for infiltration in unsaturated conditions, drainage from saturation and infiltration under ponded conditions were used to test the model along with an experiment with variable depth in a porous storage over saturated conditions. Results show that the model produces a good fit to the observed data, analytical solutions and Hydrus. Copyright © 2008 John Wiley & Sons, Ltd.     

Effects of chemical reactions on density-dependent fluid flow: on the numerical formulation and the development of instabilities

A three-dimensional, reactive numerical flow model is developed that couples chemical reactions with density-dependent mass transport and fluid flow. The model includes equilibrium reactions for the aqueous species, kinetic reactions between the solid and aqueous phases, and full coupling of porosity and permeability changes that result from precipitation and dissolution reactions in porous media. A one-step, global implicit approach is used to solve the coupled flow, transport and reaction equations with a fully implicit upstream-weighted control volume discretization. The Newton–Raphson method is applied to the discretized non-linear equations and a block ILU-preconditioned CGSTAB method is used to solve the resulting Jacobian matrix equations. This approach permits the solution of the complete set of governing equations for both concentration and pressure simultaneously affected by chemical and physical processes. A series of chemical transport simulations are conducted to investigate coupled processes of reactive chemical transport and density-dependent flow and their subsequent impact on the development of preferential flow paths in porous media. The coupled effects of the processes driving flow and the chemical reactions occurring during solute transport is studied using a carbonate system in fully saturated porous media. Results demonstrate that instability development is sensitive to the initial perturbation caused by density differences between the solute plume and the ambient groundwater. If the initial perturbation is large, then it acts as a “trigger” in the flow system that causes instabilities to develop in a planar reaction front. When permeability changes occur due to dissolution reactions occurring in the porous media, a reactive feedback loop is created by calcite dissolution and the mixed convective transport of the system. Although the feedback loop does not have a significant impact on plume shape, complex concentration distributions develop as a result of the instabilities generated in the flow system.     

Velocity variations in carbonate rocks due to dual porosity and wave-induced fluid flow

Stiffness variations in carbonates may be described as resulting from different concentrations of flat compliant pores or cracks, which can have a significant effect on the effective stiffness and acoustic properties (e.g., velocities and attenuations) of dry as well as saturated carbonates, although they carry extremely little porosity. As shown in this paper, the effects of dual porosity and wave-induced fluid flow or pore pressure communication may also play a significant role. On the basis of a previously published T-matrix approach to model the effective viscoelastic properties of cracked porous media, we illustrate the (frequency-dependent) effects of wave-induced fluid flow (mainly squirt flow) or pore pressure communication for a model structure consisting of a mixture of fluid-saturated porous grains and fluid-saturated cavities (vugs, etc.) that are embedded in a solid matrix associated with carbonates. We assume that the pores within the porous grains are decoupled from the pores in the solid matrix (and possibly saturated with different fluids) but that each pore system at the micro and/or mesoscale may or may not be connected. For each of four different connectivity models, we present numerical results for four different cases of microstructure (that emphasize the importance of cracks and flat compliant pores). Our numerical results indicate that the velocity and attenuation spectra of carbonates vary significantly, even when the crack density and all other volume concentrations are constant.     

VS2DRTI: Simulating Heat and Reactive Solute Transport in Variably Saturated Porous Media

Variably saturated groundwater flow, heat transport, and solute transport are important processes in environmental phenomena, such as the natural evolution of water chemistry of aquifers and streams, the storage of radioactive waste in a geologic repository, the contamination of water resources from acid‐rock drainage, and the geologic sequestration of carbon dioxide. Up to now, our ability to simulate these processes simultaneously with fully coupled reactive transport models has been limited to complex and often difficult‐to‐use models. To address the need for a simple and easy‐to‐use model, the VS2DRTI software package has been developed for simulating water flow, heat transport, and reactive solute transport through variably saturated porous media. The underlying numerical model, VS2DRT, was created by coupling the flow and transport capabilities of the VS2DT and VS2DH models with the equilibrium and kinetic reaction capabilities of PhreeqcRM. Flow capabilities include two‐dimensional, constant‐density, variably saturated flow; transport capabilities include both heat and multicomponent solute transport; and the reaction capabilities are a complete implementation of geochemical reactions of PHREEQC. The graphical user interface includes a preprocessor for building simulations and a postprocessor for visual display of simulation results. To demonstrate the simulation of multiple processes, the model is applied to a hypothetical example of injection of heated waste water to an aquifer with temperature‐dependent cation exchange. VS2DRTI is freely available public domain software.     

Modeling of non-isothermal multi-component reactive transport in field scale porous media flow systems

A general 2-D finite element multi-component reactive transport code, TRANQUI, was developed, using a sequential iteration approach (SIA). It is well suited to deal with complex real-world thermo-hydro-geochemical problems for single-phase variably water saturated porous media flow systems. The model considers a wide range of hydrological and thermodynamic as well as chemical processes such as aqueous complexation, acid-base, redox, mineral dissolution/precipitation, gas dissolution/ex-solution, ion exchange and adsorption via surface complexation. Under unsaturated conditions only water flow is considered, although gas pressures are allowed to vary in space in a depth-dependent manner specified by the user. In addition to the fully iterative sequential approach (SIA), a sequential non-iterative approach (SNIA), in which transport and chemistry are de-coupled, was implemented and tested. The accuracy and numerical performance of SIA and SNIA have been compared using several test cases. The accuracy of SNIA depends on space and time discretization as well as on the nature of the chemical reactions. The capability of the code to model a real case study in the field is illustrated by its application to the modeling of the hydrochemical evolution of the Llobregat Delta aquitard in northeastern Spain over the last 3500 years during when fresh-water flow from a lower aquifer displaced the native saline aquitard waters. Manzano and Custodio carried out a reactive transport model of this case study by using the PHREEQM code and considering water flow, aqueous complexation, cation exchange and calcite dissolution. Their results compare favorably well with measured porewater chemical data, except for some of the cations. Our code is not only able to reproduce the results of previous numerical models, but leads to computed concentrations which are closer to measured data mainly because our model takes into consideration redox processes in addition to the processes mentioned above. A number of sensitivity runs were performed with TRANQUI in order to analyze the effect of errors and uncertainties on cation selectivities.     

Spatial and temporal analysis of groundwater recharge with application to sampling design

A numerical experiment of flow in variably saturated porous media was performed in order to evaluate the spatial and temporal distribution of the groundwater recharge at the phreatic surface for a shallow aquifer as a function of the input rainfall process and soil heterogeneity. The study focused on the groundwater recharge which resulted from the percolation of the excess rainfall for a 90-days period of an actual precipitation record. Groundwater recharge was defined as the water flux across the moving phreatic surface. The observed spatial non-uniformity of the groundwater recharge was caused by soil heterogeneity and is particularly pronounced during the stage of recharge peak (substantial percolation stage). During that stage the recharge is associated with preferential flow paths defined as soil zones of locally higher hydraulic conductivity. For the periods of low percolation intensity the groundwater recharge was exhibiting more uniform spatial characteristics. The temporal distribution of the recharge was found to be a function of the frequency and intensity of the rainfall events. Application of sampling design demonstrates the joint influence of the spatial and temporal recharge variability on the cost-effective monitoring of groundwater potentiometric surfaces.     

Coupling multiscale finite element method for consolidation analysis of heterogeneous saturated porous media

The multiscale finite element method is developed for solving the coupling problems of consolidation of heterogeneous saturated porous media under external loading conditions. Two sets of multiscale base functions are constructed, respectively, for the pressure field of fluid flow and the displacement field of solid skeleton. The coupling problems are then solved with a multiscale numerical procedure in space and time domain. The heterogeneities induced by permeabilities and mechanical parameters of the saturated porous media are both taken into account. Numerical experiments are carried out for different cases in comparison with the standard finite element method. The numerical results show that the coupling multiscale finite element method can be successfully used for solving the complicated coupling problems. It reduces greatly the computing effort in both memory and time for transient problems.