Choice of function spaces for thermodynamic variables in mixed finite‐element methods

We study the dispersion properties of three choices for the buoyancy space in a mixed finite‐element discretization of geophysical fluid flow equations. The problem is analogous to that of the staggering of the buoyancy variable in finite‐difference discretizations. Discrete dispersion relations of the two‐dimensional linear gravity wave equations are computed. By comparison with the analytical result, the best choice for the buoyancy space basis functions is found to be the horizontally discontinuous, vertically continuous option. This is also the space used for the vertical component of the velocity. At lowest polynomial order, this arrangement mirrors the Charney–Phillips vertical staggering known to have good dispersion properties in finite‐difference models. A fully discontinuous space for the buoyancy corresponding to the Lorenz finite‐difference staggering at lowest order gives zero phase velocity for high vertical wavenumber modes. A fully continuous space, the natural choice for scalar variables in a mixed finite‐element framework, with degrees of freedom of buoyancy and vertical velocity horizontally staggered at lowest order, is found to entail zero phase velocity modes at the large horizontal wavenumber end of the spectrum. Corroborating the theoretical insights, numerical results obtained on gravity wave propagation with fully continuous buoyancy highlight the presence of a computational mode in the poorly resolved part of the spectrum that fails to propagate horizontally. The spurious signal is not removed in test runs with higher‐order polynomial basis functions. Runs at higher order also highlight additional oscillations, an issue that is shown to be mitigated by partial mass‐lumping. In light of the findings and with a view to coupling the dynamical core to physical parametrizations that often force near the horizontal grid scale, the use of the fully continuous space should be avoided in favour of the horizontally discontinuous, vertically continuous space.

