Analysis of and Solution to the Polar Numerical Noise Within the Shallow‐Water Model on the Latitude‐Longitude Grid

This study conducts an analysis of the polar numerical noise in the barotropic shallow‐water version of the Grid‐point Atmospheric Model of IAP LASG (GAMIL‐SW) and provides a good solution to the problem. GAMIL‐SW suffers from numerical noise in the polar region in some ideal test cases, which is likely to be detrimental to the full physical model. The noise is suspected to be related to the nonlinear advection term in the momentum equation. Thus, a new shallow‐water model with a vector‐invariant form of the momentum equation is developed on the latitude‐longitude grid to analyze the polar noise. It is found that the version with meridional wind component staggered on the pole is free from noise, while the version with zonal wind component staggered on the pole is still contaminated. By redefining the polar relative vorticity, the polar noise is eliminated in the latter version, and the global conservation properties are maintained. In addition, the test cases demonstrate that the new shallow‐water model maintains the properties of the original GAMIL‐SW with respect to numerical accuracy and computational stability. This study helps to identify appropriate governing equations to further develop the next generation of GAMIL dynamical core.


Introduction
Atmospheric general circulation models (AGCMs) are one of the crucial tools for operational numerical weather prediction and climate modeling and are one of the fundamental components in Earth system models (ESMs) (Wang et al., 2009). At the heart of AGCMs is the dynamical core, which is responsible for solving the governing equations of atmospheric dynamics and thermodynamics on the resolved scales (Thuburn, 2008). Along with the advances in high-performance computing, substantial investments are being made in the development of the next generation of global nonhydrostatic high-resolution numerical models at modeling centers around the world (e.g., Kühnlein et al., 2019;Satoh et al., 2014;Skamarock et al., 2012;Ullrich et al., 2017). The Grid-point Atmospheric Model of IAP LASG (GAMIL) is an AGCM based on the finite-difference scheme, developed by the State Key Laboratory of Numerical Modeling for Atmospheric Sciences and Geophysical Fluid Dynamics (LASG) in the Institute of Atmospheric Physics (IAP) of the Chinese Academy of Sciences (Wang & Ji, 2006;Wang et al., 2004). GAMIL is the atmospheric component in the coupled climate system model Flexible Global Ocean-Atmosphere-Land System model: grid-point version (FGOALS-g) (e.g., Li et al., 2013;Yu et al., 2008). However, as this model remains to solve the hydrostatic primitive equations, it is urgent to upgrade the model to a nonhydrostatic version for higher resolution global weather forecast and climate modeling. In general, in the process of developing atmospheric models, model developers tend to start from the shallow-water equations, as they mimic the important features of the horizontal aspects of the dynamics (e.g., Staniforth & Thuburn, 2012;Thuburn, 2008;. Many formulations of the shallow-water equations on a rotating sphere are available. These formulations are equivalent in the continuous case but lead to very different discretization forms. Moreover, the shallow-water equations have five basic conservative integrals, including three primary integrals (total mass, total absolute angular momentum, and total potential vorticity [PV]) and two quadratic integrals (total energy and total potential enstrophy) (Wang & Ji, 2003). If we want to simulate the continuous adiabatic frictionless governing equations faithfully, numerical methods should be applied to approximately conserve these integrals, particularly the quadratic ones, which is the necessary condition for avoiding computational instability in a nonlinear system (Arakawa & Lamb, 1977, 1981Sadourny, 1975).
GAMIL uses the IAP variable transformation method (Zhang, 1990;Zeng & Zhang, 1987;Zeng et al., 1989) to convert the standard primitive equations to the quadratic conservation form, which ensures the conservation of total effective energy in the hydrostatic model under spatial finite-difference discretization. The single-layer shallow-water version of GAMIL (GAMIL-SW) has many of the same properties as its hydrostatic version, such as antisymmetric spatial operators, and total mass and total energy conservation. Consequently, various idealized test cases for the shallow-water model (e.g., Bates & Li, 1997;Ferguson et al., 2019;Galewsk et al., 2004;McDonald & Bates, 1989;Nair et al., 2005;Shamir & Paldor, 2016;Williamson et al., 1992) can be applied to evaluate the basic numerical schemes from multiple perspectives, for example, the model grid, prognostic equations, variable placement, and numerical method. Williamson et al. (1992) proposed a standard suite of test cases for verifying numerical algorithms for solving the shallow-water equations on a sphere. We evaluated GAMIL-SW with the test case 5 of the suite and found that the polar noise is noticeable in this case. The unphysical noise is at the grid scale and propagates from the poles to the midlatitudes. The noisy solution is likely to be detrimental in three-dimensional simulations when the physical processes are coupled within the dynamical core, for example, via advection and wave propagation. Therefore, this problem is critical and needs to be solved by designing suitable numerical methods for the dynamical core.
The objective of this study is to analyze the polar noise and provide a solution to the problem. A number of approaches have been attempted without success. The major efforts in the early phase are in twofold. First, averaging in the calculation of Coriolis terms due to the use of the staggered Arakawa C-grid was suspected, which may lead to grid-scale oscillations when the Rossby deformational radius is underresolved, but this is not the case for the model resolution in this study. Second, we endeavored to add artificial diffusion or damping, a polar fast Fourier transform (FFT) filter, a digital filter or a nonlinear diffusion to remove the noise. However, none of these approaches have been found to be effective. Nevertheless, during the investigations, two features were noticed. One is that the noise would disappear if the nonlinear momentum advection terms are switched off, and the other is that the abrupt increase of the potential enstrophy may result from the noise. Therefore, the source of the noise is probably related to these two phenomena. On the basis of these observations, we turn to the vector-invariant form of the horizontal momentum equation, which bypasses the potential problematic discretization of nonlinear horizontal momentum transport (Skamarock et al., 2012). In addition, the spatial discretization can be designed to dissipate the potential enstrophy (i.e., the square of the PV) without violating the total energy conservation (e.g., Ringler et al., 2010).
The remainder of this paper is organized as follows. Section 2 introduces the shallow-water equations in vector-invariant form, and the spatial discretization of the equations on the latitude-longitude (lat-lon) grid is described in detail in section 3. Section 4 mainly discusses the polar noise in the numerical experiments, and section 5 provides a summary and conclusion.

The Shallow-Water Equations in Vector-Invariant Form
The standard nonlinear shallow-water equations, including the mass continuity equation in flux form and momentum equation in advection form, can be expressed in vector form as where h is the fluid thickness, u is the fluid velocity vector with components u and v in the longitudinal (λ) and latitudinal (φ) directions, respectively, and k is the local unit vertical vector. The other three parameters are the Coriolis parameter f ¼ 2Ω sin φ, gravity acceleration g, and bottom topography h s , in which Ω is the rotation rate of the Earth. The momentum advection term can be reformulated into a vector-invariant form by the so-called Lamb transformation (e.g., Gassmann & Herzog, 2008;Zängl et al., 2015) u where ξ=k·(∇ × u) is the relative vorticity. The first term of RHS in Equation 3 can be combined with the Coriolis term (fk × u) into the nonlinear Coriolis term If the PV is introduced, the nonlinear Coriolis term can be rewritten as the nonlinear PV flux where q=(ξ + ƒ)/h is the PV defined as the ratio between the absolute vorticity and the fluid thickness. As in Peixoto et al. (2018), the two expressions of nonlinear Coriolis term can be referred to as the nondepth-weighted form and depth-weighted form, respectively. They are equivalent in the continuous form but are different for discretization. The depth-weighted form is flexible for designing a numerical scheme of potential enstrophy conservation or dissipation by introducing the PV (e.g., Arakawa & Hsu, 1990;Arakawa & Lamb, 1981;Ringler et al., 2010;Sadourny, 1975;Takano & Wurtele, 1982). The nondepth-weighted form has been studied in some shallow-water models (e.g., Bonaventura & Ringler, 2005;Lin & Rood, 1997;Ringler & Randall, 2002;Tomita et al., 2001;Wang & Ji, 2003) and adopted in certain baroclinic models (e.g., Lin, 2004;Skamarock et al., 2012;Zängl et al., 2015). In this study, the depth-weighted form is used in the new shallow-water model. By substitution, the vector-invariant form of the momentum equation can be written as where K=u 2 /2 denotes the horizontal kinetic energy. hu and hu ⊥ =hk × u denote the normal and tangential mass flux, respectively. The second term in the momentum equation (Equation 6) involving q does not contribute to the change of total kinetic energy. Through this term, however, potential enstrophy can be damped by employing a diffusive advection scheme for PV, such as the anticipated PV method (APVM) or the scale-selective dissipation method (Chen et al., 2011;Sadourny & Basdevant, 1985). Arakawa and Lamb (1981) designed a finite-difference scheme for the shallow-water equations that simultaneously conserves mass, energy, and potential enstrophy using the vector-invariant form of momentum equations with the depth-weighted nonlinear Coriols term. They began with a general finite-difference discretization of the shallow-water equations on the staggered Arakawa C-grid. The discretization contained many unspecified coefficients, and these coefficients were determined via direct algebraic manipulation under the constrains that the scheme conserves the total energy and potential enstrophy under spatial discretization. The scheme had been extended to fourth-order accuracy for the vorticity in the nondivergent limit (Takano & Wurtele, 1982) and was generalized to a energy conserving and potential enstrophy dissipative system by Arakawa and Hsu (1990). However, as described in appendix 3 of Arakawa and Lamb (1981), some modifications of the scheme must be made if it is applied on the spherical lat-lon grid, such as the definition of kinetic energy near the poles, as well as the grid interval must be modified to satisfy some requirement of geometrical areas. Meanwhile, whether the scheme is suitable for the placement of u at the poles is unknown, since only the scheme for v at the poles on the lat-lon grid was presented in the literature. As will be described in section 3, the explicit time step with v at the poles is approximately equal to half of that with u at the poles, which is a disadvantage for v at the poles.
Cotter, 2012; Weller et al., 2012). In Thuburn et al. (2009), the linearized version of continuity equation (Equation 1) and momentum equation (Equation 6) was analyzed to derive a discretization for the Coriolis term on an arbitrarily structured orthogonal C-grid, which allows for representation of the stationary geostrophic modes in the linearized equations. The main result of Thuburn et al. (2009) is a set of weights for diagnosing the tangential velocity from the surrounding normal velocities. Ringler et al. (2010) further extended this method to the nonlinear shallow-water equations on the spherical centroidal Voronoi tessellation (SCVT) mesh and formulated a spatial discretization that conserves the total energy and PV, as well as conserves or dissipates the potential enstrophy. Weller et al. (2012) implemented the same TRiSK method for solving the shallow-water equations to investigate the computational modes and grid imprinting on five quasi-uniform spherical grids. However, although the TRiSK method is applicable to a wide variety of meshes theoretically, previous studies have utilized only quasi-uniform meshes, and the numerical accuracy of the method degrades to only first order  due to some unwanted features of those grids; for example, the dual edges do not bisect the primary edges perpendicularly, and there are several pentagon cells different from the surrounding hexagon cells (e.g., Ringler et al., 2010;Skamarock et al., 2012;Weller et al., 2012). Motivated by these considerations, in this study, the shallow-water model with equations in vector-invariant form is pursued on the regular lat-lon grid to solve the problem of polar noise and achieve second-order accuracy.

