A Compatible Finite Element Discretisation for the Moist Compressible Euler Equations

We present new discretisation of the moist compressible Euler equations, using the compatible finite element framework identified in Cotter and Shipton (2012). The discretisation strategy is described and details of the parametrisations of moist processes are presented. A procedure for establishing hydrostatic balance for moist atmospheres is introduced, and the model's performance is demonstrated through several test cases, two of which are new.

that we aim to solve in Section 2. Then Section 3 describes the spatial and temporal discretisation of the dynamical part of the model, presenting the function spaces for the variables and the strategy for solving the equations of motion. Section 4 then details the specifics of the moisture parametrisations, including a discussion of how to combine fields from different function spaces within a parametrisation. After presenting our approach to initialising states with hydrostatic balance in Section 5, we then demonstrate its use in a number of test cases in Section 6. Finally, the model presented in this paper is written as part of Gusto, a library for dynamical cores using compatible finite element discretisations. Gusto itself is based on the Firedrake software of [7], the development of which is based at Imperial College London. This software provides automated code generation for finite element methods. Firedrake also has the functionality for tensor product element and extruded mesh functionality, as described in [8], [9] and [10], which we used substantially.

Governing Equations
We solve the compressible Euler equations featuring three species of moisture: water vapour, cloud water and rain. The ice phase is neglected. Motivated by the UK Met Office's most recent dynamical core, ENDGame, [11], our prognostic variables are the density of dry air ρ d and the dry virtual potential temperature θ vd as prognostic variables, alongside the wind velocity v and the mixing ratios r v , r c and r r . Here the subscripts respectively denote water vapour, cloud water and rain, whilst the mixing ratio r i is the ratio of the density by volume of the i-th substance to that of dry air, i.e. r i := ρ i /ρ d . The total mixing ratio of water species is r t := r v + r c + r r . The dry potential temperature θ vd is defined for r v , temperature T and air pressure p by where c pd is the specific heat capacity of dry air, p R is a reference pressure, θ d is the dry potential temperature and := R d /R v is the ratio of specific gas constants of dry air and water vapour. The choice of θ vd is motivated in [11], which notes that it is the more natural choice of variable to complement ρ d , and claims that it may be smoother than the dry potential temperature θ d .
The full equation set that we use is Dv Dr v Dt , Here Π is the Exner pressure function and Φ represents the geopotential, whilst f = f k represents the Coriolis parameter multiplied by the vertical unit vector. For the specific heat capacities c pd , c vd , c vml and c pml , the specific gas constant R m and also the latent heat of vaporization L v (T ), we follow closely the values used in [12], which are also listed in the appendix. The advective derivative is given by 3) The equation of state can be written as with κ := R d /c pd .
The terms on the left hand sides of (2.2) represent the dynamics, whilst the right hand sides are considered to be the physics. The processesṙ c cond ,ṙ accr ,ṙ r evap ,ṙ accu and S are the microphysics parametrisations and are described in Section 4.
The momentum equation can also be recast in vector invariant form: The final thing to note is the extra term proportional to ∇ · v appearing on the left hand side of (2.2b). This term is neglected in many models but mentioned in [13] and [12] to be important in fully capturing convection, particularly in a saturated atmosphere. In our model it appears in the forcing step of the dynamical core.

Function Spaces
One of the main results of [2] was the suggestion of combinations of function spaces for the velocity v and the (dry) air density ρ d within the compatible finite element framework that will satisfy the properties described in [1]. These spaces mimic the Arakawa C-grid staggering that allowed for steady geostrophic modes, avoidance of spurious pressure modes and more accurate representation of the dispersion relation for Rossby and inertia-gravity waves. Compatible sets of finite element spaces for v and ρ d will form part of a discrete de Rham complex, so that ∇ · v is in the same space as ρ d . Some such families of spaces can be found in the periodic table of finite elements [14].
In this section we will expand on this to list the function spaces (V ρ , V v , V θ ) that we use in our model for velocity, density and potential temperature respectively. We shall consider both two-dimensional vertical slices with quadrilateral cells and also three-dimensional domains using hexahedral elements with quadrilateral faces. The stratification of geophysical fluids along with their high aspect ratio motivates using an extruded mesh in the vertical, meaning that the mesh has regular layers. The finite elements on such a mesh can be constructed as tensor product elements, as described in [8]. We will also present two configurations of the model, using two different sets of spaces from the same family of finite elements spaces. We label these two sets as the lowest-order k = 0 spaces and the next-to-lowest-order k = 1 spaces. More explanation on these can be found in [15]. All of these spaces are illustrated in Table 1.
The density ρ d lies in a space that is discontinuous between elements. For the k = 0 spaces, density is constant within each cell, while it is linear in each cell for the k = 1 spaces. These are the discontinuous Galerkin spaces, denoted by dQ k .
The velocity v has continuous normal components between cells but discontinuous tangential components. The spaces that we use are the Raviart-Thomas spaces, RT k , [16]. For the k = 0 configuration this means that for each component, the field is linear in the direction of that component and continuous between cells in that direction, but constant within the cell in other directions. This becomes continuous and quadratic in the direction of the component for the k = 1 spaces, but discontinuous and linear in other directions.
To replicate the Charney-Philips grid, the potential temperature θ vd is co-located with the vertical component of velocity. This is the tensor product element of the discontinuous Galerkin element dQ k over the horizontal part of the domain with the continuous Galerkin element Q k+1 on the vertical part of the domain. We therefore denote this space by dQ k ⊗ Q k+1 . Moisture variables also lie in this space, as they can be considered as a dynamical adjustment to θ vd and so as to easily facilitate latent heat transfer as the water changes phase.

