Numerical instabilities of spherical shallow water models 3 considering small equivalent depths 4

Shallow water models are often adopted as an intermediate step in the development of atmosphere and ocean models, though they are usually tested only with fluid depths relevant to barotropic fluids. Here we investigate numerical instabilities emerging in shallow water models considering small fluid depths, which are relevant for baroclinic fluids. Different numerical instabilities of similar nature are investigated. The first one is due to the adoption of the vector invariant form of the momentum equations, related to what is known as the Hollingsworth instability. We provide examples of this instability with finite volume and finite element schemes used in modern quasi-uniform spherical grid based models. The second is related to an energy conserving form of discretization of the Coriolis term in finite difference schemes on latitude-longitude global models. Simple test cases with shallow fluid depths are proposed as a means of capturing and predicting stability issues that can appear in three-dimensional models using only twodimensional shallow-water codes. 11


Introduction
In the development of an ocean or atmosphere model it is important to establish that the numerical methods used are stable. In practice, a mathematical stability analysis might not be tractable, for example if the geometry, grid structure, or background state are not simple enough; in that case stability can only be determined by empirical tests. It may also be the case that a mathematical analysis is not done because a particular type of instability was not anticipated. It is highly desirable to be able to identify any instability of some candidate numerical method before investing the effort to develop a three-dimensional model. Shallow-water models are often adopted as an intermediate step in the development of atmosphere and ocean models and can be useful in testing the stability of numerical methods. However, standard numerical test cases (e.g. Williamson et al., 1992;Galewsky et al., 2004;Shamir and Paldor, 2016) are usually run with fluid depths of thousands of metres, relevant to barotropic fluids. In this paper we point out the usefulness of spherical shallow-water models run with very shallow fluid depths, mimicking the internal modes of three-dimensional models, for testing the stability of numerical methods. Hollingsworth et al. (1983) discovered that an implementation of a well-known energy-and enstrophy-conserving scheme, originally due to Sadourny, when used in a hydrostatic primitive equation model (Burridge and Haseler, 1977), was prone to near-grid-scale instabilities, with severe consequences for high-resolution forecasts. This kind of instability is intrinsically related to the use of the vector-invariant form of the momentum equations (e.g. Vallis, 2006), which expresses the advection of momentum as where v = (u, v) defines the horizontal velocities, K = v 2 /2 is the kinetic energy, ζ = k · ∇ × v is the relative vorticity and k is a unit vector pointing normal to the horizontal surface. Expanding the differential expressions, and looking at the equations for each velocity component separately, one notices that on the left-hand side of the above relations there are no derivatives of v in the u equation, and no derivatives of u in the v equation. In contrast, we see that these derivatives exist in each of the two terms on the righthand sides, respectively, but they cancel out. The cancellation is a general feature of this formulation of the equation, and exists independently of the particular choice of coordinate system. In numerical schemes, this cancellation is not always exact, and this may lead to stability issues. This is what happens with the original energy-and enstrophy-conserving scheme of Sadourny (Burridge and Haseler, 1977), hereafter referred to as the 'EEN' scheme. Hollingsworth et al. (1983) proposed a modification of the scheme to ensure cancellation at a discrete level of the linearized equations, avoiding the instability while preserving its conservation properties. Similarly, the energy-and enstrophyconserving scheme of Arakawa and Lamb (1981), hereafter referred to as the 'AL' scheme, was also shown to be prone to such instabilities, and again a modification was proposed to avoid it.
There has been a renewed interest in these instabilities, as the vector-invariant form of the momentum equations have recently been adopted in many novel atmospheric and ocean models on quasi-uniform spherical grids (e.g. Tomita et al., 2008;Skamarock et al., 2012;Gassmann, 2013;Wan et al., 2013;Ringler et al., 2013). As a consequence, novel high-resolution ocean and atmosphere models are presenting instabilities of a similar nature (e.g. Skamarock et al., 2012;Gassmann, 2013). A problematic point is that such instabilities are only being detected once the full 3D model is implemented, as no trace of instability is usually detected in shallow-water prototypes. Bell et al. (2017), hereafter referred to as BPT, examined the instability for the two well-known energy-and enstrophyconserving schemes, the EEN and AL schemes, for the vector-invariant hydrostatic Boussinesq equations. As found by Hollingsworth et al. (1983), the schemes were shown to be linearly unstable for height coordinate models. Also, BPT showed that it is possible to detect such instabilities on shallow-water versions of the schemes, as long as the model adopts a small equivalent depth, that is, the shallow-water layer is very thin. This allows easy testing of novel and existing schemes for unstable linear modes of similar nature to the original instability detected by Hollingsworth et al. (1983). All the analysis of BPT was performed for planar quadrilateral grids, but it suggests a way of testing for the instability on more general discrete domains, such as the quasi-uniform spherical ones.
Shallow-water models considering small equivalent depths are directly related to reduced gravity layer models (e.g. Vallis, 2006). Reducing the inertia-gravity wave speed increases the importance of the nonlinear terms, such as advection terms in the shallowwater equations, and so exacerbates the effect of any numerical errors in those terms. Gassmann (2011) shows an experiment that mimics the case of relatively small equivalent depths occurring in atmospheric models using a planar triangular-hexagonal grid shallow-water model. The results allow interpretation, within a simpler 2D shallow-water framework, of how a checkerboard divergence mode is expected to interfere in a 3D model. Similarly, Peixoto (2016) uses a shallow-water experiment with reduced depth to illustrate and foresee inaccuracies that might appear due to the finite-volume discretizations of the nonlinear terms on quasi-uniform spherical grids.
The main goal of this article is to discuss possible numerical instabilities that may arise in shallow-water spherical models when small equivalent depths are adopted. Two different kinds of instabilities that appear in linear analysis for the equations linearized about a non-resting basic state with small depth are analysed. One of the instabilities is the one arising from the use of the vector-invariant form of the momentum equations, which is of similar nature to that of Hollingsworth et al. (1983) and is what motivated this study. The other instability occurs in a latitude-longitude C-staggered semi-Lagrangian semi-implicit finite-differences model, as used, for example, in the Unified Model of the UK MetOffice, also known as ENDGame (Wood et al., 2014). The scheme used in ENDGame is based on the discretization proposed in Zerroukat et al. (2009), where the Coriolis term discretization, brought forward from Thuburn and Staniforth (2004), was inspired by the work of Arakawa and Lamb (1981). Therefore, the two cases analysed are more intimately related than it would seem at first, as will be discussed in what follows.
A key concept used in this article is that the stability analysis of three-dimensional models may be done by separating the linear modes into horizontal and vertical parts, and these are connected by equivalent depths. A discussion of the equivalent depths emerging in typical three-dimensional sets of equations is presented in section 2, followed by a discussion of how different vertical coordinates may be interpreted either using depthweighted or non-depth-weighted vorticity terms in shallow-water systems. Based on these discussions, we can thereafter limit our attention only to shallow-water equations, but with conclusions that can be interpreted for three-dimensional models.
In section 3 we discuss the influence of the Coriolis force in the Hollingsworth instability, showing that the instability exists even in models without Coriolis terms. Section 4 presents the common framework of tests and methods that will be used to investigate several different spherical shallow-water models. Section 5 examines the Hollingsworth instability in quasi-uniform spherical grids, such as cubed and icosahedral spheres, for finite-volume and finite-element schemes. Section 6 analyses the instabilities detected for ENDGame, starting from an analytical examination of the linear modes on a plane, followed by numerical experiments on the sphere.

