Semi-implicit integrations of perturbation equations for all-scale atmospheric dynamics

This paper presents a generalised perturbation form of the nonhydrostatic partial differen- tial equations (PDEs) that govern dynamics of all-scale global atmospheric ﬂows. There can be many alternative perturbation forms for any given system of the governing PDEs, depending on the assumed ambient state about which perturbations are taken and on subject preferences in the numerical model design. All such forms are mathematically equivalent, yet they have different implications for the design and accuracy of effective semi-implicit numerical integrators of the governing PDEs. Practical and relevant arguments are presented in favour of perturbation forms that maximise the degree of implicitness of the associated integrators. Two optional forms are implemented in the high-performance ﬁnite-volume module (IFS-FVM) for simulating global all-scale atmospheric ﬂows Smolarkiewicz et al. (2016) [32]. Their relative performance is veriﬁed with a class of ambient states of reduced complexity. A series of numerical simulations of the planetary baroclinic instability assumes geostrophically balanced zonally uniform ambient ﬂows with signiﬁcant meridional and vertical shear and exempliﬁes accuracy gains enabled by the generalised perturbation approach. The newly developed semi-implicit integrators show the potential for numerically accurate separation of a background state from ﬁnite-amplitude perturbations of the global atmosphere and provide a base for further development. ©


Introduction
An important characteristic of the atmospheric dynamics is that it constitutes a relatively small perturbation about dominant hydrostatic and geostrophic balances established in effect of the Earth gravity, rotation, stably-stratified thermal structure of its atmosphere and the energy provided by the incoming flux of solar radiation. For illustration, consider that the total conversion of available potential energy to kinetic energy in the global atmosphere is only 0.26% of the incoming solar radiation at the top of the atmosphere; cf. Fig. 1 in [46]. Preserving this fundamental equilibrium, while accurately resolving the perturbations about it, conditions the design of atmospheric models and subjects their numerical procedures to intricate stability and accuracy requirements. In particular, extended-range (11-45 days) and seasonal weather predictions (up to 1 year ahead) are subjected to dominant initial and model biases in the simulation of the atmosphere/ocean and suffer from low signal-to-noise ratios. In seasonal forecasting, this proofs a difficult problem for the detection of anomalies, especially when there is also a model drift relative to long-term biases in the simulated quantities [16]. The conceptual idea underlying the approach presented in this paper is a numerically accurate separation and dual evolution of an arbitrary (but advantageously selected) ambient state and finite-amplitude perturbations around this state, to improve conservation properties, mean state bias and drift, and the reliability of extended-range predictions.
Given the nature of atmospheric dynamics, it is compelling to formulate governing PDEs in terms of perturbation variables about an arbitrary state of the atmosphere that already satisfies some or all of the dominant balances. In fact, this is a standard approach in small-scale atmospheric modelling based on soundproof equations, such as the classical incompressible Boussinesq system or its anelastic or pseudo-incompressible generalisations. The key role of the reference state in these equation systems is to uncouple a dominating hydrostatic balance from buoyancy driven motions, and to justify linearisations required to achieve a physically-relevant complexity reduction. Otherwise, the reference state does not have to be unique nor represent the actual mean state of the system. The primary strength of established reference states is their generally recognized theoretical and practical validity. However, such reference states are overly limited, most often consisting only of a single stably or neutrally stratified vertical profile of thermodynamic variables. To alleviate this limitation, soundproof models often account for an additional, more general balanced state, hereafter referred to as an "ambient" state [23,24].
The notion of ambient states is distinct from that of the reference state. The utilization of an ambient state is justified by expediency and optional for any system of the governing equations. The role of ambient states is to enhance the efficacy of numerical simulation-e.g. by facilitating the design of the initial and boundary conditions and/or improving the conditioning of elliptic boundary value problems-without resorting to linearisation of the system. The underlying assumption is that the ambient state satisfies the governing equations of the full problem of interest, so that subtracting its own minimal set of PDEs from the governing equations can provide a useful perturbation form of the governing system. In general, ambient states are not limited to stable or neutral stratifications [2], can be spatially and temporally varying to represent, e.g. thermally balanced large-scale steady flows in atmospheric models [24,31] or prescribe oceanic tidal motions [38]. Here, we present a generalised formalism allowing for an arbitrary ambient state that satisfies the actual governing equations. However, the developed numerical apparatus is more general, as it offers technical advantages equally applicable to approximate ambient states, which can be particularly useful for weather and climate applications. We shall return to this point while concluding the paper. All proposed theoretical/numerical developments directly pertain to two advanced all-scale modelling systems, the Eulerian-Lagrangian (EULAG) research model for multiscale flows [20,30,31] and the Finite-Volume Module (IFS-FVM) [32,11,33] in the Integrated Forecasting System (IFS) of the European Centre for Medium-Range Weather Forecasts (ECMWF).
The paper is organised as follows: Section 2 derives the generalised perturbation form of the nonhydrostatic compressible equations of all-scale atmospheric dynamics. The related numerical procedures are described in Section 3. Section 4 substantiates technical developments of the preceding sections, with a series of extended-range simulations of the planetary baroclinic instability focused on the accuracy gains relevant to global weather prediction. Section 5 concludes the paper.

Generic form
The generalised PDEs of EULAG [31] and IFS-FVM [32] assume the compressible Euler/Navier-Stokes equations under gravity on a rotating sphere as default, but include reduced soundproof equations [15,5] as options. From the perspective of numerics, the design of the semi-implicit integrators in FVM follows the same path for the compressible and the soundproof system. Hence, we focus on the most general case of the compressible Euler equations. For simplicity, we start the presentation with the physically more intuitive advective form of the inviscid governing equations formulated on a rotating sphere, and introduce the corresponding conservation forms implemented in FVM afterwards in Section 3. With these caveats, the advective forms of equations for density ρ, potential temperature θ and the physical velocity vector u can be compactly written as Here, d/dt = ∂/∂t + v · ∇, with ∇ = (∂ x , ∂ y , ∂ z ) representing a vector of spatial partial derivatives with respect to the generalised curvilinear coordinates x = (x, y, z) [19,42,9] (assumed stationary for the purpose of this paper) and v = G T u, where G denotes a 3 × 3 matrix of known metric coefficients. Furthermore, G symbolises the Jacobian-i.e. the square root of the determinant of the metric tensor-whereas f and M(u) respectively mark the vectors of Coriolis parameter and metric forcings (i.e., Christoffel terms) shown for a plain geospherical framework in Appendix A. In (1b), H symbolises a heat source/sink, including diffusion. In (1c), g is the vector of gravitational acceleration, while D(u) symbolises momentum dissipation. The pressure variable φ = c p θ 0 π renormalises the Exner-function of pressure, π := (p/p 0 ) R/c p = T /θ , upon which it satisfies the gas law