Analytical and adaptable initial conditions for dry and moist baroclinic waves in the global hydrostatic model OpenIFS (CY43R3)

This article presents a description of an analytical, stable, and flexible initial background state for both dry and moist baroclinic wave simulation on an aquaplanet in order to test the dynamical core of numerical weather prediction models and study the dynamics and evolution of extratropical cyclones. The initial background state is derived from an analytical zonal wind speed field, or jet structure, and the hydrostatic primitive equations for moist adiabatic and frictionless flow in spherical coordinates. A baroclinic wave can develop if a perturbation is added to the zonal wind speed field. This new baroclinic wave configuration has been implemented in the Open Integrated Forecasting System (OpenIFS) CY43R3, a global numerical weather prediction model developed by the European Centre for Medium-Range Weather Forecasts. In total, seven parameters can be used to control the generation of the initial background state and hence the development of the baroclinic waves in the OpenIFS configuration file: the jet's width, the jet's height, the maximum zonal mean wind speed of the jet, the horizontal mean of the surface virtual temperature, the surface relative humidity, the lapse rate, and the surface roughness. Nine dry and nine moist initial background states have been generated to test their stability without perturbations. The meteorological stability of the initial states is investigated by examining the spatial distributions of the equivalent potential temperature, the absolute vorticity, and the Brunt–Väisälä frequency. Moreover, the root mean square error (RMSE) of the zonal wind speed has been computed to assess their numerical stability. Finally, six dry and six moist initial background state have been used with an unbalanced perturbation to ensure that the baroclinic life cycles that develop are physically realistic. The resulting baroclinic wave is shown to be sensitive to the jet's width. This configuration for baroclinic wave simulations will be used to create a large ensemble of baroclinic life cycles to study how extratropical cyclones may evolve in the future.

1 Introduction

General circulation models (GCMs) are an important tool to predict the extent of global climate change as documented in the IPCC reports (IPCC, 2021). These GCMs provide numerical solutions to the governing equations of the atmosphere. They can take into consideration real-world data to predict the short-term evolution of the weather or they can be used to simulate idealised weather systems to study specific phenomena of our climate such as convection (Khairoutdinov et al., 2022) or baroclinic waves (Ullrich et al., 2015). Baroclinic waves are the synoptic-scale patterns of high- and low-pressure systems that develop in the mid-latitudes. These waves develop due to the release of baroclinic instability, and the resulting patterns are important parts of the Earth's global circulation as they transport energy polewards (Simmons and Hoskins, 1978; Thorncroft et al., 1993; Beare, 2007).

Two main reasons exist to perform baroclinic wave simulation (BWS). To further improve our weather prediction and climate models, the dynamical cores have to be tested. Most BWS experiments specify a zonally uniform solution to the hydrostatic primitive equations that is statically stable and stable to inertial and symmetric instabilities and then run a numerical model with this specified as the initial state. This type of simulation tests the ability of the numerical model to retain this exact solution in the presence of numerical errors. These initial states are baroclinically and barotropically unstable, and therefore adding a perturbation triggers the development of a baroclinic wave – to which there is no exact solution (Hoskins and Simmons, 1975; Simmons and Hoskins, 1975; Jablonowski and Williamson, 2006). The baroclinic wave development can be simulated on an f-plane or, its extension, a β-plane in both Cartesian and spherical geometries (Feldstein and Held, 1989; Staniforth and White, 2007; Ullrich et al., 2015). These models are often less expensive to run from a computational point of view by simplifying the Coriolis forces. Another approximation that is used alongside f-plane and β plane is the restriction of the size of the model domain where the zonal extent of the domain is set roughly equal to the most unstable wavelength (∼ 4000 km) (Hoskins et al., 1977; Ullrich et al., 2015). With this limitation, any upstream or downstream development is forced to occur on top of the main perturbation. This representation can efficiently display the energy propagation of the baroclinic wave (Hoskins et al., 1977). However, it becomes increasingly difficult to study dynamical and synoptic properties of the cyclones without the realistic simulation of upstream and downstream developments. Moreover, numerous descriptions and specifications of the initial states are available for Cartesian geometry and for channel models (Hoskins et al., 1977; Feldstein and Held, 1989; Wang and Polvani, 2011; Ullrich et al., 2015; Terpstra and Spengler, 2015) as well as for spherical geometry and fully global models (Polvani et al., 2004; Jablonowski and Williamson, 2006; Staniforth and White, 2011; Ullrich et al., 2014; Hughes and Jablonowski, 2023), but none proposes an initial state with tunable parameters (e.g. defining the width of the jet), which is one of the main motivations for this work.

Baroclinic wave simulations are of interest to study extratropical cyclones, extreme cases which will likely become more frequent in the future (IPCC, 2021). Extreme cyclones are characterised by strong winds, heavy precipitation, and powerful ocean waves. Consequently, these extreme events can damage infrastructure, forests, and homes, cause flooding, and result in injuries and even death. Depending on the location and the state of the large-scale background environment that the cyclone develops in, the structure and intensity of a given cyclone can vary considerably (Tang et al., 2020). Traditionally, BWSs were adiabatic and the simulations were run without physics parameterisation schemes (Simmons and Hoskins, 1978; Thorncroft et al., 1993), the main reason being that the synoptic-scale dynamics of baroclinic waves can be largely explained by the classic quasigeostrophic theories of dry baroclinic instability (Charney, 1947; Eady, 1949). Furthermore, the impact of latent heat on extratropical cyclone intensity has been heavily investigated in the last 3 decades (e.g.  Kuo et al., 1991; Stoelinga, 1996; Willison et al., 2013; Park et al., 2021). For example, diabatic processes are important to the evolution of precipitation (Kuo et al., 1991; Park et al., 2021) and smaller-scale systems (Stoelinga, 1996). Moreover, to predict how cyclones may change in the future, it is necessary to determine how sensitive baroclinic waves are to changes in temperature and moisture. This motivated the development of new BWSs which were designed to be run with physics parameterisation schemes active and moisture present (Beare, 2007; Kirshbaum et al., 2018; Tierney et al., 2018; Rantanen et al., 2019).

The sensitivity of the resultant extratropical cyclones to the jet structure has been studied (Thorncroft et al., 1993; Shapiro et al., 1999; Rupp and Birner, 2021). Popular zonal jet structures are zonal jet 1, zonal jet 2, and zonal jet 3 (denoted Z1, Z2, and Z3) resulting in baroclinic life cycles 1, 2, and 3, respectively (denoted LC1, LC2, and LC3) (Thorncroft et al., 1993; Agustí-Panareda et al., 2005). Z2 and Z3 differ from the zonally quasi-symmetric jet of Z1 by including a cyclonic (Z2) or anticyclonic (Z3) barotropic shear (Thorncroft and Hoskins, 1990; Thorncroft et al., 1993; Shapiro et al., 1999; Polvani and Esler, 2007). Depending on the barotropic shear, baroclinic life cycles have different structures and intensities (Agustí-Panareda et al., 2005). However, it is difficult to control the height and the width of the jet and test the sensitivity of the BWSs to these parameters due to the finite number of jet structures tested. Here, an initial state is developed, in which the jet's width and vertical structure can be varied in addition to the jet strength. However, setting up a balanced and flexible background state with several jet structures is difficult. The challenges lie in balancing the initial conditions at high resolutions in state-of-the-art models. Few cases are fully documented and can be difficult to reproduce. Many are based on the model developed by Hoskins and Simmons (1975) and rely on numerical integration, which is prone to truncation errors. Some are based on Cartesian geometry as presented in Kirshbaum et al. (2018) and their jet structures are obtained from the potential vorticity inversion method (Heckley and Hoskins, 1982; Olson and Colle, 2007). Having an analytical structure of the jet may allow more control on its structure and strength.

The aim of this study is to describe a balanced and flexible initial background state for a baroclinic life cycle experiment that can be entirely expressed analytically and that produces relatively realistic weather systems. The analytical solution is derived from a steady-state momentum equation for the meridional wind speed. In other words, the meridional wind speed is set to 0.0 m s⁻¹, which leads to a gradient–wind balance. Moreover, the proposed background state is also in hydrostatic balance; i.e. the hydrostatic equation was used to derive the virtual temperature anomaly from the geopotential field. The initial background state is also based on a flexibly defined jet structure and can furthermore be initialised with moisture by changing the relative humidity profile. The theoretical description and derivation of the initial state with mathematical formulae together with the method for including moisture are presented in Sect. 2. The technical implementation into the global state-of-the-art numerical weather prediction model, the Open Integrated Forecasting System (OpenIFS), is described in Sect. 3. The different experiments are described in Sect. 4 and their associated results in Sect. 5. First, the new initial states (both dry with no physics parameterisation schemes and moist with almost all physics parameterisation schemes) are run for 15 d with no perturbation to confirm that the initial states are indeed stable both numerically and from a meteorological perspective. Second, six dry and six moist initial background states have been used with an unbalanced perturbation to generate baroclinic life cycles. The evolution of the resulting baroclinic waves that develop from our default dry and moist initial state is shown. The results from the moist default simulation also show how the precipitation patterns form as the cyclones develop.

2 Initial condition for the baroclinic wave

This section presents a balanced, steady-state initial condition for a 3D hydrostatic atmospheric model in spherical coordinates with a flexible jet structure. The moisture field is defined to be consistent with the virtual temperature field. This section describes the theoretical background for the initial conditions and is divided into four parts: (1) analytical derivation of the geopotential and virtual temperature fields, (2) initialisation of moisture, (3) initialisation of the surface (both sea surface temperature – SST – and roughness), and (4) description of the unbalanced perturbation, which when added triggers the development of the baroclinic wave.

2.1 Analytical geopotential and virtual temperature fields

The derivation of the analytical initial conditions for geopotential and virtual temperature fields starts from the primitive equations for moist adiabatic and frictionless flow in spherical coordinates and normalised pressure levels for a planet with no topography (i.e. surface geopotential is zero). The geopotential and virtual temperature anomaly fields are derived from hydrostatic equations, and the derived initial states apply to both hydrostatic and nonhydrostatic models. The geopotential and virtual temperature fields are described as the horizontal mean field as a function of vertical levels plus an anomaly field which is a function of longitude, latitude, and vertical levels ($\mathit{\lambda },\mathit{\varphi },\mathit{\eta }$, respectively). In OpenIFS, the vertical η levels are defined as $\mathit{\eta }=p/p_{\mathrm{s}}=a/p_{\mathrm{s}}+b$, where p_s = 1013.25 hPa is the pressure at the surface pressure, and a and b are hybrid coefficients defined for each vertical resolution. The horizontal means proposed by Ullrich et al. (2015) are used in this background state. The analytical formula of the horizontal mean geopotential is described as

$$\begin{array}{}\text{(1)}& \langle \mathrm{\Phi }\left(\mathit{\eta }\right)\rangle ={\displaystyle \frac{T_{v,\mathrm{0}}g}{\mathit{\gamma }}}(\mathrm{1}-{\mathit{\eta }}^{\frac{R_{\mathrm{d}}\mathit{\gamma }}{g}})\end{array}$$

and the horizontal mean virtual temperature field as

$$\begin{array}{}\text{(2)}& \langle T_{\mathrm{v}}\left(\mathit{\eta }\right)\rangle =T_{v,\mathrm{0}}{\mathit{\eta }}^{\frac{R_{\mathrm{d}}\mathit{\gamma }}{g}},\end{array}$$

where γ is the specified lapse rate, T_v,0 the reference virtual temperature, R_d the gas constant for dry air, and g the gravity constant. The derivation for the horizontal mean geopotential is available in the Appendix.

To be able to derive the anomaly fields of geopotential and virtual temperature, a jet structure has been defined similar to the one proposed by Ullrich and Jablonowski (2012) and Ullrich et al. (2015) with the only difference being the power of the sine function. The power of the sine is described as 2n, allowing for a narrower jet when n increases. The chosen formula can be expressed as

$$\begin{array}{}\text{(3)}& u(\mathit{\lambda },\mathit{\varphi },\mathit{\eta })=-u_{\mathrm{0}}\ln \left(\mathit{\eta }\right)\exp \left[-{\left({\displaystyle \frac{\ln \mathit{\eta }}{b}}\right)}^{\mathrm{2}}\right]{\sin }^{\mathrm{2}n}\left(\mathrm{2}\mathit{\varphi }\right),\end{array}$$

where n is a positive integer defining the width of the jet, b is a nondimensional parameter representing the depth of the jet, and u₀ is the reference zonal wind speed and defines the zonal mean speed of the jet in the troposphere. As expressed by Eq. (3), the jet width decreases with an increase in n, and the height of the centre of the jet and the vertical width of the jet increase with increasing values of b. The jet reaches its maximum wind speed at ϕ = 45° and $\mathit{\eta }=\exp (-b/\sqrt{\mathrm{2}})$. Furthermore, the value for b needs to be positive and smaller than $-\sqrt{\mathrm{2}}\ln \left({\mathit{\eta }}_{\text{top}}\right)$, where η_top is the ratio between the top pressure level and the surface pressure level. For example, with a top pressure level of 0.01 hPa and surface pressure level of 1000 hPa, the upper limit for b is about 16. If the value of b exceeds the upper limit, the centre of the jet is located outside the model domain. The analytical geopotential and virtual temperature fields are solved for any jet structure defined by Eq. (3), i.e. for arbitrary values of n and b.

The derivation of the geopotential and virtual temperature anomaly fields starts from the primitive equations for moist adiabatic and frictionless flow (Holton and Hakim, 2013). Following the instructions given in Appendix A of Jablonowski and Williamson (2006), the geopotential anomaly field has been derived from the steady-state momentum equation for the meridional flow ($\partial v/\partial t=\mathrm{0}$) by inserting our choice of jet structure and solving for ${\mathrm{\Phi }}^{\prime }(\mathit{\lambda },\mathit{\varphi },\mathit{\eta })$:

$$\begin{array}{}\text{(4)}& {\displaystyle \frac{\mathrm{1}}{a}}{\displaystyle \frac{\partial {\mathrm{\Phi }}^{\prime }}{\partial \mathit{\varphi }}}=-u\left(\mathrm{2}\mathrm{\Omega }\sin \mathit{\varphi }+{\displaystyle \frac{u}{a}}\tan \mathit{\varphi }\right),\end{array}$$

where Ω is the angular velocity of the Earth, u the jet structure given by Eq. (3), and a the radius of the Earth. A steady-state solution leads to a gradient–wind balance, where the centrifugal, Coriolis, and pressure gradient forces are in balance. Integrating Eq. (4) analytically over ϕ results in

$$\begin{array}{}\text{(5)}& \begin{array}{rl}{\mathrm{\Phi }}^{\prime }(\mathit{\lambda },\mathit{\varphi },\mathit{\eta })=& -u_{\mathit{\eta }}\mathrm{2}a\mathrm{\Omega }{\mathrm{4}}^n\sum _{k=\mathrm{0}}^n\left(\genfrac{}{}{0ex}{}{n}{k}\right)(-\mathrm{1}{)}^k\\ & \times {\displaystyle \frac{\mathrm{1}}{\mathrm{2}(k+n)+\mathrm{1}}}{\cos }^{\mathrm{2}(k+n)+\mathrm{1}}\mathit{\varphi }\\ & -u_{\mathit{\eta }}^{\mathrm{2}}{\mathrm{16}}^n\sum _{k=\mathrm{0}}^{\mathrm{2}n-\mathrm{1}}\left(\genfrac{}{}{0ex}{}{\mathrm{2}n-\mathrm{1}}{k}\right)(-\mathrm{1}{)}^k\\ & \times {\displaystyle \frac{\mathrm{1}}{\mathrm{2}(k+\mathrm{2}n+\mathrm{1})}}{\sin }^{\mathrm{2}(k+\mathrm{2}n+\mathrm{1})}\mathrm{2}\mathit{\varphi }+{\mathrm{\Phi }}_{\mathrm{0}}\left(\mathit{\eta }\right).\end{array}\end{array}$$

