Richards Equation Research Articles

How pore structure non‐stationarity compromises flow properties representativity (REV) for soil samples: Pore‐scale modelling and stationarity analysis

AbstractClassic soil physics relies heavily on the concept of representative elementary volume (REV), which is necessary to perform upscaling from the studied soil samples and parameterize continuum scale hydrological models (e.g., based on Richards or Darcy equations). In this paper, we explore the boundaries of the classic REV concept and conventional representativity studies that claim REV for any given physical property, and that its values converge to a steady value with increasing sample volume. We chose two undisturbed soil samples of standard size from Ah and B horizons, subcropped two subvolumes within each of them and performed pore‐scale flow simulations using binarized X‐ray microtomography scans as input data. The volume of the simulation domains was 9003 voxels, with a physical volume within two orders of magnitude of the whole soil core. Based on the 3D pore geometry images and the resulting flow velocity and pressure fields, we performed REV analysis for the saturated hydraulic conductivity and porosity. Although in general subvolumes showed classical REV behaviour (convergence of porosity and saturated hydraulic conductance (Ksat)to plateau‐like behaviour), their flow properties converged to different REV values. We also evaluated the stationarity of pore structures by computing the directional correlation functions to explain the observed non‐unique behaviour. We concluded that neither of the studied samples can be considered to be representative due to their structural non‐stationarity. We extensively discussed the implications of these results for the upscaling and parameterization of continuum‐scale flow models. We argued that the plateau in Ksat (or any other physical property) is a necessary, but insufficient, condition for the REV, which requires (at least) pore structure stationarity as an additional criterion. REV as a concept is much broader than previously anticipated and the possibility of establishing REVs in structured soils will require an analysis of tensorial flow properties with correct boundary conditions, multiscale soil structure imaging with stationarity analysis, and pore‐scale simulations on fused multiscale images.Highlights Ksat analysis within subvolumes showed classical REV behaviour but resulted in non‐unique REVs for the soil samples studied Based on the stationarity analysis we argued that the conventional soil core samples studied were not REVs Plateau‐like behaviour in physical properties is a necessary but insufficient condition for REV Soil structure information is necessary to properly establish REV or upscale flow properties

Darcy, Forchheimer, Brinkman and Richards: classical hydromechanical equations and their significance in the light of the TPM

In hydromechanical applications, Darcy, Brinkman, Forchheimer and Richards equations play a central role when porous media flow under saturated and unsaturated conditions has to be investigated. While Darcy, Brinkman, Forchheimer and Richards found their equations mainly on the basis of flow observations in field and laboratory experiments, the modern Theory of Porous Media allows for a scientific view at these equations on the basis of precise continuum mechanical and thermodynamical investigations. The present article aims at commenting the classical equations and at deriving their counterparts by the use of the thermodynamical consistent Theory of Porous Media. This procedure will prove that the classical equations are valid under certain restrictions and that extended equations exist valid for arbitrary cases in their field.

Hole irrigation process simulation using a soil water dynamical model with parameter inversion method

Soil water dynamics are an important part of agricultural irrigation evaluation standards; therefore, supplying water efficiently has become the focus of agricultural water-saving research, and different soil water movement models have been developed for various water supply models. In this work, a parameter inversion model is coupled with the soil water movement equation (the Richards equation) based on the main variable of the head (h)-based Richards equation to investigate the dynamics of soil moisture under hole irrigation. A mathematical model for hole irrigation with two hole diameters (5 cm or 7 cm) and two hole numbers (one or two) was constructed after verifying the model with field-scale soil water content data from different soil profile depths. The Richards equation based on the Gauss–Newton parameter inversion method was found to accurately simulate the soil moisture dynamics under hole irrigation (R2>0.9). The difference in the water content of the soil profile due to the difference in hole diameter did not exceed 2.7%, and the soil water content under single-hole irrigation was significantly higher than that under double-hole irrigation due to lower evaporation. The infiltration distance of the wetting front under single-hole irrigation, which was mainly distributed on the hole side of the planting ridge, was greater than that under double-hole irrigation at 378 ml irrigation water per plant, and the soil moisture was unevenly distributed at a certain depth. However, the wetting fronts under double-hole irrigation intersected and began to infiltrate the soil evenly within a short time (only 6% of the entire irrigation cycle). The results for wetting front movement under the different hole settings show that double-hole irrigation does not stimulate crop root growth on both sides of the plant in the long term; however, this growth pattern can be stimulated by alternating single-hole irrigation.

