Environmental tracer transport (H and SF6) in the saturated and unsaturated zones and its use in nitrate pollution management

An important quantity in groundwater protection is the residence time of water in an aquifer. It relates to both the travel time of a pollutant to arrive at a well and the time span required for self-purification of a polluted aquifer after removal of pollutant inputs. Time scales for aquifers can be gained from artificial tracer experiments or from environmental tracer data, the latter offering the only realistic alternative if time scales of years or decades have to be taken into account.

Different tracers show different time scales due to their different transport mechanisms especially in the unsaturated zone. While solute tracers are moved advectively with the seepage water, gas tracers pass the unsaturated zone diffusively through the air phase. Depending on the properties of the unsaturated zone (hydraulic properties, thickness) this difference in behavior can be used to separate the subsurface transport process into the unsaturated and the saturated parts.

In a field study in Germany, SF₆ and ³H were used as environmental tracers. Both have a relatively well-known input function. Interpretation of data from observation wells by a box model approach led to spatially and temporally varying residence times. This was an indication that the influence of the unsaturated zone could not be neglected. While the gas tracer SF₆ shows only residence times in the saturated zone, the tracer ³H reflects the whole travel time of water including both the unsaturated and saturated zones. Using a one-dimensional plug-flow model for the unsaturated zone combined with a detailed two-dimensional flow and transport model for the saturated zone leads to a holistic and consistent interpretation of the measured tracer concentrations. The observed pattern of old water under thick loess cover and younger water under areas where the fractured basalt aquifer crops out is reproduced after adjusting only two parameters: the effective porosity of the saturated aquifer and the product of field capacity and thickness of the unsaturated zone. While the effective porosity of the saturated zone is adjusted by means of the SF₆ data, the field capacity of the loess layer is adjusted by means of the ³H observations. The thickness of the unsaturated zone is deduced from geological and pedological maps. All flow data are obtained from a calibrated flow model, which is based on geological data, observed heads and pumping tests only.

The transport model for the saturated zone was calibrated by fitting the porosity by means of gaseous tracer concentrations (SF₆). The combined saturated–unsaturated zone model was then calibrated by fitting the field capacity of the unsaturated zone by means of ³H concentrations. With this model it was possible to verify the observed NO₃ concentrations at the drinking water wells and to develop predictions for their future development under various scenarios of fertilizer input reduction in specific areas.     

Finite-element analysis of the transport of water and solutes in title-drained soils

The finite-element method based on a Galerkin technique was used to formulate the problem of simulating the two-dimensional (cross-sectional) transient movement of water and solute in saturated or partially saturated nonuniform porous media. The numerical model utilizes linear triangular elements. Nonreactive, as well as reactive solutes whose behaviour can be described by a distribution coefficient or first-order reaction term were considered. The flow portion of the model was tested by comparison of the model results with experimental and finite-difference results for transient flow in an unsaturated sand column and the solute transport portion of the model was tested by comparison with analytical solution results. The model was applied to a hypothetical case involving movement of water and solutes in tile-drained soils. The simulation results showed the development of distinct solute leaching patterns in the soil as drainage proceeded. Although applied to a tile drainage problem in this study, the model should be equally useful in the study of a wide range of two-dimensional water and solute migration problems.     

Estimation of water saturation dependence of dispersion in unsaturated porous media: experiments and modelling analysis

In the dispersion theory, a linear relationship has been verified between the coefficient of hydrodynamic dispersion and water velocity, both in saturated and in unsaturated porous media. But for unsaturated soils the variability of flow directions and microscopic velocities can be larger than in saturated soils because of the lower degree of water saturation. This leads to an increased dispersion. Therefore, relationships between water content and relative water velocity fluctuations and water content together with the coefficient of dispersivity in unsaturated porous media respectively have been investigated systematically by displacement experiments in glass beads and coarse-textured sandy soil columns. The breakthrough curves (BTCs) of chloride showed that an increase of solute mixing with a decrease of water content was caused by an increase of flow velocity fluctuations for different pathways. In order to explain the observed tailing effect in unsaturated flow, two mathematical models were used to fit theoretically derived nonlinear functions of water content dependent dispersivities for both porous media. The close agreement between the observed and computed results suggests that the theoretical model of hydrodynamic dispersion can be extended to transport in unsaturated porous media, providing that BTCs of the effluent water are used to estimate representative dispersivity parameters of soils.     