Since the deviations of Φ^′ vanish when averaging horizontally, Φ₀ is solved by inserting ${\mathrm{\Phi }}^{\prime }(\mathit{\lambda },\mathit{\varphi },\mathit{\eta })$ in the horizontal mean equation

$$\begin{array}{}\text{(6)}& {\displaystyle \frac{\mathrm{1}}{\mathrm{4}\mathit{\pi }}}\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\underset{-\mathit{\pi }/\mathrm{2}}{\overset{\mathit{\pi }/\mathrm{2}}{\int }}{\mathrm{\Phi }}^{\prime }(\mathit{\lambda },\mathit{\varphi },\mathit{\eta })\cos \mathit{\varphi }\mathrm{d}\mathit{\varphi }\mathrm{d}\mathit{\lambda }=\mathrm{0},\end{array}$$

which then gives the analytical geopotential anomaly field Φ^′ as

$$\begin{array}{}\text{(7)}& \begin{array}{rl}{\mathrm{\Phi }}^{\prime }(\mathit{\lambda },\mathit{\varphi },\mathit{\eta })=& u_{\mathit{\eta }}a\mathrm{\Omega }{\mathrm{4}}^n\left(F_{\mathrm{3}}-\mathrm{2}{F}_{\mathrm{1}}\right)\\ & +u_{\mathit{\eta }}^{\mathrm{2}}{\mathrm{16}}^n\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{F}_{\mathrm{4}}-F_{\mathrm{2}}\right),\end{array}\end{array}$$

where

$$\begin{array}{}\text{(8a)}& {\displaystyle }{F}_{\mathrm{1}}=\sum _{k=\mathrm{0}}^n\left({\displaystyle \genfrac{}{}{0ex}{}{n}{k}}\right)(-\mathrm{1}{)}^k{\displaystyle \frac{\mathrm{1}}{\mathrm{2}(k+n)+\mathrm{1}}}{\cos }^{\mathrm{2}(k+n)+\mathrm{1}}\mathit{\varphi },\text{(8b)}& {\displaystyle }\begin{array}{rl}{F}_{\mathrm{2}}=& \sum _{k=\mathrm{0}}^{\mathrm{2}n-\mathrm{1}}\left({\displaystyle \genfrac{}{}{0ex}{}{\mathrm{2}n-\mathrm{1}}{k}}\right)(-\mathrm{1}{)}^k\\ & \times {\displaystyle \frac{\mathrm{1}}{\mathrm{2}(k+\mathrm{2}n+\mathrm{1})}}{\sin }^{\mathrm{2}(k+\mathrm{2}n+\mathrm{1})}\mathrm{2}\mathit{\varphi },\end{array}\text{(8c)}& {\displaystyle }{F}_{\mathrm{3}}=\sum _{k=\mathrm{0}}^n\left({\displaystyle \genfrac{}{}{0ex}{}{n}{k}}\right)(-\mathrm{1}{)}^k{\displaystyle \frac{\mathrm{1}}{\mathrm{2}(k+n)+\mathrm{1}}}\sqrt{\mathit{\pi }}{\displaystyle \frac{\mathrm{\Gamma }(k+n+\mathrm{3}/\mathrm{2})}{\mathrm{\Gamma }(k+n+\mathrm{2})}},\text{(8d)}& {\displaystyle }\begin{array}{rl}{F}_{\mathrm{4}}=& \sum _{k=\mathrm{0}}^{\mathrm{2}n-\mathrm{1}}\left({\displaystyle \genfrac{}{}{0ex}{}{\mathrm{2}n-\mathrm{1}}{k}}\right)(-\mathrm{1}{)}^k{\displaystyle \frac{\mathrm{1}}{\mathrm{2}(\mathrm{2}n+k+\mathrm{1})}}\\ & \times {\displaystyle \frac{\mathrm{2}}{\mathrm{2}(\mathrm{2}n+k+\mathrm{1})+\mathrm{1}}},\end{array}\text{(8e)}& {\displaystyle }{u}_{\mathit{\eta }}=u_{\mathrm{0}}\ln \mathit{\eta }\exp (-[\ln \mathit{\eta }/b{]}^{\mathrm{2}}).\end{array}$$

Note that $\left(\genfrac{}{}{0ex}{}{n}{k}\right)$ is the binomial coefficient representing the k unordered outcomes from n possibilities, Γ(x) is the gamma function for positive half-integer $x=z+\mathrm{1}/\mathrm{2}$ with z a positive integer, and 2n is the power of the sine in the jet structure.

The total geopotential field is described as the sum of the mean horizontal geopotential field and the anomaly geopotential field as

$$\begin{array}{}\text{(9)}& \begin{array}{rl}\mathrm{\Phi }(\mathit{\lambda },\mathit{\varphi },\mathit{\eta })& ={\displaystyle \frac{T_{v,\mathrm{0}}g}{\mathit{\gamma }}}\left(\mathrm{1}-{\mathit{\eta }}^{\frac{R_{\mathrm{d}}\mathit{\gamma }}{g}}\right)+u_{\mathit{\eta }}a\mathrm{\Omega }{\mathrm{4}}^n\left(F_{\mathrm{3}}-\mathrm{2}{F}_{\mathrm{1}}\right)\\ & +u_{\mathit{\eta }}^{\mathrm{2}}{\mathrm{16}}^n\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{F}_{\mathrm{4}}-F_{\mathrm{2}}\right).\end{array}\end{array}$$

The virtual temperature anomaly field is then derived by inserting ${\mathrm{\Phi }}^{\prime }(\mathit{\lambda },\mathit{\varphi },\mathit{\eta })$ into the hydrostatic equation and taking the derivative of Φ^′ with respect to η,

$$\begin{array}{}\text{(10)}& T_{\mathrm{v}}^{\prime }(\mathit{\lambda },\mathit{\varphi },\mathit{\eta })=-{\displaystyle \frac{\mathit{\eta }}{R_{\mathrm{d}}}}{\displaystyle \frac{\partial {\mathrm{\Phi }}^{\prime }(\mathit{\lambda },\mathit{\varphi },\mathit{\eta })}{\partial \mathit{\eta }}},\end{array}$$

which gives the virtual temperature field:

$$\begin{array}{}\text{(11)}& \begin{array}{rl}{T}_{\mathrm{v}}(\mathit{\lambda },\mathit{\varphi },\mathit{\eta })=& \langle T_{\mathrm{v}}\left(\mathit{\eta }\right)\rangle +T_{\mathrm{v}}^{\prime }(\mathit{\lambda },\mathit{\varphi },\mathit{\eta })\\ =& T_{v,\mathrm{0}}{\mathit{\eta }}^{\frac{R_{\mathrm{d}}\mathit{\gamma }}{g}}+{\displaystyle \frac{u_{\mathrm{0}}}{R_{\mathrm{d}}}}\exp [-(\ln \mathit{\eta }/b{)}^{\mathrm{2}}]\\ & \times \left[{\displaystyle \frac{\mathrm{2}(\ln \mathit{\eta }{)}^{\mathrm{2}}}{b^{\mathrm{2}}}}-\mathrm{1}\right]\left[a\mathrm{\Omega }{\mathrm{4}}^n\right(F_{\mathrm{3}}-\mathrm{2}{F}_{\mathrm{1}})\\ & +{\mathrm{16}}^n{u}_{\mathit{\eta }}(F_{\mathrm{4}}-\mathrm{2}{F}_{\mathrm{2}})],\end{array}\end{array}$$

where $F_{\mathrm{1}},F_{\mathrm{2}},F_{\mathrm{3}},F_{\mathrm{4}}$, and u_η are as defined in Eq. (8). A detailed step-by-step derivation of the analytical geopotential and temperature anomaly fields is available in Appendix A.

Figure 1 Relative humidity profile for RH₀ = 80 % as a function of height. The left-hand y axis shows η levels and the right-hand y axis the corresponding pressure levels (hPa).

2.2 Moisture initialisation

As stated in the Introduction, the proposed background state can be used in dry and moist case studies. In order to set the latter, a relative humidity profile with respect to water RH(η), depending on the model level η and the surface relative humidity RH₀, has been defined. It is inspired from the ERA-Interim (Romps, 2014) and ERA5 (Gamage et al., 2020) average relative humidity profile. The profile, as shown in Fig. 1, has a maximum value of the RH₀ at the surface, above which it decreases to 70 % of RH₀ at η=0.8. Between 0.8 and 0.3, RH is constant and for η<0.3 it again decreases linearly to 0 at η=0.1. Above 0.1, RH is set to 0 %. The profile is given as

