Modelling non-equillibrium unsaturated flow in soils during sudden pressure head changes by solving Richards' equation with a Method of lines approach

The classical formulation of Richards' equation is relying on a unique functional relationship between water content, conductivity and pressure head. Some phenomena like hystersis effects in the water content during wetting and drying cycles and hydraulic non-equillibrium cannot be accounted for with this formulation. Therefor it has been extended in different ways in the past to be able to include these effects in the simulation. Each modification comes with its own challenges regarding implementation and numerical stability.
The Method Of Lines approach to solving the Richards' equation has already be shown to be an efficient and stable alternative to established solution methods, such as low-order finite difference and finite element methods applied to the mixed form of Richards' equation.
In this work a slightly modified Method Of Lines approach is used to solve the pressure based 1D Richards' equation. A finite differencing scheme is applied to the spatial derivative and the resulting system of ordinary differential equations is reformulated as differential-algebraic system of equations. The open-source code IDAS from the Sundials suite is used to solve the DAE system. Different extensions to Richards' equation have been incorporated into the model to address the shortcomings mentioned above. These extensions are a model able to simulate preferential flow using a coupled two domain approach, a simple hysteretic model to account for hysteresis in the water retention curve and also two models to either fully or partially calculate hydraulic non-equillibrium effects. To verify the numerical robustness of the extended model, stochastic parameterizations were generated that represent the full range of all soil types. Simulations were carried out using these parameter sets and real-world meteorological boundary conditions at 10 minutes time intervals, that exhibit drastic flux changes and poses numerical challenges for classical solution methods.

The results show that not only does the extended model converge for all parameterizations, but that numerical robustness and performance is maintained. Where applicable the results have been verified against solutions from the software Hydrus and show good agreement with those.