Scalable and locally mass-conservative simulation of basin-scale variably saturated groundwater flow

Groundwater is the largest store of liquid fresh water, yet because it is hidden beneath the surface, it is primarily understood through models of porous-media flow. Richards' equation, a physically faithful description of variably saturated groundwater flow, remains costly for basin-scale groundwater modelling. Consequently, operational models use reduced physics, either by assuming linearity or by using approximate solutions. Here, we introduce a scalable, open-source, three-dimensional solver for variably saturated groundwater flow implemented in G-ADOPT, an automated finite-element platform built on Firedrake. The implementation combines a locally conservative discontinuous Galerkin discretisation, Newton iteration with an exact Jacobian, implicit time integration through the Irksome library, used here with Backward Euler, and vertically lumped multigrid preconditioners designed for the extreme aspect ratios of production scenarios. We verify accuracy, convergence, and strict local mass conservation on two- and three-dimensional benchmarks, and apply the solver to the Lower Murrumbidgee floodplain in New South Wales, Australia. The model covers approximately 30 000 km^2, contains several hundred million degrees of freedom, and uses public stratigraphic, rainfall, and water-table data. Our scaling experiments show that the effectiveness of the preconditioner is controlled primarily by the physical regime in which the model is integrated. A single-level method (e.g., block-Jacobi ILU(0)) suffices at the short time steps of an infiltration study, but at seasonal time steps, the admissible time step decreases significantly under horizontal refinement. In a near-saturated profile without specific storage, such single-level methods do not converge at any tested step size. By contrast, a vertically lumped coarse-space preconditioner takes the full three-month step on every mesh, requiring four to ten linear iterations per Newton step and reducing wall-clock time by a factor of two to four. It also integrates the near-saturated profile on every mesh. Because basin simulations use seasonal scenarios, we report, for the first time with a discontinuous Galerkin discretisation and a vertically lumped multigrid preconditioner, a seasonal, basin-scale, fully three-dimensional solution of the full Richards' equation on an unstructured, terrain-following mesh. This study demonstrates basin-scale capability rather than predictive skill for the synthetic Murrumbidgee test case, whose hydraulic parameters require calibration. The algorithmic differentiation in the underlying Firedrake and pyadjoint frameworks supports future inverse and data-assimilation studies, including the calibration needed to assess predictive skill.