We study the dispersion properties of three choices for the buoyancy space in a mixed finite-element discretization of geophysical fluid flow equations. The problem is analogous to that of the staggering of the buoyancy variable in finite-difference discretizations. Discrete dispersion relations of the two-dimensional linear gravity wave equations are computed. By comparison with the analytical result, the best choice for the buoyancy space basis functions is found to be the horizontally discontinuous, vertically continuous option. This is also the space used for the vertical component of the velocity. At lowest polynomial order, this arrangement mirrors the Charney-Phillips vertical staggering known to have good dispersion properties in finite-difference models. A fully discontinuous space for the buoyancy corresponding to the Lorenz finite-difference staggering at lowest order gives zero phase velocity for high vertical wavenumber modes. A fully continuous space, the natural choice for scalar variables in a mixed finite-element framework, with degrees of freedom of buoyancy and vertical velocity horizontally staggered at lowest order, is found to entail zero phase velocity modes at the large horizontal wavenumber end of the spectrum. Corroborating the theoretical insights, numerical results obtained on gravity wave propagation with fully continuous buoyancy highlight the presence of a computational mode in the poorly resolved part of the spectrum that fails to propagate horizontally. The spurious signal is not removed in test runs with higher-order polynomial basis functions. Runs at higher order also highlight additional oscillations, an issue that is shown to be mitigated by partial mass-lumping. In light of the findings and with a view to coupling the dynamical core to physical parametrizations that often force near the horizontal grid scale, the use of the fully continuous space should be avoided in favour of the horizontally discontinuous, vertically continuous space. motion. In this context, selection of a particular set-up has commonly revolved around the implied presence or otherwise of so-called computational modes. These computational modes typically take the form of grid-scale zigzag patterns that are invisible to the discrete derivative operator of the scheme and hence these patterns fail to propagate. Unchecked perturbations of this kind can pollute simulations by interacting with physical modes of meteorological interest, whilst ad hoc filtering techniques aimed at reducing the impact of spurious modes run the risk of affecting the overall accuracy and efficiency of the scheme. Absence or efficient control of computational modes is seen as a fundamental requirement for grid arrangements in modern atmospheric dynamical cores (Staniforth and Thuburn, 2012). The sign of the group velocity and the dynamical response to forcing are other examples of issues related to the placement of discrete variables on the grid.
Properties of several grid arrangements in the vertical direction have been extensively studied in Thuburn and Woollings (2005) using normal mode expansions (also Thuburn et al. (2002), Cordero et al. (2002), Benacchio and Wood (2016)). Two of the most common arrangements are the Lorenz (1960) staggering, whereby horizontal velocity is vertically co-located with the potential temperature and vertically staggered with respect to the vertical velocity, and the Charney-Phillips staggering (Charney and Phillips, 1953), whereby the potential temperature is co-located with the vertical velocity and vertically staggered with respect to the horizontal velocity. Due to the computational mode of the Lorenz staggering, and notwithstanding its better conservation properties, the Charney-Phillips staggering has emerged as the preferred arrangement by many groups for dynamical cores in comparisons with the Lorenz staggering, e.g. Arakawa and Moorthi (1988), Fox-Rabinovitz (1994), Arakawa and Konor (1996) and Holdaway et al. (2013a). Examples of operational forecast models running with vertical Charney-Phillips staggering are the Global Environmental Multiscale model (Girard et al., 2014) and the Met Office's Unified Model, which also uses C-grid staggering in the horizontal direction (Cullen et al., 1997;Davies et al., 2005;Wood et al., 2014).
Mixed finite-element methods (Cotter and Shipton, 2012) have a number of properties that make them appealing for modelling geophysical flows: they are inf-sup stable and mimetic in that they carry over continuous properties, such as irrotationality of the gradient operator or solenoidality of the curl operator, at the discrete level (see also Natale et al. (2016) and references therein). Compatible finite-element methods are useful for building discretizations for numerical weather prediction since they allow (a) flexibility in the type of grids that can be used, since they do not require orthogonal grids for their fundamental properties, (b) flexibility in the choice of finite element spaces leading to discretizations on triangular meshes that maintain the 2:1 ratio of velocity to pressure degrees of freedom that is necessary to avoid spurious modes (Staniforth and Thuburn, 2012), and (c) flexibility in the choice of the consistency order of the approximation through using higher-order polynomials. The grid flexibility is particularly important since it allows the recovery of C-grid wave dispersion properties on pseudo-uniform grids such as the cubed sphere that avoid the parallel scalability issues associated with the poles of a latitude-longitude grid. It also allows for unstructured adaptive mesh refinement. Other applications of continuous and discontinuous Galerkin methods in numerical weather prediction can be found in, e.g. Kelly and Giraldo (2012), Marras et al. (2013), Guerra and Ullrich (2016), and, exclusively for the vertical direction, in Simarro and Hortal (2012) and Yi and Park (2017), as well as in the operational model of the European Centre for Medium-range Weather Forecasts (Untch and Hortal, 2004). At lowest polynomial order on a horizontal quadrilateral grid, the finite-element equivalent of the C-grid staggering for wind and pressure variables is the RT0-DG0 pair of mixed finite-element function spaces that has cellwise continuous normal components for the velocity and cellwise discontinuous pressure. However, this leaves the question of what space should be used for the buoyancy variable.
This article compares three choices for the function space for the buoyancy variable in a mixed finite-element scheme for atmospheric dynamics. At lowest polynomial order, with the RT0-DG0 choice for velocity and pressure, the first option is to place the buoyancy in a space with basis functions discontinuous at cell edges in both the horizontal and the vertical direction, corresponding to the Lorenz vertical staggering. The second option, with horizontally discontinuous, vertically continuous basis functions, corresponds to the Charney-Phillips staggering. A third option of a space with continuous basis functions in both directions corresponds to the natural choice for scalars in the so-called deRham complex of mixed finite-element function spaces (Thuburn and Cotter, 2015 and references therein). Moreover, the fully continuous option eases integration by parts of the pressure gradient term in the momentum equation of the fully compressible Euler model (Natale et al., 2016). Unlike in the first two cases, with the third option the buoyancy degrees of freedom are horizontally staggered with respect to those of the vertical velocity.
The three options are evaluated using a discrete dispersion analysis on a weak formulation of the two-dimensional linear gravity wave equations. With buoyancy in the fully continuous space, it is found that gravity waves are slowed down on a large range of horizontal wavenumbers, with the shortest two-gridpoint wave failing to propagate. The analysis also confirms the issues with the Lorenz vertical staggering on the large vertical wavenumber end of the spectrum, while the horizontally discontinuous, vertically continuous option gives the best match with the analytical dispersion relation for gravity waves. The theoretical findings are corroborated by numerical simulations of a non-hydrostatic gravity wave in a horizontal channel. Horizontally co-locating the buoyancy with the vertical velocity gives results in line with the literature, whilst a spurious zigzag pattern appears in the coarse-grid simulations using the fully continuous space. Furthermore, a simulation using one-degree higher polynomial basis functions highlights widespread oscillations when using the fully continuous space. It is observed that the oscillations can be smoothed out by employing a partial mass-lumping strategy in the spirit of Melvin et al. (2014).
The rest of the article is structured as follows. Section 2 contains the analytical model used in this study together with the analytical dispersion analysis. The weak formulation and the mixed finite-element formulation are given in section 3. The discrete system is explicitly derived for the lowest-order finite-element case in section 3.1, which also contains the discrete dispersion analysis for that case. Section 3.2 discusses the next higher-order case. Section 4 contains the numerical results, which are discussed together with concluding remarks in section 5.