A two‐dimensional model of hydraulic performance of stormwater infiltration systems

Stormwater infiltration systems are a popular method for urban stormwater control. They are often designed using an assumption of one‐dimensional saturated outflow, although this is not very accurate for many typical designs where two‐dimensional (2D) flows into unsaturated soils occur. Available 2D variably saturated flow models are not commonly used for design because of their complexity and difficulties with the required boundary conditions. A purpose‐built stormwater infiltration system model was thus developed for the simulation of 2D flow from a porous storage. The model combines a soil moisture–based model for unsaturated soils with a ponded storage model and uses a wetting front‐tracking approach for saturated flows. The model represents the main physical processes while minimizing input data requirements. The model was calibrated and validated using data from laboratory 2D stormwater infiltration trench experiments. Calibrations were undertaken using five different combinations of calibration data to examine calibration data requirements. It was found that storage water levels could be satisfactorily predicted using parameters calibrated with either data from laboratory soils tests or observed water level data, whereas the prediction of soil moistures was improved through the addition of observed soil moisture data to the calibration data set. Copyright © 2012 John Wiley & Sons, Ltd.     

Cyclic hysteretic flow in porous medium column: model, experiment, and simulations

A periodic vertical movement of the groundwater table results in a subsequent cyclic response of the water content and pressure profiles in the vadose zone. The sequence of periodic wetting and drying processes can be affected by hysteresis effects in this zone. A one-dimensional saturated/unsaturated flow model based on Richards’ equation and the Mualem (Soil Sci. 137 (1984) 283) hysteresis model is formulated which can take into account multi-cycle hysteresis effects in the relation between capillary pressure and water content. The numerical integration of the unsaturated flow equation is based on a Galerkin-type finite element method. The flow domain is discretised by finite elements with linear shape functions. Simulations start with static water content and pressure profiles, which correspond to either a boundary drying or wetting retention curve. To facilitate the numerical solution of the hysteretic case an implicit non-iterative procedure was chosen for the solution of the nonlinear differential equation. Laboratory experiments were performed with a vertical sand column by imposing a high frequency periodic pressure head at the lower end of the column. The total water volume in the column, and the periodic water content profile averaged over time were measured. The boundary drying and wetting curves of the relation between water content and capillary pressure were determined by independent experiments. The simulations of the experimental conditions show a clear effect of the hysteresis phenomenon on the water content profile. The simulations with hysteresis agree well with the measurements. Computed dimensionless water content profiles are presented for different oscillation frequencies with and without consideration of hysteresis.     

Implementation of Solute Transport in the Vadose Zone into the “HYDRUS Package for MODFLOW”

The “HYDRUS package for MODFLOW” is an existing MODFLOW package that allows MODFLOW to simultaneously evaluate transient water flow in both unsaturated and saturated zones. The package is based on incorporating parts of the HYDRUS-1D model (to simulate unsaturated water flow in the vadose zone) into MODFLOW (to simulate saturated groundwater flow). The coupled model is effective in addressing spatially variable saturated-unsaturated hydrological processes at the regional scale. However, one of the major limitations of this coupled model is that it does not have the capability to simulate solute transport along with water flow and therefore, the model cannot be employed for evaluating groundwater contamination. In this work, a modified unsaturated flow and transport package (modified HYDRUS package for MODFLOW and MT3DMS) has been developed and linked to the three-dimensional (3D) groundwater flow model MODFLOW and the 3D groundwater solute transport model MT3DMS. The new package can simulate, in addition to water flow in the vadose zone, also solute transport involving many biogeochemical processes and reactions, including first-order degradation, volatilization, linear or nonlinear sorption, one-site kinetic sorption, two-site sorption, and two-kinetic sites sorption. Due to complex interactions at the groundwater table, certain modifications of the pressure head (compared to the original coupling) and solute concentration profiles were incorporated into the modified HYDRUS package. The performance of the newly developed model is evaluated using HYDRUS (2D/3D), and the results indicate that the new model is effective in simulating the movement of water and contaminants in the saturated-unsaturated flow domains.     