Modeling Soil Water Redistribution under Gravity Irrigation with the Richards Equation

Soil water movement is important in fields such as soil mechanics, irrigation, drainage, hydrology, and agriculture. The Richards equation describes the flow of water in an unsaturated porous medium, and analytical solutions exist only for simplified cases. However, numerous practical situations require a numerical solution (1D, 2D, or 3D) depending on the complexity of the studied problem. In this paper, numerical solution of the equation describing water infiltration into soil using the finite difference method is studied. The finite difference solution is made via iterative schemes of local balance, including explicit, implicit, and intermediate methods; as a special case, the Laasonen method is shown. The found solution is applied to water transfer problems during and after gravity irrigation to observe phenomena of infiltration, evaporation, transpiration, and percolation.

Inverse modeling towards parameter estimation of the nonlinear soil hydraulic functions using developed multistep outflow procedure

Modeling and monitoring vadose zone processes have a key role to effectively identify and solve important emerging issues in hydrology and related disciplines. However, major gaps and challenges remain concerning the further conceptualization of the zone complexities and parameterization of relevant models. Consequently, the main focus of this investigation was to establish a modern developed version of the multistep outflow experiments for the vadose zone characterization and modeling. Besides, soil hydraulic properties (K-θ-h) were determined by the Mualem/Van Genuchten model based on the inverse solution of the Richards equation for the four loamy textured soils. The adjustable parameters were estimated with the least-square approach through the implementation of the generalized reduced gradient (GRG) algorithm. Then, the inverse manner was evaluated by comparing optimized hydraulic functions with those measured independently using an equilibrium condition. The results revealed the optimized attributes complied well with those determined using equilibrium analysis and steady-state situations. The reliable estimation of the parameters to inversely optimize the K-θ-h functions was possible only by measuring the cumulative outflow as a function of time. The simultaneous optimization of the unsaturated hydraulic functions carried out from a single transient experiment. This protocol appears well suited for vadose zone simulation. Furthermore, along with the tensiometer elimination and no independent measurements of Ks and θs, the developed version is generally more cost-effective, faster and more convenient than the conventional method. Finally, we conclude the combination of proposed procedure with GRG algorithm can provide a reasonable procedure to inversely estimate the hydraulic parameters.

Simulation of Preferential Flow in Snow With a 2‐D Non‐Equilibrium Richards Model and Evaluation Against Laboratory Data

AbstractRecent studies of water flow through dry porous media have shown progress in simulating preferential flow propagation. However, current methods applied to snowpacks have neglected the dynamic nature of the capillary pressure, such as conditions for capillary pressure overshoot, resulting in a rather limited representation of the water flow patterns through snowpacks observed in laboratory and field experiments. Indeed, previous snowmelt models using a water entry pressure to simulate preferential flow paths do not work for natural snowpack conditions where snow densities are less than 380 kg m−3. Because preferential flow in snowpacks greatly alters the flow velocity and the timing of delivery of meltwater to the base of a snowpack early in the melt season, a better understanding of this process would aid hydrological predictions. This study presents a 2‐D water flow through snow model that solves the non‐equilibrium Richards equation. This model, coupled with random perturbations of snow properties, can represent realistic preferential flow patterns. Using 1‐D laboratory data, two model parameters were linked to snow properties and model boundary conditions. Parameterizations of these model parameters were evaluated against 2‐D snowpack observations from a laboratory experiment, and the resulting model sensitivity to varying inputs and boundary conditions was calculated. The model advances both the physical understanding of and ability to simulate water flow through snowpacks and can be used in the future to parameterize 1‐D snowmelt models to incorporate flow variations due to preferential flow path formation.

