Numerical entropy conservation without sacrificing Charney–Phillips grid optimal wave propagation

Conservation of entropy or potential temperature and accurate representation of wave propagation without computational modes are both desirable properties for a numerical model of the atmosphere. However, they appear to require different model formulations, forcing model developers to choose between them. Here it is shown that, by a straightforward modification of the horizontal entropy fluxes, numerical entropy conservation can be achieved without sacrificing accurate wave propagation. The result is confirmed by a numerical linear normal mode analysis for a simple but suitably modified finite‐volume scheme, and by buoyant bubble and gravity‐wave test cases in a vertical slice model using a suitably modified conservative semi‐Lagrangian transport scheme.


INTRODUCTION
Conservative transport of entropy is a desirable property for the dynamical core of an atmospheric numerical model. So, too, is an accurate representation of wave propagation without computational modes. However, these two properties seem to require different model formulations, forcing model developers to choose between them. In this note, a scheme for obtaining numerical entropy conservation without sacrificing optimal wave propagation is proposed and tested.
On a global scale, numerical entropy conservation is important because the entropy budget constrains the behaviour of the climate system (Goody, 2000), and it has been suggested that spurious numerical entropy production could lead to systematic biases in numerical models (Johnson, 1997). The entropy budget is closely related to the budget of available potential energy (Lorenz, 1955;Peixoto and Oort, 1992), which controls the strength of midlatitude eddies and other aspects of the circulation. The entropy budget is also important for a variety of smaller scale phenomena such as the growth of the convective boundary layer (e.g., Stull, 1988), and in precipitating convection (Emanuel et al., 1994;Pauluis and Held, 2002;Raymond, 2013). Numerical difficulties in simulating the convective boundary layer in a single-column two-fluid model, associated with the use of a nonconservative transport scheme for entropy (Thuburn et al., 2019), were a key motivation for the work presented here.
It is relatively straightforward to achieve conservation of entropy (or some related quantity such as potential temperature) in a dynamical core by formulating and discretizing its prognostic equation as a flux-form conservation law. Moreover, the developer retains considerable flexibility in how the fluxes are chosen, allowing higher-order accuracy, upwinding, and monotonicity constraints, for example.
The vertical placement of variables on the model grid can significantly affect the properties of a numerical model, including conservation properties, the propagation of marginally resolved waves, and the ability to represent balanced flows. Two main alternatives are Charney-Phillips vertical staggering (Charney and Phillips, 1953), in which the entropy is staggered vertically relative to the density and horizontal velocity components u and v, and Lorenz vertical staggering (Lorenz, 1960), in which the entropy is located at the same vertical levels as the density and horizontal velocity. The extensions of these two grids to the fully compressible nonhydrostatic case are shown in Figure 1.
By predicting at the same grid location as , the Lorenz grid facilitates the development of schemes that conserve entropy and also energy. However, studies of the effects of grid staggering on wave propagation (e.g., Tokioka, 1978;Lesley and Purser, 1992;Fox-Rabinovitz, 1994;1996;Thuburn and Woollings, 2005;Liu, 2008;Girard et al., 2014;Thuburn, 2017b) have shown that the Charney-Phillips grid gives more accurate wave propagation; provided the pressure-gradient term is evaluated appropriately (Thuburn, 2006;Toy and Randall, 2007), the wave propagation is "optimal" in the sense that it is as good as can be achieved by any scheme based on two-point second-order centred differences. Moreover, the Lorenz grid supports a computational mode, that is, a vertical pattern in the thermodynamic variables that spuriously satisfies hydrostatic balance and so is invisible to the dynamics. The existence of the computational mode can lead to the appearance of vertical grid-scale noise, an unphysical response to forcing (e.g., Schneider, 1987), and even spurious baroclinic instability (Arakawa and Moorthi, 1988). The reduced accuracy of wave propagation on the Lorenz grid is associated with impaired adjustment towards hydrostatic and geostrophic balance for disturbances with small vertical scale (e.g., Arakawa and Konor, 1996), and with a reduction in the effective Rossby deformation radius, implying an increased susceptibility to the so-called Hollingsworth instability (Hollingsworth et al., 1983;Bell et al., 2017) for models using the vector-invariant form of the momentum equation.
On the Lorenz grid, the calculation of the buoyancy term in the vertical momentum equation requires to be averaged vertically from its native levels to w-levels. Also, if the entropy equation is written in advective form, then the vertical velocity w must be averaged vertically to -levels to calculate w ∕ z. Less obviously, if the entropy equation is written in conservative form, there is still an implied averaging of w (Appendix A). This averaging is responsible for the less accurate wave propagation, reduced effective Rossby deformation radius, and computational mode of the Lorenz grid, and the avoidance of such averaging is critical for the good wave propagation behaviour of the Charney-Phillips grid.
Is it not therefore possible to obtain conservation of entropy together with accurate wave propagation by using a Charney-Phillips staggering of variables combined with a flux-form discretization of the entropy equation? At first glance, this does not seem to be possible. All of the studies showing accurate wave propagation on the Charney-Phillips grid assume the advective form for the entropy or potential temperature equation. If the entropy equation is written in flux form, then w must be averaged vertically to compute the vertical entropy fluxes (and there is a further implied vertical averaging from the discrete product rule, Appendix A). Thus it appears impossible to avoid unwanted vertical averaging if the entropy equation is solved in flux form, even on a Charney-Phillips grid.
The key to obtaining entropy conservation without losing the optimal Charney-Phillips grid wave propagation is to take a finite-volume perspective of the behaviour of a vertical-grid-scale disturbance, and recognise that the entropy tendency in an -cell should arise primarily through horizontal rather than vertical fluxes (Section 2). This insight then suggests a straightforward modification to a finite-volume entropy-transport scheme that gives it the desired properties. A numerical normal mode analysis of the discrete linearized equations confirms that such a modified finite-volume scheme does indeed give optimal wave propagation (Section 2).
Once this key idea is recognized, it can be adapted and applied to other conservative advection schemes. In Section 3, it is applied to the conservative semi-Lagrangian SLICE scheme of Zerroukat et al. (2007;2009). The procedure is not quite straightforward, because the key idea refers to the calculation of fluxes, whereas SLICE works in terms of remapping. This modified SLICE scheme is then tested in a two-dimensional model (Section 4). The conservation properties are verified in an idealized saturated buoyant bubble test and the accuracy of wave propagation is tested by simulating gravity waves with small vertical scale.

