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.     

Operator-splitting procedures for reactive transport and comparison of mass balance errors

Operator-splitting (OS) techniques are very attractive for numerical modelling of reactive transport, but they induce some errors. Considering reactive mass transport with reversible and irreversible reactions governed by a first-order rate law, we develop analytical solutions of the mass balance for the following operator-splitting schemes: standard sequential non-iterative (SNI), Strang-splitting SNI, standard sequential iterative (SI), extrapolating SI, and symmetric SI approaches. From these analytical solutions, the operator-splitting methods are compared with respect to mass balance errors and convergence rates independently of the techniques used for solving each operator. Dimensionless times, NOS, are defined. They control mass balance errors and convergence rates. The following order in terms of decreasing efficiency is proposed: symmetric SI, Strang-splitting SNI, standard SNI, extrapolating SI and standard SI schemes. The symmetric SI scheme does not induce any operator-splitting errors, the Strang-splitting SNI appears to be O(N2OS) accurate, and the other schemes are first-order accurate.     

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.     

On efficient numerical solution of one-dimensional convection–diffusion equations in modelling atmospheric processes

A new approach is proposed to the numerical solution of one-dimensional convection–diffusion equations that arise in modelling atmospheric processes and air pollution modelling. The technique is based on upstream-type difference approximations for first-order derivatives and non-standard difference approximations for second-order derivatives of convection–diffusion equations. This approach leads to the significant qualitative improvements in the numerical solutions behaviour. The relative contribution of convection and diffusion is directly incorporated into the corresponding numerical scheme in such a way that large spatial grids can be taken without affecting solution stability. The method is compared with the contemporary computational schemes for solving problems with severe internal and boundary gradients and is shown to be stable and computationally efficient. The results of a numerical experiment are given.     

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.     

Effect of density-driven flow on the through-diffusion experiment

Diffusion is one of the main mechanisms of solute transport in pore water of geologic media. The effective diffusion coefficient of a solute in a rock is usually measured by the through-diffusion experiment. However, in this experiment, the effect of advection, induced by density difference between dense aqueous solution and pure water, has not been considered. To evaluate the effect of density-driven flow, a through-diffusion experiment using Fontainebleau sandstone was conducted for KCl and KI aqueous solutions with various densities. The measured effective diffusion coefficients were positively correlated with the density difference; the effective diffusion coefficient of a 1 M KI solution (density difference, 0.119 g/cm³) was one order of magnitude larger than that of a 0.1 M KCl solution (density difference, 0.005 g/cm³). The result is explained by a theoretical model using a diffusion–advection equation including Darcy's law. Based on the theory, a diagram to evaluate the condition at which the measured effective diffusion coefficient does not include the effect of advection is presented.     

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.     

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.     

Contaminant transport in groundwater in the presence of colloids and bacteria: Model development and verification

Colloids and bacteria (microorganisms) naturally exist in groundwater aquifers and can significantly impact contaminant migration rates. A conceptual model is first developed to account for the different physiochemical and biological processes, reaction kinetics, and different transport mechanisms of the combined system (contaminant–colloids–bacteria). All three constituents are assumed to be reactive with the reactions taking place between each constituent and the porous medium and also among the different constituents. A general linear kinetic reaction model is assumed for all reactive processes considered. The mathematical model is represented by fourteen coupled partial differential equations describing mass balance and reaction processes. Two of these equations describe colloid movement and reactions with the porous medium, four equations describe bacterial movement and reactions with colloids and the porous medium, and the remaining eight equations describe contaminant movement and its reactions with bacteria, colloids, and the porous medium. The mass balance equations are numerically solved for two-dimensional groundwater systems using a third-order, total variance-diminishing scheme (TVD) for the advection terms. Due to the complex coupling of the equations, they are solved iteratively each time step until a convergence criterion is met. The model is tested against experimental data and the results are favorable.     

First-passage-time transfer functions for groundwater tracer tests conducted in radially convergent flow

