the Creative Commons Attribution 4.0 License.
the Creative Commons Attribution 4.0 License.
A mixed finiteelement discretisation of the shallowwater equations
James Kent
Thomas Melvin
Golo Albert Wimmer
This paper introduces a mixed finiteelement shallowwater model on the sphere. The mixed finiteelement approach is used as it has been shown to be both accurate and highly scalable for parallel architecture. Key features of the model are an iterated semiimplicit timestepping scheme, a finitevolume transport scheme, and the cubed sphere grid. The model is tested on a number of standard spherical shallowwater test cases. Results show that the model produces similar results to other shallowwater models in the literature.
The works published in this journal are distributed under the Creative Commons Attribution 4.0 License. This license does not affect the Crown copyright work, which is reusable under the Open Government Licence (OGL). The Creative Commons Attribution 4.0 License and the OGL are interoperable and do not conflict with, reduce or limit each other.
© Crown copyright 2022
The dynamical core of an atmospheric model numerically approximates the solution to the governing fluid dynamics equations that determine the evolution of the atmosphere. An operational dynamical core must be both accurate, to give confidence in the forecast, and efficient, to produce a forecast in the desired time frame. Current supercomputer architectures focus on using an increasing number of processors to decrease runtime, and so future dynamical cores must be suitable for massively parallel architecture. Key aspects of the dynamical core that will greatly affect both accuracy and parallel scalability are the type of grid and the numerical methods used to discretise the equations.
Traditionally, atmospheric models used a latitude–longitude grid (Wood et al., 2014). However, the convergence of the meridians at the pole leads to a computational bottleneck, and thus the latitude–longitude grid is not suitable for future supercomputers. Hence, a number of different approaches to gridding the sphere without a pole have been developed and used in dynamical cores, such as icosahedral grids and the yin–yang grid (Staniforth and Thuburn, 2012). One such promising grid is the cubed sphere, which uses quadrilateral cells on six panels that cover the sphere (Ronchi et al., 1996).
The choice of numerical methods used in a dynamical core will affect both the accuracy and the efficiency of the model. Highorder methods generally improve accuracy but increase the computational cost. There are a variety of spatial methods (for example, finitedifference, finitevolume, semiLagrangian) and temporal methods (for example, explicit, semiimplicit, horizontally explicit, vertically implicit) used in atmospheric models by different operational centres and modelling groups (see Ullrich et al., 2017, and references within). The mixed finiteelement method of Cotter and Shipton (2012) and Cotter and Thuburn (2014) is a spatial method that allows highorder schemes to be defined whilst retaining good parallel properties (Melvin et al., 2019). Mixed finite elements retain many properties of Cgrid finitedifference and finitevolume methods, such as good numerical dispersion relations. They also have the advantage of not requiring orthogonal grids (Cotter and Shipton, 2012), better consistency of the Coriolis operator (Thuburn and Cotter, 2015), and the flexibility to increase accuracy by using higherorder elements. In practice, most current atmospheric dynamical cores aim for secondorder accuracy overall (with the transport scheme often higher order) (Wood et al., 2014). For this reason the lowestorder mixed finiteelement method will be a building block of the model described in this article.
The shallowwater equations are an important step in testing the methods used for the development of a dynamical core (Williamson et al., 1992). The shallowwater equations contain many of the properties of the full atmospheric governing equations but with a reduction in complexity. In practice, if a method is not good enough (by whatever metric) for a shallowwater model, it will not be good enough for an atmospheric dynamical core. This allows the methods to be evaluated in a simpler setting and design decisions to be made early in development.
This document presents a novel mixed finiteelement discretisation, extending that of Melvin et al. (2019) for atmospheric dynamics in Cartesian geometry, of the rotating shallowwater equations in a spherical domain. This is the natural next step from Melvin et al. (2019) in the development of a threedimensional dynamical core in spherical geometry, allowing the impact of the chosen spherical grid to be investigated in a simplified context using an evolution of the same numerics. As with the model of Melvin et al. (2019), the shallowwater discretisation presented here has been developed in the LFRic framework: see Adams et al. (2019). Key components of the model are the semiimplicit iterated timestepping scheme and the finitevolume transport. The governing shallowwater equations are given in Sect. 2. The discretisation, including the spatial and temporal aspects, is given in Sect. 3, the finitevolume transport scheme is described in Sect. 4, and the solution procedure is outlined in Sect. 5. Results from standard shallowwater test suites are given in Sect. 7, with a concluding summary in Sect. 8.
The shallowwater equations in a rotating domain are given in vectorinvariant form by
where $\mathit{u}=(u,v)$ is the velocity, Φ=gh is the geopotential (with gravity g and free surface height h), Φ_{s} is the surface geopotential, and $K\equiv \mathrm{1}/\mathrm{2}\mathit{u}{}^{\mathrm{2}}$ is the kinetic energy. The perpendicular operator is defined by $(u,v{)}^{\u27c2}=(v,u)$.
The potential vorticity (PV) q is defined as
where f is the Coriolis parameter and ${\mathrm{\nabla}}^{\u27c2}\cdot \equiv \mathit{k}\cdot \mathrm{\nabla}\times $, and taking the curl of Eq. (1) gives a conservative transport equation for PV:
The governing Eqs. (1) and (2) are discretised in time using a twotimelevel iterated implicit scheme (Sect. 3.1) and in space using a mixed finiteelement method (Sect. 3.2) for the wavedynamic terms and a highorder upwind finitevolume scheme (Sect. 4) for the advection terms. This process is an extension to the shallowwater equations on the sphere of Melvin et al. (2019), who presented a similar discretisation in a threedimensional Cartesian domain.
3.1 Temporal discretisation
To achieve secondorder temporal accuracy, a timecentred approach is used. The target discretisation, as in Melvin et al. (2019), is a twotimelevel iterated implicit scheme where terms responsible for the fastwave dynamics are treated by the iterative semiimplicit scheme and the transport terms (those involving the mass flux in Eq. 2 and the potential vorticity in Eq. 1) are computed using a highorder, upwind, explicit finitevolume scheme. Taking Eqs. (1) and (2) and applying this discretisation results in
where ${\mathit{\delta}}_{t}s\equiv \left({s}^{n+\mathrm{1}}{s}^{n}\right)/\mathrm{\Delta}t$, ${\stackrel{\mathrm{\u203e}}{s}}^{\mathit{\alpha}}\equiv \mathit{\alpha}{s}^{n+\mathrm{1}}+\left(\mathrm{1}\mathit{\alpha}\right){s}^{n}$, and 𝓕(s,u) are the fluxes computed by the transport scheme of variable s by wind field u. We use $\mathit{\alpha}=\mathrm{1}/\mathrm{2}$ to achieve the second order centred in the time scheme.
3.2 Mixed finiteelement discretisation
The mixed finiteelement formulation in two spatial dimensions requires the specification of three finiteelement function spaces: ${\mathbb{V}}_{\mathrm{0}},\phantom{\rule{0.125em}{0ex}}{\mathbb{V}}_{\mathrm{1}},\phantom{\rule{0.125em}{0ex}}{\mathbb{V}}_{\mathrm{2}}$ (cf. the four function spaces ${\mathbb{W}}_{i},\phantom{\rule{0.125em}{0ex}}i=\mathrm{0}\mathrm{\dots}\mathrm{3}$ used in Melvin et al., 2019). The scalar spaces are an H_{1} space consisting of pointwise scalars, 𝕍_{0} (zero forms), or an L_{2} space consisting of areaintegrated scalars, 𝕍_{2} (two forms). There are two choices for the vector space 𝕍_{1} corresponding to either ${\mathbb{V}}_{\mathrm{1}}^{C}$, a H_{curl} space of circulation vectors, or ${\mathbb{V}}_{\mathrm{1}}^{D}$, a H_{div} space of flux vectors. Each choice has an associated discrete de Rham complex,
for curlconforming vectors and
for divconforming vectors. For this paper only the divconforming complex Eq. (8) will be considered: ${\mathbb{V}}_{\mathrm{1}}\equiv {\mathbb{V}}_{\mathrm{1}}^{D}$, with u∈𝕍_{1} and Φ∈𝕍_{2}. This is analogous to the standard Cgrid staggering (Arakawa and Lamb, 1977) where the normal components of the velocity vector are stored at the cell edges. The divconforming complex, with the locations of the prognostic variables u, v, and Φ, is shown in Fig. 1.
The potential vorticity q∈𝕍_{2} is collocated with Φ and is computed diagnostically as follows. The velocity u is represented in the H_{curl} space as $\mathit{v}\in {\mathbb{V}}_{\mathrm{1}}^{\mathrm{C}}$ using a Galerkin projection, and then the curl is taken to give the relative vorticity $\mathit{\omega}\equiv \mathit{k}\cdot \mathrm{\nabla}\times \mathit{v}\in {\mathbb{V}}_{\mathrm{2}}$. The absolute vorticity is projected into 𝕍_{2} and then divided by the geopotential. As the lowestorder space is used, this gives the potential vorticity q∈𝕍_{2}.
Taking Eq. (5), multiplying by a test function w from the velocity space and integrating over the domain D gives
where, since the geopotential and kinetic energy are discontinuous between cells, the third term has been integrated by parts and the boundary term vanishes due to the continuity of the test functions w. Similarly, taking Eq. (6), multiplying by a test function σ from the geopotential space and integrating over the domain gives
where it is assumed that the advection scheme returns fluxes 𝓕∈𝕍_{1}.
3.3 Transforms
As in Melvin et al. (2019), the equations are transformed from a physical cell C to a reference cell $\widehat{\mathrm{C}}$ using the mapping $\mathit{\varphi}:\widehat{\mathrm{C}}\to \mathrm{C}$. The physical cell, on the cubed sphere, has coordinates χ, and the reference cell, a unit square, has coordinates $\widehat{\mathit{\chi}}$, and the transform is such that $\mathit{\chi}=\mathit{\varphi}\left(\widehat{\mathit{\chi}}\right)$. Transforming the equations to a single reference cell provides a number of computational efficiencies such as a single set of basis functions and quadrature points (Rognes et al., 2009). The Jacobian of this transformation is defined as $\mathbf{J}\equiv \partial \mathit{\varphi}\left(\widehat{\mathit{\chi}}\right)/\partial \widehat{\mathit{\chi}}$ and is used in transforming variables between the physical and reference cells. The transformations used for spaces 𝕍_{1} and 𝕍_{2} are designed to preserve fluxes through an edge (𝕍_{1}) and areaintegrated values (𝕍_{2}) respectively. The transformation for 𝕍_{1} is $\mathit{v}\left(\mathit{\chi}\right)\equiv \mathit{v}\left(\mathit{\varphi}\left[\widehat{\mathit{\chi}}\right]\right)=\mathbf{J}\widehat{\mathit{v}}\left(\widehat{\mathit{\chi}}\right)/\text{det}\mathbf{J}$. Following Melvin et al. (2019), for the 𝕍_{2} transformation rehabilitation, the approach of Bochev and Ridzal (2010) is used so that the 𝕍_{2} mapping is modified to $\mathit{\sigma}\left(\mathit{\chi}\right)\equiv \mathit{\sigma}\left(\mathit{\varphi}\left[\widehat{\mathit{\chi}}\right]\right)=\widehat{\mathit{\sigma}}\left(\widehat{\mathit{\chi}}\right)$. Applying these to Eqs. (9) and (10) gives
and
The transport scheme is an extension to the methodoflines scheme used by Melvin et al. (2019), computing fluxes 𝓕 of a scalar field s by a wind field u:
The flux 𝓕 is obtained using a methodoflines scheme where a conservative transport equation
is solved to obtain 𝓕. The temporal aspects of this scheme are handled in the same manner as Melvin et al. (2019) using an mstage Runge–Kutta scheme
The coefficients a_{i,j} and b_{k} in Eqs. (15) and (16) are given by the Butcher tableau for the scheme:
0  