FINITE-VOLUME TRANSPORT
For the rest of this article, we restrict attention to a Charney-Phillips grid in a height-based vertical coordinate. For clarity, we will also restrict attention to the two-dimensional (x, z)-plane and assume uniform horizontal and vertical grid spacing Δx and Δz respectively. In this case, all vertical averages between -levels and w-levels (indicated by an overbar) can be taken to have simple 1 2 -1 2 weights. A spatially discrete conservation law for mass may be written Here, i and are the horizontal and vertical grid indices of the cell of interest. F x i+1∕2 are the horizontal mass fluxes. The index i + 1∕2 indicates that they are evaluated at the lateral faces of the cell, and they will typically be expressed as F x i+1∕2 = u i+1∕2̂i+1∕2 , wherêis a cell face value of that must be reconstructed from the cell average values i . Similarly, F z i +1∕2 are vertical mass fluxes evaluated at the lower and upper cell faces and typically expressed as F z i +1∕2 = w i +1∕2̂i +1∕2 . There is considerable freedom in how thêare chosen. We also require a discrete conservation law for entropy, but, since and are stored at different levels, this must be a conservation law for the quantity . In order that the scheme should be able to preserve an initially uniform , this conservation law must reduce to a conservation law for that is consistent with Equation (1) when ≡ 1. The discrete conservation law for is obtained simply by taking a vertical average of Equation (1): (At the bottom boundary, level 1∕2, the conservation law applies in a layer of thickness Δz∕2 and reduces to with i 1∕2 = i 1 , F x i+1∕2 1∕2 = F x i+1∕2 1 , F z i 1 = 1 2 F z i 3∕2 , and F z i 1∕2 = 0. A similar modification is made at the top boundary.) We now seek a discrete flux-form conservation law for . Its general form must be (suitably modified at the lower and upper boundaries), where G x and G z are the horizontal and vertical entropy fluxes. Based on Equation (2), we might anticipate that G x and G z must be related to F x and F z . Let us first clarify why a naive choice for G x and G z does not give optimal wave propagation. Consider the situation shown in Figure 2a. Suppose that there is a background stratification in which increases with height, and that the mass fluxes F x and F z have an oscillation with vertical scale 2Δz. In this situation, the 2Δz structure in the fluxes should lead to a 2Δz structure in the tendencies. This behaviour is captured correctly by the advective-form equation via the w ∕ z term, which involves no averaging of w. However, we wish to use the flux-form , Equation (4). Suppose the entropy fluxes are defined by for some reconstructed -cell face valueŝ. It is clear that if the mass fluxes F x and F z have a vertical 2Δz oscillation, as in Figure 2a, then the vertically averaged fluxes F x and F z will vanish and so, too, will the entropy fluxes (Equations (5) and (6)). Consequently the tendencies must vanish, and the correct behaviour is not captured.
To make progress, let us examine the entropy budget for the -cell shown by the dotted line in Figure 2a. Because of the descent at the cell centre, the cell-average value of should increase. However, this increase cannot occur via vertical fluxes through the lower and upper cell faces, because the mass fluxes F z vanish there. The only possibility, then, is that the increase in cell-average occurs via horizontal fluxes. Now, the net mass flux at the right face of the -cell F x i+1∕2 +1∕2 is zero. However, it is made up of two nonzero but cancelling contributions: 1 2 F x i+1∕2 on the lower part of the face and 1 2 F x i+1∕2 +1 on the upper part of the face. Because of the background stratification, F x i+1∕2 +1 ought to carry a greater entropy flux into Schematics showing the implied vertical subgrid reconstruction of the profiles of F x (z) and̂(z) in the derivation of Equation (7). Open circles indicate the values of F x at -levels and filled circles indicate the values of̂at w-levels obtained by horizontal reconstruction. Note that the piecewise constant vertical reconstruction of̂between z and z +1 preserves the vertical average of in that interval equal tô+ 1∕2 . (d) Schematic showing the implied vertical subgrid reconstruction of (z) in the derivation of Equation (21). Note that the piecewise constant reconstruction ofb etween z and z +1 does not preserve the vertical average of̂in that interval equal tô+ 1∕2 the -cell than F x i+1∕2 carries out. Thus, there should be a net horizontal flux of entropy into the -cell even though the net horizontal mass flux vanishes.
The situation just described can be captured by a straightforward modification of the horizontal entropy fluxes: Here, ∕ z is an estimate for the vertical derivative of at the lateral cell faces, obtained, for example, using a finite-difference approximation.
One way to obtain the expression in Equation (7) is as follows. For the flux G x i+1∕2 +1∕2 , approximate the vertical profile of F x i+1∕2 (z) by a piecewise constant subgrid reconstruction, constant between neighbouring pairs of w-levels, and the vertical profile of̂i +1∕2 (z) by a piecewise constant subgrid reconstruction, constant across each half-interval:̂i ( Figure 2b,c). The entropy flux is then given by Evaluating the integral then gives Equation (7). 1 1 As an alternative, we may use a (discontinuous) piecewise linear reconstruction for̂: and approximate the integral using the trapezium rule, which again leads to Equation (7). If, however, we use this piecewise linear reconstruction and evaluate the integral exactly, then we obtain an expression similar to Equation (7) but with a factor of 1∕8 rather than 1∕4. However, only with a factor of 1∕4 does the implied advective form equation for reduce to Equation (12) in the relevant case, and only with a factor of 1∕4 do we obtain optimal wave dispersion in the tests shown in Figure 3. Thus, the factor 1∕4, as in Equation (7), is indeed what is required.  (5) and (6)), and (c) modified flux form (Equations (7) and (8)). In each panel, the three branches represent acoustic modes (highest frequency), inertia-gravity modes, and Rossby modes (lowest frequency), all westward-propagating; there are also eastward-propagating acoustic and gravity-mode branches (not shown). The parameters used are as follows: domain depth 10 4 m; number of vertical levels 20; gravitational acceleration 9.80616 m ⋅ s −2 ; gas constant for dry air 287.05 J ⋅ kg −1 ⋅ K −1 ; specific heat capacity at constant pressure 1,005 J ⋅ kg −1 ⋅ K −1 ; Coriolis parameter 1.031 × 10 −4 s 1 ; northward gradient of Coriolis parameter 1.619 × 10 −11 s −1 ⋅ m −1 ; background temperature 250 K; surface pressure 10 5 Pa; zonal wavenumber 2 × 10 6 m −1 ; meridional wavenumber zero In this argument, it has been assumed that a vertical mass flux divergence in any -cell is accompanied by a compensating horizontal mass flux convergence, so that the net divergence of mass flux is small. This assumption is a good approximation for gravity waves and Rossby waves that are adversely affected by vertical averaging of and w. The assumption would not be a good approximation for acoustic waves, but acoustic waves are not significantly affected by vertical averaging of or w in the equation. The validity of the assumption is confirmed by the numerical normal mode analysis discussed below.
It is useful to examine the implied advective form equation for obtained by taking Equation (4) minus i +1∕2 × Equation (2): This shows that if the fluxes have a 2Δz vertical structure and has a background stratification then, provided that stratification is captured by the estimates for ∕ z, there will be a nonzero tendency of arising through the horizontal flux terms rather than the vertical flux terms. Also, if is independent of x and the mass fluxes are approximately nondivergent (i.e., the time derivative in Equation (1) is small) then Equation (11) reduces to the expected form: Finally, Equation (11) confirms that if is uniform, and provided̂takes that same uniform value of and the estimates for ∕ z vanish, then the tendency vanishes and the uniform value of is preserved.
To investigate whether this modified flux-form conservation equation for based on Equations (7) and (8) gives the accurate wave propagation expected for the Charney-Phillips grid, a numerical linear normal mode analysis was carried out, following the methodology of Thuburn and Woollings (2005) and Thuburn (2006). The compressible Euler equations on a -plane, with , , u, v, and w as prognostic variables, were linearized about an isothermal state of rest. The system was discretized in the vertical on a Charney-Phillips grid, and solutions proportional to exp{i(kx + ly − t)} were sought. For given values of k and l, the system comprises an eigenvalue problem for the normal-mode frequencies and their vertical structures. In this linear calculation, the values of̂are given by the reference profile about which we linearize; thus, the details of the advection scheme do not matter, except for the inclusion of the modification Equation (7).
Some example results are shown in Figure 3 for three versions of the discrete equation: (a) advective form; (b) naive flux form (Equations (5) and (6)), and (c) modified flux form (Equations (7) and (8)). For the advective-form equation, the numerical frequencies agree very well with the exact frequencies for the continuous linearized equations; this result is identical to that shown by Thuburn (2006) for the same system. For the naive flux-form equation the higher internal Rossby modes are significantly retarded due to the explicit and implicit vertical averaging of w. The modified flux form, however, captures all wave modes with the same optimal accuracy as the advective form equation.
Although it is the higher internal Rossby modes that are adversely affected by the averaging in the discrete equations in this example, when the horizontal wavelength is much shorter it is the higher internal gravity modes that are affected, as discussed by Thuburn (2006). The testing in Section 4 below focuses on gravity waves.