Influence of water pressure dynamics and fluid flow on the streaming‐potential response for unsaturated conditions

Streaming‐potentials are produced by electrokinetic effects in relation to fluid flow and are used for geophysical prospecting. The aim of this study is to model streaming potential measurements for unsaturated conditions using an empirical approach. A conceptual model is applied to streaming potential measurements obtained from two drainage experiments in sand. The streaming potential data presented here show a non‐monotonous behaviour with increasing water saturation, following a pattern that cannot be predicted by existing models. A model involving quasi‐static and dynamic components is proposed to reproduce the streaming potential measurements. The dynamic component is based on the first time derivative of the driving pore pressure. The influence of this component is investigated with respect to fluid velocity, which is very different between the two experiments. The results demonstrate that the dynamic component is predominant at the onset of drainage in experiments with the slowest water flow. On the other hand, its influence appears to vanish with increasing drainage velocity. Our results suggest that fluid flow and water distribution at the pore scale have an important influence on the streaming potential response for unsaturated conditions. We propose to explain this specific streaming potential response in terms of the behaviour of both rock/water interface and water/air interfaces created during desaturation processes. The water/air interfaces are negatively charged, as also observed in the case of water/rock interfaces. Both the surface area and the flow velocity across these interfaces are thought to contribute to the non‐monotonous behaviour of the streaming potential coefficient as well as the variations in its amplitude. The non‐monotonous behaviour of air/water interfaces created during the flow was highlighted as it was measured and modelled by studies published in the literature. The streaming potential coefficient can increase to about 10 to 40 when water saturation decreases. Such an increase is possible if the amount of water/air interfaces is increased in sufficient amount, which can be the case.     

Hydrological field data from a modeller's perspective: Part 1. Diagnostic tests for model structure

Hydrological scientists develop perceptual models of the catchments they study, using field measurements and observations to build an understanding of the dominant processes controlling the hydrological response. However, conceptual and numerical models used to simulate catchment behaviour often fail to take advantage of this knowledge. It is common instead to use a pre‐defined model structure which can only be fitted to the catchment via parameter calibration. In this article, we suggest an alternative approach where different sources of field data are used to build a synthesis of dominant hydrological processes and hence provide recommendations for representing those processes in a time‐stepping simulation model. Using analysis of precipitation, flow and soil moisture data, recommendations are made for a comprehensive set of modelling decisions, including Evapotranspiration (ET) parameterization, vertical drainage threshold and behaviour, depth and water holding capacity of the active soil zone, unsaturated and saturated zone model architecture and deep groundwater flow behaviour. The second article in this two‐part series implements those recommendations and tests the capability of different model sub‐components to represent the observed hydrological processes. Copyright © 2010 John Wiley & Sons, Ltd.     

Free surface and saturated-unsaturated analyses of borehole infiltration tests above the water table

Constant head borehole infiltration tests are widely used for the in situ evaluation of saturated hydraulic conductivities of unsaturated soils above the water table. The formulae employed in analysing the results of such tests disregard the fact that some of the infiltrating water may flow under unsaturated conditions. Instead, these formulae are based on various approximations of the classical free surface theory which treats the flow region as if it were fully saturated and enclosed within a distinct envelope, the so-called ‘free surface’. A finite element model capable of solving free surface problems is used to examine the mathematical accuracy of the borehole infiltration formulae. The results show that in the hypothetical case where unsaturated flow does not exist, the approximate formulae are reasonably accurate within·a practical range of borehole conditions. To see what happens under conditions closer to those actually encountered in the field, the effect of unsaturated flow on borehole infiltration is investigated by means of two different numerical models: a mixed explicit-implicit finite element model, and a mixed explicit-implicit integrated finite difference model. Both of these models give nearly identical results; however, the integrated finite difference model is considerably faster than the finite element model. The relatively low computational efficiency of the finite element scheme is attributed to the large number of operations required in order to re-evaluate the conductivity (stiffness) matrix at each iteration in this highly non-linear saturated-unsaturated flow problem. The saturated-unsaturated analysis demonstrates that the classical free surface approach provides a distorted picture of the flow pattern in the soil. Contrary to what one would expect on the basis of this theory, only a finite region of the soil in the immediate vicinity of the borehole is saturated, whereas a significant percentage of the flow takes place under unsaturated conditions. As a consequence of disregarding unsaturated flow, the available formulae may underestimate the saturated hydraulic conductivity of fine grained soils by a factor of two, three, or more. Our saturated-unsaturated analysis leads to an improved design of borehole infiltration tests and a more accurate method for interpreting the results of such tests. The analysis also shows how one can predict the steady state rate of infiltration from data collected during the early transient period of the test.     

