Numerical issues arising in determination of interlayer conductivities in layered unsaturated soils

The one-dimensional (1D) Richards equation is used to model infiltration flow in an unsaturated two-layer soil. Imposing continuities of both flux and pressure head at the interface yields a nonlinear equation for determining interface conductivities. The authors show that multiple solutions of this nonlinear interface equation may exist if the spatial discretization is not fine enough around the interface, in particular as sharp wetting fronts pass through the interface, or for flow across highly dissimilar materials. Three hydraulic models, the Gardner model (G), the Mualem-van Genuchten model (MvG), and the Fredlund-Xing-Leong-Rahardjo model (FXLR), are investigated to demonstrate the nonuniqueness of solutions to the interface problem. For the simplest G model, a full mathematical analysis of the interface problem in terms of two parameters depending on local hydraulic conditions and mesh size is possible. For more advanced models, the interface equation can only be analyzed numerically. In all cases, the interface equation exhibits multiple solutions under suitable conditions, and the impact on the overall numerical solutions is studied. In particular, both a fixed-point iteration and a Newton iteration are used to solve the interface equation numerically.