c_{2}  a_{2,1}  
c_{3}  a_{3,1}  a_{3,2}  
⋮  ⋮  ⋱  
c_{m}  a_{m,1}  ⋮  ${a}_{m,m\mathrm{1}}$  
b_{1}  b_{2}  ⋮  b_{m−1}  b_{m} 
Here the threestage thirdorder strong stability preserving (SSP3) method of Gottleib (2005) is used, which has the Butcher tableau:
0  

1  1  
$\mathrm{1}/\mathrm{2}$  $\mathrm{1}/\mathrm{4}$  $\mathrm{1}/\mathrm{4}$  
$\mathrm{1}/\mathrm{6}$  $\mathrm{1}/\mathrm{6}$  $\mathrm{4}/\mathrm{6}$ 
At each stage i we need to compute $\mathit{F}\left({s}^{\left(i\right)},\mathit{u}\right)\equiv {\stackrel{\mathrm{\u02c7}}{s}}^{\left(i\right)}\mathit{u}$, where $\stackrel{\mathrm{\u02c7}}{s}$ is a highorder upwind reconstruction of s. We use a quadratic reconstruction for $\stackrel{\mathrm{\u02c7}}{s}$ (see Sect. 4.1), resulting in a transport scheme that is third order in both time and space. As the model uses the cubed sphere grid, the spatial reconstruction of Melvin et al. (2019) is extended to take into account the nonuniformity of the mesh by using a twodimensional horizontal reconstruction. The scheme defined in this section works on order l=0 spaces.
4.1 Reconstruction of a scalar field
The advection scheme computes a highorder upwind reconstruction $\stackrel{\mathrm{\u02c7}}{s}$ of a given scalar field s. The reconstructed field is computed at points staggered half a grid length in all directions from the original field, so for a field s∈𝕍_{2}, which is located at cell centres, the reconstructed field $\stackrel{\mathrm{\u02c7}}{s}$ is computed at the centre of each cell edge. The reconstruction is computed by fitting a polynomial through a number of cells and evaluating this polynomial at the staggered points. The reconstruction is such that the integral of the polynomial is equal to the integral of s within each cell. This is given an upwind bias by choosing evenorder polynomials for the reconstruction which require an odd number of s points and hence can be weighted to the upwind side of the point at which the reconstruction is needed; for example, for a onedimensional quadratic reconstruction at a point ${\stackrel{\mathrm{\u02c7}}{s}}_{i+\mathrm{1}/\mathrm{2}}$ with a positive wind field, two upwind points s_{i−1} and s_{i} and one downwind point s_{i+1} are used.
The horizontal spatial reconstruction is based on that used in Thuburn et al. (2014), which uses a similar method to Baldauf (2008) and Skamarock and Menchaca (2010). To summarise, a series of polynomials P_{k} of a given order n in a local Cartesian coordinate system (x,y) is defined over a stencil of n_{s} cells. The polynomial is required to fit (in a least squares sense) the discrete field being reconstructed. The integral along the cell edge of the reconstructed field $\stackrel{\mathrm{\u02c7}}{s}$ is approximated by Gaussian quadrature and is given by
where (x_{j},y_{j}) are the integration points and w_{j} the weights of the Gaussian quadrature. In practice, a n_{q}=2 point quadrature is used, and this is found to give a small improvement over singlepoint quadrature.
The weights P_{k}(x_{r},y_{r}) that multiply each value s_{k} in the stencil are obtained by evaluating a polynomial
at (x_{r},y_{r}). The coefficients ${a}_{i,j}^{k}$ of P_{k} are determined by minimising the residual r_{k}
so that the integral of P_{k}=1 in cell k and P_{k}=0 otherwise.
For an order n reconstruction, there are ${n}_{m}\equiv (n+\mathrm{1})(n+\mathrm{2})/\mathrm{2}$ coefficients ${a}_{i,j}^{k}$, and so, to avoid an underdetermined problem, this requires at least n_{m} cells in the stencil. Additionally, the stencil should be symmetric about the central cell. To ensure these properties hold, the stencils are generated in the same manner as Thuburn et al. (2014). To summarise, the following algorithm is used.