Forced-gradient groundwater tracer tests may be conducted using a variety of hydraulic schemes, so it is useful to have simple semi-analytic models available that can examine various injection/withdrawal scenarios. Models for radially convergent tracer tests are formulated here as transfer functions, which allow complex tracer test designs to be simulated by a series of simple mathematical expressions. These mathematical expressions are given in Laplace space, so that transfer functions may be placed in series by simple multiplication. Predicted breakthrough is found by numerically inverting the composite transfer function to the time-domain, using traditional computer programs or commercial mathematical software. Transport is assumed to be dictated by a radially convergent or uniform flow field, and is based upon an exact first-passage-time solution of the backward Fokker–Planck equation. These methods are demonstrated by simulating a weak-dipole tracer test conducted in a fractured granite formation, where mixing in the injection borehole is non-ideal.     

Sensitivity of mesoscale simulations of land–sea breeze to boundary layer turbulence parameterization

The sensitivity of mesoscale simulations of land and sea breeze circulation on the south east coast in the Chennai region of India to boundary layer turbulence parameterizations is studied using the community based PSU/NCAR mesoscale model MM5. High-resolution simulations are carried out with nested domains. Four widely used boundary layer turbulence parameterization schemes are selected for the study. Two of these schemes (Blackadar (BK) and medium range forecast (MRF)) are simple first-order non-local schemes and the other two (Mellor–Yamada (MY) Eta planetary boundary layer (PBL) and Gayno–Seaman (GS)) are more complex 1.5 order local schemes that include a prognostic equation for turbulence kinetic energy. The other physics used in the model are identical. The micro-meteorological tower and flux observations, GP sonde and radiosonde observations from the study region are used to compare the simulated mean variables. In spite of differences in complexity, all the schemes show similar capability in simulating the boundary layer temperature, humidity and winds. The land–sea breeze circulation and internal boundary layer formation, which are special phenomena at the coastal site, could be simulated by all the schemes. The BK, MRF schemes produced boundary layers that are more mixed than those produced with the MY and GS schemes. While all the schemes underestimated surface sensible heat fluxes during stable night conditions they are reasonably simulated during daytime by MRF and BK schemes. Among the tested schemes, the BK scheme has reasonably produced the PBL height, humidity, wind fields and proves suitable for operational forecasting. The results suggest that there is little improvement in overall accuracy of predictions with more complex turbulence parameterizations.     

Intercomparison of reactive transport models applied to UO2 oxidative dissolution and uranium migration

Oxidative dissolution of uranium dioxide (UO₂) and the subsequent migration of uranium in a subsurface environment and an underground waste disposal have been simulated with reactive transport models. In these systems, hydrogeological and chemical processes are closely entangled and their interdependency has been analyzed in detail, notably with respect to redox reactions, kinetics of mineralogical evolution and hydrodynamic migration of species of interest.Different codes, where among CASTEM, CHEMTRAP and HYTEC, have been used as an intercomparison and verification exercise. Although the agreement between codes is satisfactory, it is shown that the discretization method of the transport equation (i.e. finite elements (FE) versus mixed-hybrid FE and finite differences) and the sequential coupling scheme may lead to systematic discrepancies.     

Obtaining the porewater composition of a clay rock by modeling the in- and out-diffusion of anions and cations from an in-situ experiment

A borehole in the Callovo–Oxfordian clay rock in ANDRA's underground research facility was sampled during 1 year and chemically analyzed. Diffusion between porewater and the borehole solution resulted in concentration changes which were modeled with PHREEQC's multicomponent diffusion module. In the model, the clay rock's pore space is divided in free porewater (electrically neutral) and diffuse double layer water (devoid of anions). Diffusion is calculated separately for the two domains, and individually for all the solute species while a zero-charge flux is maintained. We explain how the finite difference formulas for radial diffusion can be translated into mixing factors for solutions. Operator splitting is used to calculate advective flow and chemical reactions such as ion exchange and calcite dissolution and precipitation. The ion exchange reaction is formulated in the form of surface complexation, which allows distributing charge over the fixed sites and the diffuse double layer. The charge distribution affects pH when calcite dissolves, and modeling of the experimental data shows that about 7% of the cation exchange capacity resides in the diffuse double layer. The model calculates the observed concentration changes very well and provides an estimate of the pristine porewater composition in the clay rock.     

