# Energy conserving SUPG methods for compatible finite element schemes in numerical weather prediction

We present an energy conserving space discretisation based on a Poisson bracket that can be used to derive the dry compressible Euler as well as thermal shallow water equations. It is formulated using the compatible finite element method, and extends the incorporation of upwinding for the shallow water equations as described in Wimmer, Cotter, and Bauer (2019). While the former is restricted to DG upwinding, an energy conserving SUPG scheme for the (partially) continuous Galerkin thermal field space is newly introduced here. The energy conserving property is validated by coupling the Poisson bracket based spatial discretisation to an energy conserving time discretisation. Further, the discretisation is demonstrated to lead to an improved temperature field development with respect to stability when upwinding is included. An approximately energy conserving full discretisation with a smaller computational cost is also presented.

## Authors

• 1 publication
• 13 publications
• 5 publications
10/04/2019

### A Compatible Finite Element Discretisation for the Moist Compressible Euler Equations

We present new discretisation of the moist compressible Euler equations,...
03/12/2022

### Energy conserving particle-in-cell methods for relativistic Vlasov–Maxwell equations of laser-plasma interaction

Energy conserving particle-in-cell schemes are constructed for a class o...
02/26/2021

### An HPC-Based Hydrothermal Finite Element Simulator for Modeling Underground Response to Community-Scale Geothermal Energy Production

Geothermal heat, as renewable energy, shows great advantage with respect...
03/12/2021

### Energy stable and accurate coupling of finite element methods and finite difference methods

We introduce a hybrid method to couple continuous Galerkin finite elemen...
05/03/2021

### High-order space-time finite element methods for the Poisson-Nernst-Planck equations: Positivity and unconditional energy stability

We present a novel class of high-order space-time finite element schemes...
12/21/2019

### A structure-preserving approximation of the discrete split rotating shallow water equations

We introduce an efficient split finite element (FE) discretization of a ...
06/19/2019

### Dynamic coupling between particle-in-cell and atomistic simulations

We propose a method to directly couple molecular dynamics, finite elemen...
##### This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

## 1 Introduction

Finite element methods have recently gained an increased interest in numerical weather prediction (NWP), as they allow for higher order discretisations and more general meshes, thus avoiding the parallel computing issues associated with grid poles. This includes spectral element methods, discontinuous Galerkin methods, and the compatible finite element method [COTTER20127076], where in the latter finite element spaces are mapped to one another via differential operators [natale2016compatible]. In the context of NWP, the compatible finite element method can be seen as an extension of the Arakawa C finite difference grid, and a dynamical core based on it is currently in development at the UK Met Office [ford2013gung], due to replace the current finite difference latitude longitude mesh discretisation.

An important aspect of discretisations in NWP, particularly for climate simulations, is conservation of quantities such as mass and energy. While the former can be conserved using a suitable discretisation of the continuity equation, the latter requires a careful discretisation of all prognostic equations, ensuring that the energy losses and gains are balanced between the discretised terms. One way to guide this process is to consider the equations in a Hamiltonian framework [morrison1998hamiltonian], where the Hamiltonian is given by the system’s total energy, and the equations are inferred by a Poisson bracket. Conservation of energy then follows directly from the bracket’s antisymmetry, and any space discretisation maintaining this property will then also conserve energy. In particular, this framework facilitates the construction of advection schemes, which may otherwise not have an easily accountable effect on the total energy development. For the dry compressible Euler equations, which form the basic equation set of dynamical cores in NWP, such Hamiltonian based discretisations already exist e.g. for hexagonal finite difference C-grids [gassmann2013global] (which includes higher order transport schemes) and the compatible finite element method [lee2020mixed].

For finite element methods, the choice of advection scheme depends on the underlying space, and for the compatible finite element method, this requires different formulations for the range of different spaces in use, including continuous and discontinuous ones. Within the Poisson bracket framework, DG upwinding for the depth field in the rotating shallow water equations as well as for the buoyancy field in the thermal shallow water equations has already been introduced in [WIMMER2019109016] and [ELDRED20191], respectively. Further, an upwinded energy conserving velocity advection scheme was presented in the context of Lie derivatives for the incompressible Euler equations in [natale2016variational], which was extended to the shallow water case in [natale2016compatible] and the shallow water case from a Poisson bracket point of view in [WIMMER2019109016]. Finally, an energy conserving SUPG advection scheme was considered for potential vorticity advection in the shallow water equations in [BAUER2018171], exploiting the fact that the Poisson bracket term corresponding to velocity advection is antisymmetric in itself.

These energy conserving upwinded advection schemes can, with the exception of the potential vorticity scheme, in principle be implemented readily for the three-dimensional dry compressible Euler equations. However, an upwinded, energy conserving potential temperature advection scheme is still missing if a continuous (or partially continuous) finite element space is used for the temperature field. In particular, this is the case for the finite element equivalent of a Charney-Phillips finite difference grid, where the space’s node locations are set to coincide with the velocity space nodes corresponding to the vertical velocity [natale2016compatible]. This choice of temperature space was recommended in a dispersion property study in [melvin2018choice] for mixed finite element methods in NWP when compared to a fully discontinuous or fully continuous finite element space, and will be used in the UK Met Office’s next dynamical core.

