An OpenFOAM®-based coupled surface-subsurface flow model for simulating watershed hydrodynamics

Watershed hydrodynamics is governed by various hydrological flow processes that occur at different spatiotemporal scales. Most hydrological models couple the surface flow solver with the standard empirical infiltration models for flood propagation modeling. However, the empirical infiltration models are not applicable for heterogeneous and anisotropic soils and shallow groundwater tables, which are most vulnerable to waterlogging problems. Hence, simultaneous and integrated modeling of the surface and subsurface flow processes is essential for the continuous monitoring of watershed hydrodynamics. A physically based unified multi-region, multi-process watershed model integrates the various hydrological flow components in different regions through unique coupling mechanisms at the interfaces. The current work presents a Finite Volume (FV) method-based watershed flow model developed using the OpenFOAM^® framework [1]. The developed model framework utilizes the ‘multi-region’ structure from the OpenFOAM^® library to integrate the OpenFOAM^®-based solvers for the individual processes of surface overland flow [2,3] and saturated-unsaturated subsurface flow [4] through the imposition of appropriate interface boundary conditions or addition of source/sink terms at the interfaces of the flow regions. The surface flow component is modeled using the diffusive wave or the zero-inertia (ZI) approximation of the two-dimensional (2D) depth-averaged shallow water equations (SWE). On the other hand, the flow through the variably saturated subsurface media is modeled using the ‘mixed form’ of the 3D modified Richards Equation. The flux exchange between the surface and subsurface regions (infiltration or exfiltration rate) is modeled using a switching algorithm to impose the boundary condition on the interface between the two regions. The algorithm changes the interface to a Dirichlet or a Neumann type boundary condition based on the rainfall intensity and the saturated hydraulic conductivity of the ground surface. A stabilized and adaptive time-stepping algorithm has been implemented to ensure smooth convergence of the iterative technique used for linearizing the nonlinear governing equations. The developed model is equipped with parallelization strategies to be run on multi-core processors, which is essential for increased computational efficiency while solving regional-scale watershed flow problems. The developed watershed model has been verified and validated against the standard benchmark problems on saturation excess and infiltration excess from the literature. Moreover, the applicability of the developed model has been extended to solve complex hydrological problems on exfiltration occurring over natural catchments, yielding satisfactory results.

References

[1] Jasak, H., A. Jemcov, Z. Tukovic. (2007). OpenFOAM: A C++ library for complex physics simulations. In Vol. 1000 of Proc., Int. Workshop on Coupled Methods in Numerical Dynamics,1–20. Dubrovnik, Croatia: Inter-University Center

[2] Dey, S., Dhar, A. (2024). Applicability of Zero-Inertia Approximation for Overland Flow Using a Generalized Mass-Conservative Implicit Finite Volume Framework. Journal of Hydrologic Engineering, 29(1), 04023042.

[3] Dey, S. (2025). zeroInertiaFlowFOAM – a OpenFOAM^®-based computationally efficient, mass-conservative, implicit zero-inertia flow model for flood inundation problems on collocated grid-systems (No. EGU25-17402). Copernicus Meetings.

[4] Dey, S., & Dhar, A. (2022). Generalized mass-conservative finite volume framework for unified saturated–unsaturated subsurface flow. Journal of Hydrology, 605, 127309.