Exact Solutions of the Generalized Richards Equation with Power-Law Nonlinearities

We consider a generalized Richards equation with power-law nonlinearities modeling filtration in porous media. Conditions are derived under which the problem can be reduced to the linear heat equation or to nonlinear equations with known solutions. The families of explicit exact solutions that can be expressed via elementary functions or Lambert’s $$W $$ -function are found. Some examples illustrating the results are given.

Computational investigation on surface water distribution and permeability of porous asphalt pavement

ABSTRACT This paper used the finite element method to analyse surface water distribution characteristics on porous asphalt layer under unsaturated seepage, and further compared the critical rainfall intensity of porous asphalt pavements with different design parameters. In this study, the Richards equation was used to describe the unsaturated seepage, and the results were validated with experimental data. The results show that when the water accumulates on the pavement surface, it would firstly emerge at the lowest level of pavement along the longitudinal direction, and then spread along transverse and longitudinal directions. The water depth has no remarkable difference along the longitudinal direction while it differs a lot among different transverse positions. The distribution of surface water depth was analysed and the correlation between the surface water coverage in the transverse direction and the rainfall intensity was established with high R-square values. For pavements with different geometries, the porous layer thickness has the most remarkable effect on the critical rainfall intensity of porous asphalt pavements.

Effects of soil heterogeneity on subsurface water movement in agricultural fields: A numerical study

This study aims to numerically examine the effects of soil heterogeneity on subsurface water movement in agricultural fields under wheat (Triticum aestivum L.) and rice (Oryza sativa L.) crop covers. The 3-D Richards equation that accounts for root water uptake was coupled with a random field generation method to produce realizations of soil hydraulic properties, and the effects of soil heterogeneity were analyzed under realistic atmospheric boundary conditions. Four soil types and three different soil heterogeneity levels were considered to investigate the response of soil heterogeneity on soil water dynamics. It was observed that the effect of soil heterogeneity varies across soil types and its effect was more pronounced in silty clay loam soils. Further, an increase in soil heterogeneity level can increase or decrease the mean soil moisture, hydraulic conductivity, and flow velocity based on land cover, upper boundary condition, and atmospheric demand, and its effect was higher during irrigation periods. A maximum difference of 30% in soil moisture values was observed near the surface soil layers due to an increase in soil heterogeneity on irrigation days. The mean soil moisture content under rice crop cover was always lower than that of wheat crop cover due to higher atmospheric demand in the former. Soil layering effect on soil moisture dynamics, which is related to capillary and hydraulic barrier formation, was found to depend on the layering configuration, available soil moisture, and land cover. The numerical experiments provided useful insights into the complex soil moisture dynamics in agricultural fields.

Can structure data obtained from CT images substitute for parameters of a preferential flow model?

Modeling preferential flow in soils is still a challenge for the scientific community working on water resources. Indeed, it is an issue to determine the functional parameters of models dedicated to water flow, that are currently obtained by fitting processes, whereas their relationships with soil structure remain poorly known. Improved models are expected from a better understanding of the links between functional and structural parameters, which can be achieved thanks to recent developments in imaging methods such as X-ray Computed Tomography (CT). The paper seeks to improve a dual-porosity model, coupling matrix flow (by Richards equation) and preferential flow (by a Kinematic Dispersive Wave), by substituting some model parameters, usually obtained by inversion of experimental data, by those assessed from CT images of the soil structure. Thus, two versions of the model are compared, the “classical” and the “advanced” one including parameters determined using the 3D images of the sample structures. To compare model versions with real situations, infiltration experiments were conducted in lab on two different soils at two initial water contents. An X-ray medical scanner allowing acquisitions of large soil volumes (≈1700 cm3) with a voxel size of 400 µm was used to image the sample structures. Then, we derived two geometrical parameters from the macroporosity network: the percolating macroporosity and a characteristic dimension of this macropore network, the mean inter-macropore distance. The sensitivity analysis conducted on the classical version of the model showed that the kinematic coefficient and the dimensional parameter of the porous medium are the two main contributors to the cumulated drainage whatever the initial condition. Although experimental data are better simulated by the classical version of the model, drainage dynamics is also well simulated by the advanced version. However, differences between the model versions that are small for both soils at field capacity become significant for the dried state (mean initial matric potential of −3.5 m). This emphasizes the crucial effect of the sink-source term and in particular the complex effect of the dimensional parameter that it contains. Indeed, difficulties to simulate properly water exchange between porosity domains are encountered for both versions of the model. We conclude that empirical parameters that were up to now fitted from experiments could be deduced from geometrical indicators computed from CT images and that owning to these first results the applied methodology is promising to achieve a better understanding and modeling of preferential flow processes and to improve model predictability.

