A global finite-element shallow-water model supporting continuous and discontinuous elements

This paper presents a novel nodal finite-element method for either continuous and discontinuous elements, as applied to the 2-D shallow-water equations on the cubed sphere. The cornerstone of this method is the construction of a robust derivative operator that can be applied to compute discrete derivatives even over a discontinuous function space. A key advantage of the robust derivative is that it can be applied to partial differential equations in either a conservative or a non-conservative form. However, it is also shown that discontinuous penalization is required to recover the correct order of accuracy for discontinuous elements. Two versions with discontinuous elements are examined, using either the g1 and g2 flux correction function for distribution of boundary fluxes and penalty across nodal points. Scalar and vector hyperviscosity (HV) operators valid for both continuous and discontinuous elements are also derived for stabilization and removal of grid-scale noise. This method is validated using four standard shallow-water test cases, including geostrophically balanced flow, a mountain-induced Rossby wave train, the Rossby–Haurwitz wave and a barotropic instability. The results show that although the discontinuous basis requires a smaller time step size than that required for continuous elements, the method exhibits better stability and accuracy properties in the absence of hyperviscosity.


Introduction
Modeling of the 2-D shallow-water equations is an important step in understanding the behavior of a numerical discretization for atmospheric modeling.In particular, the dynamical character of the global shallow-water equations is governed by features common with atmospheric motions including nonlinearity, barotropic Rossby waves and inertiagravity waves, without the added complexity of a vertical dimension.
This paper introduces a novel discrete derivative operator that is applied to the shallow-water equations on a manifold using continuous and discontinuous finite elements.This work is motivated by the flux correction methods of Huynh (2007) and Vincent et al. (2011), is an alternative to formulations with discontinuous elements that discretize the conservative equations of motion with explicit momentum fluxes (Giraldo et al., 2002;Nair et al., 2005) and generalizes both spectral element and discontinuous Galerkin methods.This approach is also quadrature free, requiring no integral computation.This paper further introduces a general variational discretization of the scalar and vector Laplacian oper-Published by Copernicus Publications on behalf of the European Geosciences Union.

P. A. Ullrich: A global finite-element shallow-water model
ator which is valid for continuous or discontinuous elements and only requires one communication per application of the Laplacian.
Discontinuous elements are potentially more desirable than continuous elements for several reasons: first, discontinuous elements only require parallel communication along coordinate axes, whereas continuous elements also require parallel communication along diagonals, a doubling of the total number of communications in 2-D.Second, discontinuous elements provide a natural mechanism to enforce stabilization via discontinuous penalization (or Riemann solvers, for equations in conservation form).Third, discontinuous elements can be used in conjunction with upwind methods, which are generally better for tracer transport and associated problems.However, discontinuous elements also have a number of disadvantages, including higher storage requirements (for the same order of accuracy), a maximum time step size which is typically smaller than that imposed for continuous elements (Ullrich, 2013) and added computational expense for many hyperbolic operations.
The outline of this paper is as follows.Section 2 presents the shallow-water equations on a manifold.The cubedsphere grid, which will be used for simulations on the sphere, is described in Sect.3. The discretizations of the dynamics and hyperviscosity (HV) are then presented in Sects.4 and 5, respectively.Results from four standard shallow-water test cases are given in Sect.6 and conclusions follow in Sect.7.

The shallow-water equations on a manifold
The 2-D shallow-water equations in on a Riemannian manifold with coordinate indices x s = {α, β} can be written as ∂u α ∂t + u s ∇ s u α + g αs ∂ ∂x s (g c H ) + f (k × u) α = 0, (1) ∂H ∂t + ∇ s (hu s ) = 0. (3) The prognostic variables are free-surface height H and vector velocity u = u α g α + u β g β , where g α = ∂x/∂α and g β = ∂x/∂β are the natural basis vectors on the manifold.The fluid height h and height of the bottom topography z are related to the free-surface height via H = h + z.Here g rs denotes the contravariant metric with covariant inverse g rs , J = √ detg rs is the metric Jacobian, g c is gravity, f is the Coriolis parameter and k is the vertical basis vector of unit length.The Einstein summation notation (implied summation) is used for repeated indices.These equations further make use of the covariant derivative ∇ s , which expands as where d sr denote the Christoffel symbols of the second kind associated with the coordinate transform (again with summation over repeated indices s and r implied).
The mass equation, Eq. ( 3), has been kept in conservative form to enforce strict mass conservation.On the other hand, Eqs. ( 1)-( 2) are given in a non-conservative form; this formulation is selected over the flux-form equations (where hu α and hu β are prognostic variables).Angular momentum and potential enstrophy are particularly relevant to atmospheric motions (Thuburn, 2008) and can be easily conserved under a non-conservative formulation of the shallow-water equations (Taylor and Fournier, 2010).Conservation of these quantities is more difficult when they are diagnosed from the flux-form prognostic variables.The non-conservative formulation also has the advantage of leading to a more accurate treatment of wave-like motion when formulated appropriately (Thuburn and Woollings, 2005).