CONTINUOUS EQUATIONS
The linear gravity wave equations in a two-dimensional Cartesian vertical slice are: for the prognostic velocity vector u ≡ (u, w), pressure field p, and buoyancy b. Moreover, subscript t denotes time derivatives, c s denotes the speed of sound, N the buoyancy frequency (both assumed constant) andẑ ≡ (0, 1) is a unit vector in the vertical direction. Equations 1-3 can be derived from the compressible equations by linearization and neglecting the effects of stratification on compressibility (Durran, 2010). The analytical dispersion relation is reviewed here as a reference for the discrete versions which will be derived in the following sections. As the equations are linear and have constant coefficients, to analyze the dispersion relation we seek harmonic solutions in the form: yielding the system: which can be satisfied if and only if the determinant of the matrix vanishes. Hence, admissible solutions satisfy the dispersion relation: which is a biquadratic equation for , with the four solutions: (6) The solution with positive sign in Equation 6 corresponds to acoustic waves, and the one with negative sign to gravity waves (Durran, 2010). Expression 6 can be simplified under the assumption N 2 ≪ k 2 c 2 s , which holds in most atmospheric applications. In that case, the dispersion relation for acoustic waves becomes: and the one for gravity waves, using the first-order expansion

MIXED FINITE-ELEMENT DISCRETIZATION
Writing the velocity equation componentwise, the weak formulation of Equations 1-3 reads, on a two-dimensional domain Ω, find u, w, b, p such that for all test functions and in the horizontal and vertical velocity spaces, in the pressure space, and in the buoyancy space.
The discrete finite-element formulation is obtained from the weak formulation 9-12 by placing the unknowns and test functions in the target finite-element spaces. In particular, we consider a compatible mixed finite-element formulation (Cotter and Shipton, 2012) that seeks discrete solutions in separate function spaces for the different variables. For r ≥ 0, velocities are sought in the discrete space such that both the u and w components are polynomial functions that are continuous and cellwise of degree r+1 in the normal direction, and discontinuous and cellwise of degree r in the tangential direction. This means that u is continuous and cellwise of degree r + 1 in the horizontal direction, discontinuous and cellwise of degree r in the vertical direction, while w is continuous and cellwise of degree r + 1 in the vertical direction, discontinuous and cellwise of degree r in the horizontal direction. Pressures will be sought in the space of scalar polynomial functions that are discontinuous functions of degree r in both directions, denoted V 2 (Cotter and Shipton, 2012). As for the buoyancy, we consider three possible choices for its function space:

