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.     

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.     

Closed-form solutions for unsaturated flow under variable flux boundary conditions

It is shown that from any solution of the linear diffusion equation, we may construct a solution of a realistic form of the Richards equation for unsaturated flow. Compared to the usual direct linearization method, our inverse approach involves a quite different sequence of transformations. This opens the possibility of exact solutions with a wider variety of continuously varying flux boundary conditions. Closed-form solutions are presented for two examples. In these, the varying water flux boundary conditions resemble (i) the passage of a peaking storm and (ii) the continuous opening of a valve preceding a steady water supply. Unlike earlier more systematic approaches to this problem, our method does not require the numerical solution of an integral equation.     

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.     

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.     

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.     

A fractal Richards’ equation to capture the non-Boltzmann scaling of water transport in unsaturated media

The traditional Richards’ equation implies that the wetting front in unsaturated soil follows Boltzmann scaling, with travel distance growing as the square root of time. This study proposes a fractal Richards’ equation (FRE), replacing the integer-order time derivative of water content by a fractal derivative, using a power law ruler in time. FRE solutions exhibit anomalous non-Boltzmann scaling, attributed to the fractal nature of heterogeneous media. Several applications are presented, fitting the FRE to water content curves from previous literature.     

Evaluation of the Ross fast solution of Richards’ equation in unfavourable conditions for standard finite element methods

Ross [Ross PJ. Modeling soil water and solute transport – fast, simplified numerical solutions. Agron J 2003;95:1352–61] developed a fast, simplified method for solving Richards’ equation. This non-iterative 1D approach, using Brooks and Corey [Brooks RH, Corey AT. Hydraulic properties of porous media. Hydrol. papers, Colorado St. Univ., Fort Collins; 1964] hydraulic functions, allows a significant reduction in computing time while maintaining the accuracy of the results. The first aim of this work is to confirm these results in a more extensive set of problems, including those that would lead to serious numerical difficulties for the standard numerical method. The second aim is to validate a generalisation of the Ross method to other mathematical representations of hydraulic functions.     

On the implementation of Prabhu’s exact solution of the stochastic reservoir equation

After 50 years of Prabhu’s paper on the exact solution of the stochastic reservoir equation for the important class of gamma inflow distributions with an integral shape parameter, a detailed implementation of the exact solution is still lacking, despite its potential usefulness from both theoretical and practical points of view. This paper explores some properties of Prabhu’s exact solution and investigates the numerical difficulties associated with its implementation. The solution is also extended to derive the distributions of deficit, spillage, yield, and actual release from the reservoir. Explicit analytical solutions for three relatively simple cases are given in detail as examples and comparisons with approximate numerical solutions are made, which reveal some shortcomings of approximate methods. The implementation of the solution in the general case reveals some numerical problems associated with large values of the shape parameter of the inflow distribution and large ratios of reservoir size to draft, mainly due to accumulation of round-off errors. A Matlab program has been developed to calculate emptying and filling probabilities over a wide range of reservoir parameters using extended precision. Comparison of Prabhu’s solution with the numerical solution of the reservoir integral equation highlights possible problems with the numerical solution, which may produce inaccurate or even invalid results for large reservoirs, large drift, and large skewness of the inflow distribution. A comparison between gamma and lognormal distributions as models of skew revealed that as the reservoir size, drift, and skewness increase, the probability of emptying of the reservoir becomes smaller for the case of gamma inflow than in the case of lognormal flow having the same skewness coefficient.     

A new saturated/unsaturated model for stormwater infiltration systems

Infiltration systems are widely used as an effective urban stormwater control measure. Most design methods and models roughly approximate the complex physical flow processes in these systems using empirical equations and fixed infiltration rates to calculate emptying times from full. Sophisticated variably saturated flow models are available, but rarely applied owing to their complexity. This paper describes the development and testing of an integrated one‐dimensional model of flow through the porous storage of a typical infiltration system and surrounding soils. The model accounts for the depth in the storage, surrounding soil moisture conditions and the interaction between the storage and surrounding soil. It is a front‐tracking model that innovatively combines a soil‐moisture‐based solution of Richard's equation for unsaturated flow with piston flow through a saturated zone as well as a reservoir equation for flow through a porous storage. This allows the use of a simple non‐iterative numerical solution that can handle ponded infiltration into dry soils. The model is more rigorous than approximate stormwater infiltration system models and could therefore be valuable in everyday practice. A range of test cases commonly used to test soil water flow models for infiltration in unsaturated conditions, drainage from saturation and infiltration under ponded conditions were used to test the model along with an experiment with variable depth in a porous storage over saturated conditions. Results show that the model produces a good fit to the observed data, analytical solutions and Hydrus. Copyright © 2008 John Wiley & Sons, Ltd.     

Influence of heterogeneous air entry pressure on large scale unsaturated flow in porous media

The paper presents numerical simulations of water infiltration in unsaturated porous media containing coarse-textured inclusions embedded in fine-textured background material. The calculations are performed using the two-phase model for water and air flow and a simplified model known as the Richards equation. It is shown that the Richards equation cannot correctly describe flow in the presence of heterogeneities. However, its performance can be improved by introducing appropriately defined effective capillary and permeability functions, representing largescale behaviour of the heterogeneous medium.     

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 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.     

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.     