Advection terms for fields in continuous or partially continuous finite element spaces can be discretised using the Streamline Upwind Petrov Galerkin (SUPG) method [brooks1982streamline], where a diffusive term is added along the direction of the flow to all test functions. In this paper, we present a discretised almost-Poisson bracket that includes SUPG upwinding for the thermal field, which, depending on the choice of Hamiltonian, leads to an energy conserving space discretisation for the dry compressible Euler equations or the thermal shallow water equations. The energy conservation property as well as the thermal field’s qualitative development are then verified in numerical test cases, using an energy conserving time discretisation as described in [cohen2011linear]. For simplicity, while the formulation is valid in three dimensions, we consider two-dimensional scenarios in this paper. This includes vertical slice test cases of the Euler equations with a Charney-Phillips type temperature space, as well as a spherical test case of the thermal shallow water equations with a fully continuous buoyancy space. Further, in view of the increased computational cost to achieve energy conservation to machine precision, we also consider a simplified discretisation, which relies on an approximation to the Poisson integrator and a small adjustment in the almost Poisson bracket. To arrive at a fully upwinded bracket, the latter additionally includes upwinding in density and velocity.

The rest of the paper is structured as follows: In Section 2, we review the Poisson bracket as well as the SUPG scheme, and present the space discretised formulation including upwinding. In Section 3, we describe the full discretisation, introduce the simplified discretised form, and present numerical results. Finally, in Section 4, we review the formulation and the corresponding numerical results, and discuss ongoing work.

## 2 Formulation

In this section, we first review the Poisson bracket, which we present in its continuous form. We then describe how to derive the resulting sets of equations, considering Hamiltonians that lead to the dry compressible Euler and the thermal shallow water equations. Next, we describe the finite element space discretisation, including the Streamline Upwind Petrov Galerkin stabilisation scheme for advection terms, as well as a mixed SUPG/DG upwinding scheme for the Charney-Phillips type finite element space, as described in [natale2016compatible]. Finally, we derive a set of upwinded equations within the Poisson bracket framework, and further discuss how to incorporate existing upwinding schemes for the velocity and density fields in the bracket.

### 2.1 Hamiltonian formulation

Many fluid dynamical equations can be formulated within a Hamiltonian framework, using the system’s Hamiltonian , i.e. the total amount of energy, and a Poisson bracket , which is an antisymmetric bilinear form that satisfies the Jacobi identity (for more details, see [shepherd1990symmetries]). The time evolution of any functional of the prognostic variables is then given by

 dFdt={F,H}. (2.1.1)

Here, we consider prognostic variables , , and in a two-dimensional domain , which denote the flow velocity, a density related field, and a thermal type field, respectively. Further, we use a bracket of form

 +⟨δHδρ,∇⋅δFδu⟩+⟨1ρδHδθ∇θ,δFδu⟩ (2.1.2) −⟨δFδρ,∇⋅δHδu⟩−⟨1ρδFδθ∇θ,δHδu⟩, (2.1.3)

for L2 inner product . denotes a vorticity type variable, given by

 q=(∇×u+2Ω)/ρ, (2.1.4)

for rotation vector

. Note that in the two-dimensional test cases considered in the numerical results section below, corresponds to a scalar vorticity field multiplied by the 2D domain’s outward unit vector, with given by

 q=(∇⊥⋅u+f)/ρ, (2.1.5)

for 2D curl and Coriolis parameter depending on the test case. The functional derivatives appearing in the bracket are defined weakly by

 ⟨δFδu,w⟩\coloneqqlimϵ→01ϵ(F(u+ϵw,D)−F(u,D))∀w∈V(Ω), (2.1.6)

for a suitable space to be defined, and similar for and . In the continuous setting, we consider continuously differentiable spaces for , which are then discretised using suitable finite element spaces. The resulting fluid dynamical equations follow from their respective Hamiltonians and functionals corresponding to weak forms of the prognostic variables, i.e. , , and , respectively, for arbitrary test functions in the corresponding spaces . For , we have

 δFδu=v,δFδρ=0,δFδθ=0, (2.1.7)

so that the Poisson bracket reduces to

 {F,H}=−⟨v,q×δHδu⟩+⟨δHδρ,∇⋅v⟩+⟨1ρδHδθ∇θ,v⟩, (2.1.8)

for any test function . Using the Poisson system (2.1.1), we then obtain a momentum equation given by

 ⟨v,ut⟩=dFdt=−⟨v,q×δHδu⟩+⟨δHδρ,∇⋅v⟩+⟨1ρδHδθ∇θ,v⟩, (2.1.9)

where the subscript in denotes differentiation with respect to . Similarly, for and , we obtain weak versions of the continuity and transport equations of form

 ⟨ψ,ρt⟩=−⟨ψ,∇⋅δHδu⟩, (2.1.10) ⟨χ,θt⟩=−⟨χ,1ρδHδu⋅∇θ⟩, (2.1.11)