The cubed-sphere grid
The Eqs. ( 1)-( 3) are now applied to a particular choice of coordinate system.The cubed-sphere grid (Sadourny, 1972;Ronchi et al., 1996) consists of a cube with six Cartesian patches arranged along each face, which is then inflated onto a tangent spherical shell, as shown in Fig. 1.The cubed sphere is a quasi-uniform spherical grid, that is, it is in the class of grids that provide an approximately uniform tiling of the sphere (see Staniforth and Thuburn (2012), for example, for a review of different options for global grids).On the equiangular cubed-sphere grid, coordinates are given as (α, β, p), with central angles α, β ∈ [− π 4 , π 4 ] and panel index p ∈ {1, 2, 3, 4, 5, 6}.By convention, we choose panels 1-4 to be along the equator and panels 5 and 6 to be centered on the northern and southern pole, respectively.With uniform grid spacing, each panel consists of a square array of n e × n e elements.
The contravariant 2-D metric on the equiangular cubed sphere of radius a is given by where δ = 1 + tan 2 α + tan 2 β.The Jacobian on the manifold, denoted by J , is then and induces the infinitesimal area element dA = J dα dβ.
The Christoffel symbols of the second kind are given by Spherical coordinates (λ, φ) for longitude λ ∈ [0, 2π ] and latitude φ ∈ [−π/2, π/2] will also be used for plotting and specification of tests.Coordinate transforms between spherical and equiangular coordinates can be found in Ullrich and Jablonowski (2012) Appendix A.

The nodal basis
A nodal finite-element method is employed (Taylor et al., 1997;Giraldo et al., 2002;Hesthaven and Warburton, 2007).The 1-D reference element is defined as the interval x ∈ [−1, 1] along with a set of test functions φ(i) (x).The test functions are defined such that test function φ(i) (x) is the unique polynomial of degree n p that is one at the ith Gauss-Lobatto-Legendre (GLL) node ,i ∈ (0, . .., n p − 1), and zero at all other GLL nodes.Each basis polynomial then has a corresponding weight, defined by where α = α 2 − α 1 and β = β 2 − β 1 .The accompanying 2-D tensor-product basis is then defined by Figure 2 provides a depiction of GLL nodes within a single element.For vector quantities (such as velocity u), test functions are instead vector fields.Uniqueness of the variational system is retained if exactly 2 degrees of freedom are allowed at each nodal location for the vector test function φ.As we shall see, the most natural choice is test functions φ