CONSERVATIVE SEMI-LAGRANGIAN TRANSPORT
Section 2 discusses how the fluxes should be discretized in a flux-form conservation equation for entropy so as to avoid losing the optimal wave propagation characteristics of the Charney-Phillips grid. Conservation equations can also be discretized in terms of remapping operators. 2 Such remapping-based schemes are attractive because they can be designed to be stable while remaining accurate even for large time steps. This section discusses how a such a remapping scheme can be modified so as to retain optimal wave propagation, using the SLICE transport scheme (Zerroukat et al., 2007;2009) for illustration.
Given a set of a set of trajectory departure points for the velocity points at the faces of the -cells, SLICE constructs the corresponding -cell departure volumes. It then effectively carries out a multidimensional remapping, via a "cascade" of one-dimensional remappings, to determine the density in the departure volumes. The mass in each departure volume is then assumed to be transported during the model time step to the corresponding arrival -cell. Figure 4a,b illustrates the idea in two dimensions. Using information about the grid and the departure points, SLICE constructs the departure volume (bounded by dashed lines in Figure 4b) for each -cell of the grid (bounded by solid lines in Figure 4b). Certain intermediate Eulerian control volumes are also constructed; their lateral boundaries are represented by vertical dashed lines in Figure 4a. In the first stage of SLICE, the density field is remapped in the x-direction, using x as the remapping coordinate, from the -cells to the intermediate Eulerian  Thuburn et al. (2010) showed that the effects of flow divergence could be captured more accurately by estimating the area of departure cells (or volume in three dimensions) using the trajectory average divergence, then using cumulative column area as the coordinate for the final remapping stage. This modification is used here too. (The freedom to choose different remapping coordinates in SLICE is discussed briefly in Appendix B.) Finally, the mass in each departure volume is assumed to be transported during the time step to the corresponding arrival -cell.
An obvious way to obtain conservative transport of entropy on a Charney-Phillips grid using SLICE would be to construct the departure volumes, intermediate Eulerian control volumes, and intermediate Lagrangian control volumes corresponding to -cells, and then to apply the above algorithm to the quantity (Figure 4c,d). However, if we use x and column integrated area as the remapping coordinates in the first and second stages, as we do for density, then the result of transporting with ≡ 1 initially is different from the result of transporting followed by averaging to w-levels. Thus will no longer be identically equal to 1 at the end of the time step; the property of preserving a constant is lost.
The property of preserving a constant can be recovered by again taking advantage of the freedom to choose alternative remapping coordinates. In this case, we must In the second stage of SLICE, is remapped quasivertically from intermediate Lagrangian control volumes to departure volumes. Finally, the entropy content in each departure volume is assumed to be transported during the time step to the corresponding arrival -cell use row-integrated mass as the remapping coordinate in the first stage and column-integrated mass as the remapping coordinate in the second stage. The density in -cells is simply ; the mass in -cell intermediate control volumes is just the vertical average of the masses in the -cell intermediate control volumes immediately above and below, and the mass in the -cell departure volumes is just the vertical average of the masses in the -cell departure volumes immediately above and below. From these densities and masses and the cell geometries, the required integrated mass coordinates can be constructed straightforwardly.
The resulting scheme is conservative and preserves a constant . Unfortunately, it suffers from essentially the same problem as the naive finite-volume scheme discussed in Section 2: the positions of the lateral boundaries of the -cell intermediate Eulerian control volumes and the upper and lower boundaries of the -cell departure volumes are determined by vertically averaged velocities, so a 2Δz pattern in the velocity is invisible to the transport; thus the scheme does not retain the Charney-Phillips grid optimal wave propagation. We refer to this as the naive SLICE scheme in the gravity-wave test of Section 4.
In order to obtain optimal wave propagation, we can include a correction to the horizontal remapping of that is analogous to the correction to the horizontal fluxes in Equation (7). The light shaded rectangular region in Figure 4c represents an -cell intermediate Eulerian control volume. In the first stage of the naive SLICE scheme, is remapped horizontally from -cells to these intermediate Eulerian control volumes. To apply the correction, observe that the difference of mass fluxes in the correction term in Equation (7), integrated over a time step, corresponds to twice the sum of the two dark shaded areas in Figure 4c (with appropriate allowance for sign) multiplied bŷ, an estimate of the density at that edge of the intermediate Eulerian control volume.̂may be estimated as a by-product of the horizontal remapping of . This difference of mass fluxes is multiplied by an estimate for ∕ z × Δz to obtain an entropy flux correction. An estimate for ∕ z × Δz at the edge of the intermediate Eulerian control volume can be estimated as half the difference between edge values of at the levels above and below; these edge values of again can be obtained as by-products of the horizontal remapping. The entropy flux correction is then applied to shift entropy conservatively between neighbouring intermediate Eulerian control volumes.