Geochemical evaluation of different groundwater–host rock systems for radioactive waste disposal

The geochemical suitability of a deep bedrock repository for radioactive waste disposal is determined by the composition of geomatrix and groundwater. Both influence radionuclide solubility, chemical buffer capacity and radionuclide retention. They also determine the chemical compatibility of waste forms, containers and backfill materials. Evaluation of different groundwater–host rock systems is performed by modeling the geochemical environments and the resulting radionuclide concentrations. In order to demonstrate the evaluation method, model calculations are applied to data sets available for various geological formations such as granite, clay and rocksalt.The saturation state of the groundwater–geomatrix system is found to be fundamental for the evaluation process. Hence, calculations are performed to determine if groundwater is in equilibrium with mineral phases of the geological formation. In addition, corrosion of waste forms in different groundwater is examined by means of reaction path modeling. The corrosion reactions change the solution compositions and pH, resulting in significant changes of radionuclide solubilities. The results demonstrate that geochemical modeling of saturation state and compatibility of the host formation environment with the radioactive waste proves to be a feasible tool for evaluation of various sites considered as deep underground repositories.     

Intercomparison of reactive transport models applied to UO2 oxidative dissolution and uranium migration

Oxidative dissolution of uranium dioxide (UO(2)) and the subsequent migration of uranium in a subsurface environment and an underground waste disposal have been simulated with reactive transport models. In these systems, hydrogeological and chemical processes are closely entangled and their interdependency has been analyzed in detail, notably with respect to redox reactions, kinetics of mineralogical evolution and hydrodynamic migration of species of interest. Different codes, where among CASTEM, CHEMTRAP and HYTEC, have been used as an intercomparison and verification exercise. Although the agreement between codes is satisfactory, it is shown that the discretization method of the transport equation (i.e. finite elements (FE) versus mixed-hybrid FE and finite differences) and the sequential coupling scheme may lead to systematic discrepancies.     

Validation of the UK Met. Office’s name model against the ETEX dataset

The ETEX 1 data set has been used to assess the performance of the UK Met Office’s long-range dispersion model NAME. In terms of emergency response modelling the model performed well, successfully predicting the overall spread and timing of the plume across Europe. However, in common with most other models, NAME overpredicted the observed concentrations. This is in contrast with other NAME validation studies which indicate either no significant bias or a tendency to underpredict concentrations. This suggests the reasons for overpredicting are specific to the ETEX situation. Explanations include inadequate vertical diffusion or transport, possible venting by convective activity, and experimental errors. An assessment of a range of advection schemes of varying complexity indicated no clear advantage, at present, in using more sophisticated random walk techniques at long range, a simple diffusion coefficient based scheme providing some of the best results. A brief look is also taken at a simulation of the more problematical ETEX 2 release.     

Aquifer washing by micellar solutions: 1: Optimization of alcohol–surfactant–solvent solutions

