Flow in unsaturated random porous media,nonlinear numerical analysis and comparison to analytical stochastic models

This work presents a rigorous numerical validation of analytical stochastic models of steady state unsaturated flow in heterogeneous porous media. It also provides a crucial link between stochastic theory based on simplifying assumptions and empirical field and simulation evidence of variably saturated flow in actual or realistic hypothetical heterogeneous porous media. Statistical properties of unsaturated hydraulic conductivity, soil water tension, and soil water flux in heterogeneous soils are investigated through high resolution Monte Carlo simulations of a wide range of steady state flow problems in a quasi-unbounded domain. In agreement with assumptions in analytical stochastic models of unsaturated flow, hydraulic conductivity and soil water tension are found to be lognormally and normally distributed, respectively. In contrast, simulations indicate that in moderate to strong variable conductivity fields, longitudinal flux is highly skewed. Transverse flux distributions are leptokurtic. the moments of the probability distributions obtained from Monte Carlo simulations are compared to modified first-order analytical models. Under moderate to strong heterogeneous soil flux conditions (σ²_y≥1), analytical solutions overestimate variability in soil water tension by up to 40% as soil heterogeneity increases, and underestimate variability of both flux components by up to a factor 5. Theoretically predicted model (cross-)covariance agree well with the numerical sample (cross-)covarianaces. Statistical moments are shown to be consistent with observed physical characteristics of unsaturated flow in heterogeneous soils.©1998 Elsevier Science Limited. All rights reserved     

A numerical model for simulating Hortonian overland flow on tropical hillslopes with vegetation elements

This paper describes a two‐dimensional hydrodynamic model that characterizes surface runoff process resulting from a varying rainfall intensity event, on an infiltrating soil surface. The soil surface has spatially varied soil physical, hydraulic and microtopographic characteristics. Infiltration process is modelled with the Philip two‐term equation and the time before ponding approximated with the time compression algorithm. Vegetation is modelled as a dynamic component with the modified Gash model. The equation is solved with a modified second order Leapfrog explicit finite difference scheme with centred time and space derivatives. The model was validated with standard analytical solutions. Evaluation with results from field campaigns in the Volta Basin of West Africa during the 2002 rainfall season indicates good agreement, with r² values ranging from 0·89 to 0·96. The developed method will be useful in studying the dynamics of surface runoff generation under complex microtopographic conditions, spatially varying soil hydraulic characteristics and temporally dynamic rainfall intensity, as found in many tropical catchments. Copyright © 2007 John Wiley & Sons, Ltd.     

On estimating the hydraulic properties of soil,Part 2. A new empirical equation for estimating hydraulic conductivity for sands

A new empirical equation to estimate hydraulic conductivity is proposed, based on a large set of measured data for hydraulic properties of soil. The equation is simpler and more accurate than the series-parallel model. Under conditions of insufficient data, the new equation provides a good estimation of hydraulic conductivity for sands. For the same class of soils, another empirical equation is proposed to estimate the power N in the Averjanov-Irmay function.     

Vulnerability of three spatially variable soil groups for solute leaching

Solute leaching in unsaturated soil is influenced by the variability in hydraulic functions (water retention and conductivity) that govern the flow process. Variability in measured soil hydraulic functions of a coarse-, medium- and fine-textured soil group was quantified with the scaling theory of similar media. Solute leaching in these soils was calculated with Monte Carlo simulation assuming, successively, hydraulic conductivity, K, volumetric water content, 0, and pressure head, h, to be constant. In addition to variability in hydraulic functions, variability in the solute retardation factor was also taken into account. To examine this effect five solutes were considered: a conservative solute (chloride), a non-retarded solute subject to decay (nitrate), a retarded solute that does not decay (cadmium) and two organic solutes which are retarded but have different sorption and decay parameters (the pesticide atrazine and a chlorinated hydrocarbon). The numerical results obtained with Monte Carlo simulation were in a number of instances verified with analytical solutions. The three soil groups distinguished showed considerable differences in vulnerability for leaching of the five solutes, emphasizing the importance of the effect of variability in soil hydraulic functions when studying solute leaching. Numerical and analytical results showed good agreement. Therefore, in relatively simple situations analytical solutions are attractive. However, in complicated situations, analytical solutions are cumbersome and numerical solutions are the only realistic alternative.     

Evaluation of transport and storage properties in the soil and groundwater zone from induced polarization measurements1

