Thank you, it is helpful for the ocean model development community to see your success with a split explicit method. This paper takes great care in the detailed presentation of the numerical methods for the new split explicit solver in FESOM2.

What order of convergence do you expect for your time-stepping method? Please test the convergence rate in time. The test case would need to have both slow and fast dynamics to be a useful test, but you could test them first separately.

You could test JUST the barotropic part with the Kelvin Wave or Inertial-Gravity wave tests here in Bishnu et al. (2024). Just use redundant layers in the vertical.  Then you can compare against an exact solution for the convergence test.

You could then use the baroclinic channel test case from Ilicak et al. (2012) and Petersen et al (2015), refine in time, and compute convergence by comparing to a short-timestep case. The duration of the simulation would need to be short (30 minutes, say) and sufficiently laminar that the results do not differ due to small differences in the turbulent flow. This test case definitely has baroclinic dynamics. It is not designed for barotropic dynamics, but should include some surface gravity waves as the SSH and barotropic eastward flow adjust to a geostrophic balance. Thanks to Ange Ishimwe (Université catholique de Louvain) for pointing out this test case during his talk at AGU Ocean Sciences. He was able to show second-order convergence for his split explicit scheme using this method in his recent paper published in December, in Ishimwe et al. (2023) figure 6. This would be a good paper to reference for comparison in the current work.

Smaller comments:

1. You subcycle the external mode to exactly one baroclinic time step, and do not use a filter, as explained in lines 30-34. I think this is important enough to put in the abstract, as it is a potential 2x speed-up for the barotropic stage compared to models that subcycle to n+2
2. On equation 8, first line, I believe the sign of gH grad eta should be positive.
3. The surface flux W was dropped in equation 8. It would be good to just comment that W was added to the baroclinic (eqn 14), and not the barotropic mode, and your reasoning for that. Because you have a function at the top with the delta in equation 14, it makes sense that this is a baroclinic addition.
4. Fig 3 caption m2 needs superscript

I appreciate your description of the reasoning behind your choices. For example, the discussion of whether to add the bottom drag term to the barotropic dynamics at line 95, and the discussion on abandoning filtering at line 110. I have considered these exact issues and it is good to hear your thoughts on these.

Thank you,

Mark R. Petersen, Los Alamos National Laboratory

Bishnu, S., Petersen, M.R., Quaife, B., Schoonover, J.A. A Verification Suite of Test Cases for the Barotropic Solver of Ocean Models. Submitted to JAMES. https://essopenarchive.org/users/566154/articles/612943-a-verification-suite-of-test-cases-for-the-barotropic-solver-of-ocean-models

Ishimwe, A.P., Deleersnijder, E., Legat, V., Lambrechts, J., 2023. A split-explicit second order Runge–Kutta method for solving 3D hydrodynamic equations. Ocean Modelling 186, 102273. https://doi.org/10.1016/j.ocemod.2023.102273

Ilicak, A.J. Adcroft, S.M. Griffies, R.W. Hallberg. Spurious dianeutral mixing and the role of momentum closure  Ocean Modell., 45 (2012), pp. 37-58, 10.1016/j.ocemod.2011.10.003

Petersen, M.R., D.W. Jacobsen, T.D. Ringler, M.W. Hecht, M.E. Maltrud, Evaluation of the arbitrary Lagrangian–Eulerian vertical coordinate method in the MPAS-Ocean model, Ocean Modelling, Volume 86, February 2015, Pages 93-113, ISSN 1463-5003, http://dx.doi.org/10.1016/j.ocemod.2014.12.004.

Review

The authors describe the implementation of a split-explicit algorithm for the FESOM2 ocean model, using for time integration of the external mode the modified FB scheme or the modified AB3-AM4 scheme proposed by Demange et al. The implementation is described for the z-star and z-tilde vertical coordinates. Numerical integration tests are carried out for the z-star split-explicit model, and compared with the current semi-implicit algorithm of FESOM2. These tests show numerical solutions that are qualitatively close to those obtained with the semi-implicit version. They also show better scalability, particularly for highly parallelized workflows.

General comments :

The paper addresses the interesting question of how to construct a split-explicit algorithm based on the FESOM time-stepping scheme, which has the peculiarity of staggering the prognostic variables in time. Original questions arise about how to phase the time integration of 3D variables, momentum and tracers, and the time integration of the barotropic mode. However, it seems to me that the paper needs to be clarified on several points and requires major revisions.

Section 2 details the proposed split-explicit algorithm. Some clarifications seem necessary to me concerning (i) the choice made for the order of integration of the variables Ubar and eta with the modified FB scheme and the modified AB3-AM4 scheme and the implications of this choice, (ii) the initialisation of the variable Ubar for the barotropic integration, (iii) the correction of the 3d variables after completion of the barotropic integration.

