An alternative Boussinesq equation considering the effect of hysteresis on coastal groundwater waves

This paper presents an alternative Boussinesq equation considering hysteresis effect via a third‐order derivative term. By introducing an improved moisture–pressure retention function, this equation describes, with reasonable precision, groundwater propagation in coastal aquifers subject to Dirichlet boundary condition of different oscillation frequencies. Test results confirmed that it is necessary to consider horizontal and vertical flows in unsaturated zone, because of their variable influences on hysteresis. Hysteresis in unsaturated zone can affect the water table wave number of groundwater wave motion, such as wave damping rate and phase lag. Oscillations with different periods exert different hysteresis effect on wave propagation. Truncation/shrinkage of unsaturated zones also affects the strength of hysteresis. These impacts can be reflected in the alternative Boussinesq equation by adjusting the parameter representing the variation rate of moisture associated with pressure change, as opposed to traditional computationally expensive hysteresis algorithms. The present Boussinesq equation is simple to use and can provide feasible basis for future coupling of groundwater and surface water models. Copyright © 2016 John Wiley & Sons, Ltd.     

Characterizing the hydraulic properties of unsaturated sand considering various porosities and drying‐wetting paths of the retention curve

Previous studies have shown that water retention curve (WRC) and the hydraulic conductivity vary because of changes of the void ratio or porosity of soil. However, limited documents pointed out the change of hydraulic properties of soil when compacted to different porosities while considering both of the drying and wetting processes of the WRC. This information is sometimes necessary for research like finger flow analysis or the occurrence of wetting and drying cycles as what would be seen in the field. Therefore, this study aims to examine the change of WRC characteristics with varied porosity considering both of the drying and wetting path in WRC by conducting a sand box experiment. Results show that the same type of sand compacted to various porosities have different hydraulic parameters. Hydraulic conductivities generally decrease with reduced porosities; shape parameter α of the van Genuchten equation (1980) linearly decreases with declining porosity and shape parameter n in a reversal manner for the sands of interest whether in the drying process or wetting process. The unsaturated properties of sand are further characterized by inspecting the variations of moisture content, matric suction and vertical displacement of soil body subject to periodic changes of the water level by another sand box experiment. The outcomes suggest that the saturated water content and residual water content are changing during the wetting–drying process, which can be an implication of the changed properties of WRC. The characteristics of volumetric deformation might be varied as well because of the observation of the dissimilar patterns of the changing vertical displacements among each wetting–drying process. Infiltration patterns of the sands also are identified through numerical modelling by introducing a constant infiltration flux from the surface followed by a no‐influx condition. Results indicate that less water accumulates in the sand near the surface for the sand compacted to higher porosity, but water can move deeper. Hydraulic conductivity is found as the prime factor dominating the evolvement of wetting fronts. However, shape parameters of water retention curves also affect the infiltration pattern to some extent. In addition, different sands with similar porosities can have quite different infiltrating characteristics. Copyright © 2012 John Wiley & Sons, Ltd.     

The influence of water retention curve hysteresis on the stability of unsaturated soil slopes

Hysteresis is a common feature exhibited in hydraulic properties of an unsaturated soil. The movement of wetting front and the hysteresis effect are important factors which impact the shear strength of the unsaturated soil and the mechanics of shallow landslides. These failures are mainly triggered by the deepening of the wetting front accompanied by a decrease in matric suction induced by infiltration. This research establishes a method for determining a stability analysis of unsaturated infinite soil slopes, integrating the influence of infiltration and the water retention curve hysteresis. Furthermore, the present stability analysis method including the infiltration model and the advanced Mohr–Coulomb failure criterion calculates the variations of the safety factor (FS) in accordance with different slope angle, depth and hydrological processes. The experimentally measured data on the effect of hysteresis are also carried out for comparison. Numerical analyses, employing both wetting and drying hydraulic behaviour of unsaturated soil, are performed to study the difference in soil‐water content as observed in the experiments. The simulating approximations also fully responded to the experimental data of sand box. The results suggest that the hysteresis behaviour affect the distribution of soil‐water content within the slope indeed. The hysteresis made the FS values a remarkable recovery during the period of non‐rainfall in a rainfall event. The appropriate hydraulic properties of soil (i.e. wetting or drying) should be used in accordance with the processes that unsaturated soil actually experience. This method will enable us to acquire more accurate matric suction head and the unsaturated soil‐shear strength as it changes with the hysteretic flow, in order to calculate into the stability analysis of shallow landslides. An advanced understanding of the process mechanism afforded by this method is critical to realizing a reliable and appropriate design for slope stabilization. It also offers some immediate reference information to the disaster reduction department of the government. Copyright © 2011 John Wiley & Sons, Ltd.     