Spectral induced polarization as well as complex electrical measurements are used to estimate, on a non-invasive basis, hydraulic permeability in aquifers. Basic laboratory measurements on a variety of shaly sands, silts and clays showed that the main feature of their conductivity spectra in the frequency range from 10^-3 to 10³ Hertz is a nearly constant phase angle. Thus, a constant-phase-angle model of electrical conductivity is applied to interpret quantitatively surface and borehole spectral induced polarization measurements. The model allows for the calculation of two independent electrical parameters from only one frequency scan and a simple separation of electrical volume and interface effects. The proposed interpretation algorithm yields the true formation factor, the cation exchange capacity and the surface-area-to-porosity ratio, which corresponds to the inverse hydraulic radius. Using a Kozeny–Carman-like equation, the estimation of hydraulic permeability is possible.     

Effect of mineral reactions on the hydraulic properties of unsaturated soils: Model development and application

The selective radius shift model was used to relate changes in mineral volume due to precipitation/dissolution reactions to changes in hydraulic properties affecting flow in porous media. The model accounts for (i) precipitation/dissolution taking place only in the water-filled part of the pore space and further that (ii) the amount of mineral precipitation/dissolution within a pore depends on the local pore volume. The pore bundle concept was used to connect pore-scale changes to macroscopic soil hydraulic properties. Precipitation/dissolution induces changes in the pore radii of water-filled pores and, consequently, in the effective porosity. In a time step of the numerical model, mineral reactions lead to a discontinuous pore-size distribution because only the water-filled pores are affected. The pore-size distribution is converted back to a soil moisture characteristic function to which a new water retention curve is fitted under physically plausible constraints. The model equations were derived for the commonly used van Genuchten/Mualem hydraulic properties. Together with a mixed-form solution of Richards’ equation for aqueous phase flow, the model was implemented into the geochemical modelling framework PHREEQC, thereby making available PHREEQC’s comprehensive geochemical reactions. Example applications include kinetic halite dissolution and calcite precipitation as a consequence of cation exchange. These applications showed marked changes in the soil’s hydraulic properties due to mineral precipitation/dissolution and the dependency of these changes on water contents. The simulations also revealed the strong influence of the degree of saturation on the development of the saturated hydraulic conductivity through its quadratic dependency on the van Genuchten parameter α. Furthermore, it was shown that the unsaturated hydraulic conductivity at fixed reduced water content can even increase during precipitation due to changes in the pore-size distribution.     

Dissolution of nonaqueous phase liquid pools in anisotropic aquifers

A two-dimensional numerical transport model is developed to determine the effect of aquifer anisotropy and heterogeneity on mass transfer from a dense nonaqueous phase liquid (DNAPL) pool. The appropriate steady state groundwater flow equation is solved implicitly whereas the equation describing the transport of a sorbing contaminant in a confined aquifer is solved by the alternating direction implicit method. Statistical anisotropy in the aquifer is introduced by two-dimensional, random log-normal hydraulic conductivity field realizations with different directional correlation lengths. Model simulations indicate that DNAPL pool dissolution is enhanced by increasing the mean log-transformed hydraulic conductivity, groundwater flow velocity, and/or anisotropy ratio. The variance of the log-transformed hydraulic conductivity distribution is shown to be inversely proportional to the average mass transfer coefficient.     

LE TRITIUM DANS LES PRÉCIPITATIONS DE L'ASIE DU SUD-EST ET LE PHÉNOMÉNE DE LA MOUSSON EN 1966 ET EN 1967

Abstract

The problem of non-steady flow of water in a soil-plant system can be described by adding a sink term to the continuity equation for soil water flow. In this paper the sink term is defined in two different ways. Firstly it is considered to be dependent on the hydraulic conductivity of the soil, on the difference in pressure head between the soil and the root-soil interface and some root effectiveness function. Secondly the sink is taken to be a prescribed function of the soil water content. The partial differential equation applying to the first problem is solved by both a finite difference (FD 1) and a finite element (FE 1) technique, that applying to the second problem by a finite difference approach (FD 2). The purpose of this paper is to verify the numerical models against field measurements, to compare the results obtained by the three numerical methods and to show how the finite element method can be applied to complex but realistic two-dimensional flow situations. Two examples are given. The first concerns one-dimensional flow and it compares numerical results with those obtained experimentally in the field from water balance studies on red cabbage (Brassica oleracea L. ‘Rode Herfst’) grown on a clay soil in the presence of a water table. The second example describes two-dimensional flow in a complex field situation in the Netherlands where flow takes place under cropped field conditions through five anisotropic layers. Water is supplied to the system by infiltration from two unlined ditches and is withdrawn from the system by evapotranspiration and by leakage to an underlying pumped aquifer.     

