One-dimensional simulation of solute transfer in saturated–unsaturated porous media using the discontinuous finite elements method

A one-dimensional transport model for simulating water flow and solute transport in homogeneous–heterogeneous, saturated–unsaturated porous media is presented. The model is composed of a combination of accurate numerical algorithms for solving the nonlinear Richard's and advection–dispersion equations (ADE). The mixed form of Richard's equation is solved using a standard finite element method (FEM) with primary variable switching. The transport equation is solved using operator splitting, with the discontinuous finite element method (DFE) for discretization of the advective term. A slope limiting procedure for DFE avoids numerical instabilities but creates very limited numerical dispersion for high Peclet numbers. An implicit finite differences scheme (FD) is used for the dispersive term.The unsaturated flow and transport model (Wamos-T) is applied to a variety of rigorous problems including transient flow, heterogeneous medium and abrupt variations of velocity in magnitude and direction due to time-varying boundary conditions. It produces accurate and mass-conservative solutions for a very large range of grid Peclet numbers. The Wamos-T model is a good and robust alternative for the simulation of mass transport in unsaturated domain.     

One-dimensional simulation of solute transfer in saturated-unsaturated porous media using the discontinuous finite elements method

A one-dimensional transport model for simulating water flow and solute transport in homogeneous-heterogeneous, saturated-unsaturated porous media is presented. The model is composed of a combination of accurate numerical algorithms for solving the nonlinear Richard's and advection-dispersion equations (ADE). The mixed form of Richard's equation is solved using a standard finite element method (FEM) with primary variable switching. The transport equation is solved using operator splitting, with the discontinuous finite element method (DFE) for discretization of the advective term. A slope limiting procedure for DFE avoids numerical instabilities but creates very limited numerical dispersion for high Peclet numbers. An implicit finite differences scheme (FD) is used for the dispersive term. The unsaturated flow and transport model (Wamos-T) is applied to a variety of rigorous problems including transient flow, heterogeneous medium and abrupt variations of velocity in magnitude and direction due to time-varying boundary conditions. It produces accurate and mass-conservative solutions for a very large range of grid Peclet numbers. The Wamos-T model is a good and robust alternative for the simulation of mass transport in unsaturated domain.     

Modeling colloid-facilitated transport of multi-species contaminants in unsaturated porous media

Colloid-facilitated transport has been recognized as a potentially important and overlooked contaminant transport process. In particular, it has been observed that conventional two phase sorption models are often unable to explain transport of highly sorbing compounds in the subsurface appropriately in the presence of colloids. In this study a one-dimensional model for colloid-facilitated transport of chemicals in unsaturated porous media is developed. The model has parts for simulating coupled flow, and colloid transport and dissolved and colloidal contaminant transport. Richards' equation is solved to model unsaturated flow, and the effect of colloid entrapment and release on porosity and hydraulic conductivity of the porous media is incorporated into the model. Both random sequential adsorption and Langmuir approaches have been implemented in the model in order to incorporate the effect of surface jamming. The concept of entrapment of colloids into the air-water interface is used for taking into account the effect of retardation caused due to existence of the air phase. A non-equilibrium sorption approach with options of linear and Langmuir sorption assumptions are implemented that can represent the competition and site saturation effects on sorption of multiple compounds both to the solid matrix and to the colloidal particles. Several demonstration calculations are performed and the conditions in which the non-equilibrium model can be approximated by an equilibrium model are also studied.     

Analytical model of solute transport by unsteady unsaturated gravitational infiltration