Measurement and modeling of moisture content above an oscillating water table: implications for beach surface moisture dynamics

This study examined the influence of tidally‐induced oscillations of the beach water table in regulating beach surface moisture dynamics. A series of laboratory experiments were conducted to assess the influence of hysteresis and transient flow effects on surface moisture variability. The experimental apparatus utilized a column of well‐sorted fine sand partially immersed in a reservoir of water. The water level in the reservoir was raised and lowered via a diaphragm‐metering pump to simulate tidally induced fluctuations of the water table, and the moisture content profile within the column was monitored using an array of Delta‐T probes. Moisture contents at specific elevations within the column were utilized as proxies to represent various ‘surface’ elevations (relative to the high water table). Results indicate that surface moisture content behaves in a distinctly hysteretic manner. Examination of water flow scanning curves illustrated that for all surface elevations considered, higher moisture contents for a given pressure head occurred during the drying cycle than during the wetting cycle. This observation is particularly evident with shallow surface elevations (i.e. water table close to the surface) where the Haines Jump phenomenon was found to have a significant influence on moisture content dynamics. Additionally, an assessment of the accuracy of hysteretic and non‐hysteretic models to predict the measured moisture contents demonstrated that hysteretic simulations consistently provide a better representation of the observed moisture contents than non‐hysteretic simulations. A time lag was found between the respective maxima and minima in water table elevation surface moisture content. At the near surface water table positions the time lag ranged between 30 and 100 minutes, and it increased to 240 minutes (four hours) with the high water table at 60 cm below the surface. Copyright © 2013 John Wiley & Sons, Ltd.     

Drainage in finite-sized unsaturated zones

After the initiation of gravity drainage, water is often assumed to be either (a) draining under unit gradient, or (b) at capillary/gravity equilibrium. Both of these simplifications can be useful, but the regimes of validity of each assumption must be delineated. Water pressures are measured versus time and distance as water drains out of a 1.6 m long sand column to determine the relative effects of capillary and gravitational forces during drainage. For medium sized sands (0.15–0.3 mm in diameter), the capillary pressure is constant in space in a large region of the column for over 12 days, and the water continues to flow under unit gradient for relatively long time scales. Similar results are seen for finer sands, but with a much faster approach to equilibrium. Numerical simulations and analytical estimates are presented and compare favorably to the measurements. Together, the experimental, theoretical and analytical results are used to calculate when capillary/gravity equilibrium is reached as a function of porous media properties and length of the unsaturated zone. The ratio of the length of the unsaturated zone to the bubbling pressure is a key parameter in determining the drainage regime, and that even for relatively short unsaturated zones the equilibrium time scale can be on the order of years.     

Analytical Solutions for Water Flow and Solute Transport in the Unsaturated Zone

Analytical solutions for the water flow and solute transport equations in the unsaturated zone are presented. We use the Broadbridge and White nonlinear model to solve the Richards’ equation for vertical flow under a constant infiltration rate. Then we extend the water flow solution and develop an exact parametric solution for the advection-dispersion equation. The method of characteristics is adopted to determine the location of a solute front in the unsaturated zone. The dispersion component is incorporated into the final solution using a singular perturbation method. The formulation of the analytical solutions is simple, and a complete solution is generated without resorting to computationally demanding numerical schemes. Indeed, the simple analytical solutions can be used as tools to verify the accuracy of numerical models of water flow and solute transport. Comparison with a finite-element numerical solution indicates that a good match for the predicted water content is achieved when the mesh grid is one-fourth the capillary length scale of the porous medium. However, when numerically solving the solute transport equation at this level of discretization, numerical dispersion and spatial oscillations were significant.     

An isotopic investigation of a near-surface groundwater system

ABSTRACT

The distribution of environmental tritium, deuterium and oxygen-18 in the unsaturated zone and the underlying sandy phreatic aquifer was studied throughout 1981 in an area of high pine forests in the Rhine valley near Heidelberg. The observed vertical distribution of isotopes in the unsaturated zone can be satisfactorily explained by the combined use of a multi-cell model for moisture transport and an evapotranspiration model. The distribution in the underlying aquifer of the tracer input at the water table obtained using this method is found by replacing the total vertical diffusion coefficient in the diffusion equation with the dispersion coefficient. In this way observed tritium profiles are satisfactorily simulated for the period 1966–1981. The stable isotope profile in the unsaturated zone however remains largely unexplained due to inadequate data on the stable isotope content of precipitation over the investigated area.     

