Numerical simulation of infiltration, evaporation and shallow groundwater levels with the Richards equation

Analysis of water and solute movement in unsaturated–saturated soil systems would greatly benefit from an accurate and efficient numerical solution of the Richards equation. Recently the mass balance problem has been solved by proper evaluation of the water capacity term. However, the Darcy fluxes as calculated by various numerical schemes still deviate significantly due to differences in nodal spacing and spatial averaging of the hydraulic conductivity K. This paper discusses a versatile, implicit, backward finite difference scheme which is relatively easy to implement. Special attention is given to the selection of a head or flux controlled top boundary condition during the iterative solution of the Richards equation. The stability of the scheme is shown for extreme events of infiltration, evaporation and rapidly fluctuating, shallow groundwater levels in case of two strongly non-linear soils. For nodal distances of 5 cm, arithmetic means of K overestimate the soil water fluxes, while geometric means of K underestimate these fluxes. At smaller nodal distances, arithmetic means of K converge faster to the theoretical solution than geometric means. In case of nodal distances of 1 cm and arithmetic averages of K, errors due to numerical discretization are small compared to errors due to hysteresis and horizontal spatial variability of the soil hydraulic functions.