Equivalent depths in 3D models
A relatively good approximation to ocean dynamics is to consider the hydrostatic, incompressible, adiabatic, Boussinesq equations. The linearized Boussinesq equations enjoy separable solutions (Gill, 1982;Vallis, 2006), and the important terms connecting the vertical and horizontal modes are the equivalent depths, which are eigenvalues of the vertical mode problem. These are named equivalent depths because they give rise to shallowwater systems with such mean fluid depth. Bell et al. (2017) analyse the Boussinesq equations considering a constant Coriolis parameter plane and two vertical coordinate systems: height and isopycnal. The numerical vertical modes are investigated for a Lorenz and a Charney-Phillips staggering, respectively for height and isopycnal coordinates, and estimates of equivalent depths for real application parameters are discussed. On a Charney-Phillips grid, considering 100 vertical layers, typical ocean parameters can result in equivalent depths of less than a metre. The Lorenz grid has equivalent depths inversely proportional to the square of the number of modes, so when many vertical levels (e.g. 100) are used, this may result in equivalent depths smaller than millimetres.
Many weather and climate models adopt the primitive equations for modelling the atmospheric dynamics (Lauritzen et al., 2011;Holton and Hakim, 2012). Similar vertical versus horizontal separation is also possible for the linearized primitive equations, and again the key connecting parameters are the equivalent depths (Tribbia and Temam, 2011). Terasaki and Tanaka (2007) investigated the equivalent depths occurring in the primitive equations. For a fully spectral analysis and standard atmospheric parameters, considering 22 vertical spectral modes, the minimum equivalent depth calculated was about 8 m. Also, Kasahara and Puri (1981) perform a full analysis of the 3D modes and calculate the equivalent depths emerging for this kind of equation set for a sigma coordinate model. An example for a model discretized vertically with finite differences with nine sigma levels shows that the smallest equivalent depth is 3 m. The smaller equivalent depth estimated in Kasahara and Puri (1981), even with only nine levels, is due to use of a vertical finite-difference scheme, whereas Terasaki and Tanaka (2007) consider a spectral vertical analysis.
High-resolution global atmospheric models often adopt the full, non-hydrostatic, compressible Euler equations. Under the shallow atmosphere approximation, the linear compressible Euler equations have normal mode solutions that separate into a product of a vertical structure function and a horizontal structure function; however, in contrast to the hydrostatic case, the vertical structure equation involves the mode frequency (Daley, 1988).
Consequently, the horizontal structure equation no longer has a single equivalent depth independent of the mode frequency. Nevertheless, for gravity modes with large vertical wavenumber, which are the most problematic modes for the instabilities of interest here, non-hydrostatic effects are rather weak. Therefore, analysis based on the hydrostatic assumption, particularly the smallest equivalent depths, should give a useful indication of model behaviour in the non-hydrostatic case too.
To summarize, in realistic atmosphere and ocean models the equivalent depth of higher internal modes is expected to be small (only a few metres or smaller), and can be much smaller depending on the vertical coordinate and grid adopted.

Form of the vorticity term, and relation to vertical coordinate systems
For the vector-invariant form of the shallow-water equations, there are two distinct ways of writing the vorticity term: where f is the Coriolis parameter, η is the fluid depth, and q = (f + ζ )/η is the potential vorticity. We will refer to these as the non-depth-weighted form and depth-weighted form, respectively. For the continuous equations the two forms are equivalent, but this is no longer the case after discretization. The depth-weighted form is attractive because it facilitates the design of numerical schemes that conserve potential-vorticity-related quantities such as potential enstrophy (Arakawa and Lamb, 1981;Ringler et al., 2010). Importantly, the presence or absence of depth weighting can affect the stability of the scheme (BPT, also see below). A similar distinction arises when the advective form of the momentum equation is used, except that now only the Coriolis term is affected. The non-depth-weighted and depth-weighted forms are respectively. When analysing the stability of a three-dimensional numerical method by separation into a vertical-structure problem and a shallow-water problem, it is important to determine which form of the vorticity term is appropriate. The form of the vorticity term in a three-dimensional model is closely related to the type of vertical coordinate it uses.
Layer-based vertical coordinates offer the flexibility to choose either a depth-weighted or a non-depth-weighted discretization of the vorticity term, and the depth-weighted option is a common choice. By 'layer-based vertical coordinate' we mean one that carries some quasi-Lagrangian information about the thickness of material layers. Common examples of layerbased coordinates include isentropic, isopycnal, and Lagrangian coordinates. Throughout much of the atmosphere, diabatic heating is weak and potential temperature θ is approximately materially conserved. Thus isentropic surfaces (surfaces of constant θ ) are approximately Lagrangian surfaces, and a coordinate system with θ as the vertical coordinate provides a quasi-Lagrangian coordinate system that has some attractive features for atmospheric modelling (e.g. Hsu and Arakawa, 1990). In an analogous way, potential density is approximately materially conserved throughout much of the ocean, providing the basis for an isopycnal vertical coordinate (e.g. Bleck and Boudra, 1981). For atmospheric modelling a Lagrangian vertical coordinate, in which a set of material surfaces are used to define the vertical coordinate, has also been used (e.g. Lin, 2004). In practice all such models require some measures to ensure that model layers do not fold over or become too thin. Nevertheless, they all carry information about the thickness of quasi-Lagrangian layers and so can use the depth-weighted vorticity or Coriolis terms.
In level-based coordinates, on the other hand, the non-depthweighted form of the vorticity term is generally used. Level-based models adopt a monotonic variable, such as geometric height, to define fixed vertical levels. For the present discussion, pressurebased coordinates and mass-based coordinates, including their terrain-following variants, should also be thought of as levelbased. Although the height of pressure levels or mass levels can vary in time, these variations are relatively small. Levelbased coordinate systems do not carry direct information on the thickness of material layers. In principle some estimate of material layer depth such as (∂θ/∂z) −1 could be calculated in order to use a depth-weighted vorticity term; however, the numerical errors in this estimate would affect small vertical scales and would need to be accounted for in any stability analysis.
In some height-, pressure-, or mass-based atmospheric models the vorticity or Coriolis terms are weighted by density or by a pseudo-density proportional to density times model layer thickness (e.g. Skamarock et al., 2012;Wood et al., 2014;Dubos et al., 2015). However, local variations in density are relatively small, so this density weighting does not have the same effect as depth weighting. Similarly, the thickness of model layers does not correspond to the thickness of material layers (Arakawa, 2000), so pseudo-density weighting does not have the same effect as depth weighting. Therefore, models based on such coordinate systems must be interpreted as using non-depth-weighted vorticity or Coriolis terms and analysed accordingly.