A dripper-TDR method for in situ determination of hydraulic conductivity and chemical transport properties of surface soils

Field determined hydraulic and chemical transport properties can be useful for the protection of groundwater resources from land-applied chemicals. Most field methods to determine flow and transport parameters are either time or energy consuming and/or they provide a single measurement for a given time period. In this study, we present a dripper-TDR field method that allows measurement of hydraulic conductivity and chemical transport parameters at multiple field locations within a short time period. Specifically, the dripper-TDR determines saturated hydraulic conductivity (K_s), macroscopic capillary length (λ_c), immobile water fraction (θ_im/θ), mass exchange coefficient (α) and dispersion coefficient (D_m). Multiple dripper lines were positioned over five crop rows in a field. Background and step solutions were applied through drippers to determine surface hydraulic conductivity parameters at 44 locations and surface transport properties at 38 locations. The hydraulic conductivity parameters (K_s, λ_c) were determined by application of three discharge rates from the drippers and measurements of the resultant steady-state flux densities at the soil surface beneath each dripper. Time domain reflectometry (TDR) was used to measure the bulk electrical conductivity of the soil during steady infiltration of a salt solution. Breakthrough curves (BTCs) for all sites were determined from the TDR measurements. The K_s and λ_c values were found to be lognormally distributed with average values of 31.4 cm h⁻¹ and 6.0 cm, respectively. BTC analysis produced chemical properties, θ_im/θ, α, and D_m with average values of 0.23, 0.0036 h⁻¹, and 1220 cm² h⁻¹, respectively. The estimated values of the flow and transport parameters were found to be within the ranges of values reported by previous studies conducted at nearby field locations. The dripper TDR method is a rapid and useful technique for in situ measurements of hydraulic conductivity and solute transport properties. The measurements reported in this study give clear evidence to the occurrence of non-equilibrium water and chemical movement in surface soil. The method allows for quantification of non-equilibrium model parameters and preferential flow. Quantifying the parameters is a necessary step toward determining the influences of surface properties on infiltration, runoff, and vadose zone transport.     

Assessing the impact of the hydraulic properties of a crusted soil on overland flow modelling at the field scale

Soil surface crusts are widely reported to favour Hortonian runoff, but are not explicitly represented in most rainfall‐runoff models. The aim of this paper is to assess the impact of soil surface crusts on infiltration and runoff modelling at two spatial scales, i.e. the local scale and the plot scale. At the local scale, two separate single ring infiltration experiments are undertaken. The first is performed on the undisturbed soil, whereas the second is done after removal of the soil surface crust. The HYDRUS 2D two‐dimensional vertical infiltration model is then used in an inverse modelling approach, first to estimate the soil hydraulic properties of the crust and the subsoil, and then the effective hydraulic properties of the soil represented as a single uniform layer. The results show that the crust hydraulic conductivity is 10 times lower than that of the subsoil, thus illustrating the limiting role the crust has on infiltration. Moving up to the plot scale, a rainfall‐runoff model coupling the Richards equation to a transfer function is used to simulate Hortonian overland flow hydrographs. The previously calculated hydraulic properties are used, and a comparison is undertaken between a single‐layer and a double‐layer representation of the crusted soil. The results of the rainfall‐runoff model show that the soil hydraulic properties calculated at the local scale give acceptable results when used to model runoff at the plot scale directly, without any numerical calibration. Also, at the plot scale, no clear improvement of the results can be seen when using a double‐layer representation of the soil in comparison with a single homogeneous layer. This is due to the hydrological characteristics of Hortonian runoff, which is triggered by a rainfall intensity exceeding the saturated hydraulic conductivity of the soil surface. Consequently, the rainfall‐runoff model is more sensitive to rainfall than to the subsoil's hydrodynamic properties. Therefore, the use of a double‐layer soil model to represent runoff on a crusted soil does not seem necessary, as the increase of precision in the soil discretization is not justified by a better performance of the model. Copyright © 2005 John Wiley & Sons, Ltd.     

Bulk density effects on soil hydrologic and thermal characteristics: A numerical investigation