A new variable-order fractional derivative with non-singular Mittag–Leffler kernel: application to variable-order fractional version of the 2D Richard equation

In this study, to overcome the limitations of non-singular fractional derivatives in the Caputo–Fabrizio and Atangana–Baleanu senses (especially in dealing with the variable-order (VO) fractional calculus), we introduce a new non-singular VO fractional derivative with Mittag–Leffler function as its kernel. Some useful results are derived from this fractional derivative. Moreover, this fractional derivative is used for introducing the VO fractional version of the 2D Richard equation. A meshless scheme based on the thin plate spline radial basis functions (RBFs) is developed for solving this equation. The validity of the formulated method is investigated through three numerical examples.

Adsorption of Dye on a Tunisian Unsaturated Layered Soil: Physical and Numerical Modeling

The main objectives of this study is to model the transport of methylene blue (MB) in homogenous and layered soils which has been studied experimentally by fixed bed column. The effect of soil stratification has been studied through numerical investigation based on the coupled solute transport model in three-layered soil. The effect of some significant parameters such as flow rate, initial concentration of MB, thickness of each layer and the numbers of layers on the breakthrough curves have been undertaken. A finite element analysis model was employed to predict the transport of MB in soils. A two dimensional model based on Richards equation and advection-dispersion equation coupled with adsorption model has been developed and an analytical model has been used to predict the dispersivity. The results shows that the soil heterogeneity has a significant effect on methylene blue (MB) adsorption through unsaturated layered media and the effect of the parameters of the upper layer of the soil is more significant than those of the lower layer but this effect can be controlled with the choice of the thicknesses of each layer.

Coupling surface flow and subsurface flow in complex soil structures using mimetic finite differences

We explore the coupling of surface and subsurface flows on fully unstructured meshes that conform to complex soil structures. To accommodate the distorted meshes that inevitably result from explicit representation of complex soil structures, we leverage the structure of the Mimetic Finite Difference (MFD) spatial discretization scheme to couple surface and subsurface flows. The MFD method achieves second-order accuracy and maintains local mass conservation on distorted meshes. We couple the diffusion wave approximation for surface flows to the Richards equation for subsurface flow, ensuring continuity of both pressure and flux between the surface and subsurface. The MFD method is particularly convenient for this coupling because it uses face-based constraints in the subsurface system that can be expressed as face-pressure unknowns. Those unknowns are coincident with surface cell-based unknowns, thus allowing the discrete surface system to be directly substituted into the subsurface system and solved implicitly as a global system. Robust representation of the transition between wet and dry surface conditions requires upwinding of the relative permeability and is facilitated by globalization in the nonlinear solver. The approach and its implementation in the Advanced Terrestrial Simulator (ATS) are evaluated by comparison to previously published benchmarks. Using runoff from soils with patchy groundcover (duff) as an example, we show that the new method converges significantly faster in mesh convergence tests than the commonly used two-point flux approximation.

A continuum model of unstable infiltration in porous media endowed with an entropy function

We propose a thermodynamic approach to modeling unsaturated flow in porous media, where the liquid saturation is understood as the state variable. The free energy functional is designed as a symmetric expansion of the traditional capillary energy density in Richards equation, therefore removing ambiguities on the interpretation of the higher-order term in the model equation. The proposed definition renders a formulation that leads naturally to an entropy function of the system, and we show that the model describes an entropy-increasing process for an isolated system.The new formulation reproduces gravity fingering during infiltration in soil. We show that the nonlinear and singular structure of the capillary pinning function in the fourth-order term plays a fundamental role in the behavior and stability of infiltration fronts, promoting front pinning and the persistence of fingered infiltration at relatively large flux ratios.

