Numerical solution of two-layer shallow water equations in terms of barotropic and baroclinic variables using fourth-order compact method

Two-dimensional shallow water equations are commonly used to study the dynamics of large-scale flows in the atmosphere and oceans that are nearly horizontal. Over the past years, the numerical solution of multi-layer shallow water systems has been widely researched. In stratified geophysical flows, two-layer shallow water equations are proper models for the simulation of certain phenomena in the atmosphere and oceans. In this model, the fluid is assumed to be composed of two shallow layers of immiscible fluids in which the superposed layers differ in velocity and density in a two-dimensional domain. The constant density of upper layer is less than the density of the lower one.
For the solution of shallow water equations, high order compact finite difference schemes have been widely used owing to their good accuracy compared with standard finite difference schemes for a large range of wave numbers and a low numerical diffusion with small dispersion errors.
In this research, fourth-order compact finite-difference method was used for spatial differencing of f-plane two-layer shallow-water equations in the vorticity-divergence formulation for a rectangular domain with periodic boundaries; the results were compared to those of conventional second order finite difference method. For time integration, a three-level semi-implicit formulation was applied with the Robert-Asselin time filter which prevents the numerical instability caused by the computational mode of the three-time-level scheme. The equations were derived in terms of barotropic and baroclinic variables such that they were split into two coupled systems (barotropic and baroclinic systems) consisting of all variables of both upper and lower layers in each system. These systems are different than standard two-layer shallow water systems.
A perturbed unstable zonal jet, in a two-layer shallow-water flow, was considered for initial value problem with baroclinic instabilities. The two-layer baroclinic initial values were extracted from a one-layer initial condition such that the potential vorticity of one-layer initial condition was equal to the upper layer potential vorticity; other initial values of the two-layer shallow water equations were determined by taking the definition of baroclinic initial condition, in the situation where the initial current in the upper layer is opposite to one in the lower layer such that the initial barotropic variables are all zero. Over time, the initial state broke up into smaller barotropic and baroclinic vortices.
To control the numerical nonlinear instability due to the interaction of nonlinear terms of the vorticity equations and subsequently produce the aliasing errors, the barotropic and baroclinic planetary vorticities were explicitly damped by adding the hyperdiffusion operator acting on the same vorticity equations.
There exist myriad methods for assessing the accuracy of numerical schemes such as the conservation of total mass and energy. The ability of each simulation to conserve the total energy on the total layer-depth (i.e. energy error) and mass between isolevels of the potential vorticity (i.e. mass error) within each layer, was measured in order to assess the numerical accuracy. Results showed that this model meets the conservation laws of mass and energy in the test problem considered here. The fourth-order compact finite-difference method entailed less energy error than the second order finite difference method during the time integration of two-layer shallow water equations; additionally, it is more accurate regarding the conservation of mass for a long time integration.