A solution to the trilemma of the moist Charney–Phillips staggering

The Charney–Phillips grid, used in many numerical models of the atmosphere, involves vertically staggering the nodes of the density variable with the nodes of the entropy‐type variable. When moisture is included in such a model, it is either co‐located with density so that moisture can be transported conservatively and consistently with dry mass, or with the entropy‐type variable so that the coupling between moisture and temperature can be represented well. Both properties are desirable, yet at first it appears difficult to obtain both simultaneously. Here, we present a framework to resolve this problem, by co‐locating the moisture mixing ratio with potential temperature but formulating its transport as that of a density on a vertically shifted mesh. Within this framework, particular choices of the operators involved provide the desired conservation and consistency properties of the moisture transport. The framework is described in the context of a finite‐element approach. We also present an explicit Runge–Kutta time‐stepping scheme that is appropriate for use within this framework. This approach is then illustrated through numerical tests, which demonstrate that it does indeed have the desired conservation and consistency properties.


The trilemma
For numerical models of the dry atmosphere, there are two common choices for the vertical placement of the density/pressure variable and the buoyancy/entropy variable.
Here, these prognostic variables are taken to be the dry density d and the (dry) potential temperature 1 . The first 1 The potential temperature is defined by ∶= T(p 0 ∕p) R d ∕c p , where T is the temperature, p is the air pressure, p 0 = 1 × 10 5 Pa is a reference choice of vertical placement is the Lorenz grid, in which d and are co-located (Lorenz, 1960). This can be exploited to facilitate the global and local conservation of entropy and energy by writing the transport of as On the other hand, the Charney-Phillips grid staggers d and vertically, by co-locating with w, the vertical pressure, R d is the gas constant of dry air, and c p is the specific heat capacity of dry air at constant pressure. component of the velocity (Charney and Phillips, 1953). A number of studies, such as Fox-Rabinovitz (1994) and Thuburn and Woollings (2005), have shown that the Charney-Phillips grid gives a better discrete representation of gravity waves than the Lorenz grid does. However, with the Charney-Phillips grid it is not easy to locally conserve entropy, creating a dilemma. This was the motivation for Thuburn (2022), who presented a conservative transport scheme for entropy in the Charney-Phillips grid without losing good wave dispersion. When moisture is added to a Charney-Phillips grid, the dilemma becomes a trilemma. Here, species X of moisture is represented through the mixing ratio m X , given by where X is the mass density of species X. Many models choose to co-locate m X with d , again to attain conservation of moisture, which is considered an important aspect of climate and numerical weather prediction models (Bush et al., 2020). However, as argued by Konor and Arakawa (2000), co-locating m X with gives a better representation of the coupling between moisture and temperature, through the calculation of the saturation vapour pressure, the latent heat that is released or absorbed as water changes phase, and the representation of inversions of temperature. It is this that has motivated the Met Office to co-locate m X and with a Charney-Phillips staggering in its Unified Model (Wood et al., 2014;Walters et al., 2017), but at the cost of losing local conservation of moist mass. It seems that any choice of staggering excludes at least one desirable property. Although any moist Charney-Phillips staggering will have this trilemma, this work is particularly motivated to find a solution for the Met Office's new dynamical core, GungHo (Melvin et al., 2019). The current ENDGame dynamical core (Wood et al., 2014) uses a latitude-longitude grid that suffers from a scalability bottleneck associated with the Poles (Lawrence et al., 2018). To avoid the scalability bottleneck, GungHo uses a cubed-sphere mesh, which is quasi-uniform across the globe. This is combined with a mixed finite-element discretisation, in which the prognostic variables are chosen from a suitable family of finite-element spaces. It was shown by Cotter and Shipton (2012) and Cotter and Thuburn (2014) that such a model can replicate the desirable properties described by Staniforth and Thuburn (2012) that are provided by ENDGame's Arakawa C-grid staggering on an orthogonal latitude-longitude mesh.
The lowest order finite elements used by GungHo on a mesh of hexahedral cells can be thought of as the finite-element analogy of the variable staggerings used by ENDGame. The density d can be described by its value in the centre of cells, and the wind v can be described through its value in the centre of the cells' facets, thus staggering it from d in each direction. It was shown by Melvin et al., (2018) that the properties of the Charney-Phillips grid could be retained in this set-up by co-locating with the vertical component of v. Throughout this work, the finite-element function spaces used by d , v, and are denoted by V , V v , and V respectively. Although the language of finite elements is used in the remainder of this work, the formulation of Sections 2-5 could also be applied to finite-difference or finite-volume discretisations, in which case V , V v , and V correspond to the appropriate representations of , v, and in those discretisations.
In the context of GungHo, the natural choices of space for the moisture mixing ratio would be V or V . To replicate the arrangement used by ENDGame, this work takes m X ∈ V , which was also the choice made by Bendall et al., (2020) in a similar finite-element discretisation of the moist atmosphere. The trilemma challenge that this article addresses is then how to achieve conservation of total water content in a manner consistent with that for dry density held in V .