$$\begin{array}{}\text{(12)}& \mathrm{RH}\left(\mathit{\eta }\right)=\left\{\begin{array}{ll}\mathrm{0.0}\mathit{\%}& \text{between }\mathit{\eta }=\mathrm{0}\text{ and }\mathrm{0.1}\\ (\mathrm{3.5}\mathit{\eta }-\mathrm{0.35}){\mathrm{RH}}_{\mathrm{0}}& \text{between }\mathit{\eta }=\mathrm{0.1}\text{ and }\mathrm{0.3}\\ \mathrm{0.7}{\mathrm{RH}}_{\mathrm{0}}& \text{between }\mathit{\eta }=\mathrm{0.3}\text{ and }\mathrm{0.8}\\ (\mathrm{1.5}\mathit{\eta }-\mathrm{0.5}){\mathrm{RH}}_{\mathrm{0}}& \text{above }\mathit{\eta }>\mathrm{0.8}.\end{array}\right.\end{array}$$

The specific humidity field $q(\mathit{\lambda },\mathit{\varphi },\mathit{\eta })$ is then computed to ensure concordance with the proposed virtual temperature and jet structure by assuming T=T_v. The specific humidity field is derived from the relative humidity (RH(η)), the saturation vapour pressure (e_s), and the saturation mixing ratio (w_s) using the Bolton approximation for the saturation vapour pressure (Bolton, 1980; Yau and Rogers, 1996) as presented in the following equations:

$$\begin{array}{}\text{(13a)}& {\displaystyle }{e}_{\mathrm{s}}(\mathit{\lambda },\mathit{\varphi },\mathit{\eta })=\mathrm{611.21}\exp {\displaystyle \frac{\mathrm{17.67}\left(T\right(\mathit{\lambda },\mathit{\varphi },\mathit{\eta })-\mathrm{273.15})}{T(\mathit{\lambda },\mathit{\varphi },\mathit{\eta })-\mathrm{29.65}}},\text{(13b)}& {\displaystyle }{w}_{\mathrm{s}}(\mathit{\lambda },\mathit{\varphi },\mathit{\eta })=\mathrm{0.622}{\displaystyle \frac{e_{\mathrm{s}}(\mathit{\lambda },\mathit{\varphi },\mathit{\eta })}{p(\mathit{\lambda },\mathit{\varphi },\mathit{\eta })-e_{\mathrm{s}}(\mathit{\lambda },\mathit{\varphi },\mathit{\eta })}},\end{array}$$

where $T(\mathit{\lambda },\mathit{\varphi },\mathit{\eta })$ is the temperature field (K), $e_{\mathrm{s}}(\mathit{\lambda },\mathit{\varphi },\mathit{\eta })$ is the saturation vapour pressure (Pa), and p is the pressure (Pa). Finally, $w_{\mathrm{s}}(\mathit{\lambda },\mathit{\varphi },\mathit{\eta })$ and RH(η) are used to infer $q(\mathit{\lambda },\mathit{\varphi },\mathit{\eta })$ as

$$\begin{array}{}\text{(14)}& q(\mathit{\lambda },\mathit{\varphi },\mathit{\eta })={\displaystyle \frac{w_{\mathrm{s}}(\mathit{\lambda },\mathit{\varphi },\mathit{\eta })\mathrm{RH}\left(\mathit{\eta }\right)}{w_{\mathrm{s}}(\mathit{\lambda },\mathit{\varphi },\mathit{\eta })\mathrm{RH}\left(\mathit{\eta }\right)+\mathrm{100}\mathit{\%}}}.\end{array}$$

The formulations of the virtual temperature and the specific humidity lead to the following expression for the temperature field:

$$\begin{array}{}\text{(15)}& T(\mathit{\lambda },\mathit{\varphi },\mathit{\eta })={\displaystyle \frac{T_{\mathrm{v}}(\mathit{\lambda },\mathit{\varphi },\mathit{\eta })}{\mathrm{1}+\mathrm{0.608}q(\mathit{\lambda },\mathit{\varphi },\mathit{\eta })}}.\end{array}$$

The process for updating the temperature and the specific humidity needs to be computed iteratively. The iteration starts by first setting T=T_v after which e_s is computed using the Bolton equation (Eq. 13a) and w_s and q are computed using Eqs. (13b) and (14), respectively. When q is updated, then $T=T(T_{\mathrm{v}},q)$ is updated. In the following iteration, the estimated T is again used to estimate a new T. Tests show that this iterative process converges quickly, and after 10 cycles, the algorithm has reached an optimum. The final result is T and q. Figure 2 shows how this iterative process works for computing T and q.

Figure 2 Flowchart showing how to compute T and q. During each iteration e_s, w_s, and q are computed from T after which the temperature T is updated using the new q and the virtual temperature T_v. The iterative process ends when the number of iterations is 10 and the final T and q are returned.

2.3 Surface initialisation: temperature and roughness

A uniform sea surface with no land has been chosen for the described background state. Thus, the experiments presented here were conducted using an aquaplanet setting, which is the traditional configuration of baroclinic wave simulation (Jablonowski and Williamson, 2006; Ullrich et al., 2015). The sea surface temperature SST is zonally uniform and is specified to equal the temperature field at η=1 (see Eqs. 11 and 15), which means negative temperatures are allowed and the zonal wind is equal to 0.0 m s⁻¹. If the SST differed from the near-surface atmospheric temperature, then in the moist cases there would be nonzero surface sensible heat fluxes, which could either heat or cool the boundary layer. Such fluxes could trigger convection, destabilising the proposed background state. The proposed SST is stated as

$$\begin{array}{}\text{(16)}& T_{\text{SST}}(\mathit{\lambda },\mathit{\varphi })={\displaystyle \frac{T_{v,\mathrm{0}}-\frac{u_{\mathrm{0}}a\mathrm{\Omega }}{R_{\mathrm{d}}}{\mathrm{4}}^n(F_{\mathrm{3}}-\mathrm{2}{F}_{\mathrm{1}})}{\mathrm{1}+\mathrm{0.608}q(\mathit{\lambda },\mathit{\varphi },\mathit{\eta }=\mathrm{1})}},\end{array}$$

where F₁ and F₃ are described in Eq. (8).

To complete the surface initialisation, the Charnock parameter is specified to control the surface roughness (Charnock, 1955). The surface roughness lengths for momentum (M), heat (H), and total water (Q) air–surface transfers are defined in OpenIFS (Eq. 3.26, ECMWF, 2017b) (and Eq. 25, Beljaars, 1995) as

$$\begin{array}{}\text{(17a)}& {\displaystyle }{z}_{\mathrm{0}M}& {\displaystyle }={\mathit{\alpha }}_M{\displaystyle \frac{\mathit{\nu }}{u_{\star }}}+{\mathit{\alpha }}_{\text{Ch}}{\displaystyle \frac{u_{\ast }^{\mathrm{2}}}{g}},\text{(17b)}& {\displaystyle }{z}_{\mathrm{0}H}& {\displaystyle }={\mathit{\alpha }}_H{\displaystyle \frac{\mathit{\nu }}{u_{\star }}},\text{(17c)}& {\displaystyle }{z}_{\mathrm{0}Q}& {\displaystyle }={\mathit{\alpha }}_Q{\displaystyle \frac{\mathit{\nu }}{u_{\star }}},\end{array}$$

where α_Ch is the Charnock parameter, u_⋆ the friction velocity, ν kinematic viscosity, and α_M, α_H, and α_Q constants set as 0.11, 0.40, and 0.62, respectively. Being able to tune the Charnock parameter allows the modification of the surface friction which previous studies have shown to influence the intensity of extratropical cyclones (Adamson et al., 2006; Sinclair et al., 2010) and the structure of warm and cold fronts (Hines and Mechoso, 1993; Sinclair and Keyser, 2015).

2.4 Initial perturbation

The baroclinic wave can be triggered by adding a localised unbalanced wind perturbation to a baroclinically unstable background state such as the one described in Sect. 2.1. A Gaussian perturbation was chosen and it was centred at $({\mathit{\lambda }}_c,{\mathit{\varphi }}_c)=(\frac{\mathit{\pi }}{\mathrm{9}},\frac{\mathrm{2}\mathit{\pi }}{\mathrm{9}})$, which corresponds to 40° N, 20° E (Jablonowski and Williamson, 2006; Ullrich et al., 2015). The equation of the perturbation is given by

$$\begin{array}{}\text{(18)}& u_{\mathit{ϵ}}(\mathit{\lambda },\mathit{\varphi },\mathit{\eta })=u_{\mathrm{p}}\exp \left[-{\left({\displaystyle \frac{r}{R}}\right)}^{\mathrm{2}}\right],\end{array}$$

where $R=\frac{a}{\mathrm{10}}$, u_p = 1.0 m s⁻¹, and r is the great circle distance given by

$$\begin{array}{}\text{(19)}& r=a\arccos (\sin {\mathit{\varphi }}_c\sin \mathit{\varphi }+\cos {\mathit{\varphi }}_c\cos \mathit{\varphi }\cos (\mathit{\lambda }-{\mathit{\lambda }}_c\left)\right).\end{array}$$

The final zonal wind field is obtained by adding u_ϵ to u at each grid point at all model levels:

$$\begin{array}{}\text{(20)}& u_{\text{total}}(\mathit{\lambda },\mathit{\varphi },\mathit{\eta })=u(\mathit{\lambda },\mathit{\varphi },\mathit{\eta })+u_{\mathit{ϵ}}(\mathit{\lambda },\mathit{\varphi },\mathit{\eta }).\end{array}$$

3 Implementation in OpenIFS

3.1 OpenIFS

The proposed background state has been implemented in the Open Integrated Forecasting System (OpenIFS) cycle 43R3v2 (CY43R3), which is based on the Integrated Forecasting System of the European Centre for Medium-Range Weather Forecasts (ECMWF) cycle 43R3 that was operational from July 2017 to June 2018 (ECMWF, 2017a). OpenIFS is a version of the Integrated Forecasting System model but does not include data assimilation capacities. Despite its name, OpenIFS is not open-source but available to universities and research institutions under license. The model is hydrostatic and spectral and has the same physics parameterisation schemes as the full version of the IFS. In terms of applications, OpenIFS is able to compute deterministic and ensemble forecasts from either real, specified, or idealised initial conditions. The project is coded in FORTRAN and in C, which is efficient for intensive and scientific computing.

3.2 Existing implementation

In OpenIFS CY43R3v2, the idealised background state implemented for the baroclinic wave test case is the one developed by Jablonowski and Williamson (2006). Originally, this background state was implemented in the full version of the IFS, and hence OpenIFS, to test the dynamical core. This background state was referred to as NTESTCASE 41 (dry case) and 42 (moist case) for the Dynamical Core Intercomparison Project (DCMIP). The original initial state of Jablonowski and Williamson (2006) has a very strong meridional temperature gradient, which means that the near-surface temperature reaches −50 °C at high latitudes. In the dry case with no physics parameterisation schemes, the surface heat fluxes are not computed, meaning that the SSTs can be specified to be much warmer (or colder) than the near-surface atmospheric temperatures without causing any problems such as destabilisation of the boundary layer or convection. In contrast, in the moist case with physics parameterisation schemes turned on, exceptionally cold conditions at high latitudes with physically realistic SSTs cause large surface heat fluxes to develop and in the extreme case can result in low-pressure centres resembling polar lows developing at high latitudes. Therefore, modifications to the Jablonowski and Williamson (2006) case are needed to enable it to be run with physics parameterisation schemes and to allow it to be used to investigate cyclone dynamics rather than the numerical accuracy of dynamical cores. Hence, the SST definition presented in Sect. 2.3 was used. Lastly, many aspects of the existing implementation are hard-coded, and the parameters used to compute the background state were not accessible via the OpenIFS namelist.

3.3 The new implementation: OpenIFS baroclinic wave v1.0

The proposed background state is implemented into OpenIFS based on the derived analytical equations for the geopotential and temperature field as detailed in Sect. 2. Both fields contain nontrivial functions – such as the gamma function (see Eq. 8) – and can be difficult to implement. In order to avoid the costly use of factorials, F₃ was expressed as a binomial coefficient fraction and all the binomial coefficients were computed once with the multiplicative method, since $\left(\genfrac{}{}{0ex}{}{z}{k+\mathrm{1}}\right)=\frac{z-k}{k+\mathrm{1}}\left(\genfrac{}{}{0ex}{}{z}{k}\right)$ with z and k being integers. By using the definition for the gamma function for positive integers z,

$$\begin{array}{}\text{(21)}& \mathrm{\Gamma }\left(z\right)=(z-\mathrm{1})\mathrm{!}\text{ and }\mathrm{\Gamma }(z+\mathrm{1}/\mathrm{2})={\displaystyle \frac{\left(\mathrm{2}z\right)\mathrm{!}}{{\mathrm{4}}^{z}z\mathrm{!}}}\sqrt{\mathit{\pi }},\end{array}$$

the gamma function can be replaced by binomial coefficients in F₃ as follows.

$$\begin{array}{ll}{\displaystyle \frac{\mathrm{\Gamma }(k+n+\mathrm{3}/\mathrm{2})}{\mathrm{\Gamma }(k+n+\mathrm{2})}}& {\displaystyle }={\displaystyle \frac{\mathrm{\Gamma }(k+n+\mathrm{1}+\mathrm{1}/\mathrm{2})}{\mathrm{\Gamma }(k+n+\mathrm{1}+\mathrm{1})}}\\ {\displaystyle }& {\displaystyle }={\displaystyle \frac{\mathrm{\Gamma }(z+\mathrm{1}/\mathrm{2})}{\mathrm{\Gamma }(z+\mathrm{1})}}\text{ where }z=k+n+\mathrm{1}\\ {\displaystyle }& {\displaystyle }={\displaystyle \frac{\left(\mathrm{2}z\right)\mathrm{!}}{{\mathrm{4}}^{z}z\mathrm{!}}}\sqrt{\mathit{\pi }}{\displaystyle \frac{\mathrm{1}}{z\mathrm{!}}}\text{ Note: }\left({\displaystyle \genfrac{}{}{0ex}{}{\mathrm{2}z}{z}}\right)={\displaystyle \frac{\left(\mathrm{2}z\right)\mathrm{!}}{z\mathrm{!}z\mathrm{!}}}\\ {\displaystyle }& {\displaystyle }=\left({\displaystyle \genfrac{}{}{0ex}{}{\mathrm{2}z}{z}}\right){\displaystyle \frac{\sqrt{\mathit{\pi }}}{{\mathrm{4}}^z}}\\ {\displaystyle }& {\displaystyle }=\left({\displaystyle \genfrac{}{}{0ex}{}{\mathrm{2}(k+n)+\mathrm{2}}{k+n+\mathrm{1}}}\right){\displaystyle \frac{\sqrt{\mathit{\pi }}}{{\mathrm{4}}^{k+n+\mathrm{1}}}}\end{array}$$

F₃ can then be rewritten as

$$\begin{array}{}\text{(22)}& \begin{array}{rl}{F}_{\mathrm{3}}=& \sum _{k=\mathrm{0}}^n\left(\genfrac{}{}{0ex}{}{n}{k}\right)(-\mathrm{1}{)}^k{\displaystyle \frac{\mathrm{1}}{\mathrm{2}(k+n)+\mathrm{1}}}\left(\genfrac{}{}{0ex}{}{\mathrm{2}(k+n)+\mathrm{2}}{k+n+\mathrm{1}}\right)\\ & \times {\displaystyle \frac{\mathit{\pi }}{{\mathrm{4}}^{k+n+\mathrm{1}}}}.\end{array}\end{array}$$

Table 1 Template of the namelist used to set the proposed background state in the moist case. All other parameters under NAEPHY were set to false.

The proposed solution has been implemented as a new idealised case (indicated by the NTESTCASE parameter in OpenIFS), where the model state variables are initialised based on the equations for geopotential, virtual temperature, the horizontal wind components, and, in the case of moist simulations, the specific humidity that were derived above. Once the initial values of the state variables are defined in the model, the OpenIFS simulations are integrated forward in time on an aquaplanet. In the moist case, most of the physics parameterisation schemes of OpenIFS are switched on, but the radiation parameterisation scheme and the wave model are deactivated. The customised SST function presented in Eq. (16) has been added to the other SST schemes already implemented in OpenIFS and is identified with the number 10 in the NAMDYNCORE namespace (variable name MSSTSCHEME in OpenIFS). In the OpenIFS namelist, NAMDYNCORE and NAEPHY are important namespaces to fill in order to ensure the correct configuration and set up of the baroclinic wave simulation. The NAMDYNCORE namespace sets up all of the idealised model configurations and NAEPHY the different physics parameterisation schemes (i.e. whether they are activated or not). Finally, the NAMCT0 namespace sets the main model control variables. N3DINI was set to 2, meaning that the meteorological values in the initial grib files were ignored and replaced by the idealised background state. The default values for the simulation in the moist case are presented in Table 1.

In this version there is no de-centring or Asselin filter. The spectral diffusion used by default is of fourth order (with the exponent of the wavenumber dependency REXPDH = 4) and is set to be rather weak, the strength of which is related to the used model time step. The coefficients for T_L319 are 2100.0 s for vorticity (HDIRVOR), divergence (HDIRDIV), temperature (HDIRT), and humidity (HDIRQ) diffusions, and the other coefficients are set to zero.

Table 2 Modifiable parameters with their default values and short description with units in NAMDIM namespace.

3.4 How to use the new implementation

Subsequent baroclinic wave simulations were run in the dry and moist case. The difference between the dry and moist case is the computation of virtual temperature (see Eq. 15) and a nonzero specific humidity. In the dry case, the computation of specific humidity is disabled and thus the virtual temperature is equal to the real temperature at all times. The dry and moist test case can be computed by setting the NTESTCASE value to 41 or 42, respectively, replacing the previous implementation de facto. The current solution allows the user to switch on or off the perturbation specified in Sect. 2.4. Moreover, the user can define the amplitude of the Gaussian hill zonal wind perturbation by changing the value of ZUP in the namelist. It would be possible to use a perturbation with a different structure, but that would require the user to modify the source code. In total, six parameters were input to create various different background states and influenced the resulting baroclinic wave; one controlled the surface roughness (α_Ch) and one triggered the initial perturbation (u_p). All parameters were included in the NAMDIM namespace in the OpenIFS namelist. A default case was defined as shown in Table 2. Of these parameters, only ZCHAR and ZUP are not used to compute the initial background state.

A stand-alone version has been developed in FORTRAN and is available on Zenodo (https://doi.org/10.5281/zenodo.10592587, Bouvier et al., 2024). This stand-alone version is divided into two parts: (1) a main programme setting all the variables to compute the zonal fields and (2) a subroutine computing the zonal fields detailed in Sect. 2.1 and 2.2. In Sect. 5, the “dry case” refers to the dry simulations without physics parameterisation schemes and the “moist case” case describes the case with physics parameterisation schemes and with moisture included.

4 Description of the experiments and diagnostics

This section is divided into four parts: (1) numerical stability of the initial background state (dry case and moist case), (2) meteorological stability and structure of the initial states, (3) temporal evolution of the default dry and moist baroclinic waves, and (4) sensitivity of the evolution to different initial background states. All simulations are run at T_L319 L137 resolution (i.e. 63 km horizontal resolution at the Equator and 137 vertical levels with a model top of 0.01 hPa; https://confluence.ecmwf.int/display/UDOC/L137+model+level+definitions, last access: 5 December 2023), with a time step of 900 s for 15 d and an output frequency of 3 h. Simulations without the perturbation were conducted to test the stability of the proposed background state: if implemented correctly and numerical errors are small, the initial state should not change in time. The jet width and height were varied by changing n and b, and all other parameters presented in Table 2 were set to their default value. In total, 18 background states without the unbalanced wind perturbation were tested for three values of n (1, 3, 6), three values of b (1, 1.5, 2), and for both the dry and moist cases.

The root mean square error (RMSE) is computed across all vertical levels for the zonal wind (Eq. 4.2, Jablonowski and Williamson, 2006) as

$$\begin{array}{}\text{(23)}& \begin{array}{rl}& \mathrm{RMSE}\left(u_{\text{za}}\right(t)-u_{\text{ideal}}(t=\mathrm{0}\left)\right)\approx \\ & {\left({\displaystyle \frac{\begin{array}{c}{\sum }_{{\mathit{\eta }}_i={\mathit{\eta }}_{\text{surface}}}^{{\mathit{\eta }}_{\text{top}}}{\sum }_{{\mathit{\varphi }}_j=-\mathrm{90}\mathit{°}}^{\mathrm{90}\mathit{°}}\left[u_{\text{za}}\right({\mathit{\varphi }}_j,{\mathit{\eta }}_i,t)\\ -u_{\text{ideal}}({\mathit{\varphi }}_j,{\mathit{\eta }}_i,t=\mathrm{0}){]}^{\mathrm{2}}\times w_{{\mathit{\varphi }}_j}\mathrm{\Delta }{\mathit{\eta }}_i\end{array}}{{\sum }_{{\mathit{\eta }}_i={\mathit{\eta }}_{\text{surface}}}^{{\mathit{\eta }}_{\text{top}}}{\sum }_{{\mathit{\varphi }}_j=-\mathrm{90}\mathit{°}}^{\mathrm{90}\mathit{°}}{w}_{{\mathit{\varphi }}_j}\mathrm{\Delta }{\mathit{\eta }}_i}}\right)}^{\frac{\mathrm{1}}{\mathrm{2}}},\end{array}\end{array}$$

where u_za is the zonal average of the zonal wind speed, u_ideal is the ideal zonal average of the zonal wind speed computed from the analytical expression for the zonal wind (Eq. 3), $w_{{\mathit{\varphi }}_j}$ represents the weights to correct the convergence of the meridians ϕ_j, and Δη_i is the thickness of the model layer η_i.

Previous studies have used several metrics to test the meteorological stability of their background states such as the temperature, geopotential height, and zonal winds (Khairoutdinov et al., 2022); sometimes potential temperature, absolute vorticity, and the Brunt–Väisälä frequency are added (Jablonowski and Williamson, 2006; Ullrich et al., 2015). The absolute vorticity, potential temperature, equivalent potential temperature (for the moist cases), zonal wind, and Brunt–Väisälä frequency were computed to test if the initial state is stable to static (gravitational), inertial, and symmetric instability. For the initial state to be absolutely stable to dry and saturated vertical displacements (static stability), equivalent potential temperature must increase with height everywhere. In the situation where equivalent potential temperature decreases with height, conditional instability is present, meaning that the atmosphere is stable to displacements of dry and unsaturated air parcels but unstable to displacements of saturated air parcels. If potential temperature decreases with height, then the atmosphere is absolutely unstable – both dry and saturated displacements are unstable. Thus, for the initial state to be absolutely stable potential temperature must increase with height and the Brunt–Väisälä frequency must be positive. Regions where equivalent potential temperature decreases with height are also stable and acceptable in the initial state as long as these regions are not saturated. For the initial state to be stable to horizontal displacements (inertial stability) the absolute vorticity must be positive (negative) in the Northern (Southern) Hemisphere. Situations can exist where the atmosphere is statically and inertially stable, but the atmosphere is unstable to slantwise displacements (symmetric instability). This exists when the potential vorticity is negative. Hence for our initial state to be stable to inertial and symmetric instability, both the absolute vorticity and potential vorticity must be positive in the Northern Hemisphere. However, although the initial state must be stable to static, inertial, and symmetric stability, the setup must be baroclinically unstable. This requires the presence of a meridional temperature gradient in the mid-latitudes and a well-defined zonal jet.

Several baroclinic waves (background states with perturbation) are generated for six values of n (1 to 6) and the default values for the remaining parameters (presented in Table 2). All baroclinic wave simulations (BWSs) have been run for 15 d. During the 15 d simulation, the first cyclone to emerge develops directly from the initial perturbation; however, upstream and downstream development also occurs, resulting in multiple cyclones and anticyclones. To objectively identify the centre of the cyclone which develops directly from the initial perturbation (the first cyclone) the TRACK software (Hodges, 1994, 1995, 1999) is used. Cyclones are identified as localised maxima in the 850 hPa relative vorticity field truncated to T42 spectral resolution. Only tracks lasting for 4 d and which travel 1000 km are retained. The weak tracks are removed by setting a threshold of 0.4 × 10⁻⁵ s⁻¹ for the T42 850 hPa relative vorticity. This threshold is weaker than the threshold usually applied when using TRACK (1 × 10⁻⁵ s⁻¹) to identify cyclones in the real world to enable the first cyclone to be detected as soon as possible.

Figure 3 RMSE as a function of time for all cases. The solid lines represent the moist cases and the dotted lines the dry cases. The colour map is from Crameri et al. (2020).

5 Results

5.1 Numerical stability of the initial background state

As shown in Fig. 3, the RMSEs for all of the dry cases with no perturbation included are below 0.005 m s⁻¹ with no noticeable increase over time. These values are lower than the RMSE in the Jablonowski and Williamson (2006) report for the same semi-Lagrangian setup and lower than the RMSE reported by Khairoutdinov et al. (2022). This means that the specified background states are stable and correctly implemented into OpenIFS. A slight increase in RMSE over time can be observed in the moist cases; however, in comparison to the dry cases the RMSE in the moist cases is at most 6 % larger by the end of the 15 d simulation. The lower the n value and the higher the b value, the higher the RMSE. This tendency can be understood as follows: the more extensive the jet, the higher the RMSE.

Figure 4 Cross-sections of the default moist initial background state (n=3 and b=2.0) for (a) zonal wind speed, (b) absolute vorticity, (c) potential temperature, (d) Brunt–Väisälä frequency, (e) equivalent potential temperature, and (f) specific humidity.

Figure 5 The zonal wind speed in m s⁻¹ (black contours), temperature in K (red contours), and dynamical tropopause at 2 PVU (blue contours) fields of the moist initial state for different values of n and b.

5.2 Meteorological stability and structure of the initial background state

The initial fields for the default values (n=3 and b=2.0) in the moist case are presented in Fig. 4. The initials fields for the default values in the dry case are presented in the Supplement (Fig. S4) as they are very similar to the moist case with the obvious exceptions of the specific humidity (zero everywhere in the dry case) and the equivalent potential temperature. In both the dry and moist cases, a strong baroclinic zone is present between 30 and 60° (Figs. 4a, c and 5d) with temperatures at the surface in the moist (dry) case reaching a maximum of 22.6 °C (26.0 °C) at the Equator and decreasing to −13.7 °C (−12.8 °C) at the poles. These temperatures are similar to the temperatures reported in the ERA5 global reanalysis, which display a surface temperature between 20 °C and 30 °C in the tropics and between −10 °C and −20 °C at the poles (Hersbach et al., 2020). The baroclinic zone is co-located with a jet stream between 20° N/S and 70° N/S with a maximum zonal wind speed of 30 m s⁻¹ (identical in the dry case) located in the centre of the jet at 45° N/S and 250 hPa as shown in Fig. 4a and c. When the width and height of the idealised jet streams that we propose in this study are compared to realistic values from ERA5 (see Lee et al., 2023), reasonable agreement is found with the climatological winter means over the northwestern Atlantic. However, the maximum zonal wind speed is in the lower range of the maximum wind speed of the ERA5 profile, which is not a problem considering that u₀ is a tunable parameter as presented in Table 2. Moreover, the December mean zonal wind speed over the years 1979–2021 was computed (Fig. S1 in the Supplement) and compared to our idealised jet structures (Fig. S2 in the Supplement). A 10° location difference of the centre of the northern jet has been found, but the proposed analytical solution is close in shape and intensity to this profile. The absolute vorticity, presented in Fig. 4b, demonstrates that the default initial state is inertially stable. The condition for symmetric stability, which is that potential vorticity is positive in the Northern Hemisphere and negative in the Southern Hemisphere, is also met. This can be inferred from the potential temperature distribution (Fig. 4c) combined with the absolute vorticity distribution (Fig. 4b). The potential temperature (Fig. 4c) increases with height everywhere in the model domain, indicating that the initial state is stable to dry static stability. Figure 4e shows that the equivalent potential temperature increases everywhere, except for the tropical regions. This means that the initial state is stable to moist static stability, as long as air parcels are not saturated in the tropics. The static stability is also indicated by the Brunt–Väisälä frequency presented in Fig. 4d. Relatively high values for the Brunt–Väisälä frequency were seen at high latitudes at most pressure levels, indicating high static stability. A lower Brunt–Väisälä frequency and thus lower static stability were observed around the Equator and near the surface in the polar regions.

Figure 5 shows the initial temperature, zonal wind speed, and dynamical tropopause (taken here to be the 2 PVU surface) for different values of n (1, 3, and 6) and b (1, 1.5, and 2), with Fig. 5d showing the default initial state which was presented in Fig. 4. The same figure has been produced for the dry cases (Fig. S5 in the Supplement), but this figure is not included here because of the similarity to Fig. 5. Increasing n causes the meridional width of the jet to decrease. In the case of n=1, the jet is very wide, extending from 15° N/S to 75° N/S. In contrast, when n=6, the jet is constrained between 30 and 60° N/S. Increasing n also causes the width of the baroclinic zone to decrease and the surface temperature gradient between the Equator and pole to decrease: in the case of n=1, the temperature difference between the Equator and pole is 57.3 °C, whereas for n=6 the temperature difference is 26.3 °C (for b=2.0). Decreasing b leads to the displacement of the centre of the jet toward the surface. For all initial states, higher values for n lead to a narrower jet and a stronger temperature gradient (Fig. 5). Moreover, high values of b lead to a stronger temperature gradient, a higher centre of the jet, and an increase in the maximum wind speed of the jet as shown by the innermost contour (e.g. Fig. 5a and c). Additionally, the centre of the jet gets closer to the surface with decreasing b. For b=1, the centre of the jet is lower than the dynamical tropopause at 2 PVU but does not resemble a realistic jet structure as found in reanalysis (Lee et al., 2023, and Figs. S1 and S2). Thus, when b=1, the jet's width is unrealistic, meaning that it is not recommended to use b=1 when studying extratropical cyclone dynamics. The latitudinal location of the centre of the jet is independent of n and b and located at 45° in both hemispheres for all the different initial cases presented. The height of the dynamical tropopause increases with increasing n, which means that the static stability increases with higher n (Held, 1982).

Based on the results of the Sect. 5.1 and 5.2, it is evident that the initial background states are stable through time and balanced from a meteorological perspective. As shown in Fig. 5, the background states display strong baroclinicity in the mid-latitudes. The next section will present the evolution of the baroclinic waves triggered in some of the presented background states (b=2.0).

Figure 6 The development of the baroclinic wave for the default dry scenario (n=3 and b=2.0) at (a) t = 144 h, (b) t = 168 h, (c) t = 192 h, (d) t = 228 h, (e) t = 264 h, and (f) t = 312 h. The black contours show mean sea level pressure (hPa), and the shading shows the temperature (°C) at 850 hPa. Note: this figure does not show the whole model domain; the x axis ranges from 40–220° E, while the y axis ranges from 30–65° N.

5.3 Temporal evolution of the dry and moist default baroclinic waves

5.3.1 Dry case

The evolution of the dry default baroclinic wave (n=3, b=2), which is also run with no physics parameterisation schemes switched on, is shown in Fig. 6. The formation of a very weak closed low-pressure system is just evident after 144 h (Fig. 6a). At this time, the minimum mean sea level pressure (MSLP) is located slightly poleward of 45° N, which was the central latitude of the initial perturbation. Also at 144 h, a small-amplitude wave in the 850 hPa temperature is evident, co-located with the developing cyclonic circulation, but pronounced fronts have not yet developed. After 168 h of development (Fig. 6b), the cyclone has a minimum MSLP of 1006.75 hPa, which is 6.5 hPa lower than the initial surface pressure. Two areas of high pressure are evident on either side of this cyclone, and another cyclone has begun developing upstream of the cyclone which developed directly from the initial unbalanced perturbation. By 192 h (Fig. 6c), a well-developed cyclone is present with a minimum MSLP of 1000.8 hPa. Cold and warm fronts are now evident in the 850 hPa temperature. The centre of the cyclone has continued to move polewards and the two anticyclones on either side have also intensified slightly and moved equatorwards. In addition to the upstream development that was evident 24 h previously, downstream development has also started to take place by 192 h; a third low-pressure centre/cyclone is now visible at 150° E, although this is the weakest of all three cyclones. Rapid intensification takes place and by 228 h (Fig. 6d), and the initial cyclone has deepened considerably (minimum MSLP of 983.3 hPa and almost 30 hPa lower than the initial surface pressure). Cold and warm fronts with large thermal gradients are present and an occlusion has begun to develop. Furthermore, at 228 h, both the upstream and downstream cyclones have continued deepening, with the upstream cyclone being deeper yet spatially smaller than the downstream cyclone. The horizontal distance between the upstream cyclone and the initial cyclone is smaller than between the downstream cyclone and the initial cyclone. By 264 h (Fig. 6e), the initial cyclone is very mature and deep, has a clear warm seclusion, and has moved farther polewards. The downstream cyclone has undergone notable deepening since 228 h and a fourth relatively small cyclone has started to develop upstream at 60° E. After 312 h of simulation, a much more chaotic picture is evident. The first upstream cyclone has merged with the trailing cold front of the initial cyclone, leading to the presence of a very large cyclonic system which has a very strong north–south pressure gradient on its southern flank. The downstream cyclone has continued intensifying and has moved polewards, leading to a very meridionally extended system that has strong fronts associated with it. The upstream cyclone has also continued to intensify but remains weaker and smaller in spatial scale than the other cyclones.

Figure 7 The development of the baroclinic wave for the default moist scenario (n=3 and b=2.0) at (a) t = 144 h, (b) t = 168 h, (c) t = 192 h, (d) t = 228 h, (e) t = 264 h, and (f) t = 312 h. The black contours show mean sea level pressure (hPa), and the shading shows the temperature (°C) at 850 hPa. Note: this figure does not show the whole model domain; the x axis ranges from 40–220° E, while the y axis ranges from 30–65° N.

Figure 8 The development of the baroclinic wave for the default moist scenario (n=3 and b=2.0) at (a) t = 144 h, (b) t = 168 h, (c) t = 192 h, (d) t = 228 h, (e) t = 264 h, and (f) t = 312 h. The black contours show mean sea level pressure (hPa), and the shading shows the 3 h precipitation (in mm). Note: this figure does not show the whole model domain; the x axis ranges from 40–220° E, while the y axis ranges from 30–65° N.

5.3.2 Moist case

The evolution of the default baroclinic wave (n=3, b=2) with moisture is shown in Fig. 7 and the associated precipitation patterns in Fig. 8. After 144 h of development (Fig. 7a), the structure of the moist baroclinic wave is very similar to the dry case at the same time; only a very weak cyclonic circulation and no precipitation have developed at this time (Fig. 8a). At 168 h (Fig. 7b), the initial cyclone has deepened to 1007.5 hPa, which is slightly weaker than the dry case. A thermal wave is also evident at 168 h and an upstream cyclone has begun to develop but no precipitation has developed at this stage (Fig. 8b). At 192 h, the general evolution of the moist baroclinic continues in a similar manner to the dry case. A well-developed cyclone is now present with a minimum MSLP below 1003 hPa, which is again slightly higher than in the dry case. Cold and warm fronts are now evident in the 850 hPa temperature and precipitation is now evident on the warm front of the initial cyclone and, to a lesser extent, in the developing warm sector of the upstream cyclone. Rapid development takes places between 192 and 228 h, and by 228 h, a strong cyclone with clear fronts is evident. Furthermore, by 228 h moderate to heavy precipitation (values exceeding 3 mm per 3 h) is present along the warm fronts of the initial cyclone and the upstream cyclone as well as on the warm side of the cold front and in the warm sector these two cyclones. In the moist case, after 264 h, three well-developed cyclones are evident, all of which have a surface pressure trough co-located with the cold front (Fig. 7e). In comparison to the dry case, the moist case has much weaker horizontal pressure gradients, particularly on the southern and northern sides of the initial cyclone. Precipitation now covers a larger area than 24 h earlier and extends further south along the cold fronts of each of these three well-developed cyclones (Fig. 8e). The second and weaker upstream cyclone also produces precipitation at this stage. After 312 h of development, the initial cyclone and the first upstream cyclone have become very mature systems. The strongest cyclone at this stage is the downstream cyclone, which has a minimum mean sea level pressure below 972 hPa (41 hPa lower than the initial surface pressure). Precipitation is still associated with all four cyclones (Fig. 8f). The downstream cyclone has the heaviest and most expansive area of precipitation, whereas the first upstream cyclone only has moderate precipitation remaining on the trailing cold front. In comparison to the dry case at 312 h, all cyclones in the moist case are weaker and less developed, and the first upstream cyclone and initial cyclone remain as separate features.

The slower rate of development in the moist case is very likely caused by the presence of surface friction (and other physics parameterisation schemes) in the moist simulations, whereas the dry cases are run with no friction and no parameterised physics schemes. This results is in agreement with Boutle et al. (2010), who, in idealised baroclinic life cycle simulations, found that a deeper cyclone developed in their dry no-boundary-layer case (i.e. no boundary layer scheme nor surface friction was included in the simulations) than in their moist simulation in which the boundary layer scheme was activated.

Overall, Figs. 6–8 demonstrate that the cyclones and anticyclones that develop from the default initial state are realistic and consequently can be used to investigate cyclone dynamics. In comparison to previous BWS studies that have used limited model domains with periodic boundaries (Feldstein and Held, 1989; Hoskins et al., 1977), the setup presented here and the evolution of the default baroclinic wave allow for studies of both upstream and downstream development.

Figure 9 The development of the baroclinic wave at t = 204 h for (a) n=1, (b) n=2, (c) n=3, (d) n=4, (e) n=5, and (f) n=6 in the dry case (for fixed b=2.0). The black contours show mean sea level pressure (hPa), and the shading shows the temperature (°C) at the 850 hPa level. Note: this figure does not show the whole model domain; the x axis ranges from 40–220° E, while the y axis ranges from 30–65° N.

Figure 10 The development of the baroclinic wave at t = 204 h for (a) n=1, (b) n=2, (c) n=3, (d) n=4, (e) n=5, and (f) n=6 in the moist case (for fixed b=2.0). The black contours show mean sea level pressure (hPa), and the shading shows the temperature (°C) at the 850 hPa level. Note: this figure does not show the whole model domain; the x axis ranges from 40–220° E, while the y axis ranges from 30–65° N.

5.4 Impact of different initial states on the baroclinic wave evolution

Figures 9 and 10 demonstrate how the structure of the dry and moist baroclinic wave depends on the value of n after 204 h. In both the dry and moist case, decreasing n increases the width of the jet stream and the baroclinic zone (Fig. 5), which means that the resultant cyclones have a greater meridional extent with smaller values of n (Figs. 9 and 10). Decreasing n leads to a reduction in the minimum surface pressure of the cyclones in both the dry and moist cases. In the dry case, for n=1, the minimum surface pressure of the first developing cyclone after 204 h is 989 hPa (Fig. 9a), whereas for n=6 the corresponding value is 1000 hPa (Fig. 9f). The decrease in surface pressure with increasing n also holds true for the cyclones which develop upstream and downstream and is likely because the 850 hPa temperature difference between the Equator and the pole, and hence available potential energy and baroclinicity, is larger with smaller n (Fig. 10). Decreasing n also causes the phase speed of the cyclones to increase and the cyclones and anticyclones travel eastward faster with larger n, despite the fact that the maximum speed of the zonal jet does not change considerably. This increase in phase speed is observed in both the dry and moist cases. A detailed investigation of how all parameters which control the structure of the initial state (as presented in Table 2) affect the intensity and structure of the resultant cyclones will be the topic of a future study.

Figure 11 Evolution of the maximum relative vorticity of the first cyclone to develop as part of the baroclinic wave as a function of time. Solid lines represent the moist cases and the dashed lines the dry cases (for fixed b=2.0). The black horizontal line shows the threshold for relative vorticity set to 0.4 × 10⁻⁵ s⁻¹ used by the tracking algorithm TRACK. The colour map is from Crameri et al. (2020).

Figure 11 shows the temporal evolution of the maximum 850 hPa relative vorticity as identified by TRACK of the cyclone which develops first in both the dry and moist experiments with different n. Note that this is not necessarily the cyclone which experiences the largest deepening rates or the largest maximum vorticity, but it is the cyclone which is directly caused by the initial unbalanced perturbation. The first cyclone in all experiments is detected by TRACK after 5–6 d. Between day 7 and day 10, a period of rapid intensification occurs for all experiments. In both the dry and moist cases, as n decreases, the rate of intensification increases, particularly between day 8 and day 10, and the eventual maximum value of relative vorticity is larger for smaller n. This is consistent with Figs. 9 and 10, which show that the minimum surface pressure was also lower with smaller n and that the 850 hPa temperature difference between the poles and the Equator was larger with smaller n. Furthermore, increasing n results in the maximum vorticity being reached at an earlier time in both the dry and moist cases. Figure 11 also clearly highlights the difference between the moist and dry simulations. The maximum vorticity of the first developing cyclone is larger and occurs later in time in the dry cases compared to the corresponding moist cases.

6 Conclusions

This article introduced idealised initial background states for baroclinic life cycle simulations. The main advantages of these background states are that they can entirely be expressed analytically and controlled through configuration files. The jet structure and strength can be tuned, as can the average virtual temperature, the surface relative humidity, the lapse rate and the surface roughness. This flexibility allows an easy generation of different background states and their related baroclinic waves. All studied initial background states are proven to be stable even in the moist case scenario. Moreover, a Gaussian perturbation of the zonal wind speed allows the development of a baroclinic wave in the dry and moist cases, which depends on the given jet structure. The presented solution is appealing for two main reasons.

First, the proposed solution is implemented in OpenIFS CY43R3, which is a popular state-of-the-art model used extensively for meteorological and climate research (Carver, 2022). The idea was to propose a solution easy to test and use it to allow wide accessibility. Moreover, the Appendix presents the exhaustive integration and derivation of the initial background state to allow easy modification of the analytical formulae of the zonal wind speed field and the easy identification of the potential difficulties of its integration. One limitation of our initial state is that in its current form it is not possible to easily add barotropic shear to the low level of the jet. Barotropic shear has proven to be a defining condition for the structure development of extratropical cyclones (Agustí-Panareda et al., 2005). Future work will address this issue by proposing an analytical solution including an analytical barotropic shear.

Second, even if the initial background states may be unstable to moist saturated air parcel displacements, for our default simulation of n=3 and b=2, no convective available potential energy (CAPE) is present in the tropics (0° N). This shows that even when the air parcel at the surface is lifted, it will stay cooler than its environment, meaning that the atmosphere is stable. In the case of n=1 and b=1, the CAPE is large due to a stronger decrease in temperature with height in the troposphere. However, as this CAPE can only be released once the surface parcel reaches saturation or is substantially lifted and this area is not affected by the baroclinic wave, this does not influence the meteorological stability of our setup. The initial background states are proven to be realistic considering ERA5 average temperature and zonal wind speed cross-sections. Moreover, this study has proven the sensitivity of the BWS to the height and width of the jet, which has been possible by the control of the different parameters through the OpenIFS namelist. The proposed solution is of interest to create a large ensemble of baroclinic life cycles which would allow the study of the sensitivity of ETC to various initial background states. Future study will investigate the dependency of several measures of cyclone intensities – such as mean sea level pressure, relative vorticity at 850 hPa, 10 m wind gust, and storm severity index (SSI) (Leckebusch et al., 2008) – on the initial conditions of cyclogenesis.

Appendix A: Detailed derivation of analytical initial conditions

The derivation of the analytical initial conditions starts from the primitive equations for moist adiabatic and frictionless flow in spherical coordinates (λ,ϕ) and vertical η levels by the equations for the following.

$$\begin{array}{}\text{(A1a)}& {\displaystyle }\begin{array}{rl}u\text{momentum}:& {\displaystyle \frac{\mathrm{d}u}{\mathrm{d}t}}-{\displaystyle \frac{uv\tan \mathit{\varphi }}{a}}=-{\displaystyle \frac{\mathrm{1}}{a\cos \mathit{\varphi }}}\\ & \times \left({\displaystyle \frac{\partial \mathrm{\Phi }}{\partial \mathit{\lambda }}}+R_{\mathrm{d}}{T}_{\mathrm{v}}{\displaystyle \frac{\partial \ln p}{\partial \mathit{\lambda }}}\right)+fv\end{array}\text{(A1b)}& {\displaystyle }\begin{array}{rl}v\text{momentum}:& {\displaystyle \frac{\mathrm{d}v}{\mathrm{d}t}}+{\displaystyle \frac{u^{\mathrm{2}}\tan \mathit{\varphi }}{a}}=-{\displaystyle \frac{\mathrm{1}}{a}}\\ & \left({\displaystyle \frac{\partial \mathrm{\Phi }}{\partial \mathit{\varphi }}}+R_{\mathrm{d}}{T}_{\mathrm{v}}{\displaystyle \frac{\partial \ln p}{\partial \mathit{\lambda }}}\right)-fu\end{array}\text{(A1c)}& {\displaystyle }\begin{array}{r}\text{hydrostatic balance}:{\displaystyle \frac{\partial \mathrm{\Phi }}{\partial \mathit{\eta }}}=-{\displaystyle \frac{R_{\mathrm{d}}{T}_{\mathrm{v}}}{p}}{\displaystyle \frac{\partial p}{\partial \mathit{\eta }}}\end{array}\text{(A1d)}& {\displaystyle }\begin{array}{rl}\text{continuity}:& {\displaystyle \frac{\partial }{\partial t}}\left({\displaystyle \frac{\partial p}{\partial \mathit{\eta }}}\right)+{\displaystyle \frac{\mathrm{1}}{a\cos \mathit{\varphi }}}{\displaystyle \frac{\partial }{\partial \mathit{\lambda }}}\left(u{\displaystyle \frac{\partial p}{\partial \mathit{\eta }}}\right)\\ & +{\displaystyle \frac{\mathrm{1}}{a\cos \mathit{\varphi }}}{\displaystyle \frac{\partial }{\partial \mathit{\varphi }}}\left(\left(v\cos \mathit{\varphi }\right){\displaystyle \frac{\partial p}{\partial \mathit{\eta }}}\right)\\ & +{\displaystyle \frac{\partial }{\partial \mathit{\eta }}}\left(\dot{\mathit{\eta }}{\displaystyle \frac{\partial p}{\partial \mathit{\eta }}}\right)=\mathrm{0},\end{array}\text{(A1e)}& {\displaystyle }\begin{array}{r}\text{thermodynamics}:{\displaystyle \frac{\mathrm{d}T}{\mathrm{d}t}}-{\displaystyle \frac{R_{\mathrm{d}}{T}_{\mathrm{v}}\mathit{\omega }}{c_{\mathrm{p}}p}}=\mathrm{0}\end{array}\text{(A1f)}& {\displaystyle }\begin{array}{r}\text{moist ideal gas law: }p=\mathit{\rho }{R}_{\mathrm{d}}{T}_{\mathrm{v}}\end{array}\end{array}$$

Here, $\mathrm{d}/\mathrm{d}t$ is the full time derivative in spherical coordinates given as

$${\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\left(\right)={\displaystyle \frac{\partial }{\partial t}}+{\displaystyle \frac{u}{a\cos \mathit{\theta }}}{\displaystyle \frac{\partial }{\partial \mathit{\lambda }}}\left(\right)+{\displaystyle \frac{v}{a}}{\displaystyle \frac{\partial }{\partial \mathit{\varphi }}}\left(\right)+\dot{\mathit{\eta }}{\displaystyle \frac{\partial }{\partial \mathit{\eta }}}\left(\right).$$

Please note that the specific heat capacity of air, c_p, in the thermodynamic equation needs to be corrected with a correction factor δ when moisture is included in the model. The correction factor is defined as

$${\displaystyle }\mathit{\delta }=\mathrm{1}+(c_{\text{p,vap}}/c_{\text{p, dry}}-\mathrm{1})q,$$

where $c_{\text{p,vap}}/c_{\text{p, dry}}$ = 1860/1004 (units: J (kg K)⁻¹) and q is the specific humidity (units: kg kg⁻¹) (ECMWF, 2017a, Eq. 2.3).

The following equalities are used in the derivations of the geopotential field.

$$\begin{array}{c}{}^{\ast }\sin \left(\mathrm{2}\mathit{\varphi }\right)=\mathrm{2}\sin \mathit{\varphi }\cos \mathit{\varphi }\\ {}^{\ast \ast }{\sin }^{\mathrm{2}}\mathit{\varphi }=\mathrm{1}-{\cos }^{\mathrm{2}}\mathit{\varphi }\\ {}^{\ast \ast \ast }(a+b{)}^n=\sum _{k=\mathrm{0}}^n\left({\displaystyle \genfrac{}{}{0ex}{}{n}{k}}\right)a^{n-k}{b}^k\\ {\displaystyle }+\tan \mathit{\varphi }={\displaystyle \frac{\sin \mathit{\varphi }}{\cos \mathit{\varphi }}}\\ {\displaystyle }\begin{array}{rl}\mathit{†}\underset{-\mathit{\pi }/\mathrm{2}}{\overset{\mathit{\pi }/\mathrm{2}}{\int }}{\cos }^{\mathrm{2}(k+n+\mathrm{1})}\mathit{\varphi }\mathrm{d}\mathit{\varphi }=& \sqrt{\mathit{\pi }}{\displaystyle \frac{\mathrm{\Gamma }(k+n+\mathrm{3}/\mathrm{2})}{\mathrm{\Gamma }(k+n+\mathrm{2})}},\\ & \text{ see e.g. WolframAlpha}\end{array}\\ {\displaystyle }\begin{array}{rl}\mathit{‡}\underset{-\mathit{\pi }/\mathrm{2}}{\overset{\mathit{\pi }/\mathrm{2}}{\int }}{\sin }^{\mathrm{2}(k+\mathrm{2}n+\mathrm{1})}\mathit{\varphi }\cos \mathit{\varphi }\mathrm{d}\mathit{\varphi }=& \underset{-\mathrm{1}}{\overset{\mathrm{1}}{\int }}{u}^{\mathrm{2}(k+\mathrm{2}n+\mathrm{1})}\mathrm{d}u\\ =& {\displaystyle \frac{\mathrm{2}}{\mathrm{2}(\mathrm{2}n+k+\mathrm{1})+\mathrm{1}}}.\end{array}\end{array}$$

Whenever an equality is used its symbol is noted above the equality sign, for example $\stackrel{\ast }{=}$.

A1 Mean geopotential field

The mean virtual temperature field is defined as

$$\begin{array}{}\text{(A2)}& \langle T_{\mathrm{v}}\rangle =T_{v,\mathrm{0}}-\mathit{\gamma }z,\end{array}$$

where T_v,0 is the reference virtual temperature, γ the lapse rate, and z height (moist version of Eq. 10,  Ullrich et al., 2015). The geopotential mean field is derived from the mean virtual temperature field and the hydrostatic balance equation as follows.

$$\begin{array}{ll}{\displaystyle \frac{\partial p}{\partial z}}& {\displaystyle }=-\mathit{\rho }g\\ {\displaystyle \frac{\partial p}{\partial z}}& {\displaystyle }=-{\displaystyle \frac{p}{R_{\mathrm{d}}\langle T_{\mathrm{v}}\rangle }}|\text{ Eq. (A2)}\\ {\displaystyle \frac{\partial p}{p}}& {\displaystyle }=-{\displaystyle \frac{g}{R_{\mathrm{d}}(T_{v,\mathrm{0}}-\mathit{\gamma }z)}}\partial z|\text{ integrate LH and RH }\\ \text{(A3)}& {\displaystyle }\underset{p_{\mathrm{s}}}{\overset{p}{\int }}{\displaystyle \frac{\partial p}{p}}& {\displaystyle }=-{\displaystyle \frac{g}{R_{\mathrm{d}}}}\underset{z_{\mathrm{0}}=\mathrm{0}}{\overset{z}{\int }}{\displaystyle \frac{\mathrm{1}}{(T_{v,\mathrm{0}}-\mathit{\gamma }z)}}\partial z\end{array}$$

The left-hand side of Eq. (A3) is solved as

$$\begin{array}{}\text{(A4)}& \underset{p_{\mathrm{s}}}{\overset{p}{\int }}{\displaystyle \frac{\partial p}{p}}=\ln {\displaystyle \frac{p}{p_{\mathrm{s}}}},\end{array}$$

and the right-hand side of Eq. (A3) is solved by substitution of $u=T_{v,\mathrm{0}}-\mathit{\gamma }z$. The right-hand side is solved as

$$\begin{array}{}\text{(A5)}& -{\displaystyle \frac{g}{R_{\mathrm{d}}}}\underset{z=\mathrm{0}}{\overset{z}{\int }}{\displaystyle \frac{\mathrm{1}}{(T_{v,\mathrm{0}}-\mathit{\gamma }z)}}\partial z={\displaystyle \frac{g}{R_{\mathrm{d}}\mathit{\gamma }}}\ln {\displaystyle \frac{Tv,\mathrm{0}-\mathit{\gamma }z}{T_{v,\mathrm{0}}}}.\end{array}$$

By solving for z, this gives us

$$\begin{array}{}\text{(A6)}& z={\displaystyle \frac{T_{v,\mathrm{0}}}{\mathit{\gamma }}}\left[\mathrm{1}-{\left({\displaystyle \frac{p}{p_{\mathrm{s}}}}\right)}^{\frac{R_{\mathrm{d}}g}{\mathit{\gamma }}}\right],\end{array}$$

and by realising that $p/p_{\mathrm{s}}$ is η, the mean geopotential field is then

$$\begin{array}{}\text{(A7)}& \langle \mathrm{\Phi }\rangle ={\displaystyle \frac{g{T}_{v,\mathrm{0}}}{\mathit{\gamma }}}\left(\mathrm{1}-{\mathit{\eta }}^{\frac{R_{\mathrm{d}}g}{\mathit{\gamma }}}\right).\end{array}$$

A2 Geopotential field

The geopotential anomaly field is derived from the steady-state momentum equation for v with the zonal flow defined by $u=-u_{\mathrm{0}}\ln \mathit{\eta }\exp (-[\ln \mathit{\eta }/b{]}^{\mathrm{2}}){\sin }^{\mathrm{2}n}\mathrm{2}\mathit{\varphi }$:

$$\begin{array}{}\text{(A8)}& \begin{array}{rl}{\displaystyle \frac{\mathrm{1}}{a}}{\displaystyle \frac{\partial {\mathrm{\Phi }}^{\prime }}{\partial \mathit{\varphi }}}=& -\left(-u_{\mathit{\eta }}{\sin }^{\mathrm{2}n}\mathrm{2}\mathit{\varphi }\right)\\ & \times \left(\mathrm{2}\mathrm{\Omega }\sin \mathit{\varphi }+{\displaystyle \frac{-u_{\mathit{\eta }}{\sin }^{\mathrm{2}n}\mathrm{2}\mathit{\varphi }}{a}}\tan \mathit{\varphi }\right),\end{array}\end{array}$$

where $u_{\mathit{\eta }}=u_{\mathrm{0}}\ln \mathit{\eta }\exp (-[\ln \mathit{\eta }/b{]}^{\mathrm{2}})$. Equation (A8) is multiplied by a and then integrated analytically over ϕ:

$$\begin{array}{}\text{(A9)}& \begin{array}{rl}{\mathrm{\Phi }}^{\prime }(\mathit{\lambda },\mathit{\varphi },\mathit{\eta })=& \int \left[u_{\mathit{\eta }}{\sin }^{\mathrm{2}n}\mathrm{2}\mathit{\varphi }\right(\mathrm{2}a\mathrm{\Omega }\sin \mathit{\varphi }\\ & -u_{\mathit{\eta }}{\sin }^{\mathrm{2}n}\mathrm{2}\mathit{\varphi }\tan \mathit{\varphi }\left)\right]\mathrm{d}\mathit{\varphi }+{\mathrm{\Phi }}_{\mathrm{0}}\left(\mathit{\eta }\right),\end{array}\end{array}$$

which is solved by dividing the integral into two parts that are solved separately.

$$\begin{array}{}\text{(A10a)}& {\displaystyle }\begin{array}{rl}\int u_{\mathit{\eta }}{\sin }^{\mathrm{2}n}\mathrm{2}\mathit{\varphi }& \mathrm{2}a\mathrm{\Omega }\sin \mathit{\varphi }\mathrm{d}\mathit{\varphi }=\\ & u_{\mathit{\eta }}\mathrm{2}a\mathrm{\Omega }\int {\sin }^{\mathrm{2}n}\mathrm{2}\mathit{\varphi }\sin \mathit{\varphi }\mathrm{d}\mathit{\varphi }\end{array}\text{(A10b)}& {\displaystyle }\begin{array}{rl}\int -u_{\mathit{\eta }}^{\mathrm{2}}{\sin }^{\mathrm{2}n}\mathrm{2}\mathit{\varphi }& {\sin }^{\mathrm{2}n}\mathrm{2}\mathit{\varphi }\tan \mathit{\varphi }\mathrm{d}\mathit{\pi }=\\ & -u_{\mathit{\eta }}^{\mathrm{2}}\int {\sin }^{\mathrm{4}n}\mathrm{2}\mathit{\varphi }\tan \mathit{\varphi }\mathrm{d}\mathit{\varphi }\end{array}\end{array}$$

Equation (A10a) is solved as

$$\begin{array}{ll}{\displaystyle }\int {\sin }^{\mathrm{2}n}\left(\mathrm{2}\mathit{\varphi }\right)\sin \mathit{\varphi }\mathrm{d}\mathit{\varphi }& {\displaystyle }\stackrel{\ast }{=}\int (\mathrm{2}\sin \mathit{\varphi }\cos \mathit{\varphi }{)}^{\mathrm{2}n}\sin \mathit{\varphi }\mathrm{d}\mathit{\varphi }\\ {\displaystyle }& {\displaystyle }={\mathrm{2}}^{\mathrm{2}n}\int {\sin }^{\mathrm{2}n}\mathit{\varphi }{\cos }^{\mathrm{2}n}\mathit{\varphi }\sin \mathit{\varphi }\mathrm{d}\mathit{\varphi }\\ {\displaystyle }& {\displaystyle }\stackrel{\ast \ast }{=}{\mathrm{4}}^n\int (\mathrm{1}-{\cos }^{\mathrm{2}}\mathit{\varphi }{)}^n{\cos }^{\mathrm{2}n}\mathit{\varphi }\sin \mathit{\varphi }\mathrm{d}\mathit{\varphi },\end{array}$$

which is integrated by substituting u=cos ϕ and $\mathrm{d}u=-\sin \mathit{\varphi }$:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathrm{4}}^n\int (\mathrm{1}-{\cos }^{\mathrm{2}}\mathit{\varphi }{)}^n{\cos }^{\mathrm{2}n}\mathit{\varphi }\sin \mathit{\varphi }\mathrm{d}\mathit{\varphi }\\ {\displaystyle }& {\displaystyle }=-{\mathrm{4}}^n\int (\mathrm{1}-u^{\mathrm{2}}{)}^n{u}^{\mathrm{2}n}\mathrm{d}u\\ {\displaystyle }& {\displaystyle }\stackrel{\ast \ast \ast }{=}-{\mathrm{4}}^n\int \left(\sum _{k=\mathrm{0}}^n\left({\displaystyle \genfrac{}{}{0ex}{}{n}{k}}\right)(-\mathrm{1}{)}^k{u}^{\mathrm{2}(k+n)}\right)\mathrm{d}u\\ {\displaystyle }& {\displaystyle }=-{\mathrm{4}}^n\sum _{k=\mathrm{0}}^n\left({\displaystyle \genfrac{}{}{0ex}{}{n}{k}}\right)(-\mathrm{1}{)}^k\int u^{\mathrm{2}(k+n)}\mathrm{d}u\\ {\displaystyle }& {\displaystyle }=-{\mathrm{4}}^n\sum _{k=\mathrm{0}}^n\left({\displaystyle \genfrac{}{}{0ex}{}{n}{k}}\right)(-\mathrm{1}{)}^k{\displaystyle \frac{\mathrm{1}}{\mathrm{2}(k+n)+\mathrm{1}}}{u}^{\mathrm{2}(k+n)+\mathrm{1}}\\ {\displaystyle }& {\displaystyle }=-{\mathrm{4}}^n\sum _{k=\mathrm{0}}^n\left({\displaystyle \genfrac{}{}{0ex}{}{n}{k}}\right)(-\mathrm{1}{)}^k{\displaystyle \frac{\mathrm{1}}{\mathrm{2}(k+n)+\mathrm{1}}}{\cos }^{\mathrm{2}(k+n)+\mathrm{1}}\mathit{\varphi }.\end{array}$$

Hence,

$$\begin{array}{}\text{(A11)}& \begin{array}{rl}& \int u_{\mathit{\eta }}{\sin }^{\mathrm{2}n}\mathrm{2}\mathit{\varphi }\mathrm{2}a\mathrm{\Omega }\sin \mathit{\varphi }\mathrm{d}\mathit{\varphi }=-u_{\mathit{\eta }}\mathrm{2}a\mathrm{\Omega }{\mathrm{4}}^n\\ & \times \sum _{k=\mathrm{0}}^n\left(\genfrac{}{}{0ex}{}{n}{k}\right)(-\mathrm{1}{)}^k{\displaystyle \frac{\mathrm{1}}{\mathrm{2}(k+n)+\mathrm{1}}}{\cos }^{\mathrm{2}(k+n)+\mathrm{1}}\mathit{\varphi }.\end{array}\end{array}$$

Then, Eq. (A10b) is solved as

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\int {\sin }^{\mathrm{4}n}\left(\mathrm{2}\mathit{\varphi }\right)\tan \mathit{\varphi }\mathrm{d}\mathit{\varphi }\stackrel{+,\ast }{=}\int (\mathrm{2}\sin \mathit{\varphi }\cos \mathit{\varphi }{)}^{\mathrm{4}n}{\displaystyle \frac{\sin \mathit{\varphi }}{\cos \mathit{\varphi }}}\mathrm{d}\mathit{\varphi }\\ {\displaystyle }& {\displaystyle }={\mathrm{2}}^{\mathrm{4}n}\int {\sin }^{\mathrm{4}n+\mathrm{1}}\mathit{\varphi }{\cos }^{\mathrm{4}n-\mathrm{1}}\mathit{\varphi }\mathrm{d}\mathit{\varphi }\\ {\displaystyle }& {\displaystyle }={\mathrm{16}}^n\int {\sin }^{\mathrm{4}n+\mathrm{1}}\mathit{\varphi }{\cos }^{\mathrm{2}(\mathrm{2}n-\mathrm{1})}\mathit{\varphi })\cos \mathit{\varphi }\mathrm{d}\mathit{\varphi }\\ {\displaystyle }& {\displaystyle }\stackrel{\ast \ast }{=}{\mathrm{16}}^n\int {\sin }^{\mathrm{4}n+\mathrm{1}}\mathit{\varphi }(\mathrm{1}-{\sin }^{\mathrm{2}}\mathit{\varphi }{)}^{\mathrm{2}n-\mathrm{1}}\cos \mathit{\varphi }\mathrm{d}\mathit{\varphi },\end{array}$$

which is integrated by substituting u=sin ϕ and du=cos ϕ. That gives us

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathrm{16}}^n\int {\sin }^{\mathrm{4}n+\mathrm{1}}\mathit{\varphi }(\mathrm{1}-{\sin }^{\mathrm{2}}\mathit{\varphi }{)}^{\mathrm{2}n-\mathrm{1}}\cos \mathit{\varphi }\mathrm{d}\mathit{\varphi }\\ {\displaystyle }& {\displaystyle }={\mathrm{16}}^n\int u^{\mathrm{4}n+\mathrm{1}}(\mathrm{1}-u^{\mathrm{2}}{)}^{\mathrm{2}n-\mathrm{1}}\mathrm{d}u\\ {\displaystyle }& {\displaystyle }\stackrel{\ast \ast \ast }{=}{\mathrm{16}}^n\int \left(u^{\mathrm{4}n+\mathrm{1}}\sum _{k=\mathrm{0}}^{\mathrm{2}n-\mathrm{1}}\left({\displaystyle \genfrac{}{}{0ex}{}{\mathrm{2}n-\mathrm{1}}{k}}\right)(-\mathrm{1}{)}^k{u}^{\mathrm{2}k}\right)\mathrm{d}u\\ {\displaystyle }& {\displaystyle }={\mathrm{16}}^n\int \left(\sum _{k=\mathrm{0}}^{\mathrm{2}n-\mathrm{1}}\left({\displaystyle \genfrac{}{}{0ex}{}{\mathrm{2}n-\mathrm{1}}{k}}\right)(-\mathrm{1}{)}^k{u}^{\mathrm{2}k+\mathrm{4}n+\mathrm{1}}\right)\mathrm{d}u\\ {\displaystyle }& {\displaystyle }={\mathrm{16}}^n\sum _{k=\mathrm{0}}^{\mathrm{2}n-\mathrm{1}}\left({\displaystyle \genfrac{}{}{0ex}{}{\mathrm{2}n-\mathrm{1}}{k}}\right)(-\mathrm{1}{)}^k\int \left(u^{\mathrm{2}k+\mathrm{4}n+\mathrm{1}}\right)\mathrm{d}u\\ {\displaystyle }& {\displaystyle }={\mathrm{16}}^n\sum _{k=\mathrm{0}}^{\mathrm{2}n-\mathrm{1}}\left({\displaystyle \genfrac{}{}{0ex}{}{\mathrm{2}n-\mathrm{1}}{k}}\right)(-\mathrm{1}{)}^k{\displaystyle \frac{\mathrm{1}}{\mathrm{2}(k+\mathrm{2}n+\mathrm{1})}}{u}^{\mathrm{2}(k+\mathrm{2}n+\mathrm{1})}\\ {\displaystyle }& {\displaystyle }={\mathrm{16}}^n\sum _{k=\mathrm{0}}^{\mathrm{2}n-\mathrm{1}}\left({\displaystyle \genfrac{}{}{0ex}{}{\mathrm{2}n-\mathrm{1}}{k}}\right)(-\mathrm{1}{)}^k{\displaystyle \frac{\mathrm{1}}{\mathrm{2}(k+\mathrm{2}n+\mathrm{1})}}{\sin }^{\mathrm{2}(k+\mathrm{2}n+\mathrm{1})}\mathrm{2}\mathit{\varphi }.\end{array}$$

Hence,

$$\begin{array}{}\text{(A12)}& \begin{array}{rl}& -u_{\mathit{\eta }}^{\mathrm{2}}\int {\sin }^{\mathrm{4}n}\mathrm{2}\mathit{\varphi }\tan \mathit{\varphi }\mathrm{d}\mathit{\varphi }=-u_{\mathit{\eta }}^{\mathrm{2}}{\mathrm{16}}^n\sum _{k=\mathrm{0}}^{\mathrm{2}n-\mathrm{1}}\left(\genfrac{}{}{0ex}{}{\mathrm{2}n-\mathrm{1}}{k}\right)\\ & \times (-\mathrm{1}{)}^k{\displaystyle \frac{\mathrm{1}}{\mathrm{2}(k+\mathrm{2}n+\mathrm{1})}}{\sin }^{\mathrm{2}(k+\mathrm{2}n+\mathrm{1})}\mathrm{2}\mathit{\varphi }.\end{array}\end{array}$$

Combining the solutions of Eqs. (A10a) and (A10b) gives us the solution for Φ^′:

$$\begin{array}{}\text{(A13)}& \begin{array}{rl}{\mathrm{\Phi }}^{\prime }(\mathit{\lambda },\mathit{\varphi },\mathit{\eta })=& -u_{\mathit{\eta }}\mathrm{2}a\mathrm{\Omega }{\mathrm{4}}^n\sum _{k=\mathrm{0}}^n\left(\genfrac{}{}{0ex}{}{n}{k}\right)\\ & \times (-\mathrm{1}{)}^k{\displaystyle \frac{\mathrm{1}}{\mathrm{2}(k+n)+\mathrm{1}}}{\cos }^{\mathrm{2}(k+n)+\mathrm{1}}\mathit{\varphi }\\ & -u_{\mathit{\eta }}^{\mathrm{2}}{\mathrm{16}}^n\sum _{k=\mathrm{0}}^{\mathrm{2}n-\mathrm{1}}\left(\genfrac{}{}{0ex}{}{\mathrm{2}n-\mathrm{1}}{k}\right)\\ & \times (-\mathrm{1}{)}^k{\displaystyle \frac{\mathrm{1}}{\mathrm{2}(k+\mathrm{2}n+\mathrm{1})}}{\sin }^{\mathrm{2}(k+\mathrm{2}n+\mathrm{1})}\mathrm{2}\mathit{\varphi }\\ & +{\mathrm{\Phi }}_{\mathrm{0}}\left(\mathit{\eta }\right).\end{array}\end{array}$$

Since the deviations of Φ^′ vanish when averaging horizontally, we solve for Φ₀ by solving the equation

$$\begin{array}{}\text{(A14)}& {\displaystyle \frac{\mathrm{1}}{\mathrm{4}\mathit{\pi }}}\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\underset{-\mathit{\pi }/\mathrm{2}}{\overset{\mathit{\pi }/\mathrm{2}}{\int }}{\mathrm{\Phi }}^{\prime }(\mathit{\lambda },\mathit{\varphi },\mathit{\eta })\cos \mathit{\varphi }\mathrm{d}\mathit{\varphi }\mathrm{d}\mathit{\lambda }=\mathrm{0}\end{array}$$

by inserting the expression for ${\mathrm{\Phi }}^{\prime }(\mathit{\lambda },\mathit{\varphi },\mathit{\eta })$:

$$\begin{array}{}\text{(A15)}& \begin{array}{rl}& {\displaystyle \frac{\mathrm{1}}{\mathrm{4}\mathit{\pi }}}\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\underset{-\mathit{\pi }/\mathrm{2}}{\overset{\mathit{\pi }/\mathrm{2}}{\int }}[-u_{\mathit{\eta }}\mathrm{2}a\mathrm{\Omega }{\mathrm{4}}^n\sum _{k=\mathrm{0}}^n\left(\genfrac{}{}{0ex}{}{n}{k}\right)(-\mathrm{1}{)}^k{\displaystyle \frac{\mathrm{1}}{\mathrm{2}(k+n)+\mathrm{1}}}{\cos }^{\mathrm{2}(k+n)+\mathrm{1}}\mathit{\varphi }\\ & -u_{\mathit{\eta }}^{\mathrm{2}}{\mathrm{16}}^n\sum _{k=\mathrm{0}}^{\mathrm{2}n-\mathrm{1}}\left(\genfrac{}{}{0ex}{}{\mathrm{2}n-\mathrm{1}}{k}\right)(-\mathrm{1}{)}^k{\displaystyle \frac{\mathrm{1}}{\mathrm{2}(k+\mathrm{2}n+\mathrm{1})}}{\sin }^{\mathrm{2}(k+\mathrm{2}n+\mathrm{1})}\mathrm{2}\mathit{\varphi }\\ & +{\mathrm{\Phi }}_{\mathrm{0}}\left(\mathit{\eta }\right)]\cos \mathit{\varphi }\mathrm{d}\mathit{\varphi }\mathrm{d}\mathit{\lambda }=\mathrm{0}.\end{array}\end{array}$$

Equation (A15) is divided into three parts as follows.

$$\begin{array}{ll}{\displaystyle }& {\displaystyle \frac{\mathrm{1}}{\mathrm{4}\mathit{\pi }}}\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\underset{-\mathit{\pi }/\mathrm{2}}{\overset{\mathit{\pi }/\mathrm{2}}{\int }}(-u_{\mathit{\eta }}\mathrm{2}a\mathrm{\Omega }{\mathrm{4}}^n\sum _{k=\mathrm{0}}^n\left({\displaystyle \genfrac{}{}{0ex}{}{n}{k}}\right)\\ \text{(A16a)}& {\displaystyle }& {\displaystyle }\times (-\mathrm{1}{)}^k{\displaystyle \frac{\mathrm{1}}{\mathrm{2}(k+n)+\mathrm{1}}}{\cos }^{\mathrm{2}(k+n)+\mathrm{1}}\mathit{\varphi })\cos \mathit{\varphi }\mathrm{d}\mathit{\varphi }\mathrm{d}\mathit{\lambda }{\displaystyle }& {\displaystyle \frac{\mathrm{1}}{\mathrm{4}\mathit{\pi }}}\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\underset{-\mathit{\pi }/\mathrm{2}}{\overset{\mathit{\pi }/\mathrm{2}}{\int }}(-u_{\mathit{\eta }}^{\mathrm{2}}{\mathrm{16}}^n\sum _{k=\mathrm{0}}^{\mathrm{2}n-\mathrm{1}}\left({\displaystyle \genfrac{}{}{0ex}{}{\mathrm{2}n-\mathrm{1}}{k}}\right)\\ \text{(A16b)}& {\displaystyle }& {\displaystyle }\times (-\mathrm{1}{)}^k{\displaystyle \frac{\mathrm{1}}{\mathrm{2}(k+\mathrm{2}n+\mathrm{1})}}{\sin }^{\mathrm{2}(k+\mathrm{2}n+\mathrm{1})}\mathit{\varphi })\cos \mathit{\varphi }\mathrm{d}\mathit{\varphi }\mathrm{d}\mathit{\lambda }\text{(A16c)}& {\displaystyle }& {\displaystyle \frac{\mathrm{1}}{\mathrm{4}\mathit{\pi }}}\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\underset{-\mathit{\pi }/\mathrm{2}}{\overset{\mathit{\pi }/\mathrm{2}}{\int }}{\mathrm{\Phi }}_{\mathrm{0}}\left(\mathit{\eta }\right)\cos \mathit{\varphi }\mathrm{d}\mathit{\varphi }\mathrm{d}\mathit{\lambda }\end{array}$$

The first part, Eq. (A16a), is solved as follows

$$\begin{array}{ll}{\displaystyle }& {\displaystyle \frac{\mathrm{1}}{\mathrm{4}\mathit{\pi }}}\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\underset{-\mathit{\pi }/\mathrm{2}}{\overset{\mathit{\pi }/\mathrm{2}}{\int }}(-u_{\mathit{\eta }}\mathrm{2}a\mathrm{\Omega }{\mathrm{4}}^n\sum _{k=\mathrm{0}}^n\left({\displaystyle \genfrac{}{}{0ex}{}{n}{k}}\right)\\ {\displaystyle }& {\displaystyle }\times (-\mathrm{1}{)}^k{\displaystyle \frac{\mathrm{1}}{\mathrm{2}(k+n)+\mathrm{1}}}{\cos }^{\mathrm{2}(k+n)+\mathrm{1}}\mathit{\varphi })\cos \mathit{\varphi }\mathrm{d}\mathit{\varphi }\mathrm{d}\mathit{\lambda }\\ {\displaystyle }& {\displaystyle }=-u_{\mathit{\eta }}\mathrm{2}a\mathrm{\Omega }{\mathrm{4}}^n{\displaystyle \frac{\mathrm{1}}{\mathrm{4}\mathit{\pi }}}\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\underset{-\mathit{\pi }/\mathrm{2}}{\overset{\mathit{\pi }/\mathrm{2}}{\int }}(\sum _{k=\mathrm{0}}^n\left({\displaystyle \genfrac{}{}{0ex}{}{n}{k}}\right)\\ {\displaystyle }& {\displaystyle }\times (-\mathrm{1}{)}^k{\displaystyle \frac{\mathrm{1}}{\mathrm{2}(k+n)+\mathrm{1}}}{\cos }^{\mathrm{2}(k+n+\mathrm{1})}\mathit{\varphi })\mathrm{d}\mathit{\varphi }\mathrm{d}\mathit{\lambda }\\ {\displaystyle }& {\displaystyle }=-u_{\mathit{\eta }}\mathrm{2}a\mathrm{\Omega }{\mathrm{4}}^n{\displaystyle \frac{\mathrm{1}}{\mathrm{4}\mathit{\pi }}}\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}(\sum _{k=\mathrm{0}}^n\left({\displaystyle \genfrac{}{}{0ex}{}{n}{k}}\right)\\ {\displaystyle }& {\displaystyle }\times (-\mathrm{1}{)}^k{\displaystyle \frac{\mathrm{1}}{\mathrm{2}(k+n)+\mathrm{1}}}\underset{-\mathit{\pi }/\mathrm{2}}{\overset{\mathit{\pi }/\mathrm{2}}{\int }}{\cos }^{\mathrm{2}(k+n+\mathrm{1})}\mathit{\varphi }\mathrm{d}\mathit{\varphi })\mathrm{d}\mathit{\lambda }\end{array}$$

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }\stackrel{\mathit{†}}{=}-u_{\mathit{\eta }}\mathrm{2}a\mathrm{\Omega }{\mathrm{4}}^n{\displaystyle \frac{\mathrm{1}}{\mathrm{4}\mathit{\pi }}}\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}(\sum _{k=\mathrm{0}}^n\left({\displaystyle \genfrac{}{}{0ex}{}{n}{k}}\right)\\ {\displaystyle }& {\displaystyle }\times (-\mathrm{1}{)}^k{\displaystyle \frac{\mathrm{1}}{\mathrm{2}(k+n)+\mathrm{1}}}\sqrt{\mathit{\pi }}{\displaystyle \frac{\mathrm{\Gamma }(k+n+\mathrm{3}/\mathrm{2})}{\mathrm{\Gamma }(k+n+\mathrm{2})}})\mathrm{d}\mathit{\lambda }\\ {\displaystyle }& {\displaystyle }=-u_{\mathit{\eta }}\mathrm{2}a\mathrm{\Omega }{\mathrm{4}}^n{\displaystyle \frac{\mathrm{1}}{\mathrm{4}\mathit{\pi }}}\mathrm{2}\mathit{\pi }(\sum _{k=\mathrm{0}}^n\left({\displaystyle \genfrac{}{}{0ex}{}{n}{k}}\right)\\ {\displaystyle }& {\displaystyle }\times (-\mathrm{1}{)}^k{\displaystyle \frac{\mathrm{1}}{\mathrm{2}(k+n)+\mathrm{1}}}\sqrt{\mathit{\pi }}{\displaystyle \frac{\mathrm{\Gamma }(k+n+\mathrm{3}/\mathrm{2})}{\mathrm{\Gamma }(k+n+\mathrm{2})}})\\ {\displaystyle }& {\displaystyle }=-u_{\mathit{\eta }}a\mathrm{\Omega }{\mathrm{4}}^n(\sum _{k=\mathrm{0}}^n\left({\displaystyle \genfrac{}{}{0ex}{}{n}{k}}\right)\\ {\displaystyle }& {\displaystyle }\times (-\mathrm{1}{)}^k{\displaystyle \frac{\mathrm{1}}{\mathrm{2}(k+n)+\mathrm{1}}}\sqrt{\mathit{\pi }}{\displaystyle \frac{\mathrm{\Gamma }(k+n+\mathrm{3}/\mathrm{2})}{\mathrm{\Gamma }(k+n+\mathrm{2})}}),\end{array}$$

