the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
| 09 Oct 2025
High-performance coupled surface-subsurface flow simulation with SERGHEI-SWE-RE
Abstract. This work presents SERGHEI-SWE-RE, a performance-portable, parallel model that couples a fully dynamic two-dimensional Shallow Water Equation (SWE) solver with a three-dimensional Richards Equation (RE) solver within the Kokkos framework to simulate surface–subsurface flow exchange. The model features a modular architecture with sequential coupling strategy, supporting both synchronous and asynchronous executions of surface and subsurface modules. The SERGHEI-SWE-RE model is validated against five benchmark problems incorporating stationary and fluctuating free-surface tests, a tilted v-catchment, a lateral-flow slope without ponding, and a heterogeneous superslab. The results demonstrate good agreement with established models. Asynchronous coupling reduces wall-clock time by up to about 60 % in the superslab case while preserving simulation accuracy. Strong and weak scaling tests on multiple Intel Xeon CPUs and NVIDIA GPUs reveal robust portability, with near-ideal RE scaling and less-satisfactory SWE scaling at high GPU counts, suggesting future improvements on differentiated meshes or more advanced domain decomposition strategies. Overall, the results presented establish SERGHEI-SWE-RE as an efficient, flexible and scalable model for integrated surface-subsurface flow simulations.
-
Notice on discussion status
The requested preprint has a corresponding peer-reviewed final revised paper. You are encouraged to refer to the final revised version.
-
Preprint
(7116 KB)
-
The requested preprint has a corresponding peer-reviewed final revised paper. You are encouraged to refer to the final revised version.
- Preprint
(7116 KB) - Metadata XML
- BibTeX
- EndNote
- Final revised paper
Journal article(s) based on this preprint
Interactive discussion
Status: closed
-
RC1: 'Comment on egusphere-2025-4246', Anonymous Referee #1, 24 Nov 2025
- AC1: 'Reply on RC1', Zhi Li, 22 Dec 2025
-
RC2: 'Comment on egusphere-2025-4246', Anonymous Referee #2, 21 Mar 2026
The paper presents SERGHEI-SWE-RE, a coupled surface–subsurface model built on the Kokkos framework that targets performance portability across CPUs and GPUs. The work is clearly within GMD's scope, and I found the asynchronous coupling strategy and the multi-platform scalability analysis to be the most interesting contributions. The benchmark suite covers a reasonable range of conditions, and the comparisons with ParFlow, CATHY, HGS and others are generally convincing. That said, I have a few concerns I would like the authors to address. Two of them are substantive and touch on aspects of the model that I think are currently under-described; the third is a straightforward request for numbers that the authors say they have but chose not to show. None of these require new simulations. My recommendation is minor revision, but I want to be clear that I expect concrete answers rather than a restatement of what is already in the paper.
Major Comments:
My main technical concern is with Section 2.2. The four-case boundary condition scheme is described clearly enough, but the manuscript is completely silent on what happens numerically when the system moves between these cases, for instance when a dry cell receives exfiltration and becomes wet, or when the top subsurface layer saturates and the pressure head approaches the surface elevation. These are precisely the transitions that are known to cause convergence problems in this class of models. Coon et al. (2020) document this explicitly for the ATS model, noting that convergence degrades near the dry-to-ponding transition, and Furman (2008) gives a broader review of the numerical difficulties associated with boundary condition switching. The present manuscript does not acknowledge any of this.
I should also point out that more elegant solutions to this problem exist in the literature and are not cited. Casulli (2017) derived a coupled surface–subsurface formulation in which the governing equations for both domains are solved simultaneously using the NCZ algorithm, the same algorithm on which SERGHEI-RE is based, entirely avoiding the switching problem. Earlier, Casulli (2009) introduced the idea of an additional computational node at the soil surface governed by the equation state H(ψ), which allows rainfall to be prescribed as a permanent Neumann boundary condition regardless of ponding state. Tubini and Rigon (2022) implemented exactly this approach in WHETGEO-1D and showed that it handles both infiltration-excess and saturation-excess ponding without any case logic or numerical oscillations. I am not asking the authors to rewrite their code. I am asking them to engage with this literature and explain, in the paper, whether their four-case scheme has been tested under conditions of rapid state transitions and what its behaviour is. If fully saturated conditions with return flow are never actually reached in any of the five benchmark cases, that should be stated explicitly, because right now the reader cannot tell.
The introduction lists modularisation and flexibility as key advantages of SERGHEI-SWE-RE. I agree the internal modularisation between the SWE and RE solvers is well demonstrated, but the paper says nothing about whether the model can be coupled with external codes. This matters because the community has largely converged on standard interfaces for this purpose, the Basic Model Interface (BMI; Hutton et al., 2020) being the most widely adopted, and models that do not address interoperability are increasingly difficult to integrate into larger workflows. WHETGEO (Tubini and Rigon, 2022) does this through OMS3. The Next-Generation Water Resources Modeling Framework uses BMI natively. The paper should state clearly whether SERGHEI-SWE-RE supports any such interface, and if not, acknowledge this as a current limitation rather than implying that flexibility is already there.
A related point: the consistent domain decomposition (same MPI partition for SWE and RE) is already shown in Section 4 to be a bottleneck at high GPU counts. This is an inherent consequence of tying the two solvers' parallelisation together, and it would also complicate any future attempt to expose an external coupling interface. The authors discuss this briefly but do not connect it to the broader architectural question of whether the design supports or hinders future extension. A paragraph addressing this would strengthen the Discussion considerably.
Lines 173–176 state that "further quantitative examination (not shown) confirms" that the C-property is satisfied. This is not acceptable in a model description paper submitted to GMD. The C-property is a fundamental numerical requirement and its verification should be documented with actual numbers, the maximum absolute deviation from the initial condition is the standard metric and takes one line to report. The authors clearly have these numbers. They should include them, either in the text or in a supplementary table. Reviewer 1 made the same point and I fully agree.
Specific and Minor Comments:
The rainfall routing assumption (Section 2.2, lines 140–143) is validated only against Case 4, which was designed with no ponding. This is not a meaningful test of the assumption. The authors should acknowledge the scenarios where routing rain through the surface first introduces a numerical lag, the obvious case being a very dry soil with high infiltration capacity.
On the asynchronous coupling: Case 5 tests two values of Δtsub,max and both give satisfactory results for that specific benchmark. But the paper offers no guidance on how to choose this parameter in general. Users running scenarios with fast-moving wetting fronts or rapidly alternating wet/dry conditions will need some basis for this choice. Even a qualitative discussion would help.
The scalability breakdown in Section 4 isolates SWE, RE, and MPI communication, but never the coupling overhead itself (halo transfer, four-case logic). Given that the paper's main claim includes efficiency, it would be worth confirming that this overhead is indeed negligible.
Figure 3 needs a clearer caption. The meaning of the tick mark spacing and the shaded tsub interval is not self-evident from the figure or the surrounding text.
Line 179: the y-range of the inlet boundary in Case 2 is not specified.
Section 2.3 should state explicitly whether the model requires a structured Cartesian grid. This is relevant context for the domain decomposition limitations discussed in Section 4.
The finding in Case 5 that full SWE produces systematically higher ponding depths than kinematic or diffusive wave approximations (lines 247–252) is more interesting than the current discussion suggests. This has direct relevance for flood hazard applications and deserves more than a brief mention. Quantifying the discrepancy and commenting on the conditions under which it is largest would strengthen the paper.
Line 246: "Fig. 14(b)" should be "Fig. 15(b)". Figures 15(a) and 19 have overlapping lines/markers that are hard to distinguish; different line styles or marker sizes would help.
Line 125: "that" should be "than".
References
Casulli, V.: Int. J. Numer. Meth. Fl., 60, 391–408, doi:10.1002/fld.1896, 2009.
Casulli, V.: Int. J. Numer. Meth. Fl., 85, 449–464, doi:10.1002/fld.4389, 2017.
Chourdakis, G., et al.: Open Res. Europe, 2, 51, doi:10.12688/openreseurope.14445.2, 2022.
Coon, E. T., et al.: Adv. Water Resour., 144, 103701, doi:10.1016/j.advwatres.2020.103701, 2020.
Furman, A.: Vadose Zone J., 7, 741–756, doi:10.2136/vzj2007.0065, 2008.
Hutton, E. W. H., Piper, M. D., and Tucker, G. E.: J. Open Source Softw., 5, 2317, doi:10.21105/joss.02317, 2020.
Peckham, S. D., Hutton, E. W., and Norris, B.: Comput. Geosci., 53, 3–12, doi:10.1016/j.cageo.2012.04.002, 2013.
Tubini, N. and Rigon, R.: Geosci. Model Dev., 15, 75–104, doi:10.5194/gmd-15-75-2022, 2022.
Citation: https://doi.org/10.5194/egusphere-2025-4246-RC2 - AC2: 'Reply on RC2', Zhi Li, 06 Apr 2026
Interactive discussion
Status: closed
-
RC1: 'Comment on egusphere-2025-4246', Anonymous Referee #1, 24 Nov 2025
This manuscript presents SERGHEI-SWE-RE, a coupled surface–subsurface model that links a fully dynamic 2D shallow-water solver with a 3D Richards solver. Overall, the work is solid and clearly within the scope of GMD. The model is well described, the benchmark cases are appropriate, and the scaling analysis across several HPC systems is particularly useful. The asynchronous coupling strategy is also a strength of the paper, and the manuscript is generally clear and well organized. I find the manuscript suitable for publication after a set of minor revisions. My comments below are intended to improve clarity and strengthen the presentation.
Major Comments:
The sequential coupling strategy is well explained, but it would be helpful to briefly discuss its potential limitations. In particular, readers would benefit from understanding situations in which asynchronous coupling may introduce small inaccuracies, or cases where a fully coupled approach might be more appropriate. A short clarification in Section 2.2 or in the Discussion would be sufficient. For example, in Case 5 the authors test two maximum RE time steps, which affects the overall runtime almost linearly. While I am not requesting additional experiments, it would strengthen the paper to comment on how sensitive the results are to the chosen coupling time window. How does varying the subsurface time step influence accuracy and computational cost? A brief discussion of this trade-off would help contextualize the asynchronous approach.
One overall question the scalability discussion does not fully address is the performance impact introduced by the coupling itself. In other words, how much does the integration of the SWE and RE solvers degrade performance compared to running the two components independently? A brief quantitative or qualitative assessment of this overhead, whether due to synchronization, data exchange, or load imbalance, would help readers better understand the true cost of coupling and the efficiency of the current implementation.
Rainfall is always applied to the SWE, even when no ponding exists. Since this is a modeling choice and may not hold in all situations, it would be good to clarify the limitations of this assumption. The Case 4 comparison shows the approach works well there, but a sentence noting scenarios where this may be less appropriate would be helpful.
Specific and Minor Comments:
Line 125, “that” to “than”
Figure 3 and the accompanying explanation are hard to follow. In the figure, if the intervals are meant to represent the time steps Δt, their lengths should be consistent. It is also not clear what the shaded/boxed “tsub” period represents. Does this interval mark the window during which the SWE and RE solvers exchange information? If so, this should be stated more explicitly in both the figure and the text.
In Section 2.3, the description of domain decomposition suggests that the model uses a structured grid without regional refinement. If this is indeed the case, it would be helpful to state this explicitly. Clarifying this will make it easier for readers to understand the limitations in the current domain decomposition strategy that you discuss later in the manuscript.
Line 174-176. The author should either report the results even in the supplement, or they should not use such a statement to support their conclusion.
Line 179, what is the range of y for the inlet boundary?
Line 246 “Fig. 14(b)” -> “Fig. 15(b)”
Lines 247–248: The explanation provided is reasonable, but I suggest elaborating a bit more here. Since simplified governing equations such as the kinematic and diffusive wave formulations are widely used in rainfall–runoff modeling, the distinction you highlight is important. This result that full SWE produces systematically larger ponding depths compared to the simplified approaches deserves to be emphasized more clearly, as it is likely to be of broad interest to the hydrologic modeling community.
For Figures 15 and 19, several curves overlap almost exactly, which makes them hard to distinguish. In Fig. 15, the blue and red lines lie directly on top of each other; using different line thicknesses (or slightly different styles) would make the overlap clearer. Similarly, in Fig. 19, the black circular markers are difficult to see because they coincide with another line. Adjusting line weights or marker sizes would improve readability in both figures.
Citation: https://doi.org/10.5194/egusphere-2025-4246-RC1 - AC1: 'Reply on RC1', Zhi Li, 22 Dec 2025
-
RC2: 'Comment on egusphere-2025-4246', Anonymous Referee #2, 21 Mar 2026
The paper presents SERGHEI-SWE-RE, a coupled surface–subsurface model built on the Kokkos framework that targets performance portability across CPUs and GPUs. The work is clearly within GMD's scope, and I found the asynchronous coupling strategy and the multi-platform scalability analysis to be the most interesting contributions. The benchmark suite covers a reasonable range of conditions, and the comparisons with ParFlow, CATHY, HGS and others are generally convincing. That said, I have a few concerns I would like the authors to address. Two of them are substantive and touch on aspects of the model that I think are currently under-described; the third is a straightforward request for numbers that the authors say they have but chose not to show. None of these require new simulations. My recommendation is minor revision, but I want to be clear that I expect concrete answers rather than a restatement of what is already in the paper.
Major Comments:
My main technical concern is with Section 2.2. The four-case boundary condition scheme is described clearly enough, but the manuscript is completely silent on what happens numerically when the system moves between these cases, for instance when a dry cell receives exfiltration and becomes wet, or when the top subsurface layer saturates and the pressure head approaches the surface elevation. These are precisely the transitions that are known to cause convergence problems in this class of models. Coon et al. (2020) document this explicitly for the ATS model, noting that convergence degrades near the dry-to-ponding transition, and Furman (2008) gives a broader review of the numerical difficulties associated with boundary condition switching. The present manuscript does not acknowledge any of this.
I should also point out that more elegant solutions to this problem exist in the literature and are not cited. Casulli (2017) derived a coupled surface–subsurface formulation in which the governing equations for both domains are solved simultaneously using the NCZ algorithm, the same algorithm on which SERGHEI-RE is based, entirely avoiding the switching problem. Earlier, Casulli (2009) introduced the idea of an additional computational node at the soil surface governed by the equation state H(ψ), which allows rainfall to be prescribed as a permanent Neumann boundary condition regardless of ponding state. Tubini and Rigon (2022) implemented exactly this approach in WHETGEO-1D and showed that it handles both infiltration-excess and saturation-excess ponding without any case logic or numerical oscillations. I am not asking the authors to rewrite their code. I am asking them to engage with this literature and explain, in the paper, whether their four-case scheme has been tested under conditions of rapid state transitions and what its behaviour is. If fully saturated conditions with return flow are never actually reached in any of the five benchmark cases, that should be stated explicitly, because right now the reader cannot tell.
The introduction lists modularisation and flexibility as key advantages of SERGHEI-SWE-RE. I agree the internal modularisation between the SWE and RE solvers is well demonstrated, but the paper says nothing about whether the model can be coupled with external codes. This matters because the community has largely converged on standard interfaces for this purpose, the Basic Model Interface (BMI; Hutton et al., 2020) being the most widely adopted, and models that do not address interoperability are increasingly difficult to integrate into larger workflows. WHETGEO (Tubini and Rigon, 2022) does this through OMS3. The Next-Generation Water Resources Modeling Framework uses BMI natively. The paper should state clearly whether SERGHEI-SWE-RE supports any such interface, and if not, acknowledge this as a current limitation rather than implying that flexibility is already there.
A related point: the consistent domain decomposition (same MPI partition for SWE and RE) is already shown in Section 4 to be a bottleneck at high GPU counts. This is an inherent consequence of tying the two solvers' parallelisation together, and it would also complicate any future attempt to expose an external coupling interface. The authors discuss this briefly but do not connect it to the broader architectural question of whether the design supports or hinders future extension. A paragraph addressing this would strengthen the Discussion considerably.
Lines 173–176 state that "further quantitative examination (not shown) confirms" that the C-property is satisfied. This is not acceptable in a model description paper submitted to GMD. The C-property is a fundamental numerical requirement and its verification should be documented with actual numbers, the maximum absolute deviation from the initial condition is the standard metric and takes one line to report. The authors clearly have these numbers. They should include them, either in the text or in a supplementary table. Reviewer 1 made the same point and I fully agree.
Specific and Minor Comments:
The rainfall routing assumption (Section 2.2, lines 140–143) is validated only against Case 4, which was designed with no ponding. This is not a meaningful test of the assumption. The authors should acknowledge the scenarios where routing rain through the surface first introduces a numerical lag, the obvious case being a very dry soil with high infiltration capacity.
On the asynchronous coupling: Case 5 tests two values of Δtsub,max and both give satisfactory results for that specific benchmark. But the paper offers no guidance on how to choose this parameter in general. Users running scenarios with fast-moving wetting fronts or rapidly alternating wet/dry conditions will need some basis for this choice. Even a qualitative discussion would help.
The scalability breakdown in Section 4 isolates SWE, RE, and MPI communication, but never the coupling overhead itself (halo transfer, four-case logic). Given that the paper's main claim includes efficiency, it would be worth confirming that this overhead is indeed negligible.
Figure 3 needs a clearer caption. The meaning of the tick mark spacing and the shaded tsub interval is not self-evident from the figure or the surrounding text.
Line 179: the y-range of the inlet boundary in Case 2 is not specified.
Section 2.3 should state explicitly whether the model requires a structured Cartesian grid. This is relevant context for the domain decomposition limitations discussed in Section 4.
The finding in Case 5 that full SWE produces systematically higher ponding depths than kinematic or diffusive wave approximations (lines 247–252) is more interesting than the current discussion suggests. This has direct relevance for flood hazard applications and deserves more than a brief mention. Quantifying the discrepancy and commenting on the conditions under which it is largest would strengthen the paper.
Line 246: "Fig. 14(b)" should be "Fig. 15(b)". Figures 15(a) and 19 have overlapping lines/markers that are hard to distinguish; different line styles or marker sizes would help.
Line 125: "that" should be "than".
References
Casulli, V.: Int. J. Numer. Meth. Fl., 60, 391–408, doi:10.1002/fld.1896, 2009.
Casulli, V.: Int. J. Numer. Meth. Fl., 85, 449–464, doi:10.1002/fld.4389, 2017.
Chourdakis, G., et al.: Open Res. Europe, 2, 51, doi:10.12688/openreseurope.14445.2, 2022.
Coon, E. T., et al.: Adv. Water Resour., 144, 103701, doi:10.1016/j.advwatres.2020.103701, 2020.
Furman, A.: Vadose Zone J., 7, 741–756, doi:10.2136/vzj2007.0065, 2008.
Hutton, E. W. H., Piper, M. D., and Tucker, G. E.: J. Open Source Softw., 5, 2317, doi:10.21105/joss.02317, 2020.
Peckham, S. D., Hutton, E. W., and Norris, B.: Comput. Geosci., 53, 3–12, doi:10.1016/j.cageo.2012.04.002, 2013.
Tubini, N. and Rigon, R.: Geosci. Model Dev., 15, 75–104, doi:10.5194/gmd-15-75-2022, 2022.
Citation: https://doi.org/10.5194/egusphere-2025-4246-RC2 - AC2: 'Reply on RC2', Zhi Li, 06 Apr 2026
Peer review completion
Journal article(s) based on this preprint
Viewed
| HTML | XML | Total | BibTeX | EndNote | |
|---|---|---|---|---|---|
| 480 | 311 | 39 | 830 | 36 | 39 |
- HTML: 480
- PDF: 311
- XML: 39
- Total: 830
- BibTeX: 36
- EndNote: 39
Viewed (geographical distribution)
| Country | # | Views | % |
|---|
| Total: | 0 |
| HTML: | 0 |
| PDF: | 0 |
| XML: | 0 |
- 1
Na Zheng
Gregor Rickert
Mario Morales-Hernández
Ilhan Özgen-Xian
Daniel Caviedes-Voullième
The requested preprint has a corresponding peer-reviewed final revised paper. You are encouraged to refer to the final revised version.
- Preprint
(7116 KB) - Metadata XML
This manuscript presents SERGHEI-SWE-RE, a coupled surface–subsurface model that links a fully dynamic 2D shallow-water solver with a 3D Richards solver. Overall, the work is solid and clearly within the scope of GMD. The model is well described, the benchmark cases are appropriate, and the scaling analysis across several HPC systems is particularly useful. The asynchronous coupling strategy is also a strength of the paper, and the manuscript is generally clear and well organized. I find the manuscript suitable for publication after a set of minor revisions. My comments below are intended to improve clarity and strengthen the presentation.
Major Comments:
The sequential coupling strategy is well explained, but it would be helpful to briefly discuss its potential limitations. In particular, readers would benefit from understanding situations in which asynchronous coupling may introduce small inaccuracies, or cases where a fully coupled approach might be more appropriate. A short clarification in Section 2.2 or in the Discussion would be sufficient. For example, in Case 5 the authors test two maximum RE time steps, which affects the overall runtime almost linearly. While I am not requesting additional experiments, it would strengthen the paper to comment on how sensitive the results are to the chosen coupling time window. How does varying the subsurface time step influence accuracy and computational cost? A brief discussion of this trade-off would help contextualize the asynchronous approach.
One overall question the scalability discussion does not fully address is the performance impact introduced by the coupling itself. In other words, how much does the integration of the SWE and RE solvers degrade performance compared to running the two components independently? A brief quantitative or qualitative assessment of this overhead, whether due to synchronization, data exchange, or load imbalance, would help readers better understand the true cost of coupling and the efficiency of the current implementation.
Rainfall is always applied to the SWE, even when no ponding exists. Since this is a modeling choice and may not hold in all situations, it would be good to clarify the limitations of this assumption. The Case 4 comparison shows the approach works well there, but a sentence noting scenarios where this may be less appropriate would be helpful.
Specific and Minor Comments:
Line 125, “that” to “than”
Figure 3 and the accompanying explanation are hard to follow. In the figure, if the intervals are meant to represent the time steps Δt, their lengths should be consistent. It is also not clear what the shaded/boxed “tsub” period represents. Does this interval mark the window during which the SWE and RE solvers exchange information? If so, this should be stated more explicitly in both the figure and the text.
In Section 2.3, the description of domain decomposition suggests that the model uses a structured grid without regional refinement. If this is indeed the case, it would be helpful to state this explicitly. Clarifying this will make it easier for readers to understand the limitations in the current domain decomposition strategy that you discuss later in the manuscript.
Line 174-176. The author should either report the results even in the supplement, or they should not use such a statement to support their conclusion.
Line 179, what is the range of y for the inlet boundary?
Line 246 “Fig. 14(b)” -> “Fig. 15(b)”
Lines 247–248: The explanation provided is reasonable, but I suggest elaborating a bit more here. Since simplified governing equations such as the kinematic and diffusive wave formulations are widely used in rainfall–runoff modeling, the distinction you highlight is important. This result that full SWE produces systematically larger ponding depths compared to the simplified approaches deserves to be emphasized more clearly, as it is likely to be of broad interest to the hydrologic modeling community.
For Figures 15 and 19, several curves overlap almost exactly, which makes them hard to distinguish. In Fig. 15, the blue and red lines lie directly on top of each other; using different line thicknesses (or slightly different styles) would make the overlap clearer. Similarly, in Fig. 19, the black circular markers are difficult to see because they coincide with another line. Adjusting line weights or marker sizes would improve readability in both figures.