PaperlayerView exact revision on arXiv

Source-linked AI summary

Mixed finite elements for numerical weather prediction

C. J. Cotter, J. Shipton

arXiv:1103.2440v3math.NAcs.NA

TL;DR

Mixed finite element discretisations offer a way to use flexible, potentially non-orthogonal grids while addressing spurious wave branches associated with degree-of-freedom ratios. The paper applies these spaces to linear rotating shallow-water equations and shows exact steady geostrophic modes on arbitrary meshes. Its supported scope also includes energy and local mass conservation properties of the formulation.

  • Problem

    Triangular and hexagonal C-grids alter the velocity-to-pressure degree-of-freedom ratio, which can introduce spurious inertia-gravity or Rossby modes.

    Source: §1 · p. 2

  • Method

    The paper constructs mixed finite element discretisations of the linear rotating shallow-water equations using spaces with div-curl structure.

    Source: §2.1 · p. 3

  • Results

    Geostrophic modes are exactly steady for the paper’s linear rotating shallow-water discretisations on arbitrary meshes.

    Source: §1 · p. 2

  • Takeaways & Limitations

    Mixed finite elements provide a flexible extension of C-grid methods that can use non-orthogonal grids and adjust the pressure-to-velocity degree-of-freedom ratio to avoid spurious mode branches.

    Source: §5 · p. 21

    Benchmarking the augmented mixed element pair for numerical weather prediction and ocean modelling remains future work.

    Source: §5 · p. 21

Abstract

from arXiv · show

We show how two-dimensional mixed finite element methods that satisfy the conditions of finite element exterior calculus can be used for the horizontal discretisation of dynamical cores for numerical weather prediction on pseudo-uniform grids. This family of mixed finite element methods can be thought of in the numerical weather prediction context as a generalisation of the popular polygonal C-grid finite difference methods. There are a few major advantages: the mixed finite element methods do not require an orthogonal grid, and they allow a degree of flexibility that can be exploited to ensure an appropriate ratio between the velocity and pressure degrees of freedom so as to avoid spurious mode branches in the numerical dispersion relation. These methods preserve several properties of the C-grid method when applied to linear barotropic wave propagation, namely: a) energy conservation, b) mass conservation, c) no spurious pressure modes, and d) steady geostrophic modes on the $f$-plane. We explain how these properties are preserved, and describe two examples that can be used on pseudo-uniform grids: the recently-developed modified RT0-Q0 element pair on quadrilaterals and the BDFM1-\pdg element pair on triangles. All of these mixed finite element methods have an exact 2:1 ratio of velocity degrees of freedom to pressure degrees of freedom. Finally we illustrate the properties with some numerical examples.

1. Introduction

Pseudo-uniform grids support direct addressing while avoiding the polar singularity, but triangular and hexagonal C-grids alter the velocity-to-pressure degree-of-freedom ratio and can support spurious modes. The paper uses mixed finite element spaces to retain key C-grid properties while choosing degrees of freedom to avoid spurious branches.

1. Introduction

  • Pseudo-uniform spherical grids avoid the polar singularity and provide direct addressing, while C-grid staggering prevents non-constant pressure modes with zero numerical gradient.Discrete curl and divergence operators satisfying div curl = 0 can also yield stationary geostrophic modes on the f-plane, a necessary condition for accurate geostrophic adjustment.Source: §1 · p. 1
  • Triangular and hexagonal C-grids depart from the quadrilateral grid’s 2:1 velocity-to-pressure DOF ratio, producing extra inertia-gravity or Rossby modes.On triangular grids, the extra inertia-gravity branch can produce checkerboard divergence patterns at small Rossby deformation radius; hexagonal grids support an extra low-frequency Rossby mode.Source: §1 · p. 2
  • Finite element spaces offer flexibility to adjust triangular-element degrees of freedom, while the BDM1 pair has 3 velocity DOF and 1 pressure DOF per triangle.The BDM1 ratio is associated with too many Rossby modes per inertia-gravity mode, as in the hexagonal C-grid.Source: §1 · p. 2
  • The paper shows that discrete Helmholtz decomposition makes geostrophic modes exactly steady for linear rotating shallow-water discretisations on arbitrary meshes.The decomposition requires a commuting diagram, and the paper identifies dim(S) = 2 dim(V) as an optimal choice satisfying necessary conditions against both spurious wave types in the periodic plane.Source: §1 · pp. 2–3

2. Mixed finite elements for geophysical fluid dynamics