where Γ(z) is the gamma function. The second part, Eq. (A16b), is integrated by substituting u=sin ϕ and du=cos ϕ as follows

$$\begin{array}{rl}& {\displaystyle \frac{\mathrm{1}}{\mathrm{4}\mathit{\pi }}}\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\underset{-\mathit{\pi }/\mathrm{2}}{\overset{\mathit{\pi }/\mathrm{2}}{\int }}\left(u_{\mathit{\eta }}^{\mathrm{2}}{\mathrm{16}}^n\sum _{k=\mathrm{0}}^{\mathrm{2}n-\mathrm{1}}\left(\genfrac{}{}{0ex}{}{\mathrm{2}n-\mathrm{1}}{k}\right)\right(-\mathrm{1}{)}^k{\displaystyle \frac{\mathrm{1}}{\mathrm{2}(k+\mathrm{2}n+\mathrm{1})}}{\sin }^{\mathrm{2}(k+\mathrm{2}n+\mathrm{1})}\mathit{\varphi })\cos \mathit{\varphi }\mathrm{d}\mathit{\varphi }\mathrm{d}\mathit{\lambda }\\ & =-u_{\mathit{\eta }}^{\mathrm{2}}{\mathrm{16}}^n{\displaystyle \frac{\mathrm{1}}{\mathrm{4}\mathit{\pi }}}\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\left(\sum _{k=\mathrm{0}}^{\mathrm{2}n-\mathrm{1}}\left(\genfrac{}{}{0ex}{}{\mathrm{2}n-\mathrm{1}}{k}\right)\right(-\mathrm{1}{)}^k{\displaystyle \frac{\mathrm{1}}{\mathrm{2}(k+\mathrm{2}n+\mathrm{1})}}\underset{-\mathit{\pi }/\mathrm{2}}{\overset{\mathit{\pi }/\mathrm{2}}{\int }}{\sin }^{\mathrm{2}(k+\mathrm{2}n+\mathrm{1})}\mathit{\varphi }\cos \mathit{\varphi }\mathrm{d}\mathit{\varphi })\mathrm{d}\mathit{\lambda }\\ & \stackrel{\mathit{‡}}{=}-u_{\mathit{\eta }}^{\mathrm{2}}{\mathrm{16}}^n{\displaystyle \frac{\mathrm{1}}{\mathrm{4}\mathit{\pi }}}\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\left(\sum _{k=\mathrm{0}}^{\mathrm{2}n-\mathrm{1}}\left(\genfrac{}{}{0ex}{}{\mathrm{2}n-\mathrm{1}}{k}\right)\right(-\mathrm{1}{)}^k{\displaystyle \frac{\mathrm{1}}{\mathrm{2}(k+\mathrm{2}n+\mathrm{1})}}{\displaystyle \frac{\mathrm{2}}{\mathrm{2}(\mathrm{2}n+k+\mathrm{1})+\mathrm{1}}})\mathrm{d}\mathit{\lambda }\\ & =-u_{\mathit{\eta }}^{\mathrm{2}}{\mathrm{16}}^n{\displaystyle \frac{\mathrm{1}}{\mathrm{4}\mathit{\pi }}}\mathrm{2}\mathit{\pi }\sum _{k=\mathrm{0}}^{\mathrm{2}n-\mathrm{1}}\left(\genfrac{}{}{0ex}{}{\mathrm{2}n-\mathrm{1}}{k}\right)(-\mathrm{1}{)}^k{\displaystyle \frac{\mathrm{1}}{\mathrm{2}(k+\mathrm{2}n+\mathrm{1})}}{\displaystyle \frac{\mathrm{2}}{\mathrm{2}(\mathrm{2}n+k+\mathrm{1})+\mathrm{1}}}\\ & =-u_{\mathit{\eta }}^{\mathrm{2}}{\displaystyle \frac{{\mathrm{16}}^n}{\mathrm{2}}}\sum _{k=\mathrm{0}}^{\mathrm{2}n-\mathrm{1}}\left(\genfrac{}{}{0ex}{}{\mathrm{2}n-\mathrm{1}}{k}\right)(-\mathrm{1}{)}^k{\displaystyle \frac{\mathrm{1}}{\mathrm{2}(\mathrm{2}n+k+\mathrm{1})}}{\displaystyle \frac{\mathrm{2}}{\mathrm{2}(\mathrm{2}n+k+\mathrm{1})+\mathrm{1}}}.\end{array}$$

