Scalable and locally mass-conservative simulation of basin-scale variably saturated groundwater flow
Abstract
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.
Related
- Simulation of Groundwater Pollutant Transport Characteristics in Mining Areas Driven by Subsurface Flow Velocity Fields.
- Anthropogenic groundwater declines may exert far-reaching climatic impacts
- The GRACE satellite mission and transformation of global groundwater science
- Global projections of groundwater sustainability through the end of the century
- Groundwater recharge sources and flow systems in the Matsumoto Basin, Japan: Insights from ⁸⁷Sr/⁸⁶Sr isotopes, δ¹⁸O, δ²H, and hydrochemistry
BibTeX
@misc{morrow2026scalable,
title = {{Scalable and locally mass-conservative simulation of basin-scale variably saturated groundwater flow}},
author = {Liam Morrow and Sia Ghelichkhan and Stephan C. Kramer},
howpublished = {EGUsphere},
year = {2026},
month = oct,
doi = {10.31223/x57221},
url = {https://doi.org/10.31223/X57221},
note = {CC BY 4.0},
}
RIS
TY - UNPB
TI - Scalable and locally mass-conservative simulation of basin-scale variably saturated groundwater flow
AU - Liam Morrow
AU - Sia Ghelichkhan
AU - Stephan C. Kramer
PY - 2026
DA - 2026/10/26
PB - EGUsphere
DO - 10.31223/x57221
UR - https://doi.org/10.31223/X57221
AB - 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.
ER -