Robust differentiation
A robust differentiation operator is now constructed for both continuous and discontinuous finite elements.Let f : (α, β) → R be defined and continuous on Z ∪ ∂Z with basis φ (i,j ) : for coefficients f (p,q) ∈ R. Further, let f : (α, β) → R be defined and continuous on ∂Z.Here f represents the evaluation of f in neighboring elements.Note that for a continuous finite-element method, f and f must satisfy f (α, β) = f (α, β) on ∂Z, whereas no such restriction is imposed for discontinuous finite elements.Following Huynh (2007), robust differentiation in the α direction is defined at GLL nodes via where the overscore denotes the co-located average of f and f , An analogous definition holds in the β direction.Here g L and g R are the local flux correction functions, which are chosen to satisfy and otherwise are chosen to approximate zero throughout Several options for g L and g R will lead to a stable discretization, including g 1 (Radau polynomials), which will lead to the discontinuous Galerkin method, and g 2 , which will lead to the mass-lumped discontinuous Galerkin method (Huynh, 2007).Hereafter discontinuous elements with the g 1 flux correction function will be referred to as discontinuous g 1 elements, whereas elements using of the g 2 flux correction function will be referred to as discontinuous g 2 elements.Observe that for continuous finite elements, the rightmost two terms in Eq. ( 15) are exactly zero.
With the definition of a robust discrete derivative Eq. ( 15), discretization of the shallow-water system, Eqs.(1)-(3), is straightforward.Note that for continuous finite elements, this discretization is identical to the approach of Taylor et al. (1997) when applied in conjunction with direct stiffness summation (DSS, that is, projection into the space of continuous functions; see Appendix A).If the conservative form of the shallow-water equations was employed, this discretization would be the same as Giraldo et al. (2002) when mass lumping is not employed (discontinuous g 1 ) and Nair et al. (2005) if mass lumping is applied (discontinuous g 2 ).To the best of the author's knowledge, no previous work has used both discontinuous elements and a non-conservative form of the shallow-water system.

Discontinuous penalization
At element boundaries, the use of one-sided derivatives will cause the discontinuity between neighboring elements to ex-hibit an error with magnitude O( x n p −1 ), an effective loss of 1 order of accuracy from the expected convergence rate.To reduce errors associated with the discontinuity, a penalization term is added in each coordinate direction.In the α direction this term reads as where λ(α, β) = |u α | + √ gh/a represents the maximum local wave speed in the α direction.An analogous term is added in the β direction.Note that with this choice of penalization, the evolution equation for H is identical to the evolution equation that would arise from a traditional conservative discontinuous Galerkin method with local Lax-Friedrichs flux.Since the penalization term is equivalent to upwinding, it is weakly diffusive and so allows the discontinuous scheme to maintain stability even in the absence of explicit viscosity.

Implementation considerations
On the cubed-sphere grid, the discontinuous method has 6 n 2 e n 2 p degrees of freedom compared to 8 + 8(n e (n p − 1) − 1)+6(n e (n p −1)−1) 2 for the continuous method.In the limit as n e → ∞, this yields a ratio of (n p −1) 2 /n 2 p degrees of freedom for the continuous formulation versus the discontinuous formulation.Note that in practice, the continuous formulation typically stores redundant degrees of freedom in order to reduce computational expense associated with indexing and so memory requirements are typically identical.
The primary computational difference between the continuous and discontinuous formulations is due to the evaluation of the penalty terms, Eqs. ( 18)-( 19).Note that although the robust differentiation operation, Eq. ( 15), does require additional computation for discontinuous methods, the cost of evaluating the discontinuous terms in this expression is roughly equivalent to the computational cost of the direct stiffness summation operation needed for continuous elements.Nonetheless, from numerical experiments the discontinuous method has an approximately 30 % overhead compared with a continuous method (when run with the same time step size).

Viscosity and hyperviscosity
Stabilization is typically needed for co-located (or unstaggered) finite-element methods, whether implicitly in the form of upwinding or explicitly in the form of a diffusive operator, to avoid high-frequency dispersive errors associated with spectral ringing.In general, it is preferred that this operator is consistent with the underlying geometry of the Riemannian manifold, which precludes, for example, the Boyd-Vandeven filter (Boyd, 1996).There has been considerable success with the use of hyperviscosity in the spectral element method (Dennis et al., 2011), which maintains geometric consistency by mimicking the natural fourth-order hyperviscosity operator.Previously, it has not been clear how to extend this operator to a discontinuous function space.However, the robust derivative, Eq. ( 15), provides a direct mechanism by which the hyperviscosity operator can be constructed.The viscosity operator for both the continuous and discontinuous function space will be discussed here.
Note that any viscosity operator will lead to a loss of energy conservation of the underlying numerical method.This loss is exhibited in two obvious ways: first, for geostrophically balanced flows the error will tend to grow over time.Second, energy conservation is lost leading to a decay in the total energy content of the system over time.