for any test functions , . Given a system of equations of form (2.1.9) - (2.1.11), conservation of energy then follows directly from the bracket framework. First, the Poisson system (2.1.1) with bracket (2.1.2) - (2.1.3) can be recovered from (2.1.9) - (2.1.11) for any functional via

 dFdt=⟨δFδu,ut⟩+⟨δFδρ,ρt⟩+⟨δFδθ,θt⟩={F,H}. (2.1.12)

In particular, this also holds for , and noting the bracket’s antisymmetry, we arrive at

 dHdt={H,H}=−{H,H}=0. (2.1.13)

#### 2.1.1 Euler equations

A Poisson bracket based formulation of the dry compressible Euler equations can be found e.g. in [gassmann2013global, lee2020mixed]

. Note that the Poisson bracket and the skew symmetric operator presented in these papers correspond to a different Poisson bracket to the one considered here, relying on the use of a different set of underlying fields (with mass weighted potential temperature

). In our case, the Hamiltonian is given by

 H(u,ρ,θ)=∫Ω(ρ2|u|2+gρz+cvρθπ)dx, (2.1.14)

for wind velocity , air density , and potential temperature . , and denote the vertical coordinate, gravitational acceleration and specific heat of air at constant volume, respectively. Further, denotes the Exner pressure, which is given by the ideal gas law as

 π1−κκ=Rp0ρθ, (2.1.15)

for reference pressure , ideal gas constant , and non-dimensional parameter , where denotes the specific heat at constant pressure. In view of the Poisson system (2.1.1), in the non-discretised case the variational derivatives are given by

 (2.1.16) δHδθ=cpρπ, (2.1.17)

where we used (2.1.15) to derive the pressure related expressions for the variations in and . The usual form of the Euler equations then follows from (2.1.9) - (2.1.11). For the velocity equation, we have

 ⟨v,ut⟩= (2.1.18) =−⟨v,q×(ρu)⟩+⟨12|u|2+gz+cpθπ,∇⋅v⟩+⟨cpπ∇θ,v⟩ (2.1.19) =−⟨v,(∇×u+2Ω)×u⟩−⟨12∇|u|2+gk+cp∇(θπ)⟩+⟨cpπ∇θ,v⟩ (2.1.20) =−⟨v,(u⋅∇)u+2Ω×u+gk+cpθ∇π⟩, (2.1.21)

for vertical unit vector , and where we have used integration by parts (with on ) and . Since this holds for any test function in , we arrive at the usual strong form of the momentum equation given by

 (2.1.22)

Similarly, we obtain the continuity and temperature transport equations of form

 ρt+∇⋅(ρu)=0, (2.1.23) θt+u⋅∇θ=0, (2.1.24)

where we used in (2.1.10) and (2.1.11), which hold for any test functions and , respectively.

#### 2.1.2 Thermal shallow water equations

For the thermal shallow water equations, we follow [ELDRED20191] and consider a Hamiltonian of form

 H(u,ρ,θ)=∫Ω(ρ2|u|2+ρθ(ρ2+b))dx, (2.1.25)

where in this context corresponds to the fluid height (and for consistency of notation, we keep instead of using the more usual or ), and to the buoyancy. Further, denotes the topographic height. The Hamiltonian variations are now given by

 (2.1.26) δHδθ=ρ(ρ2+b), (2.1.27)

and we obtain a weak momentum equation of form

 ⟨v,ut⟩= −⟨δFδu,q×δHδu⟩+⟨δHδρ,∇⋅δFδu⟩+⟨1ρδHδθ∇θ,δFδu⟩ (2.1.28) = −⟨v,q×(ρu)⟩+⟨12|u|2+θ(ρ+b),∇⋅v⟩+⟨(ρ2+b)∇θ,v⟩ (2.1.29) = −⟨v,(∇×u)×u⟩−⟨v,2Ω×u⟩ (2.1.30) (2.1.31) = −⟨v,(u⋅∇)u+(2Ω×u)+θ∇(ρ+b)+ρ2∇θ⟩, (2.1.32)

In the non-discretised case, this leads to

 ut+(u⋅∇)u+fu⊥+θ∇(ρ+b)+ρ2∇θ=0, (2.1.33)

noting that given the spherical context of the shallow water equations, we rewrote the Coriolis term as , for Coriolis parameter and perpendicular , where , , , and correspond to the sphere’s radius, rotation rate, rotational axis coordinate, and unit vertical vector, respectively. Similarly, as in the Euler case, equations for and follow from (2.1.10) and (2.1.11), and since the variational derivative is the same here, the resulting equations are given as before by

 ρt+∇⋅(ρu)=0, (2.1.34) θt+u⋅∇θ=0. (2.1.35)

### 2.2 Space discretisation

