Numerical solution of two-layer shallow water equations using mode splitting method

In the numerical models that use iterative methods to solve the momentum equations by applying the rigid-lid approximation, the number of iterations increases for high resolution, therefore processing time increases. An alternative method is applying a free surface and splitting equations to barotropic and baroclinic modes. The surface gravity waves that are faster than slow moving internal gravity waves; impose a limitation on the time steps with the CFL condition. Thus, mode splitting method is computationally efficient to handle the multiple time steps separating the barotropic and baroclinic mode equations. In this method, the barotropic mode equations are solved at small time steps consistent with the fast surface gravity wave speeds and the baroclinic mode equations are solved at larger time steps consistent with the slow internal gravity wave speeds. This is used in most of the ocean circulation models and is an unavoidable choice to high resolution models.
In this study, we considered the shallow water equations for two-layer basin with vorticity-divergence formulation using mode splitting method by a small time step of barotropic mode within a larger time step of baroclinic mode. The primary systems of equations that contain both upper and lower layer variables, are rewritten in terms of new (barotropic and baroclinic) variables without any variation or more approximation of primary systems. This procedure can be extended to multi-layer systems so that primary N-layer system of equations is changed to 1 system of barotropic mode equations and N-1 systems of baroclinic mode equations coupled together.
For numerical experiments, a fully baroclinic (non-barotropic) initial condition is considered in a constant depth rectangular domain with 64, 128 and 256 grid points in each direction and periodic boundaries. For spatial differencing, second order centered scheme with low computational cost and fourth-order compact scheme with high computational cost are used. For time integration, a semi-implicit descretization based on leapfrog scheme is implemented with the Robert-Asselin time filter for both barotropic and baroclinic systems of equations, similarly.
Mode splitting method may presents numerical instabilities on the larger baroclinic time steps, in spite of time step limitation based on CFL condition coming from each system of barotropic and baroclinic mode equations taken individually. Here, it is controlled by increasing the coefficient of time filter to some extent.
First, we solve the baroclinic mode equations to derive all baroclinic variables that are necessary to solve barotropic mode equations during a baroclinic time step. In this case, these variables can be taken to be constant up to the next baroclinic time level or determined by time interpolation between two successive baroclinic time levels.
To assess the performance of the numerical method, relative error of energy conservation is calculated. Results show that for the ratio of baroclinic time step to up to 20 times of that of barotropic time step, time evolution of the barotropic and baroclinic variables have appropriate correspondence to the basic state, in which the barotropic mode has the same time step as the baroclinic mode. When this ratio increases, the differences of errors from basic state are presented more clearly. These errors are increased on fourth order compact method insofar as it leads to numerical instabilities so the time filter coefficient had to be increased, while second order scheme is not sensitive and stays stable with small coefficient. Moreover, taking constant for baroclinic variables to solve barotropic mode equations makes the solution on fourth order compact scheme for large baroclinic time step unstable, but on the other hand time interpolation provides more stable condition and has a good performance almost on both spatial schemes.