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.     

A numerical model for transient-hysteretic flow and solute transport in unsaturated porous media

A two-dimensional flow and transport model was developed for simulating transient water flow and nonreactive solute transport in heterogeneous, unsaturated porous media containing air and water. The model is composed of a unique combination of robust and accurate numerical algorithms for solving the Richards', Darcy flux, and advection-dispersion equations. The mixed form of Richards' equation is solved using a finite-element formulation and a modified Picard iteration scheme. Mass lumping is employed to improve solution convergence and stability behavior. The flow algorithm accounts for hysteresis in the pressure head-water content relationship. Darcy fluxes are approximated with a Galerkin and Petrov-Galerkin finite-element method developed for random heterogeneous porous media. The transport equation is solved using an Eulerian-Lagrangian method. A multi-step, fourth-order Runge-Kutta, reverse particle tracking technique and a quadratic-linear interpolation scheme are shown to be superior for determining the advective concentration. A Galerkin finite-element method is used for approximating the dispersive flux. The unsaturated flow and transport model was applied to a variety of rigorous problems and was found to produce accurate, mass-conserving solutions when compared to analytical solutions and published numerical results.     

Truncation errors in finite difference models for solute transport equation with first-order reaction

The truncation errors associated with finite difference solutions of the advection-dispersion equation with first-order reaction are formulated from a Taylor analysis. The error expressions are based on a general form of the corresponding difference equation and a temporally and spatially weighted parametric approach is used for differentiating among the various finite difference schemes. The numerical truncation errors are defined using Peclet and Courant numbers and a new Sink/Source dimensionless number. It is shown that all of the finite difference schemes suffer from truncation errors. In particular it is shown that the Crank–Nicolson approximation scheme does not have second order accuracy for this case. The effects of these truncation errors on the solution of an advection–dispersion equation with a first order reaction term are demonstrated by comparison with an analytical solution. The results show that these errors are not negligible and that correcting the finite difference scheme for them results in a more accurate solution.     

Coupling geochemistry with a particle tracking transport model

Coupling geochemistry and transport appears unavoidable since it is rare that either of these two phenomena alone can account for the movement of solutes in groundwater. The chemical model is based on thermodynamic equilibrium. The method used is a Gibbs free energy minimization constrained by mass balances. The model calculates the aqueous speciation, the precipitation and the dissolution of pure minerals or solid solutions. The transport equation is solved by the random walk technique which avoids the problem of numerical dispersion for transport, but may be more time consuming than finite differences or elements if a large number of particles are necessary in order to get a sufficiently “smooth” solution. However, when the chemistry deals with a realistic number of elements (e.g., > 10), the cost of the chemistry computation largely dominates that of transport. Special techniques had to be developed in order to solve problems linked to the conditions present in some of the CEC CHEMVAL tests (boundary with fixed concentrations and very low Péclet numbers). The coupling consists of calculating the exchanges of chemical elements between two populations. The first population is sedentary, constituted by a mesh of fixed cells representing the composition of the solid phase. The other population is nomadic, represented by a set of particles which are advected by groundwater flow. A vector of real numbers is associated with each mobile particle. This vector accounts for the mass of each element dissolved in the moving liquid phase. For this reason, the transport equation is only solved once for the whole set of elements. The main assumptions that were necessary to perform the coupling in a simple way are discussed. Two applications are presented: (1) a verification compared to an analytical solution; and (2) the simulation of a percolation experiment through a sandstone core.     

Pore-scale modeling of dispersion in disordered porous media

We employ a direct pore-level model of incompressible flow that uses the modified moving particle semi-implicit (MMPS) method. The model is capable of simulating both unsteady- and steady-state flow directly in microtomography images of naturally-occurring porous media. We further develop this model to simulate solute transport in disordered porous media. The governing equations of flow and transport at the pore level, i.e., Navier-Stokes and convection-diffusion, are solved directly in the pore space mapped by microtomography techniques. Three naturally-occurring sandstones are studied in this work. We verify the accuracy of the model by comparing the computed longitudinal dispersion coefficients against the experimental data for a wide range of Peclet numbers, i.e., 5×10(-2)相似文献   

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.     

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 modelling of stable isotope fractionation of volatile organic compounds in the unsaturated zone

Analytical models were developed that simulate stable isotope ratios of volatile organic compounds (VOCs) near a point source contamination in the unsaturated zone. The models describe diffusive transport of VOCs, biodegradation and source ageing. The mass transport is governed by Fick's law for diffusion. The equation for reactive transport of VOCs in the soil gas phase was solved for different source geometries and for different boundary conditions. Model results were compared to experimental data from a one-dimensional laboratory column and a radial-symmetric field experiment. The comparison yielded a satisfying agreement. The model results clearly illustrate the significant isotope fractionation by gas phase diffusion under transient state conditions. This leads to an initial depletion of heavy isotopes with increasing distance from the source. The isotope evolution of the source is governed by the combined effects of isotope fractionation due to vaporisation, diffusion and biodegradation. The net effect can lead to an enrichment or depletion of the heavy isotope in the remaining organic phase, depending on the compound and element considered. Finally, the isotope evolution of molecules migrating away from the source and undergoing degradation is governed by a combined degradation and diffusion isotope effect. This suggests that, in the unsaturated zone, the interpretation of biodegradation of VOC based on isotopic data must always be based on a model combining gas phase diffusion and degradation.     