3.1
Lowest-order r = 0 elements For r = 0 lowest-order elements, the one-dimensional linear and constant basis functions for an interval s ∈ [ s m , s m+1 ] are, following the notation of Melvin et al. (2014), The two-dimensional basis functions can then be created as tensor products of the appropriate one-dimensional functions; these sets are: The variables, restricted to a single element [ are expanded as: Inserting Equations 22-27 into 9-12 with test functions 16-21, one gets two discrete equations for the horizontal momentum equation, two equations for the vertical momentum equation for each choice of the buoyancy space V 0 , V cp and V 2 , one equation for the pressure, and four, two, or one equations for the buoyancy variable for the choices V 0 , V cp , and V 2 , respectively. Calculating the integrals of products of basis functions within these integrals, considering a uniform grid spacing Δx ≡ x m+1 − x m and Δz = z n+1 − z n in both directions and integrating over two neighbouring elements Location of degrees of freedom, r = 0 order elements for buoyancy in (a) V 0 space, (b) V cp space, and (c) V 2 space, and for (d) horizontal and (e) vertical components of velocity and (f) the pressure field. m and n are the horizontal and vertical indices for cell (m, n) such that m + 1∕2, n + 1∕2 is located at the centre of the cell and simplifying, one gets the following expression for the horizontal momentum equation: For the vertical momentum equation integrating over two neighbouring elements spanning , with b ∈ V 0 one gets: For the vertical momentum equation, with b ∈ V cp one gets: For the vertical momentum equation, with b ∈ V 2 one gets: Integrating while integrating over four neighbouring elements spanning , the buoyancy equation in the V 0 case reads: Integrating over two neighbouring elements spanning for the buoyancy in V cp , one gets: Finally, integrating over a single element spanning for the buoyancy in V 2 one gets: For buoyancy in V cp , as can be seen in expression 34, the coefficients of w and b match. By contrast, horizontal and vertical averaging is evident in the expressions 33 and 35 for V 0 and V 2 , respectively (also Figure 2). As a consequence, in the two latter cases the scheme will fail to detect and propagate the shortest resolved waves.

Discrete dispersion analysis
Solutions are sought in the form: with k = (k, l) and x m,n = (x m,n , z m,n ). Inserting Equations 36-39 in 28-35 and using the properties of the exponential, after simplification and grouping the equations in matrix form we have: with the switches: where: , ) .
(42) It is worth noting that the Ms represent the action of mass matrices, Ss are the action of the discrete derivative and the Cs represent averaging, or equivalently projection from one finite-element function space to another. For the left-hand side of Equation 40 to vanish, the determinant of the matrix needs to vanish, i.e. expanding with respect to the last row and simplifying we have: Solving, we find the two roots for = 2 : (44) Figure 3 shows the discrete gravity wave dispersion relation (positive square root of 2 ) for Δx = Δz = 1 000 m. With the horizontally discontinuous, vertically continuous buoyancy space V cp , the analytical behaviour of − in Figure 1b is well reproduced. This is the desired behaviour, which reproduces the properties of the well-known Charney-Phillips staggering in the finite-difference context. With the fully discontinuous space V 2 , the analytical behaviour is well reproduced for small vertical wavenumbers lΔz, while the scheme fails to propagate waves with large lΔz. This is the well-known shortcoming of the Lorenz staggering with grid-scale vertical motions. The situation is reversed with the fully continuous space V 0 , where the analytical behaviour is well reproduced for very small horizontal wavenumbers kΔx only, while virtually every signal with kΔx > ∕4 is slowed down. Grid-scale short waves with kΔx approaching are not propagated by the scheme. Except for an overshoot for large wavenumbers, the discrete acoustic wave dispersion relation (positive square root of 1 ; Appendix A, Figure A1a) reproduces the analytical dispersion of acoustic waves + in Figure 1a for each choice of the finite-element space for the buoyancy, and so only the result for V 0 is shown in Figure A1.