In this section, we describe the space discretisation for the Poisson bracket introduced above. First, we discuss our choice of compatible finite element spaces for the fields , , and , as well as the additional space for the thermal field . Next, we describe SUPG stabilised schemes for the advection of and present the incorporation of these schemes into the Poisson bracket framework, which forms the core part of this paper. Finally, we discuss how to include upwinding for and in this bracket, as presented in [WIMMER2019109016].

#### 2.2.1 Choice of compatible finite element spaces

To discretise (2.1.9) - (2.1.11) and the diagnostic vorticity equation (2.1.5), we consider compatible finite element spaces for , , and . They are given by , , and such that

 Vq∇⊥⟶Vu∇⋅⟶Vρ, (2.2.1)

that is the differential operators appearing in our equations map one finite element space to another. In the case of the Euler equations, where we consider a (horizontally periodic) vertical slice of the atmosphere, the spaces are based on tensor products of continuous Galerkin (CG) and discontinuous Galerkin (DG) finite elements defined on intervals and are given by

[natale2016compatible]

 Vq=CGr+1(T)⊗CGs+1(S) (2.2.2) Vu=[DGr(T)k⊗CGs+1(S)]⊕[CGr+1(T)⊗DGs(S)k], (2.2.3) Vρ=DGr⊗DGs(S), (2.2.4)

where , denote horizontal and vertical reference intervals, respectively. In the numerical results section below, we will set and use bilinear () and biquadratic () finite elements for the density space. The corresponding compatible finite elements for are depicted in Figure 1. Note that in three dimensions, we would consider compatible horizontal reference cells instead of intervals in the tensor product. The discretisation described below is still valid in the 3D case, except that the diagnostic scalar vorticity equation would be replaced by its corresponding vector field form.

For the thermal shallow water equations, where we consider triangular spherical meshes, we follow [WIMMER2019109016] and use the degree two Brezzi-Douglas-Marini finite element space for the velocity field. The compatible depth and vorticity spaces are then given by and , i.e. the first and third (polynomial) order discontinuous and continuous triangular Galerkin spaces, respectively.

Finally, we need to specify finite element spaces for the thermal field . In the case of the thermal shallow water equations, we use the continuous space , given by . For the Euler equations, we consider a Charney-Phillips (CP) type space by choosing the finite element nodes to coincide with the vertical wind nodes of the velocity space element. The resulting elements for the thermal field space in the case of the Euler and thermal shallow water equations are given in Figure 2.

In terms of the equations derived by the Poisson bracket in Section 2.1, we obtain a space discretisation given by

 ∀w∈Vu, (2.2.5) ⟨ϕ,ρt⟩+⟨ϕ,∇⋅δHδu⟩=0 ∀ϕ∈Vρ, (2.2.6) ⟨γ,θt⟩+⟨γ,1ρδHδu⋅∇θ⟩=0 ∀γ∈Vθ, (2.2.7)

together with a discretised diagnostic scalar vorticity equation given by

 ⟨η,qρ⟩=−⟨∇⊥η,u⟩+⟨⟨η,n⊥⋅u⟩⟩+⟨η,f⟩ ∀η∈Vq, (2.2.8)

for L2 inner product on . Note that the variational derivatives of are now given by projections into the relevant finite element spaces, e.g. for the Hamiltonian corresponding to the Euler equations, we find

 δHδu=PVu(ρu),δHδρ=PVρ(12|u|2+gz+cpθπ), (2.2.9) δHδθ=PVθ(ρcpπ), (2.2.10)

where denotes projection into . Further, we note that the Poisson bracket structure still holds in the discretised case, i.e. the discretised evolution equations (2.2.5) - (2.2.7) can be used for evaluating as in (2.1.12), leading to a Poisson system (as in (2.1.1)) with a bracket whose form is equal to the continuous one, given by (2.1.2) - (2.1.3). While it is not clear if the discretised bracket still satisfies the Jacobi identity, it is still antisymmetric, implying that energy is conserved. We therefore find that the discretised bracket is at least an almost Poisson bracket, i.e. a bilinear antisymmetric form.

### 2.3 SUPG stabilisation

To review the SUPG advection scheme, we consider the thermal field transport equation in its continuous form, given by

 θt+u⋅∇θ=0. (2.3.1)

In a finite element setting, we would then space discretise using its weak form

 ⟨γ,θt⟩+⟨γ,u⋅∇θ⟩=0 ∀γ∈Vθ. (2.3.2)

However, for standard continuous Galerkin spaces, the advection term tends to produce spurious oscillations, which can be counteracted by introducing diffusion in the direction of the flow [kuzmin2010guide]. This is achieved by replacing the test functions by

 γ→γ+τu⋅∇γ\eqqcolonγu, (2.3.3)

for a suitable coefficient . We then arrive at

 ⟨γu,θt+u⋅∇θ⟩=0 ∀γ∈Vθ, (2.3.4)

and this form of upwinding is commonly referred to as the Streamline Upwind Petrov Galerkin (SUPG) method [kuzmin2010guide], and was first introduced in [brooks1982streamline]. We find that the modified test function of the advection term acts as diffusion along the velocity field , noting that for , the SUPG equation (2.3.4) leads to

 12∂∂t∥θ∥22=−⟨θ,u⋅∇θ⟩−∥√τu⋅∇θ∥22−τ⟨u⋅∇θ,θt⟩, (2.3.5)