Scalar viscosity
For stabilization of the method, diffusion is added in the form of either viscosity or hyperviscosity, which corresponds to multiple applications of the viscosity operator.A scalar viscosity operator is constructed to satisfy where ∇ 2 = ∇ • ∇ is the usual scalar Laplacian.The operator is defined implicitly via a variational construction.If f = H(ν)ψ then, multiplying Eq. ( 20) by a test function and applying integration by parts, one obtains where dS is the infinitesimal line element along the boundary of Z and dA is the infinitesimal area element.The two terms on the right-hand side of this expression correspond to the viscosity flux through element boundaries and the Laplacian within the element.Under a continuous element formulation, only the rightmost term is retained and fluxes are instead accounted for via direct stiffness summation.Under a discontinuous formulation, both terms are retained and discretized.The discrete equation satisfied by f (i,j ) that follows from Eq. ( 21) is written as where f B (i,j ) denotes the discretization of the boundary integral and f A (i,j ) denotes the discretization of the area integral.After a lengthy derivation (see Appendix B), these discretizations read as and where δ i,j is the Krönecker delta.Here the contravariant derivative of ψ has been used, defined via Note that the contravariant derivatives ∇ p ψ are multivalued along this edge, and so must be adjusted using the robust derivative operator Eq. ( 15).

Vector viscosity
Vector viscosity is used for damping of the velocity field, and takes the form Observe that if ν = ν d = ν v then this expression is exactly the standard vector Laplacian operator ∇ 2 u, with coefficient ν.By writing the vector Laplacian as Eq. ( 26), the combined operator separates into two distinct operators that affect divergence damping (with coefficient ν d ) and vorticity damping (with coefficient ν v ).This result can be quickly verified by taking the divergence and curl of Eq. ( 26), For simplicity of calculation, we treat divergence damping and vorticity damping separately.For divergence damping, the operator is constructed by taking the inner product of f = H(ν d , ν v )u with the vector test function φ, integrating over Z and applying integration by parts, www.geosci-model-dev.net/7/3017/2014/Geosci.Model Dev., 7, 3017-3035, 2014 For vorticity damping an analogous procedure leads to Note that for shallow-water flows, only the radial component of the vorticity is relevant for this calculation.The discrete value of f α (i,j ) and f β (i,j ) that arises from this calculation then has contributions from the area integral via f A,d  (i,j ) and boundary integral via f B,d  (i,j ) , Following another lengthy derivation (see Appendix B), the area integral term appears as whereas the boundary integral term is Applying an analogous procedure for test function φ The divergence and curl, which are needed for evaluation of the Laplacian, are computed via where All partial derivatives are evaluated using the robust derivative operator (15).

Hyperviscosity
For stabilization of a high-order discretization, hyperviscosity is preferred since it retains the order of accuracy of the underlying scheme.In practice, hyperviscosity is implemented by repeated application of the viscosity operator.For instance, for fourth-order hyperviscosity, the ∇ 4 operator is approximated as follows

Computational considerations
Calculation of hyperviscosity in the form presented here requires one parallel exchange per application of the Laplacian operator.For continuous elements, this communication is manifested through the application of DSS, which averages away any discontinuity that has been generated along element edges.For discontinuous elements, scalar viscosity requires pointwise updates along element edges computed from Eq. ( 24), whereas vector viscosity requires both onesided values of u, (∇ •u) and (∇ ×u) r , which are in turn used for computing nodal values of (∇ • u) and (∇ × u) r needed for Eqs. ( 32)-( 35).This constitutes a doubling of the overall bandwidth requirement relative to continuous elements.

Results
In this section selected results are provided to evaluate the relative performance of the methods described in this paper.Four test cases are evaluated: from the Williamson et al. (1992) test case suite, steady-state geostrophically balanced flow, zonal flow over an isolated mountain and the Rossby-Haurwitz wave will be analyzed, in addition to the barotropic instability test of Galewsky et al. (2004).For all test cases, time integration is handled via the strong-stability preserving three-stage, third-order Runge-Kutta method (Gottlieb et al., 2001).Note that some improvement in the maximum time step size for discontinuous elements can be obtained from the use of the five-stage, third-order Runge-Kutta method (Ruuth, 2006), which has a stability region that covers a larger portion of the negative real plane.The time step t for each test is chosen to be close to the stability limit in each case (observed empirically).The value of t has a negligible effect on the results (not shown), suggesting that spatial errors dominate in each case.Further, note that mass conservation is maintained to machine truncation for all simulations (not shown).From the shallow-water equations, the values of g c and f for the Earth are used, All simulations are performed with n p = 4.A thorough investigation of different values of n p would greatly extend the length of the manuscript, so n p was chosen in accordance with the Community Atmosphere Model spectral element dynamical core.As argued by Ullrich (2013), this choice is also optimal when considering the accurate treatment of linear waves.
When required, the standard L 2 error measure is calculated via where h T is the height field at the initial time (which is the analytical solution for steady-state test cases) and I denotes an approximation to the global integral given by where the subscript k denotes the values of x and J within the kth element.When applied, hyperviscosity uses a single coefficient for both scalar and vector hyperviscosity, This choice of scaling for the hyperviscosity coefficient is based on Takahashi et al. (2006).