Potential evaporation dynamics over saturated bare soil and an open water surface

Actual evaporation (Ea) can be calculated as a fraction of potential evaporation (PE), which refers to the evaporation rate if supply water is unlimited. Potential evaporation depends on the available energy and the underlying material, and different approaches to estimate potential evaporation exist nowadays. This study provides a detailed analysis of the evaporation dynamics over fully saturated, sandy soil (PEs) and an open water surface (PEw). Moreover, the performance of commonly used methods to estimate PE is assessed. At the basis of these analyses is a lysimeter experiment in the Guanzhong Basin, China, which allowed a precise measurement of PE with a very high temporal resolution. Temperature profiles in lysimeters and meteorological data were also measured during the experiment. A comparison of PEs and PEw was carried out for seven consecutive days (August 11th to 17th, 2016). Results show that PEw is smaller than PEs on a daily scale, with PEw rates being bigger than PEs at night but smaller during daytime. Furthermore, the temporal dynamics of PEw lags 4–5 h behind PEs. In accordance with the energy balance equation, PE dynamics are mainly governed by “available energy”. The PE rates calculated by Penman-Monteith (PM) and Priestly-Taylor (PT) based on these measurements were also evaluated. The measured PE is relatively well reproduced by PM and PT equations. Finally, the effect of using different approaches to estimate PE on calculating Ea was tested by an integrated hydrological model that calculates water flow in the unsaturated zone by solving the Richards equation. The relative differences were up to 17.5%.

The importance of unsaturated hydraulic conductivity measurements for green roof detention modelling

Characterising the unsaturated hydraulic conductivity of a green roof substrate is essential for accurately modelling runoff detention in response to rainfall events. In this paper, the unsaturated hydraulic conductivities for four representative green roof substrates were determined in an infiltration column using steady state and transient techniques. The conventional Durner-Mualem Hydraulic Conductivity Function (HCF) model, for which parameters were calibrated based on the measured Soil Water Release Curve (SWRC) data, was shown to provide a poor fit to the experimental data. A new three-segment HCF was, therefore, proposed to fit measured unsaturated hydraulic conductivity data. Detention tests were carried out on 100 mm and 200 mm deep substrates using four simulated storm events. The runoff and moisture content data collected during the detention tests was used to validate the HCFs using the Richards Equation. The new three-segment HCF resulted in simulated runoff and moisture content profiles that closely matched the measured data (with mean Rt2 = 0.754 for modelled runoff), in contrast to predictions made using the conventional Durner-Mualem model (with mean Rt2 = 0.409 for modelled runoff). It was also demonstrated that further simplification of the HCF to a function defined by moisture content at just two points – the saturated hydraulic conductivity and at an unsaturated hydraulic conductivity of 0.1 cm/min – provides a model that is fit-for-purpose for green roof runoff estimation (with mean Rt2 = 0.629 for modelled runoff).

Methane production from mesophilic fed-batch anaerobic digestion of empty fruit bunch of palm tree

In this work, mesophilic anaerobic digestion (AD) of empty fruit bunch (EFB) of the date palm tree was investigated in order to evaluate the application for this type of substrate, using a fed-batch process. Several parameters were monitored such as temperature, pH, and methane production throughout the experiments. The biodegradability was calculated by determining the volatile solids (VS) and total solids (TS) and the end of each added load. The process was started up with 8 gVS/L of inoculum, provided from a full-scale reactor from the molasses processing industry. The digester was loaded in fed-batch mode by increasing gradually the load from 1 to 3.5 gVS/L. The results showed that the methane yield coefficient and the organic loading rate (OLR) were 118.5 NmL/gVS and 0.27 gVS/L.d, respectively. Three kinetic models were tested the best fit (modified Gompertz models, modified Richard, and transference equation). The results show that the modified Gompertz model was the most suitable for fitting the EFB substrate and more adjusted to the experimental results.