Phase diagrams were used for the formulation of alcohol–surfactant–solvent and to identify the DNAPL (Dense Non Aqueous Phase Liquid) extraction zones. Four potential extraction zones of Mercier DNAPL, a mixture of heavy aliphatics, aromatics and chlorinated hydrocarbons, were identified but only one microemulsion zone showed satisfactory DNAPL recovery in sand columns. More than 90 sand column experiments were performed and demonstrate that: (1) neither surfactant in water, alcohol–surfactant solutions, nor pure solvent can effectively recover Mercier DNAPL and that only alcohol–surfactant–solvent solutions are efficient; (2) adding salts to alcohol–surfactant or to alcohol–surfactant–solvent solutions does not have a beneficial effect on DNAPL recovery; (3) washing solution formulations are site specific and must be modified if the surface properties of the solids (mineralogy) change locally, or if the interfacial behavior of liquids (type of oil) changes; (4) high solvent concentrations in washing solutions increase DNAPL extraction but also increase their cost and decrease their density dramatically; (5) maximum DNAPL recovery is observed with alcohol–surfactant–solvent formulations which correspond to the maximum solubilization in Zone C of the phase diagram; (6) replacing part of surfactant SAS by the alcohol n-butanol increases washing solution efficiency and decreases the density and the cost of solutions; (7) replacing part of n-butanol by the nonionic surfactant HOES decreases DNAPL recovery and increases the cost of solutions; (8) toluene is a better solvent than D-limonene because it increases DNAPL recovery and decreases the cost of solutions; (9) optimal alcohol–surfactant–solvent solutions contain a mixture of solvents in a mass ratio of toluene to D-limonene of one or two. Injection of 1.5 pore volumes of the optimal washing solution of n-butanol–SAS–toluene–D-limonene in water can recover up to 95% of Mercier DNAPL in sand columns. In the first pore volume of the washing solution recovered in the sand column effluent, the DNAPL is in a water-in-oil microemulsion lighter than the excess aqueous phase (Winsor Type II system), which indicates that part of the DNAPL was mobilized. In the next pore volumes, DNAPL is dissolved in a oil-in-water microemulsion phase and is mobilized in an excess oil phase lighter than the microemulsion (Winsor Type I system). The main drawback of this oil extraction process is the high concentration of ingredients necessary for DNAPL dissolution, which makes the process expensive. Because mobilization of oil seems to occur at the washing solution front, an injection strategy must be developed if there is no impermeable limit at the aquifer base. DNAPL recovery in the field could be less than observed in sand columns because of a smaller sweep efficiency related to field sand heterogeneities. The role of each component in the extraction processes in sand column as well as the Winsor system type have to be better defined for modeling purposes. Injection strategies must be developed to recover ingredients of the washing solution that can remain in the soil at the end of the washing process. ©1997 Elsevier Science B.V.     

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.     

Degradation of anthraquinone dye C.I. Reactive Blue 19 in aqueous solution by ozonation

This study investigated the degradation of anthraquinone reactive dye C.I. Reactive Blue 19 (RB-19) with initial concentration of 100 mg L⁻¹ in aqueous solution by ozone oxidation. The results of UV/VIS and FTIR spectra showed that the anthraquinone structures, nitrogen linkages and amino groups of RB-19 were destroyed under direct ozone reaction. The identification by LC–MS and GC–MS analyses indicated that some organic acids (e.g., phthalic acids) and 1,3-indanone could be the primary degradation products, respectively. The Microtox toxicity of the ozonated RB-19 solution initially increased but subsequently decreased when ozonation time increased. This detoxification accompanied biodegradability enhancement revealed by BOD/COD ratio increasing from 0.15 to 0.33 after 10 min of ozonation.     

Dichloroacetic acid degradation employing hydrogen peroxide and UV radiation

The degradation reaction of dichloroacetic acid employing H(2)O(2) and UVC radiation (253.7nm) has been studied in a well mixed reactor operating inside a recycling system. It has been shown that in an aqueous solution no stable reaction intermediates are formed and, at every time during the reaction, two mols of hydrochloric acid are formed for every mol of dichloroacetic acid that is decomposed and, in the same way, there is a paired agreement between the calculated TOC concentration corresponding to the unaltered dichloroacetic acid and the experimental values measured in the solution. On this basis and classical references from the scientific literature for the H(2)O(2) photolysis, a complete reaction scheme, apt for reaction kinetics mathematical modeling and ulterior scale-up is proposed.