A scalable and implicit meshless RBF method for the 3D unsteady nonlinear Richards equation with single and multi-zone domains

The Local Hermitian Interpolation (LHI) method is a strong-form meshless numerical technique which uses Radial Basis Function (RBF) interpolants to satisfy linear and nonlinear governing equations and boundary operators. Recent developments have shown that, for linear transport problems, applying the PDE governing equation directly to the basis functions can greatly improve the accuracy and stability of the resulting solutions. In this work, the LHI formulation with local PDE-interpolation is extended to the nonlinear gravity-driven Richards equation, in order to solve unsteady problems involving flow in unsaturated porous media. The application of the linearised PDE-operator to the basis functions incorporates information, such as the effective velocity field, directly into the solution construction. This results in a form of 'analytical upwinding' which helps to stabilise the solution. In addition, the local interpolation itself satisfies the linearised governing equation, allowing for more accurate reconstruction of partial derivatives, and hence a more accurate solution. The procedure is tested using a 3D infiltration problem with a known analytical solution. The performance of the LHI method, both with and without local PDE interpolation, is compared to the Finite Element method via the FEMWATER software. The LHI formulation with PDE data centres shows consistent improvement over the FEMWATER solutions, reducing errors by several orders of magnitude. In addition, a procedure is introduced to model strongly heterogeneous and layered soils. The physically correct matching conditions are applied over layer interfaces, i.e. continuity of pressure and mass flux. The 'double collocation' property of the Hermitian RBF method is exploited to enforce both matching conditions at the same set of locations on the layer interface. This procedure allows the accurate capture of solutions across such interfaces, replicating the required discontinuities in the first derivatives of the pressure profile. The multi-layer formulation is validated using a transient two-layer infiltration problem, with the analytical solution replicated to a high precision in a variety of configurations. Copyright (C) 2010 John Wiley & Sons, Ltd.