The upper-atmosphere extension of the ICON general circulation model (version: ua-icon-1.0)

How the upper-atmosphere branch of the circulation contributes to and interacts with the circulation of the middle and lower atmosphere is a research area with many open questions. Inertia–gravity waves, for instance, have moved in the focus of research as they are suspected to be key features in driving and shaping the circulation. Numerical atmospheric models are an important pillar for this research. We use the ICOsahedral Non-hydrostatic (ICON) general circulation model, which is a joint development of the Max Planck Institute for Meteorology (MPI-M) and the German Weather Service (DWD), and provides, e.g., local mass conservation, a flexible grid nesting option, and a non-hydrostatic dynamical core formulated on an icosahedral–triangular grid. We extended ICON to the upper atmosphere and present here the two main components of this new configuration named UA-ICON: an extension of the dynamical core from shallow- to deep-atmosphere dynamics and the implementation of an upper-atmosphere physics package. A series of idealized test cases and climatological simulations is performed in order to evaluate the upper-atmosphere extension of ICON.

1 Introduction

In climate simulations and numerical weather prediction (NWP), there are ongoing efforts to raise the upper model lid, acknowledging possible influences of middle- and upper-atmosphere dynamics on tropospheric weather and climate (e.g., Thompson et al., 2002; Scaife et al., 2012; Charlton-Perez et al., 2013). The dynamics of the large-scale flow in the middle and upper atmosphere is determined, for instance, by the interaction with small-scale gravity waves. These waves, predominantly forced in the troposphere, can propagate vertically until they become unstable and break. As a result of this, a drag is exerted on the atmospheric background flow. This is an important route of momentum flux from the lower atmosphere to the middle and upper atmosphere, which shapes the meridional circulation in the latter regions (e.g., Fritts and Alexander, 2003; Kim et al., 2003). To have a model at hand that allows to study such processes on a wide range of spatial and temporal scales was one of our central motivations for the upper-atmosphere extension of the ICON model, which we present in the following.

The ICOsahedral Non-hydrostatic (ICON) general circulation model (Zängl et al., 2015; Dipankar et al., 2015; Heinze et al., 2017; Giorgetta et al., 2018; Crueger et al., 2018), a joint development of the Max Planck Institute for Meteorology (MPI-M) and the German Weather Service (DWD), is designed as a unified modeling system to allow simulations from the turbulent scales in large-eddy simulations (LESs) up to climatological timescales. Its triangular horizontal grid provides an almost uniform horizontal resolution, and the non-hydrostatic dynamical core makes it applicable on a wide range of spatial scales. The discretization of the governing equations makes use of a combination of finite-difference and finite-volume methods, e.g., to ensure mass conservation, which is an important property for climate simulations (e.g., Staniforth and Wood, 2008). In addition to a limited-area mode, ICON offers a grid-nesting option (one- or two-way nesting) for local refinement up to potentially very high horizontal resolutions. With three distinct packages of physics parameterizations, ICON meets the different needs of climate simulations, NWP and LES. To prepare ICON for simulations with model tops in the lower thermosphere, some extensions of the dynamical core and the physics parameterizations are necessary, as presented in this work.

Among the most important approximations which are applied to the dynamical core of ICON in its standard configuration are the shallow-atmosphere approximation and traditional approximation (e.g., Phillips, 1966; White and Bromley, 1995; White et al., 2005; Staniforth and Wood, 2008). These approximations are applied to the mapping of the budget equations on spherical coordinates relative to the center of the Earth. The shallow-atmosphere approximation is basically associated with the neglect of terms related to the spherical curvature of the atmosphere as well as variations of the gravitational field by assuming the field strength to be constant (see Thuburn and White, 2013, for a detailed examination of the metrical implications of this approximation). The traditional approximation refers to neglecting the contribution of the horizontal component of the Earth's angular velocity to the Coriolis acceleration. Following the usual terminology, we will call the system of equations with and without the two approximations the shallow-atmosphere equations and deep-atmosphere equations, respectively. Both approximations are generally applied together in order to satisfy the conservation of the energy, the axial component of angular momentum, and the potential vorticity (e.g., Phillips, 1966; White and Bromley, 1995; Staniforth and Wood, 2003). However, as Tort and Dubos (2014) have shown, it is possible to extend the shallow-atmosphere equations in such a way that the full Coriolis acceleration can be retained without violating the conservation principles.

The accuracy of the shallow-atmosphere and traditional approximations can be estimated by comparing the magnitude of the terms neglected in the shallow-atmosphere equations to the magnitude of the terms that are present in both systems. Such scale analysis has been used, for instance, by White and Bromley (1995), to show that for diabatically driven flows in the tropics and planetary-scale flows the neglected terms of the Coriolis acceleration might reach magnitudes up to about 10 % of the magnitude of key terms of the shallow-atmosphere momentum budget. On the other side, normal-mode analyses, done by Thuburn et al. (2002a) for the deep-atmosphere equations, and by Kasahara (2003) with focus on a Boussinesq model featuring the full Coriolis acceleration, show that the differences in the spatial structure and the frequencies of the energetically most significant modes between the shallow- and the deep-atmosphere equations are relatively small (with differences in the frequency magnitude being typically less than about 1 %; Thuburn et al., 2002b; Kasahara, 2003). Both, the scale analysis and the normal-mode analysis are important tools to figure out the differences between the deep- and shallow-atmosphere equations. However, the applicability of the results on long-term integrations might be limited. The systematic errors introduced by the approximations, albeit small in magnitude, could accumulate over time and lead to significantly different flow patterns of the model atmosphere. This might be especially important for the large-scale circulations of the middle and upper atmosphere. Furthermore, in view of the ever-increasing computational power, some approximations that were meaningful under the restrictions of past computer architectures might nowadays lose their justification. Therefore, we decided to expand the dynamical core of ICON by a deep-atmosphere option.

Examples for other models that use a deep-atmosphere formulation, or offer the option to do so, are the Met Office's Unified Model (UM) (Davies et al., 2005; Wood et al., 2014), the Non-hydrostatic Icosahedral Atmospheric Model (NICAM) (e.g., Tomita and Satoh, 2004), the Ocean Land Atmosphere Model (OLAM) (e.g., Walko and Avissar, 2008), the MCore model by Ullrich and Jablonowski (2012), and the Finite Volume Model (FVM) of the Integrated Forecasting System (IFS) developed at the European Centre for Medium-Range Weather Forecasts (ECMWF) (e.g., Smolarkiewicz et al., 2016). An overview of some of these models can be found in Ullrich et al. (2017).

Apart from the dynamics, the physics parameterizations are the second important model pillar that has to be extended for applications including part of the upper atmosphere. The ICON model offers basically three different physics packages: one which has been largely adopted from the ECHAM model intended for climate simulations (e.g., Stevens et al., 2013; Giorgetta et al., 2018; Crueger et al., 2018), another one which is used for NWP (some aspects of which can be found in Zängl et al., 2015), and a third one for LES (Dipankar et al., 2015). The upper-atmosphere-specific physics parameterizations have been integrated into the ECHAM and NWP packages, but we will focus on the extended ECHAM physics package in the remainder of this work. To avoid confusion in the following, a side note on our terminology may be in order: if we discuss the physics parameterizations, we typically use the attribute “upper atmosphere” to denote that their effects become significant in and above, say, the upper-mesosphere–lower-thermosphere region. In contrast, the attribute “deep atmosphere” is used in the context of the modifications of the dynamical core, since they apply to the entire air column, and no more or less well-defined vertical significance threshold can be made out for them. If we address both extensions as a whole, we use again the attribute “upper atmosphere”.

While the model lid is raised and the deep-atmospheric dynamics is applied, the model cannot produce physically reasonable results without the extended physics parameterizations. In fact, our experience shows that when the model lid is raised to above ∼85 km, the model without the extended physics would suffer from strong numerical instability, which could only be suppressed by using extremely large and unphysical numerical damping. The reason is related to the characteristics of the upper atmosphere, most importantly the rarefied air and the broader spectrum of incoming solar irradiance, which give rise to some physical phenomena that are negligible in the lower atmosphere and thus not parameterized in the current model but become crucial in maintaining the upper-atmospheric dynamics and thermodynamics.

One of such physical phenomena is molecular diffusion of momentum and heat, which is of negligible magnitude in the lower atmosphere compared to turbulent diffusion. In the upper atmosphere, as turbulence dies away and the air molecules are capable of traveling a long distance, molecular diffusion becomes dominant. Besides molecules, the upper atmosphere is also abundant with atoms and radicals produced as a result of photolysis. Chemical heating is the release of heat by recombination reactions between atoms or radicals, and is of particular importance in the upper atmosphere where photolysis products can travel large distances before recombining. Moreover, the higher altitude also means that the upper atmosphere would receive and absorb more solar irradiance in higher frequencies than at lower levels and on the surface. This brings the need of parameterizing the ultraviolet radiation including the Schumann–Runge bands and continuum, and the extreme ultraviolet bands. The solar radiation also acts to ionize the atmosphere, establishing the ionosphere from about 60 km. The electrically charged ions in the ionosphere are then aligned with the magnetic field of the Earth, thus creating a drag and a heating source to the neutral mass flow. Further, the usual assumption of local thermodynamical equilibrium (LTE) does not hold in the upper atmosphere as a consequence of the low collision frequency between particles; therefore, some modification must be made to the longwave radiation parameterization. The heating and the eddy diffusion of momentum and heat generated by the breaking of gravity waves have to be taken into account near the mesopause, too.

In recent decades, several general circulation models (GCMs) have been extended to the middle and upper atmosphere. An early development in this direction is the spectral model described by Miyahara et al. (1993). With an upper boundary at an altitude of about 165 km, it extends into the lower thermosphere, and has been used to study thermal tides, for instance. The influence of gravity wave parameterizations on the circulation in the middle atmosphere has been studied with the MA/ECHAM4 model, an extension of the hydrostatic spectral ECHAM4 up to the mesopause region (Manzini et al., 1997; Manzini and McFarlane, 1998). Later versions of MA/ECHAM provided the basis for the Hamburg Model of the Neutral and Ionized Atmosphere (HAMMONIA), which extends into the thermosphere and includes a package of upper-atmosphere-specific physics parameterizations, together with an interactive chemistry module (Schmidt et al., 2006). The upper-atmosphere physics package implemented into UA-ICON has been largely adopted from HAMMONIA. A good overview of the physical processes that are typically parameterized can also be found in the description of an extension of the Canadian Middle Atmosphere Model (CMAM) from an upper boundary at an altitude of about 95 to 210 km by Fomichev et al. (2002). In addition, the authors present an examination of how much these processes contribute to the energy and momentum budgets. Another spectral hydrostatic GCM that includes the upper-mesosphere–lower-thermosphere region is the Kühlungsborn Mechanistic general Circulation Model (KMCM, e.g., Becker, 2009). It has been used recently to study secondary gravity wave generation in the mesosphere as a result of the breaking of gravity waves originating from the troposphere (Becker and Vadas, 2018). In contrast to the aforementioned spectral GCMs, WACCM (Whole Atmosphere Community Climate Model; see, e.g., Richter et al., 2008 for its version 3), which is based on the Community Atmosphere Model and Community Climate Model (CAM and CCM) of the National Center for Atmospheric Research (NCAR), offers a hydrostatic finite-volume dynamical core. The different versions of WACCM generally extend to about 150 km and share many parameterizations of its upper-atmosphere physics package with CMAM and HAMMONIA. The offspring WACCM-X (Liu et al., 2010) extends the simulations even deeper into the thermosphere up to about 500 km. Recent developments, such as self-consistent electrodynamics, the transport of O⁺, and a modification of the dynamical core to account for the variation of the specific heats and mean molecular weight in the heterosphere (above ∼100 km, say) led to improved simulations of space weather and space climate (WACCM-X 2.0, Liu et al., 2018). As a final example, we mention the middle atmosphere GCM developed by Watanabe et al. (2008) and later complemented by physics parameterizations that allow for simulations covering the lower thermosphere up to about 150 km (the Japanese Atmospheric general circulation model for Upper Atmosphere Research (JAGUAR); Watanabe and Miyahara, 2009). It has been employed, for instance, to study the interaction of resolved gravity waves with thermal tides by global simulations of a horizontal, triangularly truncated spectral resolution of T213 and a vertical layer spacing of 500 m throughout the middle atmosphere (Watanabe and Miyahara, 2009). In a subsequent study on the dependence of the resolved gravity-wave-borne vertical flux of zonal momentum on the vertical resolution, the authors were able to conduct simulations with a layer spacing down to 200 m (Watanabe et al., 2015).

The outline of the paper is as follows: in Sect. 2, we describe the modifications and extensions to the dynamical core and to the physics parameterizations, followed by a presentation of results from idealized test cases and climate simulations for the evaluation of the new implementations in Sect. 3. We close with a conclusion in Sect. 4.

2 Model extension to the upper atmosphere

2.1 Deep-atmosphere dynamics

The dynamical core of the standard configuration of ICON makes use of the shallow-atmosphere approximation, which mainly consists in simplifying the governing equations measured in a spherical coordinate system in the following way: the radial distance of an air parcel to the center of the Earth r is approximated by the radius of the Earth a, and metrical terms which result from the unit vectors of the coordinate system to be functions of position are neglected. In addition, the traditional approximation is applied, by which the acceleration due to the horizontal component of the angular velocity of the Earth is neglected (Phillips, 1966; Staniforth and Wood, 2003). For atmospheric models having a model top below, say, the mesopause region at an altitude of about 70 to 100 km, the shallow-atmosphere approximation is likely a very good approximation (e.g., Ullrich et al., 2014), albeit some adverse impacts might exist, for instance, in the tropics where the cosine of latitude is of order 1, questioning the neglect of the non-traditional part of the Coriolis acceleration to some extent (White and Bromley, 1995; White et al., 2005). If the model top is raised into the lower thermosphere, the systematic errors introduced by the shallow-atmosphere approximation might start to outweigh its benefits, especially on a climatological timescale. So we think the extension of the dynamical core from shallow- to deep-atmosphere dynamics is an important component for UA-ICON.

2.1.1 Model equations

We restrict our considerations to the dry atmosphere hereafter in order to focus on the particular aspects of the deep-atmosphere dynamics. The budget equations for the momentum, mass, and heat are rearranged to the following set of prognostic equations (White et al., 2005; Zängl et al., 2015):

$$\begin{array}{}\text{(1)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\partial \mathit{v}}{\partial t}}& +\mathit{v}\cdot \mathbf{\nabla }\mathit{v}+\mathrm{2}\mathbf{\Omega }\times \mathit{v}+\mathbf{\Omega }\times \left(\mathbf{\Omega }\times \mathit{r}\right)=\\ & -c_p\mathit{\theta }\mathbf{\nabla }\mathit{\pi }-\mathbf{\nabla }{\mathit{\varphi }}_{\mathrm{g}}+{\mathit{\rho }}^{-\mathrm{1}}\mathit{F},\end{array}\text{(2)}& {\displaystyle \frac{\partial \mathit{\rho }}{\partial t}}+\mathbf{\nabla }\cdot \left(\mathit{v}\mathit{\rho }\right)=\mathrm{0},\text{(3)}& {\displaystyle \frac{\partial \mathit{\pi }}{\partial t}}+{\displaystyle \frac{R}{c_v}}{\displaystyle \frac{\mathit{\pi }}{\mathit{\rho }\mathit{\theta }}}\mathbf{\nabla }\cdot \left(\mathit{v}\mathit{\rho }\mathit{\theta }\right)={\displaystyle \frac{\mathrm{R}}{c_v{c}_p\mathit{\rho }\mathit{\theta }}}Q,\end{array}$$

and the equation of state reads

$$\begin{array}{}\text{(4)}& \mathit{\rho }\mathit{\theta }={\displaystyle \frac{p_{\mathrm{00}}{\mathit{\pi }}^{c_v/\mathrm{R}}}{\mathrm{R}}},\end{array}$$

where $\mathit{v}={\mathit{v}}_{\mathrm{h}}+w{\mathit{e}}_{\mathrm{r}}$ is the wind vector split in horizontal and radial components, Ω is the angular velocity vector of the Earth, and r denotes the position relative to the center of the Earth. In addition, $\mathit{\pi }=(p/p_{\mathrm{00}}{)}^{\mathrm{R}/c_p}$ denotes the Exner pressure relative to the reference pressure p₀₀=1000 hPa, $\mathit{\theta }=T/\mathit{\pi }$ is the potential temperature, ϕ_g denotes the gravitational potential, and c_p=1004.64 $\mathrm{J}{\mathrm{K}}^{-\mathrm{1}}{\mathrm{kg}}^{-\mathrm{1}}$, c_v=717.6 $\mathrm{J}{\mathrm{K}}^{-\mathrm{1}}{\mathrm{kg}}^{-\mathrm{1}}$, and $\mathrm{R}=c_p-c_v$ are the specific heat capacities at constant pressure and volume, and the gas constant for dry air, respectively. Finally, F denotes velocity tendencies due to dissipative processes or other parameterized processes, and Q denotes Exner pressure tendencies due to diabatic processes.