The third part, Eq. (A16c), is solved as

$${\displaystyle \frac{\mathrm{1}}{\mathrm{4}\mathit{\pi }}}\underset{\mathrm{0}}{\overset{\mathrm{2}\mathit{\pi }}{\int }}\underset{-\mathit{\pi }/\mathrm{2}}{\overset{\mathit{\pi }/\mathrm{2}}{\int }}{\mathrm{\Phi }}_{\mathrm{0}}\left(\mathit{\eta }\right)\cos \mathit{\varphi }\mathrm{d}\mathit{\varphi }\mathrm{d}\mathit{\lambda }={\mathrm{\Phi }}_{\mathrm{0}}\left(\mathit{\eta }\right),$$

which gives us

$$\begin{array}{ll}{\displaystyle }{\mathrm{\Phi }}_{\mathrm{0}}\left(\mathit{\eta }\right)=& {\displaystyle }-\left(\text{A16a }+\text{ A16b}\right)\\ {\displaystyle }=& {\displaystyle }-(-u_{\mathit{\eta }}a\mathrm{\Omega }{\mathrm{4}}^n\sum _{k=\mathrm{0}}^n\left({\displaystyle \genfrac{}{}{0ex}{}{n}{k}}\right)(-\mathrm{1}{)}^k\\ {\displaystyle }& {\displaystyle }\times {\displaystyle \frac{\mathrm{1}}{\mathrm{2}(k+n)+\mathrm{1}}}\sqrt{\mathit{\pi }}{\displaystyle \frac{\mathrm{\Gamma }(k+n+\mathrm{3}/\mathrm{2})}{\mathrm{\Gamma }(k+n+\mathrm{2})}}\\ {\displaystyle }& {\displaystyle }-u_{\mathit{\eta }}^{\mathrm{2}}{\displaystyle \frac{{\mathrm{16}}^n}{\mathrm{2}}}\sum _{k=\mathrm{0}}^{\mathrm{2}n-\mathrm{1}}\left({\displaystyle \genfrac{}{}{0ex}{}{\mathrm{2}n-\mathrm{1}}{k}}\right)(-\mathrm{1}{)}^k\\ {\displaystyle }& {\displaystyle }\times {\displaystyle \frac{\mathrm{1}}{\mathrm{2}(\mathrm{2}n+k+\mathrm{1})}}{\displaystyle \frac{\mathrm{2}}{\mathrm{2}(\mathrm{2}n+k+\mathrm{1})+\mathrm{1}}})\\ {\displaystyle }=& {\displaystyle }{u}_{\mathit{\eta }}a\mathrm{\Omega }{\mathrm{4}}^n\sum _{k=\mathrm{0}}^n\left({\displaystyle \genfrac{}{}{0ex}{}{n}{k}}\right)(-\mathrm{1}{)}^k\\ {\displaystyle }& {\displaystyle }\times {\displaystyle \frac{\mathrm{1}}{\mathrm{2}(k+n)+\mathrm{1}}}\sqrt{\mathit{\pi }}{\displaystyle \frac{\mathrm{\Gamma }(k+n+\mathrm{3}/\mathrm{2})}{\mathrm{\Gamma }(k+n+\mathrm{2})}}\\ {\displaystyle }& {\displaystyle }+u_{\mathit{\eta }}^{\mathrm{2}}{\displaystyle \frac{{\mathrm{16}}^n}{\mathrm{2}}}\sum _{k=\mathrm{0}}^{\mathrm{2}n-\mathrm{1}}\left({\displaystyle \genfrac{}{}{0ex}{}{\mathrm{2}n-\mathrm{1}}{k}}\right)(-\mathrm{1}{)}^k\\ \text{(A17)}& {\displaystyle }& {\displaystyle }\times {\displaystyle \frac{\mathrm{1}}{\mathrm{2}(\mathrm{2}n+k+\mathrm{1})}}{\displaystyle \frac{\mathrm{2}}{\mathrm{2}(\mathrm{2}n+k+\mathrm{1})+\mathrm{1}}}.\end{array}$$