3.2
Higher-order r = 1 elements For r = 1 order elements, the one-dimensional basis functions are now quadratic (E) and linear functions (F), which, for an The two-dimensional basis functions can then be created as tensor products of the appropriate one-dimensional functions. These sets are 1,2,3,4 ≡ 1,2,3,4 .
The locations of degrees of freedom in the three buoyancy spaces under consideration as well as for the other variables are as in Figure 4. The variables, restricted to a single element w(x, z) = w m,n +,0 1 +w m+1,n −,0 2 +w m,n+1∕2 +,0 3 +w m+1,n+1∕2 −,0 4 +w m,n+1 p(x, z) = p m,n +,+ 1 + p m+1,n −,+ 2 + p m,n+1 where 0, +, − subscripts denote whether the degree of freedom is continuous (0) or discontinuous (±) in the x and z directions and located to the positive or negative side of a point, i.e. p m,n+1 is discontinuous in x and z and located at point . For all three choices of the buoyancy space, Equations 9-12 can now be written in operator form as a sum over all elements e: where the superscript (e) denotes restriction of a prognostic field to a single element e. The operators are where now each of the matrices  u,w,p,b ,  1,2 ,  is a function of the spatial wavenumber k ≡ (k, l) due to the substitutions of Equations 36-39. For a given wavenumber pair (k, l) , the system 68 can be solved to yield the frequencies j , j = 1 … 16. These solutions correspond to the sets where w 1,…,4 are the eigenvalues associated with the vertical velocity components W of Equation 68 and w + j representing the acoustic waves and w − j representing the gravity waves, j = 1 … 4. For the positive roots of + and − , there are now four solutions instead of the single solution found for both the analytical case and the r = 0 solutions. The question is what do these extra solutions represent?
As in  and Melvin et al. (2014), here it is argued that these solutions are in fact aliased waves from higher wavenumbers that are now resolvable due to the finite-element method being able to support modes that oscillate on the sub-grid scale. For example, consider the pressure field. For the lowest order (r = 0) elements, there is only one degree of freedom per element p m+1∕2,n+1∕2 and this is associated with a piecewise constant function. Therefore, the highest wavenumber oscillation that can be supported is on the grid scale where p m+1∕2,n+1∕2 = 1, p m+3∕2,n+1∕2 = −1. However, for r = 1 there are now four degrees of freedom per element ] and these are associated with piecewise linear functions. Therefore, the highest wavenumber oscillation that can be supported is of a wavenumber at twice the grid scale, . By this argument, the four positive solutions to each of + and − correspond to a solution in the range (kΔx, lΔz) ∈ (0, ) × (0, ) along with solutions aliased from regions with a subgrid frequency in x (kΔx, lΔz) ∈ ( , 2 ) × (0, ), in z (kΔx, lΔz) ∈ (0, ) × ( , 2 ), and in both x and z (kΔx, lΔz) ∈ ( , 2 ) × ( , 2 ). For each wavenumber pair (k, l) we must determine how to correctly assign each of the four roots ± j to the correct quadrant of the extended wave number space. The method used here is as follows.
Consider a solution of (68) compactly rewritten as A V = 0 where V ≡ (U, W, P, B) and A is the dispersion matrix that depends upon the frequencies . Then let Y be one of U, W, P or B. For each eigenvalue of the dispersion matrix A there exists the associated eigenvector Y such that A * Y = Y with A * being the reduced 4×4 matrix that only contains the rows and columns of A that correspond to the components of Y. If the eigenvector Y belongs to the first quadrant of the extended wavenumber region, (kΔx, lΔz) ∈ (0, ) × (0, ), then there will be no subgrid variation and hence all values of Y will have the same sign: sign (Y) = [1, 1, 1, 1]. However if the mode is aliased from one of the other quadrants of the extended wavenumber space, then there will be a subgrid variation in either or both directions and the values of Y will not all have the same sign, i.e.
Using these expressions, each of the four roots for + and − can then be placed into the correct region of the extended wavenumber space by choosing Y to be U, W, P and B in turn and recording the quadrant of the extended wavenumber space that the eigenvalue lies in. Generally the procedure returns all four choices for Y placing the eigenvalue in the same quadrant. However occasionally the method does not work as well, in which case the quadrant is chosen by taking the most popular option from the four choices of Y. The results for the three choices of the buoyancy space are presented in Figure 5 for − and Figure A1b for + . Note that in both cases, to allow easier comparison with the previous analytical and r = 0 results, the wavenumber space has been scaled to [ 0, 2 ∕Δx ] 2 . This is equivalent to comparing results with r = 0 elements with grid spacing (Δx, Δz) to r = 1 elements with grid spacing (2Δx, 2Δz), i.e. the spacing between the degrees of freedom has been kept constant at (Δx, Δz) between the low-and higher-order elements. The dispersion relation for gravity waves ( − , Figure 5) shows similar results to the r = 0 elements. In particular, the dispersion issues observed at lowest order are not removed by using higher-order polynomials. With b ∈ V 0 horizontally short waves are severely retarded while with b ∈ V 2 vertically short waves are severely retarded. There are discontinuities in the dispersion relation between the four branches of solutions. However it should be noted that in the limit l = 0, both the V cp and V 2 spaces reproduce the exact solution − = N as kΔx → ∞ (Figure 5e, f) and so there is no discontinuity present in the solution. In contrast the V 0 space exhibits the discontinuity even in this limiting case because of the different horizontal representation of the buoyancy and vertical velocity variables with b ∈ V 0 . This is the same behaviour as found in  for the pure inertia waves. In general the dispersion properties for the r = 1 elements show a small improvement in accuracy compared to the r = 0 elements. In particular for the gravity waves, the frequency drops off more slowly as l is increased. The process of automatically placing frequencies in the correct quadrant of the extended wavenumber space is difficult and the method outlined above is not guaranteed to be successful. In fact, it breaks down for gravity waves with b ∈ V 0 and vertical wavenumber lΔz ≈ ∕2 as can be seen in Figure 5a where some frequencies have been misplaced in the incorrect quadrant, leading to jumps in the dispersion relation. While discontinuities in the dispersion relation are expected between the four branches when kΔx or lΔz = , within each section the dispersion relation should be continuous. To mitigate this error, a simpler one-dimensional semi-discretization (continuous in z, discrete in x) has been performed and is presented in Appendix B, see Figures B1 and B2. These results agree qualitatively with those in Figures 5 and A1b; due to the semi-discrete nature they would not be expected to produce the same results.
The dispersion relation for acoustic waves ( + ; Figure A1b, d) is similar to that found for the r = 0 elements and also, as expected, for inertia-gravity waves in the linear shallow-water results from Melvin et al. (2014) when the reference depth Φ 0 ≈ c 2 s . As in Melvin et al. (2014), there are discontinuities in the dispersion relation between the four branches of solutions that lie along lines kΔx = and lΔz = and, as in section 3.1, there is very little difference between the choices for the different spaces used for the buoyancy variable; this is expected since N 2 ≪ k 2 c 2 s and so the buoyancy terms are negligible in this regime.