Steady-state geostrophically balanced flow
Test case 2 of Williamson et al. (1992) describes a zonally symmetric geostrophically balanced flow.This test utilizes an unstable equilibrium solution to the shallow-water equations which is expected to be exactly maintained over time.However, it is generally true that only methods that satisfy the curl-grad annihilator property ∇ ×∇φ = 0 maintain some sort of discrete equilibrium.Nonetheless, since an analytical solution is known (identical to the initial conditions), this test is effective at measuring the convergence rate of a numerical method.Further, the error fields from this test provide some indication of what effect the grid has on the errors of the underlying method.The analytical height field for this test is given by with background height h 0 and velocity amplitude u 0 chosen to be This height field also serves as the reference solution.The non-divergent velocity field is specified in latitude-longitude (φ, λ) coordinates as Figure 3 shows L 2 errors in the height field after a 5-day integration of the model at n e = 4 resolution with n p = 4. Simulations were completed for continuous elements (a) with hyperviscosity and (d) without hyperviscosity, discontinuous elements (b, e) with mass lumping (the g 2 flux correction function), (c, f) without mass lumping (the g 1 flux correction function), (b, c) with discontinuous penalization and (e, f) without discontinuous penalization.The time step is t = 2200 s for simulations (a, d), t = 800 s for simulations (b, c, e) and t = 400 s for simulation (f).Increasing the magnitude of the time step by 100 s led to simulation instability in each case.Since the addition of hyperviscosity leads to loss of energy conservation, there is a slow decay of the geostrophically balanced flow towards a uniform height state, hence leading to a nearly zonally symmetric decay in the height field towards the poles.For all configurations (both continuous and discontinuous elements) visually identical results are observed when hyperviscosity is added, and so these results are not shown.All simulations exhibit a characteristic wave number 4 mode triggered by the underlying cubed sphere, although the specific error pattern differs throughout.Simulation   but show stable error norms, discontinuous elements with penalization show smaller error norms than continuous elements but a very slow growth with time due to the upwinding effect of the discontinuous penalization and discontinuous elements without penalization show rapid growth in errors (and even instability without mass lumping).
To verify that the model exhibits the correct convergence rate, Fig. 5 shows the global error norms associated with simulations with n e ∈ {4, 8, 16, 32, 64} over a 5-day integration period.At n e = 4, the time step is t = 2200 s for continuous elements, t = 800 s for g 2 discontinuous elements and g 1 discontinuous elements with penalization, and t = 400 s for g 1 discontinuous elements without penalization.Increasing the time step by 100 s led to an unstable simulation.The time step is scaled inversely with increasing resolution.Missing simulations correspond to model instability and divergence prior to simulation completion.The use of hyperviscosity reduces the convergence rate to O( x 3.2 ), as expected from the choice of hyperviscosity coefficient, Eq. ( 44).With hyperviscosity disabled, model simulations converge at O( x 4 ) for continuous elements and discontinuous elements with penalty, and O( x 3 ) for discontinuous elements without penalty.The loss of 1 order of accuracy is due to one-sided differentiation at co-located nodes along element edges, leading to enhancement of the discontinuity.Similar results (not shown) are observed when changing n pthat is, continuous elements and discontinuous elements with penalty converge at O( x n p ), whereas unpenalized discontinuous elements converge at O( x n p −1 ).

