Richards equation

The Richards equation represents the movement of water in unsaturated soils, and is attributed to Lorenzo A. Richards who published the equation in 1931.[1] It is a quasilinear partial differential equation; its analytical solution is often limited to specific initial and boundary conditions.[2] Proof of the existence and uniqueness of solution was given only in 1983 by Alt and Luckhaus.[3] The equation is based on Darcy-Buckingham law[1] representing flow in porous media under variably saturated conditions, which is stated as

${\displaystyle {\vec {q}}=-\mathbf {K} (\theta )(\nabla h+\nabla z),}$

where

${\displaystyle {\vec {q}}}$ is the volumetric flux;

${\displaystyle \theta }$ is the volumetric water content;

${\displaystyle h}$ is the liquid pressure head, which is negative for unsaturated porous media;

${\displaystyle \mathbf {K} (h)}$ is the unsaturated hydraulic conductivity;

${\displaystyle \nabla z}$ is the geodetic head gradient, which is assumed as ${\displaystyle \nabla z=\left({\begin{smallmatrix}0\\0\\1\end{smallmatrix}}\right)}$ for three-dimensional problems.

Considering the law of mass conservation for an incompressible porous medium and constant liquid density, expressed as

${\displaystyle {\frac {\partial \theta }{\partial t}}+\nabla \cdot {\vec {q}}+S=0}$,

where

${\displaystyle S}$ is the sink term [T${\displaystyle ^{-1}}$], typically root water uptake.[4]

Then substituting the fluxes by the Darcy-Buckingham law the following mixed-form Richards equation is obtained:

${\displaystyle {\frac {\partial \theta }{\partial t}}=\nabla \cdot \mathbf {K} (h)(\nabla h+\nabla z)-S}$.

For modeling of one-dimensional infiltration this divergence form reduces to

${\displaystyle {\frac {\partial \theta }{\partial t}}={\frac {\partial }{\partial z}}\left(\mathbf {K} (\theta )\left({\frac {\partial h}{\partial z}}+1\right)\right)-S}$.

Although attributed to L. A. Richards, the equation was originally introduced 9 years earlier by Lewis Fry Richardson in 1922.[5][6]

Formulations

The Richards equation appears in many articles in the environmental literature because it describes the flow in the vadose zone between the atmosphere and the aquifer. It also appears in pure mathematical journals because it has non-trivial solutions. The above-given mixed formulation involves two unknown variables: ${\displaystyle \theta }$ and ${\displaystyle h}$. This can be easily resolved by considering constitutive relation ${\displaystyle \theta (h)}$, which is known as the water retention curve. Applying the chain rule, the Richards equation may be reformulated as either ${\displaystyle h}$-form (head based) or ${\displaystyle \theta }$-form (saturation based) Richards equation.

Head-based

By applying the chain rule on temporal derivative leads to

${\displaystyle {\frac {\partial \theta (h)}{\partial t}}={\frac {{\textrm {d}}\theta }{{\textrm {d}}h}}{\frac {\partial h}{\partial t}}}$,

where ${\displaystyle {\frac {{\textrm {d}}\theta }{{\textrm {d}}h}}}$ is known as the retention water capacity ${\displaystyle C(h)}$. The equation is then stated as

${\displaystyle C(h){\frac {\partial h}{\partial t}}=\nabla \cdot \left(\mathbf {K} (h)\nabla h+\nabla z\right)-S}$.

The head-based Richards equation is prone to the following computational issue: the discretized temporal derivative using the implicit Rothe method yields the following approximation:

${\displaystyle {\frac {\Delta \theta }{\Delta t}}\approx C(h){\frac {\Delta h}{\Delta t}},\quad {\mbox{and so}}\quad {\frac {\Delta \theta }{\Delta t}}-C(h){\frac {\Delta h}{\Delta t}}=\varepsilon .}$

This approximation produces an error ${\displaystyle \varepsilon }$ that affects the mass conservation of the numerical solution, and so special strategies for temporal derivatives treatment are necessary.[7]

Saturation-based

By applying the chain rule on the spatial derivative leads to

${\displaystyle \mathbf {K} (h)\nabla h=\mathbf {K} (h){\frac {{\textrm {d}}h}{{\textrm {d}}\theta }}\nabla \theta ,}$