The role of the Coriolis force in the Hollingsworth instability
We discussed in the introduction how the Hollingsworth instability is connected to the vector-invariant form of the momentum equations through the lack of a certain discrete cancellation. The term responsible for the lack of cancellation solely involves momentum advection, so it might be expected that the Coriolis term should play no role in the existence of such instability. Lazić et al. (1986) observed that real data runs of a 3D finitedifference ECMWF * model, which used the original energy-and enstrophy-conserving EEN scheme, collapsed after a couple of days, showing accumulation of energy in short waves with the instability naturally linked to the Hollingsworth problem. They performed a linear analysis of a planar shallow-water version of the model discretized with the EEN scheme considering a constant background velocity field and a constant Coriolis parameter. Confirming the results from Hollingsworth et al. (1983), the shallow-water model was shown to be linearly unstable. Interestingly though, the system was shown to be stable if the Coriolis parameter was set to zero.
Following the analysis of BPT it is possible to show analytically that indeed for a constant basic background velocity and null Coriolis parameter the system is neutrally stable. Since this was not explicitly proven in either BPT or Lazić et al. (1986), we describe the proof in Appendix A. Also, BPT show how the non-dimensional growth rate (ω i ) of the instability is related to the Coriolis parameter, and that it is stronger for larger grid Rossby numbers, and do not show what happens in the case of absence of rotation (f 0 = 0). BPT normalize the growth rate using the Coriolis parameter, so that the dimensional growth rate is actually f 0 ω i . An increase in the grid Rossby number, considering a fixed grid and fixed velocity, is related to the decrease of the Coriolis parameter (f 0 ). BPT shows that with increasing grid Rossby number the non-dimensional growth rate (ω i ) also increases, but this increase in ω i is small compared to the decrease in f 0 . So, in fact, the dimensional growth rate ω i f 0 reduces with a reduction of f 0 , even with an increase of the grid Rossby number. In the limit, the dimensional growth rate is zero once f 0 is also zero, in agreement with the analysis of Appendix A.
This all goes against our intuition about the known source of the instability, which is the lack of discrete cancellation of terms related to momentum advection solely. This is elucidated in what follows.
The momentum advection plus Coriolis term may be written in vector-invariant form as where we see that the Coriolis term is added to the vorticity term. The basic velocity states used in both analysis, BPT and Lazić et al. (1986), did not have background vorticity. Following the notation from BPT for the planar analysis, let us impose a non-rotating (f 0 = 0) system, but add a background basic flow with vorticity as where u 1 , v 1 are constant velocities and 2ζ 0 will be the basic flow vorticity. Assume a constant basic depth, η 0 , and let the bottom topography (b) be used to ensure steady state as where g is the gravity constant and b 0 is a constant. With this basic state, performing the linear analysis would lead to a set of stability equations that depend on the position (x, y). By assuming that the disturbances are of small scale compared with that on which the zonal flow varies, that is, that terms involving perturbation variables times position variables (x or y) are neglected, this dependence disappears. As a consequence, simple calculations show that all that changes in the linear equations derived in BPT is that ζ 0 appears instead of f 0 . Therefore, the system with f 0 = 0, but with background vorticity ζ 0 = 0, is unstable and subject to exactly the same analysis as done in BPT for the height coordinate model. We clearly see that the Coriolis term is acting as a background vorticity, which seems to be necessary to trigger the instability. Also, one could interpret the relationships between the Rossby numbers and the instability growth rates simply as Rossby numbers relative to the background vorticity of the flow, not necessarily related to the Earth rotation.
Summarizing, although the Coriolis term plays no role in the lack of cancellation that leads to the existence of the instability, it can influence the instability as a source of vorticity in the vector-invariant momentum equations. Also, in the absence of rotation (f 0 = 0), the EEN scheme is still linearly unstable when the background state has non-zero vorticity and varies slowly between grid points.
In spherical geometry one cannot define a continuous constant vector field over the whole surface. Therefore, basic states used for linear analysis naturally have some source of divergence or vorticity. Simple basic states, such as the purely zonal flows that will be defined next in section 4, should have enough vorticity to trigger the instability even without the Coriolis term.

Test case
In this section we describe a zonal balanced flow test case that mimics small equivalent depth behaviours in spherical shallowwater models. This will be used to analyse the stability properties of several schemes in later sections. Shallow-water models based on either non-depth-weighted or depth-weighted Coriolis terms will be considered, as defined in section 2.2. The test is based on existing spherical shallow-water tests (Williamson et al., 1992) and is intended to be very simple to implement in existing codes. Consider as basic state a constant fluid depth η 0 , which is used to mimic small equivalent depths and may vary from millimetres to a few metres, and a zonal flow (u = u 0 cos(φ), v = 0), with maximum velocity given by u 0 = 2π a/12 days ≈ 38.6 m s −1 , where φ is the latitude and a = 6.371 km is the Earth radius. The bottom topography is then used to ensure a steady state, where g = 9.80616 m s −2 is the gravity constant and = 7.292 × 10 −5 rad s −1 is the rotation rate of the Earth. The test should be set up identically as in test case 2 of Williamson et al. (1992), except that the topography is used to balance the flow to become steady state, and the depth is defined to be constant (η 0 ). A good starting point for the constant depth is to assume η 0 = 1 m, but some models may require tests with smaller depths. As in the Williamson et al. (1992) test cases, all the experiments performed in this study do not use any additional perturbation to the initial conditions, as the numerical errors are enough to allow possibly existing unstable modes to emerge. Schemes that exactly represent the initial steady state may need a small perturbation in the initial conditions to allow an investigation of the unstable modes.
For the analysis of the Hollingsworth instability, this test may be run without the Coriolis term ( = 0), since the instability happens due to the lack of a discrete cancellation of advection terms solely.
The theoretical linear analysis for spherical models can be rather complicated, particularly on unstructured spherical grids. Also, numerical implementations of schemes for spherical shallowwater models rarely allow the possibility of running linear shallow water only. Nevertheless, it is not difficult to numerically investigate the eigen-structure of the most unstable modes with slight modifications of a nonlinear shallow-water code. In this work, we use a variation of the power method for small perturbations, which is fully described in Appendix B.
The output of the method is the value of the largest growth rate of a possibly existing unstable mode, which can be used to infer the e-folding time (time for the numerical solution to grow by a factor of e), and its associated eigenvector, which describes the spatial structure of the unstable mode.

Instabilities on quasi-uniform spherical grids
In this section we investigate the stability of spherical shallow-water models that adopt the vector-invariant form of the momentum equations and quasi-uniform spherical grids (Staniforth and Thuburn, 2012) considering small equivalent depths.

Analysis of finite-volume schemes on icosahedral grids
For recently developed schemes that use spherical unstructured grids, it is very difficult to derive formulations in which the advection term is decomposed into vortical and kinetic energy terms in a way that satisfies the cancellation property in the discrete sense. Considering mimetic finite-volume schemes for hexagonal-pentagonal grids (Voronoi grids), following the discretizations proposed in Thuburn et al. (2009) andRingler et al. (2010), hereafter named TRSK, no discrete cancellation is expected and the model is prone to being unstable. Shallow-water experiments with TRSK on the sphere (Ringler et al., 2010;Weller et al., 2012;Peixoto, 2016) did not reveal instabilities related to the non-cancellation issue of the vector-invariant momentum equations. Nevertheless, Skamarock et al. (2012) and Gassmann (2013) show that the instability indeed appears in 3D models. We will show here that the main point for the instability not to appear in the shallow-water models was that the growth rates of the instability were too small because of the large fluid depths adopted and the run-time scales tested. TRSK is well suited for icosahedral grid-based models. Icosahedral grids can be initially built based on a triangulation of the sphere (a Delaunay grid), but usually the dual grid (its Voronoi diagram) is adopted for the finite-volume computational cells (Staniforth and Thuburn, 2012). All test cases here use a dual Voronoi (hexagonal-pentagonal) grid of level 5, which has 10 242 computational cells (12 regular pentagons and 10 230 not necessarily regular hexagons) which corresponds to an approximate grid resolution of 240 km. This resolution may not adequately resolve the Rossby radius of deformation in most experiments performed below, but it is fine enough to investigate the stability of the scheme in a non-time-consuming way. Section 5.3 gives further discussion of this point.
The grids used here adopt a Spherical Centroidal Voronoi Tessellation (SCVT), which slightly modifies the original icosahedral grid cell nodes to ensure that they are very close to the centroid of the Voronoi cell that they define (Ju et al., 2011). In what follows, all figures that show spatial distribution of a scalar field on this grid will also show the Voronoi diagram associated with the basic triangular icosahedral grid, to serve as reference of the underlying grid structure. A four-stage fourthorder explicit Runge-Kutta was used as the time-stepping scheme with a time step of 400 s.
The original TRSK scheme is proposed to ensure total energy conservation (within time truncation error) and compatibility with the potential vorticity equation. To achieve this, the vorticity terms are weighted by the layer depth. We start the analysis with this original depth-weighted (DW) formulation and later show how the scheme behaves with a non-depth-weighted (NDW) formulation, which should be more closely related to the results of Skamarock et al. (2012) and Gassmann (2013), since they adopt level-based vertical coordinate systems.