NUMERICAL EXAMPLES
Some numerical tests were carried out to confirm the conservation and wave dispersion properties of the modified SLICE scheme. The two-dimensional moist compressible Euler equations were solved using the semi-implicit, semi-Lagrangian vertical slice model described by Thuburn (2017a). The model uses a Charney-Phillips vertical grid, with specific humidity as well as entropy stored at w-points. The original formulation uses (in its default configuration) SLICE with parabolic spline subgrid reconstruction for conservative transport of density and a semi-Lagrangian scheme with cubic Lagrange interpolation for transport of entropy, specific humidity, and w. An option to use the modified SLICE scheme for transport of entropy, specific humidity, and w was implemented and compared.

Conservation
To test the conservation properties, the semi-Lagrangian and modified-SLICE versions were compared using the saturated buoyant bubble test case of Bryan and Fritsch (2002). The domain is 20 km wide and 10 km deep, discretized using 192 × 96 grid cells, giving a horizontal and vertical grid length of a little over 100 m. The time step is 10 s. Figure 5a,b compares the perturbation to equivalent potential temperature e for the two model versions after 800 s. There are some small but noticeable differences at the leading edge of the bubble, where the perturbation is slightly smaller for the modified-SLICE scheme. The overall evolution, however, is very similar for both versions. At later times (e.g., Figure 5c,d) the differences between the two versions grow. This test case, particularly the behaviour at the leading edge of the bubble, is notoriously sensitive to the details of the numerics (e.g., Duarte et al., 2014;Kurowski et al., 2014), and a similar sensitivity was found here. All of the schemes tested-semi-Lagrangian, modified SLICE, and unmodified SLICE, each with and without limiters-produced clear differences from each other at the bubble's leading edge at 1,000 s. Figure 5e,f compares the conservation of entropy (solid lines) and energy (dashed lines) in the semi-Lagrangian and modified-SLICE versions. Note the different axis scale in the two panels. See Appendix C for a discussion of how the entropy and energy changes are normalized. In the semi-Lagrangian version, there are normalized entropy and energy losses of around −0.15 over 1,000 s, while in the modified-SLICE version entropy is conserved to machine precision. The model is not formulated to conserve energy exactly, so we should not expect perfect energy conservation even with the modified-SLICE version. Figure 5f shows a very slight increase in energy between 200 and 500 s, before numerical dissipation of kinetic energy becomes significant at later times. Nevertheless, as a by-product of the entropy conservation, the energy loss in the modified-SLICE version is an order of magnitude smaller than that in the semi-Lagrangian version.
The total water content (not shown) is actually conserved to machine precision for both model versions in this test case. In the semi-Lagrangian version, this exact conservation occurs because the specific humidity is uniform. Preservation of this uniform specific humidity by semi-Lagrangian advection, together with conservation of mass by SLICE transport, implies conservation of total water. For nonuniform specific humidity, the semi-Lagrangian version would not conserve total water. The uniform specific humidity was also preserved to machine precision by the modified-SLICE version, confirming that the mass-coordinate-based remapping for entropy and specific humidity works as intended.