Soil bulk density (ρ_b) is commonly treated as static in studies of land surface dynamics. Magnitudes of errors associated with this assumption are largely unknown. Our objectives were to (a) quantify ρ_b effects on soil hydrologic and thermal properties and (b) evaluate effects of ρ_b on surface energy balance and heat and water transfer. We evaluated 6 soil properties, volumetric heat capacity, thermal conductivity, soil thermal diffusivity, water retention characteristics, hydraulic conductivity, and vapour diffusivity, over a range of ρ_b, using a combination of 6 models. Thermal conductivity, water retention, hydraulic conductivity, and vapour diffusivity were most sensitive to ρ_b, each changing by fractions greater than the associated fractional changes in ρ_b. A 10% change in ρ_b led to 10–11% change in thermal conductivity, 6–11% change in saturated and residual water content, 49–54% change in saturated hydraulic conductivity, and 80% change in vapour diffusivity. Subsequently, 3 field seasons were simulated with a numerical model (HYDRUS‐1D) for a range of ρ_b values. When ρ_b increased 25% (from 1.2 to 1.5 Mg m^?3), soil temperature variation decreased by 2.1 °C in shallow layers and increased by 1 °C in subsurface layers. Surface water content differed by 0.02 m³ m^?3 for various ρ_b values during drying events but differences mostly disappeared in the subsurface. Matric potential varied by >100 m of water. Surface energy balance showed clear trends with ρ_b. Latent heat flux decreased 6%, sensible heat flux increased 9%, and magnitude of ground heat flux varied by 18% (with a 25% ρ_b increase). Transient ρ_b impacted surface conditions and fluxes, and clearly, it warrants consideration in field and modelling investigations.     

A simplified stochastic model for infiltration into a heterogeneous soil forced by random precipitation

A simplified stochastic infiltration model is presented, aimed at representing the rainfall-runoff transformation in the presence of heterogeneity in the soil and with precipitation as a random variable with complex temporal evolution. Such a model is based on a simple water mass balance of a surface soil layer, considered as a non-linear reservoir. The explicit inclusion of spatial heterogeneity allows the model to be used in sub-grid parametrizations at a variety of scales, from the distributed modelling of the hydrological response of small watersheds to the representation of surface mass fluxes in General Circulation Models. An approximate solution procedure is developed, which allows the estimation of statistical moments of the soil effective saturation and runoff inside discrete time steps where the hydraulic saturated conductivity and the rainfall intensity are taken as random variables with known probability density functions. As a first test of the proposed model, two different simulations, relative to two soils with different hydraulic conductivity distributions, are presented and discussed. A year long record of hourly averaged rainfall intensities, as measured by a tipping bucket gauge in Central Italy, is taken as the main input. The main finding is that the non-linear nature of the soil filter is such that, for random precipitation intensity, the coefficient of variation of the runoff is always higher than that of precipitation. Such a non-linear variability enhancement, due mainly to the threshold character of the soil mass balance equation, tends to be slightly dampened by the variability of the hydraulic saturated conductivity.     

A Douglas-Jones Predictor-Corrector Program for Simulating One-Dimensional Unsaturated Flow in Soil

Abstract. A fully documented program to represent one-dimensional unsaturated flow in soil is described. The program is based on a Douglas-Jones finite-difference implicit method to solve the Richards equation. An implicit linearization scheme is used to estimate the hydraulic conductivity and specific moisture capacity functions. Predicted moisture content profiles compared with two Galerkin finite-element solutions and field observations on a Panoche clay loam soil show very good agreement.     

Effective hydraulic conductivity of nonstationary aquifers

Assuming that the ln hydraulic conductivity in an aquifer is mathematically approximated by a spatial deterministic surface, or trend, plus a stationary random noise, we treat the problem of finding what the effective hydraulic conductivity of that aquifer is. This problem is tackled by spectral methods applied to a type of diffusion equation of groundwater flow, together with suitable coordinate transformations. Analytical (exact) solutions in terms of elementary functions are presented for one- and three-dimensional finite and infinite domains. Stability criteria are obtained for the solutions, in terms of a critical parameter, that turns out to involve the product of correlation scale and trend gradient. For the case of finite and symmetrical domains, additional provisions to insure the stability of numerical calculations of effective hydraulic conductivity are provided. Effective hydraulic conductivity is an important property, with potential applications in the calibrations of groundwater and transport numerical models.     

Effective hydraulic conductivity of nonstationary aquifers

