We develop a discretisation of the semigeostrophic rotating shallow water equations, based upon their optimal transport formulation. This takes the form of a Moreau-Yoshida regularisation of the Wasserstein metric. Solutions of the optimal transport formulation provide the shallow water layer depth represented as a measure, which is itself the push forward of an evolving measure under the semigeostrophic coordinate transformation. First, we propose and study an entropy regularised version of the rotating shallow water equations. Second, we discretise the regularised problem by replacing both measures with weighted sums of Dirac measures, and approximate the (squared) L2 norm of the layer depth, which defines the potential energy. We propose an iterative method to solve the discrete optimisation problem relating the two measures, and analyse its convergence. The iterative method is demonstrated numerically and applied to the solution of the time-dependent shallow water problem in numerical examples.

Modern high-performance computers are massively parallel; for many partial differential equation applications spatial parallelism saturates long before the computer's capability is reached. Parallel-in-time methods enable further speedup beyond spatial saturation by solving multiple time steps simultaneously to expose additional parallelism. ParaDiag is a particular approach to parallel-in-time methods based on preconditioning the simultaneous time step system with a perturbation that allows block diagonalisation via a Fourier transform in time. In this article, we introduce asQ, a new library for implementing ParaDiag parallel-in-time methods, with a focus on applications in the geosciences, especially weather and climate. asQ is built on Firedrake, a library for the automated solution of finite element models, and the PETSc library of scalable linear and nonlinear solvers. This enables asQ to build ParaDiag solvers for general finite element models and provide a range of solution strategies, making testing a wide array of problems straightforward. We use a quasi-Newton formulation that encompasses a range of ParaDiag methods and expose building blocks for constructing more complex methods. The performance and flexibility of asQ is demonstrated on a hierarchy of linear and nonlinear atmospheric flow models. We show that ParaDiag can offer promising speedups and that asQ is a productive testbed for further developing these methods.

We present a particle filtering algorithm for stochastic models on infinite dimensional state space, making use of Girsanov perturbations to nudge the ensemble of particles into regions of higher likelihood. We argue that the optimal control problem needs to couple control variables for all of the particles to maintain an ensemble with good effective sample size (ESS). We provide an optimisation formulation that separates the problem into three stages, separating the nonlinearity in the ESS term in the functional with the nonlinearity due to the forward problem, and allowing independent parallel computation for each particle when calculations are performed over control variable space. The particle filter is applied to the stochastic Kuramoto-Sivashinsky equation, and compared with the temper-jitter particle filter approach. We observe that whilst the nudging filter is over spread compared to the temper-jitter filter, it responds to extreme events in the assimilated data more quickly and robustly.

A new test case is presented for evaluating the compressible dynamical cores of the atmospheric models. The test case is based on a compressible vertical slice model that can be obtained by simple modification of a standard three dimensional compressible dynamical core. On the one hand, an advantage of the test case is that is quasi-2D, so it can be run quickly on a standard workstation, enabling rapid experimentation with numerical schemes and discretisation choices. On the other hand, the test case exhibits frontogenesis, a challenging regime for numerical discretisations which usually only arises in 3D model configurations for the compressible case. Numerical results of the test case using an implicit time-stepping method with a compatible finite element discretisation are presented as a reference solution. An example comparison between advective and vector-invariant forms for the advective nonlinearity in the velocity equation demonstrates one possible use of the scheme. The comparison shows a Hollingsworth-like instability when the vector invariant form is used.

I will discuss the application of compatible finite element methods to large scale atmosphere and ocean simulation. Compatible finite element methods extend Arakawa's “C-grid” finite difference scheme to the finite element world. They are constructed from a discrete de Rham complex, which is a sequence of finite element spaces which are linked by the operators of differential calculus. The use of discrete de Rham complexes to solve partial differential equations is well established, but in this talk I focus on the specifics of dynamical cores for simulating weather, oceans and climate. The most important consequence of the discrete de Rham complex is the Hodge-Helmholtz decomposition, which has been used to exclude the possibility of several types of spurious oscillations from linear equations of geophysical flow. This means that compatible finite element spaces provide a useful framework for building dynamical cores. In this talk I will introduce the main concepts of compatible finite element spaces, and discuss their wave propagation properties. I will then cover a selection of the following topics (depending on recent advances, and interests of the audience): practical application to numerical weather prediction and ocean models, structure preserving methods, and scalable iterative solver techniques.