Combining Eqs. (A13) and (A17) leads us to the final expression for the deviation of the geopotential field:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{\mathrm{\Phi }}^{\prime }(\mathit{\lambda },\mathit{\varphi },\mathit{\eta })=\\ {\displaystyle }& {\displaystyle }{u}_{\mathit{\eta }}a\mathrm{\Omega }{\mathrm{4}}^n(\underset{:=F_{\mathrm{3}}}{\underbrace{\sum _{k=\mathrm{0}}^n\left({\displaystyle \genfrac{}{}{0ex}{}{n}{k}}\right)(-\mathrm{1}{)}^k{\displaystyle \frac{\mathrm{1}}{\mathrm{2}(k+n)+\mathrm{1}}}\sqrt{\mathit{\pi }}{\displaystyle \frac{\mathrm{\Gamma }(k+n+\mathrm{3}/\mathrm{2})}{\mathrm{\Gamma }(k+n+\mathrm{2})}}}}\\ {\displaystyle }& {\displaystyle }-\mathrm{2}\underset{:=F_{\mathrm{1}}}{\underbrace{\sum _{k=\mathrm{0}}^n\left({\displaystyle \genfrac{}{}{0ex}{}{n}{k}}\right)(-\mathrm{1}{)}^k{\displaystyle \frac{\mathrm{1}}{\mathrm{2}(k+n)+\mathrm{1}}}{\cos }^{\mathrm{2}(k+n)+\mathrm{1}}\mathit{\varphi }}})\\ {\displaystyle }& {\displaystyle }+u_{\mathit{\eta }}^{\mathrm{2}}{\mathrm{16}}^n\\ {\displaystyle }& {\displaystyle }\times ({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\underset{:=F_{\mathrm{4}}}{\underbrace{\sum _{k=\mathrm{0}}^{\mathrm{2}n-\mathrm{1}}\left({\displaystyle \genfrac{}{}{0ex}{}{\mathrm{2}n-\mathrm{1}}{k}}\right)(-\mathrm{1}{)}^k{\displaystyle \frac{\mathrm{1}}{\mathrm{2}(\mathrm{2}n+k+\mathrm{1})}}{\displaystyle \frac{\mathrm{2}}{\mathrm{2}(\mathrm{2}n+k+\mathrm{1})+\mathrm{1}}}}}\\ {\displaystyle }& {\displaystyle }-\underset{:=F_{\mathrm{2}}}{\underbrace{\sum _{k=\mathrm{0}}^{\mathrm{2}n-\mathrm{1}}\left({\displaystyle \genfrac{}{}{0ex}{}{\mathrm{2}n-\mathrm{1}}{k}}\right)(-\mathrm{1}{)}^k{\displaystyle \frac{\mathrm{1}}{\mathrm{2}(k+\mathrm{2}n+\mathrm{1})}}{\sin }^{\mathrm{2}(k+\mathrm{2}n+\mathrm{1})}\mathrm{2}\mathit{\varphi }}}).\end{array}$$

It can be rewritten as

$$\begin{array}{}\text{(A18)}& \begin{array}{rl}{\mathrm{\Phi }}^{\prime }(\mathit{\lambda },\mathit{\varphi },\mathit{\eta })=& u_{\mathit{\eta }}a\mathrm{\Omega }{\mathrm{4}}^n(F_{\mathrm{3}}-\mathrm{2}{F}_{\mathrm{1}})\\ & +u_{\mathit{\eta }}^{\mathrm{2}}{\mathrm{16}}^n\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{F}_{\mathrm{4}}-F_{\mathrm{2}}\right),\end{array}\end{array}$$