Overview
Here we present an outline of our model structure. It is based the UK Met Office ENDGame model [11] in terms of the semi-implicit formulation and physics-dynamics coupling, with the exception that Eulerian advection schemes are used with Eulerian averages in the semi-implicit time discretisation. This is in contrast to the ENDGame semi-Lagrangian approach. We illustrate this with some pseudocode that is similar to that presented in [15]. The pseudocode describes how we evolve our prognostic variables, which we denote by the single vector χ = (v, ρ d , θ vd , r v , r c , r r ).
Overall the model is a mixture of explicit and semi-implicit treatments. Terms in 2.2 are divided into forcings, F(χ), advection terms A(χ) and physics terms P(χ), with forcing being semiimplicit while advection and physics are explicit.
Given χ n , the state at the n-th time step, the first stage of the time step is to apply the explicit component of the forcing terms, returning χ * . Then we enter two loops. In the outer loop, the advecting velocity u is determined to be the average of v n and the best guess of the velocity at the time step, v n+1 p . The variables are then advected explicitly by u, returning χ p . To solve the implicit component of the forcing terms we enter the inner loop. The residual ∆χ is calculated between χ p and the implicit part of the forcing, so that the implicit part of the model is solved if ∆χ = 0. We iterate towards the solution by solve a linearised system of equations, with the right-hand-side as ∆χ. Once the outer loop is completed, the physics processes are treated explicitly in their own stage.
The pseudocode summarising the algorithm is: Table 1: An illustration of the function spaces that we use for density, velocity and potential temperature variables in vertical slice (d = 2) and three-dimensional (d = 3) simulations. We have have two configurations for each case: with vertical and horizontal degree of k = 0 and k = 1. Degrees of freedom on cell facets are shared between cells and denote the continuity of the field between them. Our model is designed for use with both the k = 0 and k = 1 sets of function spaces. The formulation is slightly different between these two cases, and here we give a summary of those differences. These are listed in Table 2.
The main area of difference is in the advection part of the model, which is briefly detailed  Table 2: A summary of the components of the different model configurations when using Gusto with the lowest order k = 0 spaces, and the next-to-lowest order k = 1 spaces.
in Section 3.4. In the k = 0 case, we use the recovered space transport schemes outlined in [15]. However when using the k = 1 spaces, a discontinuous Galerkin upwind scheme is used for the transport of ρ d , whilst the advection of θ vd and the moisture variables is performed by the embedded scheme introduced by [17]. The velocity v advection equation is written in vector invariant form, and uses the theta time stepping method. All other advection schemes use a three-step Runge-Kutta time stepping procedure (SSPRK-3) [18].
The other difference is that the limiters used with the advection schemes are different between the two set-ups.

Forcing
The 'forcing' operation constitutes the sum of the non-advective terms in the differential equations (2.2a) and (2.2b), which are written in weak form. Our goal is to find F(v) and F(θ vd ).
For both k = 0 and k = 1 configurations, the forcing that we apply to v is the solution where Ω is the domain, Γ is the set of all interior facets and the angled brackets · denote the average value on either side of a facet. The double square brackets f n = f + · n + + f − · n − denote the jump in values between either side of the facet, arbitrarily labelled with + and −. The forcing term F(v) is then taken to the solution v trial to this.
The forcing F(θ vd ) applied to θ vd is the solution θ trial , for all γ ∈ V θ to

Advection
In the advection stage, each of the variables is translated by the velocity u. We represent action this by the operator Aū. Here we briefly provide details of the advection schemes used. For a summary, see Table 2.

