the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
| 13 Aug 2025
darcyInterTransportFoam v1.0: an open-source, fully-coupled 3D solver for simulating surface water – saturated groundwater processes and exchanges
Abstract. Fully-coupled models have proven to be suitable and computationally efficient tools for studying surface water–groundwater (SW-GW) interactions. Most existing fully-coupled models use the two-dimensional, depth-averaged shallow water equations for surface flows. As a result, three-dimensional (3D) flow dynamics are ignored in the surface domain, including phenomena important for the SW-GW exchange such as turbulence. Computational Fluid Dynamics (CFD) models allow to capture 3D information on the surface turbulent flow by solving the full Navier-Stokes equations. Consequently, they are well-suited for the study of the actual exchange flows of water and solutes across the sediment-water interface. Among the available CFD software, the open-source toolbox OpenFOAM provides a flexible modelling framework to implement user-defined, fully-coupled models for the detailed investigation of SW-GW interaction processes. Based on this CFD platform, Lee et al. (2021) developed hyporheicScalarInterFoam, a fully-coupled, three-dimensional model capable of solving the flow and transport processes in both surface and subsurface domains as well as the interactions across their interface. Despite the potential of this new solver to tackle SW-GW interactions, its application to real-world hydrogeological scenarios is constrained by limitations in boundary conditions, parameter heterogeneity and key hydrodynamic and transport processes, among others, which hinder the accurate representation of the complex characteristics of natural systems. To overcome this, an updated and extended version of hyporheicScalarInterFoam, called darcyInterTransportFoam, is presented in this paper. The new fully-coupled model enhances the applicability of the original by introducing novel simulation features. These include internal solver updates, such as the definition of heterogeneous and anisotropic subsurface fields and the simulation of heat transfer in both domains, as well as newly implemented add-ons, including pre- and post-processing utilities and additional boundary conditions. A complete description of all the new features is provided in this paper. Moreover, the utility of darcyInterTransportFoam is demonstrated in a test case, where the SW-GW flow, solute transport and heat transfer processes are simulated in a highly conductive river-aquifer system.
- Preprint
(3005 KB) - Metadata XML
-
Supplement
(364 KB) - BibTeX
- EndNote
Status: final response (author comments only)
-
RC1: 'Comment on egusphere-2025-3521', Anonymous Referee #1, 10 Sep 2025
-
AC1: 'Reply on RC1', Álvaro Pardo-Álvarez, 30 Jan 2026
We would like to thank the anonymous referee for his thorough review of our manuscript and for his constructive comments. We have carefully considered all remarks and revised the manuscript accordingly. Below, we provide a point-by-point response to each comment:
-
Section 2.1 has been reduced accordingly. The description of the governing equations has been condensed and only modified or newly implemented equations are now shown in the manuscript, while the remaining formulations are referenced to previous publications. Additionally, small modifications were implemented in Section 2.2 (coupling modes) to further reduce its length. All changes are highlighted in the revised manuscript version.
-
All eleven issues identified in the manuscript have been fully addressed in the updated solver. At the same time, we acknowledge that the specific test case selected for this study does not explicitly exercise every one of these improvements. As you correctly identified, points 3, 10, and 11 are not directly demonstrated because they are not relevant to the configuration of the chosen benchmark:
- Point 3 (periodic boundary conditions): This functionality is implemented in the new solver, but it applies only to closed-domain simulations. The test case used in the manuscript is an open system, so periodic boundaries are not invoked.
- Point 10 (interface boundary conditions for one‑way or single‑domain simulations): The presented test case focuses on fully two‑way coupled surface–subsurface flow. As a result, the one‑way or single‑domain configurations covered by this improvement are not part of the scenario.
- Point 11 (runtime reinitialization of data fields): Although this capability is included in the new solver, the selected test case does not require any field reinitialization during the simulation.
To avoid any ambiguity, we have added a table in Section 2.3 which explicitly clarifies that these three points are implemented in the solver but are not triggered by the particular test case used in this study. This ensures that readers understand both the completeness of the improvements and the scope of the demonstrated results.
-
The sentence referred to was misleading at the original location and has, therefore, been relocated to line 143 in the revised manuscript, where it now more explicitly describes the pseudo-transient coupling procedure and the assumptions behind it. Specifically, the sentence clarifies that, owing to the slower response of groundwater compared to surface water, the subsurface solution is assumed to adjust quasi-instantaneously to the surface solution within each time step, thereby avoiding the need for multiple coupling iterations to ensure consistency between surface and subsurface exchanges at the same time level.
-
The volume fraction (alpha) was missing in the diffusion term of Eq. 15 and has now been added in the revised manuscript to ensure consistency with the code and the correct computation of solute diffusion restricted to the water phase.
-
Analogously to the solver’s code, the data associated with the test case was provided as part of the assets for review (https://editor.copernicus.org/EGUsphere/ms-records/egusphere-2025-3521). The link to the data repository is also included in the reference list of the manuscript as “Pardo-Álvarez, Á.: Code and data associated with Pardo-Álvarez et al., 2025, Zenodo [data set], https://doi.org/10.5281/zenodo.16283559, 2025c.”
-
Quantitative validation of variables other than temperature (e.g., flow velocities, water levels or solute concentrations) could not be performed for the selected test case due to the lack of corresponding in-situ observations. Note that the acquisition of velocity fields and the exact, spatially distributed water table constitutes a major observational challenge. Consequently, additional validation plots for these variables cannot be provided. Nevertheless, the simulated flow, temperature and concentration fields are qualitatively consistent with results reported in previous studies conducted at the same site. In addition, the simulated thermal plume compares well with that from drone-based thermal imagery acquired on-site. Together, these qualitative comparisons support the correct functioning of the newly implemented solver.
Citation: https://doi.org/10.5194/egusphere-2025-3521-AC1 -
-
AC1: 'Reply on RC1', Álvaro Pardo-Álvarez, 30 Jan 2026
-
RC2: 'Comment on egusphere-2025-3521', Anonymous Referee #2, 21 Oct 2025
This manuscript presents an updated version of an OpenFOAM-based solver for simulating coupled surface water-groundwater interactions. The work addresses real limitations of the predecessor model (hyporheicScalarInterFoam) and demonstrates improvements through a test case at the Emme River in Switzerland. Following are several comments for your reference:
- Validation and Model Performance Assessment
Only one qualitative comparison is provided (drone thermal imagery vs. model, Fig. 6), which shows notable discrepancies that are attributed to data limitations rather than rigorously investigated. Maybe the comparison with measured flow velocities, water levels, or solute concentrations are missed in the manuscript. Please provide quantitative validation metrics (RMSE, Nash-Sutcliffe efficiency, etc.) for at least the temperature field.
- The steady-state subsurface assumption is a major limitation that needs more thorough discussion. Please provide quantitative criteria for when the steady-state assumption is valid (e.g., dimensionless numbers, timescale ratios).
- The distinction between updates, extensions, and truly novel contributions is unclear. Maybe authors colud create a clear table summarizing: (1) limitations of predecessor, (2) which are addressed, (3) how, and (4) what remains.
- I didn’t see the computational cost analysis about the model despite this being a stated motivation.
- Maybe flow chart is more clearly for Figure 1.
Citation: https://doi.org/10.5194/egusphere-2025-3521-RC2 -
AC2: 'Reply on RC2', Álvaro Pardo-Álvarez, 30 Jan 2026
We would like to thank the anonymous referee for the time and effort he devoted to reviewing our manuscript. We have carefully considered all his valuable comments and suggestions and revised the manuscript accordingly. Below, we address each comment point by point:
-
A quantitative validation metric has been added in Section 4.2 by calculating the root mean square error (RMSE) of the simulated temperature field. This metric complements the qualitative comparison with drone-based thermal imagery and further supports the correct functioning of the newly implemented heat transfer module, based on the spatial heat pattern of the plume generated by the discharged groundwater. Quantitative comparisons with additional variables (e.g., flow velocities, water levels or solute concentrations) could not be performed due to the lack of corresponding field measurements for the selected test case, as already indicated in the manuscript.
-
A quantitative criterion to assess when steady-state conditions can approximate transient system response to periodic stresses has been added in Section 5, following the guidelines proposed by Haitjema (2006). This approach calculates a dimensionless aquifer response time, τ, given by:
τ = (1/4*P)*(S*L^2/K*b)
where L is the characteristic system length defined as the distance between major surface water features (m), P is the period of the forcing function (=365 days for seasonal fluctuations), K is the hydraulic conductivity (m/s), b is the saturated thickness (m) and S the is aquifer storativity (specific yield, Sy, for an unconfined aquifer).
According to the resulting value of , a steady-state approximation will be appropriate for τ > 1.0, while a transient model will be required for 0.1 < τ < 1.0 and a bounding or successive steady-state solution may be used for τ < 0.1.
An expanded discussion of this timescale classification is provided in the manuscript to clarify the conditions under which the steady-state assumption is valid.
-
A new table has been added in Section 2.3 including the suggested information, namely: (1) the limitations of the original solver, (2) which of these limitations are addressed in the new code, (3) the type of update (new or extended functionality) and (4) whether or not each new feature is applied in the test case of this study. The added table clearly distinguishes between extensions and newly introduced features, thereby clarifying the novel contributions of the upgraded solver.
-
A comprehensive and systematic computational cost analysis, including direct performance comparisons with previous versions of the code and with other solvers of the same class, is beyond the scope of the present study. The primary focus of this paper is on the description and verification of the newly introduced modelling capabilities. Nevertheless, to address practical considerations related to computational cost, Section 5 provides guidance on the computational resources required to run efficient simulations with the new solver, thereby assisting users in selecting appropriate hardware for their applications.
-
Figure 1 was designed as a flow chart to illustrate the workflow of the newly developed solver. We have reviewed the figure carefully and believe that its current form already fulfils this purpose. Could you, please, clarify this point?
Citation: https://doi.org/10.5194/egusphere-2025-3521-AC2 -
Data sets
Code and data associated with Pardo-Álvarez et al., 2025 Álvaro Pardo-Álvarez https://doi.org/10.5281/zenodo.16283559
Model code and software
darcyInterTransportFoam v1.0 Álvaro Pardo-Álvarez https://doi.org/10.5281/zenodo.15857142
MATLAB-based OpenFOAM mesh and data dictionary generator Álvaro Pardo-Álvarez https://doi.org/10.5281/zenodo.15857296
Viewed
| HTML | XML | Total | Supplement | BibTeX | EndNote | |
|---|---|---|---|---|---|---|
| 1,948 | 300 | 44 | 2,292 | 54 | 34 | 36 |
- HTML: 1,948
- PDF: 300
- XML: 44
- Total: 2,292
- Supplement: 54
- BibTeX: 34
- EndNote: 36
Viewed (geographical distribution)
| Country | # | Views | % |
|---|
| Total: | 0 |
| HTML: | 0 |
| PDF: | 0 |
| XML: | 0 |
- 1
This manuscript presents an updated version of a previous OpenFOAM solver on surface water-groundwater interactions (SW-GW). The authors described the details on what have been updated/changed/added and showed a test case. I read the manuscript with interest and found it might be a good contribution to the modeling community. I have the following comments/questions:1