Zonal flow over an isolated mountain
Test case 5 in Williamson et al. (1992) considers zonal flow with underlying topography.The wind and height fields are defined as in section 6.1, except with h 0 = 5960 m and u 0 = 20 m s −1 , respectively.A conical mountain is used for the topographic forcing, given by with The center of the mountain is at λ c = 3π/2 and φ c = π/6.
Simulation results for this test case were computed at n e = 16 and n p = 4 after 15 days of integration both with and without hyperviscosity.For discontinuous elements penalization is always used.The time step used for these runs was t = 480 s for continuous elements, t = 240 s for g 2 discontinuous elements and t = 120 s for g 1 discontinuous elements.Increasing the time step by 20 s led to an unstable simulation.These results are visually indistinguishable, so are instead compared against the continuous element run (with HV) in Fig. 6, where height field h and height field difference h − h c are plotted (where h c is the height field given in a).Simulations (b) and (c), corresponding to discontinuous elements with and without mass lumping, are very similar in structure and exhibit smooth differences from the continuous model.With no hyperviscosity applied, continuous elements (d) show significant noise which is not present for discontinuous elements (e, f).These simulations match closely with results from the literature (Nair et al., 2005;Ullrich et al., 2010).
To understand conservation of invariants over time, total energy E and potential enstrophy ξ are computed over the duration of the simulation.Since these quantities are invariant under the shallow-water equations, it would be expected that a perfect simulation would conserve these quantities exactly.They are defined via A time series of energy and potential enstrophy are plotted in Fig. 7.With hyperviscosity (a, b) all simulations exhibit nearly identical conservation properties, suggesting that both the continuous and discontinuous hyperviscosity operators (which are responsible for the loss of energy and potential enstrophy conservation) act in a nearly identical manner over the course of the simulation.Without hyperviscosity (c, d) change in energy and potential enstrophy is much smaller.Continuous elements show initiation of instability at approximately day 6, likely due to high-wave-number oscillations in the height field caused by nonlinear aliasing.On the other hand, discontinuous elements instead show a slow decay of energy and potential enstrophy driven by the weak diffusivity of the discontinuous penalization.

Rossby-Haurwitz wave
Test case 6 in Williamson et al. (1992) consists of a westward-propagating Rossby-Haurwitz wave that exactly solves the barotropic vorticity equation, but only approximately solves the nonlinear shallow-water equations.This test is particularly interesting since it is known to be sensitive to the choice of horizontal viscosity.
Results for the Rossby-Haurwitz wave are given in Figs. 8  and 9 for n e = 16 and n p = 4 horizontal resolution, respectively, after 14 days of integration.The time step used for these runs was t = 480 s for continuous elements, t = 200 s for g 2 discontinuous elements and t = 120 s for g 1 discontinuous elements.Increasing the time step by 20 s led to an unstable simulation.As expected, there are significant differences in the height field which are induced by the addition of the hyperviscosity (although both simulations appear reasonable given the coarse horizontal resolution).Except for this difference, the results are nonetheless analogous to zonal flow over an isolated mountain: stable without the addition of hyperviscosity, whereas discontinuous elements with penalization are effective at stabilizing the method for both lumped and non-lumped variants.

