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.     

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.     

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.     

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)相似文献   

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.     

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.     

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.     

Flow and transport in fractured tuff at Yucca Mountain: numerical experiments on fast preferential flow mechanisms

Recent discovery of bomb-related ³⁶Cl at depth in fractured tuff in the unsaturated zone at the Yucca Mountain candidate high-level waste (HLW) repository site has called into question the usual modeling assumptions based on the equivalent continuum model (ECM). A dual continuum model (DCM) for simulating transient flow and transport at Yucca Mountain is developed. In order to ensure properly converged flow solutions, which are used in the transport simulation, a new flow solution convergence criteria is derived. An extensive series of simulation studies is presented which indicates that rapid movement of solute through the fractures will not occur unless there are intense episodic infiltration events. Movement of solute in the environs of the repository is enhanced if the properties of the tuff layer at the repository horizon are modified from current best-estimate values. Due to a large advective–dispersive coupling between the matrix and fractures, the matrix acts as a major buffer which inhibits rapid transport along the fractures. Consequently, fast movement of solutes through the fractures to the repository depth can only be explained if the matrix–fracture coupling term is significantly reduced from a value that would be calculated on the basis of data currently available.     

Evidence of one-dimensional scale-dependent fractional advection-dispersion

A semi-analytical inverse method and the corresponding program FADEMain for parameter estimation of the fractional advection-dispersion equation (FADE) were developed in this paper. We have analyzed Huang et al.'s [Huang, K., Toride, N., van Genuchten, M.Th., 1995. Experimental investigation of solute transport in large homogeneous and heterogeneous saturated soil columns. Trans. Porous Media 18, 283-302.] laboratory experimental data of conservative solute transport in 12.5-m long homogeneous and heterogeneous soil columns to test the non-Fickian dispersion theory of FADE. The dispersion coefficient was calculated by fitting the analytical solution of FADE to the measured data at different transport scales. We found that the dispersion coefficient increased exponentially with transport scale for the homogeneous column, whereas it increased with transport scale in a power law function for the heterogeneous column. The scale effect of the dispersion coefficient in the heterogeneous soil was much more significant comparing to that in the homogeneous soil. The increasing rate of dispersion coefficient versus transport distance was smaller for FADE than that for the advection-dispersion equation (ADE). Finite difference numerical approximations of the scale-dependent FADE were established to interpret the experimental results. The numerical solutions were found to be adequate for predicting scale-dependent transport in the homogeneous column, while the prediction for the heterogeneous column was less satisfactory.     

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.     

Electrokinetic ion transport through unsaturated soil: 2. Application to a heterogeneous field site

Results of a field demonstration of electrokinetic transport of acetate through an unsaturated heterogeneous soil are compared to numerical modeling predictions. The numerical model was based on the groundwater flow and transport codes MODFLOW and MT3D modified to account for electrically induced ion transport. The 6-month field demonstration was conducted in an unsaturated layered soil profile where the soil moisture content ranged from 4% to 28% (m3 m(-3)). Specially designed ceramic-cased electrodes maintained a steady-state moisture content and electric potential field between the electrodes during the field demonstration. Acetate, a byproduct of acetic acid neutralization of the cathode electrolysis reaction, was transported from the cathode to the anode by electromigration. Field demonstration results indicated preferential transport of acetate through soil layers exhibiting higher moisture content/electrical conductivity. These field transport results agree with theoretical predictions that electromigration velocity is proportional to a power function of the effective moisture content. A numerical model using a homogeneous moisture content/electrical conductivity domain did not adequately predict the acetate field results. Numerical model predictions using a three-layer electrical conductivity/moisture content profile agreed qualitatively with the observed acetate distribution. These results suggest that field heterogeneities must be incorporated into electrokinetic models to predict ion transport at the field-scale.     

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.     

Effective rates of heavy metal release from alkaline wastes--quantified by column outflow experiments and inverse simulations