Assuming that the ln hydraulic conductivity in an aquifer is mathematically approximated by a spatial deterministic surface, or trend, plus a stationary random noise, we treat the problem of finding what the effective hydraulic conductivity of that aquifer is. This problem is tackled by spectral methods applied to a type of diffusion equation of groundwater flow, together with suitable coordinate transformations. Analytical (exact) solutions in terms of elementary functions are presented for one- and three-dimensional finite and infinite domains. Stability criteria are obtained for the solutions, in terms of a critical parameter, that turns out to involve the product of correlation scale and trend gradient. For the case of finite and symmetrical domains, additional provisions to insure the stability of numerical calculations of effective hydraulic conductivity are provided. Effective hydraulic conductivity is an important property, with potential applications in the calibrations of groundwater and transport numerical models.     

An analytical solution for non-Darcian flow in a confined aquifer using the power law function

We have developed a new method to analyze the power law based non-Darcian flow toward a well in a confined aquifer with and without wellbore storage. This method is based on a combination of the linearization approximation of the non-Darcian flow equation and the Laplace transform. Analytical solutions of steady-state and late time drawdowns are obtained. Semi-analytical solutions of the drawdowns at any distance and time are computed by using the Stehfest numerical inverse Laplace transform. The results of this study agree perfectly with previous Theis solution for an infinitesimal well and with the Papadopulos and Cooper’s solution for a finite-diameter well under the special case of Darcian flow. The Boltzmann transform, which is commonly employed for solving non-Darcian flow problems before, is problematic for studying radial non-Darcian flow. Comparison of drawdowns obtained by our proposed method and the Boltzmann transform method suggests that the Boltzmann transform method differs from the linearization method at early and moderate times, and it yields similar results as the linearization method at late times. If the power index n and the quasi hydraulic conductivity k get larger, drawdowns at late times will become less, regardless of the wellbore storage. When n is larger, flow approaches steady state earlier. The drawdown at steady state is approximately proportional to r¹⁻ⁿ, where r is the radial distance from the pumping well. The late time drawdown is a superposition of the steady-state solution and a negative time-dependent term that is proportional to t^{(1−n)/(3−n)}, where t is the time.     

The exponential decline in saturated hydraulic conductivity with depth: a novel method for exploring its effect on water flow paths and transit time distribution

The strong vertical gradient in soil and subsoil saturated hydraulic conductivity is characteristic feature of the hydrology of catchments. Despite the potential importance of these strong gradients, they have proven difficult to model using robust physically based schemes. This has hampered the testing of hypotheses about the implications of such vertical gradients for subsurface flow paths, residence times and transit time distribution. Here we present a general semi‐analytical solution for the simulation of 2D steady‐state saturated‐unsaturated flow in hillslopes with saturated hydraulic conductivity that declines exponentially with depth. The grid‐free solution satisfies mass balance exactly over the entire saturated and unsaturated zones. The new method provides continuous solutions for head, flow and velocity in both saturated and unsaturated zones without any interpolation process as is common in discrete numerical schemes. This solution efficiently generates flow pathlines and transit time distributions in hillslopes with the assumption of depth‐varying saturated hydraulic conductivity. The model outputs reveal the pronounced effect that changing the strength of the exponential decline in saturated hydraulic conductivity has on the flow pathlines, residence time and transit time distribution. This new steady‐state model may be useful to others for posing hypotheses about how different depth functions for hydraulic conductivity influence catchment hydrological response. Copyright © 2016 John Wiley & Sons, Ltd.     

Comparison of finite difference and finite element solutions to the variably saturated flow equation

Numerical solutions to the equation governing variably saturated flow are usually obtained using either the finite difference (FD) method or the finite element (FE) method. A detailed comparison of these methods shows that the main difference between them is in how the numerical schemes spatially average the variation of material properties. Further differences are also observed in the way that flux boundaries are represented in FE and FD methods. A modified finite element (MFE) algorithm is used to explore the significance of these differences. The MFE algorithm enables a direct comparison with a typical FD solution scheme, and explicitly demonstrates the differences between FE and FD methods. The MFE algorithm provides an improved approximation to the partial differential equation over the usual FD approach while being computationally simpler to implement than the standard FE solution. One of the main limitations of the MFE algorithm is that the algorithm was developed by imposing several restrictions upon the more general FE solution; however, the MFE is shown to be preferable over the usual FE and FD solutions for some of the test problems considered in this study. The comparison results show that the FE (or MFE) solution can avoid the erroneous results encountered in the FD solution for coarsely discretized problems. The improvement in the FE solution is attributed to the broader hydraulic conductivity averaging and differences in the representation of flux type boundaries.