We have added the centrifugal acceleration $\mathbf{\Omega }\times (\mathbf{\Omega }\times \mathit{r})$ in Eq. (1), since it is of importance for a test case presented in Sect. 3.1.1. In the standard configuration of ICON, where this term is absent, we have apparent gravity $-\mathbf{\nabla }{\mathit{\varphi }}_{\mathrm{g}+\mathrm{c}}$ (true gravity plus centrifugal acceleration) on the right-hand side of Eq. (1), to which the spherical–geopotential approximation is applied: in the case of Eq. (1), the vertical unit vector of the coordinate system, in which the atmospheric flow is measured, is defined as the normal vector of the spherical equipotential surfaces ${\mathit{e}}_{\mathrm{r}}=\mathbf{\nabla }{\mathit{\varphi }}_{\mathrm{g}}/\left|\mathbf{\nabla }{\mathit{\varphi }}_{\mathrm{g}}\right|$. In the case of apparent gravity, we continue to write ${\mathit{e}}_{\mathrm{r}}=\mathbf{\nabla }{\mathit{\varphi }}_{\mathrm{g}+\mathrm{c}}/\left|\mathbf{\nabla }{\mathit{\varphi }}_{\mathrm{g}+\mathrm{c}}\right|$, but the equipotential surfaces are still approximated by spherical surfaces, neglecting the oblateness caused by the centrifugal acceleration (Gill, 1982; White et al., 2005; Vallis, 2006; Staniforth and Wood, 2008).

The advection term in Eq. (1) is expressed in the so-called 2-D vector-invariant formulation (e.g., Phillips, 1966; Sadourny, 1972; Vallis, 2006):

$$\begin{array}{}\text{(5)}& {\displaystyle }{\displaystyle }\mathit{v}\cdot \mathbf{\nabla }\mathit{v}={\mathit{\omega }}_{\mathrm{h}}\times {\mathit{v}}_{\mathrm{h}}+{\mathbf{\nabla }}_{\mathrm{h}}{\displaystyle \frac{{\left|{\mathit{v}}_{\mathrm{h}}\right|}^{\mathrm{2}}}{\mathrm{2}}}+w{\displaystyle \frac{\partial {\mathit{v}}_{\mathrm{h}}}{\partial r}}+\mathit{v}\cdot \mathbf{\nabla }\left(w{\mathit{e}}_{\mathrm{r}}\right),\end{array}$$

where ${\mathit{\omega }}_{\mathrm{h}}={\mathbf{\nabla }}_{\mathrm{h}}\times {\mathit{v}}_{\mathrm{h}}=\mathit{\zeta }{\mathit{e}}_{\mathrm{r}}$ is the vertical component of the relative vorticity. During the development of ICON, it turned out that it is advantageous to separate the Exner pressure into a hydrostatically balanced part and a deviation from it $\mathit{\pi }={\mathit{\pi }}_{\mathrm{0}}\left(r\right)+{\mathit{\pi }}^{\prime }$, where $-c_p{\mathit{\theta }}_{\mathrm{0}}\mathbf{\nabla }{\mathit{\pi }}_{\mathrm{0}}-\mathbf{\nabla }{\mathit{\varphi }}_{\mathrm{g}}=\mathrm{0}$. Thus, the right-hand side of Eq. (1) reads

$$\begin{array}{}\text{(6)}& {\displaystyle }{\displaystyle }-c_p\mathit{\theta }\mathbf{\nabla }\mathit{\pi }-\mathbf{\nabla }{\mathit{\varphi }}_{\mathrm{g}}=-c_p\mathit{\theta }\mathbf{\nabla }{\mathit{\pi }}^{\prime }-c_p{\mathit{\theta }}^{\prime }{\displaystyle \frac{\mathrm{d}{\mathit{\pi }}_{\mathrm{0}}}{\mathrm{d}r}},\end{array}$$

with ${\mathit{\theta }}^{\prime }=\mathit{\theta }-{\mathit{\theta }}_{\mathrm{0}}\left(r\right)$. The background temperature profile of the hydrostatically balanced part is defined by (Zängl, 2012)

$$\begin{array}{}\text{(7)}& {\displaystyle }{\displaystyle }{T}_{\mathrm{0}}\left(z_{\mathrm{g}}\right)=T_{\mathrm{str}}+\left(T_{\mathrm{sl}}-T_{\mathrm{str}}\right)\exp \left(-{\displaystyle \frac{z_{\mathrm{g}}}{H_{\mathrm{scal}}}}\right),\end{array}$$

where $z_{\mathrm{g}}=z/(\mathrm{1}+z/a)$ denotes the geopotential height, and T_sl=288.15 K as well as T_str=213.15 K are characteristic temperature values at sea level (z_g=0) and in the stratosphere, respectively. H_scal=10 000 m is a characteristic (geopotential) temperature scale height. The geopotential height follows as “natural” height measure from the gravitational acceleration (Gill, 1982; Wood and Staniforth, 2003; Wood et al., 2014; Ullrich et al., 2017):

$$\begin{array}{}\text{(8)}& \begin{array}{rl}\mathbf{\nabla }{\mathit{\varphi }}_{\mathrm{g}}\cdot & \mathrm{d}\mathit{r}=g{\left({\displaystyle \frac{a}{r}}\right)}^{\mathrm{2}}\mathrm{d}r=g{\left({\displaystyle \frac{\mathrm{1}}{\mathrm{1}+\frac{z}{a}}}\right)}^{\mathrm{2}}\mathrm{d}z\\ & =g\mathrm{d}\left({\displaystyle \frac{z}{\mathrm{1}+\frac{z}{a}}}\right)=g\mathrm{d}{z}_{\mathrm{g}},\end{array}\end{array}$$

with the geopotential ${\mathit{\varphi }}_{\mathrm{g}}=-g{a}^{\mathrm{2}}/r+{\mathit{\varphi }}_{\mathrm{g},\mathrm{0}}$, where ϕ_g,0 denotes some constant, and g=9.80665 m s⁻² is the mean gravitational acceleration at sea level (this is actually the magnitude of apparent gravity, but since the contribution of the centrifugal acceleration is relatively small, we use this value also for the true gravity). The values of z_g range from ${lim}_{z\to -a}{z}_{\mathrm{g}}=-\mathrm{\infty }$ to ${lim}_{z\to \mathrm{\infty }}{z}_{\mathrm{g}}=a$. The profile (7) allows to integrate the hydrostatic equilibrium analytically. Note that if the centrifugal acceleration is treated in its explicit form $\mathbf{\Omega }\times (\mathbf{\Omega }\times \mathit{r})$, it is not taken into account in the hydrostatic equilibrium outlined above.

In ICON, governing Eqs. (1) to (3) are measured in a spherical coordinate system with local unit vectors e_𝔱(r), e_𝔫(r), and e_r(r) (forming a right-handed system in this order), where e_𝔱 and e_𝔫 are the horizontal unit vectors tangential and normal to an edge separating two adjacent triangular cells, and e_r is defined in the center of the bottom and top surfaces of a cell (see Fig. 1)¹. In order to find expressions for the projection of the momentum equation (Eq. 1) onto the unit vectors of the local edge coordinate system and e_r, we can make use of the projections in the standard geographical coordinate system, with its horizontal unit vectors in meridional (south–north) and zonal (west–east) directions, e_φ and e_λ (e.g., Zdunkowski and Bott, 2003). The Coriolis and centrifugal accelerations are excluded from the following consideration, since their projections onto the unit vectors of the local coordinate systems of the triangular grid can be determined directly. Each edge is part of a great circle of the Earth, which can be regarded as an imaginary Equator. Rotating the geographical coordinate system accordingly and in such a way that e_𝔱 is parallel to ${\mathit{e}}_{{\mathit{\lambda }}^{\prime }}$ at the considered location (where the prime indicates this coordinate transformation), we would find e_𝔫 being parallel to ${\mathit{e}}_{{\mathit{\phi }}^{\prime }}$ (the radial unit vectors of both systems are parallel anyway). The governing equations for the velocity components $\mathit{v}=u{\mathit{e}}_{\mathit{\lambda }}+v{\mathit{e}}_{\mathit{\phi }}+w{\mathit{e}}_{\mathrm{r}}$ (or here, $\mathit{v}=u^{\prime }{\mathit{e}}_{{\mathit{\lambda }}^{\prime }}+v^{\prime }{\mathit{e}}_{{\mathit{\phi }}^{\prime }}+w{\mathit{e}}_{\mathrm{r}}$) can be found in many textbooks (e.g., Gill, 1982; Zdunkowski and Bott, 2003; Holton, 2004; Vallis, 2006). Evaluating them at the Equator (${\mathit{\phi }}^{(\prime )}=\mathrm{0}$) yields almost immediately the components of the governing equations for $\mathit{v}=v_{\mathfrak{t}}{\mathit{e}}_{\mathfrak{t}}+v_{\mathfrak{n}}{\mathit{e}}_{\mathfrak{n}}+w{\mathit{e}}_{\mathrm{r}}$ in the local edge coordinate systems if we identify v_𝔱 with u^′ and v_𝔫 with v^′. The equation for the vertical velocity component w in the local coordinate system defined in the center of the vertical cell interfaces can be found in a similar way. Together, they read

$$\begin{array}{}\text{(9)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\partial v_{\mathfrak{n}}}{\partial t}}& +\left\{{\displaystyle \frac{a}{r}}\right\}{\displaystyle \frac{\partial K_{\mathrm{h}}}{\partial \mathfrak{n}}}+\left(\mathit{\zeta }+f_{\mathrm{r}}\right)v_{\mathfrak{t}}+w\left({\displaystyle \frac{\partial v_{\mathfrak{n}}}{\partial r}}\underset{\mathrm{¯}}{+{\displaystyle \frac{v_{\mathfrak{n}}}{r}}-f_{\mathfrak{t}}}\right)\\ & \underset{\mathrm{¯}}{+{\mathrm{\Omega }}^{\mathrm{2}}r\sin \left(\mathit{\phi }\right)\cos \left(\mathit{\phi }\right){\mathit{e}}_{\mathit{\phi }}\cdot {\mathit{e}}_{\mathfrak{n}}}\\ & =-c_p\mathit{\theta }\left\{{\displaystyle \frac{a}{r}}\right\}{\displaystyle \frac{\partial \mathit{\pi }}{\partial \mathfrak{n}}}+{\displaystyle \frac{F_{\mathfrak{n}}}{\mathit{\rho }}},\end{array}\text{(10)}& {\displaystyle }\begin{array}{rl}{\displaystyle \frac{\partial w}{\partial t}}& +v_{\mathfrak{n}}\left(\left\{{\displaystyle \frac{a}{r}}\right\}{\displaystyle \frac{\partial w}{\partial \mathfrak{n}}}\underset{\mathrm{¯}}{-{\displaystyle \frac{v_{\mathfrak{n}}}{r}}+f_{\mathfrak{t}}}\right)\\ & +v_{\mathfrak{t}}\left(\left\{{\displaystyle \frac{a}{r}}\right\}{\displaystyle \frac{\partial w}{\partial \mathfrak{t}}}\underset{\mathrm{¯}}{-{\displaystyle \frac{v_{\mathfrak{t}}}{r}}-f_{\mathfrak{n}}}\right)+w{\displaystyle \frac{\partial w}{\partial r}}\underset{\mathrm{¯}}{-{\mathrm{\Omega }}^{\mathrm{2}}r{\cos }^{\mathrm{2}}\left(\mathit{\phi }\right)}\\ & =-c_p\left[\left({\mathit{\theta }}_{\mathrm{0}}+{\mathit{\theta }}^{\prime }\right){\displaystyle \frac{\partial {\mathit{\pi }}^{\prime }}{\partial r}}+{\mathit{\theta }}^{\prime }{\displaystyle \frac{\mathrm{d}{\mathit{\pi }}_{\mathrm{0}}}{\mathrm{d}r}}\right]+{\displaystyle \frac{F_r}{\mathit{\rho }}},\end{array}\end{array}$$

where $K_{\mathrm{h}}=(v_{\mathfrak{t}}^{\mathrm{2}}+v_{\mathfrak{n}}^{\mathrm{2}})/\mathrm{2}$ is the horizontal mass-specific kinetic energy, and f_r=2Ωsin (φ), $f_{\mathfrak{t},\mathfrak{n}}=\mathrm{2}\mathrm{\Omega }\cos \left(\mathit{\phi }\right){\mathit{e}}_{\mathit{\phi }}\cdot {\mathit{e}}_{\mathfrak{t},\mathfrak{n}}$ are the Coriolis parameters. The values of the projections ${\mathit{e}}_{\mathit{\phi }}\cdot {\mathit{e}}_{\mathfrak{t},\mathfrak{n}}$ are already provided by the standard shallow-atmosphere configuration of ICON, so they pose no additional problem. Here and in the following, we will formulate the deep-atmosphere equations as a modification of the shallow-atmosphere equations. This simplifies the comparison and corresponds actually to the way we have implemented the deep-atmosphere dynamics in ICON. For instance, we made use of the expansion of the gradient in the local coordinate systems $\mathbf{\nabla }={\mathit{e}}_{\mathfrak{t}}\mathit{\{}a/r\mathit{\}}\partial /\partial \mathfrak{t}+{\mathit{e}}_{\mathfrak{n}}\mathit{\{}a/r\mathit{\}}\partial /\partial \mathfrak{n}+{\mathit{e}}_{\mathrm{r}}\partial /\partial r$ if $\mathbf{\nabla }={\mathit{e}}_{\mathfrak{t}}\partial /\partial \mathfrak{t}+{\mathit{e}}_{\mathfrak{n}}\partial /\partial \mathfrak{n}+{\mathit{e}}_z\partial /\partial z$ is its expansion under shallow-atmosphere approximation (with $\partial /\partial \mathfrak{t}$ and $\partial /\partial \mathfrak{n}$ denoting the horizontal derivatives along great circle arcs under shallow-atmosphere approximation, i.e., $r=a+z\approx a$, and e_z=e_r, $\partial /\partial z=\partial /\partial r$). So the factors {⋅} become 1 under the shallow-atmosphere approximation. In addition, the underlined terms in Eqs. (9) and (10) are neglected in the shallow-atmosphere equations.

Figure 1 Schematic illustration of a global horizontal triangular grid as used by ICON (a) and of a grid cell (b). The dashed black lines illustrate a grid of type R2B0 (two-fold root division, zero bisections). C-grid staggering is used for the prognostic variables; i.e., all scalars like density ρ, potential temperature θ, and Exner pressure π are defined in the cell center (red dot), tangential and normal wind components v_𝔱 and v_𝔫 on the side faces (corresponding to the three primal edges of a triangle on the horizontal grid), and vertical wind component w on the bottom and top surfaces of the cell (assuming $v_{\mathfrak{t}},v_{\mathfrak{n}},w>\mathrm{0}$ for the cell drawing).

The velocity component v_𝔱 is reconstructed diagnostically from the adjacent values of v_𝔫 by a radial basis function interpolation and subsequent projection. More details on this approach can be found in Wan (2009, Appendix D).

The dissipative and diabatic processes summarized in the terms ρ⁻¹F and RQ∕(c_vc_pρθ) in Eqs. (1) and (3) are not modified for the deep atmosphere. Their modification for spherical coordinates would have introduced numerous additional metric terms into the governing equations (compare, e.g., Baldauf and Brdar, 2016), whose rigorous implementation into the dynamical core would cause considerable additional computational costs, which are not affordable for us for the time being.

2.1.2 Numerical implementation

A thorough description of the spatial and temporal discretization of the governing equations used in ICON can be found in Zängl et al. (2015), and for development stages of ICON in Bonaventura and Ringler (2005), Wan (2009) and Wan et al. (2013). Here, we will focus only on those elements of the spatial discretization that are affected by the deep-atmosphere modifications.

First, we would like to point to an important simplification, which we have done in our implementation of the deep-atmosphere dynamics into ICON. Since ICON is used for operational NWP at DWD, a key criterion for new developments to be integrated into the dynamical core is computational efficiency (of course, this criterion applies to climate models as well). For this reason, priority is given to efficiency over accuracy, where we assume this to be acceptable. Where possible the deep-atmosphere dynamics are realized by modification factors to the existing terms of the shallow-atmosphere dynamics. In addition, topography is not taken into account by the deep-atmosphere modification, which, for instance, means that the dynamics in a cell of the first grid layer above the Himalaya experience the same deep-atmosphere modification as the dynamics in a cell above the ocean surface. This simplification, albeit the most severe one, has many advantages. First, we avoid the complicated calculations of the geometric measures of a cell (e.g., its volume and surface areas) in spherical geometry if they are distorted by topography and lose the center of Earth as their center of curvature. Second, the abovementioned deep-atmosphere modification factors depend only on height, which saves a considerable amount of memory and computational cost. Since topography imprints on the grid layers only up to a certain height in ICON (typically up to 16 km), the errors introduced by this measure are assumed to be relatively small. For instance, the difference of a typical grid cell volume at mean sea level and at 16 km is of the order of a few per mille in the spherical geometry. Furthermore, the mass-conserving property of the dynamical core of ICON is not affected by our simplification.