Penetration of reactive solute into a soil during a cycle of water infiltration and redistribution is investigated by deriving analytical closed form solutions for fluid flux, moisture content and contaminant concentration. The solution is developed for gravitational flow and advective transport and is applied to two scenarios of solute applications encountered in the applications: a finite pulse of solute dissolved in irrigation water and an instantaneous pulse broadcasted onto the soil surface. Through comparison to simulations of Richards' flow, capillary suction is shown to have contrasting effects on the upper and lower boundaries of the fluid pulse, speeding penetration of the wetting front and reducing the rate of drying. This leads to agreement between the analytical and numerical solutions for typical field and experimental conditions. The analytical solution is further incorporated into a stochastic column model of flow and transport to compute mean solute concentration in a heterogeneous field. An unusual phenomenon of plume contraction is observed at long times of solute propagation during the drying stage. The mean concentration profiles match those of the Monte-Carlo simulations for capillary length scales typical of sandy soils.     

A novel two-dimensional model for colloid transport in physically and geochemically heterogeneous porous media

A two-dimensional model for colloid transport in geochemically and physically heterogeneous porous media is presented. The model considers patchwise geochemical heterogeneity, which is suitable to describe the chemical variability of many surficial aquifers with ferric oxyhydroxide-coated porous matrix, as well as spatial variability of hydraulic conductivity, which results in heterogeneous flow field. The model is comprised of a transient fluid flow equation, a transient colloid transport equation, and an equation for the dynamics of colloid deposition and release. Numerical simulations were carried out with the model to investigate the colloid transport behavior in layered and randomly heterogeneous porous media. Results demonstrate that physical and geochemical heterogeneities markedly affect the colloid transport behavior. Layered physical or geochemical heterogeneity can result in distinct preferential flow paths of colloidal particles. Furthermore, the combined effect of layered physical and geochemical heterogeneity may result in enhanced or reduced preferential flow of colloids. Random distribution of physical heterogeneity (hydraulic conductivity) results in a random flow field and an irregularly distributed colloid concentration profile in the porous medium. Contrary to random physical heterogeneity, the effect of random patchwise geochemical heterogeneity on colloid transport behavior is not significant. It is mostly the mean value of geochemical heterogeneity rather than its distribution that governs the colloid transport behavior.     

Characterization of flow and transport processes within the unsaturated zone of Yucca Mountain, Nevada, under current and future climates

This paper presents a large-scale modeling study characterizing fluid flow and tracer transport in the unsaturated zone of Yucca Mountain, Nevada, a potential repository site for storing high-level radioactive waste. The study has been conducted using a three-dimensional numerical model, which incorporates a wide variety of field data and takes into account the coupled processes of flow and transport in the highly heterogeneous, unsaturated fractured porous rock. The modeling approach is based on a dual-continuum formulation of coupled multiphase fluid and tracer transport through fractured porous rock. Various scenarios of current and future climate conditions and their effects on the unsaturated zone are evaluated to aid in the assessment of the proposed repository's system performance using different conceptual models. These models are calibrated against field-measured data. Model-predicted flow and transport processes under current and future climates are discussed.     

Semi-analytical solution of two-dimensional sharp interface LNAPL transport models

In this paper, we present semi-analytical solutions for two-dimensional equations governing transport of Light Non-Aqueous Phase Liquids (LNAPL) in unconfined aquifers. The proposed model is based on sharp interface displacement and steady groundwater flow assumptions, where both the water–LNAPL interface and the LNAPL–air interface are represented as sharp interfaces. In the case of steady groundwater flow, these equations can be reduced to a two-dimensional nonlinear solute transport equation, with the LNAPL thickness in the free product lens being the primary unknown variable. The linearized form of this solute transport equation falls into the category of two-dimensional transport equation with time-dependent dispersion coefficients. This equation can be solved analytically for an infinite domain region. In this paper, the general form of the analytical solution for the transport equation, as well as the solutions for some specific cases are presented. To demonstrate the utility of the proposed solution, numerical results obtained for two example problems are discussed and presented comparatively with a finite-element solution and other more restrictive solutions available in the literature. Although the solutions discussed in this paper have some simplifying assumptions, such as sharp-interfaces between fluid phases, steady groundwater flow and homogeneous aquifer properties, the semi-analytical solutions presented in this study may be used effectively as bench mark solutions in evaluating LNAPL migration in the subsurface. These solutions are simple and cost effective to implement and may be used in the calibration of other more complex numerical solutions that can be found in the literature.     

