The discontinuous Galerkin coastal and estuarine modelling system (DGCEMS v1.0.0): a three-dimensional, nonsplit-mode, implicit–explicit Runge–Kutta hydrostatic model

Chen, Zereng; Zhang, Qinghe; Ran, Guoquan; Nie, Yang

doi:https://doi.org/10.5194/gmd-18-7199-2025

Articles | Volume 18, issue 19

https://doi.org/10.5194/gmd-18-7199-2025

© Author(s) 2025. This work is distributed under
the Creative Commons Attribution 4.0 License.

https://doi.org/10.5194/gmd-18-7199-2025

© Author(s) 2025. This work is distributed under
the Creative Commons Attribution 4.0 License.

Articles | Volume 18, issue 19

Development and technical paper

|

14 Oct 2025

Development and technical paper |

| 14 Oct 2025

The discontinuous Galerkin coastal and estuarine modelling system (DGCEMS v1.0.0): a three-dimensional, nonsplit-mode, implicit–explicit Runge–Kutta hydrostatic model

Zereng Chen, Qinghe Zhang, Guoquan Ran, and Yang Nie

Download

Final revised paper (published on 14 Oct 2025)
Preprint (discussion started on 12 Feb 2025)

Interactive discussion

Status: closed

RC1:
'Comment on gmd-2024-240', Anonymous Referee #1, 28 Feb 2025

In this paper, the authors describe their development of a novel 3D coastal and estuarine modelling system called DGCEMS based on the nodal discontinuous Galerkin method. Through some tests, it has been demonstrated that the model has second-order convergence, low spurious mixing, and capability to simulate salt-freshwater interactions in the presence of wetting and drying boundaries. The subject of the paper is well presented, and definitely of interest to the modeling community. I’d recommend the paper for publication, after addressing the following comments.
1.In the governing equations, no specific vertical stratification was given, but Figure 1 shows 2 layers, while 10(Line 230) and 20 layers(Line 254) were used in Section 3.1 and 3.2, respectively. How is vertical stratification determined?
2.When presenting the model algorithm, it is necessary to highlight the innovative points of the solution, which can help readers better understand.
3.In model validation, the sources of analytical and experimental solutions should be provided first. In other words, from cases in Section 3.1 to 3.4, which ones are referenced from others and which ones are used for the first time, there should be more specific explanations.

Citation: https://doi.org/10.5194/gmd-2024-240-RC1
- AC1: 'Reply on RC1', Qinghe Zhang, 28 Feb 2025
  
  Thank you for the suggestions and evaluation.
  Answer1: The vertical layering of the model can theoretically be arbitrary. The more layers there are, the higher the accuracy of vertical flow velocity will be, but it will not increase the convergence order of the solution and will also increase the calculation time of the model. The two layers shown in Figure 1 are intended to clearly express the spatial distribution of interpolation nodes in the vertical direction. We will add necessary explanations in the revised version.
  Answer2: Thank you. We will emphasize the significance of using the wet dry treatment and limiters in the mode nonsplit model algorithm. The 3D limiters are always applied after achieving the physical field to prevent pathological solutions, and then we can obtain the vertically averaged physical field. The WD treatment is carried out after obtaining the vertically averaged physical field to ensure the conservation of water elevation and depth-average momentum.
  Answer3: Thanks for the suggestion. The artificial analytical solution in Section 3.1 is a re-derivation of the analytical expression in the sigma coordinate system based on Kärnä et al. (2018). The case in Section 3.4 is inspired by the examples used in Chen et al. (2022) and conducted research using our own designed grid.
  
  Citation: https://doi.org/10.5194/gmd-2024-240-AC1
CC1:
'Comment on gmd-2024-240', Ange Ishimwe, 22 Apr 2025

The paper presents an interesting coastal and estuarine modelling system based on discontinuous Galerkin methods. The model is three-dimensional, hydrostatic, and employs a mode-nonsplit, implicit-explicit (IMEX) Runge–Kutta time integration. It is built using a quadrature-free nodal formulation and is claimed to be second-order accurate in both space and time. The mathematical formulation and notation are clearly presented, and the numerical methods are well described.
My main concern lies in the claim of second-order accuracy in time. The IMEX approach treats the vertical diffusion term implicitly, while all other terms are handled explicitly. When these two components are considered separately, the explicit scheme corresponds to a second-order Runge–Kutta method:
y^1 = y^n + dt f(y^n)
y^n+1 = y^n + dt/2 ( f(y^n) + f(y^1))
In contrast, the implicit part used for vertical diffusion is essentially a first-order implicit Euler method:
y^n+1 = y^n + dt f(y^n+1)
This means that when vertical diffusion is active, even if it is not dominant, the scheme loses second-order temporal accuracy and becomes effectively first-order in time for those terms.
Unfortunately, the numerical experiments do not adequately support the claim of second-order accuracy in time. The first test, based on a manufactured solution, excludes diffusion and only assesses spatial convergence. No temporal convergence is shown, so the second-order accuracy in time is not validated in this case. Moreover, since diffusion is excluded, the test cannot address the concern regarding the implicit treatment of vertical diffusion. The second test case, a lock-exchange with constant viscosity, primarily aims to evaluate spurious horizontal mixing and shows coherent evolution of the reference potential energy. However, it does not assess the temporal accuracy of the scheme either.
Could the authors clarify how second-order accuracy in time is ensured when vertical diffusion is included? Maybe additional numerical evidence demonstrating this, or further discussion on the temporal discretization strategy, would greatly strengthen the manuscript.