Modelling the migration of contaminants through variably saturated dual-porosity, dual-permeability chalk

In the Hesbaye region in Belgium, tracer tests performed in variably saturated fissured chalk rocks presented very contrasting results in terms of transit times, according to artificially controlled water recharge conditions prevailing during the experiments. Under intense recharge conditions, tracers migrated across the partially or fully saturated fissure network, at high velocity in accordance with the high hydraulic conductivity and low effective porosity (fracture porosity). At the same time, a portion of the tracer was temporarily retarded in the almost immobile water located in the matrix. Under natural infiltration conditions, the fissure network remained inactive. Tracers migrated downward through the matrix, at low velocity in relation with the low hydraulic conductivity and the large porosity of the matrix. Based on these observations, Brouyère et al. (2004a) [Brouyère, S., Dassargues, A., Hallet, V., 2004a. Migration of contaminants through the unsaturated zone overlying the Hesbaye chalky aquifer in Belgium: a field investigation, J. Contam. Hydrol., 72 (1-4), 135-164, doi: 10.1016/j.conhyd.2003.10.009] proposed a conceptual model in order to explain the migration of solutes in variably saturated, dual-porosity, dual-permeability chalk. Here, mathematical and numerical modelling of tracer and contaminant migration in variably saturated fissured chalk is presented, considering the aforementioned conceptual model. A new mathematical formulation is proposed to represent the unsaturated properties of the fissured chalk in a more dynamic and appropriate way. At the same time, the rock water content is partitioned between mobile and immobile water phases, as a function of the water saturation of the chalk rock. The groundwater flow and contaminant transport in the variably saturated chalk is solved using the control volume finite element method. Modelling the field tracer experiments performed in the variably saturated chalk shows the adequacy and usefulness of the new conceptual, mathematical and numerical model.     

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.     

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.     

Sensitivity of the solution of the Elder problem to density, velocity and numerical perturbations

In this paper the Elder problem is studied with the purpose of evaluating the inherent instabilities associated with the numerical solution of this problem. Our focus is first on the question of the existence of a unique numerical solution for this problem, and second on the grid density and fluid density requirements necessary for a unique numerical solution. In particular we have investigated the instability issues associated with the numerical solution of the Elder problem from the following perspectives: (i) physical instability issues associated with density differences; (ii) sensitivity of the numerical solution to idealization irregularities; and, (iii) the importance of a precise velocity field calculation and the association of this process with the grid density levels that is necessary to solve the Elder problem accurately. In the study discussed here we have used a finite element Galerkin model we have developed for solving density-dependent flow and transport problems, which will be identified as TechFlow. In our study, the numerical results of Frolkovic and de Schepper [Frolkovic, P. and H. de Schepper, 2001. Numerical modeling of convection dominated transport coupled with density-driven flow in porous media, Adv. Water Resour., 24, 63-72.] were replicated using the grid density employed in their work. We were also successful in duplicating the same result with a less dense grid but with more computational effort based on a global velocity estimation process we have adopted. Our results indicate that the global velocity estimation approach recommended by Yeh [Yeh, G.-T., 1981. On the computation of Darcian velocity and mass balance in finite element modelling of groundwater flow, Water Resour. Res., 17(5), 1529-1534.] allows the use of less dense grids while obtaining the same accuracy that can be achieved with denser grids. We have also observed that the regularity of the elements in the discretization of the solution domain does make a difference in obtaining a unique stationary solution for this problem. The results of our study also indicate that the density differences are critical in the solution of the Elder problem and that high density differences lead to the physical instability that is inherent with this problem. Other than the physical instability associated with the level of density differences used in the Elder problem, the following two points should be considered in solving the Elder problem in a consistent manner: (i) strict attention should be paid to the vertical grid Peclet number in developing the criteria for convergent grid selection; and, (ii) with a globally continuous velocity calculation stable solutions can be obtained at lower grid densities.     

Modeling radium and radon transport through soil and vegetation

A one-dimensional flow and transport model was developed to describe the movement of two fluid phases, gas and water, within a porous medium and the transport of 226Ra and 222Rn within and between these two phases. Included in this model is the vegetative uptake of water and aqueous 226Ra and 222Rn that can be extracted from the soil via the transpiration stream. The mathematical model is formulated through a set of phase balance equations and a set of species balance equations. Mass exchange, sink terms and the dependence of physical properties upon phase composition couple the two sets of equations. Numerical solution of each set, with iteration between the sets, is carried out leading to a set-iterative compositional model. The Petrov-Galerkin finite element approach is used to allow for upstream weighting if required for a given simulation. Mass lumping improves solution convergence and stability behavior. The resulting numerical model was applied to four problems and was found to produce accurate, mass conservative solutions when compared to published experimental and numerical results and theoretical column experiments. Preliminary results suggest that the model can be used as an investigative tool to determine the feasibility of phytoremediating radium and radon-contaminated soil.     

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.     