We start with the modification of the gravitational acceleration (Wood and Staniforth, 2003):

$$\begin{array}{}\text{(11)}& {\displaystyle }{\displaystyle }g\to g{\left\{{\displaystyle \frac{a}{r}}\right\}}^{\mathrm{2}},\end{array}$$

with the modification factor parenthesized by braces. The formula (Eq. 11) enters Eq. (10) only implicitly via the hydrostatic background state. For its computation, Eq. (11) is used only once during model initialization. Next, we consider the gradient in edge-normal and tangential directions:

$$\begin{array}{}\text{(12)}& {\displaystyle }{\displaystyle }{\mathrm{grad}}_{\mathfrak{t},\mathfrak{n}}\left(\mathit{\chi }\right)={\displaystyle \frac{\mathrm{\Delta }\mathit{\chi }}{\mathrm{\Delta }\mathfrak{t},\mathfrak{n}}}\left\{{\displaystyle \frac{a}{r}}\right\},\end{array}$$

where χ stands for the mass-specific horizontal kinetic energy K_h and the Exner pressure π in Eq. (9), and the vertical velocity w in Eq. (10). $\mathrm{\Delta }\mathit{\chi }/\mathrm{\Delta }\mathfrak{t},\mathfrak{n}$ denotes the gradient as it is computed under shallow-atmosphere approximation. The factor a∕r derives from the modification of the length of horizontal distances $l\to l\mathit{\{}r/a\mathit{\}}$. Like all other modification factors, it is pre-computed during model initialization for the vertical positions, where the gradients have to be evaluated during runtime. Another quantity the computation of which has to be modified for the deep atmosphere is the vertical vorticity component ζ on the left-hand side of Eq. (9). ζ is computed at the vertices of the triangular cells and then interpolated to the edge centers, where v_𝔫 is defined. This allows to compute the vorticity from the v_𝔫 using the Stokes theorem:

$$\begin{array}{}\text{(13)}& {\displaystyle }{\displaystyle }\mathit{\zeta }=\left({\displaystyle \frac{\mathrm{1}}{A_{\mathrm{v}}}}\sum _{i=\mathrm{1}}^{n_{\mathrm{e}}}{v}_{\mathfrak{n},i}{l}_{\mathrm{d},i}{f}_{\mathrm{o},i}\right)\left\{{\displaystyle \frac{a}{r}}\right\},\end{array}$$

where A_v is the area of the hexagonal or pentagonal cells centered at the vertices, which form the dual cells to the primal triangular cells by connecting the mass points of the triangular cells around a vertex; see Fig. 2 in Wan et al. (2013). In addition, n_e is the number of the dual edges of length l_d around the vertex, crossing the primal edges of the triangular cells perpendicularly (n_e=6, or 5 for the 12 pentagon points of the grid, respectively). Finally, f_o is an orientation factor, which is 1 or −1 according to whether e_𝔫 is parallel or antiparallel to the cyclonic direction of the integration path. Again, the first factor on the right-hand side of Eq. (13) is the vorticity under shallow-atmosphere approximation; thus, the geometric measures A_v and l_d are those at r=a. The modification factor a∕r results from the quotient l_d∕A_v in spherical geometry and is identical to the modification factor for the horizontal gradient (Eq. 12). Equations. (2) and (3) are written in flux form, which guarantees mass conservation up to machine precision and the conservation of heat as far as advective fluxes are concerned. The deep-atmosphere modifications of the flux divergences result from the quotient of area and volume A∕V for the five surfaces of a triangular cell. The modification of a cell volume itself reads

$$\begin{array}{}\text{(14)}& {\displaystyle }{\displaystyle }V=A_{\mathrm{c}}(z_{\text{top}}-z_{\text{bot}})\left\{{\displaystyle \frac{r_{\text{bot}}^{\mathrm{2}}+r_{\text{bot}}{r}_{\text{top}}+r_{\text{top}}^{\mathrm{2}}}{\mathrm{3}{a}^{\mathrm{2}}}}\right\},\end{array}$$

where $V=A_{\mathrm{c}}(z_{\text{top}}-z_{\text{bot}})$ is the cell volume under shallow-atmosphere approximation, with the cell surface area A_c, and $r_{\text{bot}}=a+z_{\text{bot}}$, $r_{\text{top}}=a+z_{\text{top}}$ are the radii of the cell's bottom surface and top surface, respectively. This modification is required, for instance, if global integrals of mass and other quantities are computed for diagnostic purposes. Together with the area modification for the three side faces and the bottom and top surfaces,

$$\begin{array}{}\text{(15)}& {\displaystyle }& {\displaystyle }{A}_{\mathrm{s}}=l_{\mathrm{p}}\left(z_{\text{top}}-z_{\text{bot}}\right)\left\{{\displaystyle \frac{r_{\text{bot}}+r_{\text{top}}}{\mathrm{2}a}}\right\},\text{(16)}& {\displaystyle }& {\displaystyle }{A}_{\text{bot},\text{top}}=A_{\mathrm{c}}{\left\{{\displaystyle \frac{r_{\text{bot},\text{top}}}{a}}\right\}}^{\mathrm{2}},\end{array}$$

where l_p is the length of the primal edges, we find for the flux divergences

$$\begin{array}{}\text{(17)}& \begin{array}{rl}& \mathrm{div}\left(\mathit{\xi }\right)=\left({\displaystyle \frac{\mathrm{1}}{A_{\mathrm{c}}}}\sum _{i=\mathrm{1}}^{\mathrm{3}}{\mathit{\xi }}_i{l}_{\mathrm{p},i}{f}_{\mathrm{o},i}\right)\left\{{\displaystyle \frac{\mathrm{3}}{\mathrm{4}}}{\displaystyle \frac{\mathrm{2}a}{r_{\text{bot}}+r_{\text{top}}}}{\left[\mathrm{1}-{\displaystyle \frac{r_{\text{bot}}{r}_{\text{top}}}{{\left(r_{\text{bot}}+r_{\text{top}}\right)}^{\mathrm{2}}}}\right]}^{-\mathrm{1}}\right\}\\ & +{\displaystyle \frac{{\mathit{\xi }}_{\text{top}}\left\{\mathrm{3}{\left[\mathrm{1}+\frac{r_{\text{bot}}}{r_{\text{top}}}+{\left(\frac{r_{\text{bot}}}{r_{\text{top}}}\right)}^{\mathrm{2}}\right]}^{-\mathrm{1}}\right\}-{\mathit{\xi }}_{\text{bot}}\left\{\mathrm{3}{\left[\mathrm{1}+\frac{r_{\text{top}}}{r_{\text{bot}}}+{\left(\frac{r_{\text{top}}}{r_{\text{bot}}}\right)}^{\mathrm{2}}\right]}^{-\mathrm{1}}\right\}}{z_{\text{top}}-z_{\text{bot}}}}.\end{array}\end{array}$$

Here, ${\mathit{\xi }}_{i,\text{top},\text{bot}}$ denotes the scalar product of the flux density and e_𝔫, or e_r at the respective cell face. The first term on the right-hand side of Eq. (17) sums the fluxes across the side faces, where the orientation factor f_o is either 1 or −1 according to whether e_𝔫 is parallel or antiparallel to the normal vector of the side faces from the point of view of the cell. The second term sums the fluxes across the top and bottom surfaces.

The underlined terms in Eqs. (9) and (10) are new implementations in the dynamical core. Their discrete formulation closely follows the formulation of structurally similar terms of the discrete shallow-atmosphere equations in Zängl et al. (2015).

The discretization employed in ICON necessitates the spatial interpolation of variables between cell center, edge midpoint, vertex, and the center of the vertical cell interfaces (bottom and top surface; Zängl et al., 2015), which could be affected potentially when changing from a shallow- to a deep-atmosphere formulation. In the horizontal, linear, bilinear, and area-weighted interpolations are used primarily. These interpolations are assumed to take place on coordinate surfaces z=const., i.e., planes in the case of the shallow atmosphere (at least practically, but see Thuburn and White, 2013, for the geometric implications of the shallow-atmosphere approximation) and spherical shells in the case of the deep atmosphere. This means that no explicit deep-atmosphere modification is necessary, because it would be the same for the respective geometric measures (distances, areas) on a coordinate surface, which cancel each other in the interpolation. For the vertical interpolation between cell centers and the cell interfaces, a linear interpolation is used in the case of the shallow atmosphere. We interpret this as an interpolation along the coordinate lines $\mathit{\lambda },\mathit{\phi }=\text{const.}$ in the case of the deep atmosphere, so that we can avoid deep-atmosphere modifications in this case, too. If we interpret the vertical interpolation, e.g., as a volume-weighted average, a modification would be necessary. Apart from the aforementioned standard interpolations, ICON makes use of an upwind-biased reconstruction of density and potential temperature from cell centers to the centers of the (horizontal and vertical) interfaces. These reconstructions enter, for instance, the divergence of mass and heat flux (Zängl et al., 2015) and make use of a Taylor expansion up to the first order along backward trajectories, which requires modifying the corresponding gradient according to Eq. (12).

For the time integration, a two-time-level predictor–corrector scheme is used, which is not directly affected by deep-atmosphere modifications. However, certain terms on the right-hand side of the governing equations, which are associated with vertical sound wave propagation, are treated implicitly in the temporal discretization in order to allow reasonable time steps. This directly affects the discretized forms of Eqs. (3) and (10), which can be combined to yield a diagnostic, elliptic equation for the vertical wind w for each vertical grid cell column of the model atmosphere (Zängl et al., 2015). Thus, some of the coefficients of this equation have to be complemented by the corresponding deep-atmosphere modification factors.

2.1.3 Model initialization