Prediction of diffusion coefficients in porous media using tortuosity factors based on interfacial areas

Determination of aqueous phase diffusion coefficients of solutes through porous media is essential for understanding and modeling contaminant transport. Prediction of diffusion coefficients in both saturated and unsaturated zones requires knowledge of tortuosity and constrictivity factors. No methods are available for the direct measurement of these factors, which are empirical in their definition. In this paper, a new definition for the tortuosity factor is proposed, as the real to ideal interfacial area ratio. We define the tortuosity factor for saturated porous media (tau5) as the ratio S/S(o) (specific surface of real porous medium to that of an idealized capillary bundle). For unsaturated media, tortuosity factor (tau(a)) is defined as a(aw)/a(aw),o (ratio of the specific air-water interfacial area of real and the corresponding idealized porous medium). This tortuosity factor is suitably measured using sorptive tracers (e.g., nitrogen adsorption method) for saturated media and interfacial tracers for unsaturated media. A model based on this new definition of tortuosity factors, termed the interfacial area ratio (IAR) model, is presented for the prediction of diffusion coefficients as a function of the degree of water saturation. Diffusion coefficients and diffusive resistances measured in a number of saturated and unsaturated granular porous media, for solutes in dilute aqueous solutions, agree well with the predictions of the IAR model. A comparison of permeability of saturated sands estimated based on tau(s) and the same based on the Kozeny-Carman equation confirm the usefulness of the tau(s) parameter as a measure of tortuosity.     

Non-linear quickflow response as indicators of runoff generation mechanisms

Linking quickflow response to subsurface state can improve our understanding of runoff processes that drive emergent catchment behaviour. We investigated the formation of non-linear quickflows in three forested headwater catchments and also explored unsaturated and saturated storage dynamics, and likely runoff generation mechanisms that contributed to threshold formation. Our analyses focused on two reference watersheds at the Coweeta Hydrologic Laboratory (CHL) in western North Carolina, USA, and one reference watershed at the Susquehanna Shale Hills Critical Zone Observatory (SHW) in Central Pennsylvania, USA, with available hourly soil moisture, groundwater, streamflow, and precipitation time series over several years. Our study objectives were to characterise (a) non-linear runoff response as a function of storm characteristics and antecedent conditions, (b) the critical levels of shallow unsaturated and saturated storage that lead to hourly flow response, and (c) runoff mechanisms contributing to rapidly increasing quickflow using measurements of soil moisture and groundwater. We found that maximum hourly rainfall did not significantly contribute to quickflow production in our sites, in contrast to prior studies, due to highly conductive forest soils. Soil moisture and groundwater dynamics measured in hydrologically representative areas of the hillslope showed that variable subsurface states could contribute to non-linear runoff behaviour. Quickflow generation in watersheds at CHL were dominated by both saturated and unsaturated pathways, but the relative contributions of each pathway varied between catchments. In contrast, quickflow was almost entirely related to groundwater fluctuations at SHW. We showed that co-located measurements of soil moisture and groundwater supplement threshold analyses providing stronger prediction and understanding of quickflow generation and indicate dominant runoff processes.     

Testing a conceptual hillslope recession model based on the storage–discharge relationship with the Richards equation

