Why modified Picard works: Mass-conservative and higher-order time integration for the Richardson–Richards equation

We show that the mass-conservative modified Picard scheme of Celia et al. (1990) for the Richardson–Richards equation is exactly the Newton–Schur reduction of a differential–algebraic (DAE) formulation in which moisture content is retained and linked to pressure head by a constitutive constraint. The reduced DAE formulation yields a pressure head-only residual and generalizes to higher-order multistep schemes (e.g., backward differentiation formulas) while preserving discrete mass balance. Benchmarks on one-dimensional infiltration problems compare the classical mixed scheme with the proposed higher-order DAE and Schur-reduced formulations, with higher-order mass-conservative formulations achieving lower errors at a given time step size.