with the following:

$$\begin{array}{ll}\text{(A19a)}& {\displaystyle }{F}_{\mathrm{1}}=& {\displaystyle }\sum _{k=\mathrm{0}}^n\left({\displaystyle \genfrac{}{}{0ex}{}{n}{k}}\right)(-\mathrm{1}{)}^k{\displaystyle \frac{\mathrm{1}}{\mathrm{2}(k+n)+\mathrm{1}}}{\cos }^{\mathrm{2}(k+n)+\mathrm{1}}\mathit{\varphi },{\displaystyle }{F}_{\mathrm{2}}=& {\displaystyle }\sum _{k=\mathrm{0}}^{\mathrm{2}n-\mathrm{1}}\left({\displaystyle \genfrac{}{}{0ex}{}{\mathrm{2}n-\mathrm{1}}{k}}\right)(-\mathrm{1}{)}^k{\displaystyle \frac{\mathrm{1}}{\mathrm{2}(k+\mathrm{2}n+\mathrm{1})}}\\ \text{(A19b)}& {\displaystyle }& {\displaystyle }\times {\sin }^{\mathrm{2}(k+\mathrm{2}n+\mathrm{1})}\mathrm{2}\mathit{\varphi },{\displaystyle }{F}_{\mathrm{3}}=& {\displaystyle }\sum _{k=\mathrm{0}}^n\left({\displaystyle \genfrac{}{}{0ex}{}{n}{k}}\right)(-\mathrm{1}{)}^k{\displaystyle \frac{\mathrm{1}}{\mathrm{2}(k+n)+\mathrm{1}}}\\ \text{(A19c)}& {\displaystyle }& {\displaystyle }\times \sqrt{\mathit{\pi }}{\displaystyle \frac{\mathrm{\Gamma }(k+n+\mathrm{3}/\mathrm{2})}{\mathrm{\Gamma }(k+n+\mathrm{2})}},{\displaystyle }{F}_{\mathrm{4}}=& {\displaystyle }\sum _{k=\mathrm{0}}^{\mathrm{2}n-\mathrm{1}}\left({\displaystyle \genfrac{}{}{0ex}{}{\mathrm{2}n-\mathrm{1}}{k}}\right)(-\mathrm{1}{)}^k{\displaystyle \frac{\mathrm{1}}{\mathrm{2}(\mathrm{2}n+k+\mathrm{1})}}\\ \text{(A19d)}& {\displaystyle }& {\displaystyle }\times {\displaystyle \frac{\mathrm{2}}{\mathrm{2}(\mathrm{2}n+k+\mathrm{1})+\mathrm{1}}},\text{(A19e)}& {\displaystyle }{u}_{\mathit{\eta }}=& {\displaystyle }{u}_{\mathrm{0}}\ln \mathit{\eta }\exp (-[\ln \mathit{\eta }/b{]}^{\mathrm{2}}).\end{array}$$