Zonal balanced flow with no Coriolis force for depth-weighted TRSK
We start the analysis with the zonal balanced flow described in section 4, but without the Coriolis force ( = 0), just to illustrate that indeed the Hollingsworth instability can be triggered independently of a Coriolis force. Using the original TRSK scheme (depth-weighted) and a constant fluid depth of 1 m, the model is unstable and the linear analysis shows that the most unstable mode has a growth rate with e-folding time of 25.7 days, and eigenvector shown in Figure 1. We see that the mode is related to the underlying grid structure. For a depth of 0.1 m, the e-folding time reduces to 19.5 days, and for a depth of 0.01 m the e-folding time drops to 5.2 days. For large depths, such as 100 m or larger, the model does not blow up in run-times of up to 1 year. This example shows how the Hollingsworth instability exists independently of the Coriolis force. Nevertheless, since most atmosphere and ocean models do have the Coriolis term included, we will continue with further analysis using only the test with Coriolis force included.

Zonal balanced flow with Coriolis force for depth-weighted versions of TRSK
We will now consider the test of balanced zonal flow with rotation, as described in section 4. Adopting a constant fluid depth of 1 m, we show in Figure 2(a) the error in the depth field for the TRSK scheme after 14.5 days, which is a few time steps before the model blows up. The model is clearly unstable and Figure 2(b) shows the dominant eigenvector. Both the error and the eigenvector patterns seem to be related to geometric properties of the grid, which agrees with previous knowledge of grid imprinting often observed in this kind of model (Weller et al., 2012;Peixoto and Barros, 2013), but they do not match each other exactly. The method used for the linear analysis is sensitive to possibly very close eigenmodes, which may be hard to separate. Also, nonlinear effects may influence the patterns in the full model run. Nevertheless, both patterns seem to indicate that the unstable mode is stronger near edges of the original icosahedral dual grid.
The blue line in Figure 3, referring to the TRSK scheme, shows the e-folding times of the instability for various depths. The instability is noticeably stronger for smaller depths, almost reaching an e-folding time of approximately 1 day. For the depths adopted in standard shallow-water test cases (Williamson et al., 1992), which are of about 10 km, the growth rate is so small that one would not observe instabilities even in long runs (of many years). For the high-resolution 3D models discussed in section 2, which have very small equivalent depths, we are very likely to be in the region where the e-folding time is near to 1 day. This can have significant impact in medium-and long-range weather forecasts, and also in long climate simulations for either the ocean or the atmosphere. Gassmann (2013) analysed the instability problem on planar regular hexagons and did not find a way to modify the kinetic energy discretization to ensure exact discrete cancellation. Nevertheless, a discretization that minimizes the non-cancellation effect is suggested. In 3D models this scheme was shown to reduce the effects of the instability in standard baroclinic wave test cases (e.g. Skamarock et al., 2012;Gassmann, 2013). The stable results observed are of course very important to enable practical usage of such schemes, but no warranty of stability is established. In fact, it might be the case that other experiments, or longer run-time periods, reveal the instability.
We examined the scheme proposed by Gassmann (2013) (hereafter denoted as the GASS scheme) using the the zonal flow test case with Coriolis force. Considering 1 m constant depth, it blows up shortly after 5 days. Figure 4(a) shows the error in the depth field a few steps before blowing up. The errors dominate at the centre of the original pentagons of the icosahedral grid, so they are clearly connected to the grid structure. The linear analysis of this scheme, confirms the lack of stability, and the pattern of the dominant eigenvector is given in Figure 4(b) and also shows larger values near the centre of the original pentagons. We notice that the blow-up time in this case was sooner than that of the original TRSK scheme. This is confirmed by the calculation of the e-folding time, which is smaller for the GASS scheme ( Figure 3 for 1 m depth).
The above result is somewhat counter-intuitive compared with the previous results discussed by Skamarock et al. (2012) and Gassmann (2013) for 3D models, since in this scenario the modification seems to have worsened the stability of the scheme. There are two main points here. First, these 3D models can be more naturally classified as having non-depth-weighted vorticity terms, since the vorticity is weighted by density, not the layer depth. A non-depth-weighted version of the TRSK scheme will  be investigated later in this article. Also, both models use a Lorenz vertical grid, which, as discussed in section 2, is related to very small equivalent depths (much smaller than the 1 m tested here). For smaller equivalent depths, the GASS scheme has e-folding time larger than the original TRSK scheme, as one can see in Figure 3, noting that the two lines cross at about half a metre. For example, for a depth of 0.1 m, TRSK has e-folding of 1.56 days and GASS has 1.8 (the growth rate is approximately 15% slower with the GASS scheme).
The modified kinetic energy proposed by Gassmann (2013) is based on a linear combination of the original kinetic energy used in TRSK, defined at cell centres and which we will denote as K c , and a kinetic energy K v , calculated based on the kinetic energies obtained from the triangles surrounding the cell, which results in the total kinetic energy K = αK c + (1 − α)K v . The choice of α = 1 gives back the original TRSK scheme, and a standard choice of α = 0.75 is suggested to be stable, as it minimizes the cancellation error in the advection decomposition on the plane. The experiments shown so far all adopted this standard parameter of α = 0.75. Now we investigate the influence of this parameter choice with respect to the growth rates of the most unstable modes. Figure 5(a) shows how the e-folding times vary with α for four choices of equivalent depths. For a depth of 1 m, we notice that the choice of parameter of α = 0.75 seems to give faster growth rates (smaller e-folding times) than the original scheme (as already observed before), but the figure also points out that it is possible to choose the parameter in a range that reduces the growth rates (α near 0.9-0.95). For a smaller depth, of 0.01 m, the scheme using α = 0.75 indeed gives smaller growth rates, but the optimal parameter in this case would be close to α = 0.625, which is in fact what is adopted in the Model for Prediction Across Scales (MPAS; Skamarock et al., 2012). Experiments with even smaller depths seems to indicate that the choice of α = 0.625 is optimal for very small equivalent depths. Peixoto (2016) noticed that the gradient of the kinetic energy discretization used in the TRSK scheme was very inaccurate on unstructured spherical grids; in fact, it was proven to be inconsistent (0th order accurate), a problem also shared with the discretization of other terms of the TRSK scheme. Peixoto (2016) suggests modifications of the TRSK scheme to ensure at least overall first-order accuracy, at the cost of losing some mimetic properties. More accurate kinetic energy and vorticity discretizations could, in theory, reduce the mis-cancellation gap that occurs in the advective term. The cancellation in this case would be ensured asymptotically for sufficiently smooth fields. Nevertheless, stability issues are usually related to near-grid-scale features, so asymptotic cancellation might not be enough to avoid the instability. We investigated the stability of the consistent scheme proposed in Peixoto (2016), which we hereafter refer to as PXT. The test case using 1 m depth shows that the scheme blows up shortly before 3 days (Figure 6(a) for the pattern of the error in the depth field at 2.8 days). The growth rate is in fact somewhat larger than that of the original TRSK scheme, as illustrated by the smaller e-folding times in Figure 3 with label PXT.