The conceptual recession model based on the storage–discharge relationship was proposed to account for the unsaturated–saturated water storage interaction. The recession model was formulated by combining the constitutive storage–discharge relationship with the integral balance equation for unsaturated and saturated water storage. The functional form of the constitutive storage–discharge relationship was determined from the spatial integration of the Richards equation. The performance of the recession model was tested by comparing with the solution of the Richards equation for different simulation geometric shapes and soil types. The conceptual recession model incorporating the unsaturated–saturated water storage interaction was in good agreement with the recession response of the Richards equation. However, the recession model that neglected the unsaturated–saturated water storage interaction was comparable to the Richards equation only for soils with the weak interaction between unsaturated water storage and saturated water storage. This result suggests the important role of the unsaturated–saturated water storage interaction in the formulation of the recession process when the derivative of the functional relationship between the unsaturated water storage and saturated water storage becomes significant. Copyright © 2007 John Wiley & Sons, Ltd.     

Changes in water table level influence solute transport in uniform porous media

Changes in the water table level result in variable water saturation and variable hydrological fluxes at the interface between the unsaturated and saturated zone. This may influence the transport and fate of contaminants in the subsurface. The objective of this study was to examine the impact of a decreasing and an increasing water table on solute transport. We conducted tracer experiments at downward flow conditions in laboratory columns filled with two different uniform porous media under static and transient flow conditions either increasing or decreasing the water table. Tracer breakthrough curves were simulated using a mobile–immobile transport model. The resulting transport parameters were compared to identify dominant transport processes. Changes in the water table level affected dispersivities and mobile water fractions depending on the direction of water table movement and the grain size of the porous media. In fine glass beads, the water flow velocity was similar to the decline rate of the water table, and the mobile water fraction was decreased compared with steady‐state saturated conditions. However, immobile water was negligible. In coarse glass beads, water flow was faster because of fingered flow in the unsaturated part, and the mobile water fraction was smaller than in the fine material. Here, a rising water table led to an even smaller mobile water fraction and increased solute spreading because of diffusive interaction with immobile water. We conclude that changes of the water table need to be considered to correctly simulate transport in the subsurface at the transition of the unsaturated–saturated zone. Copyright © 2014 John Wiley & Sons, Ltd.     

Modeling Variably Saturated Subsurface Solute Transport with MODFLOW‐UZF and MT3DMS

The MT3DMS groundwater solute transport model was modified to simulate solute transport in the unsaturated zone by incorporating the unsaturated‐zone flow (UZF1) package developed for MODFLOW. The modified MT3DMS code uses a volume‐averaged approach in which Lagrangian‐based UZF1 fluid fluxes and storage changes are mapped onto a fixed grid. Referred to as UZF‐MT3DMS, the linked model was tested against published benchmarks solved analytically as well as against other published codes, most frequently the U.S. Geological Survey's Variably‐Saturated Two‐Dimensional Flow and Transport Model. Results from a suite of test cases demonstrate that the modified code accurately simulates solute advection, dispersion, and reaction in the unsaturated zone. Two‐ and three‐dimensional simulations also were investigated to ensure unsaturated‐saturated zone interaction was simulated correctly. Because the UZF1 solution is analytical, large‐scale flow and transport investigations can be performed free from the computational and data burdens required by numerical solutions to Richards' equation. Results demonstrate that significant simulation runtime savings can be achieved with UZF‐MT3DMS, an important development when hundreds or thousands of model runs are required during parameter estimation and uncertainty analysis. Three‐dimensional variably saturated flow and transport simulations revealed UZF‐MT3DMS to have runtimes that are less than one tenth of the time required by models that rely on Richards' equation. Given its accuracy and efficiency, and the wide‐spread use of both MODFLOW and MT3DMS, the added capability of unsaturated‐zone transport in this familiar modeling framework stands to benefit a broad user‐ship.     

Experimental study and mathematical modelling of soluble chemical transfer from unsaturated/saturated soil to surface runoff