Barotropic instability
The barotropic instability test case of Galewsky et al. (2004) consists of a zonal jet with compact support at a latitude of 45 • , with a latitudinal profile roughly analogous to a much stronger version of test case 3 of Williamson et al. (1992).A small height perturbation is added atop the jet which leads to the controlled formation of an instability in the flow.The relative vorticity of the flow field at day 6 can then be visually compared against a high-resolution numerically computed solution (Galewsky et al., 2004;St-Cyr et al., 2008).
Simulation results for this test case were computed at n e = 32 and n p = 4 after 12 days of integration with hyperviscosity enabled.The time step used for these runs was t = 150 s for continuous elements, t = 75 s for g 2 discontinuous elements and t = 50 s for g 1 discontinuous elements.Increasing the time step by 10 s led to an unstable simulation.Simulations are again compared against the continuous element run (with HV) in Fig. 10, where the relative vorticity field ζ and relative vorticity field difference ζ − ζ c is plotted (where ζ c is the height field given in a).Due to the presence of sharp frontal activity in this test case and the strong resolution dependence of this problem (Ullrich et al., 2010), differences in ζ are of the same magnitude as the original field.In particular, the simulations without hyperviscosity (d, e, f) all show enhancement near wave fronts which is not apparent in the simulations with hyperviscosity (b, c).Although most differences can be found near sharp fronts, there is also a clear enhancement in the differences near 120 E associated with a trailing instability.For continuous elements without hyperviscosity (c), there is also apparent grid-scale noise which is missing from the other simulations, suggesting that this method is under-diffused.
Normalized total energy and potential enstrophy are plotted for the barotropic instability in Fig. 11 for a 12-day integration with n e = 16 and n p = 4.With hyperviscosity (a, b) there are small but visible differences in the results associated with changes in the type of elements.Without hyperviscosity (c, d) the simulation with continuous elements exhibit instability around day 6, leading to rapid growth of energy and potential enstrophy.On the other hand, with discontinuous elements there is a steady loss of energy and potential enstrophy over time due to diffusivity from discontinuous penalization.Prior to wave breaking (which occurs around day www.geosci-model-dev.net/7/3017/2014/Geosci.Model Dev., 7, 3017-3035, 2014 4), energy and potential enstrophy loss are significantly reduced compared to the simulations with hyperviscosity.After wave breaking, energy and potential enstrophy loss are of the same order of magnitude for simulations with and without hyperviscosity, associated with the fact that diffusivity is enhanced near the barotropic fronts where discontinuities are large.

Conclusions
Following Huynh (2007), a novel nodal finite-element method for continuous and discontinuous elements has been constructed using a robust derivative operator and discontinuous penalization.The resulting methodology can be used for straightforward discretization of partial differential equations in either a conservative or a non-conservative form.A hyperviscosity operator valid for both continuous and discontinuous elements was also presented that would provide a mechanism for numerical stabilization and the removal of grid-scale noise.Two versions with discontinuous elements were studied, using either the g 1 or g 2 flux correction function for distribution of boundary fluxes and penalty across nodal points.The resulting method was then applied to the 2-D shallow-water equations in cubed-sphere geometry and tested on a suite of test problems.
From the Williamson et al. (1992) test case suite, steadystate geostrophically balanced flow, zonal flow over an isolated mountain and the Rossby-Haurwitz wave were examined, in addition to the barotropic instability test of Galewsky et al. (2004).The method was shown to be stable and accurate for both continuous and discontinuous elements, with fourth-order convergence being verified for cubic basis functions.Discontinuous penalization was shown to be necessary for stability and for maintaining the correct order of accuracy of the discontinuous method.Overall the discontinuous elements required a smaller time step than for continuous elements, although all methods led to similar error norms when hyperviscosity was active.When hyperviscosity was deactivated, the discontinuous method exhibited better error norms than the continuous approach, and discontinuous penalization was shown to be sufficient for stability of the method even for complex flows.Nonetheless, differences between all three approaches appeared minor, and all methods performed well for this suite of tests.
The non-conservative discontinuous element formulation has been shown to be a potential candidate for future atmospheric modeling.It has the advantage of requiring fewer parallel communications than continuous methods, and exhibits stability even when hyperviscosity is not used for explicit stabilization.However, with the reduced time step size it remains unclear whether the discontinuous formulation would be computationally competitive with continuous element methods.
The method discussed in this paper has been implemented in the Tempest atmospheric model, available from https: //github.com/paullric/tempestmodel.where we have used g ββ = J 2 g αα .Repeating for all edges and using Eq.(B2) then yields Eq. ( 24).

B2 Vector viscosity
The area integral that appears on the left-hand side of Eqs. ( 29) and ( 30) takes the form

↵Fig. 2 .Figure 2 .
Fig. 2. A depiction of the nodal grid for a reference element on GLL nodes for n p = 4. Boundary nod which are connected to neighboring elements, are shaded.

Figure 3 .
Figure 3. L 2 errors in Williamson et al. (1992) test case 2, steady-state geostrophically balanced flow for n e = 4 and n p = 4 after a 5-day integration period.Contour spacing for plot (a) is 1 m.Contour spacing for all other plots is 0.5 m.The zero line is enhanced.Long dashed lines show the cubed-sphere grid.