The semigeostrophic equations are a frontogenesis model in atmospheric science. Existence of solutions both from the theoretical and numerical point of view is given under a change of variable involving the interpretation of the pressure gradient as an optimal transport map between the density of the fluid and its push forward. Thanks to recent advances in numerical optimal transportation, the computation of large scale discrete approximations can be envisioned. We study here the use of entropic optimal transport and its Sinkhorn Algorithm companion.

Accurate transport algorithms are crucial for computational fluid dynamics and more accurate and efficient schemes are always in development. One dimensional limiting is a commonly employed technique used to suppress nonphysical oscillations. However, the application of such limiters can reduce accuracy. It is important to identify the weakest set of sufficient conditions required on the limiter as to allow the development of successful numerical algorithms.The main goal of this paper is to identify new less restrictive sufficient conditions for flux form in-compressible advection to remain monotonic. First, we identify conditions in which the Spekreijse limiter region can fail to be monotonic for incompressible flux form advection and demonstrate this numerically. Then a convex combination argument is used to derive new sufficient conditions that are less restrictive than the Sweby region for a discrete maximum principle. This allows the introduction of two new more general limiter regions suitable for flux form incompressible advection.

In this study, we explore data assimilation for the Stochastic Camassa-Holm equation through the application of the particle filtering framework. Specifically, our approach integrates adaptive tempering, jittering, and nudging techniques to construct an advanced particle filtering system. All filtering processes are executed utilizing ensemble parallelism. We conduct extensive numerical experiments across various scenarios of the Stochastic Camassa-Holm model with transport noise and viscosity to examine the impact of different filtering procedures on the performance of the data assimilation process. Our analysis focuses on how observational data and the data assimilation step influence the accuracy and uncertainty of the obtained results.

We study parameterisation-independent closed planar curve matching as a Bayesian inverse problem. The motion of the curve is modelled via a curve on the diffeomorphism group acting on the ambient space, leading to a large deformation diffeomorphic metric mapping (LDDMM) functional penalising the kinetic energy of the deformation. We solve Hamilton's equations for the curve matching problem using the Wu-Xu element (Wu and Xu (2019) [12]) which provides mesh-independent Lipschitz constants for the forward motion of the curve, and solve the inverse problem for the momentum using Bayesian inversion. Since this element is not affine-equivalent we provide a pullback theory which expedites the implementation and efficiency of the forward map. We adopt ensemble Kalman inversion (EKI) using a negative Sobolev norm mismatch penalty to measure the discrepancy between the target and the ensemble mean shape. We provide several numerical examples to validate the approach.

The reformulation of the Met Office's dynamical core for weather and climate prediction previously described by the authors is extended to spherical domains using a cubed-sphere mesh. This paper updates the semi-implicit mixed finite-element formulation to be suitable for spherical domains. In particular the finite-volume transport scheme is extended to take account of non-uniform, non-orthogonal meshes and uses an advective-then-flux formulation so that increment from the transport scheme is linear in the divergence. The resulting model is then applied to a standard set of dry dynamical core tests and compared to the existing semi-implicit semi-Lagrangian dynamical core currently used in the Met Office's operational model.