where ${\displaystyle \mathbf {K} (h){\frac {{\textrm {d}}h}{{\textrm {d}}\theta }}}$, which could be further formulated as ${\displaystyle {\frac {\mathbf {K} (\theta )}{C(\theta )}}}$, is known as the soil water diffusivity ${\displaystyle \mathbf {D} (\theta )}$. The equation is then stated as

${\displaystyle {\frac {\partial \theta }{\partial t}}=\nabla \cdot \mathbf {D} (\theta )\nabla \theta -S.}$

The saturation-based Richards equation is prone to the following computational issues. Since the limits ${\displaystyle \lim _{\theta \to \theta _{s}}||\mathbf {D} (\theta )||=\infty }$ and ${\displaystyle \lim _{\theta \to \theta _{r}}||\mathbf {D} (\theta )||=\infty }$, where ${\displaystyle \theta _{s}}$ is the saturated (maximal) water content and ${\displaystyle \theta _{r}}$ is the residual (minimal) water content a successful numerical solution is restricted just for ranges of water content satisfactory below the full saturation (the saturation should be even lower than air entry value) as well as satisfactory above the residual water content.[8]

Parametrization

The Richards equation in any of its forms involves soil hydraulic properties, which is a set of five parameters representing soil type. The soil hydraulic properties typically consist of water retention curve parameters by van Genuchten:[9] (${\displaystyle \alpha ,\,n,\,m,\,\theta _{s},\theta _{r}}$), where ${\displaystyle \alpha }$ is the inverse of air entry value [L^−1], ${\displaystyle n}$ is the pore size distribution parameter [-], and ${\displaystyle m}$ is usually assumed as ${\displaystyle m=1-{\frac {1}{n}}}$. Further the saturated hydraulic conductivity ${\displaystyle \mathbf {K} _{s}}$ (which is for non isotropic environment a tensor of second order) should also be provided. Identification of these parameters is often non-trivial and was a subject of numerous publications over several decades.[10][11][12][13][14][15]

Limitations

The numerical solution of the Richards equation is one of the most challenging problems in earth science.[16] Richards' equation has been criticized for being computationally expensive and unpredictable [17][18] because there is no guarantee that a solver will converge for a particular set of soil constitutive relations. Advanced computational and software solutions are required here to over-come this obstacle. The method has also been criticized for over-emphasizing the role of capillarity,[19] and for being in some ways 'overly simplistic' [20] In one dimensional simulations of rainfall infiltration into dry soils, fine spatial discretization less than one cm is required near the land surface,[21] which is due to the small size of the representative elementary volume for multiphase flow in porous media. In three-dimensional applications the numerical solution of the Richards equation is subject to aspect ratio constraints where the ratio of horizontal to vertical resolution in the solution domain should be less than about 7.^[citation needed]

References

