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.     

Kinematic wave approximation to solute transport along preferred flow paths in soils

A routing procedure is introduced which accounts for the loss of a conservative solute tracer from preferred paths during macropore flow. Water flow is treated as a series of kinematic waves from which the tracer is lost due to mixing previously stored soil water, and an expression for solute loss is added to a previously developed model. The model parameters are estimated through experiments at three different input rates applied to a column of a macroporous forest soil.The results of seven experimental runs indicate that solute losses are consistently highest at the early stages of infiltration and drainage flow. An empirical relationship is proposed which links the frequency distribution of the flow parameter with that for solute loss from the preferred path during transient water flow and solute transport.     

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.     

Approximate evaluation of contaminant transport through vertical barriers

Containment of groundwater contamination using physical barriers can be an important element of a subsurface remediation program. This work presents simple analytical tools for predicting the performance of barriers in terms of the steady-state contaminant flux across the barrier, the duration of the transient period following barrier installation, and the time-dependent contaminant concentration distribution within the barrier. The analytical expressions are developed from approximate boundary layer (BL) solutions to the advective–dispersive equation subject to conservative fixed concentration boundary conditions. Critical ranges of important dimensionless quantities are identified for use in barrier performance assessment, for both steady-state and transient conditions. Comparative calculations made with the BL equations and more exact semi-analytical solutions are used to characterize the accuracy and applicability of the BL approach.     

Impact of varying soil structure on transport processes in different diagnostic horizons of three soil types

When soil structure varies in different soil types and the horizons of these soil types, it has a significant impact on water flow and contaminant transport in soils. This paper focuses on the effect of soil structure variations on the transport of pesticides in the soil above the water table. Transport of a pesticide (chlorotoluron) initially applied on soil columns taken from various horizons of three different soil types (Haplic Luvisol, Greyic Phaeozem and Haplic Cambisol) was studied using two scenarios of ponding infiltration. The highest infiltration rate and pesticide mobility were observed for the Bt₁ horizon of Haplic Luvisol that exhibited a well-developed prismatic structure. The lowest infiltration rate was measured for the Bw horizon of Haplic Cambisol, which had a poorly developed soil structure and a low fraction of large capillary pores and gravitational pores. Water infiltration rates were reduced during the experiments by a soil structure breakdown, swelling of clay and/or air entrapped in soil samples. The largest soil structure breakdown and infiltration decrease was observed for the Ap horizon of Haplic Luvisol due to the low aggregate stability of the initially well-aggregated soil. Single-porosity and dual-permeability (with matrix and macropore domains) flow models in HYDRUS-1D were used to estimate soil hydraulic parameters via numerical inversion using data from the first infiltration experiment. A fraction of the macropore domain in the dual-permeability model was estimated using the micro-morphological images. Final soil hydraulic parameters determined using the single-porosity and dual-permeability models were subsequently used to optimize solute transport parameters. To improve numerical inversion results, the two-site sorption model was also applied. Although structural changes observed during the experiment affected water flow and solute transport, the dual-permeability model together with the two-site sorption model proved to be able to approximate experimental data.     

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.     

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.     

Degradation of deicing chemicals affects the natural redox system in airfield soils

During winter operations at airports, large amounts of organic deicing chemicals (DIC) accumulate beside the runways and infiltrate into the soil during spring. To study the transport and degradation of DIC in the unsaturated zone, eight undisturbed soil cores were retrieved at Oslo airport, Norway, and installed as lysimeters at a nearby field site. Before snowmelt in 2010 and 2011, snow amended with a mix of the DICs propylene glycol (PG) and formate as well as bromide as conservative tracer was applied. Water samples were collected and analyzed until summer 2012. Water flow and solute transport varied considerably among the lysimeters but also temporally between 2010 and 2011. High infiltration rates during snowmelt resulted in the discharge of up to 51 and 82 % PG in 2010 and 2011, respectively. The discharge of formate remained comparatively low, indicating its favored degradation even at freezing temperatures compared with PG. Manganese (Mn) and iron (Fe) were observed in the drainage in autumn owing to the anaerobic degradation of residual PG during summer. Our findings suggest that upper boundary conditions, i.e., snow cover and infiltration rate, and the extent of preferential flowpaths, control water flow and solute transport of bromide and PG during snowmelt. PG may therefore locally reach deeper soil regions where it may pose a risk for groundwater. In the long term, the use of DIC furthermore causes the depletion of potential electron acceptors and the transport of considerable amounts of Fe and Mn. To avoid an overload of the unsaturated zone with DIC and to maintain the natural redox system, the development of suitable remediation techniques is required.     