Column outflow experiments operated at steady state flow conditions do not allow the identification of rate limited release processes. This requires an alternative experimental methodology. In this study, the aim was to apply such a methodology in order to identify and quantify effective release rates of heavy metals from granular wastes. Column experiments were conducted with demolition waste and municipal waste incineration (MSWI) bottom ash using different flow velocities and multiple flow interruptions. The effluent was analyzed for heavy metals, DOC, electrical conductivity and pH. The breakthrough-curves were inversely modeled with a numerical code based on the advection–dispersion equation with first order mass-transfer and nonlinear interaction terms. Chromium, Copper, Nickel and Arsenic are usually released under non-equilibrium conditions. DOC might play a role as carrier for those trace metals. By inverse simulations, generally good model fits are derived. Although some parameters are correlated and some model deficiencies can be revealed, we are able to deduce physically reasonable release-mass-transfer time scales. Applying forward simulations, the parameter space with equifinal parameter sets was delineated. The results demonstrate that the presented experimental design is capable of identifying and quantifying non-equilibrium conditions. They show also that the possibility of rate limited release must not be neglected in release and transport studies involving inorganic contaminants.     

An application of the theory of kinematics of mixing to the study of tropospheric dispersion

The most common technique used for numerical simulations of tracer mixing is that of the numerical solution of the advection–diffusion equation with the unresolved fluxes parameterized using the similarity theory. Despite correct predictions of the overall directions of transport, models based on a numerical solution of the advection–diffusion equation lack sufficient accuracy to correctly reproduce the coupling of mixing with small scale processes which are sensitive to the microstructure of the tracer distribution. The objective of this paper is to revisit the basic formalism employed in numerical models used to investigate atmospheric tracers. The main mathematical method proposed here is the theory of kinematics of mixing which could be applied effectively for simulations of atmospheric transport processes. At the beginning of the paper, we introduce simple mathematical transformations in order to demonstrate how complex topological structures are created by mixing processes. These idealistic flow systems are essential to explain transport properties of much more complex three-dimensional geophysical flows. An example of the application of the kinematics of mixing to the analysis of tracer transport on a planetary scale is presented in the following sections. The complex filamentary structures simulated in the numerical experiment are evaluated using some commonly applied statistical measures in order to compare the results with the data published in the literature. The results of the experiment are also analysed with the help of simple conceptual models of fluid filaments. The microstructure of the tracer distribution introduced in the paper is essential to increase our understanding of atmospheric transport and to develop more realistic parameterizations of small-scale mixing. The presented results could also be used to improve calculations of the coupling between microphysical processes and tracer mixing.     

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.     

Single- and dual-porosity modelling of flow in reclaimed mine soil cores with embedded lignitic fragments

Lignitic mine soils represent a typical two-scale dual-porosity medium consisting of a technogenic mixture of overburden sediments that include lignitic components as dust and as porous fragments embedded within a mostly coarse-textured matrix. Flow and transport processes in such soils are not sufficiently understood to predict the course of soil reclamation or of mine drainage. The objective of this contribution is to identify the most appropriate conceptual model for describing small-scale heterogeneity effects on flow on the basis of the physical structure of the system. Multistep flow experiments on soil cores are analyzed using either mobile–immobile or mobile–mobile type 1D dual-porosity models, and a 3D numerical model that considers a local-scale distribution of fragments. Simulations are compared with time series' of upward infiltration and matric potential heads measured at two depths using miniature tensiometers. The 3D and the 1D dual-permeability models yielded comparable results as long as pressure heads are in local equilibrium; however, could describe either the upward infiltration or the matric potential curves but not both at the same time. The mobile–immobile type dual-porosity model failed to describe the data. A simultaneous match with pressure heads and upward infiltration data could only be obtained with the 1D dual-permeability model (i.e., mobile–mobile) by assuming an additional restriction of the inter-domain water transfer. These results indicate that for unsaturated flow conditions at higher matric potential heads (i.e., here >− 40 hPa), water in a restricted part of the fragment domain must be more mobile as compared to water in the sandy matrix domain. Closer inspections of the pore system and first neutron radiographic imaging support the hypothesis that a more continuous pore region exists at these pressure heads in the vicinity of the lignitic fragments possibly formed by fragment contacts and a lignitic dust interface-region between the two domains. The results suggest that the small-scale structure is too complex as to be represented by weighted contributions of individual components alone.     

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.     

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.