Section 3 details the amplitude and phase errors of the modified FB scheme, the modified AB3-AM4 scheme and the semi-implicit FESOM scheme in the context of the SW equations. This analysis, although largely based on published results by Demange et al., could nicely illustrate the contrast between the errors produced by explicit schemes and the semi-implicit FESOM scheme. The presentation and the figure should however be improved. It seems to me that sections 2 and 3 should then be swapped, so that the SW-based results of the current section 3 introduce/motivate the implementation of the split-explicit algorithm of current section 2.

Section 4 reports the results of numerical experiments. It is somewhat disturbing that the solutions for the idealised case show little difference between the split-explicit algorithm and the semi-implicit algorithm. This is probably due to the fact that the dynamics of this idealised case is essentially baroclinic and little affected by the representation of the external dynamics. It would have been desirable to consider an idealised case where barotropic dynamics play a more important role (for example, think of an idealised case of internal tide generation on a bathymetry by barotropic current oscillation). However, the results show a certain viability of the split-explicit algorithm and a gain in efficiency compared to the semi-implicit algorithm.

Specific comments :

l37 : I suppose ‘temporal interpolations’ should be replaced by ‘temporal integrations’?

l47 : A schematic describing the time-staggering of variables and the structure of the algorithm could be referred to here, upstream of the equations, to help the reader.

l53 : Should the variable h in the expression for the pressure gradient be indexed with k?

l62 : Should the transport U be indexed with k?

l100 : It seems to me that the abbreviation SE isn't very well chosen. It would be clearer to use another which suggests that the basic time-stepping is the modified FB. For the same reason,

l101-102 : Here, with the SE scheme, the choice is made to integrate U_bar before eta. With the SESM scheme, the choice will be reversed. What are the reasons for these different choices? What are the implications? In the SE case, this choice is associated with a semi-implicit discretisation of the Coriolis term, which probably requires specific numerical processing ; if so, these should be explained. On the other hand, still in the SE case, this choice leads to the expression of the mean barotropic transport (12), which would be different with a reverse order of integration.

l103 : The quoted value of theta has been obtained under rather special conditions: it is the value that, for a constant stratification of N=10-2 s-1, a water column of H=4000 m and a splitting ratio of 20, allows the stabilisation of split-explicit algorithms whose stability is limited by large-scale barotropic waves. It is not clear that the algorithm proposed in this paper is constrained by these waves and that this value is relevant here. This is plausible, however, because the FESOM scheme is a variant of the FB scheme, which is used by Demange et al. in their study.

l110 : It might be clearer to change ‘at the end of the barotropic step’ into ‘during the barotropic step’.

l114 : The reference to Shchepetkin and McWilliams on the line 112 could be moved to just after the ‘traditional notation’ line 114.

l120 : It is mentioned that it is not clear how to initialise Ubar (due, I suppose, to the time-staggering between U and Ubar), and it is noted that "We return to this topic below". But it's not clear to me where in the paper this is explained.

l125 : In this expression we see the correction by the mean barotropic transport <> of the 3D transports U_k^{n+1/2}, which is then used to transport tracers. There doesn't seem to be any place in your implementation for the other transport correction that classically appears in split-explicit algorithms based on synchronous schemes, i.e the correction by the barotropic transport of the 3D transports U_k^{n+1} (see for example p 384 in Shchepetkin and McWilliams). I'm a little embarrassed by the lack of an equivalent to this correction. Isn't there a cog missing to ensure consistency between the 3D and 2D integrations?

l156 : Should ‘temporal interpolations’ be replaced by ‘temporal integrations’?

l157-159 : It seems to me that this sentence is confusing. What is done in this section is not really the analysis of the external mode in the context of split-explicit algorithms (which is done, for example, in Demange et al. section 4.2). Rather, it is the analysis of the modified FB, modified AB3-AM4 and semi-implicit FESOM schemes in the context of SW equations. 

l190 : The solutions of the continuous problem are e^{+i c} and e^{-i c}.

l 214 : The sentence ‘This CFL number …mesh size.’ is unnecessary here because the analysis is exact in space and does not depend on Delta x.

fig. 1: For the two panels on the left, why not show the amplitude error of the explicit schemes over larger ranges on each of the two axes? This would show the damping for large values of c, and the loss of stability of the schemes. This would make it a little easier to compare the amplitude error of the explicit schemes with that of the semi-implicit scheme. For the two panels on the right, the phase error could be plotted as the ratio of the discrete phase velocity to the exact phase velocity. why not use the same horizontal axis as for the left panels, with c values varying between 0 and 30? This would make the figure easier to read.

fig. 3 : the caption refers to diffusivities that are not traced.