Discontinuous Galerkin Upwind Advection
First, we define a single forward-Euler step of discontinuous Galerkin upwinding, which we describe as the operation Lū upon a scalar quantity q. For the advective form of the transport equation this involves finding the solution q trial , for all p ∈ V q where the jump is p + = p + − p − . The outward normal on the + side of the facet is denoted by n + and q * is the upwind value identified by This equation can also be cast in continuity form, as it will be for the transport of ρ d . The SSPRK-3 time stepping scheme presented in [18] involves composing these steps as follows: This is the scheme used for transport for ρ d in the k = 1 configuration, and is the basis for the embedded DG and recovered advection schemes.

Embedded DG Advection
The embedded DG transport scheme introduced in [17] is used for advecting θ vd and the moisture variables in the k = 1 configuration. The fields are first injected into fully discontinuous analogues of their spaces, where the advection takes place. The fields are then projected back to their original spaces. See [17] for more details.

Recovery Operator and Recovered Advection
The advection schemes used for the k = 0 spaces were presented in [15]. This scheme extends the embedded DG scheme described in the previous section: the fields are 'recovered' from the original space to spaces with higher-degrees of polynomial. The transport then occurs in a discontinuous higher-degree space. As in the embedded scheme, this field is projected back to the original space to complete the advection. The spaces used as part of this scheme were described in [15].
The 'recovery' operator is also used in physics-dynamics coupling, which is described in Section 4. The basic operation recovers a field from some function space V 0 to V 1 by pointwise evaluation of the field at degrees of freedom. When V 0 is discontinuous between elements but V 1 is continuous, the new value of the field is the average over the neighbouring elements. At the exterior boundaries of the domain, this process will not have second-order numerical accuracy and a further extrapolation step is required to accurately represent gradients.

Vector invariant advection
The most major difference in the advection stage of the model between using the k = 0 and k = 1 spaces is the transport of the velocity. We write the advection equation in vector-invariant form, discretising ∂v ∂t The action of the advection operator Lū gives the with v * as the upwind value of v and · + taking the same definition as in Section 3.4.1.
The theta time stepping scheme is used, so that the velocity at the (n + 1)-th time step is found by solving We take ϑ = 1/2.

Limiting
The advection schemes that we use do not preserve the monotonicity property of the continuous transport equation. This can result in unphysical negative concentrations when advecting moisture variables. A way of preventing this is to apply a slope limiter.
We provide the option of adding a limiter to the advection of variables in V θ . The vertexbased limiter of [19] is designed for use on fields that are piecewise-linear and discontinuous between cells. It prevents the formation of new maxima and minima by separating the field into a constant mean part and adjusting the perturbed part in each cell. It is applied to the field before the advection operator Vū first acts each time step and after each stage of the three-step Runge Kutta method.
For advection of moisture variables in the k = 0 set of spaces, we use the recovered space scheme and the advection can be limited by the vertex-based limiter of [19]. When projecting back from the advecting space V 1 , as described in [15] we split the projection into two steps: projecting into the broken space V 0 before recovering any continuity.
With the k = 1 spaces, the vertex-based limiter of [19] is again used. However this field is quadratic in the vertical and the values of the field at vertical midpoint of the cell are unchanged by the vertex-based limiter. We bound these values by restricting them to the average of the values at the adjacent vertices if it falls outside of the range spanned by them. In other words, if r i+1/2 is the value of the field at a degree of freedom halfway up the i-th cell, then the new value is 3.5 Solve Figure 1: The final fields, following one rotation around the domain, of functions in the twodimensional dQ 1 ⊗ Q 2 space that were initialised with the slotted-cylinder, hump, cone initial condition of [20]. (Left) with the limiter applied, and (right) without the limiter. Overshoots and undershoots are highlighted in yellow.
To illustrate the effectiveness of this limiter, we performed the rotation of the slotted-cylinder, hump, cone of [20] in a 2D vertical slice. Figure 1 shows the fields in the V θ space for the k = 1 configuration, rotated once around the domain by a velocity defined by the stream function The domain used was the unit square, with ∆x = ∆z = 0.01, whilst the time step was ∆t = 10 −4 . The final fields were recorded at t = 1.