Wave dispersion
To test the wave dispersion properties, the model was initialized with a packet of gravity waves of small vertical wavelength, and the frequency of the waves in the model was compared with the analytical frequency and the theoretical optimal numerical frequency. As above, the domain size was 20 km × 10 km with resolution 192 × 96 grid cells. A background resting hydrostatically balanced state with surface pressure 10 5 Pa and uniform temperature T = 270 K was set up. Specific humidity was set to zero. A gravity-wave packet disturbance was then superposed, with the following distributions of buoyancy b and mass stream function : Here where (x c , z c ) = (10 4 m, 5 × 10 3 m) is the centre of the wave packet, with x r = 7 × 10 3 m and z r = 3.5 × 10 3 m defining the size of the wave packet, w 0 = 0.01 m ⋅ s −1 is the vertical velocity amplitude, k and m are the horizontal and vertical wavenumbers of the wave, N 2 = g 2 ∕c p T is the background buoyancy frequency squared, with g = 9.81 m ⋅ s −2 the gravitational acceleration and c p = 1,004 J ⋅ kg −1 ⋅ K −1 the specific heat capacity at constant pressure, and c ≈ 0.7 kg ⋅ m −3 is the background density at the centre of the wave packet. The frequency for a monochromatic gravity wave is given, to a good approximation, by the Boussinesq frequency, The initial density and entropy are adjusted to give the buoyancy field specified by Equation (13) without perturbing the pressure. The mass stream function is used to construct the initial velocity field: The resulting disturbance evolves as a packet of nearly monochromatic gravity waves, with phase propagation towards the lower left perpendicular to the phase lines and group propagation towards the upper left, approximately parallel to the phase lines, in agreement with the theory of idealized gravity waves (e.g., section 7.3 of Vallis, 2017). For well-resolved waves, the numerical frequency should be close to the analytical Boussinesq frequency (Equation (17)). However, for waves that are less well-resolved in space, even with optimal grid staggering, the inexact approximation of derivatives changes the frequency. For second-order centred-difference derivatives on a Charney-Phillips C-grid, as used here, the effect of the numerical errors can be quantified; the effect and is to replace the exact frequency (Equation 17) by the (optimal) numerical frequency, wherê are the effective wavenumbers seen by staggered centred-difference derivatives. On a Lorenz grid, or on a Charney-Phillips grid with naive flux-form transport of entropy or with a suboptimal form of the pressure-gradient term, the vertical averaging of w and/or introduces a further factor cos 2 (mΔz∕2) in Equation (19), severely slowing those waves that are marginally resolved in the vertical (Thuburn, 2006). For the experiments discussed here, the time step Δt = 10 s is much shorter than the wave period, so time discretization errors are negligible. The horizontal wavenumber was fixed at k = × 10 −3 m −1 . For a range of vertical wavenumbers, the TA B L E 1 Exact, optimal numerical, and empirical gravity wave periods for a range of vertical wavenumbers m gravity-wave packet was simulated and the empirical period of the wave was estimated from time series of u, w, and at the centre of the domain. (This estimate incurs some small errors, because attention is restricted to periods that are multiples of the time step, and because the wave packet propagates away, so that wave amplitude at the domain centre decays over time.) The empirical period was then compared with the theoretical Boussinesq and optimal numerical periods. The results are shown in Table 1.
For both semi-Lagrangian and modified SLICE transport of , the wave periods agree well with the theoretical wave periods for an optimal scheme. The same set of gravity-wave simulations was carried out with transport given by the naive SLICE scheme discussed in Section 3. The results are shown in the final column of Table 1. They confirm that this naive application of SLICE does not lead to optimal wave propagation. They also confirm that this test case can indeed discriminate between optimal and suboptimal wave propagation. In fact, these periods are longer than the optimal periods by the theoretical factor 1∕ cos(mΔz∕2).