Citation: https://doi.org/10.5194/gmd-2024-240-CC1
- AC2:
  'Reply on CC1', Qinghe Zhang, 23 Apr 2025
  
  Thank you for the thoughtful comments. We now respond to the concern regarding the temporal accuracy of our IMEX time discretization, particularly in the presence of vertical diffusion treated implicitly.
  As correctly pointed out, the vertical diffusion term is discretized using a backward Euler scheme, which is first-order accurate in time. However, our time integration strategy is based on a class of additive Runge‒Kutta methods known as IMEXRK schemes, which allow for a second-order explicit scheme to be coupled with a first-order implicit treatment, while still maintaining overall second-order accuracy. In this second-order IMEX framework, the implicit terms do not need to be of the same order as the explicit terms. As long as the overall scheme satisfies the conditions for second-order accuracy, or if the numerical error is dominated by the explicit terms, the method can still achieve second-order accuracy.
  In another word, you may find in Eq. (27) that the implicit step does not use y^n, but y^(2) instead. Here the final physical fields results are calculated by the values from the second prediction step. Besides, the coefficients in each stage must jointly satisfy the Butcher tableau in order for the coupling of the second-order explicit and first-order implicit time discretizations to achieve overall second-order accuracy. The present study follows the second-order IMEX Runge–Kutta scheme used in the Thetis model (Kärnä et al, 2018), which provides justification for the second-order accuracy in time.
  We agree with the reviewer that verifying the temporal order of accuracy through convergence tests is a meaningful effort. However, designing such test cases requires carefully eliminating the influence of spatial discretization errors and involves consideration of multiple factors. We plan to conduct a more in-depth investigation on this in future work.
  
  Citation: https://doi.org/10.5194/gmd-2024-240-AC2
  - CC2: 'Reply on AC2', Ange Ishimwe, 28 Apr 2025
    
    Thank you very much for your thoughtful reply.
    While I appreciate the clarification regarding the IMEX Runge–Kutta framework and understand that it can, in principle, achieve second-order accuracy even when the implicit part alone is first-order, I believe that such an important claim should not be left without some form of justification in the paper.
    Currently, the manuscript does not provide a direct reference that formally proves the second-order accuracy of the specific scheme employed, nor does it clearly state and verify the conditions on the coefficients that would ensure the correct order. I also consulted the Thetis paper (Kärnä et al., 2018) and, unless I overlooked something, a complete demonstration of the full second-order accuracy (including the implicit contribution) is not provided there either.
    I fully understand that setting up a dedicated numerical test to verify temporal convergence with vertical diffusion can be heavy and time-consuming. As an alternative, it could perhaps be helpful to include a brief theoretical analysis in the Appendix, expanding the discrete scheme and comparing it to a Taylor series expansion. This would allow one to verify that the scheme satisfies:
    y^n+1 = y^n + (f_ex^n + f_im^n) dt + (f_ex^n + f_im^n)*(f_ex^n + f_im^n)' dt/2
    Such a theoretical check could be relatively lightweight to include and would fully confirm the second-order accuracy claim without requiring heavy numerical experiments.
    Thank you again for your efforts and for considering this suggestion.
    
    Citation: https://doi.org/10.5194/gmd-2024-240-CC2
    
    AC3: 'Reply on CC2', Qinghe Zhang, 29 Apr 2025
    
    Thank you for your suggestion and understanding. We will include a brief theoretical analysis of the second-order temporal accuracy in the appendix. Details are shown in the supplement.
    
    Citation: https://doi.org/10.5194/gmd-2024-240-AC3
    
    CC3: 'Reply on AC3', Ange Ishimwe, 01 May 2025
    
    Thank you again for your effort in providing a theoretical analysis in the supplementary material.
    However, after reviewing it carefully, I still have some concerns regarding the validity of the second-order accuracy claim. The final formula appears to be incomplete and misses cross terms required to match the full second-order Taylor expansion.
    Specifically, in your expansion you arrive at:
    y^n+1 = y^n + (f_ex^n + f_im^n) dt + ( (f_ex^n + (f_ex^n')) + ( f_im^n + (f_im^n)') ) dt^2/2
    However, this does not match the expected second-order expansion of a solution to y′ = f_ex(y) + f_im(y) , which is:
    y^n+1 = y^n + (f_ex^n + f_im^n) dt + (f_ex^n + f_im^n)*(f_ex^n + f_im^n)' dt^2/2
    
    = y^n + (f_ex^n + f_im^n) dt + ( ( f_ex^n + (f_ex^n)' )+( f_ex^n + (f_im^n)' )+( f_im^n + (f_ex^n)' )+( f_im^n + (f_im^n) ') ) dt^2/2
    These cross terms are essential to correctly capture the interaction between the explicit and implicit parts of the system, and they are missing in your derivation.
    Additionally, in Eq. (B2), you express y^2 = y^n + (f_ex^n ) dt + (f_ex^n )*(f_ex^n )' dt^2/2 + O(dt^3) . However, in Eq. (B5), when evaluating f_im (y^2) , you write:
    f_im (y^2) = f_im(y^n) + ( f_im^n)*( f_im^n)' dt/2 + O(dt^2)
    This is inconsistent. If you substitute the expression for y^2 from B2 into a first-order expansion of f_im(y^2), you should obtain:
    f_im (y^2) = f_im(y^n) + ( f_ex^n)*( f_im^n)' dt + O(dt^2)
    That is, the derivative term should be multiplied by f_ex(y^n), and the factor 1/2 should not appear. Unfortunately, even with this correction, the final formula would remain incorrect.
    In conclusion, even with the revised steps, the derivation does not convincingly establish second-order accuracy in time for the full IMEX scheme. The missing terms are not simply technicalities. They are structurally required to validate the claim.
    That said, I appreciate the value of the work presented in the paper, and I recognize that this specific point is not central to the overall contribution. In light of the current state of the analysis and the absence of a temporal convergence test, it might be more appropriate to soften the second-order accuracy claim in the main text.
    Thank you again for the constructive exchange and for your contributions to the modeling community.
    
    Citation: https://doi.org/10.5194/gmd-2024-240-CC3
    
    AC4: 'Reply on CC3', Qinghe Zhang, 02 May 2025
    
    Thank you for your guidance on this manuscript. We have re-examined the previous derivations. We think that Eq. (24) to (25) represent a complete time step with a step size of dt, while Eq. (26) employs an intermediate step (essentially rewinding to 0.5dt) to achieve better approximation. Thus, the expansion of Equation B5 can correspond to half a time step. Nevertheless, when cross-terms are retained, it is indeed theoretically impossible to rigorously prove that the temporal accuracy reaches second order after incorporating vertical diffusion.
    Since the implicit part is computed after two explicit updates, simultaneous Taylor expansions for proving second-order temporal accuracy may not hold. The Thetis (2018) paper considered that this method achieves second-order accuracy in both time and space, which is why similar descriptions were adopted in our original text. As you noted, although the simulation results exhibit second-order convergence, a strict theoretical analysis of temporal and spatial second-order accuracy will be addressed in future studies. Accordingly, we will soften the second-order accuracy claim in the manuscript.
    Thank you again for the constructive exchange, and it has been immensely valuable and greatly enriched our understanding.
    
    Citation: https://doi.org/10.5194/gmd-2024-240-AC4
RC2:
'Comment on gmd-2024-240', Anonymous Referee #2, 25 Jul 2025

This paper presents a 3D modeling system based on the nodal discontinuous Galerkin (DG) method, utilizing a non-split mode framework and σ-coordinates, with implicit-explicit Runge-Kutta time stepping. The authors conduct numerical tests demonstrating second-order convergence. The study is well-structured and of clear scientific significance. I recommend publication after the authors address the following minor concern:
1. The conclusion section could be strengthened by including a few sentences discussing the potential improvements and possible future extensions or implementations of the proposed method.

Citation: https://doi.org/10.5194/gmd-2024-240-RC2
- AC5: 'Reply on RC2', Qinghe Zhang, 26 Jul 2025
  
  Thank you for your recognition and positive evaluation of this work. We will incorporate future plans for model improvement and potential extensions into the conclusion section. The model will be further optimized and enhanced in terms of parallel computing strategies and the overall solution framework. Future extensions will include additional functionalities and modules, such as wave–current interaction, sediment transport, and biogeochemical processes, with broader applications to realistic estuarine and coastal environments.
  
  Citation: https://doi.org/10.5194/gmd-2024-240-AC5

Peer review completion

AR: Author's response | RR: Referee report | ED: Editor decision | EF: Editorial file upload

AR by Qinghe Zhang on behalf of the Authors (30 Jul 2025) Author's response Author's tracked changes Manuscript

ED: Publish as is (30 Jul 2025) by Simone Marras

AR by Qinghe Zhang on behalf of the Authors (30 Jul 2025) Author's response Manuscript

Short summary

This study presents the development of a novel three-dimensional discontinuous Galerkin coastal and estuarine modelling system named DGCEMS. The model has low spurious mixing and second-order convergence of surface water elevation, horizontal velocity, and the tracer field. It has the capability to simulate salt–freshwater interactions in the presence of wetting and drying boundaries.