PFigure 4 .
Fig.4.L 2 error time series for geostrophically balanced flow on the cubed-sphere for n e = 16 and n p = 4 over a 5 day integration period for all continuous and discontinuous schemes.
Figure3shows L 2 errors in the height field after a 5-day integration of the model at n e = 4 resolution with n p = 4. Simulations were completed for continuous elements (a) with hyperviscosity and (d) without hyperviscosity, discontinuous elements (b, e) with mass lumping (the g 2 flux correction function), (c, f) without mass lumping (the g 1 flux correction function), (b, c) with discontinuous penalization and (e, f) without discontinuous penalization.The time step is t = 2200 s for simulations (a, d), t = 800 s for simulations (b, c, e) and t = 400 s for simulation (f).Increasing the magnitude of the time step by 100 s led to simulation instability in each case.Since the addition of hyperviscosity leads to loss of energy conservation, there is a slow decay of the geostrophically balanced flow towards a uniform height state, hence leading to a nearly zonally symmetric decay in the height field towards the poles.For all configurations (both continuous and discontinuous elements) visually identical results are observed when hyperviscosity is added, and so these results are not shown.All simulations exhibit a characteristic wave number 4 mode triggered by the underlying cubed sphere, although the specific error pattern differs throughout.Simulation (d) is exactly mimetic and leads to exact maintenance of geostrophic balance.Simulations (b)

Figure 6 .Figure 7 .
Figure 6.Height field with n e = 16 and n p = 4 at day 15 for zonal flow over an isolated mountain with (a) continuous elements and hyperviscosity (reference solution).Height difference plot from reference solution with n e = 16 at day 15 for (b) discontinuous g 2 elements with hyperviscosity, (c) discontinuous g 1 elements with hyperviscosity, (d) continuous elements without hyperviscosity, (e) discontinuous g 2 elements without hyperviscosity and (f) discontinuous g 1 elements without hyperviscosity.The time step used for these runs was (a, d) t = 480 s, (b, e) t = 240 s and (c, f) t = 120 s.Discontinuous penalization was used for both discontinuous schemes.Contour spacing is 1 m in all difference plots with the zero line removed.Long dashed lines show the cubed-sphere grid.

Figure 8 .Fig. 9 .Figure 9 .
Figure 8. Height field with n e = 16 and n p = 4 at day 14 for the Rossby-Haurwitz wave with (a) continuous elements and hyperviscosity (reference solution).Height difference plot from reference solution with n e = 16 at day 14 for (b) discontinuous g 2 elements with hyperviscosity, (c) discontinuous g 1 elements with hyperviscosity, (d) continuous elements without hyperviscosity, (e) discontinuous g 2 elements without hyperviscosity and (f) discontinuous g 1 elements without hyperviscosity.The time step used for these runs was (a, d) t = 480 s, (b, e) t = 200 s and (c, f) t = 120 s.Discontinuous penalization was used for both discontinuous schemes.Contour spacing is 1 m in plots (b) and (c) and 20 m in plots (d), (e) and (f).Long dashed lines show the cubed-sphere grid.

PFigure 10 .
Fig.10.Relative vorticity field with n e = 32 and n p = 4 at day 6 for the barotropic instability test with (a) continuous elements and hyperviscosity (reference solution).Relative vorticity difference plot from reference solution with n e = 16 at day 6 for (b) discontinuous g 2 elements with hyperviscosity, (c) discontinuous g 1 elements with hyperviscosity, (d) continuous elements without hyperviscosity, (e) discontinuous g 2 elements without hyperviscosity and (f) discontinuous g 1 elements without hyperviscosity.The time step used for these runs was (a,d) ∆t = 150 s, (b,e) ∆t = 75 s and (c,f) ∆t = 50 s.Discontinuous penalization was used for both discontinuous schemes.Contour spacing in all plots is 2 × 10 −5 s −1 with the zero line removed.Long dashed lines show the cubed-sphere grid.

Figure 11 .
Fig. 11.Normalized total energy and enstrophy change for the barotropic instability test with n e = 16 and n p = 4 over a 12 day simulation.In (c) and (d) the continuous element simulation fails after approximately 6 days, leading to unbounded growth in energy and enstrophy.The time step used for these runs was (a,d) ∆t = 300 s, (b,e) ∆t = 150 s and (c,f) ∆t = 75 s.Discontinuous penalization was used for both discontinuous schemes.46