Analysis of non-depth-weighted versions of TRSK
So far, all analysis done for the TRSK-based schemes has considered the depth-weighted form of the vorticity term, as in the original shallow-water formulation, and appropriate for a three-dimensional layer-based coordinate model. It is possible to remove the depth weighting of the vorticity, corresponding to a three-dimensional level-based coordinate, although some mimetic properties will be lost.
We analysed the non-depth-weighted version of TRSK and found it to be unstable. The error pattern of the depth field observed a few steps before blow-up assuming a mean depth of 1 m is shown in Figure 7. The calculated growth rate indicates an e-folding time of 1.85 days, which is smaller than the e-folding time of the layer version of TRSK (which has an e-folding time of approximately 6.5 days), so the instability grows faster in the non-depth-weighted version for this mean depth.
The GASS scheme can also be used without depth weighting of the vorticity term. In fact, the method as derived in Gassmann (2013) considered only the relative vorticity, not the potential vorticity. Also, in three-dimensional models such as MPAS (Skamarock et al., 2012), the vorticity is weighted by density rather than the fluid layer depth, and the behaviour is best captured by the non-depth-weighted shallow-water case. The non-depth-weighted version of the GASS scheme, using the original scheme parameter choice of 0.75, is still unstable, with e-folding times shown in Figure 5   GASS. Therefore the GASS scheme seems to be slowing down the growth speed of the instability with respect to TRSK. We show in Figure 5(b) how the parameter choice of the GASS scheme affects the e-folding time of the most unstable mode. In all our tests with the non-depth-weighted scheme, using the GASS scheme with the original parameter choice of 0.75 was always beneficial to reduce the growth speed of the instability. We also note that the modification of the GASS scheme seems to be more effective (larger e-folding times) for the non-depth-weighted scheme than for the depth-weighted scheme when compared to the original TRSK scheme (the TRSK scheme e-folding times can be observed in Figure 5(b) looking at the parameter choice of 1).
As a result from the shallow-water equation analysis for both depth-weighted and non-depth-weighted models, we conclude that the modification proposed by Gassmann (2013) does not seem to be enough to eliminate the instability. Nevertheless, it can delay its interference in the model by reducing its growth rate. The two 3D models for which the instability was recently reported (Skamarock et al., 2012;Gassmann, 2013) showed results based on a baroclinic wave test suggested by Jablonowski and Williamson (2006). The tests show a clear unstable mode at day 8 or 9 of integration for the original TRSK scheme. At day 9, Gassmann's scheme does not reveal traces of the unstable mode yet. From our analysis, apparently the instability has not yet grown significantly, due to the smaller growth rate, but it is prone to appear at later times. Further investigations to see if indeed the instability appears later in the baroclinic wave test case would be required to confirm the extension of these results for 3D models, but this is beyond the scope of this article.

Analysis of a finite-element mimetic scheme
Thuburn and Cotter (2015) proposed a finite-element numerical scheme for the shallow-water equations on a rotating sphere. It uses compound elements, which provide a generalization of the lowest-order Raviart-Thomas finite elements to arbitrary polygonal grids and give a finite-element analogue of the C-grid placement of variables. The finite-element scheme is closely related to the finite-volume scheme of Thuburn et al. (2014), which in turn is derived from TRSK. It has the same mimetic properties, and it uses the same semi-implicit time    integration scheme and the same finite-volume advection scheme for advection of mass and potential vorticity. As in the original TRSK scheme, this finite-element scheme also adopts a depthweighted vorticity in order to obtain compatibility between the linear shallow-water equation scheme and the discrete potential vorticity equation. The finite-element scheme has greater accuracy than TRSK, which comes at the price of inverting certain mass matrices and related operators at each time step. The discretizations are built starting from the vector-invariant form of the equations, and so are susceptible to non-cancellation effects of the momentum advection term decomposition. No numerical instabilities were detected within the experiments performed in the original article. Closely related numerical methods are under consideration for the next generation threedimensional atmospheric dynamical core at the UK Met Office. Therefore, insights about whether the Hollingsworth instabilities will happen in this case are highly desirable. We used our suggested test case of a very thin layer of shallow water (with Coriolis force) to investigate potential instabilities. Two grid possibilities were used: (i) a cubed-sphere grid, where we used a resolution with 13 824 quadrilateral cells, and (ii) an icosahedral based grid, formed by hexagons and 12 pentagons, where we used a resolution with 10 242 Voronoi cells. The icosahedral grid was modified with the approach suggested by Heikes and Randall (1995), and has an orthogonal dual grid, whereas the cubed-sphere grid used the same modification as adopted in Thuburn et al. (2014), and is non-orthogonal. The time step adopted was 400 s.
Our experiments showed that the model indeed reveals itself to be unstable under small enough mean fluid depth. For a mean depth of 1 m, the cubed-sphere model blows up shortly after 5.17 days (Figure 9(a)). For this same depth, the hexagonal grid did not blow up for runs of over 1 year, which shows that it is potentially less susceptible to the instability. Nevertheless, for a mean depth of 0.10 m, it blows up shortly after 5.95 days (Figure 9(a)).
An attempt to use the power method to estimate the growth rates was performed. We observed lack of convergence of the algorithm in the experiments performed, even with very small mean depths. The dominating vector seems to be strongly caught by grid structures during the iterations, and does not seem to reflect the actual dominant computational mode of the instability, so these results are not shown.

Equivalent depths and the Rossby radius of deformation
In the experiments discussed above, reducing the mean depth has the effect of reducing the Rossby radius of deformation    λ = √ gη 0 /f . It is well-known that C-grid schemes can behave poorly when the Rossby radius is not well resolved. For example, taking η 0 = 1 m, the Rossby radius at midlatitudes is approximately λ = √ gη 0 /f ≈ 30 km, which is far from being resolved on the grids used above with horizontal resolution of approximately 240 km. It is therefore legitimate to ask whether the poor resolution of the Rossby radius might be contributing to the observed instability, and whether the instability might be ameliorated by using much higher horizontal resolution typical of that in modern weather and climate models.
For the non-depth-weighted EEN and AL schemes on a plane, the analysis of BPT (their Figures 6 and 7) shows that the growth rate (non-dimensionalized using f ) in fact increases as resolution is refined, that is, for increasing grid Rossby number R u keeping Froude number F u fixed. Consistent with this, Ducousso et al. (2017) found that, in a version of the Nucleus for European Modelling of the Ocean (NEMO) model using the EEN scheme, the instability became significantly worse at finer horizontal resolution.
To examine this issue in our experiments, the test cases for depth-weighted TRSK without and with Coriolis force at η 0 = 1 m may be compared. In the case without Coriolis force, the Rossby radius, based on the the midlatitude absolute vorticity rather than f , is λ ≈ 360 km and so is marginally resolved. In the case with Coriolis force, on the other hand, the Rossby radius λ ≈ 30 km is badly under-resolved. The respective e-folding times are 25.7 days and 6.5 days. Thus, the absolute growth rate of the instability is greater for the case of badly under-resolved Rossby radius. However, the background absolute vorticity is 13 times as strong in the case with Coriolis force. Therefore, the growth rate non-dimensionalized by the background absolute vorticity is in fact greater for the case of marginally resolved Rossby radius.
We also repeated the test for depth-weighted TRSK, with η 0 = 1 m and including Coriolis force, on grids with 120 km, 60 km, and 30 km resolution. In all cases the model crashed after about 14.5 days, similar to the original 240 km run. The linear analysis indicates that the e-folding time at 240 km resolution is about 6.5 days. Increasing the grid resolution to 120 km and 60 km results in e-folding times respectively of 7.4 and 9.5 days. In this set of experiments the growth rate of the instability (both absolute and non-dimensionalized) shows a modest decrease as the resolution of the Rossby radius improves.
From this limited set of results it appears that the dependence of the instability on resolution of the Rossby radius may be quite complicated, and might be affected by the specific scheme used, depth-weighting, spherical geometry, and whether the background vorticity is provided by the flow or the planetary rotation. Nevertheless, it is clear that increasing the horizontal resolution does not, in general, suppress the instability, and may make it worse. The results also imply that a relatively coarse and computationally cheap horizontal resolution, as used in most of our tests, is adequate for diagnosing the presence of an instability.