Transport of a non-volatile contaminant in unsaturated porous media in the presence of colloids

Stable colloidal particles can travel long distances in subsurface environments and carry particle-reactive contaminants with them to locations further than predicted by the conventional advective-dispersive transport equation. When such carriers exist in a saturated porous medium, the system can be idealized as consisting of three phases: an aqueous phase, a carrier phase, and a stationary solid matrix phase. However, when colloids are present in an unsaturated porous medium, the system representation should include one more phase, i.e. the air phase. In the work reported, a mathematical model was developed to describe the transport and fate of the colloidal particles and a non-volatile contaminant in unsaturated porous media. The model is based on mass balance equations in a four-phase porous medium. Colloid mass transfer mechanisms among aqueous, solid matrix, and air phases, and contaminant mass transfer between aqueous and colloid phases are represented by kinetic expressions. Governing equations are non-dimensionalized and solved to investigate colloid and contaminant transport in an unsaturated porous medium. A sensitivity analysis of the transport model was utilized to assess the effects of several parameters on model behavior. The colloid transport model matches successfully with experimental data of Wan and Wilson. The presence of air-water interface retards the colloid transport significantly counterbalancing the facilitating effect of colloids. However, the retardation of contaminant transport by colloids is highly dependent on the properties of the contaminant and the colloidal surface.     

Multicomponent simulation of wastewater-derived nitrogen and carbon in shallow unconfined aquifers. I. Model formulation and performance

One of the most common methods to dispose of domestic wastewater involves the release of septic effluent from drains located in the unsaturated zone. Nitrogen from such systems is currently of concern because of nitrate contamination of drinking water supplies and eutrophication of coastal waters. The objectives of this study are to develop and assess the performance of a mechanistic flow and reactive transport model which couples the most relevant physical, geochemical and biochemical processes involved in wastewater plume evolution in sandy aquifers. The numerical model solves for variably saturated groundwater flow and reactive transport of multiple carbon- and nitrogen-containing species in a three-dimensional porous medium. The reactive transport equations are solved using the Strang splitting method which is shown to be accurate for Monod and first- and second-order kinetic reactions, and two to four times more efficient than sequential iterative splitting. The reaction system is formulated as a fully kinetic chemistry problem, which allows for the use of several special-purpose ordinary differential equation (ODE) solvers. For reaction systems containing both fast and slow kinetic reactions, such as the combined nitrogen-carbon system, it is found that a specialized stiff explicit solver fails to obtain a solution. An implicit solver is more robust and its computational performance is improved by scaling of the fastest reaction rates. The model is used to simulate wastewater migration in a 1-m-long unsaturated column and the results show significant oxidation of dissolved organic carbon (DOC), the generation of nitrate by nitrification, and a slight decrease in pH.     

On the selection of primary variables in numerical formulation for modeling multiphase flow in porous media