NUMERICAL RESULTS
In order to corroborate the theoretical findings of the previous section, a mixed finite-element numerical scheme is run on a wave propagation benchmark. To this end, we consider the non-hydrostatic gravity wave experiment of Skamarock and Klemp (1994), which tests the handling of horizontally propagating non-orographic gravity waves in a two-dimensional Cartesian domain [−150, 150] × [0, 10] km. The initial condition is a balanced, thermally stratified reference state with a superposed buoyancy perturbation. For the linear system of equations studied here, this translates into a motionless basic state with a superposed buoyancy perturbation of the form: with b 0 = 0.01 m s −2 , H = 10 km, A = 5 km, and centred on x c = 0. The governing Equations 1-3 are discretized using a finite element method as set out in section 3 and a (possibly off-centred with an off-centring parameter oc ) implicit temporal discretization. Defining f ≡ f n+1 − f n for a field f , expressions 1-3 become In order to evaluate the impact of different grid spacings, the model is run with two different spatial resolutions. Specifically, at lowest finite-element order r = 0 we consider the spacings Δx = 1 km and Δx = 3 km. The vertical grid spacing is kept fixed at Δz = 1 km, the focus being on horizontal wave propagation. For all the runs we set oc = 0.5. The time step is Δt = 10 s, and the final time is T = 3 000 s. Figure 6 displays the numerical results for lowest order r = 0. The initial perturbation spreads into buoyancy-driven inertia-gravity waves. The wave propagation is well reproduced for buoyancy in V cp and V 2 , with only minor differences observed in the smallest amplitude waves. For fully continuous buoyancy b ∈ V 0 , the situation is different. For Δx = 3 km (Figure 6b) a large oscillatory pattern is evident around x = 0, highlighted in the one-dimensional cut through mid-height z = 5 km in Figure 8a. At resolution Δx = 1 km, the expected behaviour is found for all choices for the buoyancy space. As observed in the previous section and Figure 3a,d, in the fully continuous case the scheme is unable to propagate modes with large kΔx. Figure 7 displays the numerical results for the next higher order r = 1. The horizontal grid spacings considered are Δx = 2 km and Δx = 6 km, so as to keep the number of degrees of freedom the same as in the case r = 0. The vertical grid spacing is Δz = 2 km. The wave pattern is well reproduced in all cases, except for Figure 7b for b ∈ V 0 with Δx = 6 km. Figure  7b is different in two respects. As in the lowest-order case, a computational mode is evident in the centre. In addition, noise appears to affect the waves away from the centre. This is evident in the most external crests in the two-dimensional plots as well as in the one-dimensional cut in Figure 8b. At first glance, the result is somehow counterintuitive. By increasing the finite-element order, one would expect a less prominent stationary mode. In terms of root mean square error with respect to a fine-grid solution (not shown), the r = 1 result is no less accurate than the r = 0 result. Nonetheless, the propagating noise remains a disturbing feature. Viewing the dispersion relation in Figure 5a,d with this result in mind, a possible cause is the discontinuity in the dispersion relation which exists for all vertical wave numbers l with the V 0 space, in contrast to the V cp and V 2 spaces. In the situation presented in Figure 7b,e, the vertical motion is fairly well resolved (lΔz is small) but the horizontal motion is poorly resolved (kΔx is large). Hence, it is likely that the initial perturbation is projecting onto wavenumbers around the discontinuity which in turn results in unexpected behaviour that is absent from both the runs with the lower order V 0 space (where there is no discontinuity) and with the higher order V cp and V 2 spaces (where the discontinuity vanishes in the limit lΔz = 0). A proposed solution to this discontinuity problem is detailed in  and involves partially lumping the mass matrix of the momentum equation. This has been adapted here to a partial lumping of the mass matrix of the V 0 where now the partially lumped test functions are given by * and similarly for * 2,3,…,9 . The new one-dimensional basis functions G (cf. Equations 45-47) are and the coefficient = 1∕10. These basis functions along with the unmodified set are shown in Figure 9. This method of partial lumping (through the modified test functions) is designed to reproduce both the one-dimensional partial lumping result of  and the two-dimensional result of Melvin et al. (2014). The resulting dispersion relation is shown in Figure 10. The curves in Figure 10 show that the partial lumping has the effect of closing the gap in the limit of lΔz → 0, as can be seen in Figure 10d. Rerunning the test from Figure 7b, the partial lumping reduces the noise seen in Figure 7b and, see Figure 11a, which now closely approximates the lower order r = 0 solution. However, note that, importantly, the unexpected signal in the centre remains. Also note that, although the discontinuity at kΔx = for the gravity wave has been reduced, the discontinuity at lΔz = remains, as do the discontinuities for the acoustic waves, Figure 10a,c. To fully remove all discontinuities, the mass matrices for the velocity components M u and M w in Equation 66, along with Q where it appears in the vertical momentum_equation, also need a partial lumping applied, as shown in Melvin et al. (2014).
The partial mass lumping does not affect the overall accuracy of the scheme since, although it reduces the accuracy in the horizontal of the temporal derivative in the buoyancy equation from r + 2 to r + 1 order accuracy (third order to second order), the overall accuracy of the scheme will be limited by the piecewise linear representation of terms such as the pressure which will be r + 1 (second) order accurate, and so overall the scheme is still r + 1 (second) order accurate.