The paper applies mixed finite elements to the linear rotating shallow-water equations as a simplified setting for studying horizontal discretisation and wave propagation. The formulation selects velocity and elevation spaces with div-curl structure, supporting the properties required for numerical weather prediction.

2. Mixed finite elements for geophysical fluid dynamics

  • The linear rotating shallow-water equations provide a simplified model for horizontal dynamical-core discretisation, with attention restricted to wave propagation on the f-plane, β-plane, or sphere.The formulation targets global energy and local mass conservation, absence of spurious pressure modes, and steady geostrophic states.Source: §2 · p. 3

2.1. Spatial discretisation of the linear rotating shallow-water equations

  • The weak finite element formulation uses velocity space S ⊂ H(div) and elevation space V ⊂ L2 after integrating the pressure-gradient term by parts.The equations describe horizontal velocity and layer elevation in the linear rotating shallow-water system.Source: §2.1 · p. 3
  • The divergence maps S onto V, while a streamfunction space E ⊂ H1 maps by curl onto the divergence-free kernel in S.This structure makes E functions continuous, S vector fields continuous in their normal components, and V functions discontinuous.Source: §2.1 · p. 4

2.2. Energy conservation

For the linearised equations, the mixed finite element Galerkin formulation conserves global energy. The paper motivates this property by its role in avoiding spurious sources of unbalanced fast waves and as a precursor to energy-conserving nonlinear discretisation.

2.2. Energy conservation

  • Global energy conservation follows immediately from the mixed finite element Galerkin formulation for the linearised equations.The derivation substitutes the discrete solutions and chooses the velocity and elevation solutions as test functions.Source: §2.2 · p. 4
  • Energy conservation helps prevent numerical sources of unbalanced fast waves in numerical weather prediction models.The paper also identifies it as a precursor to energy-conserving discretisation of the nonlinear equations in vector-invariant form.Source: §2.2 · p. 4
  • The conserved energy expression includes velocity and elevation contributions, with the elevation term weighted by c2.The expression is given for the linearised equations.Source: §2.2 · p. 4

2.3. Local mass conservation

Local mass conservation follows from including element indicators in the discontinuous pressure space and using them as test functions. Continuous normal velocity components then ensure continuous element-boundary fluxes.

2.3. Local mass conservation

  • Local mass conservation is obtained by testing with each element indicator, which belongs to the discontinuous elevation space V.The paper identifies local mass conservation as preventing spurious mass sources and sinks.Source: §2.3 · p. 4
  • Using an element indicator as the test function gives a local mass balance for that element.This follows from the finite element equation for the linear continuity equation.Source: §2.3 · p. 4
  • Continuous normal velocity components make the flux continuous across element boundaries, so elevation is locally conserved.Source: §2.3 · p. 4

2.4. Absence of spurious pressure modes and stability of discrete Poisson equation

The section defines spurious pressure modes and explains how mixed finite element stability conditions rule them out while supporting a stable discrete Poisson problem.

2.4. Absence of spurious pressure modes and stability of discrete Poisson equation

  • Centered differences on the collocated A-grid admit an oscillatory checkerboard pressure mode with zero numerical gradient, which can pollute solutions under forcing, boundaries, or nonlinearities.The mixed finite element framework instead characterizes spurious modes through their discrete gradients and a mesh-independent stability bound.Source: §2.4 · p. 5
  • A spurious pressure mode has zero discrete gradient on uniform grids, or a gradient tending to zero while its true gradient stays bounded on varying meshes, threatening optimal convergence.The space is defined as free of such modes when a mesh-independent γ2 > 0 bound holds for every pressure function.Source: §2.4 · p. 5
  • The mixed Poisson discretisation is stable when condition (5) and a second condition with mesh-independent γ1 both hold; condition (6) follows when divergence maps into V.If divergence is orthogonal to V, it must then vanish.Source: §2.4 · pp. 5–6
  • A bounded projection ΠS that commutes with divergence and the L2 projection ΠV provides the sufficient condition (5) used in the discrete Poisson stability argument.The projection uses L2 projections of normal components on element edges and orthogonality conditions within each element.Source: §2.4 · pp. 5–6

2.5. Discrete Helmholtz decomposition

The discrete Helmholtz decomposition splits each discrete velocity into curl, discrete-gradient, and harmonic components. Its dimension relations connect the velocity, pressure, and streamfunction spaces.