Effect of soil moisture dynamics on dense nonaqueous phase liquid (DNAPL) spill zone architecture in heterogeneous porous media

The amount, location, and form of NAPL in contaminated vadose zones are controlled by the spatial distribution of water saturation and soil permeability, the NAPL spill scenario, water infiltration events, and vapor transport. To evaluate the effects of these processes, we used the three-phase flow simulator STOMP, which includes a new permeability-liquid saturation-capillary pressure (k-S-P) constitutive model. This new constitutive model considers three NAPL forms: free, residual, and trapped. A 2-D vertical cross-section with five stratigraphic layers was assumed, and simulations were performed for seven cases. The conceptual model of the soil heterogeneity was based upon the stratigraphy at the Hanford carbon tetrachloride (CT) spill site. Some cases considered co-disposal of NAPL with large volumes of wastewater, as also occurred at the Hanford CT site. In these cases, the form and location of NAPL were most strongly influenced by high water discharge rates and NAPL evaporation to the atmosphere. In order to investigate the impact of heterogeneity, the hydraulic conductivity within the lower permeability layer was modeled as a realization of a random field having three different classes. For six extreme cases of 100 realizations, the CT mass that reached the water table varied by a factor of two, and was primarily controlled by the degree of lateral connectivity of the low conductivity class within the lowest permeability layer. The grid size at the top boundary had a dramatic impact on NAPL diffusive flux just after the spill event when the NAPL was present near the ground surface. NAPL evaporation with a fine grid spacing at the top boundary decreased CT mass that reached the water table by 74%, compared to the case with a coarse grid spacing, while barometric pumping had a marginal effect for the case of a continuous NAPL spill scenario considered in this work. For low water infiltration rate scenarios, the distribution of water content prior to a NAPL spill event decreased CT mass that reached the water table by 98% and had a significant impact on the formation of trapped NAPL. For all cases simulated, use of the new constitutive model that allows the formation of residual NAPL increased the amount of NAPL retained in the vadose zone. Density-driven advective gas flow from the ground surface controlled vapor migration in strongly anisotropic layers, causing NAPL mass flux to the lower layer to be reduced. These simulations indicate that consideration of the formation of residual and trapped NAPLs and dynamic boundary conditions (e.g., areas, rates, and periods of different NAPL and water discharge and fluctuations of atmospheric pressure) in the context of full three-phase flow are needed, especially for NAPL spill events at the ground surface. In addition, NAPL evaporation, density-driven gas advection, and NAPL vertical movement enhanced by water flow must be considered in order to predict NAPL distribution and migration in the vadose zone.     

Linking a one-dimensional pesticide fate model to a three-dimensional groundwater model to simulate pollution risks of shallow and deep groundwater underlying fractured till

The one-dimensional pesticide fate model MACRO was loose-linked to the three-dimensional discrete fracture/matrix diffusion model FRAC3DVS to describe transport of the pesticide mecoprop in a fractured moraine till and local sand aquifer (5-5.5 m depth) overlying a regional limestone aquifer (16 m depth) at Havdrup, Denmark. Alternative approaches to describe the upper boundary in the groundwater model were examined. Field-scale simulations were run to compare a uniform upper boundary condition with a spatially variable upper boundary derived from Monte-Carlo simulations with MACRO. Plot-scale simulations were run to investigate the influence of the temporal resolution of the upper boundary conditions for fluxes in the groundwater model and the effects of different assumptions concerning the macropore/fracture connectivity between the two models. The influence of within-field variability of leaching on simulated mecoprop concentrations in the local aquifer was relatively small. A fully transient simulation with FRAC3DVS gave 20 times larger leaching to the regional aquifer compared to the case with steady-state water flow, assuming full connectivity with respect to macropores/fractures across the boundary between the two models. For fully transient simulations 'disconnecting' the macropores/fractures at the interface between the two models reduced leaching by a factor 24. A fully connected, transient simulation with FRAC3DVS, with spatially uniform upper boundary fluxes derived from a MACRO simulation with 'effective' parameters is therefore recommended for assessing leaching risks to the regional aquifer, at this, and similar sites.     

Electrical resistivity tomography as monitoring tool for unsaturated zone transport: an example of preferential transport of deicing chemicals