When estimating quantities and fields that are difficult to measure directly, such as the fluidity of ice, from point data sources, such as satellite altimetry, it is important to solve a numerical inverse problem that is formulated with Bayesian consistency. Otherwise, the resultant probability density function for the difficult to measure quantity or field will not be appropriately clustered around the truth. In particular, the inverse problem should be formulated by evaluating the numerical solution at the true point locations for direct comparison with the point data source. If the data are first fitted to a gridded or meshed field on the computational grid or mesh, and the inverse problem formulated by comparing the numerical solution to the fitted field, the benefits of additional point data values below the grid density will be lost. We demonstrate, with examples in the fields of groundwater hydrology and glaciology, that a consistent formulation can increase the accuracy of results and aid discourse between modellers and observationalists. To do this, we bring point data into the finite element method ecosystem as discontinuous fields on meshes of disconnected vertices. Point evaluation can then be formulated as a finite element interpolation operation (dual-evaluation). This new abstraction is well-suited to automation, including automatic differentiation. We demonstrate this through implementation in Firedrake, which generates highly optimised code for solving Partial Differential Equations (PDEs) with the finite element method. Our solution integrates with dolfin-adjoint/pyadjoint, allowing PDE-constrained optimisation problems, such as data assimilation, to be solved through forward and adjoint mode automatic differentiation.

We describe a proof‐of‐concept development and application of a phase‐averaging technique to the nonlinear rotating shallow‐water equations on the sphere, discretised using compatible finite‐element methods. Phase averaging consists of averaging the nonlinearity over phase shifts in the exponential of the linear wave operator. Phase averaging aims to capture the slow dynamics in a solution that is smoother in time (in transformed variables), so that larger timesteps may be taken. We overcome the two key technical challenges that stand in the way of studying the phase averaging and advancing its implementation: (1) we have developed a stable matrix exponential specific to finite elements and (2) we have developed a parallel finite averaging procedure. Following recent studies, we consider finite‐width phase‐averaging windows, since the equations have a finite timescale separation. In our numerical implementation, the averaging integral is replaced by a Riemann sum, where each term can be evaluated in parallel. This creates an opportunity for parallelism in the timestepping method, which we use here to compute our solutions. Here, we focus on the stability and accuracy of the numerical solution. We confirm that there is an optimal averaging window, in agreement with theory. Critically, we observe that the combined time discretisation and averaging error is much smaller than the time discretisation error in a semi‐implicit method applied to the same spatial discretisation. An evaluation of the parallel aspects will follow in later work.

We derive a linearized rotating shallow water system modeling tides, which can be discretized by mixed finite elements. Unlike previous models, this model allows for multiple layers stratified by density. Like the single-layer case~\cite{kirby2021preconditioning} a weighted-norm preconditioner gives a (nearly) parameter-robust method for solving the resulting linear system at each time step, but the all-to-all coupling between the layers in the model poses a significant challenge to efficiency. Neglecting the inter-layer coupling gives a preconditioner that degrades rapidly as the number of layers increases. By a careful analysis of the matrix that couples the layers, we derive a robust method that requires solving a reformulated system that only involves coupling between adjacent layers. Numerical results obtained using Firedrake confirm the theory.

We present a compatible finite element discretisation for the vertical slice compressible Euler equations, at next-to-lowest order (i.e., the pressure space is bilinear discontinuous functions). The equations are numerically integrated in time using a fully implicit timestepping scheme which is solved using monolithic GMRES preconditioned by a linesmoother. The linesmoother only involves local operations and is thus suitable for domain decomposition in parallel. It allows for arbitrarily large timesteps but with iteration counts scaling linearly with Courant number in the limit of large Courant number. This solver approach is implemented using Firedrake, and the additive Schwarz preconditioner framework of PETSc. We demonstrate the robustness of the scheme using a standard set of testcases that may be compared with other approaches.