where the indefinite SUPG term (i.e. the last term) is kept small by a suitable choice of . Note that since all test functions where modified to include an upwind term, the resulting weak equation (2.3.4) is still clearly consistent, in the sense that an exact solution to the continuous equation (2.3.1) also solves the weak equation.

The upwinded advection equation (2.3.4) can be used for fully continuous spaces , as is the case for our choice of buoyancy space for the thermal shallow water equations. For spaces that are not fully continuous, such as the Charney-Phillips type finite element space defined above, further care needs to be taken for the discontinuous components. In our case, we follow [natale2016compatible], and apply the standard DG upwinding scheme [kuzmin2010guide] in the horizontal (i.e. discontinuous) direction, and restrict the SUPG setup to the vertical direction. To this end, we first integrate the weak form (2.3.2) by parts and arrive at

 ⟨γ,θt⟩−⟨∇⋅(uγ),θ⟩+∫Γv[[uγ]]~θdS=0 ∀γ∈Vθ, (2.3.6)

where denotes the set of all vertical interior facets of the underlying mesh. The jump operator (for vectors and scalars , respectively) and upwind value are defined by

 [[v]]=v+⋅n++v−⋅n−,[[ψ]]=ψ+−ψ−,~θ={θ+% if u⋅n+<0,θ−otherwise, (2.3.7)

noting that the two sides of each mesh facet are arbitrarily denoted by + and - (and hence ). Before incorporating the upwind modification (2.3.3) to the test functions in (2.3.6), we need to integrate by parts again to avoid applying the differential operator to the upwinded part , as the double differentiation may not be well-defined for . Further, we restrict the SUPG scheme to the vertical direction by using modified test functions of form

 γ→γ+τ(k⋅u)(k⋅∇γ)\eqqcolonγ(k⋅u)k, (2.3.8)

for vertical unit vector . The resulting, fully upwinded advection equation is then given by

 ⟨γ(k⋅u)k,θt⟩+⟨γ(k⋅u)k,u⋅∇θ⟩+∫Γv([[uγ(k⋅u)k]]~θ−[[uγ(k⋅u)kθ]])dS=0 ∀γ∈Vθ. (2.3.9)

In the next section, we will incorporate the stabilised advection schemes (2.3.4) and (2.3.9) in the energy conserving bracket (2.1.2) - (2.1.3). To simplify notation, we write the upwinded test functions (2.3.3) and (2.3.8) as well as advection terms (2.3.4) and (2.3.9) as operators given in the definition below.

###### Definition 1

The SUPG contribution to is given by

 ~γu\coloneqq{τu⋅∇γCG type Vθ,τ(k⋅u)(k⋅∇γ)CP type Vθ. (2.3.10)

Further, the discrete thermal field transport operator is given by

 (2.3.11)

Using this definition, we can rewrite the advection equations (2.3.4) and (2.3.9) as a single equation of form

 ⟨γ+~γu,θt⟩+Lθu(θ;γu)=0 ∀γ∈Vθ. (2.3.12)

Finally, we note that the lowered at corresponds to the advecting velocity, while the lowered at corresponds to the velocity used for upwinding. Depending on the time discretisation, we will find that these may be distinct in the almost Poisson bracket setup.

### 2.4 SUPG upwinding in the discretised bracket

It remains to incorporate SUPG upwinding, as described in Section 2.3, into the discretised almost Poisson bracket. Without upwinding, the bracket is given by (2.1.2) - (2.1.3), such that , , , and the variations are functions of the corresponding finite element spaces. We aim to replace the corresponding thermal field advection equation (2.2.7) by the SUPG version (2.3.12), i.e.

 ⟨γ,θt⟩+⟨γ,1ρδHδu⋅∇θ⟩=0→⟨γ+~γu,θt⟩+LθδHδu/ρ(θ;γu)=0, (2.4.1)

where as for the non-upwinded version, we anticipate an advecting velocity of form when upwinding is included. Note that in the Poisson system (2.1.1), the time derivative term corresponds to , while the advection term is derived from the bracket . Considering the upwinded equation, we find that choosing

 F=⟨γ+~γu,θ⟩ (2.4.2)

in the Poisson system is impractical, since then the time derivative of will include a time derivative in the upwinding velocity . To avoid this, we move the corresponding upwinding part to the right-hand side:

 ⟨γ,θt⟩=−LθδHδρ/u(θ;γu)−⟨~γu,θt⟩, (2.4.3)

that is we will consider the upwind part of the time derivative term as part of the Poisson bracket. Including this term in the bracket, however, leads to further complications. First, a corresponding antisymmetric term needs to be added, which will appear in the momentum equation’s forcing part. Second, in view of the energy conserving time integration scheme to be used below, the bracket should be independent of time derivatives. This means that in the bracket should be replaced by another, time-derivative free term. Considering the upwinded equation (2.4.1), we find that depends on the advection operator as well as the upwinded test function contribution , motivating Definition 2 below.

###### Definition 2

The thermal field advection operator is given by

 ⟨γ,θa(w)⟩=−Lθw/ρ(θ;γu) ∀γ∈Vθ. (2.4.4)

Further, the SUPG recovery operator is given by

 ⟨γ,Su(τ)⟩=⟨γ+~γu,τ⟩∀γ∈Vθ. (2.4.5)

Note that both and implicitly also depend on the state . Using these two operators, we can rewrite the upwinded thermal field advection equation (2.4.1) as

 ⟨γ,Su(θt)⟩=⟨γ,θa(δHδu)⟩ ∀γ∈Vθ, (2.4.6)

for any test function . Since both and map into , this implies point-wise, and therefore , provided exists. This motivates a modified, upwinded variation part of the bracket according to

 +⟨1ρδHδθ∇θ,δFδu⟩ +LθδFδu/ρ(θ;(δHδθ)u)+⟨(˜δHδθ)u,S−1u(θa(δFδu))⟩ (2.4.7) → −⟨1ρδFδθ∇θ,δHδu⟩ −LθδHδu/ρ(θ;(δFδθ)u)−⟨(˜δFδθ)u,S−1u(θa(δHδu))⟩, (2.4.8)

where (2.4.8) corresponds to the thermal field advection equation, while (2.4.7) does to the corresponding antisymmetric bracket term, which appears as a forcing term in the momentum equation. Note that is bilinear in and (considering its definition (2.3.11)) and further that since both and are linear, the corresponding upwind contribution bracket term is also bilinear. Since (2.4.7) and (2.4.8) form an antisymmetric pair by construction, we therefore find that the modified, upwinded bracket given in Definition 3 below is still an antisymmetric bilinear form. In Proposition 1, we confirm that this setup does indeed lead to an SUPG equation of form (2.4.1).

###### Definition 3

The thermal field upwinded almost Poisson bracket is given by

 −⟨δFδu,q×δHδu⟩ (2.4.9) (2.4.10)

with operators , and given by Definitions 1 and 2.

###### Proposition 1

Consider the space discretised almost Poisson bracket given by Definition 3. Then the thermal field transport equation based on the Poisson system (2.1.1) with this bracket is given by

 ⟨γ+~γu,θt⟩+LθδHδu/ρ(θ;γu)=0 ∀γ∈Vθ. (2.4.11)

Proof. We test the upwinded bracket using the usual functional of form . In this case, as before all parts of the bracket except for the ones corresponding to temperature advection (i.e. (2.4.8)) are zero, and we are left with

 dFdt=⟨γ,θt⟩=−LθδHδu/ρ(θ;γu)−⟨~γu,S−1u(θa(δHδu))⟩. (2.4.12)

By definition of , the advection operator part of the right-hand side is equal to , which in turn can be expressed using as

 ⟨γ+~γu,S−1u(θa(δHδu))⟩. (2.4.13)

We can then cancel with the upwind contribution term (i.e. the last term in (2.4.12)), so that

 (2.4.14)

which holds for any . Hence, the term is identically equal to , and we may replace it in (2.4.12). Moving the upwind contribution term and the advection operator to the left-hand side then gives the required upwinded form.

The second implication of including a term in the bracket lies in the corresponding antisymmetric term (i.e. the second term in (2.4.7)), given by

 ⟨(˜δHδθ)u,S−1u(θa(δFδu))⟩. (2.4.15)

The test function appears as an argument of , which we would be required to solve for for every test function in the variational solver for the velocity. To avoid this, we use the transpose of , allowing us to move the operator away from the test function.

###### Proposition 2

Consider the space discretised almost Poisson bracket given by Definition 3. Then the momentum equation based on the Poisson system (2.1.1) with this bracket is given by

 ∀w∈Vu, (2.4.16)

where is given by

 ⟨s,v⟩=−Lθv/ρ(θ;σu) ∀v∈Vu, (2.4.17)

for such that

 ⟨η,σu⟩=⟨η,(˜δHδθ)u⟩ ∀η∈Vθ. (2.4.18)

Proof. The first three terms (excluding the time derivative term) follow directly from the bracket, using in the Poisson system. To obtain the last term, we evaluate the the corresponding bracket term (2.4.15) as

 ⟨(θa)∗(S−1u)∗((˜δHδθ)u),δFδu⟩, (2.4.19)

for transpose operator . First, for , we find directly that its transpose is given by

 ⟨(θa)∗(γ),v⟩=−Lθv/ρ(θ;γu) ∀v∈Vu, (2.4.20)

while for the upwind operator we note that , with clearly following from the definition of as

 S∗u(γ)=PVθ(γ+~γu). (2.4.21)

To find the inverse, we search for a such that, given “an upwinded function” ,

 (2.4.22)

Finally, we solve this for with

 χ=(˜δHδθ)u, (2.4.23)

which yields as defined above. then follows applying (2.4.20) on , as required.

Using Propositions 1 and 2, we obtain from the -upwinded discretised bracket in Definition 3 a system of equations of form

 ∀w∈Vu, (2.4.24) ⟨ϕ,ρt⟩+⟨ϕ,∇⋅δHδu⟩=0 ∀ϕ∈Vρ, (2.4.25) ⟨γ+~γu,θt⟩+LθδHδu/ρ(θ;γu)=0 ∀γ∈Vθ. (2.4.26)

Finally, in a similar fashion to the proofs above, we find that the Poisson system can be recovered from these equations from the chain rule for

as in (2.1.12). Since the -upwinded discretised bracket is antisymmetric, we therefore find that equations (2.4.24) - (2.4.26) are still energy conserving.

###### Remark 1

To also incorporate energy conserving upwinding schemes for and in the equation set (2.4.24) - (2.4.26), we can use the framework provided in [WIMMER2019109016], where a bracket for the rotating shallow water equations was considered. The bracket is identical to the one introduced above, up to excluding the two thermal terms (i.e. setting in (2.4.9) and in (2.4.10)), and we can directly apply the upwinding schemes here. The shallow water upwinded bracket is given by

 {F,H}\coloneqq (2.4.27) (2.4.28) (2.4.29)

where (2.4.27) corresponds to velocity advection, (2.4.28) to the forcing terms in the momentum equation, and (2.4.29) to density advection. For the purpose of density upwinding, a velocity recovery operator of form

 (2.4.30)

was introduced (for details, see [WIMMER2019109016]). Note that corresponds to a discrete division by , and recalling that is given by the discrete flux , we find that corresponds to the advecting velocity for and advection. The bracket then leads to an upwinded density equation of form

 ⟨ϕ,ρt⟩=⟨ρU(ρ,δHδu),∇ϕ⟩−∫Γ[[ϕU(ρ,δHδu)]]~ρdS ∀ϕ∈Vρ, (2.4.31)

which corresponds to the standard DG upwind scheme. For the velocity equation, we obtain

 ⟨w,ut⟩ (2.4.32) −⟨ρU(ρ,w),fU(ρ,δHδu)⊥⟩ (2.4.33) −⟨ρU(ρ,w),∇δHδρ⟩+∫Γ[[δHδρU(ρ,w)]]~ρdS∀w∈Vu, (2.4.34)

where upwinding was applied to the rotational velocity advection part, given by . Note that the velocity recovery operator is applied to all test functions on the right hand side of the velocity equation. For a consistent use of test functions, when extending this bracket to the thermal case, we then also apply this operator to the test functions of the corresponding thermal bracket term (2.4.7). To maintain antisymmetry, we then also need to apply it in (2.4.8), and we arrive at a fully upwinded bracket given by (2.4.27) - (2.4.29), together with

 (2.4.35) (2.4.36)

and now defined by

 ⟨γ,θa(v)⟩=−Lθv(θ;γu) ∀γ∈Vθ. (2.4.37)

The flux recovered velocity then also serves as the advecting velocity for the thermal field transport equation, given by

 ⟨γ+~γu,θt⟩+LθU(ρ,δHδu)(θ;γu) ∀γ∈Vθ. (2.4.38)

Finally, we note that the underlying 2D velocity advection scheme and Coriolis term in (2.4.27) and (2.4.28) above can be extended readily to the three-dimensional case [natale2016compatible], which in our case leads to

 (2.4.39) −⟨ρU(ρ,δFδu),2Ω×U(ρ,δHδu)⟩, (2.4.40)

where

 {{n×w}}=n+×w++n−×w−. (2.4.41)

## 3 Numerical results

In this section, we confirm numerically the upwind stabilised and energy conserving properties of the newly introduced almost Poisson bracket (2.4.9) - (2.4.10). First, we review the energy conserving time discretisation used for this purpose, as well as the resulting non-linear set of equations. Next, we discuss the operational applicability of the bracket, including simplifications in view of computational costs and the extension to a fully upwinded bracket. Finally, we present test cases, consisting of a perturbed thermogeostrophic balance scenario for the thermal shallow water equations, as well as a cold and hot air bubble scenario for the Euler equations.

### 3.1 Time discretisation and nonlinear system

To confirm the energy conserving property of the space discretised equations (2.4.24) - (2.4.26), we apply an energy conserving time discretisation, thus expecting energy conservation up to machine precision. It is given by a Poisson integrator as introduced in [cohen2011linear], and can be applied to the framework used here as detailed in [WIMMER2019109016]. For the prognostic fields , we have

 (3.1.1)

where the skew symmetric transformation is related to the almost Poisson bracket via

 {F,H}=⟨δFδz,J(z)δHδz⟩, (3.1.2)

and the time averaged Hamiltonians are given by

 ¯¯¯¯¯¯¯¯¯¯δHδu\coloneqq∫10δδuH(zn+s(zn+1−zn))ds, (3.1.3)

and similarly for the variations in and . The expressions can be integrated exactly for the thermal shallow water Hamiltonian, leading to

 (3.1.4)

and similarly for the other two variations. However, for the Euler Hamiltonian, we find that the internal energy is a non-polynomial function in , , thus requiring an approximate integration of (3.1.3) for the variations in and . For the numerical results below, a fourth order Gaussian quadrature was used.

The resulting nonlinear system of equations is given by

 ⟨w,un+1−un⟩+Δt (⟨w,¯q×¯¯¯¯¯¯¯¯¯¯δHδu⟩−⟨¯¯¯¯¯¯¯¯¯¯δHδρ,∇⋅w⟩ (3.1.5) ∀w∈Vu, (3.1.6) ⟨ϕ,ρn+1−ρn⟩+Δt ⟨ϕ,∇⋅¯¯¯¯¯¯¯¯¯¯δHδu⟩=0 ∀ϕ∈Vρ, (3.1.7) ⟨γ¯u,θn+1−θn⟩+Δt Lθ¯¯¯¯¯¯¯δHδu/¯ρ(¯θ;γ¯u)=0 ∀γ∈Vθ, (3.1.8)

for midpoint time average , and similar for and , with and given by the diagnostic scalar vorticity relation (2.2.8), using the prognostic fields at time levels and , respectively. Note that all occurrences of the state in and are also discretised at the half-time level. Further, we remark that the derivation of the fully discretised thermal field equation (3.1.8) follows the proof of Proposition 1 in a time discretised form, with replaced by .

To solve for , given the fully discretised residual above (left-hand sides of (3.1.5) - (3.1.8)), we apply the same procedure as detailed in [WIMMER2019109016] and revert to a Picard iteration scheme. For update

with unknown next time step estimate

, we set

 δR′δz(δz)=−R(zn+1,k), (3.1.9)

and . The left hand side corresponds to the Jacobian of a linearised version of without Hamiltonian projections, and is given by

 δR′δz(δz)=⎛⎜ ⎜ ⎜⎝[l]⟨δu,w⟩+Δt2(2Ω×δu,w⟩−⟨gδρ+cp(¯¯θδπ+δθ¯¯π),∇⋅w⟩+cp⟨¯¯π∇δθ,w⟩)⟨δρ,ϕ⟩+Δt2⟨¯¯ρ∇⋅δu,ϕ⟩⟨δθ,γ⟩⎞⎟ ⎟ ⎟⎠, (3.1.10)

where double-barred entries correspond to background fields, with , , and . Note that in view of the test cases to follow, we assumed the background thermal field to be constant (equal to for the thermal shallow water equations, and an isentropic background potential temperature for the Euler equations). This leads to a vanishing term in the velocity and thermal field equation, respectively, thus uncoupling the latter from the density and velocity equations. Given the right-hand side we can then first solve for , followed by a mixed solve for the density and velocity updates. Note that since this paper focuses on the underlying space discretisation, the nonlinear solve procedure is kept simple. The discretisation can in principle be applied equally to fully three-dimensional problems and scenarios with a non-constant potential background temperature (by approximately eliminating in the linearised momentum equation). A suitable nonlinear solver strategy for this case can be found in [bendall2019compatible].

The mesh, finite element discretisation, and solver were implemented using the automated finite element toolkit Firedrake111for further details, see [bercea2016structure, Gibson:aa, homolya2018tsfc, Luporini2017] or http://firedrakeproject.org [rathgeber2016firedrake], with a hybridised solver (see e.g. [shipton2018higher]) used to solve for the mixed system in .

### 3.2 Fully upwinded, approximately energy conserving formulation

While both the space and time discretisations conserve energy, the resulting Picard iteration scheme does not. In the numerical tests below, we find that the number of Picard iterations required to achieve energy conservation up to machine precision is much higher than is usually considered in forecasting models. This suggests the use of simplifications to the fully discretised scheme, which lead to an additional energy error smaller than the error due to a small number of Picard iterations. Further, the temperature upwinded bracket presented in Definition 3 does not include upwinding in the velocity and density fields, leading to a non-satisfactory quality of the field development. In this section, we present a fully upwinded fully discretised formulation based on the bracket as described in Remark 1 and a small modification to the Poisson time discretisation.

In the fully upwinded bracket, we find that in the temperature advection equation the advecting and upwind velocities are given by and , respectively, which time-discretise within the Poisson integrator as (given by (3.1.4)) and , respectively. The equations derived from the bracket can be simplified by ensuring that these two velocities are equal after time-discretisation, which can be achieved by time discretising the Hamiltonian variation in as

 δHδu=PVu(ρu)→PVu(¯ρ¯u)\eqqcolon¯¯¯¯¯¯¯¯¯¯δHδu′. (3.2.1)

In this case, we find that the velocity recovery operator is evaluated as

 ⟨¯ρv,U(¯ρ,¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯¯δHδu′)⟩=⟨v,PVu(¯ρ¯u)⟩=⟨v,¯ρ¯u⟩ ∀v∈Vu, (3.2.2)

so that

 U(¯ρ,¯¯¯¯¯¯¯¯¯¯δHδu′)