Solve
We now present the strategy for the linear solve step of the model. When initialising the model, we set mean ρ d and θ vd states which we denote respectively as ρ d and θ vd . The residual ∆χ between the predicted explicit χ p and the latest guess χ n+1 p of the prognostic variables is computed first. Reducing this residual to zero solves the implicit part of the model. This is done by solving a linear set of equations, with the residual on the right hand side, to calculate increments that are added to χ n+1 p , updating it.
The linearised problem that we are attempting to solve is where the primes represent the perturbations to be found.
These equations will be solved in weak form. Howver, we simplify them first by eliminating θ vd . This allows us to create a hybridised mixed system to be solved for v trial and ρ trial . Taking k as the unit upward normal, we introduce The Exner pressure perturbation is then approximated as (3.13) The hybridised mixed system involves breaking the continuity of normal components of v between elements, introducing the broken space V v . In conjunction, a trace field trial is introduced approximating the average values of Π perturbations on the trace space V trace : the discontinuous space of functions that live on the boundary of V v . This variable acts as Lagrange multipliers to provide the continuity constraint of the pre-hybridised system. The resulting mixed system for (3.14) where w = k · ψ is the vertical component of the velocity test function and ∂Ω is the external boundary of the domain. This mixed system can be statically condensed into a single system for the Lagrange multipliers, giving a single system to be solved. Once this system has been solved, the calculated value of is used to find ρ d and then the broken v . The recovery operator (averaging values between cells) is then used to restore the continuous velocity v . Finally, the value of θ vd is found as the θ trial that solves for all More on this solver and its performance will be detailed in future work.

Physics Parametrisations
In this section we consider the discretisation of the 'physics' terms in equation set (2.2). Each term is evaluated individually in a separate 'physics' routine. The term labelledṙ c cond represents 4.1 Combining Fields From Different Function Spaces the condensation of water vapour into cloud water and the evaporation of cloud water into water vapour. We use the rate given in [21] and [12], with the saturation mixing ratio stemming from Tetens' empirical formula [22], giving where C sat 0 , C sat 1 and C sat 2 are constants whose values are given in the appendix. The pressure p will be explained in the following section.
The similar termṙ r evap describes the evaporation of rain water into water vapour, with the value taken from [23].
The coalescence of cloud water into rain droplets is described byṙ accr andṙ accu : the accretion and auto-accumulation processes respectively. The formulas that we use are also those described in [23], with values from [24].
The sedimentation of rain is given by S. Our approach to determining S is similar to the single moment scheme described in [25]. We assume that the number n r (D) of raindrops of diameter D forms a spectrum described by a Gamma distribution. This is used to determine the terminal velocity of raindrops averaged over D. The rain field r r is then advected by this averaged velocity using the same advection scheme used by the dynamics for r r described in Section 3.

Combining Fields From Different Function Spaces
Many diagnostic fields, such as temperature T or saturation mixing ratio of water vapour r sat must be determined from prognostic fields in different function spaces, usually ρ d , θ vd and r v . Such diagnostic fields are important in the equations used for the physics parametrisations and can also be important in setting the initial state of the model.
As ρ d and θ vd lie in different function spaces, when determining some diagnostic variable q there are a number of different choices that can be taken. With our emphasis on parametrisations of moist processes, we take only the case of q ∈ V θ , so that it is also in the same function space as the moisture variables. Our approach to combining these fields is to 'recover' ρ d is into V θ using the recovery operator outlined in Section 3.4.3 to give ρ d . Then q can be calculated algebraically within V θ . This approach can be used with both the k = 0 and the k = 1 sets of spaces, and with the k = 0 spaces it has second-order numerical accuracy (including at the boundaries of the domain). For example, the temperature T and pressure p used in parametrisations are calculated as At present we are performing simple explicit first-order integration in time for the physics routines. In other words for a processṙ affecting a variable r, the new value r new will be related to the old r old by r new = r old + ∆tṙ(r old ).
The exception here is the treatment of the sedimentation of rainfall, which was described in the previous section. The value ofṙ comes from the state of the model just before this physics routine is called. This is the state of the model after either the dynamics or the preceding physics routine has been completed. The physics routines are executed consecutively, (this is commonly known as "sequential splitting"), in the following order: 1. accretion of cloud water, 2. auto-accumulation of rain water, 3. sedimentation of rain, 4. evaporation of rain water, 5. evaporation of cloud water/condensation of water vapour.
We choose to do the evaporation/condensation step last so as to prevent any supersaturation at the end of the time step.

Hydrostatic Balance Routines
For many test cases the background or initial state of the model will be in hydrostatic balance. The procedure that we use to obtain a discrete hydrostatic balance is based upon that presented in [4]. Given a boundary condition for the pressure and a θ d field, this procedure finds the ρ d which gives rise to zero vertical velocity. In this section we describe two developments to this: the extension of the routine of [4] to use a hybridised formulation and the treatment of cases in which thermodynamic and moist variables and the prognostic variables need to be found.