Spatial Discretization of Equations
The Arakawa C-grid has good dispersion properties of the inertia-gravity wave and gives an accurate representation of the geostrophic adjustment process provided the Rossby deformational radius is well resolved (Arakawa & Lamb, 1977). The prognostic variables are discretized with the C-grid staggering as illustrated in Figure 1. The horizontal velocity normal to the cell edge (u and v) as a pointwise value is prognosed at the cell edges. The fluid thickness h i,j as a cell-averaged value is prognosed at the primal cell center (i,j). All vorticity-related variables, including relative vorticity ( ξ i − 1 , and PV (q i − 1 2 ; j − 1 2 ), are diagnosed or mapped from prognostic variables on the primal cell. In view of computational efficiency, all the time-independent values, including area, length, and distance, are precalculated and stored before the time integration begins. The spatial discretization is a mixed finite-volume/finite-difference scheme, and the main calculations are described in the following subsections.

Divergence Operator
The continuity equation (Equation 1) is in flux form that guarantees mass conservation in the discrete system of equations; the mass-flux divergence can be given by using the Gauss divergence theorem over a primal cell where v represents either u or v component of velocity at the edge point and positive flux is outward, A j and l e denote the area and the edge length of the primal cell, respectively. The thickness at the cell edge ĥ e needs to be interpolated from the primal cell centers. The velocity divergence operator has second-order accuracy on the lat-lon grid because the velocity is centered on the primal cell edge, which is not always the case for the Voronoi mesh. As a result, the accuracy of the mass flux divergence operator is determined by the interpolation approximation of the thickness at the cell edge. Midpoint interpolation ĥ e ¼ h i þ h i þ 1 2 is used on a Voronoi grid in Ringler et al. (2010); that is, the edge weight is 1/2. It has second-order accuracy because the edges bisect the lines between the Voronoi generating points.
However, the midpoint scheme decreases the interpolation accuracy for the edge thickness on the lat-lon grid. The edge weight used for the interpolation is not equal in the zonal and meridional directions. In the zonal direction, the two triangles at the west and east side of one longitude line are symmetric and equal to each other; for example, A e þ j þ 1 is equal to A e − j þ 1 . However, in the meridional direction, the area of the rectangular mesh changes with respect to latitude; thus, the two triangular areas at the north and south side of one latitude line are not equal; for example, . Given the triangular area as weight, 10.1029/2020MS002047

Journal of Advances in Modeling Earth Systems
LI ET AL. Weller et al. (2012) proposed an alternative interpolation for the non-Voronoi grid, which ensures the second-order accuracy of conservative mapping between primal and dual meshes. For example, in Figure 1, where A e þ . In addition, it should be noted that the two triangles in the zonal direction within one primal cell are exactly spherical triangles because each edge of the triangle is one part of a great circle line, for example, A e − j þ 1 and A e þ j þ 1 , but the two triangles in the meridional direction are not spherical triangles because one edge of the triangle is a part of the latitude line that is not a great circle except the Equator, for example, A e − j þ 1 2 and A e þ j þ 1 2 . Thus, care must be taken to calculate the triangular areas in the spherical geometry on the lat-lon grid.

Discretization of Normal Gradient
In the momentum equation (Equation 6), to calculate ∂u/∂t, one needs to discretize the gradient of kinetic energy and geopotential on the primal cell edge. On the lat-lon grid, the edge of the primal cell is perpendicular to and bisected by the edge of the dual cell; hence, the second-order central finite-difference method can be implemented straightforwardly. For example, on the edge points i þ 1 2 ; j and i; j þ 1 2 , the gradients of geopotential in the zonal and meridional direction are calculated as where n λ and n φ denote the unit normal vectors in the zonal and meridional direction, respectively. d e j is the zonal distance between h i,j and h i+1,j ; d e j þ 1 2 is the meridional distance between h i,j and h i,j+1 . This simple two-point central finite difference on the lat-lon grid ensures that the gradient has second-order accuracy and is curl free around vertices.
is the overlapping area between the primal and dual cells, and A e j þ 1 (A e j þ 1 2 ) is the unique area associated with each edge e, which is the sum of two triangles on either side of the edge. All areas A and edge lengths l e are calculated in the spherical geometry, and the distance between neighboring primal cells is set by d e =2A e /l e .

Discretization of PV
In the momentum equation (Equation 6), to calculate ∂u/∂t, one also needs to interpolate q,h, and u ⊥ onto the u grid. This subsection describes the interpolation of q and h; how to calculate the tangential u ⊥ is to be described in section 3.6. PV is defined as absolute vorticity divided by the thickness within the shallow-water model. On the staggered C-grid, it can be calculated at either the primal cell center where thickness is located or at the dual cell center where relative vorticity is naturally calculated by applying the circulation theorem to the surrounding u and v. Lin and Rood (1997) defined PV on the primal cell to ensure that PV is accompanied by a valid thickness equation. Sadourny (1975), Arakawa and Lamb (1981), and Ringler et al. (2010) defined PV on the dual cell instead of the primal cell to obtain the PV properties of compatibility and consistency. To avoid the creation of a null space in the divergence field (Ringler et al., 2010;Skamarock, 2008), the PV is defined on the dual cell center in this study. First, the relative vorticity is discretized by applying Stokes' circulation theorem over a dual cell. For example, in Figure 1, the relative vorticity for dual cell center ði þ 1 2 ; j þ 1 2 Þ can be calculated by  Figure 1, On the lat-lon grid, the overlapping area between the primal and dual cells is equal in the zonal direction; thus, the same area weight A v þ j is used for both h i−1,j and h i,j , and the same area weight A v − j þ 1 is used for both h i−1,j+1 and h i,j+1 . Moreover, the area of the primal cell, dual cell, and the intersection are accurate calculations rather than approximate ones in spherical geometry. Consequently, the sum of four overlapping areas in one primal cell is exactly equal to the area of the primal cell, and the sum of four overlapping areas in one dual cell is exactly equal to the area of the dual cell. As a result, this interpolation scheme ensures that the divergence field on the dual mesh is compatible with the divergence field on the primal cell. In addition, on the lat-lon grid, the thickness over the primal and dual cells is located at the center of the primal and dual cells. These constraints on the lat-lon grid ensure that the interpolation operator has second-order accuracy .
After the calculation of PV on a dual cell, the next step is the interpolation of PV to the edge point from the dual cell for the nonlinear Coriolis force term. Two interpolation algorithms are used in this study. One is the midpoint scheme, which allows for the conservation of potential enstrophy. Another is APVM (Ringler et al., 2010;Sadourny & Basdevant, 1985), which leads to potential enstrophy dissipation as it is an upwind-biased estimate of the edge PV and provides an enstrophy sink.
where q v1 and q v2 are the PV at two ends of one primal edge, ∇ e q ¼ ð ∂q ∂x ; ∂q ∂y Þ is the gradient of q at edge point, and δt is the time step. On the lat-lon grid, for the u(v) point, ∂q ∂y ð ∂q ∂x Þ can be directly calculated using the two endpoint PVs of one edge with the finite-difference method, while ∂q ∂x ð ∂q ∂y Þ needs to be interpolated from surrounding values that are calculated on the edge point. The interpolation method is the same as that of tangential velocity described in section 3.6. With respect to other kinds of APVM, refer to Weller

10.1029/2020MS002047
Journal of Advances in Modeling Earth Systems (2012). In this study, the midpoint and APVM methods are applied to obtain potential enstrophy conservation and dissipation, respectively.

Calculation of Kinetic Energy
In the momentum equation (Equation 6), if the kinetic energy is defined on the primal cell centers, the calculation of the kinetic gradient has second-order accuracy using the central finite-difference method. The expression of kinetic energy in the discrete system is determined by the constraint of conservative exchange between potential and kinetic energy. Therefore, based on the calculation of thickness on the edge (Equation 8), where the triangular area at either side of each edge is the weight, the kinetic energy in terms of the normal velocities is derived (see Appendix A). For example, in Figure 1, This definition is the same as Weller et al. (2012) and more suitable for non-Voronoi grids. There are two points to note for the calculation of kinetic energy on the primal cell. First, the distance between cells is set by d e =2A e /l e instead of direct calculation in spherical geometry. In other words, the area associated with edge e and the edge length l e are calculated in spherical geometry before the distance d e is determined. Second, to ensure energetic consistency, the same area weight must be used in the interpolation for thickness from cell centers to edge points.

Discretization of the Nonlinear Coriolis Term
One of the requisites for total energy conservation is that the Coriolis force neither creates nor destroys kinetic energy as the Coriolis force is always orthogonal to the velocity vector. Under this constraint, the nonlinear Coriolis force term (or PV flux) is constructed between the target and surrounding edges (Ringler et al., 2010). On the lat-lon grid, for the PV flux on i þ 1 2 ; j

Journal of Advances in Modeling Earth Systems
where d e is the distance between two neighboring primal cells that share edge e, w ee ′ is the interpolation weight to be specified in section 3.6. l e ′ denotes the edge length of surrounding four edges. For example, as shown in Figure 2, in order to calculate the PV flux qhu at i; j þ 1 2 , one must first calculate the PV values at this target point and the surrounding four u points. One symmetric formulation of e q ee ′ is ðe q e þ e q e ′ Þ= 2, where e q e and e q e ′ are the PVs at the target point and the surrounding points, respectively. These PV values are remapped from PVs at the primal vertices (see section 3.3). In this way, the symmetry of the two e q terms along with the remapping weight w ee ′ following Thuburn et al. (2009) ensures that the Coriolis term is energy conserving. An alternative non-symmetric formulation of e q ee ′ is q e , which is directly defined on the edge using one interpolation scheme for the PV (see section 3.3). This formulation does not guarantee the nonlinear Coriolis force is energetically-neutral, but it would give an accurate treatment of PV from a potential enstrophy perspective (Ringler et al., 2010;Thuburn et al., 2014). In summary, there are two formulations of symmetry and nonsymmetry for calculating the nonlinear Coriolis term as well as two methods of midpoint (Equation 13) and APVM schemes (Equation 14) for calculating the edge PV. In this case, four combinations can be made for different purposes as given in Table 1. In the following analysis of the polar noise, we will be concerned with the energy-conserving and potential enstrophy-dissipating scheme along with the potential enstrophy-conserving scheme.

Calculation of Weights and Tangential Velocity
The weight w ee ′ used in the previous subsection is a key parameter because it is used to calculate the PV flux as well as the tangential velocities (vector quantities along edges). The weight w ee ′ is derived from the sum of area fractions for each cell vertex where A iv is the overlapping area between the dual cell around vertex v and the primal cell i, and the vs are the vertices in a walk between edge e and e′. The weights are proportional to the overlap area. For the detailed computational method, refer to Thuburn et al. (2009) andWeller et al. (2012). As Thuburn et al. (2009) noted, the weight can also be set to 1/4 on the lat-lon grid. The two options are provided in the new shallow-water model, and the experimental results are similar; thus, 1/4 is adopted at present.
The tangential velocity can be reconstructed from neighboring normal components at the edges of the cells.
The interpolation needs to ensure that the divergence of vector field on the dual cell is a convex combination of the divergence on the primal cell: where u e ′ denotes the normal velocity and e′∈ECP(e) denotes the four edges belonging to the two primal cells that share edge e. Specifying the interpolation operator on the lat-lon grid, Figure 2 illustrates how to interpolate neighboring normal wind components (u and v) to the tangential wind components (u ⊥ and v ⊥ ). With respect to the velocity flux and spherical coordinates, the tangential velocity flux is averaged by the surrounding four normal velocity fluxes, for example, as given by Equations 20 and 21. By derivation, u ⊥ is the 4-point arithmetic mean from surrounding us, but for v ⊥ , the interpolation needs to contain geometrical factors (cosφ) due to grid intervals varying with latitude. The tangential velocity is calculated as follows

Journal of Advances in Modeling Earth Systems
and can be simplified as where the overline with superscript stands for an equally weighted 2-point average in the spatial direction (λ or φ) on the lat-lon grid (e.g., Thuburn & Staniforth, 2004;Thuburn et al., 2009).

Calculation of PV on the Pole
As the staggered Arakawa C-grid is applied on the lat-lon grid, two different placements of prognostic variables at the poles are available (e.g., Thuburn, 2016;Thuburn & Staniforth, 2004), that is, u at pole and v at pole, as shown in Figures 3a and 3b. Because of the two singular points on the spherical polar coordinate system, the row of grid cells associated with each pole are treated triangles. These two placements result in different calculation of PV near the poles. The first placement needs to calculate one circle of PV at half a grid length to the pole, while each relative vorticity is calculated in one small sector with the u and v wind components using Stokes' theorem. The calculation is similar to the first generation of European Centre for Medium-Range Weather Forecasts (ECMWF) finite-difference gridpoint model with the same vector-invariant form of the momentum equation (e.g., Burridge, 1980;Simmons et al., 1989). However, the calculation of nonlinear Coriolis term in their discretization needs the "mirror value" of PV beyond the poles and the mass zonal flux at the poles. These auxiliary calculations are necessary to conserve energy and potential enstrophy or conserve potential enstrophy but not energy. The second placement of v at pole requires calculation of only one value of relative vorticity at the pole with a ring of u values closest to the pole using Stokes' theorem. The calculation is similar to Arakawa and Lamb (1981) with respect to the vorticity at the pole, but the kinetic energy and the grid meridional distance near the poles must be redefined in Arakawa and Lamb (1981) to conserve both potential enstrophy and energy.
One important issue associated with the two variable placements at the poles is the size of integration time step, because the maximum numerically stable time step for the finite-difference scheme is mostly limited by the zonal grid spacings on the lat-lon grid. Assuming the same Courant-Friedrichs-Lewy (CFL) number, CFL=cΔt/Δx, the time step Δt will become smaller with decreasing the zonal grid spacing Δx under the condition of the same wave speed c. Therefore, with respect to the time step, the arrangement with v located at the pole is approximately half of that with u at pole, which is a disadvantage of v at pole. For this reason, the u at pole is preferred in general. In the following tests, however, it is found that u at pole exhibits polar noises as in GAMIL-SW; therefore, another PV calculation method is casted. The circle of relative vorticity adjacent to the poles is calculated using only the circle of u nearest the poles, as shown in Figure 3c. The fluid thickness h is remapped from the surrounding thicknesses to calculate the PV at the primal vertex. Thus, the circle of PV next to the poles is different from each other only due to the different h eventually. In other words, the v nearest the poles is excluded from the calculation of PV, which imitates v at pole. To further check this modification, two sensitivity experiments are conducted in section 4.1 for flow over mountain test case. Another major difference between u at pole and v at pole is the mismatch in the degrees of freedom (DOFs) of mass and velocity. For u at pole, there is one single mass point located at the pole correspond to one circle of v wind component (Figure 3a), which results in much more than twice as many velocity DOFs as mass and will therefore suffer from spurious computational modes (Staniforth & Thuburn, 2012). Nevertheless, for the triangular cells next to the poles with v at pole (Figure 3b), each triangle has 1 mass DOF and 1 1 2 velocity DOFs. Although this triangular C-grid does not have enough velocity DOFs, it is less severe than that for u at pole. From the test results below, the new calculation method of polar PV with u at pole seems to control the spurious modes effectively.

Analysis of Polar Noise and Assessment of the New Shallow-Water Model
The spatial discretization of the vector-invariant form of the equations on the lat-lon grid is described in the previous section. Another important part is the temporal discretization scheme in the evolution equations.
In the following numerical experiments, the explicit second-order predictor-corrector method (Wang & Ji, 2006) is applied. The full discretized equations can be written in the following form: where F= (u,v,h) T is the prognostic variables in the vector form and L denotes the discrete spatial operator; thus, LF denotes the explicit time tendency. To solve the above equation, the predictor-corrector method needs three iterative calculations, including one predictor substep and two corrector substeps The first two substeps in the predictor-corrector integrator update the state from time level n to an intermediate time level, and the last substep obtains the state at the next time level n+1. The detailed derivation of this time integration scheme is presented in Appendix B.
In addition, as in Williamson et al. (1992), numerical errors are estimated quantitatively by using two norm and infinity norm defined by where λ and φ are the longitude and latitude of the grid points, respectively, h is the model output, h t is the true solution if there is an analytic solution or a reference solution, and I is a discrete approximation to the global integral In addition, when an analytic solution is not available, the simulation calculated by the National Center for Atmospheric Research (NCAR) Spectral Transform Shallow Water Model Worley & Toonen, 1995) is used for comparison.