Modeling of two-phase flow in membranes and porous media in microgravity as applied to plant irrigation in space

In traditional applications in soil physics it is convention to scale porous media properties, such as hydraulic conductivity, soil water diffusivity, and capillary head, with the gravitational acceleration. In addition, the Richards equation for water flux in partially saturated porous media also contains a gravity term. With the plans to develop plant habitats in space, such as in the International Space Station, it becomes necessary to evaluate these properties and this equation under conditions of microgravitational acceleration. This article develops models for microgravity steady state two-phase flow, as found in irrigation systems, that addresses critical design issues. Conventional dimensionless groups in two-phase mathematical models are scaled with gravity, which must be assigned a value of zero for microgravity modeling. The use of these conventional solutions in microgravity, therefore, is not possible. This article therefore introduces new dimensionless groups for two-phase models. The microgravity models introduced here determined that in addition to porous media properties, important design factors for microgravity systems include applied water potential and the ratio of inner to outer radii for cylindrical and spherical porous media systems.     

Water and solute movement in a coarse-textured water-repellent field soil

Unstable water flow in water-repellent unsaturated soils can significantly affect the processes of infiltration and soil water redistribution. A field experiment was carried out to study the effect of water-repellency on water and bromide movement in a coarse-textured soil in the southwestern part of The Netherlands. The field data were analyzed using a relatively simple numerical model based on the standard Richards' equation for unsaturated water flow and the Fickian-based convection-dispersion equation for solute transport. Water-repellency was accounted for by multiplying the water content and the unsaturated hydraulic conductivity of the soil with F, a factor equal to the volumetric fraction of soil occupied by preferential flow paths resulting from the unstable flow process. The good comparison of simulated and measured bromide concentrations suggests that the model provides a viable method for simulating unstable water flow in water-repellent soils.     

北京大学陆面过程模式PKULM(Peking University Land Model)介绍及检验

陆面过程模式是气候模式和天气模式的核心组成部分之一.在土壤—植被—大气耦合模式(Soil-PlantAtmosphere Model,SPAM)的基础上,发展了新一代北京大学陆面过程模式PKULM(Peking University Land Model).本文首先介绍了PKULM的辐射传输、湍流输送、光合作用、土壤水热输送等过程的参数化方案;采用隐式迭代计算框架,发展并应用了一个快速的线性方程组求解算法,提高了模式计算稳定性;提出并使用了二分搜索算法计算气孔阻抗,避免了CLM(Community Land Model)等使用的迭代方法在干旱区不稳定的情况,提高了模式的适用性;采用水势为基础的土壤水分扩散方程,使模式能够模拟土壤饱和区的水分输送过程,为进一步与水文过程模式耦合奠定了基础;还发展了一个地表积水与径流过程的机理模型,提高了模式对地表水分平衡过程的模拟能力;最后,使用"中国西北干旱区陆—气相互作用观测试验"平凉站的资料对模式进行了检验并与NOAH(National Center for Environmental Prediction,Oregon State University,Air Force,and Hydrology Lab model)陆面过程模式的模拟结果进行了比较,结果表明PKULM能够较好地模拟西北半干旱区农田下垫面地气交换过程.     

Remarks upon numerical solutions of infiltration into a soil profile

Abstract

A methodology of time-step estimation for numerically solving the Richards equation is discussed. Its importance in simulating water movement in unsaturated—saturated soils is shown for infiltration into a soil profile by applying various time-step estimations and boundary conditions for different soils. In order to test the results of the computations, infiltration theory was applied. According to infiltration theory, the pressure head in the initially unsaturated part will not take positive values as long as the moisture front has not reached the phreatic level, or, in the case of a profile with a free-draining lower boundary, it is not saturated at the base. In other cases, the appearance of positive values of the pressure head produces incorrect values for the inflow rate q.     

On the primary variable switching technique for simulating unsaturated–saturated flows

Primary variable switching appears as a promising numerical technique for variably saturated flows. While the standard pressure-based form of the Richards equation can suffer from poor mass balance accuracy, the mixed form with its improved conservative properties can possess convergence difficulties for dry initial conditions. On the other hand, variable switching can overcome most of the stated numerical problems. The paper deals with variable switching for finite elements in two and three dimensions. The technique is incorporated in both an adaptive error-controlled predictor–corrector one-step Newton (PCOSN) iteration strategy and a target-based full Newton (TBFN) iteration scheme. Both schemes provide different behaviors with respect to accuracy and solution effort. Additionally, a simplified upstream weighting technique is used. Compared with conventional approaches the primary variable switching technique represents a fast and robust strategy for unsaturated problems with dry initial conditions. The impact of the primary variable switching technique is studied over a wide range of mostly 2D and partly difficult-to-solve problems (infiltration, drainage, perched water table, capillary barrier), where comparable results are available. It is shown that the TBFN iteration is an effective but error-prone procedure. TBFN sacrifices temporal accuracy in favor of accelerated convergence if aggressive time step sizes are chosen.