Hybridised Hydrostatic Balance
The hybridised hydrostatic balance system uses a similar approach to that described in Section 3.5. As in that case, the variables are expressed in discontinuous spaces and Lagrange multipliers are introduced in a trace space on the horizontal facets to provide continuity constraints. We also introduceV v , the subspace of the broken velocity space V v with zero flow normal to the boundaries. The problem is then to find the (v trial , ρ trial , trial ) ∈ V v , V ρ , V trace , such that for where ∂Ω 0 is the boundary upon which the condition that Π = Π 0 is to be satisfied and Π is a linearisation of Π. Equation (5.1a) can be statically condensed to an elliptic equation for , which can be efficiently inverted in each column of the mesh. In fact it turns out that is an approximation for c pd θ vd Π/(1 + r t ) on the horizontal facets of the mesh. Once is determined, v and ρ d can be solved for via local back-substitution in each cell of the mesh. Continuity is restored in v again by using the recovery operator. By hybridising, we avoid solving a larger, more illconditioned mixed system and instead inverted a condensed elliptic problem. The implementation of this procedure is described by [26], with the application taking the form of a custom Pythonbase preconditioner conforming to standard PETSc library options [27,28], as described in [26]. The advantage being that solver options for solving the condensed system arising from reducing (5.1a) can easily be updated, even when nested inside a non-linear method.

Saturated Conditions
This set-up involves initial conditions like those in the moist benchmark of [12]. The problem is to find ρ d , θ vd and r v given θ e , r t and a boundary condition on the pressure. Assuming the absence of rain, r v = r sat , and r c = r t − r v . The wet-equivalent potential temperature, θ e , is a conserved quantity in reversible, moist adiabatic processes, i.e. Dθ e /Dt = 0. Following [29], for our equation set (2.2), θ e can be written as which with H = 1 for a saturated atmosphere is the same as that used by [12] and that derived in the appendix of [30] from the second law of thermodynamics. The saturation mixing ratio is given by (4.1).
The challenge is to obtain the θ vd and ρ d fields that satisfy the specified θ e field whilst ensuring that r v = r sat . We use an initial guess for θ vd and feed it to equation (5.1a), to generate a guess for ρ d . This density is then converted into V θ , before a nested fixed-point iteration-style procedure is used to invert θ e (θ vd , ρ d , r v ) and r sat (θ vd , ρ d , r v ), to obtain θ vd and r v .
Let l, m and n count the number of iterations to find ρ d , θ vd and r v respectively. These form nested loops, such that the outer loop uses the latest approximations of θ vd and r v with equation (5.1a) to obtain ρ (l+1/2) d . We then combine this with the previous value to update the approximation to ρ d , doing ρ

Unsaturated Conditions
where δ = 0.8. The next loop finds θ vd using whilst the inner loop computes where θ e without arguments denotes the specified value. These loops are iterated until θ e (θ vd , ρ d , r v ) converges to the specified value to some tolerance.

Unsaturated Conditions
We now discuss how to find the prognostic thermodynamic variables given θ d and the relative humidity H, such as the case in Section 6.3 or [31]. The relative humidity is related to r v by There is only one inner loop in this case, so that m = n. The new value of r v is found by rearranging 5.6 so that where H is the specified value of the relative humidity. The final step is to use the specified value of θ d to get θ (m) This iterative process continues until H θ has converged to its specified value to some tolerance.

Test Cases
In this section we demonstrate the discretisation detailed in previous sections through a series of test cases, with some comparison of the k = 0 and k = 1 configurations of the model. Two new test cases are presented, featuring a gravity wave in a saturated atmosphere and a threedimensional rising thermal in a saturated atmosphere.