Non-invasive spatially resolved monitoring techniques may hold the key to observe heterogeneous flow and transport behavior of contaminants in soils. In this study, time-lapse electrical resistivity tomography (ERT) was employed during an infiltration experiment with deicing chemical in a small field lysimeter. Deicing chemicals like potassium formate, which frequently impact soils on airport sites, were infiltrated during snow melt. Chemical composition of seepage water and the electrical response was recorded over the spring period 2010. Time-lapse electrical resistivity tomographs are able to show the infiltration of the melt water loaded with ionic constituents of deicing chemicals and their degradation product hydrogen carbonate. The tomographs indicate early breakthrough behavior in parts of the profile. Groundtruthing with pore fluid conductivity and water content variations shows disagreement between expected and observed bulk conductivity. This was attributed to the different sampling volume of traditional methods and ERT due to a considerable fraction of immobile water in the soil. The results show that ERT can be used as a soil monitoring tool on airport sites if assisted by common soil monitoring techniques.     

Analytical solutions for reactive transport under an infiltration-redistribution cycle

Transport of reactive solute in unsaturated soils under an infiltration-redistribution cycle is investigated. The study is based on the model of vertical flow and transport in the unsaturated zone proposed by Indelman et al. [J. Contam. Hydrol. 32 (1998) 77], and generalizes it by accounting for linear nonequilibrium kinetics. An exact analytical solution is derived for an irreversible desorption reaction. The transport of solute obeying linear kinetics is modeled by assuming equilibrium during the redistribution stage. The model which accounts for nonequilibrium during the infiltration and assumes equilibrium at the redistribution stage is termed partial equilibrium infiltration-redistribution model (PEIRM). It allows to derive approximate closed form solutions for transport in one-dimensional homogeneous soils. These solutions are further applied to computing the field-scale concentration by adopting the Dagan and Bresler [Soil Sci. Soc. Am. J. 43 (1979) 461] column model. The effect of soil heterogeneity on the solute spread is investigated by modeling the hydraulic saturated conductivity as a random function of horizontal coordinates. The quality of the PEIRM is illustrated by calculating the critical values of the Damk?hler number which provide the achievable accuracy in estimating the solute mass in the mobile phase. The distinguishing feature of transport during the infiltration-redistribution cycle as compared to that of infiltration only is the finite depth of solute penetration. For irreversible desorption, the maximum solute penetration W/theta(r) is determined by the amount of applied water W and the residual water content theta(r). For sorption-desorption kinetics, the maximum depth of penetration z(r)(e, infinity ) also depends on the ratio between the rate of application and the column-saturated conductivity. It is shown that z(r)(e, infinity ) is bounded between the depths W/(theta(r)+K(d)) and W/theta(r) corresponding to the maximum solute penetration for equilibrium transport and for irreversible desorption, respectively. This feature of solute penetration explains the unusual phenomena of plume contraction after an initial period of spreading [Lessoff, S.C., Indelman, P., Dagan, G., 2002. Solute transport in infiltration-redistribution cycles in heterogeneous soils. In Raats, P.A.C., Smiles, D.,Warrick, A.W. (Eds), Environmental Mechanics: Water, Mass and Energy Transport in the Biosphere. American Geophysical Union, pp. 133-144]. Unlike transport under equilibrium conditions, when the solute is completely concentrated at the front, the solute under nonequilibrium conditions is spread out behind the front. Heterogeneity leads to additional spreading of the plume.     

Second-order modeling of arsenite transport in soils

Rate limited processes including kinetic adsorption-desorption can greatly impact the fate and behavior of toxic arsenic compounds in heterogeneous soils. In this study, miscible displacement column experiments were carried out to investigate the extent of reactivity during transport of arsenite in soils. Arsenite breakthrough curves (BTCs) of Olivier and Windsor soils exhibited strong retardation with diffusive effluent fronts followed by slow release or tailing during leaching. Such behavior is indicative of the dominance of kinetic retention reactions for arsenite transport in the soil columns. Sharp decrease or increase in arsenite concentration in response to flow interruptions (stop-flow) further verified that non-equilibrium conditions are dominant. After some 40-60 pore volumes of continued leaching, 30-70% of the applied arsenite was retained by the soil in the columns. Furthermore, continued arsenite slow release for months was evident by the high levels of residual arsenite concentrations observed during leaching. In contrast, arsenite transport in a reference sand material exhibited no retention where complete mass recovery in the effluent solution was attained. A second-order model (SOM) which accounts for equilibrium, reversible, and irreversible retention mechanisms was utilized to describe arsenite transport results from the soil columns. Based on inverse and predictive modeling results, the SOM model successfully depicted arsenite BTCs from several soil columns. Based on inverse and predictive modeling results, a second-order model which accounts for kinetic reversible and irreversible reactions is recommended for describing arsenite transport in soils.     

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.     