Motion of wetting fronts moving into partially pre-wet soil

We study the motion of wetting fronts for vertical infiltration problems as modeled by Richards’ equation. Parlange and others have shown that wetting fronts in infiltration flows can be described by traveling wave solutions. If the soil layer is not initially dry, but has an initial distribution of water content then the motion of the wetting front will change due to the interaction of the infiltrating flow with the pre-existing soil conditions. Using traveling wave profiles, we construct simple approximate solutions of initial-boundary value problems for Richards’ equation that accurately describe the position and moisture distribution of the wetting front. We show that the influences of surface boundary conditions and initial conditions produce shifts to the position of the wetting front. The shifts can be calculated by examining the cumulative infiltration, and are validated numerically for several problems for Richards’ equation and the linear advection–diffusion equation.     

Importance of unsaturated zone flow for simulating recharge in a humid climate

Transient recharge to the water table is often not well understood or quantified. Two approaches for simulating transient recharge in a ground water flow model were investigated using the Trout Lake watershed in north-central Wisconsin: (1) a traditional approach of adding recharge directly to the water table and (2) routing the same volume of water through an unsaturated zone column to the water table. Areas with thin (less than 1 m) unsaturated zones showed little difference in timing of recharge between the two approaches; when water was routed through the unsaturated zone, however, less recharge was delivered to the water table and more discharge occurred to the surface because recharge direction and magnitude changed when the water table rose to the land surface. Areas with a thick (15 to 26 m) unsaturated zone were characterized by multimonth lags between infiltration and recharge, and, in some cases, wetting fronts from precipitation events during the fall overtook and mixed with infiltration from the previous spring snowmelt. Thus, in thicker unsaturated zones, the volume of water infiltrated was properly simulated using the traditional approach, but the timing was different from simulations that included unsaturated zone flow. Routing of rejected recharge and ground water discharge at land surface to surface water features also provided a better simulation of the observed flow regime in a stream at the basin outlet. These results demonstrate that consideration of flow through the unsaturated zone may be important when simulating transient ground water flow in humid climates with shallow water tables.     

A quasi three-dimensional model for flow and transport in unsaturated and saturated zones: 1. Implementation of the quasi two-dimensional case

A quasi three-dimensional (QUASI 3-D) model is presented for simulating the subsurface water flow and solute transport in the unsaturated and in the saturated zones of soil. The model is based on the assumptions of vertical flow in the unsaturated zone and essentially horizontal groundwater flow. The 1-D Richards equation for the unsaturated zone is coupled at the phreatic surface with the 2-D flow equation for the saturated zone. The latter was obtained by averaging 3-D flow equation in the saturated zone over the aquifer thickness. Unlike the Boussinesq equation for a leaky-phreatic aquifer, the developed model does not contain a storage term with specific yield and a source term for natural replenishment. Instead it includes a water flux term at the phreatic surface through which the Richards equation is linked with the groundwater flow equation. The vertical water flux in the saturated zone is evaluated on the basis of the fluid mass balance equation while the horizontal fluxes, in that equation, are prescribed by Darcy law. A 3-D transport equation is used to simulate the solute migration. A numerical algorithm to solve the problem for the general quasi 3-D case was developed. The developed methodology was exemplified for the quasi 2-D cross-sectional case (QUASI2D). Simulations for three synthetic problems demonstrate good agreement between the results obtained by QUASI2D and two fully 2-D flow and transport codes (SUTRA and 2DSOIL). Yet, simulations with the QUASI2D code were several times faster than those by the SUTRA and the 2DSOIL codes.     

The soil-moisture zone in a physically-based hydrologic model

The present paper describes an approach to modelling the unsaturated soil-moisture zone in the framework of an integrated physically-based hydrologic response model. It is supposed that the subsurface flow regime may be viewed as two separate entities — a saturated flow system which may be modelled by standard two-dimensional regional techniques, and a single overlying unsaturated zone in which the flow is essentially vertical. Coupling takes place via the definition of saturation at the lower boundary of the unsaturated zone, and via a conservative water balance. Attention is focused on the computational procedure for the unsaturated zone as a self-contained module. The major difficulties are the definition of the interface between the saturated and unsaturated zones, the nonlinear character of the equation used to describe unsaturated flow, the inclusion of realistic atmospheric boundary conditions, and, the interaction between water uptake by plants and available soil-moisture. Each of these points is discussed, in turn, with the emphasis on mathematically formulating the problem in such a way that the most important physical features are reproduced with a minimal amount of computational effort. The text concludes with a few illustrative examples.     