Bryan and Fritsch Moist Benchmark
The first demonstration of our discretisation is through the moist benchmark test case of [12], which simulates a rising thermal through a cloudy atmosphere. The domain is a vertical slice of width L = 20 km and height H = 10 km. Periodic boundary conditions are applied at the vertical boundaries, but the top and bottom boundaries are rigid, so that v · n = 0 along them. As in [12], we include no rain microphysics and no Coriolis force.
The initial conditions defined in [12] specify a background state with constant r t = 0.02 kg kg −1 and constant wet-equivalent potential temperature θ e = 320 K, which is defined in (5.2). Along with these, the background state is given by the requirements of hydrostatic balance, r v = r sat and p = 10 5 Pa at the bottom boundary. The procedure described in Section 5.2 allows us to obtain the prognostic variables θ vd , ρ d , r v and r c that approximately satisfy these conditions.
A perturbation is then applied to θ vd . With (x, z) as the horizontal and vertical coordinates, the perturbed field is θ vd = ∆Θ cos 2 πr 2rc , r < r c , 0, otherwise, where ∆Θ = 2 K, r c = 2 km and with x c = L/2 and z c = 2 km we define The initial θ vd field is given in terms of the background field θ vd : In the test case of [12], the pressure field is unchanged by the perturbation. To replicate this with our prognostic variables, we find ρ d such that for all ζ ∈ V ρ , Ω ζρ d θ vd dx = Ω ζρ d θ vd dx, (6.4) where ρ d and θ vd are the background states for ρ d and θ vd . The system is returned to saturation by finding the r v that solves, for all φ ∈ V θ , where r sat ( ρ d , θ vd , r v ) is the expression for saturation mixing ratio in terms of the initial ρ d recovered into V θ and the initial θ vd field that have already been found, and also the unknown r v to be solved for. Finally, r c is found from applying r c = r t − r v . The initial velocity field is zero in each component. Figures 2 and 3 show the θ e and vertical velocity w fields at t = 1000 s. Figure 2 shows these fields for the configuration using the k = 0 lowest-order spaces, with the k = 1 spaces shown in Figure 3. Both simulations used ∆x = ∆z = 100 m and ∆t = 1 s. These final states are visibly different: whilst the k = 0 solutions resemble those of [12], the k = 1 solution displays an extra plume forming at the top of the rising thermal. We believe this to be a manifestation of a physical instability that is damped by numerical diffusion in the k = 0 case. The k = 1 solution appears highly sensitive to the choice of mesh, as at higher resolution the θ e field does not appear to converge to a single solution. Indeed if the domain is spanned horizontally by an odd number of cells, rather than a secondary plume emerging, the top of the primary plume appears to collapse. This behaviour is also observed in the absence of moisture.  [12] representing a thermal rising through a saturated atmosphere. The 320 K contour has been omitted for clarity in the θ e field. This simulation is with the k = 0 lowest-order set of spaces, with grid spacing ∆x = ∆z = 100 m and a time step of ∆t = 1 s. These solutions are visibly similar to those presented in [12]. Figure 3: Outputted fields from the k = 1 next-to-lowest degree space simulation at t = 1000 s of the moist benchmark from [12]. (Left) θ e with contours spaced by 0.25 K and (right) vertical velocity w contoured every 2 m s −1 . The simulation used grid spacing ∆x = ∆z = 100 m and a time step of ∆t = 1 s. The 320 K contour has been omitted for clarity in the θ e field. A second plume can be seen forming at the top of the primary plume.
We present here a new test case, a moist version of the non-hydrostatic gravity wave test of [32], but in a saturated atmosphere like that of the moist benchmark in [12]. The final state of this test is spatially smooth, making this test appropriate for convergence tests. No rain physics is used in this test case, and there is also no Coriolis force.
The problem is set up in a two-dimensional vertical slice of length L = 300 km and height H = 10 km. The dry gravity wave setup used in [32] applies a perturbation to θ d , which has a stratified background profile and is in hydrostatic balance. Our variation on this is to apply the perturbation to the stratified background profile of θ e in hydrostatic balance. Using (x, z) as the horizontal and vertical coordinates, the specified θ e profile is where Θ 0 = 300 K and N 2 = 10 −4 s −2 . With r t = 0.02 kg kg −1 everywhere and the boundary condition of p = 10 5 Pa at z = 0, we use the hydrostatic balance procedure laid out in Section 5.2 to find the θ vd , ρ d , r v and r c fields that correspond to these initial conditions with the requirements of hydrostatic balance and that r v = r sat everywhere. The initial velocity applied is v = (U, 0) with U = 20 m s −1 describing uniform flow in the x-direction. This defines each of the mean fields.
A perturbation is then added, which is specified as with a = 5×10 3 m and ∆Θ = 0.01 K. The perturbed initial condition is then given by θ e = θ e +θ e . Setting the new requirements that both r t and the pressure are unchanged by the addition of the perturbation, and that we still have r v = r sat , defines the problem necessary to solve to find the initial prognostic fields. We do this via a nested iterative process related to that described in Section 5.2. In the outer loop we find ρ h d such that for all ζ ∈ V ρ Ω ζρ h d θ n vd dx = Ω ζρ d θ vd dx, (6.8) which is combined with the previous best estimate of ρ n d to give ρ n+1 d = (1 − δ)ρ n d + δρ h d , where δ = 0.8. Nested inside this process are more damped iterations to find θ vd and r v , exactly as in Section 5.2. Figure 4 shows the perturbation to the final diagnostic θ e field at t = 3600 s for simulations with the k = 0 lowest-degree spaces with ∆x = ∆z = 500 m and the k = 1 set of spaces, where ∆x = ∆z = 1000 m, both with ∆t = 1.2 s. These different cases are not visibly different from one another, and closely resemble the final state of the dry case from [32].