Analysis of instabilities in ENDGame
The discretization of the atmosphere dynamical core of the Unified Model of the UK MetOffice, also known as ENDGame (Wood et al., 2014), is based on the shallow-water discretization discussed in Zerroukat et al. (2009). In this section, we will first perform a linear stability analysis of the Zerroukat et al. (2009) scheme and two alternatives, under a similar approach as adopted in BPT, for the planar version of the scheme. This shows that depth weighting of the Coriolis term on its own can give rise to instabilities, depending how it is applied. Then, we will show numerical results for the spherical shallow-water model. Implications for the stability of the three-dimensional ENDGame will be discussed at the end.

Original ENDGame scheme
Following the notation of BPT, consider the shallow-water equations written in advective form for a planar domain with velocity components (u, v), water height η and bottom topography (bathymetry) b, where the subscripts refer to partial derivatives, g is the gravity constant, f 0 is the Coriolis parameter (constant) and the free surface height is given by η + b.
The scheme proposed in Zerroukat et al. (2009) adopts a semi-implicit semi-Lagrangian approach, but we are interested in the instabilities that arise due to spatial discretizations related to geostrophic and inertia-gravity linear modes, so we will adopt a simplified advection scheme. The scheme is constructed on a usual horizontal C-staggered grid. The semi-discrete version of these equations for a planar domain of constant Coriolis parameter can then be written as where the over-lines indicate a centred averaging in the direction of the superscript and the δs indicate a centred differencing in the direction of the subscript (BPT give details on the notation).
The key point of this discretization is that it uses depth-  energy budget and have good Rossby mode dispersion properties (Thuburn and Staniforth, 2004).
The linearization will be taken with respect to a geostrophically balanced non-resting steady state. To ensure the flow is balanced, one could either think of adding a forcing to the equations, or else, and as we will adopt, use the bottom topography b to enforce the steady state. To do so, consider the constant basic velocity and mean water depth (u 1 , v 1 , η 0 ), and bottom topography defined as with b 0 constant. Clearly the bottom topography terms in the right-hand side of the discrete equations reduce to f 0 v 1 and −f 0 u 1 for the u and v equations, respectively. The linearized version of these equations may then be written for the perturbations variables (u where the over-lines with double superscript indicates averaging in both directions indicated. Assume that the perturbations are of a wave-like form, where ω is a non-dimensional frequency normalized using f 0 , and κ and λ are non-dimensional horizontal wavenumbers for the xand y-directions normalized using the grid spacings x and y respectively. Define the following convenient non-dimensional quantities respectively as the Froude numbers for the basic flows u 1 and v 1 , R c as twice the ratio of the Rossby radius (c/f 0 ) and the grid spacing y, where The grid-scale Rossby numbers R u and R v for the flows u 1 and v 1 can be constructed using the above parameters Substituting the wave-like form into the discrete linearized equations (11), and using the above relations together with the definitions, one obtains a matrix form of the stability problem, Although it is possible to derive analytical solutions for special cases, such as those aligned with the grid (κ = 0), these do not capture the most unstable modes, so direct numerical evaluations of the eigenvalues (ω) of the stability matrices were performed. We considered only the imaginary part of the eigenvalues, ω i , which is the part that leads to the instability, and present their maximum absolute values for a few parameter settings in Figures 10 and  11. These show that for a zonal flow (v 1 = 0), there are unstable modes that are stronger for larger horizontal wavenumbers (κ). The maximum growth rate increases as equivalent depth decreases (increasing Froude number) and also as horizontal resolution is refined (increasing grid Rossby number). Figure 11(b) shows that the most unstable modes are not necessarily aligned with the grid. Therefore, the planar version of the ENDGame formulation of the shallow-water equations may suffer from grid-scale instabilities. However, in comparison to BPT, their dispersion relation is not the same as that for Hollingsworth instabilities; in fact, as will be shown below, it is solely related to the discretization of the Coriolis term.

Alternative simplified scheme
The instability detected in the previous section is directly related to the form of discretization of the Coriolis term. This can be seen by noting that if the terms E u and E v in (17) were set to zero then the matrix of the system could be written as ( I + H), where H is Hermitian; would then be real-valued and the problem is neutrally stable. The terms E u and E v arise from the the last terms on the left-hand sides of the u and v equations in Eq. (11), which in turn arise from the depth weighting of the Coriolis term.
With this insight, it is simple to derive a stable discretization for this problem. For example, using simple uniform averaging results in E u = 0, E v = 0, and the scheme is then stable for this problem. This is not a new scheme; in fact it is the standard kind of discretization one would first apply in finite difference Cgrid models. It was frequently used in the early days of weather forecasting (e.g. McDonald and Bates, 1989) and is still considered in some models (e.g. Barros and Garcia, 2004). Nevertheless, this approach lacks some important properties of the scheme adopted in ENDGame. For example, in this scheme the Coriolis term is not neutral in the energy budget when used in a layer-based model. On an f -plane, Eq. (19) is obtained from Eq. (9) simply by removing the depth weighting.

Alternative energy-conserving scheme
Alternatively, one can build a depth-weighted discretization that still retains an energy conserving Coriolis term and is stable for the planar analysis: This is identical to the energy conserving scheme of Sadourny (1975), with the exception that here only the planetary vorticity (Coriolis parameter) is used instead of the absolute or relative fluid vorticity, since the momentum equations are not in vectorinvariant form. With this discretization, again it is possible to show that E u = 0 and E v = 0, so the scheme is neutrally stable. One intellectually satisfying aspect of this discretization is that f is evaluated at vorticity points, which feels more natural and is closer to the approach of Arakawa and Lamb (1981) than the original scheme, which evaluates f at mass/pressure points.

Analysis on the sphere 6.2.1. Description of the schemes
We have described three numerical schemes so far. First the original ENDGame scheme, described in Eq. (9), which we will denote as 'ORIG'. Second, the simplified scheme, described in Eq. (19), which we will denote as 'SIMP'. And third, an alternative energy-conserving scheme, analogous to Sadourny's energy-conserving scheme, described in Eq. (20), which will be denoted as 'ALTEC'.
On the sphere, considering variable Coriolis parameter (f = 2 sin φ) and the spherical metric terms, ORIG is written as where φ refers to latitude and λ to longitude coordinates on the sphere. This scheme ensures that the Coriolis term does not contribute to the energy budget and also has steady geostrophic modes on the f -sphere. Thuburn (2007) shows that it also has accurate representation of Rossby waves (β-effect), as the Coriolis parameter f is calculated at mass points. Unfortunately, the previous section shows that it is numerically unstable on an f -plane.
SIMP is calculated on the sphere as On the calculation of the energy budget, it is possible to show that no exact cancellation is expected, and therefore the Coriolis term might undesirably contribute to the energy budget. Based on Thuburn (2007), it is possible to show that it also has steady geostrophic modes and accurate representation of Rossby waves. Also, it is numerically stable on an f -plane.
(a) ENDGAME-ORIG Height deviation at 14.  ALTEC considers the Coriolis parameter calculated at vorticity points, and may be written on the sphere as Again this matches the energy-conserving scheme of Sadourny (1975). It is energy-conserving and has steady geostrophic modes but, due to the position of the calculation of the Coriolis parameter at vorticity points, the β-effect may not be as accurate as ORIG.

Numerical experiments
The numerical experiments will be based on the zonal flow test case with Coriolis force, as described in section 4. The experiments will be performed with the original configurations of the shallowwater ENDGame scheme, with semi-Lagrangian semi-implicit discretization with 256 × 128 grid points and time-step size of 1200 s. A simulation run with this configuration assuming η 0 = 1m reveals that ORIG is unstable. Figure 12(a) shows the error in the height field at 14.5 days, a few time steps before the model blew up (the figure shows the step where the model attained maximum error in the velocity greater than 10 m s −1 ). The spherical linear analysis shows that the dominant eigenvector has larger values concentrated over midlatitudes (Figure 12(b)), as happens with the error pattern observed a few time steps before blow-up.
The e-folding times of the instability for different mean fluid depths (η 0 ) are shown in Figure 14. For depths larger than 4 m the linear analysis method did not converge, indicating that either the model is stable, or that the methodology is not sufficiently accurate to capture the instability. Indeed, for such larger layer depths we did not observe instability within reasonable run time of the model (up to 1 year).
The most unstable mode lies around 60 • latitude, where, for a resolution of 256 × 128 in latitude-longitude coordinates, f 0 ≈ 7 × 10 −5 , x ≈ 78 km and y ≈ 156 km, so the ratio of grid spacings is X ≈ 0.5. Also, the velocity field is zonal with speed approximately 19 m s −1 at this latitude. Therefore, we may estimate the local planar Froude number F u ≈ 6 and the Rossby number R u ≈ 7. The linear analysis on the plane provides an estimate for these parameters of a maximum growth rate of approximately ωf 0 ≈ 6.5 × 10 −5 s −1 , which gives an e-folding time of about 4 h. This is pessimistic compared to our numerical estimates of e-folding time, which, for this scenario on the sphere, are of about 1 day (Figure 14), but not so different given the great deal of approximation. Also, this scheme has implicit diffusion within the semi-Lagrangian advection, which might be responsible for the slower growth of the instability. The simple averaging scheme for the Coriolis term (SIMP) was found to be stable in all experiments performed, as predicted by the planar linear analysis. The error of the scheme remains small throughout the experiment for long periods of time.
The alternative energy conserving scheme (ALTEC) was shown to be stable on the f -plane. Nevertheless, on the sphere, the variable velocity and combination of averagings involving the depth variable and cos(φ) factors in the Coriolis term breaks the exact cancellation observed on the plane. A simulation run assuming η 0 = 1 m reveals that the alternative scheme is also unstable on the sphere, but it takes much longer for the model to crash. Figure 13(a) shows the error in the depth field at 60 days, only a few steps before blow-up. The unstable pattern emerging before blowing up seems to be very close to the poles.
The dominant eigenvector (Figure 13(b)) of the most unstable mode also shows that the problem is mainly closer to the poles. A closer look at how the growth rates change with the mean depth ( Figure 14) reveals that indeed at 1 m height the e-folding time of the alternative scheme is far (about four times) larger than that of the original scheme, which is approximately the extra time required for this model to blow up compared to the original scheme. Figure 14 also shows that for larger depths (larger than 4 m), the original scheme is in fact more stable than the alternative. A similar situation may happen if the depth is too small (smaller than 0.1 m). So, except for the interval of [0.01, 4] mean depths, the original scheme has growth rates smaller than the alternative scheme.

Implications for 3D ENDGame
The three-dimensional operational ENDGame model (Wood et al., 2014) used at the UK Met Office uses ORIG but in a heightbased vertical coordinate with density times model layer depth as weights for the Coriolis term. In a height-based coordinate, the model layer depths are fixed and the local variations in density are relatively small. Therefore, as discussed in section 2.2, the appropriate shallow-water model to analyse the stability of 3D ENDGame is one with non-depth-weighted Coriolis terms, similar to Eq. (22) but with some cos φ weights. To date, no traces of this instability have been detected in the operational model, and indeed, the above analysis suggests the model should be stable. A version of 3D ENDGame with a Lagrangian vertical coordinate has been developed by (Kavčič and Thuburn, 2017;personal communication). Early versions of that model used a depth-weighted Coriolis term, but were found to suffer from an instability with short vertical and meridional scales that appeared in regions of strong zonal wind. Removing the depth weighting from the Coriolis term eliminated the instability. All of these symptoms are consistent with the analysis above for the ORIG scheme.

Discussion and concluding remarks
We have proposed the use of shallow-water models with very small layer depth as a useful means of investigating certain modes of instability of three-dimensional numerical models. We have suggested a shallow-water test case that is simple to implement, and also a straightforward variant of the power method to identify the most unstable mode. In using this approach it is important to identify which shallow-water scheme corresponds to a given three-dimensional model since stability can depend crucially on the details. This is especially true for depth weighting of vorticity or Coriolis terms, whose effects can be either stabilizing or destabilizing.
The results obtained for all models investigated here are summarized in Table 1.
The finite-difference energy-enstrophy conserving schemes (EEN and AL) were only investigated by BPT on the f -plane, and empirical results (e.g. Lazić et al., 1986) seem to corroborate that the schemes are unstable for level-based models with nondepth-weighted vorticity terms. The planar analysis indicates that these schemes are stable for depth-weighted vorticity terms, and no evidence of instabilities in layer-based models exists to these authors' knowledge. Further investigations would be required to confirm the stability properties of these schemes for spherical shallow-water models with depth-weighted vorticity terms, which may be done following the lines discussed in this article.
A key point shown here, complementing BPT, is that the Hollingsworth instability exists even without the presence of Coriolis terms (Earth rotation). The Hollingsworth instability is due to a lack of discrete cancellation of advection terms and seems to exist only in the presence of some background vorticity. That said, the Coriolis force can influence the Hollingsworth instability by simply acting as a source of vorticity.
The finite-volume schemes investigated, based on the TRSK scheme, all seem to be unstable, either with or without depth weighting of vorticity terms. These analysis have direct impact, for example, on models such as MPAS (Skamarock et al., 2012) or DYNAMICO (Dubos et al., 2015), that use this scheme. The modifications proposed by Gassmann (2013) seem to slow down the growth rate of the instability, and is particularly beneficial for non-depth-weighted (level-based) models that run in configurations that imply very small equivalent depths (e.g. use a Lorenz grid). For the scheme with the modifications proposed by Peixoto (2016), we see that ensuring more accurate calculations of the terms involved in the existence of the instability is not a sufficient condition to improve the stability properties of the scheme.
The finite-element scheme tested in this work shows strong instabilities, particularly on the cubed-sphere grid, with slower growth rates on the hexagonal grid. Further analysis on the stability properties of these and other finite-element schemes are currently under investigation and will be reported elsewhere.
An interesting point is how the different weighting of the vorticity term influences the stability of models based on the vector-invariant momentum equations. For a depth-weighted scheme, the depth field appears in the vorticity term but does not appear in the kinetic energy term. This suggests that the non-cancellation would be more pronounced in this case when compared to a non-depth-weighted formulation. Indeed, our results for TRSK-based schemes using very small equivalent depths (< 1 m) show that the e-folding times of the depthweighted formulations are smaller than the non-depth-weighted versions ( Figure 5), indicating that the instability grows faster in depth-weighted models. Intriguingly though, BTP shows that for FD schemes on an f-plane the opposite happens: the depthweighted model is stable, whereas the non-depth-weighted is not. The reason is that the stability of the depth-weighted scheme shown in BTP is not directly connected to a perfect cancellation, but related to the conservation of uniform potential vorticity (Appendix F of BPT). Further analytical inspection of the unstable modes on hexagons would be required to make a detailed connection between the results of BTP for the EEN scheme and the results shown in this article for TRSK.
The instability detected for the ENDGame scheme was surprising at first, since it is not in principle related to the Hollingsworth instability, as it does not use the vector-invariant momentum equations. However, it shows that the suggested test cases seem to be applicable in a more general sense. In this case the suggested test case revealed an instability that had not been anticipated, and led us subsequently to perform the analysis presented in section 6.1.
As a final point, we note that the investigation of numerical instabilities arising in shallow-water systems considering small equivalent depths seems to be not only of theoretical interest, but of practical importance. Non-idealized weather and ocean models indeed possess vertical modes corresponding to very small equivalent depths and may be subject to these instabilities. The main purpose of this article was not to show solutions to stability issues, but more to enlighten the investigation path with tools and better understanding of such instabilities.

A. Linear analysis for the EEN scheme with no Coriolis force
In this appendix we follow the analysis of BPT to show that, with a constant velocity basic state and null Coriolis force (f 0 = 0), the original EEN scheme is neutrally stable on a plane.
Using the same notation as in BPT, one may write the linearized form of the shallow-water equations for the height coordinate model, ignoring the Coriolis term, as (Eqs (61-62) of BPT) ∂v ∂t +gδ y η +u 1 δ x v xν E +v 1 δ y v y +u 1 δ y u x −u xν E = 0, ∂η ∂t +u 1 δ x η x +v 1 δ y η y +η 0 δ x u +δ y v = 0, where the basic state has constant velocities (u 1 , v 1 ) and constant depth η 0 and the system refers to the perturbed variables for velocity and depth (primed variables). Let us assume wave-like solutions as where ω is the frequency and (κ, λ) are non-dimensional horizontal wavenumbers for the x-and y-directions normalized using the grid spacings x and y respectively. Substituting the wave-like forms in the perturbation equations, and using that, for a quantity ψ, the x-and y-averaging operators and the differencing operators δ x and δ y give, respectively, one obtains the linear system ⎡ where E 11 = 2u 1 x s κ c κ + 2v 1 y s λ c λ μ E , As in BPT, and c, c κ , c λ , s κ and s λ are defined in Eqs (14) and (16). It will be convenient to introduce Using these definitions with (A6) one sees that E 11 = E 33 − T v s λ , E 22 = E 33 − T u s κ , ( A 9 ) and the matrix (A5) becomes ⎡ ⎢ ⎣ in which is the Doppler-shifted non-dimensional frequency of the perturbation and X = x/ y. By direct calculation of the determinant D M of the matrix in (A10), one finds that where P = T u s κ + T v s λ , The solutions of D M = 0, = ±Q and = −P, are real-valued, so the system is neutrally stable.

B. Linear analysis on the sphere
In this appendix we describe an algorithm to evaluate the eigen-structure of the most unstable modes using a nonlinear shallow-water code modified to apply a variation of the power method to small perturbations. Let x (k) be the vector of values representing the model state at time step k and let G represent the action of the nonlinear shallow-water model integration scheme over one time step, so that the unmodified shallow-water model obeys We are interested in the evolution of perturbations to some basic statex, which should be steady. In practice, for a given numerical method it may be difficult or impossible to find a state close to the desired basic state that satisfiesx = G(x). Therefore a constant forcing term F =x − G(x) is introduced to compensate for the numerical drift of the desired basic state; F is easily computed by taking one model time step from the statex. The model governing equation (Eq. (B1)) is then replaced by Clearlyx is now a steady solution for the model (B2). Note that F must be a constant forcing, not a relaxation back towards the basic state; such a relaxation would damp perturbations and so affect the diagnosed eigenmode growth rates. Now consider the evolution of a small perturbation y (k) to the basic state, so that x (k) =x + y (k) . Linearizing (B2) gives where G (x) is the Jacobian matrix of the model evolution operator evaluated for the statex, and O(2) denotes higher-(second-)order terms in y (k) . It is the dominant eigenvalue and eigenvector of G that we wish to determine. Provided y remains small, taking repeated model time steps according to (B2) will cause y to evolve according to (B3) with O(2) negligible, and y should evolve towards the dominant eigenvector of G , as in the power method. However, with this method, it is likely that y fails to remain small and so the linearization breaks down before the dominant eigenvector emerges. Therefore, another modification is needed to the model to rescale perturbations to ensure they remain small. Thus we take a preliminary step forward using the model with constant forcing x * = G(x (k) ) + F, ( B 4 ) diagnose the perturbation then rescale the perturbation before computing the state at the new time level Here α k+1 = / r (k+1) , for > 0 a small constant. Iterating (B4)-(B6) should then determine the same eigenvector as the power method applied to G . If the method converges, then the absolute value of the dominant eigenvalue (λ) may be obtained using the converged value of α k , which we will denote simply as α, as λ = 1/α.
Provided the basic state is fluid-dynamically stable, the shallowwater equations should not have any growing modes. Then, for a stable numerical scheme, it is expected that λ = 1 for any parameter choice. The unstable modes are detected when λ > 1 (α < 1), and, in this case, the associated eigenvector will be given by the converged vector y.
The growth rate (ν) may be then obtained observing that λ = e ν t , where t is the time step used in the calculation of G. Consequently, the e-folding time may be also calculated directly from ν.
In our experiments we adopted = 1 × 10 −5 and the 2 norm of the velocity field for calculation of α k . A small local perturbation was added to the initial height field to trigger any unstable modes. This approach can usually be easily incorporated in standard shallow-water model codes.
Some important points must be made with respect to the application of this scheme and the interpretation of its results. First, we are assuming that the forcing F is small. For the basic states used in this paper, which are all in analytical balance, F represents the numerical error between the numerical adjustment with respect to the analytical steady state, and therefore should be limited to local truncation errors, and thus small. Second, the eigenvalues might not necessarily be well separated, which means that the eigenvector obtained might not always match the dominant pattern appearing in a full model run. This lack of matching might also happen due to nonlinear effects influencing the full model run. Therefore, if the method converges, we have as a result one of the possible dominant eigenmodes, which is enough to show that the model is linearly unstable and will provide estimates of the instability growth rate.