Zonal Flow Over an Isolated Mountain
Test case 5 of Williamson et al. (1992) describes a zonal flow impinging on an isolated mountain of conical shape. The surface or mountain height h s is given by h s =h s0 (1−r/R), where h s0 = 2,000 m, R=π/9, and 10.1029/2020MS002047

Journal of Advances in Modeling Earth Systems
The center of the mountain is located at λ c =3π/2, φ c =π/6. The wind velocity and height field are similar to the steady-state geostrophic flow test case in section 4.4, except α = 0, h 0 = 5,960 m, and u 0 = 20 ms −1 . The initial steady shear-free westerly flow is in geostrophic balance with the geopotential height. As the flow impinges on the mountain, an imbalance between the Coriolis force and the pressure gradient is induced, generating large-amplitude inertia-gravity waves and Rossby waves. After 15 days of simulation, these waves spread around the globe (including the poles). The interaction between the sole forcing orography and the zonal flow leads to strong nonlinearity (e.g., Ringler et al., 2011), which is particularly appropriate for assessing the effectiveness of the numerical method in conserving integral invariants, such as total mass, total energy, and total potential enstrophy (e.g., Nair et al., 2005).
First, the geopotential height and the zonal velocity component simulated by GAMIL-SW at day 20 are shown in Figure 4. Observing the animation of the u, the grid-scale noise propagates out from the North Pole. Moreover, the propagation is in the meridional direction, which could not be filtered by the FFT filter. There is, however, no noise in the geopotential height field during the same simulation. This evidence suggests that the wind field is more sensitive than the height field, which is probably due to their different time tendencies. If the nonlinear momentum advection terms are switched off, the noise would disappear in the wind field. This implies that the nonlinear momentum advection terms have an important role in the presence of noise. Second, as shown in Figure 5, the potential enstrophy increases abruptly when the waves arrive at the poles and the grid-scale noise appears, even if the global total energy is exactly conserved over the 20 simulation days. This demonstrates that energy conservation alone does not prevent the buildup of grid-scale oscillations, though does help to suppress nonlinear instability (Thuburn, 2008). The potential enstrophy is approximately conserved before the presence of the noise, which raises the question that whether the noise can be suppressed by conserving or dissipating the potential enstrophy. The equations of GAMIL-SW does not include PV explicitly, and the model is not designed to control the potential enstrophy either. Hence, for the two considerations, we decided to apply the vector-invariant form of shallow-water equations. This suite of equations does not explicitly contain the nonlinear momentum advection terms and can be flexibly designed to dissipate potential enstrophy, as described in section 3. Replacing the antisymmetric equations in the original GAMIL-SW with this suite of equations, a new shallow-water model is developed on the same lat-lon grid.
The simulation results of the new shallow-water model with two variable placements of u at pole and v at pole and different schemes for the same test case are shown in Figure 6. For v at pole, there is no noise at all in the simulation by using any scheme, as shown in Figures 6a and 6b,