Convergence
As the solution and evolution of this example is spatially smooth, we can use it to form a convergence test upon our model. Although the final state (at t = 3600 s) does not have an analytic solution, we use the fields from a high resolution simulation as the true solution. To measure the spatial accuracy of our model, we ran this test case at different resolutions, with each using a time step of ∆t = 1.2 s, which is small enough to avoid the Courant number breaching its critical value in the highest resolution cases. The error is measured looking at the θ e diagnostic in V θ at t = 3600 s. The θ e fields are interpolated onto the finest mesh, which has ∆x = 100 m for the k = 0 case but ∆x = 200 m for the k = 1 case. The error between the high resolution solution for θ e and those run at coarser resolutions is measured in the L 2 norm. Results for our model are plotted in Figure 5, which indicates that in both the k = 0 and k = 1 cases the model has second-order spatial accuracy, with the error proportional to (∆x) 2 .

Rising Thermal with Rain
This test case is based upon one described in [31]. This involves a thermal rising in an unsaturated atmosphere, forming a cloud which rains out. This is another two dimensional vertical slice test, this time with domain of height H = 2.4 km and length L = 3.6 km, again with periodic conditions at the vertical sides. The Coriolis force is neglected.
In contrast to the saturated atmosphere initial conditions of Sections 6.1, the initial state is defined by the dry potential temperature θ d and a relative humidity field H. The background fields are H = 0.2 everywhere and θ d = Θe Sz , (6.9) where Θ is the dry potential temperature corresponding to T surf = 283 K and p = 8.5 × 10 4 Pa, which also provides the pressure condition at the boundary. The stratification is given by S = 1.3 × 10 −5 m −1 . We then use the procedure outlined in Section 5.3 to find the background θ vd , ρ d and r v fields that satisfy hydrostatic balance. The initial r c and r r fields are zero.
The perturbation is then applied to the relative humidity field H, with a circular bubble that is just saturated, with an outer disk smoothing the perturbation into the background state. This initial relative humidity field is given by where H = 0.2 and Fields are displayed in Figures 6 and 7 for θ vd and r r at t = 300 and 600 s. Both simulations use the k = 0 lowest-order space set-up, with ∆x = 20 m and ∆t = 1 s. Figure 6 shows the results with no limiter applied to the advected moisture fields, whilst Figure 7 shows the same set-up but with a limiter applied to all the moisture variables.
The length scales of the simulation are small enough for it to be highly turbulent, with the final state dependent on the resolution in the absence of a turbulence parametrisation. Indeed, the lack of turbulence parametrisation in our model explains why these results look significantly different to those of [31]. Comparing Figures 6 and 7 demonstrates the effect of limiting the transport of moisture species. In the absence of the limiter, the mass of water vapour depreciates less and so more cloud is formed, associated with a greater release of latent heat and a stronger updraught. We see that rain forms earlier in the absence of a limiter. However, some negative moisture values do form, which are absent from the limited case.

Three-Dimensional Thermal in a Saturated Atmosphere
We now demonstrate the use of our discretisation upon small-scale dynamics in three dimensions. This test case is a three-dimensional version of the moist benchmark of [12] that was described in Section 6.1. Rain and the effects of planetary rotation are not included.
6.4 Three-Dimensional Thermal in a Saturated Atmosphere The domain is now periodic in the horizontal directions, with length, width and height 10 km. The background state set-up is the same as that in Section 6.1, with θ e = 320 K, r t = 0.02 kg kg −1 and the pressure as p = 10 5 Pa on the bottom surface. Using the initialisation procedure that was outlined in Section 5.2 generates the values of the prognostic variables such that the model is in hydrostatic balance and saturated with respect to water vapour.
Using x c = y c = 5 km and z c = 2 km, and defining we apply the perturbation θ vd = ∆Θ cos 2 πr 2rc , r < r c , 0, otherwise, (6.13) Figure 7: The same results as in Figure 6, but with a limiter applied to the moisture species.
with ∆Θ = 1 K and r c = 2 km. As in Section 6.1, the perturbation is applied using the background field θ vd θ vd = θ vd 1 + θ vd 300 K . (6.14) The same routine as used in Section 6.1 is then applied to obtain the initial ρ d , r v and r c fields, ensuring that the atmosphere is exactly saturated and that the initial pressure field is equal to the background pressure field.
Cross-sections of the θ e and vertical velocity w fields at t = 1000 s and y = 5 km are shown in Figures 8 and 9, for both the set-ups using k = 0 lowest-order spaces (which had ∆x = ∆y = ∆z = 100 m) and the k = 1 spaces (which had ∆x = ∆y = ∆z = 200 m). Both simulations had ∆t = 1 s. As in Section 6.1, we see a secondary plume beginning to form at the top of the rising thermal in the k = 1 case.