A one‐dimensional, two‐layer solute transport model is developed to simulate chemical transport process in an initially unsaturated soil with ponding water on the soil surface before surface runoff starts. The developed mathematical model is tested against a laboratory experiment. The infiltration and diffusion processes are mathematically lumped together and described by incomplete mixing parameters. Based on mass conservation and water balance equations, the model is developed to describe solute transport in a two‐zone layer, a ponding runoff zone and a soil mixing zone. The two‐zone layer is treated as one system to avoid describing the complicated chemical transport processes near the soil surface in the mixing zone. The proposed model was analytically solved, and the solutions agreed well with the experimental data. The developed experimental method and mathematical model were used to study the effect of the soil initial moisture saturation on chemical concentration in surface runoff. The study results indicated that, when the soil was initially saturated, chemical concentration in surface runoff was significantly (two orders of magnitude) higher than that with initially unsaturated soil, while the initial chemical concentrations at the two cases were of the same magnitude. The soil mixing depth for the initially unsaturated soil was much larger than that for the initially saturated soil, and the incomplete runoff mixing parameter was larger for the initially unsaturated soil. The higher the infiltration rate of the soil, the greater the infiltration‐related incomplete mixing parameter. According to the quantitative analysis, the soil mixing depth was found to be sensitive for both initially unsaturated and saturated soils, and the incomplete runoff mixing parameter was sensitive for initially saturated soil but not for the initially unsaturated soil; the incomplete infiltration mixing parameter behaved just the opposite. Some suggestions are made for reducing chemical loss from runoff. Copyright © 2010 John Wiley & Sons, Ltd.     

Ground-Water Tracing with Injected Helium

Helium has several characteristics that make it attractive for use as a tracer in hydrological studies. Two types of experiments were conducted to investigate applicability of helium as a tracer of ground-water movement. The experiments included studies using laboratory sand and soil columns and field ground-water tracing in a basaltic aquifer. A water helium analyzer comprised of a thin quartz glass membrane and diode ion pump (making use of the preferential permeation of helium through the quartz glass into an evacuated space) was developed and used for the experiments. Results of our studies demonstrated that breakthrough curves of specific conductance and helium were similar under saturated conditions. In the unsaturated sand/soil columns, breakthrough curves of helium were retarded relative to specific conductance reducing the usefulness of helium as a tracer.     

Numerical studies of two-dimensional saturated/ unsaturated drainage models

The partial differential equation which governs the seepage of water in unsaturated and saturated porous media is solved numerically by a generalized Newton iteration technique for two models, one ditch drainage model and one earth dam model. For each model, which is two-dimensional, a few hypothetical soils with different moisture retention curves are considered. In both models only drainage from an initially saturated soil occurs; thus, the problem of hysteresis is avoided. The results of the computations are compared with those of corresponding saturated (pure groundwater) models; solutions obtained earlier by this author and others.

Computational instability phenomena appear when the slope of the retention curves is made steep, i.e., for poorly-graded soils.     

Analysis of geomorphologic and hydrological characteristics in watershed saturated areas using topographic‐index threshold and geomorphology‐based runoff model

The objective of this paper is to investigate the variation of geomorphology and runoff characteristics in saturated areas under different partial contributing area (PCA) conditions. Geomorphologic information and hydrologic records from two mid‐size watersheds in northern Taiwan were selected for analysis. The PCA ratio in the watershed during a storm was assumed equal to the ratio of the surface‐flow volume to the direct runoff volume from measured hydrologic data. The extents of PCA regions were then determined by using a topographic‐index threshold. Consequently, the geomorphologic factors in saturated and unsaturated areas could be calculated using a digital elevation model, and these factors could then be linked to a geomorphology‐based IUH model for runoff simulation, which can consider both the surface‐ and subsurface‐flow processes in saturated and unsaturated areas, respectively. The results show that geomorphologic characteristics in the saturated areas vary significantly with different PCA ratios especially for higher order streams. A large PCA ratio results in a sharp hydrograph because the quick surface flow dominates the runoff process, whereas the hydrologic response in a low PCA case is dominated by the delayed subsurface flow. Copyright © 2007 John Wiley & Sons, Ltd.     

MRCON: COMPUTER PROGRAM FOR ESTIMATING SOIL MOISTURE RETENTION AND UNSATURATED HYDRAULIC CONDUCTIVITY

Abstract. MRCON is an interactive computer program designed to approximate soil moisture retention ( Θ-Ψ ) and/or unsaturated hydraulic conductivity (K- Θ ) for soils. MRCON, which is written in BASIC, uses empirical methods obtained from the literature to calculate K- Θ and a modified method to calculate Θ-Ψ . Input data for the program consist of a saturated moisture content and a minimum of three values of Θ-Ψ . A measured value of saturated hydraulic conductivity can be used as input to better approximate the K -Θ curves to field conditions.