Nonequilibrium statistical mechanical derivation of a nonlocal Darcy's Law for unsaturated/saturated flow

As illustrated variously by wetting and drying scanning curves, flow in unsaturated porous media is inherently nonlocal. This nonlocality is also manifest in hysteresis in the classical Darcy conductivity. It is the authors' belief that most current theories of unsaturated/saturated flow are often inadequate, as they do not account for spatial nonlocality and memory. Here we provide a fundamental theory in which nonlocality of the flow constitutive theory is a natural consequence of force balances. The results are derived from general principles in statistical physics and under appropriate limiting conditions, the classical Darcy's Law is recovered for saturated flow. A notable departure in this theory from other nonlocal flow theories is that a classical Darcy type equation on a small scale need not exist.     

Influence of soil moisture hysteresis on the functioning of capillary barriers

Previous experimental studies of capillary barriers have identified highly hysteretic soil moisture retention characteristics in the materials used. In this study, numerical modelling is used to analyse the role of soil moisture hysteresis in capillary barrier functioning. Comparisons between published experimental results and model simulations indicate that soil moisture hysteresis was a necessary inclusion in the modelling approach to adequately reproduce pore water pressure distributions and the timing of breakthrough occurrences. Under hypothetical intermittent infiltration and evaporation conditions, the predicted volumetric water content in the moisture retention layer was significantly different for hysteretic and non‐hysteretic models. The hysteresis effect was found to be dependent on the nature of infiltration–evaporation cycling, although the predicted volume of flow through the hysteretic barrier was lower than that of the non‐hysteretic case, regardless of the nature of the cyclic upper boundary conditions. For practical engineering designs, where the water leakage through the barrier is the primary concern, the inclusion of soil moisture hysteresis in numerical modelling is needed. Copyright © 2009 John Wiley & Sons, Ltd.     

A wetting and drying algorithm with a combined pressure/free-surface formulation for non-hydrostatic models

A wetting and drying method for free-surface problems for the three-dimensional, non-hydrostatic Navier–Stokes equations is proposed. The key idea is to use a horizontally fixed mesh and to apply different boundary conditions on the free-surface in wet and dry zones. In wet areas a combined pressure/free-surface kinematic boundary condition is applied, while in dry areas a positive water level and a no-normal flow boundary condition are enforced. In addition, vertical mesh movement is performed to accurately represent the free-surface motion. Non-physical flow in the remaining thin layer in dry areas is naturally prevented if a Manning–Strickler bottom drag is used. The treatment of the wetting and drying processes applied through the boundary condition yields great flexibility to the discretisation used. Specifically, a fully unstructured mesh with any finite element choice and implicit time discretisation method can be applied. The resulting method is mass conservative, stable and accurate. It is implemented within Fluidity-ICOM [1] and verified against several idealized test cases and a laboratory experiment of the Okushiri tsunami.     

Non-isothermal soil water transport and evaporation

A detailed model was formulated to describe the non-isothermal transport of water in the unsaturated soil zone. The model consists of the coupled equations of mass conservation for the liquid phase, gas phase and water vapor and the energy conservation equation. The water transport mechanisms considered are convection in the liquid phase, and convection, diffusion and dispersion of vapor in the gas phase. The boundary conditions at the soil–atmosphere interface include dynamical mass flux and energy flux that accounts for radiation transport. Comparison of numerical simulations results with published experimental data demonstrated that the present model is able to describe water and energy transport dynamics, including situations of low and moderate soil moisture contents. Analysis of field studies on soil drying suggests that that dispersion flux of the water vapor near the soil surface, which is seldom considered in soil drying models, can make a significant contribution to the total water flux.     

Measuring the effect of structural connectivity on the water dynamics in heterogeneous porous media using speedy neutron tomography

