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

Abstract. The faithful numerical representation of the coupled dynamics between groundwater (GW), the unsaturated zone (UZ) and surface water (SW) remains one of the more persistent challenges of regional-scale integrated hydrological modelling. Physically-based models that solve the Richards equation in three dimensions provide a rigorous description of vertical fluxes, but their computational cost, sensitivity to parameter uncertainty, and tendency to over-emphasise capillary forces at the expense of preferential and non-equilibrium flow regimes limit their suitability as self-contained UZ modules for regional integrated hydrological models, ensemble simulations, and multi-decadal climate-impact studies. Conversely, the conceptual and water-balance bucket models routinely adopted at regional scales rarely retain enough vertical resolution to track the position of a moving groundwater head, to represent capillary rise from the GW into a depleted root zone, or to impose a physically consistent evapotranspiration (ET) stress that responds to the local soil-moisture state. This paper introduces GWSWEX (Groundwater–Surface Water EXchange), a vertically resolved process-based modelling package designed to occupy the middle ground between these two extremes and to act as a UZ coupler between an external GW model and an external SW model within an integrated modelling chain. GWSWEX represents each model element as a layered 1-D soil column, tracks GW elevation, per-layer UZ storage, and SW ponding as prognostic states, and exposes two interchangeable numerical solvers behind a unified Python API: an explicit operator-split bucket-sequence solver with CFL-adaptive sub-stepping, and an implicit mixed-form Richards solver with Picard linearisation and Thomas-algorithm tridiagonal inversion. Both solvers share the same spatial discretisation, the same Mualem–van Genuchten retention and conductivity relations, the same vertically integrated drainable-volume function for the moving groundwater head, and the same atmospheric and lateral boundary conditions. The model is implemented as a modern Fortran 2008 kernel with OpenMP parallelism over elements, exposed to Python via f2py, and configured through a Pydantic-validated user-facing API. Verification against HYDRUS-1D across six soil profiles and two contrasting forcing scenarios, together with an iterated one-at-a-time sensitivity analysis of the empirical parameters of both solvers and a parallel-ensemble computational performance benchmark, demonstrates sub-centimetre to near-centimetre GWH accuracy across smooth and intensive forcing regimes, characterises the physical limits of the layered-bucket abstraction, and examines the computational cost of the model at the ensemble sizes that regional integrated modelling requires.