GWSWEX v1.0: a dual-solver 1D unsaturated zone model for mass-conservative groundwater recharge and runoff computation in distributed hydrological modelling

Abstract. Faithful numerical representation of the coupled groundwater (GW)–unsaturated-zone (UZ)–surface-water (SW) dynamics at regional scales is constrained by the structural extremes currently available: fully three-dimensional Richards-equation solvers resolve variably saturated flow rigorously but at a per-evaluation cost that forecloses ensemble calibration and multi-decadal projection, whereas conceptual bucket schemes cannot track a moving groundwater head, represent capillary rise, or impose a soil-moisture-dependent evapotranspiration (ET) stress. This study introduces GWSWEX (Groundwater–Surface Water EXchange), a vertically resolved, process-based UZ model designed to sit between an external GW model and an external SW model within an integrated modelling chain. Its principal innovations are two interchangeable solvers behind a single unified API: an explicit operator-split bucket-sequence solver with CFL-adaptive sub-stepping and an implicit mixed-form Richards solver with Picard linearisation, sharing the same Mualem–van Genuchten constitutive relations and piecewise-linear drainable-volume function for a freely moving groundwater head; the per-step storage budget closes to floating-point tolerance in the explicit solver and to a small, well-characterised residual in the implicit solver, enabling a coupled GW–UZ–SW iteration to converge without under-relaxation. Verification against HYDRUS-1D across six soil profiles under two contrasting forcing scenarios shows sub-centimetre to near-centimetre groundwater-head accuracy under smooth daily forcing, with eleven of twelve basic-scenario cases attaining an RMSE at or below the 2 cm discretisation step and a wet-phase mean RMSE of approximately 4 cm (implicit) and 8 cm (explicit) under intensive hourly forcing; an iterated one-at-a-time sensitivity analysis characterises the empirical parameters and the structural limits of the layered-bucket abstraction, while a parallel-ensemble benchmark shows both solvers completing a thousand-column ensemble six to twelve times faster than parallelised HYDRUS-1D, at sub-linear run-time scaling exponents of approximately 0.92–0.93 that render ensemble-grade calibration and multi-decadal projection tractable on mainstream consumer-grade hardware.