Journal of Advances in Modeling Earth Systems
where the energy-conserving scheme and potential enstrophy-conserving scheme are applied, respectively. Meanwhile, Figure 7 shows the time series of the relative change in the total energy and total potential enstrophy with different schemes. For v at pole, the total energy is conservative within temporal truncation error. The total potential enstrophy shows slightly dissipation due to the nonlinear cascade from resolved to unresolved scales, and the rate of dissipation would decrease with increasing the spatial resolution (not shown). In contrast, for u at pole, whether implementing the energy-conserving and potential enstrophy-dissipating (APVM, Equation 14) scheme or the potential enstrophy-conserving scheme, the simulated u wind component (Figures 6c and 6d) shows noises propagating from the pole when the waves arrive at the pole, which is similar to that of GAMIL-SW. In other words, the conservation of energy or potential enstrophy does not suppress the presence of polar noise, and the dissipation of potential enstrophy does not remove the noise either. Moreover, the flow is weakly nonlinear during the first 15days, and the total potential enstrophy should be approximately conserved (Thuburn et al., 2014). However, in this study, the conservation properties of total energy and total potential enstrophy are disturbed once the noise begins to appear at about model day 13, which can be seen in Figure 7. This phenomena is similar with GAMIL-SW that shown in Figure 5. Thus, it indicates that the vector-invariant form of equations does not exterminate the source of polar noise; other more effective treatments must be sought. In addition, it is also found that the integration would blow up suddenly at about model day 40 with the potential enstrophy-conserving scheme, whereas the energy-conserving scheme does not. Moreover, it is also noted that the rate of energy dissipation in u-at-pole model is slightly larger than that in v-at-pole model, which is caused by the larger time step size used in u at pole (30s) than the one used in v at pole (15s). As shown in Figure 7, if the time step size in u-at-pole model is reduced to 15s, the rate of energy dissipation would be almost identical to that of v at pole. As a result, it can be concluded that the energy-conserving scheme conserves total energy within time truncation error.
The numerical schemes implemented in the simulations for u at pole and v at pole are identical except variable placement at the pole. The simulated difference reveals that the PV or relative vorticity must be well defined at the pole nonsingularly (e.g., Lin & Rood, 1997). Figures 6e and 6f show the results simulated with u at pole but the circle of relative vorticity next to the pole are calculated using Stokes' theorem ( Figure 3c) with the energy-conserving scheme and potential enstrophy-conserving scheme, respectively. In the calculation, the circle of v wind component nearest the poles is excluded from the calculation of polar PV. It can be clearly observed that the noises are eliminated by the modification (Figures 6e and 6f) and the conservation properties of energy and potential enstrophy are guaranteed much well (Figure 7). This minor modification to the calculation of PVs closest to the pole for u at pole results in an unexpected benefit of