Moist baroclinic wave
The final test case that we present is the moist baroclinic wave outlined in [33]. This is the only test featuring the Coriolis force, although rain is again neglected. A three dimensional channel of  This test case uses initial conditions that are analytically in thermal wind balance. In [33], the vertical coordinate used is a pressure coordinate η. This is used to define the background zonal wind u, the geopotential Φ, the virtual temperature T v and the specific humidity q according to the following equations: The constants take the values Γ = 0.005 K m −1 , f 0 = 2.00×10 −6 s −1 , a = 6.37×10 6 m, T 0 = 288 K, u 0 = 35 m s −1 , b = 2, ∆y w = 3.2 × 10 6 m and η w = 0.3. We use the slightly lower value of q 0 = 0.016 than [33] to prevent our model being initially too close to saturation.
This initial condition is converted to our prognostic variables using the Newton iteration procedure suggested in [33] to find η. Using the requirement that Φ = gz and taking η ∈ V θ , the procedure uses: This η field is then used to compute u, T v and q. In order to convert T v and q into θ vd and r v , we use the definitions of θ vd , η = p/p s for p s = 10 5 Pa and the equations and These θ vd and r v fields are used to compute the ρ d field that provides hydrostatic balance for the background state using the procedure outlined in Section 5.
The perturbation is added to the x-component of the background velocity, which we will denote as u, so that u = u + u . With x c = 2 × 10 6 m, y c = 2.5 × 10 6 , L p = 6 × 10 5 , u p = 1 m s −1 we use r = (x − x c ) 2 + (y − y c ) 2 , (6.19) to get u = u p exp − r L p 2 .
(6.20) Figures 10 and 11 show fields from this test case at t = 12 days for the k = 0 and k = 1 spaces respectively. For the k = 0 configuration, ∆x = ∆y = 200 km and ∆z = 1 km, whilst for k = 1 these were ∆x = ∆y = 250 km and ∆z = 1.5 km. For both simulations, ∆t = 300 s. In the k = 1 simulation, the baroclinic wave becomes much stronger than in the k = 0 simulation. As the wave develops, the maxima and minima in the temperatures in the k = 1 case are higher than for the k = 0 lowest-order spaces. When these minima coincide with regions close to water vapour saturation, more condensation occurs. This releases latent heat and strengthens the baroclinic wave, reinforcing the behaviour.  : Cross-sections of fields at t = 12 days from the moist baroclinic wave test case using the k = 1 spaces. The grid sizes used were ∆x = ∆y = 250 km and ∆z = 1.5 km, with ∆t = 300 s. Plots shown are the same as in Figure 10, but note the different scale in the zonal wind perturbation plot (bottom right).

Conclusions
We have presented a discretisation of the moist compressible Euler equations that uses a compatible finite element framework. Building upon the work of [2], [4], [5] and [15], we have presented two sets of compatible spaces that can be used. The model configurations for each of the sets of spaces has been described, detailing the discrete equations that are solved and the transport schemes used, including slope limiters that can be used with moisture variables. A discussion is given of the parametrisation of the moist processes and how these are coupled to the description of the resolved flow. There is also a presentation of how a discrete hydrostatic balance can be set-up for the model under moist conditions. The performance of the model is displayed through the presentation of five test cases, including two new ones which have introduced.
Future work planned includes testing the choices made in these models and comparing the qualities of the model using each of the different sets of spaces. More physics parametrisations will be added and the procedure used in the linear solve stage of the model will be thoroughly

Notation, thermodynamic variables and constants
Specific heat capacity of dry air at const. volume c vd 717 J kg −1 K −1 Specific heat capacity of dry air at const. pressure c pd 1004.5 J kg −1 K −1 Specific heat capacity of water vapour at const. volume c vv 1424 J kg −1 K −1 Specific heat capacity of water vapour at const. pressure c pv 1885 J kg −1 K −1 Specific heat capacity of liquid water at const. pressure c pl 4186 J kg −1 K −1 Specific heat capacity of moist air at const. volume c vml c vd + r v c vv + (r c + r r )c pl Specific heat capacity of moist air at const. pressure c pml c pd + r v c pv + (r c + r r )c pl Specific gas constant for dry air R d 287 J kg −1 K −1 Specific gas constant for water vapour R v 461 J kg −1 K −1 Specific gas constant for moist air R m R d + r v R v Reference latent heat of vaporisation of water at T R L vR 2.5×10 6 J kg −1 Latent heat of vaporisation of water L v L vR − (c pl − c pv )(T − T R ) Reference temperature T R 273. 15