Multi-process herbicide transport in structured soil columns: experiments and model analysis

Model predictions of pesticide transport in structured soils are complicated by multiple processes acting concurrently. In this study, the hydraulic, physical, and chemical nonequilibrium (HNE, PNE, and CNE, respectively) processes governing herbicide transport under variably saturated flow conditions were studied. Bromide (Br-), isoproturon (IPU, 3-(4-isoprpylphenyl)-1,1-dimethylurea) and terbuthylazine (TER, N2-tert-butyl-6-chloro-N4-ethyl-1,3,5-triazine-2,4-diamine) were applied to two soil columns. An aggregated Ap soil column and a macroporous, aggregated Ah soil column were irrigated at a rate of 1 cm h(-1) for 3 h. Two more irrigations at the same rate and duration followed in weekly intervals. Nonlinear (Freundlich) equilibrium and two-site kinetic sorption parameters were determined for IPU and TER using batch experiments. The observed water flow and Br- transport were inversely simulated using mobile-immobile (MIM), dual-permeability (DPM), and combined triple-porosity (DP-MIM) numerical models implemented in HYDRUS-1D, with improving correspondence between empirical data and model results. Using the estimated HNE and PNE parameters together with batch-test derived equilibrium sorption parameters, the preferential breakthrough of the weakly adsorbed IPU in the Ah soil could be reasonably well predicted with the DPM approach, whereas leaching of the strongly adsorbed TER was predicted less well. The transport of IPU and TER through the aggregated Ap soil could be described consistently only when HNE, PNE, and CNE were simultaneously accounted for using the DPM. Inverse parameter estimation suggested that two-site kinetic sorption in inter-aggregate flow paths was reduced as compared to within aggregates, and that large values for the first-order degradation rate were an artifact caused by irreversible sorption. Overall, our results should be helpful to enhance the understanding and modeling of multi-process pesticide transport through structured soils during variably saturated water flow.     

Importance of considering intraborehole flow in solute transport modeling under highly dynamic flow conditions

Correct interpretation of tracer test data is critical for understanding transport processes in the subsurface. This task can be greatly complicated by the presence of intraborehole flows in a highly dynamic flow environment. At a new tracer test site (Hanford IFRC) a dynamic flow field created by changes in the stage of the adjacent Columbia River, coupled with a heterogeneous hydraulic conductivity distribution, leads to considerable variations in vertical hydraulic gradients. These variations, in turn, create intraborehole flows in fully-screened (6.5m) observation wells with frequently alternating upward and downward movement. This phenomenon, in conjunction with a highly permeable aquifer formation and small horizontal hydraulic gradients, makes modeling analysis and model calibration a formidable challenge. Groundwater head data alone were insufficient to define the flow model boundary conditions, and the movement of the tracer was highly sensitive to the dynamics of the flow field. This study shows that model calibration can be significantly improved by explicitly considering (a) dynamic flow model boundary conditions and (b) intraborehole flow. The findings from this study underscore the difficulties in interpreting tracer tests and understanding solute transport under highly dynamic flow conditions.     

Comparison of a two-dimensional numerical dust transport model with experimental dust emissions from soil surfaces in a wind tunnel

Near-surface wind-tunnel fugitive dust concentration profiles arising from soil surfaces beds were compared to a finite difference numerical dust transport model. Comparisons of the type shown in this study were previously non-existent in the literature due to the lack of experimental wind-tunnel data for near-surface concentrations over a soil bed. However, in a previous study by the authors, near-surface steady-state concentration profiles were measured in order to obtain fugitive dust emission rates, thus allowing the comparison to models shown in this paper. The novel aspects of the current study include: comparison of concentration profiles of dust obtained experimentally in the wind tunnel with those calculated numerically; comparison of the calculated numerical fetch effect on dust emissions with that obtained in the wind tunnel; and comparison of the emission rates calculated numerically with those obtained experimentally in the wind tunnel. Initial comparisons with the model indicate good agreement implying that the physical mechanism of advection–diffusion is reasonably modeled with the choice of equations for the simple “steady-state” process near the surface. Furthermore, the numerical solutions presented in this paper provide a means to systematically explore the relative impact of varied surface boundary conditions upon the emission process and provide a potential link between wind-tunnel simulations and field scale models.