SUMMARY AND DISCUSSION
On a Charney-Phillips vertical grid, for which entropy is staggered vertically relative to density , conservation of entropy can be obtained by integrating a flux-form conservation equation for the quantity . A naive discretization of that conservation equation involves explicit and implicit vertical averaging, so that optimal wave propagation is lost, despite the use of a Charney-Phillips grid. This article presents a straightforward and general method for modifying the horizontal fluxes in the entropy conservation equation so as to restore optimal wave propagation. An analogous modification can be made in a conservative semi-Lagrangian scheme based on remapping.
The proposed idea has been tested, and the predicted behaviour confirmed, by computing numerical linear normal modes for an idealized basic state, and by simulating a saturated buoyant bubble and marginally resolved gravity waves in a two-dimensional vertical slice model.
For clarity of presentation, the idea has been presented in the two-dimensional context and for a vertically uniform grid. However, it extends straightforwardly to three dimensions and to vertically nonuniform grids.
At an early stage of this work, an alternative modification of the horizontal fluxes was considered: In this scheme, is first vertically averaged to -levels, then used to construct values at the lateral faces of -cellŝ i+1∕2 and hence entropy fluxes at the lateral faces of -cells F x i+1∕2̂i +1∕2 , which are then averaged back to w-levels to give G x i+1∕2 +1∕2 . Interestingly, this scheme can be obtained by a slight modification of the derivation in Equations (9) and (10) using a different piecewise constant subgrid reconstruction of̂i +1∕2 (z): ( Figure 2d). This scheme is conservative and has the optimal wave dispersion property. However, it is less accurate than Equation (7). This is most clear if we consider the case of F x independent of z, whereupon The double vertical average of̂is only a second-order approximation tôi +1∕2 +1∕2 , so the scheme is, at best, second-order accurate for advection. This reduction in accuracy is noticeable in advection tests with the SLICE analogue of the scheme. On the other hand, the modification described in Section 3, which is the SLICE analogue of Equation (7), remains as accurate as the unmodified SLICE scheme. The modification described in this article permits considerable flexibility in the choice of the cell-edge valueŝ, or in the choice of subgrid reconstruction in the case of the SLICE scheme. In particular, schemes with a high order of accuracy are possible, and so are flux limiters that prevent the numerical generation of overshoots and undershoots in . The results shown in Section 4 all use a parabolic spline subgrid reconstruction for SLICE and include a limiter for the transport of entropy and water (Zerroukat et al., 2006). The results using semi-Lagrangian advection shown in Table 1 use two-dimensional cubic Lagrange interpolation with a simple monotonicity limiter. Switching off the limiters makes negligible difference to the results in Table 1.
In summary, the modified SLICE scheme presented here, applied to the transport of entropy on a Charney-Phillips vertical grid, achieves several desirable properties: high overall accuracy, stability at large time steps, conservation, preservation of a constant, and prevention of overshoots and undershoots. It does all this without sacrificing the optimal wave propagation permitted by the Charney-Phillips grid.
An analogous choice for remapping density in the second stage of SLICE would be s = z and q = x, where x(z) is the width of the column being remapped. However, as discussed by Thuburn et al. (2010), the effects of flow divergence can be captured more accurately by using a volume-based coordinate (in two dimensions an area-based coordinate) that incorporates the column width together with q = . This area-based coordinate is used to remap density in the second stage of SLICE (Figure 4b).
To ensure consistency with the transport of mass, and hence ensure preservation of a constant , the transport of uses column-integrated mass as the coordinate; for the first stage, and for the second stage, together with q = . This choice ensures that, when ≡ 1, the entropy content in a remapped cell agrees with the mass content: Q =s +1∕2 −s −1∕2 = ∫ cell dA.
It is useful to note that the first stage of the density remapping may be re-interpreted as using row-integrated area as the remapping coordinate: These choices of remapping coordinate then suggest the following general rule of thumb.
• For remapping density, in order to capture the effects of velocity divergence accurately, use integrated area (integrated volume in three dimensions) as the remapping coordinate in all stages of SLICE. The same applies when remapping the velocity divergence itself in order to compute the departure cell areas (or volumes) (Thuburn et al., 2010).
• For remapping tracer density, for example , in order to ensure compatibility with the density remapping and ensure preservation of constant , use integrated mass as the remapping coordinate in all stages of SLICE.

APPENDIX C. NORMALIZATION OF ENTROPY AND ENERGY CHANGES
A natural way to normalize the entropy and energy changes in the buoyant bubble test might appear to be to compute the fractional changes in these quantities. However, arbitrary constants may be added to the definitions of specific internal energy, potential energy, and specific entropy without changing any of the essential physics (e.g., Feistel et al., 2008). Since the diagnosed fractional change in entropy and energy will depend on the choices for these constants, the fractional change is not a unique and objective measure. Instead, we normalize the energy change by KE max ≈ 1.3 × 10 9 J, the maximum domain-integrated kinetic energy during the 1,000 s run, and the entropy change by KE max ∕T max , where T max ≈ 289.6 K is the temperature near the surface.