Journal of Advances in Modeling Earth Systems
controlling the noise. Although the modification may affect the accuracy near and at the pole, the influence is relatively small, which will be verified in the following test case of steady-state geostrophic flow. Furthermore, the global conservation properties are maintained when the noises are suppressed. For example, as shown in Figure 7, the modified u-at-pole model with the same time step size (15 s) as the one used in v-at-pole model, the rates of energy and potential enstrophy dissipation are basically identical to that of v-at-pole model. Moreover, the rate of potential enstrophy dissipation in the modified u-at-pole model would be reduced when the spatial resolution is increased from 2°to 1.5°, which is similar to the situation in the v-at-pole model.
In order to further validate the effectiveness of this modification of the u at pole, two sensitivity experiments are conducted. The control experiment is the simulation with energy-conserving scheme which presents polar noise. The first sensitivity experiment runs the model with the polar PV modification but restarts from the control experiment at model day 15 when the polar noise is apparent in the wind field. Compared with the control experiment (Figure 8a), the sensitivity experiment shows that the modification suppresses the evolution of polar noise after 5days of integration as shown in Figure 8b. The second sensitivity experiment is to demonstrate that the circle of v wind component nearest the poles is connected with the noise. After every time step, the v wind component nearest the poles are restored by weighted averaging with the v wind component at the interior grid, for example, v(:,1)=0.2 * v(:,1)+0.8 * v(:,2) near the South Pole. The result shown in Figure 8c demonstrates that the noise at model day 20 is eliminated to a large extent. Therefore, the noise probably belongs to one of computational modes near the poles on the lat-lon grid. The mismatch in the number of DOFs in the wind and mass fields with u at pole may include extra branches of computational modes (Staniforth & Thuburn, 2012;. Further investigation will be needed to fully determine the mechanism of this kind of noise.

Cross-Polar Rotating High-Low
To further investigate whether the previous result is case dependent, the test case proposed by McDonald and Bates (1989) is evaluated, which simulates cross-polar flow with a geostrophically balanced initial state. This test has also been used by Giraldo et al. (2002), Nair et al. (2005) and Jablonowski et al. (2009). The initial condition consists of a height field and the wind field (u, v) derived from the height field via geostrophic relationship. They are given by where gh 0 =5.768 × 10 4 m 2 s −2 and v 0 =20 m s −1 . a and Ω are the radius and rotation rate of the Earth, 10.1029/2020MS002047

Journal of Advances in Modeling Earth Systems
respectively. It consists of a low and a high, which are symmetrically located on the west and east side of the pole in the Northern Hemisphere (their positions are reversed in the Southern Hemisphere). The low and high rotate in a clockwise direction around the North Pole and deform slightly. They exchange their positions after 5 days, and the slightly deformed pattern almost returns to the initial location after 10 days of integration. The maximum wind speed and a strong gradient of wind are near the poles; thereby, this test is well suited for the analysis of polar noise, although this test case has no analytical solution.
The numerical simulations are illustrated in Figure 9, which shows the u wind component at day 10 with different numerical schemes. It can be observed that the GAMIL-SW simulation shows grid-scale noise in the polar domain of about 30° (Figure 9a). Similarly, the new shallow-water model with u at pole also shows noise in the polar domain whether by using the energy-conserving and potential enstrophy-dissipating scheme (Figure 9b) or the potential enstrophy-conserving scheme (figure not shown). Moreover, the APVM dissipation of potential enstrophy is applied on the u at pole, but the result is nearly identical to that in Figure 9a without little improvement. However, the polar noises are eliminated by using Stokes' theorem is the result of u at pole with energy-conserving and potential enstrophy-dissipating scheme, (c) is the result of u at pole with energy-conserving scheme but relative vorticity next to the pole is calculated using Stokes' theorem, and (d) is the result of v at pole with energy-conserving scheme. Contour intervals are 4ms −1 . As in Figure 6, the blank circle shown at the pole in (d) is caused by the undefined u wind component at that grid point.

Journal of Advances in Modeling Earth Systems
for the ring of relative vorticity nearest the pole; Figure 9c shows the result with energy-conserving scheme.
With the modification, the simulation with u at pole is almost indistinguishable from the result of v at pole (Figure 9d). There are no distortions in the flow pattern with the latter two numerical schemes, and the strong gradient near the pole can be well simulated at the same time. The result of the v wind field shows the same situation and hence is not shown here. Therefore, this qualitatively comparative analysis verifies the conclusion from the previous test case.
On the basis of the analysis of polar noise in the two previous test cases, the following conclusions are obtained. There are two variable placements on the pole for both GAMIL-SW and the new shallow-water model: u at pole and v at pole. Both of the two shallow-water models with the u-at-pole configuration exhibit polar noise in the same test cases. The new shallow-water model with v at pole shows no noise but exhibits polar noise with u at pole, and the noise would be eliminated if the relative vorticity nearest the poles is calculated with one minor modification, which demonstrates that the numerical treatment of relative vorticity or PV on the pole is important for controlling polar noise (e.g., Arakawa & Lamb, 1981;Lin & Rood, 1997). The noise in the new shallow-water model with u at pole was not alleviated by conserving or dissipating potential enstrophy, which demonstrates that the noise is probably produced by some computational modes.
The new shallow-water model with v at pole exhibits no noise, while GAMIL-SW with v at pole exhibits noise (figure not shown), which implies that the noise is relevant to the form of equations, numerical schemes, pole problem, and so on. Overall, arranging the v velocity component located at the pole may be a good choice, although at the cost of a relatively small time step size.

Rossby-Haurwitz Wave
The previous two test cases analyzed the numerical polar noise. In the following, the performance of the new shallow-water model is shown in several test cases by comparison with GAMIL-SW or the spectral model. Test case 6 of Williamson et al. (1992) consists of a Rossby-Haurwitz wave of zonal wavenumber 4. This type of wave is an analytic solution for the fully nonlinear nondivergent barotropic vorticity equation on a sphere and has also been widely used to test shallow-water models. Nevertheless, the Rossby-Haurwitz wave is actually unstable as a solution of the shallow-water equations which was analyzed by Thuburn and Li (2000), because small random perturbations in the initial conditions could result in long-term disruption. This is shown to be the case for a wide range of numerical models, including the spectral model, which usually apply diffusion or damping to simulate more stably (e.g., Lauritzen et al., 2011;Whitehead et al., 2011). Moreover, if the model chooses a grid that is not symmetrical across the Equator, the disruption would come faster. However, the lat-lon grid naturally has the advantage of symmetry across the Equator.
A complete description of the test case is described in detail by Williamson et al. (1992). The initial height field is chosen to be in balance with the velocity field such that the initial velocity field is nondivergent. The minimum fluid height is 8,000 m, which occurs at the poles, and the mean fluid height is 9,523 m. In addition, the angular velocity of this pattern moving from west to east can be calculated by where R=4 for the wavenumber, ω=7.848 × 10 −6 s −1 and Ω ¼ 2π 86400 ¼ 7:272 × 10 −5 s −1 . With these parameters, the period is approximately 29.36 days for zonal wavenumber 4 of the Rossby-Haurwitz wave. Jakob et al. (1993) questioned how long the initial solution could be expected to remain stable. To address this question, they integrated a T42 model with the time step of 600 s for 60 days and concluded that a viable numerical method should be able to maintain the wavenumber 4 structure for a minimum of 14 days (Jakob-Chien et al., 1995). However, some studies selected different days for their assessment, such as 60, 14, or 10 days (e.g., Bonaventura & Ringler, 2005;Ii & Xiao, 2010;Li et al., 2008;Lin & Rood, 1997;Ringler et al., 2010). Figure 10 presents the height field after simulating three periods (88 days) using the spectral model, GAMIL-SW, and the new shallow-water model, respectively. The solution of the spectral transform method can be considered as the reference solution because the initial flow field can be exactly represented by the basic functions that are spherical harmonics (Lin & Rood, 1997). In fact, GAMIL-SW is able to maintain the basic pattern for as long as 100 days, and the new shallow-water model is also able to reach this level under certain conditions. Here, no attempt is made to obtain the longest simulation time because the spectral model is the best one without doubt. Compared to the reference solution, both the new shallow-water model and GAMIL-SW simulate steadily for three periods, and the results are remarkably similar to the reference solution, both in phase and amplitude, which demonstrates that the new shallow-water model maintains the computational stability of GAMIL-SW in this test case. However, when the same equation and numerical schemes are applied on a SCVT grid symmetrical with the Equator, the model is only able to accomplish the simulation without deformation for about 40 days ( Figure 10d). Therefore, in terms of this test case, the lat-lon grid has a natural advantage for simulating the zonally balance flow. This is also one of the main reasons why we still choose the lat-lon grid to design new dynamical core.

Steady-State Geostrophic Flow
Test case 2 of the standard shallow-water suite of Williamson et al. (1992) consists of a steady-state nonlinear zonal geostrophic flow. This test case measures the ability of the numerical scheme to maintain a large-scale geostrophic balance, which is an important property of any numerical model for the atmosphere or ocean.
Since the steady-state analytical solution is equal to the initial geostrophically balanced flow field, any deviation from the analytical solution is due to discretization errors, and the mesh convergence rate can be calculated by applying different spatial resolutions. The initial conditions of the fluid depth and the eastward and northward components of velocity at longitude λ and latitude φ are where gh 0 ¼ 2:94 × 10 4 m 2 s −2 and u 0 =2πa /86,400 × 12 m s −1 . The Coriolis parameter is The parameter α is the angle between the axis of the flow orientation and the polar axis of the Earth sphere. We analyze here only the results with α = 0 for the zonal geostrophic flow. The ℓ 2 and ℓ ∞ error norms calculated with the thickness field of the simulation after 5 days by GAMIL-SW and the new shallow-water model with energy-conserving scheme are shown in Figure 11. The spatial resolutions are 2°, 1°, 0.5°, and 0.25°, and the corresponding time steps for the simulation are 240, 120, 60, and 30 s, respectively. In this test case with α = 0, the APVM scheme is equal to the midpoint scheme because the second term in the APVM (Equation 14) will vanish because u e and ∇ e q are perpendicular (Ringler et al., 2010). It can be observed that for resolution higher than 1°, both ℓ 2 and ℓ ∞ error norms of the two models converge with second-order accuracy, and the convergence rate increases with increasing resolution. Comparing the error norms at the same resolution, the new shallow-water model generates smaller errors than GAMIL-SW. Moreover, with respect to the convergence rate, the vector-invariant equations on the lat-lon grid are better than the similar equations and numerical method implemented on other unstructured grid, such as SCVT (Ringler et al., 2010). This suggests that the new shallow-water model meets the aim of achieving second-order accuracy. In addition, the low height field error also shows that the new shallow-water model maintains large-scale balance and has steady geostrophic modes without grid-scale oscillations, which is another proof that the lat-lon grid with physically perfect zonal symmetry is in favor of being computationally free of grid-scale oscillations for this kind of zonal flow.

Barotropically Unstable Zonal Jet
The last test case to be discussed concerns the growth of rapid barotropic instability from a midlatitude jet, in which the initial wind is zonally symmetric (see Galewsk et al., 2004, for detailed numerical description).
Simulating the barotropically unstable jet with an initial small-amplitude perturbation is a challenging task for a quasi-uniform grid. As the jet is fast and narrow, numerical truncation errors arising from the misalignment of the jet with the grid would lead to perturbations that also produce instability in a similar manner to the initial perturbation (e.g., Ringler et al., 2011;Weller et al., 2012). Therefore, the shallow-water model with a quasi-uniform grid generally needs higher resolution than the previous test cases to reduce the numerical truncation errors (e.g., Ii & Xiao, 2010;Thuburn et al., 2014). For the lat-lon grid, the numerical truncation errors are relatively small because the jet is essentially aligned with the grid. However, since the velocity of gravity wave calculated by the geopotential ( ffiffiffiffiffi gh p ) is about 310 ms −1 in this test case, thus the explicit time integration step is restricted by the pole problem of the lat-lon grid. In addition, there is no polar filter or other treatments at present; only the simulations by low Figure 11. Convergence rate of steady-state geostrophic flow as measured by the ℓ 2 and ℓ ∞ norms based on the geopotential height. Norms of solid and dashed lines are computed by the new shallow-water model and GAMIL-SW with respect to the analytic solution at day 5, respectively. Lines of slope = −1 (−2) represent the firstand second-order convergence rates, respectively.

Journal of Advances in Modeling Earth Systems
spatial resolution 1°are shown. The vorticity field at different stages simulated by GAMIL-SW and the new shallow-water model with the same resolution are shown in Figure 12. Compared with fig. 4 in Galewsk et al. (2004) that simulated using the spectral model without diffusion, the envelope of the growing barotropic instability of the two models performs essentially identical to the reference. The dominant ridge-troughridge pattern with the same amplitude and phase is present in the two simulations.

Summary and Conclusion
In this study, the polar noise found in the barotropic GAMIL-SW is analyzed and a solution to the problem is provided. A global mixed finite-volume/finite-difference shallow-water model with the vector-invariant equations on the same lat-lon grid is developed and evaluated, which is designed to maintain the conservation of total mass and total energy like GAMIL-SW; in addition, the new model is able to conserve the PV and allows for the conservation or dissipation of potential enstrophy.
First, focusing on the polar numerical noise, several approaches are attempted and tested. In the zonal flow over a mountain test case, compared with the GAMIL-SW, the new model with u-at-pole configuration presents polar noise, which demonstrates that the vector-invariant form of equations is not able to suppress the polar noise. The new model with v at pole exhibits no polar noise at all, which implies that the calculation of PV or relative vorticity near the pole is crucial. Based on these facts, the noise with u at pole is prohibited completely if the ring of relative vorticity closest to the pole is calculated using Stokes' theorem. Moreover, the global conservation properties are maintained with this modification, which is beneficial to the model. By contrast, neither potential enstrophy-conserving nor enstrophy-dissipating scheme is able to suppress the noise. These analysis results are confirmed by another cross-polar flow test case. Thus, it is reasonable to conclude that the polar noise is related to the equation form and numerical treatment of the singular pole on the lat-lon grid and probably belongs to an unphysical mode due to the severe mismatch in the relative number of mass and velocity DOFs. Selecting the vector-invariant momentum equation and vat-pole configuration on the lat-lon grid is a better choice although at the cost of a smaller explicit time step for the finite-difference scheme.
Second, the new model performs as well or better than GAMIL-SW on another three test cases. Rossby-Haurwitz waves can be simulated steadily for three circulations, and the difference from the spectral model is indistinguishable, which is as well as GAMIL-SW and much better than the model with the same scheme applied on the quasi-uniform meshes. As expected, the same as GAMIL-SW, the new model converges to second-order accuracy in the geostrophic flow test case, and the error of the new model is less than that of GAMIL-SW in terms of absolute accuracy at the same resolution. The simulation results for a barotropically unstable zonal jet are nearly identical to the results of GAMIL-SW and the spectral model.
Finally, our future goal is to develop a high-resolution atmospheric nonhydrostatic model. The development and assessment of the shallow-water model in this study have helped us to identify an appropriate horizontal momentum equation. Moreover, from the view point of computational efficiency associated with the convergence of meridians toward the pole on the lat-lon grid, a practical method for this issue is being developed, which can enlarge the explicit time step size by several times without zonal FFT filter.
To ensure that the second and third terms can be canceled, the discrete kinetic energy is derived as where A i denotes the area of primal cell and A ie is the area weight of normal velocity (u e ) within one primal cell.

Appendix B: The Predictor-Corrector Time Integration Scheme
The implicit midpoint method for the evolution equation dF dt ¼ LF can be written as where F= (u,v,h) T is the prognostic variables in the vector form for the shallow-water equations and L denotes the spatial operator. Although this method is stable theoretically, it cannot be solved for F n+1 directly in general. Thus, the predictor-corrector algorithm is one kind of iterative solver with fixed iterative number. Usually, one time of predictor and two times of corrector can be chosen to end the iteration, such that the implicit algorithm becomes explicit one. First, calculate an initial guess value via Euler forward; next, improve the initial guess through two times of iteration process Þ the first corrector Þ the second corrector : 8 > > > > > > < > > > > > > : Arranging the LHS of Equation B1, it can be rewritten as Define F n þ 1 2 ¼ F n þ 1 þ F n 2 ; the equation above can be reformulated as It is also an implicit method equivalent to Equation B1, and perform the predictor-corrector method (Equation B2); similarly, one obtains Thus, Equation 25 is derived through some arrangement, and the prognostic variable F n+1 at next time step can be obtained after ending the second corrector with the variable F n at current time step.