Add the central cell to the stencil (if n=0, stop).

Loop until the number of cells in the stencil n_{s} is at least the number of monomials n_{m}.

Find the set S of all neighbouring cells of cells currently in the stencil.

Either add all cells in S that are not already in the stencil and are neighbours of two cells already in the stencil or, if no cell in S is a neighbour of two cells in the stencil, add all cells in S that are not already in the stencil.
An example of the type of stencil this generates around a corner of the cubed sphere is shown in Fig. 2.
For example, with a quadratic reconstruction n=2 there are n_{m}=6 monomials, and the stencil algorithm will generate a stencil with n_{s}=9 cells in general and n_{s}=8 cells near the corners of the cubed sphere. As in Thuburn et al. (2014) the central cell (k=1) in the stencil is fitted exactly (r_{1}=0 in Eq. 19), and the others are fitted in a least squares sense.
The local Cartesian coordinates (x,y) are computed as in Thuburn et al. (2014), the origin of the coordinate system x_{0} is taken to be the centre of the cell in the centre of the stencil, and the direction of the x axis is then taken as the direction from x_{0} to an arbitrary neighbour. Then it is straightforward to reconstruct any point ${\mathit{x}}_{i}\equiv \left({X}_{i},{Y}_{i},{Z}_{i}\right)$ in the stencil in terms of the local coordinates (x_{i},y_{i}); see Thuburn et al. (2014) for more details.
4.2 Flux computation
Once the scalar field s has been reconstructed at the locations of the 𝕍_{1} degrees of freedom to give $\stackrel{\mathrm{\u02c7}}{s}$, the flux of that field on a cell face is obtained by multiplying the reconstructed scalar by the normal component of the velocity field at that point:
where $\stackrel{\mathrm{\u02c7}}{\mathbf{s}}$ is the vector containing all $\stackrel{\mathrm{\u02c7}}{s}$ at cell edges.
4.3 Predictors
The temporal discretisation is designed to mimic a semiimplicit semiLagrangian discretisation such as that used in Wood et al. (2014). For the transport scheme this includes advecting predictors (Φ^{p},q^{p}) for the geopotential and potential vorticity fields instead of the values at the start of the time step (Φ^{n},q^{n}). This is motivated by considering a semiLagrangian discretisation of a prototypical equation
with $\mathrm{D}/\mathrm{D}t\equiv \partial /\partial t+\mathit{u}\cdot \mathrm{\nabla}$. Discretising across a trajectory from x_{D} at time level n to x_{A} at time level n+1 gives
with subscripts “A” and “D” denoting evaluation at arrival x_{A} and departure x_{D} points respectively. Evaluation of a function F at a departure point can be expressed as
where 𝒜_{sl} is the operation of the semiLagrangian advection operator. Applying this approximation to Eq. (21) gives
where the subscript “A” has been dropped for convenience. Applying this idea to Eq. (2), the geopotential predictor to be advected is then
This can also be applied to potential vorticity using Eq. (4). These predictors come from considering the continuity and PV equations in advective form $\mathrm{D}s/\mathrm{D}t+s\mathrm{\nabla}\cdot \mathit{u}=\mathrm{0}$, where s=Φ or qΦ, as would be used in a semiLagrangian model.
The procedure for the solution of the shallowwater equations is as follows. The semiimplicit governing equations, Eqs. (5) and (6), are repeated here for clarity:
An iterated implicit scheme is used to solve Eqs. (27) and (28). At each stage (k) of the iterative scheme, the timelevel n+1 terms are lagged and included in the residuals such that
where k is the estimate for the n+1 terms after k iterations of the iterative scheme. Increments ${\mathit{u}}^{\prime}\equiv {\mathit{u}}^{\left(k+\mathrm{1}\right)}{\mathit{u}}^{\left(k\right)}$ and ${\mathrm{\Phi}}^{\prime}\equiv {\mathrm{\Phi}}^{\left(k+\mathrm{1}\right)}{\mathrm{\Phi}}^{\left(k\right)}$ to Eqs. (29)–(30) are sought such that the fastwave terms are handled implicitly:
where Φ^{*} is a reference state used to obtain the linearisation and τ is a relaxation parameter (usually chosen to be $\mathit{\tau}=\mathrm{1}/\mathrm{2}$). In practice, we use ${\mathrm{\Phi}}^{*}={\mathrm{\Phi}}^{n}$ as the reference state. Applying the mixed finiteelement discretisation, this becomes the system
where 𝓡_{u} and ℛ_{Φ} are the finiteelement discretisations of Eqs. (29) and (30) respectively given by Eqs. (11) and (12). The matrices are defined as
At each iteration (k) the system (33) is solved using an iterative Krylov subspace method (in this case the generalised minimal residual method – GMRES). For the tests presented in this paper (see Sect. 7), GMRES generally converged in two to three iterations to a tolerance of 10^{−4}.
An equiangular cubed sphere mesh is used to grid the sphere. We use the notation Cn to describe a cubed sphere with n×n cells per panel. The mesh is parameterised using a finiteelement representation of the sphere within a cell with polynomials of order m. Note therefore that a point within a cell does not necessarily lie on the sphere, with the error depending on the order of the elements used. For example, a linear element lies on the sphere at the vertices of the element and uses a linear approximation to the surface of the sphere within the element. Representing the sphere with these elements removes the need for analytical transforms to the reference grid, allowing the use of an arbitrary grid if required (in this paper we use the standard equiangular cubed sphere). To compute the error of the mesh parameterisation, the geocentric coordinates $(X,Y,Z)$ at a point in the cell are computed using the finiteelement representation. The error is the difference between the true radius of the sphere and the radius using $({X}^{\mathrm{2}}+{Y}^{\mathrm{2}}+{Z}^{\mathrm{2}}{)}^{\mathrm{1}/\mathrm{2}}$.
Figure 3 shows the error within a cell when using linear and quadratic elements on a C96 grid with the Earth's radius. The error for the linear element is largest at the cell centre, with a maximum error of 426.39 m. The quadratic element reduces the error by 5 orders of magnitude, with a maximum error of 0.0018 m. Figure 3c shows the convergence of the maximum error within a cell. The linear element error converges at second order, with the quadratic element converging at fourth order. For this reason the quadratic element is used to create the mesh for the shallowwater model. Note that the quadratic element is for the representation of the mesh, whereas the lowestorder finite elements are used for the spatial discretisation of the model.
This section shows the results of the model runs using a standard set of spherical shallowwater test cases. The full initial conditions for each of the tests are given in the references under each test case.
7.1 Williamson 2: steady state
The first test case is the steadystate test described in Williamson et al. (1992). The steady flow means that the initial conditions are the analytical solution at any time, and thus error norms can be calculated for the runs at different resolutions.
The normalised ℓ_{2} error norms of the geopotential field at 15 d are given in the top two rows of Table 1 for the C24 (Δt=3600), C48 (Δt=1800), and C96 (Δt=900) resolutions. The error norms can be used to determine the convergence rate and hence the empirical order of accuracy of the model. The convergence rates of the ℓ_{2} and ℓ_{∞} error norms are approximately second order. This is expected as the finiteelement discretisation and the timestepping scheme are both secondorder methods. The error fields after 15 d for the C96 resolution are shown in Fig. 4: these errors show a wave number 4 pattern coming from the underlying cubed sphere mesh; however, the errors are large scale and are not particularly clustered around the edges and corners of the cubed sphere, indicating an acceptable level of grid imprinting.
7.2 Williamson 5: mountain test
The mountain test case of Williamson et al. (1992) is used to show the performance of the model when orography is present (i.e. Φ_{s}≠0). The initial zonal flow is over a mountain centred at 90^{∘} longitude and 30^{∘} latitude. The same grid resolutions and time steps used for the Williamson 2 test are used for the mountain test.
The geopotential at day 15, shown in Fig. 5a for the C96 grid, is comparable to other shallowwater models at similar resolutions (see for example Thuburn et al., 2010, 2014; Ullrich, 2014). We investigate the errors in the total geopotential by comparing them with a highresolution solution generated using the semiimplicit semiLagrangian (SISL) shallowwater model of Zerroukat et al. (2009). This reference solution uses a time step of Δt=60 s and a highresolution latitude–longitude grid of 1536×768 points. The highresolution solution is then interpolated, using bicubic interpolation, to the corresponding cubed sphere grid where the difference is taken to produce the error plots and diagnostics.
The error plots for the total geopotential are shown in Fig. 6 for the C24 and C48 resolutions. The error for the C48 solution is visually similar to that of the finiteelement cubed sphere model of Thuburn and Cotter (2015) (their Fig. 7), with a large negative error near the mountain and the wave train error in the Southern Hemisphere. The errors decrease significantly as the resolution increases. The normalised ℓ_{2} and absolute ℓ_{∞} error norms, calculated from the highresolution SISL model, are given in the bottom of Table 1.
The C48 maximum errors (given in Fig. 6 and in Table 1) are smaller in magnitude than those for the 13 824cell finitevolume cubed sphere model of Thuburn et al. (2014) and have a similar magnitude to the SISL model (Zerroukat et al., 2009) at 160×80 resolution (note that Thuburn et al., 2014, use total height, not total geopotential, in their results, and so their results need to be multiplied by gravity to give an equivalent value). The absolute ℓ_{∞} error norm values converge at a rate of 1.6, similarly to those given in Thuburn et al. (2014) (note that they use a smaller time step than presented in our results).
The normalised ℓ_{2} errors of the geopotential, given in the bottom of Table 1, have a convergence rate of 2.42 between C24 and C48. However, at higher resolution the error convergence stalls, similarly to that found in Thuburn et al. (2014) for comparable time steps, with a rate of only 0.53 between C48 and C96. Even with this stalling, the total rate between C24 and C96 is 1.48.
To demonstrate the conservation properties of the model, the mass, total energy and potential enstrophy are computed at each time step. The total energy is defined as
and the potential enstrophy is defined as
where the integrals are over the whole domain.
Mass, energy, and potential enstrophy are conserved in the continuous equations. As the model uses a finitevolume transport scheme, the mass is conserved to machine precision; however, energy and potential enstrophy are not conserved by the discretisation presented here, and in fact, in the discrete case, potential enstrophy cascades downscale from resolved to unresolved scales and thus is expected to decrease with time.
The normalised total energy and potential enstrophy are plotted against time in Fig. 7 for different grid resolutions (C24, C48, and C96) up to day 50. As the resolution increases, the dissipation in the model decreases, and the total energy and potential enstrophy curves are closer to conservation (a horizontal line). The flow is initially weakly nonlinear, and so for the first 15 d there are no significant cascades to unresolved scales. After 15 d the percentage loss in total energy is 0.0355 % for C24, 0.0062 % for C48, and 0.001 % for C96. For potential enstrophy the percentage loss is 0.3648 % for C24, 0.076 % for C48, and 0.014 % for C96. Extending to 50 d, the flow becomes more nonlinear, and this results in more dissipation of energy and potential enstrophy. The percentage losses are 0.221 % for C24, 0.063 % for C48, and 0.014 % for C96 in the total energy and 3.33 % for C24, 2.19 % for C48, and 1.45 % for C96 in the potential enstrophy.
The potential vorticity at day 50 is shown in Fig. 5b for the C96 resolution. As the flow has become nonlinear, the potential vorticity contours start to wrap up. The model captures the potential vorticity filaments without producing noise. The results for this test indicate that the model is correctly simulating flow over orography.
7.3 Galewsky instability test
For the Galewsky instability test (Galewsky et al., 2004), a perturbation is added to a balanced jet to create a barotropic instability. As the instability progresses, many smallscale potential vorticity filaments are produced.
The potential vorticity solutions at day 6 are shown in Fig. 8 for the C48, C96, C192, and C384 resolution runs (with Δt=900, 450, 225, and 112.5 s respectively). At C48 resolution the grid imprinting from the cubed sphere is evident. A wave number 4 pattern appears on the jet, dominating the development of the instability. Increasing the resolution to C96, C192, or C384 reduces the impact of the grid imprinting; however, the instability has still developed more than in the reference solution of Galewsky et al. (2004) at T341 resolution. This is consistent with the finiteelement cubed sphere model of Thuburn and Cotter (2015), which uses a similar mixed finiteelement discretisation to the one used here. In Thuburn and Cotter (2015) it is stated that the grid imprinting may be exaggerated by the highly unstable initial state of this test.
At the higher resolution many smallscale features are resolved by the model without producing gridscale noise. This is due to the implicit diffusion from the transport scheme damping gridscale features. These solutions are comparable to the results shown in Thuburn et al. (2014) and Thuburn and Cotter (2015).
7.4 Vortices test
The final test case is the field of vortices test from Kent et al. (2016). A field of vortices (shown by the potential vorticity in Fig. 9a) is set to evolve over a number of days. The vortices interact, leading to smallscale vorticity filaments and fingers as the vortices are stretched out as well as the merger of vortices. A highresolution run (using C192 resolution, with potential vorticity shown after 15 d in Fig. 9b) is used as a reference solution.
Figure 10 shows the day 15 potential vorticity for the C24 and C96 runs (with Δt=1800 and 450 s respectively). The C24 run is unable to resolve many of the features of the vortices, but the implicit diffusion from the transport scheme is sufficient to prevent gridscale noise. For the higherresolution C96 run, the smallscale features of the potential vorticity are well represented, and the solution matches the reference solution well. This demonstrates the model's ability to represent smallscale features without producing gridscale noise.
This article presents a new shallowwater model on the sphere. The model is comprised of a mixed finiteelement spatial discretisation, a highorder finitevolume transport scheme, and an iterated semiimplicit time scheme and makes use of the cubed sphere grid. The finiteelement discretisation is chosen as it has been shown to be both accurate and scalable on many processors (Cotter and Shipton, 2012). The lowestorder finiteelement method is used to give secondorder spatial accuracy, and the finiteelement spaces are such that they are analogous to a Cgrid staggering. The semiimplicit time stepping provides stability for the fast gravity waves along with secondorder temporal accuracy. The methodoflines advection scheme, using thirdorder Runge–Kutta time stepping with a quadratic finitevolume reconstruction, gives high accuracy and conservation for transport.
The model achieves computational efficiency through the use of the mixed finiteelement method and the cubed sphere grid. The cubed sphere cells are more uniform in size than the traditional latitude–longitude grid, and so fewer cells are required for the same resolution. For example, a 1^{∘} × 1^{∘}resolution grid requires $\mathrm{360}\times \mathrm{180}=\mathrm{64}\phantom{\rule{0.125em}{0ex}}\mathrm{800}$ cells on a latitude–longitude grid but only 48 600 cells on a C90 cubed sphere grid (with both grids having 360 cells around the Equator). The cubed sphere grid also removes the pole problem of the latitude–longitude grid and thus is more scalable on massively parallel architecture (Staniforth and Thuburn, 2012).
The model has been tested on a standard set of shallowwater test cases. The Williamson steadystate test shows that the model converges at secondorder accuracy. The mountain test demonstrates the model's ability to capture flow over orography. Grid imprinting from the Galewsky instability test is evident at coarse spatial resolutions but is significantly reduced at higher resolutions, in line with results from other models with similar discretisations and cubed sphere meshes. The vortex field test highlights that the model can resolve smallscale features without producing gridscale noise. The results presented from this testing are comparable to other models in the literature. This indicates that this model has a similar level of accuracy to these other wellknown shallowwater models.
This shallowwater model has been developed alongside the mixed finiteelement Cartesian model for atmospheric dynamics of Melvin et al. (2019). The vectorinvariant model presented here is a building block towards a mixed finiteelement spherical geometry dynamical core for the atmosphere.
The shallowwater model code, at LFRic revision r39707, can be found on Zenodo with DOI https://doi.org/10.5281/zenodo.7446738 (Kent et al., 2022). This repository includes the configuration files used for the tests in this article and instructions for running the model. In addition, the code at the same revision is also available on GitHub at https://github.com/thomasmelvin/gunghoswe (last access: 15 February 2023).
The data can be accessed by running the shallowwater model code using the configurations specified in lfric/miniapps/shallow_water/configurations at https://doi.org/10.5281/zenodo.7446738 (Kent et al., 2022).
JK and TM developed the shallowwater code with contributions from GW. JK performed the model simulations. JK and TM prepared the manuscript with contributions from GW.
The contact author has declared that none of the authors has any competing interests.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
We would like to thank the two reviewers for comments that improved the manuscript. We would also like to thank John Thuburn for providing the semiimplicit semiLagrangian comparison model.
This paper was edited by James Kelly and reviewed by two anonymous referees.
Adams, S., Ford, R., Hambley, M., Hobson, J., Kavčič, I., Maynard, C., Melvin, T., Müller, E., Mullerworth, S., Porter, A., Rezny, M., Shipway, B., and Wong, R.: LFRic: Meeting the challenges of scalability and performance portability in Weather and Climate models, J. Parall. Distr. Com., 132, 383–396, https://doi.org/10.1016/j.jpdc.2019.02.007, 2019. a
Arakawa, A. and Lamb, V. R.: Computational design of the basic dynamical processes of the UCLA general circulation model, Methods Comp. Phys., 17, 174–265, 1977. a
Baldauf, M.: Stability analysis for linear discretisations of the advection equation with RungeKutta time integration, J. Comp. Phys., 227, 6638–6659, 2008. a
Bochev, P. B. and Ridzal, D.: Rehabilitation of the lowestorder RaviartThomas element on quadrilateral grids, SIAM J. Num. Anal., 47, 487–507, 2010. a
Cotter, C. J. and Shipton, J.: Mixed finite elements for numerical weather prediction, J. Comput. Phys., 231, 7076–7091, 2012. a, b, c
Cotter, C. J. and Thuburn, J.: A finite element exterior calculus framework for the rotating shallowwater equations, J. Comp. Phys., 257, 1506–1526, 2014. a
Galewsky, J., Scott, R. K., and Polvani, L. M.: An initialvalue problem for testing numerical models of the global shallowwater equations, Tellus A, 56, 429–440, 2004. a, b
Gottleib, S.: On high order strong stability preserving RungeKutta and multi step time discretizations, J. Sci. Comput., 25, 105–128, 2005. a
Kent, J., Jablonowski, C., Thuburn, J., and Wood, N.: An energyconserving restoration scheme for the shallowwater equations, Q. J. Roy. Meteor. Soc., 142, 1100–1110, https://doi.org/10.1002/qj.2713, 2016. a
Kent, J., Melvin, T., and Wimmer, G.: thomasmelvin/gunghoswe: LFric SWE (v1.1), Zenodo [code], https://doi.org/10.5281/zenodo.7446738, 2022. a, b
Melvin, T., Benacchio, T., Shipway, B., Wood, N., Thuburn, J., and Cotter, C.: A mixed finiteelement, finitevolume, semiimplicit discretization for atmospheric dynamics: Cartesian geometry, Q. J. Roy. Meteor. Soc., 145, 2835–2853, https://doi.org/10.1002/qj.3501, 2019. a, b, c, d, e, f, g, h, i, j, k, l, m
Rognes, M. E., Kirby, R. C., and Logg, A.: Efficient assembly of H(div) and H(curl) conforming finite elements, SIAM Journal on Scientific Computing, 31, 4130–4151, 2009. a
Ronchi, C., Iacono, R., and Paolucci, P.: The cubed sphere: A new method for the solution of partial differential equations in spherical geometry, J. Comput. Phys., 124, 93–114, 1996. a
Skamarock, W. and Menchaca, M.: Conservative transport schemes for spherical geodesic grids: Highorder reconstructions for forwardintime schemes, Mon. Weather Rev., 138, 4497–4508, 2010. a
Staniforth, A. and Thuburn, J.: Horizontal grids for global weather prediction and climate models: a review, Q. J. Roy. Meteor. Soc., 138, 1–26, 2012. a, b
Thuburn, J. and Cotter, C. J.: A PrimalDual Mimetic Finite Element Scheme for the Rotating Shallow Water Equations on Polygonal Spherical Meshes, J. Comput. Phys., 290, 274–297, https://doi.org/10.1016/j.jcp.2015.02.045, 2015. a, b, c, d, e
Thuburn, J., Zerroukat, M., Wood, N., and Staniforth, A.: Coupling a massconserving semiLagrangian scheme (SLICE) to a semiimplicit discretization of the shallowwater equations: Minimizing the dependence on a reference atmosphere, Q. J. Roy. Meteor. Soc., 136, 146–154, https://doi.org/10.1002/qj.517, 2010. a
Thuburn, J., Cotter, C. J., and Dubos, T.: A mimetic, semiimplicit, forwardintime, finite volume shallow water model: comparison of hexagonal–icosahedral and cubedsphere grids, Geosci. Model Dev., 7, 909–929, https://doi.org/10.5194/gmd79092014, 2014. a, b, c, d, e, f, g, h, i, j, k
Ullrich, P. A.: A global finiteelement shallowwater model supporting continuous and discontinuous elements, Geosci. Model Dev., 7, 3017–3035, https://doi.org/10.5194/gmd730172014, 2014. a
Ullrich, P. A., Jablonowski, C., Kent, J., Lauritzen, P. H., Nair, R., Reed, K. A., Zarzycki, C. M., Hall, D. M., Dazlich, D., Heikes, R., Konor, C., Randall, D., Dubos, T., Meurdesoif, Y., Chen, X., Harris, L., Kühnlein, C., Lee, V., Qaddouri, A., Girard, C., Giorgetta, M., Reinert, D., Klemp, J., Park, S.H., Skamarock, W., Miura, H., Ohno, T., Yoshida, R., Walko, R., Reinecke, A., and Viner, K.: DCMIP2016: a review of nonhydrostatic dynamical core design and intercomparison of participating models, Geosci. Model Dev., 10, 4477–4509, https://doi.org/10.5194/gmd1044772017, 2017. a
Williamson, D. L., Drake, J. B., Hack, J. J., Jakob, R., and Swarztrauber, P. N.: A standard test set for numerical approximations to the shallowwater equations in spherical geometry, J. Comp. Phys., 102, 211–224, 1992. a, b, c
Wood, N., Staniforth, A., White, A., Allen, T., Diamantakis, M., Gross, M., Melvin, T., Smith, C., Vosper, S., Zerroukat, M., and Thuburn, J.: An inherently massconserving semiimplicit semiLagrangian discretization of the deepatmosphere global nonhydrostatic equations, Q. J. Roy. Meteor. Soc., 140, 1505–1520, https://doi.org/10.1002/qj.2235, 2014. a, b, c
Zerroukat, M., Wood, N., Staniforth, A., White, A. A., and Thuburn, J.: An inherently massconserving semiimplicit semiLagrangian discretisation of the shallowwater equations on the sphere, Q. J. Roy. Meteor. Soc., 135, 1104–1116, https://doi.org/10.1002/qj.458, 2009. a, b