1. Richards, L.A. (1931). "Capillary conduction of liquids through porous mediums". Physics. 1 (5): 318–333. Bibcode:1931Physi...1..318R. doi:10.1063/1.1745010.
2. Tracy, F. T. (August 2006). "Clean two- and three-dimensional analytical solutions of Richards' equation for testing numerical solvers: TECHNICAL NOTE". Water Resources Research. 42 (8). doi:10.1029/2005WR004638. S2CID 119938184.
3. Wilhelm Alt, Hans; Luckhaus, Stephan (1 September 1983). "Quasilinear elliptic-parabolic differential equations". Mathematische Zeitschrift. 183 (3): 311–341. doi:10.1007/BF01176474. ISSN 1432-1823. S2CID 120607569.
4. Feddes, R. A.; Zaradny, H. (1 May 1978). "Model for simulating soil-water content considering evapotranspiration — Comments". Journal of Hydrology. 37 (3): 393–397. Bibcode:1978JHyd...37..393F. doi:10.1016/0022-1694(78)90030-6. ISSN 0022-1694.
5. Knight, John; Raats, Peter. "The contributions of Lewis Fry Richardson to drainage theory, soil physics, and the soil-plant-atmosphere continuum" (PDF). EGU General Assembly 2016.
6. Richardson, Lewis Fry (1922). Weather prediction by numerical process. Cambridge, The University press. pp. 262.
7. Celia, Michael A.; Bouloutas, Efthimios T.; Zarba, Rebecca L. (July 1990). "A general mass-conservative numerical solution for the unsaturated flow equation". Water Resources Research. 26 (7): 1483–1496. Bibcode:1990WRR....26.1483C. doi:10.1029/WR026i007p01483.
8. Kuráž, Michal; Mayer, Petr; Lepš, Matěj; Trpkošová, Dagmar (2010-04-15). "An adaptive time discretization of the classical and the dual porosity model of Richards' equation". Journal of Computational and Applied Mathematics. Finite Element Methods in Engineering and Science (FEMTEC 2009). 233 (12): 3167–3177. Bibcode:2010JCoAM.233.3167K. doi:10.1016/j.cam.2009.11.056. ISSN 0377-0427.
9. van Genuchten, M. Th. (September 1980). "A Closed-form Equation for Predicting the Hydraulic Conductivity of Unsaturated Soils". Soil Science Society of America Journal. 44 (5): 892–898. Bibcode:1980SSASJ..44..892V. doi:10.2136/sssaj1980.03615995004400050002x.
10. Inoue, M.; Šimůnek, J.; Shiozawa, S.; Hopmans, J.W. (June 2000). "Simultaneous estimation of soil hydraulic and solute transport parameters from transient infiltration experiments". Advances in Water Resources. 23 (7): 677–688. Bibcode:2000AdWR...23..677I. doi:10.1016/S0309-1708(00)00011-7.
11. Fodor, Nándor; Sándor, Renáta; Orfanus, Tomas; Lichner, Lubomir; Rajkai, Kálmán (October 2011). "Evaluation method dependency of measured saturated hydraulic conductivity". Geoderma. 165 (1): 60–68. Bibcode:2011Geode.165...60F. doi:10.1016/j.geoderma.2011.07.004.
12. Angulo-Jaramillo, Rafael; Vandervaere, Jean-Pierre; Roulier, Stéphanie; Thony, Jean-Louis; Gaudet, Jean-Paul; Vauclin, Michel (May 2000). "Field measurement of soil surface hydraulic properties by disc and ring infiltrometers". Soil and Tillage Research. 55 (1–2): 1–29. doi:10.1016/S0167-1987(00)00098-2.
13. Köhne, J. Maximilian; Mohanty, Binayak P.; Šimůnek, Jirka (January 2006). "Inverse Dual-Permeability Modeling of Preferential Water Flow in a Soil Column and Implications for Field-Scale Solute Transport". Vadose Zone Journal. 5 (1): 59–76. Bibcode:2006VZJ.....5...59K. doi:10.2136/vzj2005.0008. ISSN 1539-1663. S2CID 781417.
14. Younes, Anis; Mara, Thierry; Fahs, Marwan; Grunberger, Olivier; Ackerer, Philippe (3 May 2017). "Hydraulic and transport parameter assessment using column infiltration experiments". Hydrology and Earth System Sciences. 21 (5): 2263–2275. Bibcode:2017HESS...21.2263Y. doi:10.5194/hess-21-2263-2017. ISSN 1607-7938.
15. Kuraz, Michal; Jačka, Lukáš; Ruth Blöcher, Johanna; Lepš, Matěj (1 November 2022). "Automated calibration methodology to avoid convergence issues during inverse identification of soil hydraulic properties". Advances in Engineering Software. 173: 103278. doi:10.1016/j.advengsoft.2022.103278. ISSN 0965-9978. S2CID 252508220.
16. Farthing, Matthew W., and Fred L. Ogden, (2017). Numerical solution of Richards’ Equation: a review of advances and challenges. Soil Science Society of America Journal, 81(6), pp.1257-1269.
17. Short, D., W.R. Dawes, and I. White, (1995). The practicability of using Richards' equation for general purpose soil-water dynamics models. Envir. Int'l. 21(5):723-730.
18. Tocci, M. D., C. T. Kelley, and C. T. Miller (1997), Accurate and economical solution of the pressure-head form of Richards' equation by the method of lines, Adv. Wat. Resour., 20(1), 1–14.
19. Germann, P. (2010), Comment on “Theory for source-responsive and free-surface film modeling of unsaturated flow”, Vadose Zone J. 9(4), 1000-1101.
20. Gray, W. G., and S. Hassanizadeh (1991), Paradoxes and realities in unsaturated flow theory, Water Resour. Res., 27(8), 1847-1854.
21. Downer, Charles W., and Fred L. Ogden (2003), Hydrol. Proc.,18, pp. 1-22. DOI:10.1002/hyp.1306.

See also

• Infiltration (hydrology)
• Water retention curve
• Finite water-content vadose zone flow method
• Soil Moisture Velocity Equation