## 1 Introduction

In long-term simulations of the geophysical flows, where millions of iterations are involved, numerical accuracy alone cannot ensure an accurate representation of certain essential statistics about the flows. Low-order numerical schemes armed with better conservative properties and dispersive wave representations can sometimes reproduce more realistic dynamics than higher order numerical schemes do (Arakawa1977-og; Arakawa1981-dy). Many authors have been inspired to search for numerical schemes that are not necessarily more accurate, but possess more desirable properties in conservation of key quantities, such as mass, energy, enstrophy, and/or representation of the dispersive wave relations (Heikes1995-ky; Ringler2002-jx; Chen2013-fa; Ringler2010-sm; Gassmann2012-xb; Thuburn2014-ns). It is in the same spirit that we develop, in Part I. of this project, a new numerical scheme that conserves both the total energy and potential enstrophy and possess the optimal dispersive wave relations on unstructured meshes over bounded or unbounded domains. This article constitutes Part II of the project, and its goal is to numerically evaluate the scheme.

A focal question in the evaluations of the new scheme is what advantages it has over conventional numerical schemes that lack some of its properties. More specifically, we would like to understand how the energy and enstrophy conservations, and the representation of the dispersive wave relations affect the long-term dynamics of the simulated flows. To fully address these important questions requires extensive tests using real-world applications, which are out of the reach and scope of the current project, and are left for future endeavors. In this project, we content ourselves with a preliminary study of these questions in a somewhat idealized setting. For the effect of the energy and enstrophy conservations, Arakawa and Lamb (Arakawa1981-dy) make use of an idealized test case with a meridional ridge within a zonal reentrant channel, and shows that, a few weeks out into the simulations, their total energy and potential enstrophy conserving scheme, compared with a non-conserving scheme, is better at eliminating the spurious inverse energy cascade, and leads to a less noisy and more organized wind field. Here, we take a similar approach, and perform, among all the tests, a comparative study using an external numerical model, the MPAS-sw shallow model, and an existing test case (the shallow water standard test case #5 with mountain topography). On a global sphere, the MPAS-sw model operates on exactly the same kind of orthogonal meshes as our model does, and therefore is an ideal choice for comparison. In this study, we not only examine the vorticity dynamics, but also compute and compare the energy and enstrophy spectra for the two models at different points during the simulations. As for the representation of the dispersive wave relations, it is commonly accepted it has profound impact on the geostrophic adjustment process (Arakawa1977-og; Randall1994-vu), but we are not aware of any numerical studies of its effect in long-term real applications. In this project, we perform a second comparative study, comparing our model with the MPAS-sw model in terms of the control of the divergence variable. The rationale for this study is simple: a geostrophically balanced flow should stay close to being non-divergent.

With the focus on the potential impacts of the new scheme in real applications, this numerical study is conducted around the scheme’s properties, rather than around the established test cases, as is usually done in the literature (Chen2013-fa; Ringler2010-sm; Chen2014-cy; Thuburn2015-hi). Specifically, this study considers six properties of the numerical scheme, which are deemed essential for the success of long-term simulations: the accuracy of the individual operators, the accuracy of the whole model, the conservation of key quantities, the control of the divergence variables, the representation of the energy and enstrophy spectra, and the simulation of the barotropic instabilities. Each property is evaluated using well selected or specially designed test cases. The rest of the article is organized as follows. Section 2 recalls the conservative numerical scheme. Section 3 presents the numerical results. Some remarks are provided in Section 4.

## 2 The conservative schemes

The nonlinear rotating shallow water equations read, in the vorticity-divergence form,

(1) |

Here, stands for the fluid top-to-bottom thickness, the horizontal velocity field of the fluids, the relative vorticity, and the divergence.

When treated as a single variable (Salmon2004-tt; Salmon2005-pd; Salmon2007-dm), the mass flux has a Helmholtz decomposition (Girault1986-sr)

(2) |

where and are the streamfunction and the velocity potential respectively, and

the skewed gradient operator. In the classical Helmholtz decomposition,

is already assumed to satisfy the homogeneous Dirichlet boundary condition (Girault1986-sr); to enforce no-flux boundary condition on the flow, one only need to set the normal derivative of the velocity potential to zero on the boundary, and thus the boundary conditions on and are(3a) | ||||

(3b) |

The relation between and the vorticity and divergence can be easily derived,

(4) |

Under the boundary conditions just given in (3) and with a strictly positive thickness field , the system (4) is a coupled, self-adjoint, and strictly elliptic system for .

By substituting the mass flux (2) into (5), we can eliminate the velocity variable from the shallow water system entirely,

(5) |

As it turns out, this is the form that the to-be-derived numerical schemes bear close resemblance to.

The inviscid shallow water system (5) is a Hamiltonian system (Salmon1998-eg), and can be put in the canonical form,

(6) |

Here, represents a functional associated with the shallow water system, and the Hamiltonian, which is also a functional and is, for the shallow water system (5), given by

(7) |

The Poisson bracket has three components,

(8) |

and each component is defined as follows,

(9a) | ||||

(9b) | ||||

(9c) | ||||

In the above, , etc., are short-hands for the functional derivatives , etc. The potential vorticity for the SWEs is given by

(10) |

and is the Jacobian operator defined as

(11) |

for any two scalar functions and . The Jacobian operator is skew-symmetric w.r.t. its two argument functions.

As a consequence of the skew-symmetry of the Jacobian operator and the permutations present in its third ( component, the Poisson bracket (8) is skew-symmetric. Therefore, when is replaced by in (6), the right-hand side vanishes, implying that the total energy is conserved. It is also easy to check that quantities in the form of

is a singularity of the Poisson bracket, and therefore quantities in this form, called Casimirs, such as mass, total vorticity, potential enstrophy, etc., are also conserved in the shallow water system.

The discrete numerical schemes are derived in two phases: the Hamiltonian phase and the Poisson phase. In the Hamiltonian phase, one obtains a discretization of the Hamiltonian , defined in (7), by the discrete Hamiltonian, still denoted as ,

(12) |

Here, the subscripted variables, such as , are the discrete version of the corresponding continuous variables, such as . As in CJT1, we adopt the notational convention that subscripted variables (by , , etc.) are discrete, while un-subscripted variables are continuous; un-accented discrete variables are defined at cell centers, while on the top designates edge-defined variables, and designates vertex-defined variables. The symbols and represent the finite difference approximations of the gradient operator and the skew-gradient operator , respectively. The details are omitted, but can be found in Appendix B of CJT1. The factor of in the kinetic energy in (7

) has disappeared in the discrete version here due to the fact that only one component is used in the inner product of the vector fields. One also notes that the skew-symmetry of the Jacobian term is preserved in the approximation here, which leads to a symmetric elliptic system later.

In the Poisson phase, one obtains approximations to the Poisson brackets. Just like its continuous counterpart, the discrete Poisson bracket also has three components,

(13) |

Each component is given by

(14) | ||||

(15) | ||||

(16) |

In the above, where is the potential vorticity at cell edges, which is a remapping of the PV at cell centers, and is defined (see (10)) as

(17) |

The discrete geopotential , vorticity , and divergence are defined as

(18) |

The discrete Jacobian operator is a skew-symmetric approximation to its continuous counterpart,

(19) |

On a global sphere, no boundary conditions are needed, but both the second and last equation of (18) contains redundancy, and the solutions are not unique, for any solution plus some constants will still be a solution to the system. To ensure uniqueness, one can replace one equation from each set with an equation that sets or to a fixed value at a certain grid point. Here, without loss of generality, we pick cell , and set

(20) |

On a bounded domain, the discrete streamfunction satisfies the homogeneous Dirichlet boundary conditions, while the homogeneous Neumann boundary conditions for the velocity potential are only implicitly enforced through the specification of the divergence operator, and the velocity potential is set to zero at cell to ensure unique solvability of the system. Specifically, the boundary conditions for the coupled elliptic system on a bounded domains are

(21) |

With this Poisson bracket (13), a discrete Hamiltonian system can be constructed,

(22) |

An energy-conserving numerical scheme can be obtained by sequentially setting in the above to , , and . Details of the derivation can be found in CJT1, and here we list the energy-conserving numerical scheme for the reference purpose,

(23) |

The terms preceded by in the equation for only appear for boundary cells ().

This scheme just given is guaranteed to conserve the total energy, thanks to the skew-symmetry of the discrete Poisson bracket (13), but it does not conserve the potential enstrophy. The reason is this quantity, given as

is not a singularity of the discrete Poisson bracket. More specifically, it fails to nullify the first component of the Poisson bracket. To remedy this defect, we follow Salmon (Salmon2009-xr), and replace this component by the trilinear Nambu bracket,

(24) |

This trilinear discrete bracket is skew-symmetric with respect to any two of its three arguments, and as a consequence, it vanishes when is set to .

## 3 Numerical evaluations

This section present numerical results from the model. Rather than going through various well-established test cases for the shallow water model, as has usually been done in the literature, here, we focus on a few pre-defined properties of numerical models for large-scale geophysical flows, and examine our model with regard to each one of these properties. Specifically, we will examine our model with regard to six properties: (1) the accuracy of individual operators in the model; (2) the accuracy of the entire model; (3) the conservative properties; (4) the control of the divergence variable; (5) the representation of the energy and enstrophy spectra; (6) the simulation of the barotropic instabilities. We believe that these six properties as a whole are essential for the success of long-term simulations of large-scale geophysical flows.

For this numerical study, a selected set of well-established test cases are used, including the famous shallow water standard test cases (SWSTC) #2 (a steady-state zonal flow) and #5 (zonal flow over a mountain topography) from Williamson1992-cq, and the barotropic instability test case from Galewsky2004-hr. Following Bauer2018-uv, we also use an adapted version of the SWSTC #2 on the northern hemishpere. Finally, we also make use of a new test case of our own (Chen2018-bu) involving a freely evolving gyre in the northern Atlantic.

An external model, the MPAS-sw model, is used for comparison in the examination of certain numerical properties listed above. The MPAS-sw is a C-grid momentum-based shallow water model. This model operates on the same kind of unstructured meshes as our model does, has a similar convergence rate (presumed for our model at this point), but it is designed to conserve the total energy only, not the potential enstrophy. Therefore, it is expected both models will have similar performances in certain aspects, but differ in others. What these differences mean for long-term simulations of geophysical flows is of course an interesting and important quesiton.

### 3.1 Accuracy of individual differential operators

Discrete | Analytic | |

, | ||

, | ||

, | ||

List in Table 1 are all the discrete differential operators that appear in either the energy-conserving scheme (23) or the energy and enstrophy-conserving scheme (26). The analytical differential operators, of which these discrete differential operators are, formally, approximations of, are also listed in the rightmost column of the table. This table contains four groups of differential operators, as divided by the solid lines. Discrete differential operators on the same row are algebraically, as well as numerically, equivalent, while numerical evidence shows that the operators in the same group are very close in accuracy (see below). There are a total of eight discrete differential operators approximating the analytical differential operator . These eight operators are divided into two groups. The first group of three operators (the third group from the table) all have the vertex-to-cell mapping as the last operation, while the other group (the fourth group from the table) has the same mapping but as the very first operation.

All the operators given in Table 1, except the discrete Laplace operator , are combinations of three to four individual discrete operators, e.g. the discrete gradient, the discrete divergence, the cell-to-vertex mapping, and the vertex-to-cell mapping. It is reasonable to expect that these composite operators will remain second-order accurate on a planar centroidal Voronoi mesh with perfect hexagons. However, on a regular centroidal Voronoi mesh which consists of mostly hexagons and a few pentagons, the accuracy will likely be degraded. Rigorous numerical analysis of these operators on unstructured meshes will be carried out elsewhere. Here, we examine the accuracy of these operators numerically.

Due to the prominent roles that the remapping (cell-to-vertex, vertex-to-cell, and edge-to-cell) operators play throughout the schemes, we examine the accuracies of these operators first. The cell-to-edge mapping is not included in this study, because the accuracy of this mapping is well understood to be second-order accurate, due to the fact that every two neighboring cell centers are equi-distant to the common Voronoi edge that separates them. We measure the accuracy of the remapping operatros by comparing the remapped values from one location to another with the actual values of a pre-chosen scalar analytic function. Specifically, we use a scalar analytic function with mild variations in both latitude on longitude,

where and stand for the latitude and longitude, respectively, and the earth radius. Then accuracy of the cell-to-vertex remapping operator, for example, is then computed as follows,

In the above, , and stand for the restriction operators at the cell centers and cell vertices, respectively, and the cell-to-vertex remapping, as in the specification of the schemes.

For this test, we use a series of quasi-uniform spherical centroidal Voronoi tessellations (SCVT) with resolutions ranging from km up to km. The norm and the norm of the errors from all three remapping operators are calculated, and plotted in Figure 1. The norms (left panel) of the errors in the edge-to-cell and vertex-to-cell remapping operators all converge towards zero consistently at the second order, while those of the cell-to-vertex operator converge only at the first order. This degradation in the accuracy of the cell-to-vertex remapping operator is likely due to the fact that it is a upscaling remapping from cells () to vertices (). The errors in all three remapping operators start out similarly as the errors, up to the resolution of km. But, onto higher resolutions, the curves for the vertex-to-cell and edge-to-cell operators start to tilt up, and approach the first-order convergence rate after the km resolution. The curve for the cell-to-vertex operator, which is already at the a lower convergence rate, exhibits a slight degradation at the highest resolutions. The degradation in the accuracies of all three operators are likely indicators of irregularities in the mesh at the higher resolutions.

We study the accuracies of these discrete differential operators of Table 1 by measuring their normalized truncation errors. To do so, we take

where, again, and stand for the latitude and longitude, respectively, and the earth radius. The magnitudes of these variables are chosen so that the magnitude of matches the planetary potential vorticity, while that of matches streamfunction. The discrete variables and are defined by applying the restriction operator to the analytic variable and , respectively. The normalized truncation error of , for example, is defined as

The accuracies of all the discrete operators in Table 1 are studied, using a series of mesh resolutions ranging from up to . Numerical results show that the normalized truncation errors for operators in the same group are very close, agreeing up to at least the third digits. For this reason, we only plot the results from the first representative member in each group in Figure 2. On the left are the normalized -errors, and it is seen that the errors from all four representative operators converge at a rate between the first and the second orders, with the representative operator from the third group being the most accurate (at the 2nd order) and the representative operator from the fourth the least accurate (at the 1st order). The performances of the operators from the first and second groups are in between these two extremes. The difference in the performances of the third and fourth groups are caused by the applications of the remapping operators. The operators in the fourth group all start with the cell-to-vertex mapping, while the operators in the third group all finish with the vertex-to-cell mapping. Both mappings are area weighted, but, as shown in Figure 1, the vertex-to-cell mapping is more regular, as it is a downscaling mapping from vertices () to cells (), while cell-to-vertex mapping is more singular, as it is a upscaling mapping from cells () to verteices (). The convergences for the operators from all groups are consistent across the whole range of mesh resolutions, except the operator at the highest resolution . We suspect that this flat-out is caused by an anomaly in the mesh.

On the right of Figure 2 are the normalized errors. The errors for most of the operators are not converging, which is expected. It is therefore surprising to see that the errors for the representative operator from the third group actually converges at the first order. This operator also has the best and most consistent converging rate.

The above results concerning the individual operators of Table 1, except the third group, largely agree with those of (Eldred2015-yg, Section 6.2) on the SCVT grids in that the operators are between first and second order accurate in the norm, and non-consistent in the norm. In addition to the SCVT, other types of optimized grids are also considered in the work just referenced, and one particular type, the so-called tweaked iscosahedra, is shown to lead to consistent operators in both the and norms. Operators in the third group on Table 1 are not needed, and hence not studied in Eldred2015-yg.

### 3.2 Accuracy of the entire model

#### 3.2.1 A stationary zonal flow on a global sphere (SWSTC#2) and on a hemisphere

The SWSTC #2 specifies an initial zonal flow that is in perfect geostrophic balance. Thus, this test case not only provides a rare means of testing the accuracy of numerical scheme with an analytic solution to the full SWEs, it is also an excellent test on the capability of a numerical scheme to maintain the geostrophic balance. In order to test the impact of the presence of a boundary on the accuracy of the numerical scheme as well as its capability to maintain geostrophic balance, we follow Bauer2018-uv and run the same test, with the same initial conditions, but on the northern hemisphere, with the equator serving as the southern boundary. Both the EC scheme and the EEC scheme are used for this test case, and the results are very similar. Hence we only present and discuss the results from the EEC scheme.

The exact configurations of the SWSTC #2 can be found in the original Williamson et al paper (Williamson1992-cq) and many other numerical papers that follow. We run the test for 5 days, as recommended in the original paper, on a suite of quasi-uniform meshes on the global sphere, with resolutions ranging from km up to km. Then we compare the solutions on day 5 with the intial state, and compute the normalized and errors. The same test is then carried out again, on a suite of quasi-uniform mesh on the northern hemisphere.

The results on the global sphere are plotted against the grid resolutions in Figure 3. together with the first- and second-order reference convergence curves. It is seen that the errors for both the thickness and vorticity variables converge at the the second order, with the vorticity errors larger but converging more consistently, all the way up to the 30km resolution. Since the flow starts out as non-divergent, and hence the area-weighted and -normalized absolute errors for the divergence variables are used. The errors for the divergence variable stay below for all resolutions. The relative errors of neither the thickness or the vorticity show any sustained converging trend, and the error of the divergence variables shows a steady growing trend. At 30km, the maximum magnitude of the divergence variables reaches about .

For comparision, the thickness variable under the extended Z-grid scheme of Eldred Eldred2015-yg on this particular test case converges at approximatley the first order in both the and norms, on the so-called tweaked icosahedral grids. Under our scheme, the same variable converges slightly faster in the norm on general SCVT grids, but does not converge in the norm, consistent with the convergence behaviors of the individual operators. The convergence behavior of the vorticity variable is not discussed in the work just referenced.

The result from the hemisphere case are plotted against the grid resolutions in Figure 4, together with the first- and second-order convergence reference curves. In this case, the relative errors in the thickness converge at an approximately first-order rate. The relative errors for the vorticity variable are larger, but converge faster, at an approximately second-order rate. The errors in neither variable converge consistently. As we already pointed out, the flow starts out non-divergent, and the magnitude of the divergence variable can only go up. The figure shows that magnitude in the divergence variable also go up as the grid refines, with the area-weighted and -normalized errors reaching at approximately , and the reaching at the 30km resolution.

Compared with the SWSTC #2 on the global sphere, the errors in the current case are larger by a factor of 2 to 5, across all variables and for both norms, and the convergences are also slower. The downgraded accuracies can be attributed to the presence of boundaries and the downgraded uniformity in the grid resolution. The centroidal Voronoi tessellations on the global sphere are generated using the isohedra as starting points, and the resulting meshes are highly uniform, with a ratio of between the highest and lowest resolutions, as represented by the cell-to-cell distances. The centroidal Voronoi tessellations on bounded domains start with a set of prescribled points on the boundary and a set of largely random points in the interior. With iterations, the meshes can achieve a ratio of between the highest and lowest resolutions, which qualifies as quasi-uniform meshes, but certainly worse than the meshes on a global sphere.

#### 3.2.2 Zonal flow over a mountain topography

This classical test case starts with an initially zonal flow similar to that of the SWSTC #2. The flow impinges on a mountain topography centered at latitude and longitude , and gradually evolves into a turbulent flow around day 25. The exact configurations of this test case, again, can be found in the reference Williamson1992-cq, and will not be repeated here. This test case was originally used by Takacs (Takacs1988-tx) to study the conservation of integral invariants by a posteriori methods. Subsequently, this test case is often used by authors to qualitatively study the capability of numerical schemes to simulate the nonlinear dynamics in geophysical flows (Ringler2010-sm; Weller2012-md; Chen2013-fa; Li2015-nb).

For this study, we again use a suite of quasi-uniform meshes over the global sphere, with resolutions ranging from km up to km. The relative errors in the thickness field on day 15, under both the and the norms, are computed by comparing our solutions to the high resolution spectral solution H213 from (Jakob-Chien1995-ci). The relative errors are then plotted against the grid resolutions in Figure 5, along with the reference first- and second-order convergence curves. It is seen that the errros in both norms converge consistently at a rate between the first and second orders. It should be pointed out that the reference solution is computed with diffusions, while our solutions are computed without any diffusions.

### 3.3 Conservative properties

It has been proven in CJT1 that the numerical scheme conserves the first order moments, such as the mass and the absolute vorticity, up to the machine errors, and the total energy and enstrophy, which are second-order moments, up to the time truncation errors. These results are confirmed in the numerical study. Here we present the results from two tests. The first test is on a bounded domain, and the results from this test illustrate the impact of the time truncation errors in the conservation of the total energy and the potential enstrophy. The second test is on a global sphere, and involves comparison to a third-party model that is designed to conserve the total energy but not the potential enstrophy.

In numerical schemes, the boundary is often a source of errors. Inconsistent treatment of the boundary conditions can even result in unstable simulations. In order to numerically verify the conservative properties of our numerical scheme and examine its treatment of the boundary conditions, we design a test case with a freely evolving gyre in a bounded domain, with no external forcing or diffusion. Physically, both the total energy and the potential enstrophy should be conserved, alongside the mass and total vorticity, and this should also be the case in numerical simulations with schemes that are designed to conserve these quantities. For the physical domain, we consider one section of the mid-latitude northern Atlantic ocean. The initial state of the flow is given by

In the above, and are the latitude and longitude of the center point of the Atlantic section, , . The function is used to ensure that the initial stream function is flat (with one constant value) near and along the boundary of the domain. A plot of the domain and of the stream function is shown in Figure 6. This stream function produces a circular clockwise velocity field in the middle of the ocean, with the maximum speed of about m/s near the center, which is considered fast for oceanic flows. No external forcing or diffusion is applied, and the flow is allowed to evolve freely. The circular pattern will breakdown eventually, because of the irregular shape of the domain and because of the non-uniform Coriolis force.

We cover this bounded domain with a non-uniform 3088-cell SCVT mesh with resolutions ranging from about km to km. The EEC scheme is employed, together with the 4th order Runge-Kutta time stepping scheme, and the simulation is run for one calendar year (360 days). The total energy and potential enstrophy are computed for each day, and these statistics are then plotted against time in Figure 7. In order to illustrate the impact of the time truncation errors, we run the simulation twice, with two time step sizes, s (left) and s (right). At the end (day 360) of the first run with the s time step size, the system loses of the total energy and of the potential enstrophy (Figure 7 (a)), with the latter almost in the range of roundoff errors for a simulation of this length. At the end of the second run with the s time step size, the system loses of the total energy, a reduction compared to the first run by a factor of about 8, and it loses about of the potential enstrophy, a reduction by a factor of about 4. These result confirm that the EEC numerical scheme conserve the total energy and potential enstrophy up to the time truncation errors, and with reduced time step sizes, the truncation errors do come down.

Our second test involves an initially zonal flow over a mountain topography over the global sphere(SWSTC #5, Williamson1992-cq). The purposes of this test are two folds. First, it seeks to verify the conservative properties of the EEC scheme in different setting (flow over topography in a global domain), and second, it seeks to contrast the results from our EEC scheme with those of the MPAS-sw model, which operates on the same kind of global SCVT meshes, has similar convergence rates, but is not designed to conserve the potential enstrophy. For this test, both models use a global quasi-uniform SCVT mesh with 40962 cells (res. 120km) and the RK4 time stepping scheme. The simulations are run for days, and the total energy and potential enstrophy are computed for each day, and then plotted against time in Figure 9. In the EEC scheme (left panel) the flow experiences a weak decay in both the total energy and the potential enstrophy, likely due to the numerical diffusion of the Runge-Kutter 4th-order time stepping scheme. By the end of the simulation, the flow loses about of the total energy and about of the potential enstrophy. In the MPAS-sw model, the evolution curve of the total energy is similar to that of our EEC scheme, but the potential enstrophy increases by about 40% at the end of the simulation. This shows that numerical accuracy alone cannot guarantee the conservation of all dynamical quantities. Our EEC scheme and the MPAS-sw model share a similar convergence rate, but have dramatically different behaviors in terms of potential enstrophy.

It is worth noting that, our EC scheme, while not proven to conserve the potential enstrophy, manages to only accumulate about of it (results not plotted here), highlighting the promising prospect of the Hamiltonian approach for numerical scheme design.

### 3.4 Control of the divergence variable

Large-scale geophysical flows constantly evolve around an approximate geostrophical balance, i.e. a balance between the Coriolis force and the pressure gradient. Therefore, it is crucial for numerical scheme to maintain this geostrophic balance in order to be able to accurately simulate the dynamics. It has long been recognized that an accurate representation of the dispersive wave relations is crucial to the success in this task, because a proper representation of the dispersive wave relations allow energies to propagate properly across scales rather unphysically accumulate at one scale, and eventually destroy geostrophic balance (Arakawa1977-og; Randall1994-vu). Besides the excellent conservative properties, our EEC scheme also possess the optimal dispersive wave relations. How does this last property contributes to its capability to simulate and maintain geostrophically balanced flows? We need an observable indicator. When flows evolve around the geostrophic balance, the flow is close to being, but not exactly, non-divergent. Therefore, the divergence variable can be an indicator to the capability of a numerical scheme to maintain the geostrophic balance. Ideally, the divergence variable should remain small throughout the simulation. In the left panel of Figure 9, we plot the evolution of the area-normalized norm and the norm of divergence variable, again on the km quasi-uniform global mesh, up to day 250. It is seen that the norm of the divergence variable remains largely flat, exhibiting no noticeable growth. The norm is larger, as expected, but remains largely below . Again, no noticeable and sustained growth is observed. For comparison, we also plot the results by the MPAS-sw model on the same scale, on the right panel of the figure. MPAS-sw is a C-grid scheme, with decent but less accurate (vs Z-grid) representations of the dispersive wave relations. In the plot, the norm of the divergence exhibits a slow but steady growth. The norm of the divergence frequently exceeds the upper limit of the panel, and, despite being highly oscillatory, also exhibits a noticeable growth over time.

### 3.5 Representation of the energy and enstrophy spectra

Our EEC scheme excels in combining the desirable conservative and dispersive properties into one scheme. Dynamically, what real advantages does this scheme have compared to other numerical schemes that are currently in operation? Answers to this big question can only be sought by applying the scheme to real-world applications. These important questions will be pursued in the future. Here we examine one crucial dynamical aspect of the numerical scheme, namely its representation of the spectra of the kinetic energy (KE) and potential enstrophy (PEns), and the results are compared with those from the MPAS-sw model, to see what advantages, if any, the EEC model have. According to the phenomenological theories of two-dimensional turbulence (Batchelor1969-fo; Kraichnan1967-un; Lilly1969-wl), the KE spectrum of a turbulent two-dimensional has a slope of in the inertial range, while the spectrum of the PEns has a slope of . These theories and predictions have been largely confirmed by numerical studies; see e.g. Chen2011-jh; Maltrud1991-wa. While a correct representation of the KE and PEns spectra does not guarantee a faithful simulation of the nonlinear dynamics of turbulent flows, it is a prerequisite.

For this study, we again use the SWSTC #5 with a free zonal flow over a mountain topography. We run both EEC and MPAS-sw over the same global quasi-uniform km mesh for 250 days, and we examine and plot the vorticity field and the KE and PEns spectra on day 50 and day 250. In Figure 10 are snapshots of the vorticity field on day 50 from the EEC and the MPAS-sw models. Placed below the snapshots are the plots of the KE and PEns spectra at the same instantaneous moment. The vorticity fields from these two models are qualitatively similar, while the vorticity field of the MPAS-sw appears smoother, due to the fact that, in MPAS-sw, the vorticity is natively defined at the vertices, and is mapped to the cell centers for the visualization purpose. For the EEC model, the KE spectrum demonstrates a slope of from wave number 10 to 20, while the PEns spectrum demonstrates a slope of from wave number 10 to 30, approximately verifying the phenomenological theories. The KE and PEns from the MPAS-sw model exhibit similar trends. For both models, the KE and PEns spectra decay quickly after wave number 30, apparently due to the fact that fluid motions still largely concentrate at large scales, and “cascade” is not complete yet. This completely changes at day 250. For the EEC model, the KE spectrum now display a slope of from wave number 10 up to wave number 40, indicating that a significant portion of the fluid motions has shifted from the large scales to small scales. Due to this shifting, and due to the lack of diffusion, the PEns spectrum is now distorted with the PEns piling up at high wave numbers, and with no discernible part of the curve conforming to the slope. The results from the MPAS-sw are worse. With this model, the KE spectrum is completely distorted with no discernible part conforming to the slope, and so is the PEns spectrum. For the PEns spectrum, not only is the slope is gone, the peak density at high wave number is much larger than the peak density of the PEns from the EEC model, which is most likely due to the spurious source of PEns that we saw in Figure 9, right panel. Snapshots of the vorticity fields corroberate what the KE and PEns spectra tell us. Both snapshots are noisy, but the one from the EEC model clearly still possess some large-scale structures, while the one from the MPAS-sw is completely noisy, with no discernible structures.

### 3.6 Simulation of barotropic instabilities

Instabilities, both barotropic and baroclinic, in the ocean and atmosphere are responsible for the transient and diverse weather phenomena. Therefore, in weather forecast, where accuracy is the primary concern, it is vitally important to accurately simulate the instabilities. But instabilities often happen erratically, if not randomly, and accurately simulating the instabilities in the ocean and atmosphere has been a challenge. The one-layer shallow water model, by nature, is not capable of simulating the baroclinic instabilities. Here, we only exam the capability of our model to simulate the barotropic instability. We will use the familiar Galewsky et al barotropic instability test case (Galewsky2004-hr). This test case starts with a perfectly balanced fast zonal jet at north. Due to the barotropic instability, this jet will quickly start to wobble, and eventually give rise to a stream of eddies alongside the jet around day 5, in line with the time scale observed in our weather system. The exact setup of this experiment can be found in the original reference.

As just mentioned, initially, the thickness and zonal velocity field are in a perfect geostrophic balance. A model with unlimited precision can maintain this balance indefinitely. As pointed in the original article Galewsky2004-hr, maintaining this balance is trivial for spectral models using spherical harmonics, because, with all the modes initially horizontal, it will be very difficult for non-horizontal modes to grow. This task is, on the other hand, challenging for low-order models based on grids that are not aligned with the latitudes, especially unstructured meshes such as ours. In these models, the discretization errors can quickly get into the system, and distort the geostrophic balance. In our first test, we run the test case using the unperturbed thickness field up to day 5, using a sequence of quasi-uniform SCVT meshes, with mesh resolutions ranging from 480km up to 30km. The goals are (i) to see the effect of discretization errors at various resolutions, as measured by the departure from the exact geostrophic balance, and (ii) to see what resolution is required to maintain a decent geostrophic balance up to day 5. In figure 12, we plot the vorticity field on day 5 from the simulations. In the top two panels for 480km and 240km resolutions, the vorticity exhibit a significant depart from the purely zonal flow, with a Rossby wave as wide as in latitude. At 120km (third panel), which is traditionally placed in the mesoscale range, the Rossby wave width is significantly reduced. At 60km (fourth panel), the maintenance of the geostrophic balance and zonal flow structure is decent.

In order to qualitatively compare the dynamics simulated by our model with previous results, we pick the quasi-uniform km mesh, as it is shown in the above that the effect of the truncation errors on this mesh on the overall dynamics is minimal. In Figure 13, we plot the vorticity field on day 6, for both the inviscid and viscous cases, with the perturbed thickness field. In the inviscid case, the instability is obviously more pronounced. In both cases, the results match those of the previous studies (Galewsky2004-hr; Chen2008-gz; Weller2012-jw).

Left unexplored in the work is the question of whether, in the inviscid case, the solution actually converges to a regular solution free of singularities (shocks). Common wisdom say that shocks are inevitable in a compressible system, in the absence of diffusions. The shallow water equations is a slightly compressible system, and the Coriolis force is known as a regularizing effect. It is not clear whether this regularizing effect can help to suppress the emergence of the shocks. No theoretical results exsit concerning one way or the other. In Galewsky2004-hr, the convergence of the solution in the viscous case is established numerically. The issue in the inviscid case was not explored, due to the prohibitive computational costs. Such costs are also beyond our reach at the moment.

## 4 Discussions

This overarching goal of this work is to evaluate the newly developed conservative numerical scheme (EEC), so as to preliminarily assess its potentials in real-world applications. To this end, this work focuses on, and is organized around, a set of pre-defined properties such as accuracy and conservation of key quantities. The evaluations make use of a suite of old and new test cases.

With regard to accuracy, the EEC scheme has a convergence rate between the first and second orders, as evidenced in the SWSTC #2 and #5. On bounded domains, the convergence rate is slightly lower, but remains within that range. On the SWSTC #2, our scheme exhibits a slightly faster convergence rate, second vs first, than the extended Z-grid scheme of Eldred (Eldred2015-yg). This could be attributed to the fact that we use the more accurate area-weighted mapping throughout the scheme, while the extended Z-grid scheme uses less accurate remapping operators for the sake of conservations. Both our work and that of Eldred demonstrate that the accuracy of the numerical scheme is sensitive to the mesh quality.

The EEC scheme is proven to conserve the mass and vorticity up to the roundoff errors, and to conserve the total energy and potential enstrophy up to the time truncation errors. These results are confirmed by the numerical experiments conducted in this study, including one on a bounded domain. A more intriguing question is what advantage(s) a comprehensively conservtaive numerical scheme has in a long-term simulations. Arakawa and Lamb (Arakawa1981-dy) argues, via a test case with a channel flow over a ridge, argues that a doubly conservative scheme can prevent unphysical energy cascade towards the finest resolvable scales. Our study, via a different test case and through a different perspective, confirms this point. We showed that, in the SWSTC #5, the EEC scheme, compared withan energy-conserving only scheme (MPAS-sw), can deliver a more structured vorticity field at the end of a long simulation (250 days), with more physically realistic energy and enstrophy spectra.

The EEC scheme studied here has been shown to possess the optimal dispersive wave relations among the second-order accurate numerical schemes. It is well known that properly represented dispersive wave relations with positive group speeds are essential to the geostrophic adjustment process, through which large-scale geophysical flows evolve among geostrophically balanced structures. Geostrophic balance is hard to verify directly without the expensive reconstruction of the velocity and pressure fields, but geostrophically balanced flows are close to being non-divergent. Therefore, the divergence variable is a natural choice as an indicator of the maintenance of geostrophically balanced structures in the flows. A flow that starts out geostrophically balanced with a small divergence variable should remain so throughout the simulation period. This is what is observed in the SWSTC #5 with the EEC scheme: the divergence variable remains small and exhibits no discernible growth over a period of 250 days. In contrast, the divergence variable from the MPAS-sw model, based on a C-grid scheme with a decent representation of the dispersive wave relations, shows a slight but unmistakenable growth over the same period.

Effort has also been made to assess the capability of the scheme to simulate flow dynamics close to those in the real world. Snapshots of the vorticity field from the Galewsky et al test case demonstrate that while the scheme has a low order of accuracy, it is capable of simulating the barotrophic instabilities at mid- to high-resolution, even over unstructured meshes. Unfortunately, the potentials of this scheme in real-world applications cannot be fully assessed in its current implementation, due to the limitations posed by the shallow water model, chief among them the lack of vertical variations and convections. Hence, more work remains to be done in terms of development and testing. On the development side, the scheme needs to be extended and implemented for more complex models that incorporate vertical variations and convections, such as the primitive equations. On the testing side, real-world applications, or well-designed test cases that are closer to the real-world applications, need to be utitlized on the numerical scheme. These important undertakings will be reported elsewhere.

## Acknowledgements

If you’d like to thank anyone, place your comments here and remove the percent signs.

Comments

There are no comments yet.