The temporal and spatial distribution of water within a porous medium is affected by the medium’s structure, i.e., the spatial arrangement of its constituents. To analyze structural effects on the fluid dynamics, we measured the 3D water content distribution in a heterogeneous sand column during two drainage-wetting cycles using neutron transmission tomography. The sample with a volume of 105 cm³ contained 101 cubes of fine and 49 cubes of coarse sand with particles ranging from 0.01 to 0.05 and 0.03 to 0.09 cm, respectively. The pressure at the lower boundary was determined by the water reservoir positioned between 7 and 39 cm below the top of the column. The duration of one complete 3D scanning with a spatial resolution of 127 μm was 56 s. The signal to noise ratio of the measurements was low due to the short exposure time in the neutron beam, but it was possible to quantify the water content in the individual cubes and hence the effect of structure on macroscopic water distribution. Continuous structures of coarse sand drained faster than coarse sand without connection to the upper boundary. During the initial wetting phase, cubes of coarse sand material completely embedded in the fine material remained water unsaturated due to air entrapment. The effect of the coarse sand connectivity was analyzed in two-dimensional numerical simulations based on Richards equation. In contrast to the measurements, no effect of structure connectivity was found. The coarse sand cubes embedded within the fine matrix drain as quickly as the coarse sand cubes arranged in a continuous channel due to the model assumption of a continuous air phase.     

Numerical analysis of one-dimensional unsaturated flow in layered soils

The paper deals with numerical solutions to the Richards equation to simulate one-dimensional flow processes in the unsaturated zone of layered soil profiles. The equation is expressed in the pressure-based form and a finite-difference algorithm is developed for accurately estimating the values of the hydraulic conductivity between two neighboring nodes positioned in different soil layers, often referred to as the interlayer hydraulic conductivity. The algorithm is based upon flux conservation and continuity of pressure potential at the interface between two consecutive layers, and does not add significantly to simulation run time. The validity of the model is established for a number of test problems by comparing numerical results with the analytical solutions developed by Srivastava and Yeh²⁹ which hold for vertical infiltration towards the water table in a two-layer soil profile. The results show a significant reduction in relative mass balance errors when using the proposed model. Some specific insights into its numerical performance are also gained by comparisons with a numerical model in which the more common geometric averaging operator acts on the interlayer conductivities.     

Inverse modeling of a radial multistep outflow experiment for determining unsaturated hydraulic properties

Modeling flow and solute transport in the unsaturated zone on the basis of the Richards equation requires specifying values for unsaturated hydraulic conductivity and water potential as a function of saturation. The objectives of the paper are to evaluate the design of a transient, radial, multi-step outflow experiment, and to determine unsaturated hydraulic parameters using inverse modeling. We conducted numerical simulations, sensitivity analyses, and synthetic data inversions to assess the suitability of the proposed experiment for concurrently estimating the parameters of interest. We calibrated different conceptual models against transient flow and pressure data from a multi-step, radial desaturation experiment to obtain estimates of absolute permeability, as well as the parameters of the relative permeability and capillary pressure functions. We discuss the differences in the estimated parameter values and illustrate the impact of the underlying model on the estimates. We demonstrate that a small error in absolute permeability, if determined in an independent experiment, leads to biased estimates of unsaturated hydraulic properties. Therefore, we perform a joint inversion of pressure and flow rate data for the simultaneous determination of permeability and retention parameters, and analyze the correlations between these parameters. We conclude that the proposed combination of a radial desaturation experiment and inverse modeling is suitable for simultaneously determining the unsaturated hydraulic properties of a single soil sample, and that the inverse modeling technique provides the opportunity to analyze data from nonstandard experimental designs.     

Non-passive transport of volatile organic compounds in the unsaturated zone

A detailed model was formulated to describe the non-passive transport of water-soluble chemicals in the unsaturated zone and used to illustrate one-dimensional infiltration and redistribution of alcohol–water mixtures. The model includes the dependence of density, viscosity, surface tension, molecular diffusion coefficient in the liquid-phase, and gas–liquid partition coefficient on the aqueous mixture composition. It also takes into account the decrease in the gas–liquid partition coefficient at high capillary pressures, in accordance with Kelvin’s equation for multi-component mixtures. Simulation of butanol–water mixtures infiltration in sand was in agreement with the experimental data and simulations reported in the literature. Simulation of methanol infiltration and redistribution in two different soils showed that methanol concentration significantly affects volumetric liquid content and concentration profiles, as well as the normalized volatilization and evaporation fluxes. Dispersion in the liquid-phase was the predominant mechanism in the transport of methanol when dispersivity at saturation was set to 7.8 cm. Liquid flow was mainly due to capillary pressure gradients induced by changes in volumetric liquid content. However, for dispersivity at saturation set to 0.2 cm, changes in surface tension due to variation in composition induced important liquid flow and convection in the liquid-phase was the most active transport mechanism. When the Kelvin effect was ignored within the soil, the gas-phase diffusion was significantly lower, leading to lower evaporation flux of water and higher volumetric liquid contents near the soil surface.