A triple-continuum approach for modeling flow and transport processes in fractured rock

This paper presents a triple-continuum conceptual model for simulating flow and transport processes in fractured rock. Field data collected from the unsaturated zone of Yucca Mountain, a repository site of high-level nuclear waste, show a large number of small-scale fractures. The effect of these small fractures has not been considered in previous modeling investigations within the context of a continuum approach. A new triple-continuum model (consisting of matrix, small-fracture, and large-fracture continua) has been developed to investigate the effect of these small fractures. This paper derives the model formulation and discusses the basic triple-continuum behavior of flow and transport processes under different conditions, using both analytical solutions and numerical approaches. The simulation results from the site-scale model of the unsaturated zone of Yucca Mountain indicate that these small fractures may have an important effect on radionuclide transport within the mountain.     

Operator-splitting errors in coupled reactive transport codes for transient variably saturated flow and contaminant transport in layered soil profiles

One possible way of integrating subsurface flow and transport processes with (bio)geochemical reactions is to couple by means of an operator-splitting approach two completely separate codes, one for variably-saturated flow and solute transport and one for equilibrium and kinetic biogeochemical reactions. This paper evaluates the accuracy of the operator-splitting approach for multicomponent systems for typical soil environmental problems involving transient atmospheric boundary conditions (precipitation, evapotranspiration) and layered soil profiles. The recently developed HP1 code was used to solve the coupled transport and chemical equations. For steady-state flow conditions, the accuracy was found to be mainly a function of the adopted spatial discretization and to a lesser extent of the temporal discretization. For transient flow situations, the accuracy depended in a complex manner on grid discretization, time stepping and the main flow conditions (infiltration versus evaporation). Whereas a finer grid size reduced the numerical errors during steady-state flow or the main infiltration periods, the errors sometimes slightly increased (generally less than 50%) when a finer grid size was used during periods with a high evapotranspiration demand (leading to high pressure head gradients near the soil surface). This indicates that operator-splitting errors are most significant during periods with high evaporative boundary conditions. The operator-splitting errors could be decreased by constraining the time step using the performance index (the product of the grid Peclet and Courant numbers) during infiltration, or the maximum time step during evapotranspiration. Several test problems were used to provide guidance for optimal spatial and temporal discretization.     

Field investigation into unsaturated flow and transport in a fault: model analyses

Results of a fault test performed in the unsaturated zone of Yucca Mountain, Nevada, were analyzed using a three-dimensional numerical model. The fault was explicitly represented as a discrete feature and the surrounding rock was treated as a dual-continuum (fracture-matrix) system. Model calibration against seepage and water-travel-velocity data suggests that lithophysal cavities connected to fractures can considerably enhance the effective fracture porosity and therefore retard water flow in fractures. Comparisons between simulation results and tracer concentration data also indicate that matrix diffusion is an important mechanism for solute transport in unsaturated fractured rock. We found that an increased fault-matrix and fracture-matrix interface areas were needed to match the observed tracer data, which is consistent with previous studies. The study results suggest that the current site-scale model for the unsaturated zone of Yucca Mountain may underestimate radionuclide transport time within the unsaturated zone, because an increased fracture-matrix interface area and the increased effective fracture porosity arising from lithophysal cavities are not considered in the current site-scale model.     

Analytical solution of two-dimensional solute transport in an aquifer–aquitard system

This study deals with two-dimensional solute transport in an aquifer–aquitard system by maintaining rigorous mass conservation at the aquifer–aquitard interface. Advection, longitudinal dispersion, and transverse vertical dispersion are considered in the aquifer. Vertical advection and diffusion are considered in the aquitards. The first-type and the third-type boundary conditions are considered in the aquifer. This study differs from the commonly used averaged approximation (AA) method that treats the mass flux between the aquifer and aquitard as an averaged volumetric source/sink term in the governing equation of transport in the aquifer. Analytical solutions of concentrations in the aquitards and aquifer and mass transported between the aquifer and upper or lower aquitard are obtained in the Laplace domain, and are subsequently inverted numerically to yield results in the real time domain (the Zhan method). The breakthrough curves (BTCs) and distribution profiles in the aquifer obtained in this study are drastically different from those obtained using the AA method. Comparison of the numerical simulation using the model MT3DMS and the Zhan method indicates that the numerical result differs from that of the Zhan method for an asymmetric case when aquitard advections are at the same direction. The AA method overestimates the mass transported into the upper aquitard when an upward advection exists in the upper aquitard. The mass transported between the aquifer and the aquitard is sensitive to the aquitard Peclet number, but less sensitive to the aquitard diffusion coefficient.