Moisture transport
In the absence of sources and sinks of the moisture species, the density X evolves according to the conservative form of the transport equation: Integrating Equation (3) over the volume of the domain shows that, in the absence of inflow and outflow, the mass of this moisture species is conserved. As the dry density also obeys the conservative transport equation, the mixing ratio obeys the advective form of the transport equation: We are seeking a transport scheme that is locally conservative, discretising some form of Equation (3). At the same time the transport of X must also be consistent with the transport of d . Here, as in Lauritzen et al., (2011), consistency is defined to be that a spatially constant m X field always remains constant for all transporting winds. A corollary of this is that if X is initially equal to d then the two fields will evolve in the same way.
Other studies, notably Lauritzen et al. (2011, section 8.7; and Zängl et al., (2015), describe the use of the dry flux to obtain consistent and conservative transport of tracers on the Lorenz mesh (this method is also sometimes called "free-stream preserving"). As shown by Thuburn (2022), a transport scheme can achieve these properties on a Charney-Phillips staggering by introducing a vertically shifted mesh, to which the dry fluxes are mapped. The key is that when the dry fluxes are mapped to the vertically shifted mesh the dry mass budget should be preserved, analogously to the method used by Ringler et al., (2010) to obtain consistent evolution of potential vorticity, mass, and velocity. This work builds upon the approach of Thuburn (2022) by presenting a general framework for the consistent and conservative transport of moisture to resolve the trilemma.
The remainder of this article is laid out as follows. Section 2 describes a vertically shifted mesh similar to that used by Thuburn (2022) to facilitate conservative transport of variables co-located with . The general framework for consistent and conservative transport using this shifted mesh is then presented in Section 3, which places some requirements on the operators for mapping variables to the shifted mesh and also on the scheme used for computing the moisture fluxes. Specific operators satisfying these requirements are presented in Section 4, and a flux scheme using Runge-Kutta time stepping is presented in Section 5. The framework is then demonstrated through some idealised tests in Section 6, before a summary in Section 7.

A VERTICALLY SHIFTED MESH FOR CONSERVATIVE TRANSPORT
To achieve local conservation of mass of the moisture species, the conservative form of transport equation, Equation (3), is discretised, while still using the mixing ratio m X as the prognostic variable for the physical parametrisations. Key to this is the discrete mapping between X and m X , involving d via Equation (2). As discussed in the previous section, m X is co-located with for accurate representation of latent heating and phase changes.
Where should X be located within the Charney-Phillips setup? Many conservative transport schemes for d rely on computing mass fluxes for a control volume surrounding each d node. In a finite-element context, these nodes are known as degrees of freedom (DoFs). In the Charney-Phillips grid, this control volume is naturally bounded in the vertical by the w-levels, which are co-located with , and the mass fluxes are computed at these points and are therefore staggered relative to d . In a finite-element model this is analogous to a discrete relationship between V v and V , such that for all v ∈ V v the divergence ⋅ v ∈ V . This property is assumed throughout. In contrast, the DoFs of V are co-located with the DoFs of V v that represent the vertical component of the velocity. This means that, in the standard Charney-Phillips staggering, illustrated in Figure 1, V has one more DoF per column than V does. This presents a difficulty in choosing the space of X . If X were to be computed in V , there would not be a one-to-one relationship between m X and X and information would be lost at each transport step. Additionally, there is not such a neat relationship between ⋅ v and V as there is between ⋅ v and V , making design of a conservative transport scheme in V more challenging.
Instead, we follow Thuburn (2022) in introducing a second vertical mesh, whose levels are vertically shifted relative to the primary mesh. The shifted mesh w-levels are located between the w-levels on the primary mesh. The heights of the top and bottom w-levels of the shifted mesh still match those of the primary mesh, meaning that the layers at the top and bottom of the model are of smaller depths. This same shifted mesh can be used in a lowest order finite-element model like that of GungHo, where is represented by a continuous, piecewise linear function in the vertical. This is equivalent to a finite-difference or finite-volume method, with DoFs lying in the centres of the top and bottom faces of a cell. In this article, quantities on the shifted mesh are denoted with a tilde ∼ ⋅. The levels of the two meshes can then be summarised as follows. Let d have N DoFs in a column on the primary mesh, and let have N + 1 DoFs. In a lowest order finite-element model, the facets at the top/bottom of cells coincide with the w-DoFs, which have vertical coordinates { w k } for k ∈ [1, N + 1]. The vertical coordinates of the cells' tops/bottoms on the shifted mesh are {̃w k } for k ∈ [1, N + 2], which are related to those of the primary mesh bỹ where the A k ∈ (0, 1) describe the displacement of the levels between the two meshes. For the lowest order finite-element model, a density on the shifted mesh takes a constant value in each of the cells, but the corresponding finite-volume model has -levels on the shifted mesh with vertical coordinates of {̃k }. There is flexibility in how {̃k } can be chosen, though the equivalent values to thẽ V DoF locations in the lowest-order finite element model used in Section 6 arẽk = 1 2 (̃w k +̃w k+1 ) for k ∈ [1, N + 1]. These staggerings are illustrated in Figure 1. It can be seen here thatṼ has the same number of degrees of freedom as V , although these two spaces are not the same. F I G U R E 1 An illustration of the vertically shifted mesh and the Charney-Phillips staggering. In the Charney-Phillips staggering, the dry density d is staggered vertically relative to the potential temperature . The top/bottom facets of cells for a lowest order finite-element model are shown by solid black lines, with the midpoints shown by grey dashed lines. The moisture mixing ratio m X and vertical velocity w are co-located with , and the horizontal velocity u is on the same level as d . The dry mass flux F d is co-located with the velocity (the horizontal part with u and the vertical part with w) on the primary mesh, while the horizontal and vertical parts of the fluxesF d andF X are respectively represented byũ andw on the shifted mesh. The vertical coordinates of the interfaces between cells are given by w . The vertically shifted mesh has one more layer than the primary mesh, with the levels shifted relative to those of the primary mesh (although the top and bottom of the two meshes are co-located). In the conservative transport formulation of Section 2, the moisture density X is described on the vertically shifted mesh and is transported by a shifted velocity Having introduced this mesh, the moist densitỹX is computed inṼ from d ∈ V and m x ∈ V through an operator  ∶ V , V →Ṽ . The transporting velocity or mass flux can also be expressed in the spaceṼ v on the shifted mesh, through an operator The basic premise for computing moisture conservatively is to computẽX at the start of the transport step from the respective values of d and m X , then to transport X conservatively, before evaluating m X ∈ V by inverting  with respect to m X , using the value of d at the end of the transport step. This operator is called  −1 ∶Ṽ , V → V . Using superscript n to denote values before a transport step, and n + 1 to denote values after the transport step, this is summarised as̃n There are two requirements about : Requirement 1. The operator  ∶ V , V →Ṽ is linear in both its arguments, so that for any constants , , A, B, for any m (1) Requirement 2. The operators  ∶ V , V →Ṽ and  −1 ∶Ṽ , V → V are inverses when their second arguments use the same dry density d , such that for all m X ∈ V , for all positive, non-zero d ∈ V and all̃X ∈Ṽ : Requirement 2 ensures that the operators  and  −1 are reversible so that, in the absence of transport, m X is unchanged. Likewise, if there are no changes to m X and d in between transport steps theñX at the start of a transport step is equal to its value at the end of the previous step. The reversibility of  for a specific d is possible becausẽ V and V have the same number of DoFs.
Finally, in this approach, the mass that is conserved is defined using̃X :

A CONSISTENT AND CONSERVATIVE TRANSPORT FRAMEWORK
The previous section describes the transformation of a mixing ratio m X to a densitỹX on a vertically-shifted mesh, so that moisture can be transported conservatively. This section presents a framework to also transport moisture consistently with dry air, by ensuring that a constant mixing ratio is preserved.
Defining the dry mass flux as F d ∶= d v and the moisture flux as F X ∶= X v, the continuous transport equations for d and X can be written as Using Equation (2), the moisture flux can also be expressed as F X = m X F d , so the transport equation for X becomes The crucial ingredient for transporting X consistently with d is to transport X via Equation (11) and to use the same flux F d in the transport of both densities. This is the approach discussed by Lauritzen et al. (2011Lauritzen et al. ( , 2014, Zängl et al., (2015), and Thuburn (2022).
To discretise Equation (10), we consider conservative, explicit, two-time-step transport schemes for d ∈ V and̃X ∈Ṽ that can be generally expressed via the dis- where  ∶ V v → V and ∶Ṽ v →Ṽ are the discrete divergence operators on the primary and shifted meshes respectively and, as in Section 2, the superscript n denotes a field at the nth time level. The fluxes F d andF X are computed from flux operators As in Thuburn (2022), this requires F d to be evaluated on the shifted mesh, which is denoted byF d ∈Ṽ v . This shifted mass flux is obtained through the operator To be as general as possible here, we avoid specifying any details about the discrete divergence operators, or of the spatial or temporal discretisations involved in  d and X .
Although  d and X could use the same temporal and spatial discretisations (albeit on different meshes) they do not necessarily have to do so. Now consider a mixing ratio that is spatially constant, so that m n X = C. The transport of moisture is consistent if m n+1 X = C, irrespective of the transporting wind v and the dry density d . Since m n+1 d ], this requires the process of representing d on the shifted mesh to commute with transporting it. This is illustrated by the commutativity diagram in Figure 2.
For this consistency to be achieved, in addition to Requirements 1 and 2, there are two further requirements about the operators involved in the framework: Requirement 4. For any constant C ∈ V and for all Requirement 3 is represented by the commutativity diagram in Figure 3. That these requirements lead to a consistent and conservative transport scheme can be shown by substituting Equation (12) Using Equations (6a) and (13) gives ] .
(17c) F I G U R E 2 A diagram representing the commutativity that is required by the framework presented for a transport scheme to be consistent. In the case that the mixing ratio is a constant C, transport should commute with the transformation of a density in V to a density inṼ by the operator . The algorithmic implementation of the transport of̃X is described by the left and upper arrows, which must (in combination) be equivalent to the lower and right arrows to provide consistency F I G U R E 3 A representation of the required commutativity for the operators  and  that respectively map fluxes and their divergences to the shifted mesh in the presence of a constant mixing ratio of value 1. This is described in Requirement 3 F I G U R E 4 A schematic illustrating the different variables, function spaces, and operators involved in the consistent and conservative transport framework presented in Section 3. Each row of the figure corresponds to a particular function space, with fields on the row matching their function space. Each arrow represents the action of an operator, beginning at the fields that are the arguments of that operator and ending at the field that is the result of the operator [Colour figure can be viewed at wileyonlinelibrary.com] Inserting m n X = C, this becomes m n+1 Then using the four requirements about the scheme and Equation (14), , (17e) and the transport of m X is both consistent and conservative.
To summarise, in this framework, conservation is achieved by converting m X ∈ V tõX ∈Ṽ , which is then transported. The densitỹX lies on a vertically shifted mesh to facilitate both the flux form transport scheme and the one-to-one relationship between m X and̃X . The consistency of the scheme comes from using F d to evaluate the moisture mass flux. The transported mixing ratio is then calculated from the transported values of d and̃X . The action of the different operators used in this framework is summarised by the schematic in Figure 4.

CONSISTENT MAPPING OPERATORS
This section presents specific choices for the operators that map fields from the primary mesh to the shifted mesh, satisfying Requirements 1-3. A flux operator satisfying Requirement 4 will be considered in Section 5.
First, consider an operator [m X , d ] of the form so that the value returned by  is the product of the values returned by two operators  ∶ V →Ṽ and  ∶ V → V , each of which directly maps a variable toṼ . The multiplication is pointwise at each DoF ofṼ . If  and  are linear in their own arguments then  is linear in both arguments, satisfying Requirement 1. For [m X ], we simply identify the values of m X at the DoFs ofṼ . In the top and bottom half-layers this involves linear interpolation to theṼ DoFs, whereas in the interior of the domain this is the identity operator, so that the kth value in a column is given bỹ using the vertical coordinates of the levels from Equation (5). Note that [C] = C for a constant field C. As V andṼ have the same number of DoFs, this choice of  can easily be inverted to give the operator  −1 ∶Ṽ → V . For the top and bottom values of a column,  −1 extrapolates instead of interpolates. Furthermore, by using the identity operator on the interior of the domain the generation of new maxima or minima in m X is avoided, which is generally an important property of schemes for transporting moisture.
Taking  of the form in Equation (18), the operator  −1 that satisfies Requirement 2 is given by From the definition of  it follows that if m X is a constant C, [C] = C, and so Then, Equation (15) of Requirement 3 reduces tõ and thus  and  must be chosen together to ensure the consistency of the scheme.

A specific choice for the operators
To find a form of the operators  and  satisfying Equation (22), it will be helpful to establish some more notation. Consider first an element in the primary mesh, such that the kth element in the ith column is indicated by e k i . This element has lateral surfaces that are denoted by s k i,l . Unless the element is on the exterior boundary of the domain, this surface will connect to a neighbouring element, e k i,l . Each column also contains N + 1 horizontal surfaces (associated with vertical fluxes), which form the top and bottom faces of each cell. These surfaces are denoted by s k i,b and s k i,t , with s k i,t = s k+1 i,b . All these entities are denoted with a tilde on the shifted mesh.
Before describing operators  and , it is convenient to now express the action of . If, for a flux F ∈ V v , [F] = Δ ∈ V , then let Δ k i be the value of Δ in e k i . The outward flux normal to s k i,l is F k i,l ; and similarly, F k i,t and F k i,b are the outward fluxes through the top and bottom surfaces, with F k i,t = −F k+1 i,b . Then, the action of  can be expressed as which is a discrete form of the divergence theorem. Note that F is taken here to have the same units as v.
We look for operators  and  that are linear in their arguments (and so independent of F): and k i,b are all coefficients to be found, with Equation (24a) representing , and Equations (24b) and (24c) together representing . To justify this, consider Figure 1. Each element on the shifted mesh overlaps with two from the primary mesh. Each lateral face on the shifted mesh overlaps with two from the primary mesh, and each "horizontal" face (with a vertical flux) lies between two horizontal faces on the primary mesh. The bottom layer of the shifted mesh only overlaps with a single layer of the primary mesh, so that 1 i = 0, 1 i,l = 0, and 1 i,b = 0. The same reasoning applies in the top layer of the shifted mesh, so that N+1 i = 0, N+1 i,l = 0, and N+2 i,b = 0. In Equation (24c), to simplify the notation we are taking Combining these aspects, the matrices representing these operators are of bidiagonal form within each column. It is stressed that these may not be the only forms of  and  that will satisfy Equation (22).
Assuming  and  are of the form in Equations (24a) and (24b) and  is of the form Equation (23b), substituting these into Equation (22) gives This can be rearranged to Since this must be true for any F, it is obtained that Further, mass must be preserved as density is mapped from V toṼ . This implies that and using that 1 i = 0 and N+1 i = 0, the two integrals can be incorporated into the same sum, so this becomes This will then be satisfied if which also ensures that the first term on the second line of Equation (26) also vanishes. The forms of the operators are then as follows: The presence of k−1 i before the second integral of Equation (30c) indicates that vertical fluxes are obtained on the shifted mesh from interpolation, in contrast to the method used to map the density and the horizontal fluxes.
How can the k i coefficients be obtained? It appears logical that they should relate to features of the mesh. However, since the same constant is used to map each lateral flux of a given cell, the constant cannot depend on the areas of the lateral surfaces. Similarly, a surface typically joins two elements together, which may in general not have the same volume, so the k i cannot relate to the volumes of elements. Since ∫ F dA is continuous between cells, it must be true that k i,l is equal to a corresponding value for the neighbouring column. As k i,l = k i , this means that k i must be independent of the column i. This leaves only that the coefficients can depend on k; and since the vertical coordinate is common between neighbouring cells, then if the vertical coordinates of the mesh levels are given by Equation (5)

A CONSISTENT RUNGE-KUTTA SCHEME
With the choices made in Section 4, all requirements except Requirement 4 have been met. In contrast to the other requirements, this depends on the flux operator  X , requiring that, if the mixing ratio is a constant field everywhere,  X should return a moisture flux F X that is simply the dry flux F d scaled by that constant. In this section, a "method of lines" scheme is introduced using explicit Runge-Kutta time stepping that obeys this requirement.

Standard Runge-Kutta schemes
First, consider a general N-stage explicit Runge-Kutta scheme for discretising the advective form of the transport equation for a mixing ratio m X ∈ V by a velocity v ∈ V v . Let the weights of the Runge-Kutta scheme be a ,k and the th stage of the scheme be denoted by m ( ) X . The operator  ∶ V , V v → V discretises (v ⋅ )m X , and it is required, for numerical consistency, that [C, v] = 0 for any constant field C and any v ∈ V v . Following, for instance, Shu and Osher (1988), for an explicit Runge-Kutta scheme the value of m X at the (n + 1)th time step, m n+1 X , can be written as To describe a comparable scheme for the conservative form of the transport equation for some density ∈ V , we introduce the operator  ∶ V , V v → V v , which acts as the operator for a single step in the Runge-Kutta scheme. 2 The operator  performs a flux reconstruction, creating a discrete form of v, and is assumed to satisfy [C, v] = Cv for any constant field C and any v ∈ V v . Then, the Runge-Kutta scheme for the conservative form of the transport equation is

5.2
Steps of the moisture transport scheme A Runge-Kutta scheme discretising Equation (11), but with a flux operator satisfying Requirement 4, can be made by combining Equations (31a-c) and (32a-c). In this scheme, the solution is built up by using the advective form of the transport equation of m X ∈ V for the first (N − 1) stages of the scheme, using the conservative form for̃X ∈ V only for the final stage. This is analogous to the flux form semi-Lagrangian schemes of Leonard et al., (1996) and Lin and Rood (1996), where, to ensure a scheme preserves a constant, the advective form is used in all stages except the final one. Take ∶Ṽ ,Ṽ v →Ṽ v as a flux operator on the shifted mesh. The overall scheme then consists of the following steps:

NUMERICAL RESULTS
This section demonstrates the consistent moisture transport framework of Section 3 through test cases in two-dimensional vertical slice domains, with coordinates (x, z). The shifted mesh is defined using Equation (5), taking A i = 1∕2 for all the levels. For the consistent transport of moisture, the operators of Section 4 and "method of lines" scheme of Section 5 are used. The Runge-Kutta weights are those of the three-stage strong stability-preserving scheme (SSPRK3); see Cockburn and Shu (2001) for more details. The advective operator  and the flux operator  use a third-order upwind finite-volume reconstruction. More information on this reconstruction can be found in Melvin et al., (2019). The dry density d is transported using the standard SSPRK3 scheme of Equation (32a-c). The new transport scheme is compared with an advective-form transport scheme for the moisture mixing ratio. This comparison scheme uses the standard SSPRK3 method in Equation (31a-c) with the same advective operator  as is used for the new scheme.

Idealised transport tests
To demonstrate the convergence and consistency of the framework, a mixing ratio field and a dry density field are transported subject to a prescribed deformational, divergent flow in a vertical slice of length L x and height H z , with x ∈ (−L x ∕2, L x ∕2) and z ∈ (0, H z ). Inspired by the deformational, divergent flow of Nair and Lauritzen (2010) and Lauritzen et al., (2012), the transporting velocity v = (u, w) is where L x = 2 km, H z = 2 km, and = 2000 s is the length of the simulation. The domain is periodic in the x-direction and the flow includes the contribution of uniform flow of speed U = L x ∕ in the x-direction, so taking x ′ = x − L x ∕2 − Ut the true solution at time is equal to the initial condition. The other speed is W, which can be varied relative to U to adjust the amount of deformation; but for the results presented here, W = U∕10. The flow is illustrated by the arrows in Figure 6. F I G U R E 5 The d field used in the convergence configuration of the transport test case in Section 6.1. Left: The initial condition. Right: A computed state at t = . Both are shown with a contour spacing of 0.1 kg⋅m −3 , which varies linearly in the vertical direction. Centre: A computed state of d at t = ∕2, shown with a contour spacing of 0.25 kg⋅m −3 . As the flow is divergent, the initially smooth field is deformed, before returning to close to its original state [Colour figure can be viewed at wileyonlinelibrary.com] This flow is used in two configurations, to separately show the convergence and consistency properties of the moisture transport. To describe the initial conditions, first define the distance from a point (x c , z c ) as (39) Then a Gaussian hill of width c centred at (x c , z c ) is given by For the convergence configuration, the initial conditions use two Gaussian hills for m X and a smooth d varying uniformly with height: with , and t = 0.5 kg⋅m −3 . These initial states are displayed in Figures 5 and 6, alongside numerical computations of the fields at t = ∕2.
To demonstrate the convergence of the scheme, the L 2 error norm between the initial m X field and its computed value at t = was measured for a range of resolutions. The horizontal and vertical grid spacings Δx and Δz were varied together, with uniform Δx = Δz for all simulations. The values used for Δx and Δz were those corresponding to subdividing the domain into 120, 140, 160, 180, and 200 cells in each direction. The time step Δt = 2 s was kept constant for all simulations. The convergence is displayed in the left-hand part of Figure 7, demonstrating that the consistent transport has the same order of accuracy as the advective transport scheme. However, as shown in the left-hand part of Figure 9, the consistent framework conserves mass to machine precision throughout the simulation whereas the advective transport does not.
The consistency configuration of the test case shows that the framework preserves a constant mixing ratio, even in the presence of a divergent, deformational flow. Now the dry density field is initially two Gaussian hills, so that the fields at t = 0 are with (x 1 , z 1 ) = (L x ∕8, H z ∕2), (x 2 , z 2 ) = (−L x ∕8, H z ∕2), c = 2L x ∕25, and m 0 = 0.02 kg⋅kg −1 as before, but now with f 0 = 0.5 kg⋅m −3 and b = 0.5 kg⋅m −3 . The initial d field and a computed state at t = ∕2 are shown in Figure 8. In the consistency configuration, there was just a single simulation using the flow (38), taking Δx = Δz = 20 m and Δt = 2 s. The right-hand side of Figure 7 plots the L 2 norms of d and m X as a function of time, which allows these two fields to be compared. Even though d is deformed by the flow, the mixing ratio field remains constant.

Moist rising bubble
The next test is to use the consistent and conservative moisture transport in the context of a non-hydrostatic dynamical core. This is the GungHo model of Melvin et al., (2019), which uses a semi-implicit time-stepping structure with an outer loop that involves a step to solve for the

F I G U R E 7
The convergence and consistency properties of the new framework of Section 3. Left: A convergence plot showing the change as a function of resolution of the L 2 error norm between the initial and final mixing ratio fields for the convergence configuration. The legend displays the gradients of lines of best fit through the points, approximating the order of accuracy. The consistent scheme shows a similar order of accuracy to an advective transport scheme. Right: A time series of the L 2 norms of the d and m X field from the consistency configuration. These have been normalised by subtracting and then dividing by the initial values. The legend displays the range (maximum minus minimum) for the whole time series of each line, showing that the mixing ratio field remains constant up to machine precision even while the d field is highly deformed [Colour figure can be viewed at wileyonlinelibrary.com] transport terms of the prognostic variables and an inner loop to solve a linearised problem for the implicit wave-like terms. This section presents results from a moist rising bubble test, based on the benchmark case of Bryan and Fritsch (2002). This test uses two the moisture species water vapour and cloud liquid (whose mixing ratios are represented by m v and m c respectively), which can convert between one another via condensation and evaporation.
The test describes a saturated atmosphere in a vertical slice that is initially in hydrostatic balance, with a prescribed mixing ratio of total moisture m t and a profile that is neutrally stable when the effects of latent heating are included. A perturbation is added to this state, triggering a rising thermal. There is no Coriolis force or precipitation.
The set-up is broadly the same as that of Bryan and Fritsch (2002): the domain is of length L x = 20 km and height H z = 10 km, and the balanced background state is defined by a pressure p = 10 5 Pa at the surface, along with the wet equivalent potential temperature that is e = 320 K everywhere. One difference is that the total moisture mixing ratio m t varies with height. This is to avoid fortuitous conservation: when m t is a constant field then lack of conservation in the transport of each of the two moisture species can exactly compensate for one another, resulting in a total moisture mass that is conserved. The total moisture profile is then given by with m 0 = 0.0175 kg⋅kg −1 and Δm = 5 × 10 −3 kg⋅kg −1 . The water vapour mixing ratio is set to its saturation value (with the remainder of the moisture being cloud liquid). However, the GungHo dynamical core has some thermodynamic simplifications compared with that in Bryan and Fritsch (2002). It uses a latent heat L v that is constant with respect to temperature T, and the heat capacities c v and c p do not include contributions from the moist component of the air. With these assumptions, a different definition of the wet equivalent potential temperature is more appropriate to define the initial conditions: The same perturbation of Bryan and Fritsch (2002) is then applied, but to this definition of e . To find the initial values of the prognostic variables from these requirements, the routine of Bendall et al., (2020) is used. Condensation and evaporation are applied through a saturation adjustment step at the end of each dynamical core time step; the details of this routine are described by Bendall (2019), although here it is used with constant L v .
This test is also used to compare the dynamical core using an advective form of the moisture transport with the new consistent and conservative approach. For these simulations, the grid spacing was Δx = Δz = 100 m and Δt = 2 s. The inner loop at the heart of the dynamical core changes d but not m X , and so if̃X is recomputed at the end of the inner loop the original mass of moisture will generally not be conserved. However, the dynamical core will be numerically conservative over the whole time step when the solver for the implicit terms has converged (and when using a conservative transport scheme). Therefore, to demonstrate the conservation with this test, the dynamical core was run with 16 × 1 outer and inner iterations, rather than the usual 2 × 2 that is used in GungHo. The time series of the evolution of the total mass of moisture for these cases are shown in Figure 9, showing that the new framework preserves the mass of moisture when the simple advective form does not. The diagnosed e fields at t = 800 s are shown in Figure 10, with the fields for both transport schemes appearing visually similar.

SUMMARY
The trilemma of the moist Charney-Phillips grid is that any choice of variable staggering appears to eliminate a desirable property of the whole discretisation. To address this, we have presented a framework to transport moisture conservatively and consistently with dry density, while still co-locating the mixing ratio with the potential temperature. This framework uses a vertically shifted mesh on which the moisture is expressed as a moisture density, which is transported using the same mass flux that was used to transport the dry density. The consistency between the transport of the dry and moist densities relies on a few simple properties of the operators that are involved in the scheme.

F I G U R E 9
The conservation of moisture mass using the new framework presented in Section 3, compared with a transport scheme using a standard advective form for the mixing ratio fields. Left: Time series of the mass M of the moisture species from the highest resolution simulation of the convergence configuration of Section 6.1, with Δx = Δz = 10 m. Right: Time series of the total mass of moisture from the rising bubble test of Section 6.2. Both plots show that the mass of moisture is conserved much better by the consistent, conservative scheme than the advective transport scheme. The legends display the ranges (maximum minus minimum) of the time series. The masses have been normalised by subtracting and then dividing by the initial value F I G U R E 10 The diagnostic e fields at t = 800 s from the moist rising bubble test of Section 6.2. Left: Tthe final field from a simulation using the advective form of the transport equation for the moisture species. Right: The final field from a simulation using the consistent transport framework for the moisture species. Contours are spaced every 1 K, and the 320 K contour has been omitted. The two fields are visually similar [Colour figure can be viewed at wileyonlinelibrary.com] The results in Section 6 demonstrate that the framework presented in this article does indeed provide the desired properties, without compromising on the benefits of co-locating moisture with . This was done in the context of a lowest order finite-element model, but can also be straightforwardly applied to a finite-difference or finite-volume discretisation. In further work we intend to address how this framework can be applied to a finite-element discretisation of higher polynomial order in the vertical direction. The challenge here is that thẽ V space on the shifted mesh would not necessarily have the same number of DoFs as V when the finite elements have a higher polynomial order. This problem could be bypassed if the dynamical core used a higher resolution grid of lowest order polynomials for its transport stage. An alternative approach would be to use a finite-element-based transport scheme, formulating Equation (3) in weak form and solving for the mixing ratio. Further work will also aim to show the benefits of solving the trilemma in the context of a full numerical weather prediction and climate model.

AUTHOR CONTRIBUTIONS
Thomas M. Bendall: formal analysis; investigation; software; writing -original draft; writing -review and editing. Nigel Wood: conceptualization; formal analysis; supervision; writing -review and editing. John Thuburn: conceptualization; formal analysis; writing -review and editing. Colin J. Cotter: conceptualization; formal analysis; writing -review and editing.