Selecting the proper primary variables is a critical step in efficiently modeling the highly nonlinear problem of multiphase subsurface flow in a heterogeneous porous-fractured media. Current simulation and ground modeling techniques consist of (1) spatial discretization of mass and/or heat conservation equations using finite difference or finite element methods; (2) fully implicit time discretization; (3) solving the nonlinear, discrete algebraic equations using a Newton iterative scheme. Previous modeling efforts indicate that the choice of primary variables for a Newton iteration not only impacts computational performance of a numerical code, but may also determine the feasibility of a numerical modeling study in many field applications. This paper presents an analysis and general recommendations for selecting primary variables in simulating multiphase, subsurface flow for one-active phase (Richards' equation), two-phase (gas and liquid) and three-phase (gas, water and nonaqueous phase liquid or NAPL) conditions. In many cases, a dynamic variable switching or variable substitution scheme may have to be used in order to achieve optimal numerical performance and robustness. The selection of primary variables depends in general on the sensitivity of the system of equations to the variables selected at given phase and flow conditions. We will present a series of numerical tests and large-scale field simulation examples, including modeling one (active)-phase, two-phase and three-phase flow problems in multi-dimensional, porous-fractured subsurface systems.     

Numerical methods for describing chemical transport in the unsaturated zone of the subsurface

Finite-difference and finite-element methods of approximation have been extended to solve the one-dimensional nonlinear partial differential equations that describe the simultaneous transport of heat, moisture and chemical in the unsaturated zone. Especially for chemical transport, nodal spacing criteria are required to minimize numerical dispersion and oscillatory behavior in the solution vector for chemical concentration. Conservative criteria for nodal spacing for saturated flow can be used to set nodal spacing for unsaturated zone transport. When nodal spacing criteria are satisfied, for the same set of transport and boundary conditions, chemical concentration profiles calculated by the two numerical methods will be almost the same. A situation that is simulated very well with one-dimensional models, is the application of chemicals to land surfaces. To compare and contrast the characteristics of solutions given by the two numerical methods, moisture content, temperature and chemical concentration profiles for a 75-day period after application in the unsaturated zone are calculated for two representative types of organic chemicals. In the first, the chemical is very slowly degraded in the subsurface environment but strongly sorbed to soil surfaces. In the second, the chemical is rapidly degraded but weakly sorbed to soil surfaces. Because of differences in sorption coefficients and mechanisms of degradation, for the same set of hydrodynamic properties of the subsurface, the weakly sorbed chemical is more widely distributed throughout the unsaturated zone, whereas the strongly sorbed chemical stays very close to where it is put initially with little penetration into the subsurface. Satisfying nodal spacing criteria minimizes the impact of the method of approximation on the calculated solutions of the transport equations. For better model predictive performance, however, there are needs for more fundamental information on processes governing transport in the subsurface.     

Development and analysis of an adaptive transport scheme

An approach to solving the advection dominated atmospheric mass transport problem which adaptively constructs and updates the computational mesh is implemented and analyzed. The formulation of the mesh adaptation algorithm allows for information from other model processes as well as transport to influence mesh refinement. Comparisons to other methods are limited to a pure advection test problem. The scheme is based on a Petrov–Galerkin finite element method using quadratic interpolating polynomials over triangular elements. The temporal component of the transport equation is discretized using the Crank–Nicolson method. A single time step is applied to the entire domain regardless of mesh refinement levels. Two means of upwind biasing the solution via modification of the finite element weighting functions are investigated. Criteria for refining and unrefining the mesh based on advective flux and an estimate of the error are proposed and compared. The effects of varying the number of time steps taken between mesh refinements are analyzed. Results for adapted mesh solutions of a test problem are shown to be superior to corresponding uniform mesh finite element solutions. Comparisons to several other popular uniform mesh advection schemes indicate either equal or improved accuracy for the adapted mesh solutions.     

Effects of pore volume-transmissivity correlation on transport phenomena

The relevant velocity that describes transport phenomena in a porous medium is the pore velocity. For this reason, one needs not only to describe the variability of transmissivity, which fully determines the Darcy velocity field for given source terms and boundary conditions, but also any variability of the pore volume. We demonstrate that hydraulically equivalent media with exactly the same transmissivity field can produce dramatic differences in the displacement of a solute if they have different pore volume distributions. In particular, we demonstrate that correlation between pore volume and transmissivity leads to a much smoother and more homogeneous solute distribution. This was observed in a laboratory experiment performed in artificial fractures made of two plexiglass plates into which a space-dependent aperture distribution was milled. Using visualization by a light transmission technique, we observe that the solute behaviour is much smoother and more regular after the fractures are filled with glass powder, which plays the role of a homogeneous fault gouge material. This is due to a perfect correlation between pore volume and transmissivity that causes pore velocity to be not directly dependent on the transmissivity, but only indirectly through the hydraulic gradient, which is a much smoother function due to the diffusive behaviour of the flow equation acting as a filter. This smoothing property of the pore volume-transmissivity correlation is also supported by numerical simulations of tracer tests in a dipole flow field. Three different conceptual models are used: an empty fracture, a rough-walled fracture filled with a homogeneous material and a parallel-plate fracture with a heterogeneous fault gouge. All three models are hydraulically equivalent, yet they have a different pore volume distribution. Even if piezometric heads and specific flow rates are exactly the same at any point of the domain, the transport process differs dramatically. These differences make it important to discriminate in situ among different conceptual models in order to simulate correctly the transport phenomena. For this reason, we study the solute breakthrough and recovery curves at the extraction wells. Our numerical case studies show that discrimination on the basis of such data might be impossible except under very favourable conditions, i.e. the integral scale of the transmissivity field has to be known and small compared to the dipole size. If the latter conditions are satisfied, discrimination between the rough-walled fracture filled with a homogeneous material and the other two models becomes possible, whereas the parallel-plate fracture with a heterogeneous fault gouge and the empty fracture still show identifiability problems. The latter may be solved by inspection of aperture and pressure testing.     

Transport of a decay chain in homogenous porous media: analytical solutions

With the aid of integral transforms, analytical solutions for the transport of a decay chain in homogenous porous media are derived. Unidirectional steady-state flow and radial steady-state flow in single and multiple porosity media are considered. At least in Laplace domain, all solutions can be written in closed analytical formulae. Partly, the solutions can also be inverted analytically. If not, analytical calculation of the steady-state concentration distributions, evaluation of temporal moments and numerical inversion are still possible. Formulae for several simple boundary conditions are given and visualized in this paper. The derived novel solutions are widely applicable and are very useful for the validation of numerical transport codes.     

Infiltration and redistribution of LNAPL into unsaturated layered porous media

Enhanced understanding of light non-aqueous phase liquid (LNAPL) infiltration into heterogeneous porous media is important for the effective design of remediation strategies. We used a 2-D experimental facility that allows for visual observation of LNAPL contours in order to study LNAPL redistribution in a layered porous medium. The layers are situated in the unsaturated zone near the watertable and they are inclined to be able to observe the effect of discontinuities in capillary forces and relative permeabilities. Two experiments were performed. The first experiment consisted of LNAPL infiltration into a fine sand matrix with a coarse sand layer, and the second experiment consisted of a coarse sand matrix and a fine sand layer. The numerical multi-phase flow model STOMP was validated with regard to the experimental results. This model is able to adequately reproduce the experimental LNAPL contours. Numerical sensitivity analysis was also performed. The capillarity contrast between sands was found to be the main controlling factor determining the final LNAPL distribution.     

Using the multiple cell balance method to solve the problem of two-dimensional groundwater flow and contaminant transport with nonequilibrium sorption

The problem of large-scale contamination of groundwater by relatively low levels of organic contaminants is most frequently addressed by extracting and treating the impacted groundwater. This pump-and-treat strategy is often unsuccessful because of difficulties encountered in recovering the contaminants from relatively immobile zones within the porous medium. These zones can exist at the particle scale, as intraparticle or intra-aggregate porosity, and at the larger scales, as low-permeability layers or lenses interspersed in substantially more permeable layers. This work focuses on achieving an efficient numerical solution to a system of groundwater flow and contaminant transport equations that sufficiently captures the dynamics of slow desorption in a two-dimensional porous medium. The conceptual model and governing equations are presented. A numerical method for solving the governing equations, the upstream-weighted, multiple cell balance (UMCB) method, is proposed. The UMCB algorithm has been employed previously for the case of solute transport with equilibrium sorption, and is extended here to the nonequilibrium case. The approach employs a finite-element basis function and a finite-difference local mass balance, and is designed to reduce computational and storage requirements, while minimizing the mass balance error. The computational grid is formed by division of the flow domain into triangular elements. An invented node at the center of each element divides the element into three subtriangular regions. By linking the center of each triangular element and the mid-point of each elemental side, a multiangular region, referred to as an exclusive subdomain, is defined. The discretized system of governing equations is derived from the integral form that describes the mass balance in the exclusive subdomain of each node. The paper details the application of the numerical method, and demonstrates that the method is reasonably accurate and computationally efficient for a two-dimensional domain subject to nonequilibrium sorption.     

Stochastic-convective transport with nonlinear reactions and mixing: finite streamtube ensemble formulation for multicomponent reaction systems with intra-streamtube dispersion

An effective streamtube ensemble method is developed to upscale convective-dispersive transport with multicomponent nonlinear reactions in steady nonuniform flow. The transport is cast in terms of a finite ensemble of independent discrete streamtubes that approximate convective transport along macroscopically averaged pathlines and dispersive transport longitudinally as microscopic mixing within streamtubes. The representation of fate and transport via a finite ensemble of effective linear streamtubes, allows the treatment of arbitrarily complex reaction systems involving both homogeneous and heterogeneous reactions, and longitudinal dispersive/diffusive mixing within streamtubes. This allows the use of reactive-transport codes designed to solve such problems in an Eulerian framework, as opposed to reliance on closed-form (convolutional or canonical) expressions for reactive transport in exclusively convective streamtubes. The approach requires both reactive-transport solutions for a representative ensemble of one-dimensional convective-dispersive-reactive streamtubes and the distribution of flux over the streamtube ensemble variants, and it does not allow for lateral mixing between streamtubes. Here, the only ensemble variant is travel time. The discussion details the way that the conventional Eulerian fate and transport model is converted first into an ensemble of transports along three-dimensional streamtubes of unknown geometry, and then to approximate one-dimensional streamtubes that are designed to honor the important global properties of the transport. Conditions under which such an 'equivalent' ensemble of one-dimensional streamtubes are described. The breakthrough curve of a nonreactive tracer in the ensemble is expressed as a combined Volterra-Fredholm integral equation, which serves as the basis for estimation of the distribution of flux over the variant of the ensemble, travel time. Transient convective speed and the effects of errors in flux distributions are described, and the method is applied to a demonstration problem involving nonlinear multicomponent reaction kinetics and strongly nonuniform flow.     

Variable-density flow in heterogeneous porous media — Laboratory experiments and numerical simulations

Konz, M., Ackerer, P., Younes, A., Huggenberger, P., Zechner, E., 2009a. 2D Stable Layered Laboratory-scale Experiments for Testing Density-coupled Flow Models. Water Resources Research, 45. doi:10.1029/2008WR007118., a series of laboratory-scale 2D tank experiments were conducted and accurately simulated for density driven flow problems on homogeneous porous media. In the present work, we extended the numerical and experimental studies to heterogeneous problems. The heterogeneous porous medium was constructed with a low permeability zone in the centre of the tank and had well-defined parameters and boundary conditions. Concentration distributions were measured in high resolution using a photometric method and an image analysis technique. The numerical model used for the simulations was based on efficient advanced approximations for both spatial and temporal discretizations. The Method Of Lines (MOL) was used to allow higher-order temporal discretization. Three different boundary conditions, corresponding to different localizations of the inflow and the outflow openings at the opposite edges of the tank, were applied to investigate different flow scenarios in the heterogeneous porous medium flow tank. Simulation results of all three density coupled experiments revealed a density-dependent behavior of dispersion. Thus, a reduction of dispersivites was required to obtain a good matching of the experimental data. The high quality of the experiments enabled a detailed testing of numerical variable-density flow codes under heterogeneous conditions. Therefore, the experiments were considered to be reliable benchmark tests.