Simultaneous analysis of slope instabilities on a small catchment-scale using coupled surface and subsurface flows

High-velocity runoff generated in hillslopes during heavy rainfall caused by typhoon increases chances of instability of the embankment slopes at the exit of the gully. Such effects of high-velocity runoff are usually neglected in conventional rainfall-induced slope failure analysis. In order to consider the effects of runoff on the slope instability, this study attempts to simulate the runoff, infiltration, seepage, and slope instabilities on a small catchment-scale simultaneously. For this purpose, this study firstly proposes a coupled model of surface flow, subsurface flow, and soil mechanics based on shallow water equations, Richards's equation, Green-Ampt infiltration capacity model, and local factor of safety (LFS) approach. Next, to make the proposed coupled model effective in the practical analysis of runoff, a diffusion wave approximation of shallow water equations is validated by numerical simulations, and then it is used to replace shallow water equations in the proposed coupled model. Finally, the proposed coupled model is verified by Abdul and Gillham system and applied to a natural slope in Hokkaido, Japan inside a small catchment with an area of 0.4675 km2. The numerical results highlight the influences of runoff from upstream on the embankment slope failure at the exit of the gully. Furthermore, the small catchment-scale slope instabilities assessment approach proposed in this study provides an effective approach for simulating heavy rainfall induced runoff and slope instabilities. The distribution map of the factor of safety (FOS) has significant implications for precisely determining the dangerous spots (instead of areas) on a small catchment-scale and accurately releasing early warning information to these dangerous spots.

Численное моделирование процесса просачивания в трещиновато-пористый грунт в условиях Крайнего Севера

Приводится математическое моделирование комплексной мультифизичной задачи, актуальной для территорий Крайнего Севера и Арктики. Актуальность данной задачи характеризируется важностью процесса просачивания при формировании и оттаивании слоя вечной мерзлоты. Современные прикладные задачи по большей части требуют учета сложных геометрий, а также большого количества разных процессов и их взаимовлияния. Мультифизичная модель состоит из уравнения Ричардса для описания процесса просачивания, модели двойной пористости для описания естественной трещиноватости грунта, задачи Стефана для описания температурного режима грунта в условиях криолитозоны. Вычислительный алгоритм базируется на конечно-элементной аппроксимации по пространству на триангулированных сетках Делоне и использовании схемы расщепления по времени с применением линеаризации с предыдущего временного слоя. The article provides mathematical modeling of the complex multiphysical problem relevant for the territories of the Far North and the Arctic. The relevance of this task is characterized by importance of the seepage process in the formation and thawing of the permafrost layer. Modern applications for the most part require consideration of complex geometries, as well as a large number of different processes and their mutual relationship. The multiphysical model consists of the Richards equation to describe the seepage process, the double porosity model to describe natural soil fracturing, the Stefan task to describe the temperature regime of the soil in permafrost zone conditions. The computational algorithm is based on finite-element space approximation on triangulated Delone meshes and using of a time splitting scheme using linearization from a previous time layer.

Modeling Soil Moisture Redistribution and Infiltration Dynamics in Urban Drainage Systems

Hydrological modeling has become a common approach for the design of stormwater management strategies which, at the same time, increasingly rely on permeable and decentralized systems to mitigate the impacts of urbanization. In most applications, little attention is however paid to the temporal variability of infiltration fluxes resulting from soil moisture redistribution between rain events. In this study, a conceptual infiltration-redistribution model is introduced to investigate the importance of the description of infiltration fluxes for the modeling of sustainable urban drainage systems (SUDS). Model verification against numerical solution of Richards equation is first conducted. Performance indicators simulated with the model for a large range of infiltration scenarios are later compared to those obtained with simpler approaches, commonly used in urban hydrology. The model is shown to replicate at a low computational cost numerical solutions of Richards equation. Regarding SUDS modeling, results indicate that a correct description of the temporal variability of infiltration fluxes (1) may be needed for some configurations to assess long-term volume-reduction efficiencies and (2) is more generally required when examining frequency-based performance indicators.