2.5. Discrete Helmholtz decomposition

  • Because condition (5) gives the discrete gradient no nontrivial kernel, the curl and discrete-gradient images form orthogonal subspaces of S.Together with the discrete harmonic space H, they give uh = ∇⊥ψh + Dφh + hh.Source: §2.5 · p. 6
  • The discrete harmonic space Hh has the same dimension as the continuous harmonic-velocity space H, although the spaces need not coincide.On the periodic plane these fields are constant velocities; the sphere has none.Source: §2.5 · pp. 6–7
  • The decomposition gives dim(E) = dim(V) on the periodic plane and dim(E) = dim(V) + 2 on the sphere, relating velocity-space size to wave-propagation behavior.The section notes that whether dim(S) exceeds 2 dim(V) becomes important when examining wave propagation.Source: §2.5 · p. 7

2.6. Vorticity and divergence

The section describes recovering velocity from discrete vorticity and divergence through elliptic problems for a streamfunction and velocity potential.

2.6. Vorticity and divergence

  • Given discrete vorticity and divergence, velocity can be recovered by solving separate elliptic problems for the streamfunction and velocity potential.The streamfunction problem follows by rewriting the Helmholtz decomposition as a finite element Poisson equation.Source: §2.6 · p. 7
  • The streamfunction is obtained with the usual finite element discretisation of the Poisson equation.This solve recovers the streamfunction from the discrete vorticity.Source: §2.6 · p. 7
  • The coupled system for the velocity potential is well-posed when conditions (5) and (6) hold.It is the mixed finite element approximation to the Poisson equation.Source: §2.6 · p. 7

2.7. Steady geostrophic modes

On the f-plane, geostrophic balanced states are steady in the linear rotating shallow-water equations. The section shows this by initializing velocity from a streamfunction and height from geostrophic balance.

2.7. Steady geostrophic modes

  • Geostrophic balanced states on the f-plane are steady because their velocity has zero divergence.In the quasi-geostrophic limit, such states become slowly evolving balanced states that characterize large-scale weather systems.Source: §2.7 · p. 7
  • The discrete balanced state is initialized with uh = ∇⊥ψh and ηh set from the geostrophic balance relation.The proof checks that vanishing divergence yields a steady height tendency and then a steady velocity.Source: §2.7 · p. 8
  • The resulting height tendency is zero, so the initialized discrete geostrophic balanced state is steady.The argument uses the divergence-free velocity uh = ∇⊥ψh.Source: §2.7 · p. 8

2.8. Numerical dispersion relations

The discrete Helmholtz decomposition frames the dispersion analysis on the f-plane and the β-plane in the quasi-geostrophic limit. Projection kernels can produce spurious inertia-gravity or Rossby branches, motivating dimension conditions on the finite element spaces.

2.8. Numerical dispersion relations

  • The first row of the f-plane discrete eigenvalue problem converges to the continuous problem’s eigenvalues at the optimal rate.The row discretizes the rotating shallow-water eigenvalue problem using mixed finite element spaces V and S.Source: §2.8 · p. 9
  • dim(V) ≤ dim(E) is necessary to avoid spurious divergent inertia-gravity branches because the projection from V to E and back can have a kernel.This occurs when V has larger dimension than E, as for the lowest-order Raviart–Thomas element on triangles.Source: §2.8 · p. 10
  • In the β-plane quasi-geostrophic limit, the leading-order solution matches the f-plane steady geostrophic-state formula.The Rossby-wave eigenvalue problem’s continuous finite element approximation in E has convergent eigenvalues.Source: §2.8 · p. 10
  • dim(V) = dim(E) is necessary to avoid both spurious divergent and spurious irrotational modes, but is not sufficient.A non-trivial kernel in either projection can still produce spurious branches.Source: §2.8 · p. 11
  • When P_E P_V has a non-trivial kernel, the BDM1 triangular pair can produce spurious Rossby-wave branches; dim(V) = dim(E) is therefore necessary to avoid both branch types.For BDM1 on triangles, the passage gives 2 dim(V) = dim(E) + 2 − dim(H).Source: §2.8 · p. 11

3. Examples

The examples pair mixed finite element spaces with pseudo-uniform grids while maintaining twice as many velocity as pressure degrees of freedom, preventing spurious mode branches. They illustrate quadrilateral and triangular constructions, including a divergence correction for general quadrilaterals and an augmented triangular velocity space.

3. Examples

  • Both examples have twice as many velocity as pressure degrees of freedom, a ratio that prevents spurious mode branches.The examples are the modified Raviart-Thomas element on quadrilaterals and the BDFM element on triangles.Source: §3 · p. 11

