Numerical wave propagation for the triangular P1DGP1_{DG}-P2P2 finite element pair

Inertia-gravity mode and Rossby mode dispersion properties are examined for discretisations of the linearized rotating shallow-water equations using the $P1_{DG}$-$P2$ finite element pair on arbitrary triangulations in planar geometry. A discrete Helmholtz decomposition of the functions in the velocity space based on potentials taken from the pressure space is used to provide a complete description of the numerical wave propagation for the discretised equations. In the $f$-plane case, this decomposition is used to obtain decoupled equations for the geostrophic modes, the inertia-gravity modes, and the inertial oscillations. As has been noticed previously, the geostrophic modes are steady. The Helmholtz decomposition is used to show that the resulting inertia-gravity wave equation is third-order accurate in space. In general the \pdgp finite element pair is second-order accurate, so this leads to very accurate wave propagation. It is further shown that the only spurious modes supported by this discretisation are spurious inertial oscillations which have frequency $f$, and which do not propagate. The Helmholtz decomposition also allows a simple derivation of the quasi-geostrophic limit of the discretised $P1_{DG}$-$P2$ equations in the $\beta$-plane case, resulting in a Rossby wave equation which is also third-order accurate.Comment: Revised version prior to final journal submissio

Mixed finite elements for numerical weather prediction

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.Comment: Revision after referee comment

Two-dimensional evaluation of atham-fluidity, a nonhydrostatic atmospheric model using mixed continuous/discontinuous finite elements and anisotropic grid optimization

AbstractThis paper presents the first attempt to apply the compressible nonhydrostatic Active Tracer High-Resolution Atmospheric Model–Fluidity (ATHAM-Fluidity) solver to a series of idealized atmospheric test cases. ATHAM-Fluidity uses a hybrid finite-element discretization where pressure is solved on a continuous second-order grid while momentum and scalars are computed on a first-order discontinuous grid (also known as ). ATHAM-Fluidity operates on two- and three-dimensional unstructured meshes, using triangular or tetrahedral elements, respectively, with the possibility to employ an anisotropic mesh optimization algorithm for automatic grid refinement and coarsening during run time. The solver is evaluated using two-dimensional-only dry idealized test cases covering a wide range of atmospheric applications. The first three cases, representative of atmospheric convection, reveal the ability of ATHAM-Fluidity to accurately simulate the evolution of large-scale flow features in neutral atmospheres at rest. Grid convergence without adaptivity as well as the performances of the Hermite–Weighted Essentially Nonoscillatory (Hermite-WENO) slope limiter are discussed. These cases are also used to test the grid optimization algorithm implemented in ATHAM-Fluidity. Adaptivity can result in up to a sixfold decrease in computational time and a fivefold decrease in total element number for the same finest resolution. However, substantial discrepancies are found between the uniform and adapted grid results, thus suggesting the necessity to improve the reliability of the approach. In the last three cases, corresponding to atmospheric gravity waves with and without orography, the model ability to capture the amplitude and propagation of weak stationary waves is demonstrated. This work constitutes the first step toward the development of a new comprehensive limited area atmospheric model.This research has received funding from the European Union Seventh Framework Program (FP7/2007-2013) under Grant agreement 603663 for the research project PEARL (Preparing for Extreme And Rare events in coastaL regions). The EPSRC multiphase program grant MEMPHIS is also acknowledged.This is the author accepted manuscript. The final version is available from the American Meteorological Society via http://dx.doi.org/10.1175/MWR-D-15-0398.

A correction to the enhanced bottom drag parameterisation of tidal turbines

Hydrodynamic modelling is an important tool for the development of tidal stream energy projects. Many hydrodynamic models incorporate the effect of tidal turbines through an enhanced bottom drag. In this paper we show that although for coarse grid resolutions (kilometre scale) the resulting force exerted on the flow agrees well with the theoretical value, the force starts decreasing with decreasing grid sizes when these become smaller than the length scale of the wake recovery. This is because the assumption that the upstream velocity can be approximated by the local model velocity, is no longer valid. Using linear momentum actuator disc theory however, we derive a relationship between these two velocities and formulate a correction to the enhanced bottom drag formulation that consistently applies a force that remains closed to the theoretical value, for all grid sizes down to the turbine scale. In addition, a better understanding of the relation between the model, upstream, and actual turbine velocity, as predicted by actuator disc theory, leads to an improved estimate of the usefully extractable energy. We show how the corrections can be applied (demonstrated here for the models MIKE 21 and Fluidity) by a simple modification of the drag coefficient

Geostrophic balance preserving interpolation in mesh adaptive shallow-water ocean modelling

The accurate representation of geostrophic balance is an essential requirement for numerical modelling of geophysical flows. Significant effort is often put into the selection of accurate or optimal balance representation by the discretisation of the fundamental equations. The issue of accurate balance representation is particularly challenging when applying dynamic mesh adaptivity, where there is potential for additional imbalance injection when interpolating to new, optimised meshes. In the context of shallow-water modelling, we present a new method for preservation of geostrophic balance when applying dynamic mesh adaptivity. This approach is based upon interpolation of the Helmholtz decomposition of the Coriolis acceleration. We apply this in combination with a discretisation for which states in geostrophic balance are exactly steady solutions of the linearised equations on an f-plane; this method guarantees that a balanced and steady flow on a donor mesh remains balanced and steady after interpolation onto an arbitrary target mesh, to within machine precision. We further demonstrate the utility of this interpolant for states close to geostrophic balance, and show that it prevents pollution of the resulting solutions by imbalanced perturbations introduced by the interpolation

Compatible finite element methods for numerical weather prediction

This article takes the form of a tutorial on the use of a particular class of mixed finite element methods, which can be thought of as the finite element extension of the C-grid staggered finite difference method. The class is often referred to as compatible finite elements, mimetic finite elements, discrete differential forms or finite element exterior calculus. We provide an elementary introduction in the case of the one-dimensional wave equation, before summarising recent results in applications to the rotating shallow water equations on the sphere, before taking an outlook towards applications in three-dimensional compressible dynamical cores.Comment: To appear in ECMWF Seminar proceedings 201

Towards a fully unstructured ocean model for ice shelf cavity environments: Model development and verification using the Firedrake finite element framework

Numerical studies of ice flow have consistently identified the grounding zone of outlet glaciers and ice streams (the region where ice starts to float) as crucial for predicting the rate of grounded ice loss to the ocean. Owing to the extreme environments and difficulty of access to ocean cavities beneath ice shelves, field observations are rare. Estimates of melt rates derived from satellites are also difficult to make near grounding zones with confidence. Therefore, numerical ocean models are important tools to investigate these critical and remote regions. The relative inflexibility of structured grid models means, however, that they can struggle to resolve these processes in irregular cavity geometries near grounding zones. To help solve this issue, we present a new nonhydrostatic unstructured mesh model for flow under ice shelves built using the Firedrake finite element framework. We demonstrate our ability to simulate full ice shelf cavity domains using the community standard ISOMIP+ Ocean0 test case and compare our results against those obtained with the popular MITgcm model. Good agreement is found between the two models, despite their use of different discretisation schemes and the sensitivity of the melt rate parameterisation to grid resolution. Verification tests based on the Method of Manufactured Solutions (MMS) show that the new model discretisation is sound and second-order accurate. A main driver behind using Firedrake is the availability of an automatically generated adjoint model. Our first adjoint calculations, of sensitivities of melt rate with respect to different inputs in an idealised grounding zone domain, are promising and point to the ability to address a number of important questions on ocean influence on ice shelf vulnerability in the future