DISCUSSION AND CONCLUSION
The discrete dispersion properties of three choices for the buoyancy space in a mixed finite-element discretization of geophysical fluid flow were investigated. For the lowest-order analysis on the two-dimensional gravity wave equations, velocity and pressure were placed in the RT0-P0 pair of function spaces. The best match with the analytical dispersion of gravity waves was achieved with a horizontally discontinuous, vertically continuous space for the buoyancy, which corresponds to the Charney-Phillips vertical staggering. The second option of a fully discontinuous buoyancy space, corresponding to the Lorenz vertical staggering, gave a dispersion relation with vanishing phase velocity for high vertical wavenumbers. A third, fully continuous option, corresponding to the natural space for scalars in the deRham complex of mixed finite-element spaces, was shown to slow down wave frequencies on a sizeable range of horizontal wavenumbers. These findings also broadly extended to higher-order r = 1 elements, notably the fully continuous space was shown to still slow down wave frequencies over a smaller but still sizable range of horizontal wavenumbers. Additionally, as found in earlier studies of the shallow-water equations, unphysical spectral gaps exist between the extra branches of solutions now permitted by the higher-order elements. If the horizontal representation of the buoyancy space matches that of the vertical momentum space, i.e. for b ∈ V cp , V 2 , then in the limit lΔz → 0 the exact representation of the gravity waves is obtained. By contrast, for b ∈ V 0 , only the well-resolved half of the spectrum is exact, whilst for high wavenumbers the frequency is spuriously retarded to zero as in the case for r = 0 elements with b ∈ V 0 . The outcome of the theoretical analysis was complemented by a numerical comparison of the three choices on simulations of gravity wave propagation in a vertical slice. For a coarse grid run at lowest order with the fully continuous space, a spurious zigzag grid-size pattern appeared in the central part of the domain which became decoupled from the dynamics as the scheme was unable to propagate the large horizontal wavenumber end of the spectrum. By contrast, results in line with the literature were obtained using the horizontally discontinuous spaces. The pattern observed in Figure 6b would have potentially serious effects in a numerical weather prediction model. Left unchecked, the spurious signal would be bound to interact with forcing from the physics parametrizations at the same horizontal scale in an uncontrollable fashion, inevitably affecting the overall accuracy of coupled simulations, even for motions well-resolved in the vertical direction.
Numerical results using one degree higher polynomials confirmed the presence of the computational mode for the coarse grid run in the fully continuous space. This was in some way unexpected as, at least at an intuitive level, the higher the polynomial order, the less the horizontal distance between degrees of freedom of buoyancy and vertical velocity. In the same simulation, small oscillations were shown to affect the wave dynamics at all scales. A partial mass-lumping strategy as used in Melvin et al. (2014) was shown to be effective in closing the gap and in smoothing out the oscillations. Note that, although using higher-order elements results in more terms in the basis function expansions for a variable, cf. Equations 22 and 56, so that an increased number of operations must be performed to evaluate a given term (particularly for nonlinear terms), this does not necessarily entail degrading computational performance in large-scale problems. In fact, due to the balance of costs between computation and communication, higher-order methods can perform extremely well on large-scale problems, e.g. Nair et al. (2009).
The present study has considered only the case of Cartesian grids of quadrilaterals in two spatial dimensions and used linear equations. Nonetheless, the gathered findings appear sufficient to recommend relinquishing the fully continuous option for the discretization of the temperature field in compatible finite-element discretizations for geophysical flow. The horizontally discontinuous Charney-Phillips and Lorenz-like options should be favoured instead, based upon the case studied here, with the Charney-Phillips-like option providing the most accurate dispersion properties. In a nonlinear setting, the horizontally discontinuous choice implies a technical overhead in that it requires the computation of boundary integral terms on horizontal facets that arise from the integration by parts of the pressure gradient term in the momentum equation and of the advective term in the potential temperature equation in a fully compressible model (Natale et al., 2016).
The results here extend the earlier insights in the dispersion studies of  and Melvin et al. (2014) and are part of the current effort in the development of a scalable dynamical core for the next generation of atmospheric models at the Met Office . Thanks to its grid-agnostic mimetic properties, the mixed finite-element methodology is being employed to generalize to non-orthogonal grids the accuracy properties of current models. Analysis of results with a three-dimensional compressible implementation on Cartesian and spherical test cases will be the subject of future publications.