3.1. Modified Raviart-Thomas element on quadrilaterals

  • For general quadrilaterals, a divergent basis correction makes the velocity divergence constant without changing its streamfunction embedding.The correction has vanishing normal components on the element boundary; the streamfunction space is continuous bilinear Q1.Source: §3.1 · p. 12
  • The modified RT0-Q0 quadrilateral pair has one pressure degree of freedom per element and one velocity degree of freedom per edge, yielding the required 2:1 ratio.With periodic boundaries or on the sphere, each edge is shared by two elements.Source: §3.1 · p. 12

3.2. Brezzi-Douglas-Fortin-Marini element on triangles

  • On triangles, BDFM1 paired with discontinuous P1 pressure satisfies 2 dim(V) = dim(S), unlike RT0-P0 and BDM1.The extra velocity degrees of freedom are useful because they give dim(V) = dim(E); RT and BDM use fewer degrees of freedom for the same order of convergence after suitable post-processing.Source: §3.2 · p. 12
  • The augmented triangular pair has six velocity and three pressure degrees of freedom per element, with no pressure continuity requirement across edges.Its streamfunction space uses continuous quadratic functions augmented by cubic bubbles, whose curls map into the velocity space.Source: §3.2 · p. 14

4. Numerical results

The numerical tests examine steady geostrophic states, Kelvin and Rossby wave propagation, and solid rotation on the sphere. They report steady states to machine precision in f-plane tests, no visible unbalanced motions in wave tests, and near-round-off error for spherical solid rotation.

4. Numerical results

  • BDFM1 tests used the implicit midpoint rule and hybridisation, reducing the solve to a symmetric positive definite system for Lagrange multipliers.Velocity and free-surface variables were reconstructed element-by-element; the approach permits fully implicit treatment of the Coriolis term.Source: §4 · p. 16

4.1. Steady states for the f-plane

  • The f-plane steady tests kept layer thickness and velocity constant to machine precision, both on an unstructured planar mesh and an icosahedral sphere mesh.The tests initialized velocity from a streamfunction and obtained the free surface by solving a discrete elliptic system.Source: §4.1 · p. 16

4.2. Kelvin waves in a circular basin

  • Kelvin-wave tests showed no spurious gravity waves, indicating that BDFM1 maintained normal geostrophic balance and the Kelvin-wave structure.The circular-basin test used Ro = 0.1, Fr = 1, and ran to 10 dimensionless time units.Source: §4.2 · p. 16

4.3. Rossby waves

  • Rossby-wave errors converged at O(∆x^3) until finite-Rossby-number error caused saturation.The exact Rossby-wave solution is only an asymptotic solution of the linearised rotating shallow-water equations, with O(Ro^2) error; third-order convergence is attributed to the continuous piecewise quadratic functions in E, subject to a perturbation tending to zero, for which no estimates are available.Source: §4.3 · p. 18
  • On the cylindrical β-tube mesh, no unbalanced motions were visible in the Rossby-wave plots.The test used a cylinder of unit dimensionless radius and height 2 with a spatially varying Coriolis parameter.Source: §4.3 · p. 19

4.4. Solid rotation on the sphere

  • On a level 4 icosahedral sphere mesh (edges subdivided by 8), after 10 days with a 3600s timestep, free-surface elevation errors from initial condition were almost indistinguishable from round-off error.The mapping yields ψh = z and fh = 2|Ω|z, both exactly representable in their respective finite-element spaces.Source: §4.4 · pp. 19–20

5. Summary and outlook

Mixed finite elements offer flexible discretisations for pseudo-uniform grids, including non-orthogonal grids and adjustable degree-of-freedom ratios to avoid spurious mode branches. The paper outlines low-order examples and future work on benchmarking and higher-order elements.

5. Summary and outlook

  • Mixed finite element spaces extend mimetic C-grid discretisations to non-orthogonal grids and rapid mesh-resolution changes, while adjustable pressure/velocity DOF ratios target spurious mode branches.The examples are the modified RT0-Q0 pair for the non-orthogonal cubed sphere and BDFM1-P1DG for triangles.Source: §5 · p. 21
  • The augmented pair's discontinuous Galerkin continuity equation is locally conservative, has minimal dispersion and diffusion errors, and can be made TVB with appropriate slope limiters.Conservative and consistent tracer transport must also use the pressure space, making tracer transport discontinuous Galerkin.Source: §5 · p. 21
  • Although BDFM(k)-PkDG spaces lack a 2:1 velocity-to-pressure DOF ratio, higher-order variants can achieve it by augmenting BDM(k) spaces with components that vanish on element boundaries.For k > 1, the pressure P(k) space must also be augmented to obtain stable pairs with twice as many velocity DOFs as pressure DOFs per triangle.Source: §5 · p. 21
Loading 1103.2440v3…