The geopotential field is now described as

$$\begin{array}{}\text{(A20)}& {\displaystyle }\mathrm{\Phi }(\mathit{\lambda },\mathit{\varphi },\mathit{\eta })& {\displaystyle }=\langle \mathrm{\Phi }\left(\mathit{\eta }\right)\rangle +{\mathrm{\Phi }}^{\prime }(\mathit{\lambda },\mathit{\varphi },\mathit{\eta })\text{(A21)}& {\displaystyle }& {\displaystyle }={\displaystyle \frac{T_{v,\mathrm{0}}g}{\mathit{\gamma }}}(\mathrm{1}-{\mathit{\eta }}^{\frac{R_{\mathrm{d}}\mathit{\gamma }}{g}})+{\mathrm{\Phi }}^{\prime }(\mathit{\lambda },\mathit{\varphi },\mathit{\eta }),\end{array}$$

where the expression for 〈Φ(η)〉 is the expression derived in Sect. A1 and the same as Eq. (7) (Ullrich et al., 2015).

A3 Virtual temperature field

We now have an expression for ${\mathrm{\Phi }}^{\prime }(\mathit{\lambda },\mathit{\varphi },\mathit{\eta })$, and we continue by deriving the virtual temperature anomaly field T^′ by starting from the hydrostatic balance:

$$\begin{array}{}\text{(A22)}& T_{\mathrm{v}}^{\prime }(\mathit{\lambda },\mathit{\varphi },\mathit{\eta })=-{\displaystyle \frac{\mathit{\eta }}{R_{\mathrm{d}}}}{\displaystyle \frac{\partial {\mathrm{\Phi }}^{\prime }(\mathit{\lambda },\mathit{\varphi },\mathit{\eta })}{\partial \mathit{\eta }}}.\end{array}$$

We start by inserting the expression for Φ^′ (Eq. A18) and u_η into Eq. (A22) and then take the derivative with respect to η:

$$\begin{array}{ll}{\displaystyle }& {\displaystyle }{T}_{\mathrm{v}}^{\prime }(\mathit{\lambda },\mathit{\varphi },\mathit{\eta })=-{\displaystyle \frac{\mathit{\eta }}{R_{\mathrm{d}}}}\\ {\displaystyle }& {\displaystyle }\times {\displaystyle \frac{\partial }{\partial \mathit{\eta }}}\left[u_{\mathit{\eta }}\underset{A_{\mathrm{1}}}{\underbrace{a\mathrm{\Omega }{\mathrm{4}}^n(F_{\mathrm{3}}-\mathrm{2}{F}_{\mathrm{1}})}}+u_{\mathit{\eta }}^{\mathrm{2}}\underset{A_{\mathrm{2}}}{\underbrace{{\mathrm{16}}^n({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{F}_{\mathrm{4}}-F_{\mathrm{2}}}})\right]\\ {\displaystyle }& {\displaystyle }=-{\displaystyle \frac{\mathit{\eta }}{R_{\mathrm{d}}}}\left[{\displaystyle \frac{\partial u_{\mathit{\eta }}}{\partial \mathit{\eta }}}{A}_{\mathrm{1}}+\mathrm{2}{u}_{\mathit{\eta }}{\displaystyle \frac{\partial u_{\mathit{\eta }}}{\partial \mathit{\eta }}}{A}_{\mathrm{2}}\right]\\ {\displaystyle }& {\displaystyle }=-{\displaystyle \frac{\mathit{\eta }}{R_{\mathrm{d}}}}\left[{\displaystyle \frac{\partial u_{\mathit{\eta }}}{\partial \mathit{\eta }}}(A_{\mathrm{1}}+\mathrm{2}{u}_{\mathit{\eta }}{A}_{\mathrm{2}})\right].\end{array}$$

As the terms $F_{\mathrm{1}},F_{\mathrm{2}},F_{\mathrm{3}}$, and F₄ do not depend on η, we only need to calculate

$$\begin{array}{}\text{(A23)}& \begin{array}{rl}{\displaystyle \frac{\partial u_{\mathit{\eta }}}{\partial \mathit{\eta }}}=& {\displaystyle \frac{\partial }{\partial \mathit{\eta }}}\left(u_{\mathrm{0}}\ln \mathit{\eta }\exp (-[\ln \mathit{\eta }/b{]}^{\mathrm{2}})\right)\\ =& u_{\mathrm{0}}\left[{\displaystyle \frac{\mathrm{1}}{\mathit{\eta }}}\exp \right(-[\ln \mathit{\eta }/b{]}^{\mathrm{2}})\\ & +\ln \mathit{\eta }\exp (-[\ln \mathit{\eta }/b{]}^{\mathrm{2}}\left){\displaystyle \frac{-\mathrm{2}}{b^{\mathrm{2}}}}\ln \mathit{\eta }{\displaystyle \frac{\mathrm{1}}{\mathit{\eta }}}\right]\\ =& u_{\mathrm{0}}{\displaystyle \frac{\mathrm{1}}{\mathit{\eta }}}\exp (-[\ln \mathit{\eta }/b{]}^{\mathrm{2}})\left(\mathrm{1}-{\displaystyle \frac{\mathrm{2}(\ln \mathit{\eta }{)}^{\mathrm{2}}}{b^{\mathrm{2}}}}\right).\end{array}\end{array}$$

Inserting Eq. (A23) into Eq. (A22) gives us

$$\begin{array}{}\text{(A24)}& \begin{array}{rl}& T_{\mathrm{v}}^{\prime }(\mathit{\lambda },\mathit{\varphi },\mathit{\eta })={\displaystyle \frac{u_{\mathrm{0}}}{R_{\mathrm{d}}}}\exp (-[\ln \mathit{\eta }/b{]}^{\mathrm{2}})\left({\displaystyle \frac{\mathrm{2}(\ln \mathit{\eta }{)}^{\mathrm{2}}}{b^{\mathrm{2}}}}-\mathrm{1}\right)\\ & \times \left[a\mathrm{\Omega }{\mathrm{4}}^n\left(F_{\mathrm{3}}-\mathrm{2}{F}_{\mathrm{1}}\right)+{\mathrm{16}}^n{u}_{\mathit{\eta }}\left(F_{\mathrm{4}}-\mathrm{2}{F}_{\mathrm{2}}\right)\right].\end{array}\end{array}$$

The virtual temperature field is now described by

$$\begin{array}{ll}\text{(A25)}& {\displaystyle }{T}_{\mathrm{v}}(\mathit{\lambda },\mathit{\varphi },\mathit{\eta })=& {\displaystyle }\langle T_{\mathrm{v}}\left(\mathit{\eta }\right)\rangle +T_{\mathrm{v}}^{\prime }(\mathit{\lambda },\mathit{\varphi },\mathit{\eta }){\displaystyle }=& {\displaystyle }{T}_{v,\mathrm{0}}{\mathit{\eta }}^{\frac{R_{\mathrm{d}}\mathit{\gamma }}{g}}+{\displaystyle \frac{u_{\mathrm{0}}}{R_{\mathrm{d}}}}\exp (-[\ln \mathit{\eta }/b{]}^{\mathrm{2}})\\ {\displaystyle }& {\displaystyle }\times \left({\displaystyle \frac{\mathrm{2}(\ln \mathit{\eta }{)}^{\mathrm{2}}}{b^{\mathrm{2}}}}-\mathrm{1}\right)\left[a\mathrm{\Omega }{\mathrm{4}}^n\right(F_{\mathrm{3}}-\mathrm{2}{F}_{\mathrm{1}})\\ \text{(A26)}& {\displaystyle }& {\displaystyle }+{\mathrm{16}}^n{u}_{\mathit{\eta }}\left(F_{\mathrm{4}}-\mathrm{2}{F}_{\mathrm{2}}\right)],\end{array}$$

where 〈T_v(η)〉 is the moist version of Eq. (10) (Ullrich et al., 2015).