APPENDIX A: ACOUSTIC WAVE DISPERSION RELATION
The discrete dispersion relation for the acoustic wave branch of the dispersion relation is shown in Figure A1 for (a,c) r = 0 elements and (b,d) r = 1 elements. In both cases only the result for the V 0 space is shown as results for the other spaces are visually identical.

APPENDIX B: SEMI-DISCRETE DISPERSION ANALYSIS
The noise that appears in the r = 1 elements when b ∈ V 0 at coarse resolution, Figure 7b, can be investigated in a simpler setting by using a semi-discrete approach. Since the problem is theorized to be a result of the different representations of the buoyancy space in the horizontal direction, then a one-dimensional dispersion analysis that is semi-discrete in x will reveal the problem. Substituting solutions that are harmonic in the z direction of the form (u, w, p, b) into Equations 1-3 yields the system which is a system of a single dependent variable x. In one dimension the finite element function spaces used in section 3.2 reduce to either a quadratic continuous space (Ũ) or a linear discontinuous space (P,W) and the spaces used for the buoyancy variableB are the quadratic continuous space (b ∈ V 0 ) or the linear discontinuous space (b ∈ V cp , V 2 ). It should be noted that these spaces correspond to those studied in  and Melvin et al. (2014). ) (x) = (U, W, P, B) exp (ikx) (but where all the vertical variations are neglected) to the set of Equations B1-B4 leads to a system of eight equations for the unknowns (U, W, P, B) which corresponds to two roots for each of ± + and ± − . Again, one root corresponds to a branch that has been aliased from higher wavenumbers but, in contrast to the two-dimensional case, it is straightforward to separate the two branches. The positive roots of + and − are shown in Figure B1 for a vertical wavenumber of l = ∕(4Δz). For the gravity waves when b ∈ V cp , V 2 there is an accurate representation of the dispersion relation (in fact is shown as dashed, and the semi-discrete solution as solid the discrete relation overlies the analytical results), whilst for b ∈ V 0 , although half the spectrum (for kΔx < ) is again very accurate, the other half of the spectrum (for kΔx > ) is poorly represented and there is a discontinuity between the two branches. To remove this discontinuity the partial lumping of section 4 is successfully applied, Figure B2a, at the cost of slightly degrading the accuracy of the discrete dispersion for kΔx < ∕2. The discontinuity for the acoustic waves, Figure B2b, is unaffected, but this could be removed by using the same partial lumping approach to the mass matrix in the u equation as shown in .