One motivation to implement the deep-atmosphere dynamics in ICON is to increase the accuracy of simulations with a model top ≳100 km. However, initial data are usually only available up to altitudes of about 70 to 80 km. For instance, the model top is at 75 km in the operational NWP with ICON at DWD, and the IFS model of the ECMWF, whose operational analysis data can be used to initialize ICON, has its model top at 0.01 hPa (see https://www.ecmwf.int/en/forecasts/documentation-and-support, last access: October 2018). To our knowledge, no regular and reliable measurements with global coverage are available for the (lower) thermosphere, which could be used for some kind of data assimilation. Therefore, climatological tables appear as a possible second choice. However, the momentary state at a point in the atmosphere generally deviates more or less strongly from the state given by any climatology. The model thermosphere initialized with climatological data will consequently undergo an adjustment process during some spin-up phase. With this uncertainty in mind, we decided to begin with a very simple approach. The state of the lower and middle atmosphere is initialized by an IFS analysis. In order to complement this by a state for the upper atmosphere, a mean vertical temperature profile (neglecting horizontal and temporal variations) obtained from Bates (1959), Hedin (1983), and Fleming et al. (1988) is used. It provides the basis to obtain first the pressure from the hydrostatic balance, followed by a computation of the horizontal wind from the geostrophic balance. The technical details can be found in Appendix A. Later on, this simple approach could be improved on by using a climatology generated by UA-ICON itself.

2.2 Upper-atmosphere physics

A new physics package, which parameterizes processes specific to the upper atmosphere has been developed for UA-ICON. This package, referred to as the UA package, can be called in combination with either the NWP package or the ECHAM package. The processes taken into consideration in the UA package are summarized in Table 1.

Huang et al. (1998)Banks and Kockarts (1973)Gill (1982)Hong and Lindzen (1976)Hines (1997a, 1997b)Strobel (1978)Richards et al. (1994)Fomichev and Blanchet (1995, 1998)Ogibalov and Fomichev (2003)Kockarts (1980)

Table 1 Physical parameterizations implemented in the UA package of UA-ICON shown with references and heights, above which their computation starts.

The processes are categorized into three groups (kinetics, radiation, and chemical heating), as described below. Most of the parameterizations are adopted from the HAMMONIA model, a spectral model based on ECHAM5 (Roeckner et al., 2006), covering the atmosphere up to the thermosphere ($\mathrm{1.7}\times {\mathrm{10}}^{-\mathrm{7}}$ hPa, ∼250 km). A detailed description of the physics parameterizations used in HAMMONIA can be found in Schmidt et al. (2006); thus, here we keep the description brief, only noting important and differing treatments. An overall difference is that in UA-ICON all the parameterizations are implemented such that the computation only starts at a certain altitude above which the forcings are expected to become relevant (see Table 1). This increases computational efficiency significantly.

2.2.1 Kinetics

Above the mesopause, molecular diffusion, which is negligible at lower altitude, becomes significant. In fact, there the downward transport of heat by molecular diffusion appears as a strong cooling in balance with the strong solar heating. Hence, it is of primary importance to parameterize molecular processes in this region of the upper atmosphere. Molecular transport of heat, momentum, and tracers is parameterized in the UA package following Huang et al. (1998) and Banks and Kockarts (1973), as in HAMMONIA. The computation starts at 75 km. Besides direct transport of heat by molecular diffusion, the momentum transport also leads to energy deposition in the form of heat, known as frictional heating. In the UA package, this is parameterized following Gill (1982). The computation starts at 85 km.

In the mesosphere and lower thermosphere, unlike in lower layers, a larger number of air particles are ionized and thus aligned with the electromagnetic field of the Earth. This produces a force on the neutral mass flow, its tangential and normal components known as the ion drag and the Lorenz force, respectively. Joule heating is produced by the ion drag as well. In UA-ICON, as in HAMMONIA, this effect is parameterized following the simple Hong and Lindzen (1976) approach. The computation starts at 80 km.

A large portion of gravity wave momentum energy is deposited and transferred to turbulent energy near the mesopause, where turbulence is otherwise very weak. The turbulent mixing effects induced by gravity waves can be estimated using the Hines (1997a, 1997b) parameterization included in the ECHAM package. This option is switched off in standard ICON simulations. In UA-ICON simulations, however, we enable the calculation and pass the computed turbulent diffusion coefficient to the turbulent mixing subroutine to account for gravity-wave-induced turbulent mixing.

2.2.2 Radiation

In standard ICON, the PSrad radiation package (Pincus and Stevens, 2013) is employed, which is itself an extension to the Rapid Radiative Transfer Model for GCMs (RRTMG) (Mlawer et al., 1997; Iacono et al., 2008). The shortwave component of the radiation package covers wavelengths of the solar spectrum longer than 200 nm. This is a sufficient bandwidth only up to about the mesopause, beyond which radiative heating in the ultraviolet and extreme ultraviolet frequencies become dominant. In the UA package, starting from 50 km, ultraviolet solar forcing for the O₂ Schumann–Runge bands (SRBs; 175 to 205 nm) and continuum (SRC; 125 to 175 nm) is calculated based on the model of Strobel (1978). Efficiency factors multiplied by the SRBC (SRB plus SRC) heating rates account for the loss of internal energy due to airglow processes and are taken from Mlynczak and Solomon (1993). For the extreme ultraviolet (EUV; 5 to 105 nm) solar forcing, starting above 90 km, a model based on Richards et al. (1994) taken from HAMMONIA is currently used. Efficiency factors multiplied by the EUV heating rates are based on Roble (1995). Their values also account for the energy loss due to radiative cooling in the 5.3 µm NO band. Since this process is explicitly calculated in our model (see below), a factor of 1.33 is multiplied by these efficiency factors (see Richards et al., 1982). Some additional adjustments to the PSrad/RRTMG shortwave radiation are necessary in UA-ICON due to the introduction of chemical heating. More details on this are given in Sect. 2.2.3.

The longwave component of PSrad/RRTMG covers terrestrial wavelengths shorter than 1 mm. This bandwidth is still valid at thermospheric heights, yet a few important additions have to be made:

1. The usual assumption of local thermodynamical equilibrium (LTE) does not hold above the mesopause; thus, non-LTE effects must be taken into account. As in HAMMONIA, non-LTE infrared cooling by O₃ and CO₂ is calculated from the parameterization of Fomichev and Blanchet (1995) with the modifications of Fomichev and Blanchet (1998). The calculation starts at 65 km, and the calculated values are multiplied by a scaling factor α equaling 0 at 65 km and linearly growing to 1 at 75 km. Correspondingly, the longwave radiation computed by PSrad/RRTMG is scaled with the factor 1-α, effectively discarding it above 75 km.

2. As in HAMMONIA, a parameterization of CO₂ non-LTE absorption in the near-infrared following Ogibalov and Fomichev (2003) is employed. The computed values are ignored below 19.25 km and fully considered above 24.5 km.

3. NO cooling at 5.3 µm is calculated utilizing the parameterization from Kockarts (1980). The computation starts at 60 km.

A noteworthy difference between ICON and HAMMONIA is that ICON is a non-hydrostatic model on height levels, whereas HAMMONIA is hydrostatic on hybrid pressure levels. In HAMMONIA, for the use in the radiation computation, number densities of radiatively active tracers are calculated based on the mass of air in a given layer which is derived from pressure differences between the upper and lower surfaces of a layer. This approach is only valid under the assumption of hydrostatic balance, since otherwise pressure is not guaranteed to decrease strictly monotonically with increasing altitude. Therefore, in the UA package, the computation of number density for the radiation parameterization is utilizing the mass of air, which is a globally conserved quantity in ICON.

Moreover, HAMMONIA has the upper boundary of its top pressure level at 0 hPa, effectively covering the whole atmosphere. The height levels of ICON, on the other hand, cover a finite range and unavoidably ignore the atmospheric air mass above the model lid. The effect of this missing amount of air on radiative fluxes is ignored in UA-ICON.

2.2.3 Chemical heating

Chemical heating, i.e., the release of heat due to recombination reactions between atoms or radicals produced as a result of photolysis, becomes important in maintaining the upper-atmospheric thermodynamic balance, where the photolysis products can travel a long distance before recombining. HAMMONIA employs a condensed version of the MOZART3 chemistry model (Kinnison et al., 2007) to explicitly compute chemical heating online. In our case, however, given the finer target resolution (from 160 km down to a few tens of kilometers) and the central goal of studying gravity waves, using a coupled chemistry model is overly expensive. Therefore, we deploy a simpler strategy of prescribing monthly zonal-mean climatological chemical heating rates from a 35-year HAMMONIA simulation with constant present-day boundary conditions. As an example, the chemical heating rates for January are shown in Fig. 2. Technically, below 70 km, all heating is calculated in the PSrad/RRTMG shortwave radiation code and no chemical heating rates are prescribed, whereas above 80 km full chemical heating rates are used and the radiative heating provided by the PSrad/RRTMG scheme is reduced by 23 % in order not to count solar energy twice. This approximation accounts for the energy used to break the chemical bond of ozone. Additionally, height-dependent efficiency factors from Mlynczak and Solomon (1993) are applied at these altitudes. This approach was also used in HAMMONIA (Schmidt et al., 2006). Between 70 and 80 km, the two heating sources are linearly merged.²

Figure 2 Zonal-mean chemical heating rates (K d⁻¹) averaged for the month of January from HAMMONIA simulations. Such monthly and zonal means are prescribed in UA-ICON.

3 Model evaluation

3.1 Idealized test cases

To test the deep-atmosphere implementation in the dynamical core, we used two test cases. In the first test case, the propagation of a sound wave is considered, for which an analytical solution (of the linearized equations) is available. It is aimed at testing especially the accuracy of the spherical geometry in its imprint on the grid cells and the corresponding modification factors described in Sect. 2.1.2, and at testing the metric terms and the complete Coriolis acceleration in the components of the momentum Eqs. (9) and (10). The second test case is the Jablonowski–Williamson baroclinic instability test case (Jablonowski and Williamson, 2006) in its extension for deep-atmosphere dynamical cores by Ullrich et al. (2014). It reveals if the height dependence of gravity is properly implemented and maintains the hydrostatic background state of the test case atmosphere (especially if gravity enters the momentum equation only implicitly, as in Eq. 10). In addition, the performance of the entire dynamical core is tested when it comes to reproducing the development of the baroclinic wave. Both test cases make use of the small-Earth approach of Wedi and Smolarkiewicz (2009) to pronounce deep-atmosphere effects.

3.1.1 Sound wave test case

The particular motivation for this first test case is that an analytical solution is available to which the numerical solution can be compared. We have developed this test case with a method originally proposed by Läuter et al. (2005) for the shallow-water equations on the sphere. The method was developed further for the shallow- and deep-atmosphere equations, e.g., by Staniforth and White (2008); Baldauf et al. (2014). An atmosphere at rest in the absolute frame is considered. If a non-trivial analytical solution is known for this case, it can be transformed into a rotating frame (e.g., regarded as a rotating Earth slipping through the air without exchange of tangential momentum). Depending on the solution in the absolute frame being either stationary or time dependent, potentially all aspects of a dynamical core can be tested. However, a disadvantage of this method is that the centrifugal acceleration has to be taken into account explicitly in the dynamical core (see Eq. 1). Some aspects of this transformation method are shown in Appendix B1, and a thorough mathematical description can be found in the literature cited above.

Given a solution in the absolute frame, it appears to be advected with $\mathit{v}=-{\mathit{v}}_F=-\mathbf{\Omega }\times (\mathit{X}-\mathit{A})$ from the perspective of the rotating frame (with the center of Earth A and an arbitrary point X not coincident with A). In practice, this means that the solution has merely to be rotated by an angle −Ωt about an axis being parallel to Ω and crossing A. Therefore, we will direct our attention to the solution to be rotated in the following.

Baldauf et al. (2014) derived analytical gravity and sound wave solutions for the linearized deep-atmosphere equations. However, certain terms of the equations had to be omitted in order to allow the solution to be expanded in a system of orthonormal basis functions. Although these terms were shown to be only of second order, their effects are not controllable in a dynamical core, where the omission of most of the corresponding discretized terms is unfeasible without greater effort. Therefore, it would be desirable to find solutions for the linearized equations with the omitted terms restored, or alternatively with only those omissions retained which can be realized in a dynamical core without greater effort. The fact that the gravity $-g(a/r{)}^{\mathrm{2}}(\mathit{X}-\mathit{A})/r$ imprints a spherical symmetry on the equations with its distinct point A turns out to be a severe obstacle to that aim. So we decided to switch off gravity (which in model practice means to set the constant parameter g to a very small value >0, since divisions by g are used throughout the model, especially in the physics parameterizations). This greatly simplifies the problem but has the consequence that the test case makes no statement about the implementation of gravity. Under these circumstances, the atmosphere, enclosed between the spherical boundaries of the model bottom and top, is isotropic, with the constant pressure p₀ and temperature T₀. As long as the sound waves, propagating with the speed of sound $c_{\mathrm{s}}=\sqrt{(c_p/c_v)\mathrm{R}{T}_{\mathrm{0}}}$, do not interact with the boundaries, they are not “aware” of the spherical shape of the atmosphere as a whole. Therefore, the challenge for the model is to properly simulate the sound wave propagation on the anisotropic spherically curved grid. For this test case, we consider a spherically symmetric acoustic wave, shown schematically in Fig. 3, which consists of an outward-propagating part, only. The derivation is shown in Appendix B2, and the solution for the pressure perturbation $p^{\prime }=(c_p/\mathrm{R})(p_{\mathrm{0}}/{\mathit{\pi }}_{\mathrm{0}}){\mathit{\pi }}^{\prime }$ associated with the sound wave reads

$$\begin{array}{}\text{(18)}& \begin{array}{rl}& p^{\prime }(x,t)=\\ & \mathit{\delta }p\mathit{\{}{\displaystyle \frac{x-c_{\mathrm{s}}t}{x}}\sin \left(\mathit{\pi }{\displaystyle \frac{x-b_{\mathrm{1}}-c_{\mathrm{s}}t}{b_{\mathrm{2}}-b_{\mathrm{1}}}}\right)\sin \left(\mathrm{2}\mathit{\pi }n{\displaystyle \frac{x-b_{\mathrm{1}}-c_{\mathrm{s}}t}{b_{\mathrm{2}}-b_{\mathrm{1}}}}\right)\\ & +{\displaystyle \frac{b_{\mathrm{2}}-b_{\mathrm{1}}}{x}}\left[{\displaystyle \frac{\sin \left(\mathit{\pi }\left(\mathrm{2}n-\mathrm{1}\right)\frac{x-b_{\mathrm{1}}-c_{\mathrm{s}}t}{b_{\mathrm{2}}-b_{\mathrm{1}}}\right)}{\mathrm{2}\mathit{\pi }\left(\mathrm{2}n-\mathrm{1}\right)}}\right.\\ & \left.-{\displaystyle \frac{\sin \left(\mathit{\pi }\left(\mathrm{2}n+\mathrm{1}\right)\frac{x-b_{\mathrm{1}}-c_{\mathrm{s}}t}{b_{\mathrm{2}}-b_{\mathrm{1}}}\right)}{\mathrm{2}\mathit{\pi }\left(\mathrm{2}n+\mathrm{1}\right)}}\right]\mathit{\}}\left[\mathrm{\Theta }(x-b_{\mathrm{1}}-c_{\mathrm{s}}t)\right.\\ & \left.-\mathrm{\Theta }(x-b_{\mathrm{2}}-c_{\mathrm{s}}t)\right],\end{array}\end{array}$$

where $\mathit{\delta }p=(c_p/\mathrm{R})(\mathit{\delta }T/T_{\mathrm{0}})p_{\mathrm{0}}$ denotes the pressure amplitude of the wave, determined by a temperature amplitude δT in our implementation, $x=|\mathit{X}-\mathit{B}|$ is the distance from the center of the spherical wave B, $x=b_{\mathrm{1}}>\mathrm{0}$, and $x=b_{\mathrm{2}}>b_{\mathrm{1}}$ are the radial boundaries within which the wave has a non-vanishing amplitude at t=0 (see Fig. 3), and n is the number of wave crests (in this work, we consider only n=1). In addition, Θ denotes the Heaviside step function defined as (e.g., Bronstein et al., 2001)

$$\begin{array}{}\text{(19)}& {\displaystyle }{\displaystyle }\mathrm{\Theta }\left(\mathit{\xi }\right)=\left\{\begin{array}{ll}\mathrm{0},& \text{for}\mathit{\xi }<\mathrm{0}\\ \mathrm{1},& \text{for}\mathit{\xi }\ge \mathrm{0}\end{array}\right..\end{array}$$

The solution (Eq. 18) is only valid until the first reflection occurs. Of course, sound waves have been thoroughly investigated in the literature and solutions to the sound wave equation are all but new (e.g., Kirchhoff, 1876), but since their propagation in combination with the method of Läuter et al. (2005) involves potentially all parts of a dynamical core (except for gravity), we found this test useful, also in view of the relative scarcity of test cases dedicated to deep-atmosphere dynamical cores in the literature. In order to highlight the effect of the spherical curvature (on the model grid), the radius of Earth can be rescaled a→η₁a, (η₁<1). This is the small-Earth approach proposed by Wedi and Smolarkiewicz (2009) (see also Baldauf et al., 2014; Ullrich et al., 2014). The model time step is rescaled accordingly in order to account for the correspondingly smaller mesh size of the horizontal grid. Furthermore, it might be advantageous to rescale the angular velocity of the Earth, Ω→η₂Ω, in order to control the velocity $\mathit{v}=-{\mathit{v}}_F$ with which the sound wave is advected. Further details on the implementation can be found in Appendix B2. Apart from that, we followed closely the guidelines given by Baldauf et al. (2014) for the implementation into UA-ICON.

Figure 3 Schematic illustration of the initial state of the spherical sound wave (where the blue and red color shading indicates the positive and negative values of the pressure perturbation associated with the sound wave). The small Earth (without topography) is depicted in gray, and the black circles represent the model bottom and top.

Figure 4 Pressure perturbation associated with the spherical sound wave for a height–longitude cross section at the Equator and at time t=60 s. (a, b) The numerical solution from a simulation with UA-ICON on a R2B6L180 grid is in color (where L denotes the number of vertical grid layers), isolines depict the analytical solution (solid lines denote positive values, dashed lines negative values; the zero contour is omitted). The parameters of the sound wave are listed in Table 2. (a, c) Configuration without rotation (Ω=0). (b, d) Configuration with rotation ($\mathrm{\Omega }={\mathit{\eta }}_{\mathrm{2}}\cdot \mathrm{7.29}\times {\mathrm{10}}^{-\mathrm{5}}$ rad s⁻¹ $=\mathrm{6.84}\times {\mathrm{10}}^{-\mathrm{4}}$ rad s⁻¹). (c, d) Pressure difference from subtracting the numerical from the analytical solution for the two respective configurations in the upper row. (Due to the compression of the circular sector into the rectangular plot, the sound wave appears to have an oval shape.)

Table 2 Parameters used for the sound wave test case with UA-ICON. (Where necessary, a comma separates different values for the two configurations differing in their angular velocity Ω.)

We envisaged two concrete test configurations: the first without rotation (η₂=0), to simulate the sound wave propagation in the absolute frame, and the second with rotation (η₂>0), to test if the dynamical core is able to maintain the balance of the background velocity in the rotating frame $\mathrm{d}{\mathit{v}}_F/\mathrm{d}t+\mathrm{2}\mathbf{\Omega }\times {\mathit{v}}_F+\mathbf{\Omega }\times [\mathbf{\Omega }\times (\mathit{X}-\mathit{A}\left)\right]=\mathrm{0}$, so as to advect the sound wave in a shape-conserving way. Further parameter settings are listed in Table 2. The temperature amplitude of the sound wave δT was chosen small enough that its non-linear dynamics as computed by the dynamical core is negligibly small (an initial amplitude of δT=0.1 K appears large when compared, e.g., to typical values used in Baldauf et al. (2014), but even for larger amplitudes, ∼1 K, the numerical solution of UA-ICON was in relatively good agreement with the analytical solution of the linearized equations).

A height–longitude cross section at the Equator of the numerical solution from UA-ICON, the analytical solution, as well as the difference between the two for both configurations are shown in Fig. 4, shortly before the periphery of the sound wave would impinge on the bottom and top boundaries. The angular velocity of the second configuration was chosen such that the center of the sound wave, B, would be advected by the zonal background wind $u\left(\mathit{B}\right)=-u_F\left(\mathit{B}\right)=-\mathrm{\Omega }r\cos \left(\mathit{\phi }\right){|}_{\mathit{B}}=-\mathrm{100}$ m s⁻¹ in the rotating frame. This is close to the value used by Baldauf et al. (2014) in their test scenario (C). We use a time step of $\mathrm{\Delta }t={\mathit{\eta }}_{\mathrm{1}}\cdot \mathrm{13.2}$ s for both configurations. It should satisfy the Courant–Friedrichs–Lewy (CFL) criterion in both configurations. The maximum propagation velocity in the first configuration is $|v{|}_{max}=c_{\mathrm{s}}=\mathrm{317}$ m s⁻¹, whereas in the second configuration it is $|v{|}_{max}=|u_F{|}_{max}+c_{\mathrm{s}}=\mathrm{134}$ m s⁻¹ +317 m s⁻¹ =451 m s⁻¹.

Figure 5 The L² norm and L^∞ norm of the difference between the numerical solution and the analytical solution of the pressure field at the height–longitude cross section at the Equator and at time t=60 s. Left: configuration without rotation (Ω=0). Right: configuration with rotation ($\mathrm{\Omega }={\mathit{\eta }}_{\mathrm{2}}\cdot \mathrm{7.29}\times {\mathrm{10}}^{-\mathrm{5}}$ rad s⁻¹ $=\mathrm{6.84}\times {\mathrm{10}}^{-\mathrm{4}}$ rad s⁻¹). The norms are plotted for the grids: R2B5L90 ($\mathrm{\Delta }x={\mathit{\eta }}_{\mathrm{1}}\cdot \mathrm{78.9}$ km, Δz=1111.1 m, $\mathrm{\Delta }t={\mathit{\eta }}_{\mathrm{1}}\cdot \mathrm{26.4}$ s), R3B5L136 ($\mathrm{\Delta }x={\mathit{\eta }}_{\mathrm{1}}\cdot \mathrm{52.6}$ km, Δz=735.3 m, $\mathrm{\Delta }t={\mathit{\eta }}_{\mathrm{1}}\cdot \mathrm{17.6}$ s), R2B6L180 ($\mathrm{\Delta }x={\mathit{\eta }}_{\mathrm{1}}\cdot \mathrm{39.5}$ km, Δz=555.6 m, $\mathrm{\Delta }t={\mathit{\eta }}_{\mathrm{1}}\cdot \mathrm{13.2}$ s), R3B6L277 ($\mathrm{\Delta }x={\mathit{\eta }}_{\mathrm{1}}\cdot \mathrm{26.3}$ km, Δz=361.0 m, $\mathrm{\Delta }t={\mathit{\eta }}_{\mathrm{1}}\cdot \mathrm{8.8}$ s), and R2B7L360 ($\mathrm{\Delta }x={\mathit{\eta }}_{\mathrm{1}}\cdot \mathrm{19.7}$ km, Δz=277.8 m, $\mathrm{\Delta }t={\mathit{\eta }}_{\mathrm{1}}\cdot \mathrm{6.6}$ s), where L denotes the number of vertical grid layers; Δx, Δz, and Δt are the mean horizontal mesh size, the grid layer thickness, and the time step, respectively. The dashed and solid lines indicate 𝒪(Δx) and 𝒪(Δx²) behavior, respectively.

Shape and amplitude are relatively well captured by the numerical solution in both configurations, with the difference in amplitude to the analytical solution being about 1 order of magnitude smaller than the magnitude of the wave's pressure perturbation itself. However, in the second configuration, where the sound wave is advected westward while radially propagating, the magnitude of the error has increased slightly, and the symmetry of the pressure difference with respect to a vertical axis crossing the center of the sound wave is lost due to the horizontal advection. The amount by which symmetry is lost is a measure for the phase error of the horizontal advection implementation (e.g., Skamarock and Klemp, 2008). The pressure difference in the first configuration not being radially symmetric with respect to the center of the sound wave is probably due to at least three anisotropies between the vertical and the horizontal: first, the horizontal and vertical mesh sizes (Δx=597.8 m and Δz=555.6 m) slightly differ. Second, the extension of a grid cell increases with height in the spherical geometry, and third, the horizontally explicit–vertically implicit scheme employed in the dynamical core of ICON introduces an anisotropy as well.

We repeated the simulation for different grid resolutions and computed the L² norm and L^∞ norm of the pressure difference between the numerical and analytical solutions on the entire circumequatorial height–longitude cross section, of which a part is plotted in Fig. 4, and at time t=60 s, according or in analogy to the formula employed by Baldauf and Brdar (2014, p. 1983). All pressure values entering the computation of the two norms are weighted equally, i.e., no weighting with the cell volume is applied (in which case the two norms would not be with respect to the pressure difference Δp but with respect to a work-like quantity ∝ΔpV). The results and some further information on the employed grids are shown in Fig. 5. In the first configuration, without rotation, the convergence rate is dominated by a second-order behavior, although a relatively small first-order component seems to be present, especially in the case of the L^∞ norm. In the second configuration, with rotation, the convergence rate seems to start with a second-order behavior for the lower grid resolutions and changes into a first-order behavior for the higher resolutions. This is in agreement with the findings of Baldauf et al. (2014) for their test scenario (C) (compare their Fig. 7). The reason for the first-order convergence in the presence of a background wind is still unknown. Nevertheless, we regard the agreement between the analytical and numerical solutions in both configurations as satisfactory, as we assume a critical deficiency in the deep-atmosphere modification of the dynamical core to leave a much more distinct fingerprint in the numerical solution.

3.1.2 Jablonowski–Williamson baroclinic instability test case

The previous test case focused on one particular emergent structure of the atmosphere. However, if we turn our focus to the atmospheric features on the synoptic scale, other structures, such as baroclinic waves, are much more important than sound waves. The Jablonowski–Williamson baroclinic instability test case (Jablonowski and Williamson, 2006) is a standard test to investigate the performance of atmospheric models in representing a key feature of midlatitude dynamics. It consists of a baroclinically unstable atmosphere in hydrostatic and geostrophic balance to which a perturbation is added which triggers the instability. This test case reveals, on the one hand, if the model is able to maintain the hydrostatically and geostrophically balanced background state during the first days of the wave evolution, when its amplitude is still relatively small, and on the other hand, how the model performs in reproducing the amplitude growth of the wave and its shape. However, a disadvantage of this test is that no analytical solution for the problem is known, so the evaluation has to be based on a model intercomparison. Ullrich et al. (2014) have extended this test case for deep-atmosphere dynamical cores and introduced some further improvements to the original formulation. The approach of the small Earth is employed to highlight the differences between the shallow- and the deep-atmosphere dynamics. The rescale factors are ${\mathit{\eta }}_{\mathrm{1}}=\mathrm{1}/\mathrm{20}$ and η₂=20, for the Earth radius and angular velocity, respectively. For the test with UA-ICON, we used a R3B4 grid which provides a horizontal mesh size of $\mathrm{\Delta }\mathit{\phi }={\mathrm{0.95}}^{\circ }$. This is close to the value used for the production of the numerical benchmark solution in Ullrich et al. (2014). The vertical grid is stretched, with layer thicknesses increasing from the model bottom to the model top at 30 km. Following Ullrich et al. (2014), we use n_lev=30 levels; however, the vertical stretching of the ICON grid differs from their formula (28):

$$\begin{array}{}\text{(20)}& {\displaystyle }{z}_j{\displaystyle }=H_{\text{top}}{\displaystyle \frac{\sqrt{\mathit{\mu }{\left(\frac{n_{\text{lev}}-j+\mathrm{1}}{n_{\text{lev}}}\right)}^{\mathrm{2}}+\mathrm{1}}-\mathrm{1}}{\sqrt{\mathit{\mu }+\mathrm{1}}-\mathrm{1}}},\end{array}$$

where z_j denotes the height of the jth interface separating the layer j−1 from the layer j. The index j counts from the model top to the model bottom, with z₁=H_top and $z_{n_{\text{lev}}+\mathrm{1}}=\mathrm{0}$. With a value of μ=15 for the flattening parameter, the lowermost and uppermost layers have a thickness of $z_{n_{\text{lev}}}\approx \mathrm{82}$ m and $H_{\text{top}}-z_{\mathrm{2}}\approx \mathrm{1249}$ m, respectively. For ICON, the formula reads

$$\begin{array}{}\text{(21)}& {\displaystyle }{z}_j{\displaystyle }=H_{\text{top}}{\left\{{\displaystyle \frac{\mathrm{2}}{\mathit{\pi }}}\arccos \left[{\left({\displaystyle \frac{j-\mathrm{1}}{n_{\text{lev}}}}\right)}^{\mathit{\sigma }}\right]\right\}}^{\mathit{\lambda }}.\end{array}$$

A value of 1 is used for the stretching parameter σ, and a value of 3.16 follows for the exponent λ from Eq. (21) if we require a thickness of $z_{n_{\text{lev}}}=\mathrm{100}$ m for the lowermost layer. This setting results in a thickness of $H_{\text{top}}-z_{\mathrm{2}}\approx \mathrm{1969}$ m for the uppermost layer. We assume that the differences between the vertical grids (20) and (21) are negligible, since tests with n_lev=60 and n_lev=120 revealed that the numerical solution is largely converged on the vertical grid of ICON for n_lev=30 (not shown). The results for days 8 and 10 of the simulation with UA-ICON are shown in Fig. 6a, b. In order to study the convergence behavior with regard to the horizontal grid resolution, we doubled the same twice (see Fig. 6c–f). First of all, we can state that UA-ICON is able to maintain the hydrostatic and geostrophic balances of the background state in the first days of the simulation relatively well. This indicates, for instance, that the vertical variation of the gravitational acceleration is adequately implemented. Second, the amplitude and shape of the baroclinic wave, as they show up in the surface pressure in Fig. 6, compare relatively well to the benchmark solution of Ullrich et al. (2014, their Fig. 9), and also to the solution of Wood et al. (2014, their Figs. 4 and 5). However, some differences can be recognized, especially in the tail of the baroclinic wave. The convergence tests revealed that the numerical solution is largely converged with regard to the horizontal resolution in the zonal range from 120 to 240^∘ E (120^∘ W), say. As mentioned before, the tail of the baroclinic wave in the zonal range from about 60 to 120^∘ E shows a greater variation between the different horizontal resolutions (Fig. 6). To see if the resolution R3B6 is converged in that regard, we tested the R3B7 grid as well. However, it developed a numerical instability in the tail region of the baroclinic wave around day 8. The reason for the instability has not yet been clearly identified, but we assume that the non-traditional part of the Coriolis acceleration (see Eqs. 9 and 10) may play a role in its development, as tests, in which this part was switched off, indicate. This might result from the fact that our discretization does not satisfy $[\mathit{v}\cdot (\mathrm{2}\mathbf{\Omega }\times \mathit{v}){]}_{\text{discretized formulation}}=\mathrm{0}$; i.e., it violates momentum and energy conservation. A formulation of the Coriolis acceleration that satisfies these conservation principles is not possible within the discrete formulation of Eq. (1) used in ICON for reasons discussed by Gassmann (2011, 2018). In simulations of the standard Earth, the violation of the conservation principles is typically so small that it does not pose a problem. In contrast, the small-Earth approach seems to amplify the violation to such an extent that numerical instabilities may be triggered under certain conditions (e.g., a relatively high horizontal grid resolution in this case). However, a clear statement about this is difficult, since the stationary analytical atmospheric background state prescribed for this test case follows from balances that include the non-traditional part of the Coriolis acceleration. As a consequence, the model atmosphere undergoes an adjustment process right away from the beginning of the simulation if this part is switched off. This makes it difficult to evaluate if the instability is absent in the considered time period because of the disabled Coriolis acceleration or because of the atmospheric state being changed by the adjustment process. Nevertheless, given that the instability is absent in standard Earth simulations and measures against it could at most treat the symptoms, not the cause (leaving aside an extensive reformulation of the dynamical core that satisfies the aforementioned conservation principles), we decided to postpone further investigations. Apart from that, the comparison of Fig. 6 with the benchmark from Ullrich et al. (2014) makes us confident that our deep-atmosphere implementation is satisfactory for our purposes.

Figure 6 Surface pressure (in hPa) at days 8 (a, c, e) and 10 (b, d, f) from UA-ICON simulations of the Jablonowski–Williamson baroclinic instability test for deep-atmosphere dynamical cores. Results from three different horizontal grid resolutions are shown: R3B4 ($\mathrm{\Delta }\mathit{\phi }=\mathrm{0.95}{}^{\circ }$, $\mathrm{\Delta }t={\mathit{\eta }}_{\mathrm{1}}\cdot \mathrm{96}$ s), R3B5 ($\mathrm{\Delta }\mathit{\phi }=\mathrm{0.48}{}^{\circ }$, $\mathrm{\Delta }t={\mathit{\eta }}_{\mathrm{1}}\cdot \mathrm{48}$ s), and R3B6 ($\mathrm{\Delta }\mathit{\phi }=\mathrm{0.24}{}^{\circ }$, $\mathrm{\Delta }t={\mathit{\eta }}_{\mathrm{1}}\cdot \mathrm{24}$ s), where Δφ denotes the mean horizontal mesh size and Δt is the dynamical time step. The vertical resolution is the same for all simulations: 30 levels (L30) up to a model top at 30 km. The grid is vertically stretched, from a thickness of $\mathrm{\Delta }{z}_{min}=\mathrm{100}$ m for the lowermost level up to $\mathrm{\Delta }{z}_{max}=\mathrm{1969}$ m for the uppermost level.

3.2 Climatological test cases

3.2.1 Simulation setup

For the evaluation of the model climatology, a UA-ICON simulation with the upper-atmosphere physics coupled to the ECHAM physics package has been performed. The deep-atmosphere dynamics is also switched on. The model was integrated for 20 years with climatological boundary conditions: sea surface temperature and sea ice concentration are averaged for each calendar month from the Program for Climate Model Diagnosis and Intercomparison (PCMDI) Atmospheric Model Intercomparison Project (AMIP) dataset (Taylor et al., 1998) version 1.1.2 over 1979–2014; concentrations of radiatively active gases, namely O₃, CO₂, O₂, O, NO, CH₄, and N₂O, are averaged in the same manner from a 35-year HAMMONIA simulation with fixed present-day boundary conditions; concentrations of CFC-11 and CFC-12 are fixed at 214.5 and 371.1 pptv, respectively; the 1865 condition of the tropospheric background aerosol from the MAC-v1 dataset (Kinne et al., 2013) is used; no volcanic or anthropogenic aerosols are used; land-surface parameters for the parameterization of the effects of subgrid-scale orography and for the embedded version of the Jena Scheme for Biosphere Atmosphere Coupling in Hamburg (JSBACH) land-surface model (v4; Giorgetta et al., 2018) are fixed as described by Giorgetta et al. (2018). The total solar irradiance is held constant at 1361.371 W m⁻², and the F10.7 index for the calculation of EUV heating rates is fixed at 150 sfu ($\mathrm{1}\mathrm{sfu}=\mathrm{1}\times {\mathrm{10}}^{-\mathrm{22}}$ $\mathrm{W}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{Hz}}^{-\mathrm{1}}$).

The simulation uses the R2B4 grid, which has a horizontal mesh size of about 160 km. In the vertical, the model uses 120 layers for the altitude range from the surface up to 150 km. Rayleigh damping (Klemp et al., 2008; Zängl et al., 2015) is applied above 120 km with a maximum damping coefficient of 10 s⁻¹ at the top. Such strong damping is necessary to allow for a reasonable computational time step despite the occasionally very large vertical velocities in the thermosphere. The model was integrated with a (physical) time step of 4 min and five dynamical substeps each physical time step. Radiation parameterizations – i.e., the PSrad radiation scheme of ICON, the shortwave radiation in the SRBC and EUV, the non-LTE longwave radiation, the NO heating – are evaluated once every hour, whereas all other parameterizations are evaluated every time step. For all physics parameterizations, we apply the “all-fast” treatment described by Giorgetta et al. (2018). For non-orographic gravity wave drag, a cutoff maximum vertical wavelength of 12 km is applied, thus prohibiting long gravity waves. Disabling these long gravity waves is physically sensible, as they are believed to strongly propagate horizontally and be subject to internal reflection before reaching the mesopause (Hines, 1997b).

Companion simulations have been performed with two ICON configurations using a standard model lid at 80 km, Rayleigh damping (maximum damping coefficient 1 s⁻¹) applied above 50 km, and 100 vertical levels exactly following the lower part of the vertical grid applied in the UA-ICON simulations. In the first configuration (referred to as ICON in the following), the deep-atmosphere dynamics and the upper-atmosphere physics are disabled. All other numerical and physical settings are identical to the UA-ICON run. The second configuration (referred to as ICON(UA)) additionally has the deep-atmosphere dynamics and upper-atmosphere physics enabled. With the help of these two configurations, we can estimate which of the differences between ICON and UA-ICON are due to the vertical extension and which are related to the application of extended physics and dynamics also below 80 km. The most important differences between these test cases are listed in Table 3.

Table 3 Major differences between the climatological test cases denoted UA-ICON, ICON, ICON(UA), and UAphys-ICON. (The abbreviations – “SA” is shallow atmosphere and “DA” is deep atmosphere – are used for the row “dynamics”.)

Finally, a third configuration, denoted UAphys-ICON, has been simulated, which differs from the UA-ICON setup only in the deep-atmosphere modification of the dynamics being switched off (see Table 3). We regard the comparison of UA-ICON and UAphys-ICON as a possibility to quantify the difference between the shallow-atmosphere and the deep-atmosphere dynamics.

When comparing the experiments ICON and ICON(UA), with their model top at 80 km, to UA-ICON, with its model top at 150 km, one has to take into account that the sponge layer affects the dynamics in ICON and ICON(UA) at an altitude range where there is no such impact in the UA-ICON configuration. As demonstrated by Shepherd et al. (1996), the distortion of the model dynamics by a sponge layer can extend even significantly below this layer by about two density scale heights. We have to accept this, since the use of a sponge layer is necessary to alleviate the adverse effects of a rigid model lid, such as wave reflection. Nevertheless, we assume that the model physics and dynamics of interest dominate to the extent that at least qualitative conclusions can be drawn from our comparison.

For the evaluation of the simulation results, a 15-year (2002 to 2016) temperature climatology from observations of the Sounding of the Atmosphere using Broadband Emission Radiometry (SABER) instrument on NASA's Thermosphere Ionosphere Mesosphere Energetics Dynamics (TIMED) satellite is used (v 2.0; Dawkins et al., 2018). A monthly zonal-mean zonal wind climatology is taken from the Upper Atmosphere Research Satellite Reference Atmosphere Project (URAP; Swinbank and Ortland, 2003).

3.2.2 Comparison of simulated and observed climatologies

Figure 7 shows multi-year zonal-mean temperatures for January and July from the UA-ICON and ICON simulations and from SABER. The observed temperature patterns are reasonably reproduced in the simulations for large parts of the stratosphere and mesosphere. UA-ICON simulates the low summer mesopause temperatures and its altitude well. The same is true for the weaker temperature minimum over the respective winter pole (the altitude range enclosed by the 200 K isoline near 100 km, at about 80^∘ N in Fig. 7a and c, and at about 80^∘ S in Fig. 7d and f). UA-ICON and ICON both simulate the stratopause and tropopause fairly accurately, although the winter stratosphere is slightly warmer than suggested by the observations. In addition, the tropical tropopause layer shows an increased thickness in the simulations as compared to the observations (compare, e.g., the equatorial region with temperatures between 180 and 200 K at an altitude of about 15 km between Fig. 7d and f). This may at least partly be related to the absence of stratospheric aerosol or to the prescribed climatological ozone.

Figure 7 Climatological zonal-mean temperature for (a, b, c) January and (d, e, f) July averages from the 20-year simulations with (a, d) UA-ICON, (b, e) ICON, and (c, f) SABER satellite retrievals averaged over the years 2002 to 2016.

Figure 8 Climatological zonal-mean zonal wind for (a, b, c) January and (d, e, f) July averages from the 20-year simulations with (a, d) UA-ICON, (b, e) ICON, and (c, f) the Upper Atmosphere Research Satellite Reference Atmosphere Project (URAP) climatology. The URAP data are provided on log pressure levels, which have been transformed to (geopotential) height levels assuming a scale height of 7 km.

Figure 9 Climatological zonal-mean temperature differences for (a, b) January and (c, d) July averages from the 20-year simulations between (a, c) ICON(UA) (i.e., with extended dynamics and physics, but lid at 80 km) and ICON, and (b, d) UA-ICON (i.e., with model lid at 150 km) and ICON. Only differences statistically significant at the 95 % level using a t test are shown.

Figure 10 As in Fig. 9 but for the differences between UAphys-ICON (using shallow-atmosphere dynamics) and UA-ICON (using deep-atmosphere dynamics).

Zonal-mean zonal wind climatologies are presented in Fig. 8. Patterns simulated by UA-ICON and ICON agree qualitatively with the observation-based URAP climatology. The sign reversals of zonal wind in both hemispheres near the mesopause are simulated in UA-ICON but are in general too strong; i.e., the lower-thermospheric jets are too strong and peak at too-low altitudes. We are currently investigating if this can be adjusted by tuning the non-orographic gravity wave parameterization (without losing too much performance in lower parts of the atmosphere). Concerning stratospheric winds, the winter westerlies are too weak in ICON and UA-ICON, an issue already mentioned in the ICON evaluation by Crueger et al. (2018). While UA-ICON and ICON show very similar biases for the boreal winter jet, the problem is reduced in austral winter. It is no surprise that our UA-ICON simulation, performed with the same settings for the subgrid-scale orography (SSO) parameterization as used by Crueger et al. (2018), shows similar issues. A reduction of the orographic gravity wave sources would reduce this issue in particular in the Northern Hemisphere but has not been implemented by Crueger et al. (2018) as it would deteriorate near-surface winds. Retuning of orographic and non-orographic gravity wave parameters is planned for future model versions.

In order to quantify how the vertical model extension, the upper-atmosphere physics, and the deep-atmosphere dynamics affect the state of the model atmosphere on climatological scales, we examine the zonal-mean temperature difference as one possible measure for this purpose. Figure 9 shows the differences between UA-ICON and ICON, on the one hand, and between ICON(UA) and ICON on the other hand, again for the months of January and July. The differences below about 60 km are very similar in the left and the right panels, indicating that in this region they are mostly related to the extension of the dynamical and physical processes and not to the vertical extension. In most areas, the process extension leads to higher temperatures with the strongest signals of up to about 5 K in the summer middle stratosphere and even stronger in the winter lower mesosphere. Above about 60 km, the patterns of the differences are again similar but the magnitude is stronger for the difference between UA-ICON and ICON, meaning that here also the vertical extension adds a warming in comparison to the standard configuration of ICON. At the uppermost level of comparison, the temperature differences reach several tens of Kelvin. Despite the sponge layer gradually increasing in magnitude above 50 km in the two configurations, ICON and ICON(UA), the comparisons suggest that a vertical model extension beyond 80 km as implemented in UA-ICON even influences simulated climatological means down to at least 60 km.

The difference between UAphys-ICON and UA-ICON (i.e., the difference between shallow-atmosphere and deep-atmosphere dynamics) is shown in Fig. 10. The application of the deep-atmosphere dynamics results in significantly lower temperatures above about 90 km as compared to the shallow-atmosphere dynamics. This difference increases with height up to several tens of Kelvin at the model top and shows only relatively moderate meridional and temporal variation. A possible explanation for this might read as follows: in UA-ICON, the magnitude of the gravitational acceleration decreases with height (see Eq. 11), whereas it is constant in UAphys-ICON. So, at a given height, we would expect a larger air density on average in UA-ICON than in UAphys-ICON. In addition, the heating rates at a certain altitude that result from the parameterizations of radiative processes depend on the optical thickness of the air column above (for shortwave radiation), which in turn depends more or less directly on the air mass in the column above. That is, if we would use the air mass below the model top in a grid cell column as a vertical coordinate, and compare the mean temperatures from UAphys-ICON with UA-ICON in this coordinate, we would expect significantly smaller differences as compared to the comparison in the geometric height coordinate. To test this hypothesis, we computed the global horizontal mean of the temperature, denoted $\overline{T}\left(z\right)$, and the mean of the air mass in a grid cell column between the height z and the model top at height z_top, which we denote $\overline{m}\left(z\right)$. This was done for two simulation variants that are identical to the UA-ICON and UAphys-ICON simulations, except for using a shorter period. Each simulation variant was initialized with an IFS analysis for 1 November and integrated for 3 months, of which the data for January went into our analysis. We did this for the six winters from 2013/2014 to 2018/2019 and averaged over the January months. In addition, we combined the UA physics with the NWP physics instead of the ECHAM physics, given the relatively short integration period. For a more convenient comparison, we translated $\overline{m}$ into the “mass height” $z_m=-H\ln [\overline{m}/{\overline{m}}_{\text{tot}}+(\mathrm{1}-\overline{m}/{\overline{m}}_{\text{tot}}\left)\exp \right(-z_{\text{top}}/H\left)\right]$, which results from the hydrostatic balance (for the shallow-atmosphere dynamics) $\partial \mathit{\rho }/\partial z=-\mathit{\rho }/H$. Here, ${\overline{m}}_{\text{tot}}$ is the global mean total air mass in a grid cell column and we used a constant scale height of $H=\mathrm{R}T/g=\mathrm{7.6}$ km as a representative value for the observed range of $\overline{T}$. The results for $\overline{T}\left(z\right)$ and $\overline{T}\left(z_m\right)$ are shown in Fig. 11. As expected, the difference $\left|\mathit{\right\{}\overline{T}\left(z_m\right){\mathit{\}}}_{\text{SA}}-\mathit{\left\{}\overline{T}\right(z_m\left){\mathit{\}}}_{\text{DA}}\right|$ is significantly smaller than $\left|\mathit{\right\{}\overline{T}\left(z\right){\mathit{\}}}_{\text{SA}}-\mathit{\left\{}\overline{T}\right(z\left){\mathit{\}}}_{\text{DA}}\right|$ above the mesopause region, in the geometric altitude range between about 95 km and about 130 km (between about 110 km and about 130 km in terms of the mass height). Here, “SA” and “DA” stand for the shallow atmosphere and deep atmosphere (i.e., UAphys-ICON and UA-ICON), respectively. Below and above this range, at least one of our assumptions seems to hold no longer. Below 95 km, in the mesosphere, the interplay of the radiative processes with the other physics may change. Above about 130 km, the proximity of the model top might be one reason why the hypothesis becomes invalid. As the air mass above the model top is neglected, the optical thickness in zenith direction approaches zero there. In addition to the comparison of $\overline{T}\left(z\right)$ and $\overline{T}\left(z_m\right)$, we inspected the heating rates from the shortwave radiative processes within the geometric altitude range between 95 and 130 km at the beginning of the simulations, when the model atmosphere is still in its hydrostatically balanced initial state. They are systematically stronger in the UAphys-ICON setup than in the UA-ICON setup (not shown). As mentioned above, this might be attributed to different optical thicknesses of the overlying air. To conclude, a significant part of the temperature difference shown in Fig. 10 might be explained with the aforementioned hypothesis. However, more sensitivity simulations and further diagnostics would be necessary to specify the contributions from the dynamics and the individual physics parameterizations. We think that the differences visible in Fig. 10 are significant enough to provide an argument for the use of the deep-atmosphere dynamics in a non-hydrostatic, geometric-height-based general circulation model that extends into the thermosphere.

Figure 11 (a) Global horizontal mean of temperature for January plotted against the geometric height from UA-ICON (“DA” indicates deep-atmosphere dynamics; solid black curve) and UAphys-ICON (“SA” indicates shallow-atmosphere dynamics; dashed black curve). Panel (b) is the same as (a) but plotted against the mass height as defined in the text. (c) Absolute temperature difference between UAphys-ICON and UA-ICON computed on equal geometric heights (solid black) and on equal mass heights (dashed black).

The examination of zonal climatological means of the temperature and the zonal wind component provides the first part of a comprehensive evaluation of UA-ICON. It served us as one guideline for the general performance of UA-ICON and the identification of deficiencies during the main development phase. This development has now reached a stage which we regard as mature enough to form the first completed version of UA-ICON. Nevertheless, we are still facing significant biases in certain regions of the middle and upper atmosphere, which are especially apparent in the zonal wind component, as mentioned earlier. We assume that a retuning of some parameters of the parameterizations of orographically and non-orographically forced unresolved gravity waves is a promising, although not straightforward, route to reduce the biases.

4 Conclusions

An upper-atmosphere extension of the ICON general circulation model has been presented. This includes the extension of the dynamical core from a shallow-atmosphere to a deep-atmosphere formulation in order to account for the spherical shape of the atmosphere and the gravitational field as well as to account for the non-traditional part of the Coriolis acceleration. In addition, the physics parameterizations have been complemented by processes which become relevant in the rarified air of the upper mesosphere and lower thermosphere. For instance, molecular diffusion takes over the lead role for mixing from turbulence, and various processes are linked to the relatively strong high-frequency solar irradiance in the upper atmosphere: its absorption is a source of heat; air compounds being ionized by the radiation align with the electromagnetic field of the Earth, forming currents that in turn interact with the neutral air flow via ion drag and Joule heating; atoms and radicals, as the products of photolysis, undergo a recombination reaction that is accompanied by chemical heating.

The new implementations, subsumed under the configuration name UA-ICON, have been validated by means of idealized test cases in terms of the dynamical core, and by climate simulations the results of which were compared to satellite observations. Two test cases were performed. The first one follows a method proposed by Läuter et al. (2005) and considers the propagation of a spherically symmetric sound wave in an atmosphere at rest in an absolute frame and so in a state of solid-body rotation relative to the rotating Earth. This requires accounting for the centrifugal acceleration explicitly in the dynamical core. In addition, gravity is switched off in order to simplify the derivation of an analytical solution, the comparison to which allows to quantify, e.g., how well the isotropic sound wave propagation is reproduced by the dynamical core on the anisotropic spherically curved grid, and how well it maintains the solid-body rotation of the air. The second test case is the Jablonowski–Williamson baroclinic instability test case in its formulation for the deep atmosphere by Ullrich et al. (2014). Its focus is on testing the representation of synoptic-scale flows through the simulation of the life cycle of a baroclinic wave on a hydrostatically and geostrophically balanced, zonally symmetric background state of the atmosphere. As with the first test case, the small-Earth approach of Wedi and Smolarkiewicz (2009) is applied to amplify the effects of the spherical curvature. No analytical solution is available for the second test case, so it relies on a model intercomparison. In either test case, UA-ICON showed satisfying performance, and no indication of a severe deficiency of the deep-atmosphere modification of the dynamical core was found.

For the evaluation of the upper-atmosphere physics, four AMIP-type simulations, each spanning 20 years, were performed and compared to temperature and zonal wind climatologies obtained from measurements of the SABER and TIMED satellite instruments and the URAP project, respectively. The first simulation uses a setup which we regard as typical for our envisaged first applications of UA-ICON. The model top is located at a height of 150 km, and the newly implemented upper-atmosphere physics are enabled in addition to the standard physics of ICON. The deep-atmosphere modifications of the dynamical core are applied as well. In the second simulation, the model top is lowered to 80 km, and the sponge layer is adjusted accordingly. Apart from that, the settings are identical to those of the first simulation. The third simulation equals the second simulation but is a standard ICON simulation with the upper-atmosphere physics and the deep-atmosphere modifications of the dynamics switched off. Finally, a fourth simulation is identical to the simulation with a model top at 150 km, except for the deep-atmosphere modification of the dynamical core being switched off. By comparing the four simulation variants, we try to quantify the effects of the vertical extension of the model domain into the lower thermosphere and of the upper-atmosphere physics on the middle and lower atmosphere. In addition, the difference between the deep-atmosphere and the shallow-atmosphere dynamics is quantified.

The temperature climatologies for January and July from the first type of simulation with UA-ICON are generally in good agreement with the observations. For instance, the atmospheric temperature minimum at the summer mesopause is relatively well captured. Some biases remain, however; a slightly too-warm winter stratosphere and too-thick tropical tropopause layer might be due to the absence of stratospheric aerosol in the simulation or result from the employed ozone climatology. When it comes to the climatologies of the zonal wind for January and July, we find again that UA-ICON reproduces the qualitative structure of the wind in that part of the atmosphere observed by URAP. One remaining issue is that the westerly jets in the lower thermosphere are stronger in UA-ICON, roughly by 50 %. The comparison of UA-ICON simulations and two different model configurations with a lower model top at 80 km has shown that the addition of upper-atmosphere physics and dynamics also affects stratospheric temperatures (increasing them by up to 5 K), and that the vertical extension has noticeable effects at least down to about 60 km. The difference between the shallow-atmosphere and the deep-atmosphere dynamics manifests itself in a significant decrease in temperature in the upper-mesosphere–lower-thermosphere region if the latter is applied.

As an important application for UA-ICON, we have in view the investigation of the impact of gravity waves on the global atmospheric circulation with a special focus on the feedback of the middle and upper atmosphere on the tropospheric weather and climate. The high scalability of ICON on massively parallel computers will allow us to employ model configurations with much higher grid resolutions and eventually resolve large parts of the atmospheric gravity wave spectrum instead of parameterizing it.

Appendix A: Model initialization

We assume the data used to define the initial state of the model atmosphere of UA-ICON to originate from hydrostatic models (this is true for the IFS data, as well as for the climatology employed here). The altitude assigned to such a data point is actually the geopotential height. This is no issue if the shallow-atmosphere approximation is made and the geopotential height coincides with the geometric height employed in the dynamical core of ICON. However, in a deep-atmosphere model, the geometric and geopotential heights differ. Therefore, the data initialization of UA-ICON takes place on geopotential heights, to and from which the geometric heights of the grid levels are transformed using

$$\begin{array}{}\text{(A1)}& {\displaystyle }{\displaystyle }{z}_{\mathrm{g}}={\displaystyle \frac{z}{\mathrm{1}+\frac{z}{a}}},z={\displaystyle \frac{z_{\mathrm{g}}}{\mathrm{1}-\frac{z_{\mathrm{g}}}{a}}}.\end{array}$$

The geopotential height below which initial data, e.g., from the IFS model, are available, will be denoted 𝔷_g in the following (we use a value of 𝔷_g=70 km). Climatological temperatures are taken from Fleming et al. (1988), who offer tables with zonally averaged, monthly temperature values, denoted T_F, from mean sea level to a geopotential height of 120 in 5 km intervals. For our current simple approach, these datasets are averaged temporally and meridionally to obtain a single vertical temperature profile (see Fig. A1). Temperature values within the 5 km intervals are computed by a cubic spline interpolation (e.g., Bronstein et al., 2001). Above 120 km, an analytical temperature profile from Bates (1959) (see also Hedin, 1983) is used, which is formally identical to Eq. (7):

Figure A1 Vertical climatological temperature profile from Bates (1959), Hedin (1983), and Fleming et al. (1988), employed to define the initial state of the model atmosphere above the geopotential height $z_{\mathrm{g}}={\mathfrak{z}}_{\mathrm{g}}=\mathrm{70}$ km. $z_{\mathrm{g}}=z_{\mathrm{g},\text{top}}=\mathrm{146.55}$ km is the geopotential height of the model top at a geometric height of $z=z_{\text{top}}=\mathrm{150}$ km, and T_∞=400 K is the temperature for z→∞ (this value is considerably lower than the mean exospheric temperature of about 1035 K, but higher values could be challenging for the numerical stability at the beginning of the simulation).

$$\begin{array}{}\text{(A2)}& {\displaystyle }{\displaystyle }{T}_B\left(z_{\mathrm{g}}\right)=T_{\mathrm{\infty }}+\left[T_{\mathrm{120}\mathrm{km}}-T_{\mathrm{\infty }}\right]\exp \left[-{\displaystyle \frac{z_{\mathrm{g}}-\mathrm{120}\mathrm{km}}{H_B}}\right],\end{array}$$

where $T_{\mathrm{120}\mathrm{km}}=T_F(z_{\mathrm{g}}=\mathrm{120}\mathrm{km})$, and T_∞ is approximately the temperature for the limit z→∞. This limit corresponds to a geopotential height of z_g=a, which follows from Eq. (A1) by multiplying the right-hand side of the first equation with $\mathrm{1}=(\mathrm{1}/z)/(\mathrm{1}/z)$ and applying the limit. The value of T_∞ could be set, for instance, to the mean exospheric temperature of about 1035 K (in practice, we use a value of 400 K, since higher values could challenge the numerical stability in the initial phase of the simulation). For a steady transition between T_F and T_B, the scale height is set to $H_B=(T_{\mathrm{\infty }}-T_{\mathrm{120}\mathrm{km}})/(\mathrm{d}{T}_F/\mathrm{d}{z}_{\mathrm{g}}{)}_{z_{\mathrm{g}}=\mathrm{120}\mathrm{km}}$. In addition, the temperature blending requires extrapolating the temperature data from below 𝔷_g to above, which is done by a simple linear extrapolation:

$$\begin{array}{}\text{(A3)}& {\displaystyle }{T}_{\text{IFS}}(\mathit{\lambda },\mathit{\phi },z_{\mathrm{g}})=T_{\text{IFS}}(\mathit{\lambda },\mathit{\phi },z_{\mathrm{g}}={\mathfrak{z}}_{\mathrm{g}})+\mathit{\gamma }(\mathit{\lambda },\mathit{\phi })\left(z_{\mathrm{g}}-{\mathfrak{z}}_{\mathrm{g}}\right),\end{array}$$

with $\mathit{\gamma }=\mathrm{d}{T}_{\text{IFS}}/\mathrm{d}{z}_{\mathrm{g}}{|}_{z_{\mathrm{g}}={\mathfrak{z}}_{\mathrm{g}}}$. To obtain a statically stable stratification, γ is limited by the dry adiabatic lapse rate: $\mathit{\gamma }\ge -{\mathrm{\Gamma }}_{\mathrm{d}}=-g/c_p$ (e.g., Holton, 2004). The blending reads

$$\begin{array}{}\text{(A4)}& {\displaystyle }\begin{array}{rl}& T(\mathit{\lambda },\mathit{\phi },z_{\mathrm{g}})=T_{\text{IFS}}(\mathit{\lambda },\mathit{\phi },z_{\mathrm{g}})\mathit{\alpha }\left(z_{\mathrm{g}}\right)\\ & +T_{\text{clim}}\left(z_{\mathrm{g}}\right)\left[\mathrm{1}-\mathit{\alpha }\left(z_{\mathrm{g}}\right)\right],\end{array}\text{(A5)}& {\displaystyle }\begin{array}{rl}& \mathit{\alpha }\left(z_{\mathrm{g}}\right)=\\ & \left\{\begin{array}{ll}\mathrm{1},& \text{for}{z}_{\mathrm{g}}<{\mathfrak{z}}_{\mathrm{g}}\\ \frac{\mathrm{1}}{\mathrm{2}}\left[\mathrm{1}+\cos \left(\frac{z_{\mathrm{g}}-{\mathfrak{z}}_{\mathrm{g}}}{H_{\text{blend}}}\mathit{\pi }\right)\right],& \text{for}{\mathfrak{z}}_{\mathrm{g}}\le z_{\mathrm{g}}\le {\mathfrak{z}}_{\mathrm{g}}+H_{\text{blend}}\\ \mathrm{0},& \text{for}{\mathfrak{z}}_{\mathrm{g}}+H_{\text{blend}}<z_{\mathrm{g}}\end{array}\right..\end{array}\end{array}$$

where T_clim is T_F or T_B, respectively, and H_blend is a tunable blending scale height, which allows to control over what distance the transition between T_IFS and T_clim takes place. Currently, we use a value of H_blend=10 km to avoid negative absolute temperatures which could result from Eq. (A3), although somewhat larger values for H_blend do likely satisfy this as well. The blending factor α satisfies $\mathrm{d}\mathit{\alpha }/\mathrm{d}{z}_{\mathrm{g}}{|}_{z_{\mathrm{g}}={\mathfrak{z}}_{\mathrm{g}},{\mathfrak{z}}_{\mathrm{g}}+H_{\text{blend}}}=\mathrm{0}$ in order to guarantee a steady transition at z_g=𝔷_g and $z_{\mathrm{g}}={\mathfrak{z}}_{\mathrm{g}}+H_{\text{blend}}$. Given the temperature field (A4), the other variables above 𝔷_g are determined by the hydrostatic and geostrophic balances. On the one hand, this provides a relatively simple method for their computation, and on the other hand it helps to reduce the magnitude of the dynamic tendencies and therefore the strength of the adjustment process during the first time steps of the numerical integration. The pressure is computed from a numerical integration of the discretized form of $\partial p/\partial z_{\mathrm{g}}+gp/\left(RT\right)=\mathrm{0}$, starting at 𝔷_g with $p(\mathit{\lambda },\mathit{\phi },z_{\mathrm{g}}={\mathfrak{z}}_{\mathrm{g}})=p_{\text{IFS}}(\mathit{\lambda },\mathit{\phi },z_{\mathrm{g}}={\mathfrak{z}}_{\mathrm{g}})$, where the deep-atmosphere-specific terms, underlined in Eq. (10), are neglected to simplify matters. Once temperature and pressure are known, ρ, π, and θ can be diagnosed. The horizontal wind is determined from a blending formally identical to the temperature blending (Eq. A4) (using the same α). The IFS part for the blending above 𝔷_g is a simple linear extrapolation of the horizontal velocity, formally identical to Eq. (A3). The “climatological” part shall satisfy $f_{\mathrm{r}}{\mathit{e}}_{\mathrm{r}}\times {\mathit{v}}_{\mathrm{h},\text{clim}}=-\mathit{\beta }\left(\mathit{\phi }\right)\left({\mathbf{\nabla }}_{\mathrm{h}}p\right)/\mathit{\rho }$, where p and ρ are the hydrostatic pressure and density, respectively. Associated with the thermal wind balance (Zdunkowski and Bott, 2003; Holton, 2004), relatively strong horizontal temperature gradients between 𝔷_g and 𝔷_g+H_blend can cause the magnitude of v_h,clim to increase with height and reach values which violate the CFL stability criterion (Zdunkowski and Bott, 2003; Holton, 2004). To avoid this, v_h,clim is multiplied by a factor $[\mathrm{1}+(z_{\mathrm{g}}-{\mathfrak{z}}_{\mathrm{g}})/H_{v_{\mathrm{h}}}]\exp [-(z_{\mathrm{g}}-{\mathfrak{z}}_{\mathrm{g}})/H_{v_{\mathrm{h}}}]$, with a tuneable decay scale height $H_{v_{\mathrm{h}}}$ (we use a value of $H_{v_{\mathrm{h}}}=\mathrm{10}$ km). The value of this factor and its vertical derivative is 1 and 0 at z_g=𝔷_g, respectively, so as not to affect the continuity of the extrapolated wind at that altitude. Of course this factor causes v_h,clim to violate the geostrophic balance to some extent, but it turned out that this is less severe for the spin-up phase than an initial wind field violating the CFL criterion locally. In addition, the horizontal pressure gradient is multiplied by the factor

$$\begin{array}{}\text{(A6)}& {\displaystyle }{\displaystyle }\mathit{\beta }\left(\mathit{\phi }\right)=\left\{\begin{array}{ll}\mathrm{1},& \text{for}{\mathit{\phi }}_{\text{trop}}<\left|\mathit{\phi }\right|\\ \frac{\mathrm{1}}{\mathrm{2}}\left[\mathrm{1}-\cos \left(\frac{\left|\mathit{\phi }\right|}{{\mathit{\phi }}_{\text{trop}}}\mathit{\pi }\right)\right],& \text{for}\left|\mathit{\phi }\right|\le {\mathit{\phi }}_{\text{trop}}\end{array}\right.,\end{array}$$

in order to reduce its magnitude smoothly to zero towards the Equator, where the geostrophic balance does not apply. ${\mathit{\phi }}_{\text{trop}}>\mathrm{0}{}^{\circ }$ is a tuneable tropical latitude (we use ${\mathit{\phi }}_{\text{trop}}=\mathrm{10}{}^{\circ }$). Finally, the vertical wind above 𝔷_g is computed from a blending, again formally identical to Eq. (A4), of a linearly extrapolated w_IFS and a w_clim=0, in accordance with the hydrostatic balance and the boundary condition of a vanishing vertical wind at the model top.

Appendix B: Sound wave test case

B1 Formulation in the rotating frame

Here, we consider the transformation of a solution to the deep-atmosphere equations, linearized about an atmosphere at rest in the absolute frame, to a frame rotating with angular velocity Ω. For the atmosphere at rest in the absolute frame, ${\mathit{v}}_{\mathrm{a}}={\mathit{v}}_F+\mathit{v}=\mathrm{0}$ holds, where ${\mathit{v}}_F=\mathbf{\Omega }\times \mathit{r}$, v being the velocity observed in the rotating frame, and $\mathit{r}=\mathit{X}-\mathit{A}$ (with the center of the Earth A and an arbitrary point X not coincident with A). Therefore, in the rotating frame, the motion of the air follows

$$\begin{array}{}\text{(B1)}& {\displaystyle }\begin{array}{rl}& \mathit{\rho }\left[{\displaystyle \frac{\mathrm{d}\mathit{v}}{\mathrm{d}t}}+\mathrm{2}\mathbf{\Omega }\times \mathit{v}+\mathbf{\Omega }\times \left(\mathbf{\Omega }\times \mathit{r}\right)\right]=\mathrm{0}\\ \iff & \mathit{\rho }\left[{\displaystyle \frac{\mathrm{d}}{\mathrm{d}t}}\left(-\mathbf{\Omega }\times \mathit{r}\right)-\mathrm{2}\mathbf{\Omega }\times \left(\mathbf{\Omega }\times \mathit{r}\right)+\mathbf{\Omega }\times \left(\mathbf{\Omega }\times \mathit{r}\right)\right]=\mathrm{0}\\ \iff & -\mathit{\rho }\mathbf{\Omega }\times \left({\displaystyle \frac{\mathrm{d}\mathit{r}}{\mathrm{d}t}}+\mathbf{\Omega }\times \mathit{r}\right)=\mathrm{0}\end{array}\text{(B2)}& {\displaystyle }\Rightarrow {\displaystyle \frac{\mathrm{d}\mathit{r}}{\mathrm{d}t}}=-\mathbf{\Omega }\times \mathit{r}=-\mathbf{\Omega }\times \left(\mathbf{1}\cdot \mathit{r}\right)=-\left(\mathbf{\Omega }\times \mathbf{1}\right)\cdot \mathit{r}=-\mathbf{W}\cdot \mathit{r},\end{array}$$

where Eqs. (B1) and (B2) are equivalent. In the last step, we have introduced the identity tensor 1 and the antisymmetric Coriolis tensor $\mathbf{W}=\mathbf{\Omega }\times \mathbf{1}$ to simplify the following considerations (compare Wilson, 1929; Zdunkowski and Bott, 2003). The solution to Eq. (B2) is

$$\begin{array}{}\text{(B3)}& {\displaystyle }{\displaystyle }\mathit{r}\left(t\right)=\exp (-\mathbf{W}t)\cdot {\mathit{r}}_{\mathrm{0}},\end{array}$$

where ${\mathit{r}}_{\mathrm{0}}=\mathit{r}(t=\mathrm{0})$. Introducing $\widehat{\mathbf{\Omega }}=\mathbf{\Omega }/\mathrm{\Omega }$, and $\widehat{\mathbf{W}}=\mathbf{W}/\mathrm{\Omega }$, with $\mathrm{\Omega }=\left|\mathbf{\Omega }\right|$ and using

$$\begin{array}{}\text{(B4)}& \begin{array}{rl}\exp (-\mathrm{\Omega }\widehat{\mathbf{W}}t)=& \mathbf{1}-{\displaystyle \frac{\mathrm{1}}{\mathrm{1}\mathrm{!}}}\widehat{\mathbf{W}}\mathrm{\Omega }t+{\displaystyle \frac{\mathrm{1}}{\mathrm{2}\mathrm{!}}}{\widehat{\mathbf{W}}}^{\mathrm{2}}{\mathrm{\Omega }}^{\mathrm{2}}{t}^{\mathrm{2}}\\ & -{\displaystyle \frac{\mathrm{1}}{\mathrm{3}\mathrm{!}}}{\widehat{\mathbf{W}}}^{\mathrm{3}}{\mathrm{\Omega }}^{\mathrm{3}}{t}^{\mathrm{3}}\pm \mathrm{\cdots },\end{array}\end{array}$$

where ${\widehat{\mathbf{W}}}^{\mathrm{2}}=\widehat{\mathbf{W}}\cdot \widehat{\mathbf{W}}$, for instance, in combination with

$$\begin{array}{}\text{(B5)}& {\displaystyle }{\displaystyle }{\widehat{\mathbf{W}}}^n=\left\{\begin{array}{ll}{\left(-\mathrm{1}\right)}^{\left(n-\mathrm{1}\right)/\mathrm{2}}\widehat{\mathbf{W}},& \text{for odd }n\\ {\left(-\mathrm{1}\right)}^{n/\mathrm{2}}\widehat{\mathbf{1}},& \text{for even }n\end{array}\right.,\end{array}$$

where $\widehat{\mathbf{1}}=\mathbf{1}-\widehat{\mathbf{\Omega }}\widehat{\mathbf{\Omega }}$ and $\widehat{\mathbf{\Omega }}\widehat{\mathbf{\Omega }}$ denote the dyadic product of $\widehat{\mathbf{\Omega }}$ with itself, the exponential factor in Eq. (B3) can be expanded as

$$\begin{array}{}\text{(B6)}& {\displaystyle }{\displaystyle }\exp (-\mathrm{\Omega }\widehat{\mathbf{W}}t)=\widehat{\mathbf{\Omega }}\widehat{\mathbf{\Omega }}+\cos \left(\mathrm{\Omega }t\right)\widehat{\mathbf{1}}-\sin \left(\mathrm{\Omega }t\right)\widehat{\mathbf{W}}.\end{array}$$

So Eq. (B3) can be rewritten in a form which is easier to evaluate (compare Wilson, 1929; Bronstein et al., 2001; Zdunkowski and Bott, 2003):

$$\begin{array}{}\text{(B7)}& \begin{array}{rl}\mathit{r}\left(t\right)=& \widehat{\mathbf{\Omega }}\widehat{\mathbf{\Omega }}\cdot {\mathit{r}}_{\mathrm{0}}+\cos \left(\mathrm{\Omega }t\right)\widehat{\mathbf{1}}\cdot {\mathit{r}}_{\mathrm{0}}-\sin \left(\mathrm{\Omega }t\right)\widehat{\mathbf{W}}\cdot {\mathit{r}}_{\mathrm{0}}\\ & =\cos \left(\mathrm{\Omega }t\right){\mathit{r}}_{\mathrm{0}}-\sin \left(\mathrm{\Omega }t\right)\widehat{\mathbf{\Omega }}\times {\mathit{r}}_{\mathrm{0}}\\ & +\left\{\left[\mathrm{1}-\cos \left(\mathrm{\Omega }t\right)\right]\widehat{\mathbf{\Omega }}\cdot {\mathit{r}}_{\mathrm{0}}\right\}\widehat{\mathbf{\Omega }}.\end{array}\end{array}$$

Next, given, for instance, an analytical solution in the absolute frame, which is spherically symmetric with respect to the point B, therefore depending only on $|\mathit{X}-\mathit{B}|=\left|\right(\mathit{X}-\mathit{A})-(\mathit{B}-\mathit{A}\left)\right|=\mathit{\left\{}\right[(\mathit{X}-\mathit{A})-(\mathit{B}-\mathit{A})]\cdot [(\mathit{X}-\mathit{A})-(\mathit{B}-\mathit{A})]{\mathit{\}}}^{\mathrm{1}/\mathrm{2}}$, we would find the solution unaltered in the rotating frame relative to $\exp (-\mathbf{W}t)\cdot (\mathit{B}-\mathit{A})$ since, e.g.,

$$\begin{array}{}\text{(B8)}& \begin{array}{rl}& \left[\exp (-\mathbf{W}t)\cdot \left(\mathit{X}-\mathit{A}\right)\right]\cdot \left[\exp (-\mathbf{W}t)\cdot \left(\mathit{B}-\mathit{A}\right)\right]\\ & =\left[\left(\mathit{X}-\mathit{A}\right)\cdot \exp (-\mathbf{W}t{)}^{\mathit{⊺}}\right]\cdot \left[\exp (-\mathbf{W}t)\cdot \left(\mathit{B}-\mathit{A}\right)\right]\\ & =\left(\mathit{X}-\mathit{A}\right)\cdot \left[\exp \left(\mathbf{W}t\right)\cdot \exp (-\mathbf{W}t)\right]\cdot \left(\mathit{B}-\mathit{A}\right)\\ & =\left(\mathit{X}-\mathit{A}\right)\cdot \left(\mathit{B}-\mathit{A}\right).\end{array}\end{array}$$

Here, (⋅)^⊺ denotes the transpose, and we have used the identity $\exp (-\mathbf{W}t{)}^{\mathit{⊺}}=\exp (\mathbf{W}t)$ which can be derived from Eq. (B6) and $\widehat{\mathbf{\Omega }}{\widehat{\mathbf{\Omega }}}^{\mathit{⊺}}=\widehat{\mathbf{\Omega }}\widehat{\mathbf{\Omega }}$, 1^⊺=1 and ${\widehat{\mathbf{W}}}^{\mathit{⊺}}=-\widehat{\mathbf{W}}$. In addition, $\exp \left(\mathbf{W}t\right)\cdot \exp (-\mathbf{W}t)=\mathbf{1}$ can be shown to hold, using Eqs. (B6), (B5), and the definition of $\widehat{\mathbf{W}}$. The above is not a mathematically rigorous proof for the general case, but we hope it can help to illustrate the basic idea behind the method. For detailed mathematical examinations, we refer the reader to Läuter et al. (2005); Staniforth and White (2008). In the following section, such a spherically symmetric solution to be rotated is derived.

B2 Derivation of the sound wave solution

Here, we derive the analytical solution for the expansion of sound waves for deep-atmosphere Eqs. (1)–(4) reduced to the reversible processes (i.e., F=0 and Q=0), without gravity (g=0), linearized about an isothermal atmosphere of temperature T₀=const., at rest in the absolute frame (i.e., ${\mathit{v}}_{\mathbf{a}}={\mathit{v}}_F+{\mathit{v}}_{\mathrm{0}}=\mathbf{\Omega }\times \mathit{r}+{\mathit{v}}_{\mathrm{0}}=\mathrm{0}$). In particular, the neglect of gravity leads to a constant base state for pressure π₀=const. (from Eq. 1) and density ρ₀=const. (from Eq. 4), too. Otherwise, the derivation of an analytic solution even for the linearized equations would be very difficult.

Linearization of Eqs. (1) to (4) about this base state leads to

$$\begin{array}{}\text{(B9)}& {\displaystyle }& {\displaystyle \frac{\partial {\mathit{v}}^{\prime }}{\partial t}}=-c_p{\mathit{\theta }}_{\mathrm{0}}\mathbf{\nabla }{\mathit{\pi }}^{\prime },\text{(B10)}& {\displaystyle }& {\displaystyle \frac{\partial {\mathit{\rho }}^{\prime }}{\partial t}}+{\mathit{\rho }}_{\mathrm{0}}\mathbf{\nabla }\cdot {\mathit{v}}^{\prime }=\mathrm{0},\text{(B11)}& {\displaystyle }& {\displaystyle \frac{\partial {\mathit{\pi }}^{\prime }}{\partial t}}+{\displaystyle \frac{\mathrm{R}}{c_v}}{\mathit{\pi }}_{\mathrm{0}}\mathbf{\nabla }\cdot {\mathit{v}}^{\prime }=\mathrm{0},\end{array}$$

resulting in the sound wave equation:

$$\begin{array}{}\text{(B12)}& {\displaystyle }{\displaystyle \frac{{\partial }^{\mathrm{2}}{\mathit{\pi }}^{\prime }}{\partial t^{\mathrm{2}}}}-c_{\mathrm{s}}^{\mathrm{2}}{\mathrm{\nabla }}^{\mathrm{2}}{\mathit{\pi }}^{\prime }=\mathrm{0},\end{array}$$

with the squared speed of sound $c_{\mathrm{s}}^{\mathrm{2}}=(c_p/c_v)\mathrm{R}{T}_{\mathrm{0}}$.

Without gravity and disregarding the spherical boundaries of the atmosphere, there is no longer a distinct direction (or rather a distinct point). Therefore, we seek a solution which is spherically symmetric with respect to an arbitrary point B within the boundaries. Using the expansion of the Laplace operator in spherical coordinates (e.g., Bronstein et al., 2001), we can write

$$\begin{array}{}\text{(B13)}& {\displaystyle }{\displaystyle \frac{{\partial }^{\mathrm{2}}{\mathit{\pi }}^{\prime }}{\partial t^{\mathrm{2}}}}-c_{\mathrm{s}}^{\mathrm{2}}\left({\displaystyle \frac{\mathrm{2}}{x}}{\displaystyle \frac{\partial {\mathit{\pi }}^{\prime }}{\partial x}}+{\displaystyle \frac{{\partial }^{\mathrm{2}}{\mathit{\pi }}^{\prime }}{\partial x^{\mathrm{2}}}}\right)=\mathrm{0},\end{array}$$

with the radial coordinate $x=|\mathit{X}-\mathit{B}|$ (to avoid any confusion with the distance r from the center of the Earth). This equation is valid for x>0. Using the transform

$$\begin{array}{}\text{(B14)}& \stackrel{\mathrm{̃}}{\mathit{\pi }}=x{\mathit{\pi }}^{\prime },\end{array}$$

one finds for x⋅ (Eq. B13) (Nolting, 2004)

$$\begin{array}{}\text{(B15)}& {\displaystyle }{\displaystyle \frac{{\partial }^{\mathrm{2}}\stackrel{\mathrm{̃}}{\mathit{\pi }}}{\partial t^{\mathrm{2}}}}-c_{\mathrm{s}}^{\mathrm{2}}{\displaystyle \frac{{\partial }^{\mathrm{2}}\stackrel{\mathrm{̃}}{\mathit{\pi }}}{\partial x^{\mathrm{2}}}}=\mathrm{0}.\end{array}$$

In addition, we can derive from Eq. (B11)

$$\begin{array}{}\text{(B16)}& \begin{array}{rl}& {\displaystyle \frac{\partial {\mathit{\pi }}^{\prime }}{\partial t}}=\left.-{\displaystyle \frac{\mathrm{R}}{c_v}}{\mathit{\pi }}_{\mathrm{0}}\left({\displaystyle \frac{\mathrm{2}}{x}}{v}^{\prime }+{\displaystyle \frac{\partial v^{\prime }}{\partial x}}\right)\right|x\cdot \\ \Rightarrow & {\displaystyle \frac{\partial \stackrel{\mathrm{̃}}{\mathit{\pi }}}{\partial t}}=-{\displaystyle \frac{\mathrm{R}}{c_v}}{\mathit{\pi }}_{\mathrm{0}}{\displaystyle \frac{\mathrm{1}}{x}}{\displaystyle \frac{\partial x\stackrel{\mathrm{̃}}{v}}{\partial x}}\end{array},\end{array}$$

where ${\mathit{v}}^{\prime }=v^{\prime }{\mathit{e}}_x$, ${\mathit{e}}_x=(\mathit{X}-\mathit{B})/|\mathit{X}-\mathit{B}|$, and $\stackrel{\mathrm{̃}}{v}=x{v}^{\prime }$. Equation (B16) is required for the specification of the initial conditions. Given 1-D wave Eq. (B15), a general solution is of the form (Nolting, 2004)

$$\begin{array}{}\text{(B17)}& {\displaystyle }{\displaystyle }\stackrel{\mathrm{̃}}{\mathit{\pi }}(x,t)=f_{\mathrm{1}}(x+c_{\mathrm{s}}t)+f_{\mathrm{2}}(x-c_{\mathrm{s}}t),\end{array}$$

where f₁ and f₂ denote incoming and outgoing waves of arbitrary shape, respectively. Note that the transform (Eq. B14) induces a boundary condition at x=0, because finite values of π^′ require $\stackrel{\mathrm{̃}}{\mathit{\pi }}(x=\mathrm{0},t)=\mathrm{0}$. Therefore, the outgoing wave for later times is determined by the incoming wave at the origin:

$$f_{\mathrm{2}}(-c_{\mathrm{s}}t)=-f_{\mathrm{1}}(+c_{\mathrm{s}}t)\text{ for }t\ge \mathrm{0}.$$

Since wave Eq. (B15) is of second order in time, initial conditions for $\stackrel{\mathrm{̃}}{\mathit{\pi }}(x,t=\mathrm{0})$ and for its derivative

$$\begin{array}{}\text{(B18)}& {\displaystyle }{\displaystyle }{\left.{\displaystyle \frac{\partial \stackrel{\mathrm{̃}}{\mathit{\pi }}}{\partial t}}\right|}_{t=\mathrm{0}}=c_{\mathrm{s}}\left({\displaystyle \frac{\mathrm{d}{f}_{\mathrm{1}}}{\mathrm{d}x}}-{\displaystyle \frac{\mathrm{d}{f}_{\mathrm{2}}}{\mathrm{d}x}}\right)\end{array}$$

are required. Instead, we require that the incoming wave vanishes (f₁=0) and prescribe the initial wind field:

$$\begin{array}{}\text{(B19)}& \begin{array}{rl}{v}^{\prime }\left(x\right)=& \mathit{\delta }v\sin \left(\mathit{\pi }{\displaystyle \frac{x-b_{\mathrm{1}}}{b_{\mathrm{2}}-b_{\mathrm{1}}}}\right)\sin \left(\mathrm{2}\mathit{\pi }n{\displaystyle \frac{x-b_{\mathrm{1}}}{b_{\mathrm{2}}-b_{\mathrm{1}}}}\right)\\ & \left[\mathrm{\Theta }(x-b_{\mathrm{1}})-\mathrm{\Theta }(x-b_{\mathrm{2}})\right],\end{array}\end{array}$$

where δv is a constant wind amplitude, n denotes the number of wave crests, and $x=b_{\mathrm{1}}>\mathrm{0}$ and $x=b_{\mathrm{2}}>b_{\mathrm{1}}$ are the radii, between which the sound wave has a non-vanishing amplitude at t=0. The Heaviside step function is given by (Bronstein et al., 2001)

$$\begin{array}{}\text{(B20)}& {\displaystyle }{\displaystyle }\mathrm{\Theta }\left(\mathit{\xi }\right)=\left\{\begin{array}{ll}\mathrm{0},& \text{for}\mathit{\xi }<\mathrm{0}\\ \mathrm{1},& \text{for}\mathit{\xi }\ge \mathrm{0}\end{array}\right..\end{array}$$

The function in Eq. (B19) was chosen not least because it satisfies $v^{\prime }(x=b_{\mathrm{1}},b_{\mathrm{2}})=\mathrm{0}$ and $\mathrm{d}{v}^{\prime }/\mathrm{d}x{|}_{x=b_{\mathrm{1}},b_{\mathrm{2}}}=\mathrm{0}$, and guarantees that the derived fields π^′, θ^′, and ρ^′ do likewise, so that discontinuities pose no extra challenge to the numerical simulation. For the implementation of the test case, we preferred to specify the wave amplitude by means of a temperature amplitude δT, instead of the wind amplitude δv, and use the relation $\mathit{\delta }v=(c_v/\mathrm{R})(\mathit{\delta }T/T_{\mathrm{0}})c_{\mathrm{s}}$ for this purpose. From the initial wind field (B19), we can calculate the derivative df₂(x)∕dx by Eqs. (B16) and (B18). Integration of df₂∕dx yields the Exner pressure perturbation (compare Bronstein et al., 2001):

$$\begin{array}{}\text{(B21)}& \begin{array}{rl}& {\mathit{\pi }}^{\prime }(x,t)\\ & ={\displaystyle \frac{\mathrm{R}{\mathit{\pi }}_{\mathrm{0}}\mathit{\delta }v}{c_v{c}_{\mathrm{s}}}}\left\{{\displaystyle \frac{x-c_{\mathrm{s}}t}{x}}\sin \left(\mathit{\pi }{\displaystyle \frac{x-b_{\mathrm{1}}-c_{\mathrm{s}}t}{b_{\mathrm{2}}-b_{\mathrm{1}}}}\right)\sin \left(\mathrm{2}\mathit{\pi }n{\displaystyle \frac{x-b_{\mathrm{1}}-c_{\mathrm{s}}t}{b_{\mathrm{2}}-b_{\mathrm{1}}}}\right)\right.\\ & \left.+{\displaystyle \frac{b_{\mathrm{2}}-b_{\mathrm{1}}}{x}}\left[{\displaystyle \frac{\sin \left(\mathit{\pi }\left(\mathrm{2}n-\mathrm{1}\right)\frac{x-b_{\mathrm{1}}-c_{\mathrm{s}}t}{b_{\mathrm{2}}-b_{\mathrm{1}}}\right)}{\mathrm{2}\mathit{\pi }\left(\mathrm{2}n-\mathrm{1}\right)}}-{\displaystyle \frac{\sin \left(\mathit{\pi }\left(\mathrm{2}n+\mathrm{1}\right)\frac{x-b_{\mathrm{1}}-c_{\mathrm{s}}t}{b_{\mathrm{2}}-b_{\mathrm{1}}}\right)}{\mathrm{2}\mathit{\pi }\left(\mathrm{2}n+\mathrm{1}\right)}}\right]\right\}\\ & \left[\mathrm{\Theta }(x-b_{\mathrm{1}}-c_{\mathrm{s}}t)-\mathrm{\Theta }(x-b_{\mathrm{2}}-c_{\mathrm{s}}t)\right],\end{array}\end{array}$$

or, alternatively, expressed as pressure perturbation:

$$\begin{array}{}\text{(B22)}& {\displaystyle }{\displaystyle }{p}^{\prime }={\displaystyle \frac{c_p}{\mathrm{R}}}{\displaystyle \frac{p_{\mathrm{0}}}{{\mathit{\pi }}_{\mathrm{0}}}}{\mathit{\pi }}^{\prime }.\end{array}$$

The initialization of the remaining thermodynamic fields (here, e.g., ρ and θ) can be chosen arbitrarily as long as the linearized ideal gas law is fulfilled. From Eq. (4), it follows

$$\begin{array}{}\text{(B23)}& {\displaystyle }{\displaystyle \frac{{\mathit{\rho }}^{\prime }}{{\mathit{\rho }}_{\mathrm{0}}}}+{\displaystyle \frac{{\mathit{\theta }}^{\prime }}{{\mathit{\theta }}_{\mathrm{0}}}}={\displaystyle \frac{c_v}{\mathrm{R}}}{\displaystyle \frac{{\mathit{\pi }}^{\prime }}{{\mathit{\pi }}_{\mathrm{0}}}}.\end{array}$$

In our implementation, we used ${\mathit{\rho }}^{\prime }(x,t=\mathrm{0})=c_v{\mathit{\rho }}_{\mathrm{0}}$ (Eq. B21) ∕(Rπ₀) and determined ${\mathit{\theta }}^{\prime }(x,t=\mathrm{0})$ from Eq. (B23).

Of course, the solution outlined above holds only until the wave front, initially at x=b₂, impinges on the solid boundaries, either the model bottom or top, and is reflected. Therefore, the integration time for this test case is limited.

The above solution is what an observer would see in the absolute frame. To transform it into the rotating frame, every point X (including B) has to be rotated according to the formula

$$\begin{array}{}\text{(B24)}& {\displaystyle }{\displaystyle }\mathit{X}\left(t\right)=\mathit{A}+\exp \left(-\mathbf{W}t\right)\cdot \left[\mathit{X}(t=\mathrm{0})-\mathit{A}\right],\end{array}$$

which follows from Eq. (B3). Alternatively, the formulation in Eq. (B7) could be used. Equation (B24) is the coordinate-independent formulation of the rotation. It has to be measured in the coordinate system of choice for the application in practice. As shown in Eq. (B8), the shape of the sound wave relative to B(t) is unaffected by the rotation.

Finally, we will give a short outline for the implementation of this test case into a model. First of all, the gravitational acceleration g has to be set to zero, or at least to a relatively small value. If rotation is switched on (Ω≠0), it is essential to take into account the centrifugal acceleration $\mathbf{\Omega }\times [\mathbf{\Omega }\times (\mathit{X}-\mathit{A}\left)\right]$ explicitly in the dynamical core, as shown in Eq. (B1). If desired, the radius and angular velocity of the Earth are rescaled, a→η₁a and Ω→η₂Ω, and the topography is set to zero. If possible, the slip boundary condition should be applied at model bottom and top. All physics parameterizations should be switched off. After specifying the background state of the atmosphere (at rest in the absolute frame) by T₀ and p₀, and the sound wave parameters B, b₁, b₂, δv or δT, and n, the dynamic fields of the (dry) atmosphere are initialized according to

$$\begin{array}{}\text{(B25)}& {\displaystyle }& {\displaystyle }\mathit{v}\left(\mathit{X}\right)=-{\mathit{\eta }}_{\mathrm{2}}\mathbf{\Omega }\times \left(\mathit{X}-\mathit{A}\right)+v^{\prime }\left(x\right){\displaystyle \frac{\left(\mathit{X}-\mathit{B}\right)}{x}},\text{(B26)}& {\displaystyle }& {\displaystyle }\mathit{\pi }\left(\mathit{X}\right)={\mathit{\pi }}_{\mathrm{0}}+{\mathit{\pi }}^{\prime }\left(x\right),\text{(B27)}& {\displaystyle }& {\displaystyle }\mathit{\theta }\left(\mathit{X}\right)={\mathit{\theta }}_{\mathrm{0}}+{\mathit{\theta }}^{\prime }\left(x\right),\text{(B28)}& {\displaystyle }& {\displaystyle }\mathit{\rho }\left(\mathit{X}\right)={\mathit{\rho }}_{\mathrm{0}}+{\mathit{\rho }}^{\prime }\left(x\right).\end{array}$$

In order to keep the non-linear dynamics of the sound wave as simulated by the dynamical core negligibly small, the non-dimensional wind amplitude $\mathit{\delta }v/c_{\mathrm{s}}=(c_v/\mathrm{R})(\mathit{\delta }T/T_{\mathrm{0}})$ should be small enough. Apart from that, we adhered closely to the guidelines given by Baldauf et al. (2014) in our implementation of the test case into UA-ICON.