Line Intensity Mapping (LIM) is a new observational technique that uses low-resolution observations of line emission to efficiently trace the large-scale structure of the Universe out to high redshift. Common mm/sub-mm emission lines are accessible from ground-based observatories, and the requirements on the detectors for LIM at mm-wavelengths are well matched to the capabilities of large-format arrays of superconducting sensors. We describe the development of an $ R$ = $\lambda / \Delta \lambda = 300$ on-chip superconducting filter-bank spectrometer covering the 120–180 GHz band for future mm-LIM experiments, focusing on SPT-SLIM, a pathfinder LIM instrument for the South Pole Telescope. Radiation is coupled from the telescope optical system to the spectrometer chip via an array of feedhorn-coupled orthomode transducers. Superconducting microstrip transmission lines then carry the signal to an array of channelizing half-wavelength resonators, and the output of each spectral channel is sensed by a lumped element kinetic inductance detector (leKID). Key areas of development include incorporating new low-loss dielectrics to improve both the achievable spectral resolution and optical efficiency and development of a robust fabrication process to create a galvanic connection between ultra-pure superconducting thin-films to realize multi-material (hybrid) leKIDs. We provide an overview of the spectrometer design, fabrication process, and prototype devices.

This paper introduces the r-Camassa–Holm (r-CH) equation, which describes a geodesic flow on the manifold of diffeomorphisms acting on the real line induced by the W1,r metric. The conserved energy for the problem is given by the full W1,r norm. For r = 2, we recover the Camassa–Holm equation. We compute the Lie symmetries for r-CH and study various symmetry reductions. We introduce singular weak solutions of the r-CH equation for r⩾2 and demonstrates their robustness in numerical simulations of their nonlinear interactions in both overtaking and head-on collisions. Several open questions are formulated about the unexplored properties of the r-CH weak singular solutions, including the question of whether they would emerge from smooth initial conditions.

<p>Parallel-in-time algorithms provide a route to increased parallelism for weather and climate models, addressing the issue of how to make efficient use of future supercomputers. In this talk I will present an overview of the approaches implemented in Gusto, the compatible finite element dynamical core toolkit build on top of the Firedrake finite element library. Compatible finite element methods are of interest for weather and climate modelling due to their conservation and wave propagation properties on non-orthogonal meshes such as the cubed-sphere. These non-orthogonal meshes allow for better scaling from spatial domain decomposition than meshes based on the latitude-longitude grid which have grid points clustered at the poles. However, the sequential nature of classical timestepping algorithms is a bottleneck to increased parallelisation. Numerical weather prediction is a challenging application for time-parallel schemes due to the hyperbolic nature of the partial differential equations that make up the dynamical core. Several different time-parallel schemes are under investigation in Gusto: parallel exponential integrators using a rational approximation (REXI); asymptotic parareal, which uses averaged equations to construct the coarse approximation; and schemes based on deferred correction. I will give an overview of these methods and present the latest results and challenges.</p>

Compatible finite-element discretisations for the atmospheric equations of motion have recently attracted considerable interest. Semi-implicit timestepping methods require the repeated solution of a large saddle-point system of linear equations. Preconditioning this system is challenging, since the velocity mass matrix is nondiagonal, leading to a dense Schur complement. Hybridisable discretisations overcome this issue: weakly enforcing continuity of the velocity field with Lagrange multipliers leads to a sparse system of equations, which has a similar structure to the pressure Schur complement in traditional approaches. We describe how the hybridised sparse system can be preconditioned with a non-nested two-level preconditioner. To solve the coarse system, we use the multigrid pressure solver that is employed in the approximate Schur complement method previously proposed by the some of the authors. Our approach significantly reduces the number of solver iterations. The method shows excellent performance and scales to large numbers of cores in the Met Office next-generation climate and weather prediction model LFRic.

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

This paper introduces the r-Camassa-Holm (r-CH) equation, which describes a geodesic flow on the manifold of diffeomorphisms acting on the real line induced by the W1,r metric. The conserved energy is for the problem is given by the full W1,r norm and the for r = 2, we recover the Camassa-Holm equation. We compute the Lie symmetries for r-CH and study various symmetry reductions. We introduce singular weak solutions of the r-CH equation for r >= 2 and demonstrates their robustness in numerical simulations of their nonlinear interactions in both overtaking and head-on collisions. Several open questions are formulated about the unexplored properties of the r-CH weak singular solutions, including the question of whether they would emerge from smooth initial conditions.