The upper-atmosphere extension of the ICON general circulation model

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 5 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 shallowto deep-atmosphere dynamics, and the implementation of an upper-atmosphere physics package. A series of test cases and climatological simulations show that UA-ICON performs satisfactorily and is in good agreement with the observed global atmospheric circulation. 10


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, by which, e.g., momentum carried along their path is deposited to 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.
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 electronically 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 thermodyamical 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.
The outline of the paper is as follows: in section 2 we describe the modifications and extensions to the dynamical core and to the physics parameterizations, followed by a presentation of results from test cases and climate simulations for the evaluation of the new implementations in section 3. We close with a conclusion in section 4.
2 Model extension to the upper atmosphere

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 one, questioning the neglect of the non-traditional part of the Coriolis accleration 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 shallowatmosphere approximation might start to outweigh its benefits, especially on a climatological time scale.Just to give a simple example concerning the mass of air: assuming an isothermal atmosphere in hydrostatic balance, the mass in a grid cell in a deep model atmosphere is m(z) ≈ m 0 (z)(r/a) 2 exp[(z − z g )/H], where m 0 (z) denotes the mass contained by a cell under shallow-atmosphere approximation, z is the altitude of the cell center above a spherical shell with a radius equal to the mean radius of the Earth a = 6371229 m, r = a + z is the distance from the center of the sphere, H is the temperature scale height and z g is the geopotential height motivated in Eq. ( 8) below.The comparison assumes that H and the air density at z = 0 are the same in a shallow and deep model atmosphere.The factor after m 0 is always > 1 above sea level, and already at an altitude of z ∼ 100 km there is about 20% more mass in a grid cell according to this estimate (assuming H ≈ 10 km).This is due to the fact that the gravitational acceleration decreases with height in the deep atmosphere, so that the atmosphere can extend further radially, which finds its expression in the factor exp[(z − z g )/H].In addition, the cell volume in the deep atmosphere increases with height, so that it contains more mass at a given density (the factor (r/a) 2 , compare section 2.1.2).The relevance of this for the dynamics is certainly a separate question, nevertheless, we think the extension of the dynamical core from shallow-to deep-atmosphere dynamics is an important component for UA-ICON.

Model equations
We restrict our considerations to the dry atmosphere hereafter, 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) and the equation of state reads where v = v h + we 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, π = (p/p 00 ) R/cp denotes the Exner pressure relative to the reference pressure p 00 = 1000 hPa, θ = T /π is the potential temperature, φ g denotes the gravitational potential, and c p = 1004.64J K −1 kg −1 , c v = 717.6J K −1 kg −1 , 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 momentum sources, and Q denotes Exner pressure tendencies due to diabatic processes.
We have added the centrifugal acceleration Ω×(Ω×r) in Eq. ( 1), since it is of importance for a test case presented in section 3.1.1.In the standard configuration of ICON, where this term is absent, we have apparent gravity −∇φ g+c (true gravity plus centrifugal acceleration) on the right-hand side of Eq. ( 1), to which the spherical-geopotential approximation is applied: in 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 e r = ∇φ g /|∇φ g |.In case of apparent gravity we continue to write e r = ∇φ g+c /|∇φ g+c |, 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).Equation (3) follows from the budget equation of the heat density ∂(ρs)/∂t + ∇ • (vρs) = Q/T , by formulation in terms of the potential temperature ∂(ρθ)/∂t + ∇ • (vρθ) = Q/(c p π) and multiplication with Rπ/(c v ρθ), where s = s 0 + c p ln(θ/θ 0 ) denotes the mass-specific heat (or entropy) and s 0 , θ 0 are some reference values (e.g., Gassmann and Herzog, 2015).The advection term in Eq. ( 1) is expressed in the so-called 2d-vector-invariant formulation (e.g., Phillips, 1966;Sadourny, 1972;Vallis, 2006) where The background temperature profile of the hydrostatically balanced part is defined by (Zängl, 2012) where z g = z/(1 + z/a) denotes the geopotential height, and T sl = 288.15K as well as T str = 213.15K are characteristic temperature values at sea level (z g = 0) and in the stratosphere, respectively.H scal = 10000 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) with the geopotential φ g = −ga 2 /r + φ g,0 , where φ g,0 denotes some constant, and g = 9.80665 m s −2 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→−a z g = −∞ to The profile (7) allows to integrate the hydrostatic equilibrium analytically.Note, if the centrifugal acceleration is treated in its explicit form Ω × (Ω × r), it is not taken into account in the hydrostatic equilibrium outlined above.
In ICON the governing equations (1) to (3) are measured in a spherical coordinate system with local unit vectors e t (r), e n (r) and e r (r) (forming a right-handed system in this order), where e t and e n 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) 1 .In order to find expressions for the projection of the momentum equation (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,   1 It would be more correct to refer to the horizontal unit vectors as et(r j ) and en(r j ), since they are defined only at the grid triangle edges j.In contrast to the zonal and meridional unit vectors of the geographical coordinate system, e λ (r) and eϕ(r), et and en do not converge to a differentiable continuum for the number of triangles going to infinity, since the jumps in orientation from one edge to another adjacent edge remain, whatever the horizontal grid resolution.
Nevertheless, if computations require spatial derivatives of vectorial quantities, they can be performed in geographical coordinates and the result can projected onto the et(r j ) and en(r j ).So we regard et and en as an "indirectly differentiable" continuum and write et(r) and en(r), in order to simplify matters.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 t is parallel to e λ at the considered location (where the prime indicates this coordinate transformation), we would find e n being parallel to e ϕ (the radial unit vectors of both systems are parallel anyway).The governing equations for the velocity components v = ue λ + ve ϕ + we r (or here v = u e λ + v e ϕ + we r ) can be found in many textbooks (e.g., Gill, 1982;Zdunkowski and Bott, 2003;Holton, 2004;Vallis, 2006).Evaluating them at the equator (ϕ ( ) = 0) yields almost immediately the components of the governing equations for v = v t e t +v n e n +we r in the local edge coordinate systems, if we identify v t with u and v n 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 where is the horizontal mass-specific kinetic energy, and f r = 2Ω sin(ϕ), f t,n = 2Ω cos(ϕ)e ϕ • e t,n are the Coriolis parameters.The values of the projections e ϕ • e t,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  9) and ( 10) are neglected in the shallow-atmosphere equations.
The velocity component v t is reconstructed diagnostically from v n (e.g., through an interpolation by radial basis functions).
At the end of this section, we have a brief look at the expressions Q and F in Eqs. ( 1) and (3).Although this might be a matter of taste, they can be separated into two parts: which the first part is still counted as contribution to the dynamics, whereas the second part summarizes the tendencies from the physics parameterizations, which we will not consider further in this section (if a prognostic equation for turbulent kinetic energy is part of the model, the considerations about Q and F need special care, to avoid "double countings", but we explicitly exclude that case in the following considerations).
Expressions for Q dyn and F dyn typically follow from the budget (and state) equations of the infinitesimal air parcels of the continuum view-point, to which a mass-weighted Reynolds-average (Hesselberg, 1925;Favre, 1969) is applied (together with the assumption that the air flow behaves quasi-ergodically on the averaging scales, e.g., Herbert (1975);Sievers (1982)) where χ is some mass-specific extensive variable with its molecular flux density J χ and source density σ χ .In addition, (•) denotes the Reynolds-average, (•) = (ρ •)/ρ the mass-weighted average, and The turbulent flux density is typically parameterized by some gradient ansatz where the diffusion coefficient tensor K χ may be some function of the state variables.
In the lower atmosphere the magnitude of the turbulent flux density is typically orders of magnitude larger than the magnitude of the molecular flux density, for what reason Jχ might be neglected there.For Q dyn this ansatz usually means some expression similar to if the budget equation for the heat density ρŝ is written in terms of the potential temperature θ.We dropped again the bars and hats in Eq. ( 13), as this has been done right away from the beginning for Eqs.
(1) to (3), which actually also represent averaged equations of the form (11).Although the ansatz (13) can lead to violations of the second law of thermodynamics (Herbert, 1975;Gassmann and Herzog, 2015;Gassmann, 2018), it is probably the most common ansatz and also applied in ICON.
The expression for F dyn has to parameterize the momentum transfer between neighbouring air volumes due to relative motions and typically reads (e.g., Pichler, 1984;Zdunkowski and Bott, 2003) where F denotes the turbulent momentum flux density tensor, 1 is the identity tensor and (•) denotes the tensor transpose.
The three terms on the right-hand side of Eq. ( 14) correspond to the three forms of (linear) relative motion, into which ∇v can be expanded.These are the volume-conserving shear motion (or deformation), the shape-conserving compression/stretching, and the volume-and shape-conserving rotation (Pichler, 1984;Zdunkowski and Bott, 2003), respectively.The coefficients κ 1 , κ 2 , and κ 3 may be functions of the state variables.If F parameterizes the symmetric Reynolds stress tensor v ρv the antisymmetric rotational part of ( 14) with coefficient κ 3 has to be dropped.So by means of the Reynolds average ansatz any flux between the external angular momentum density of an air volume l ext = r × ρv and its internal angular momentum density l int = Θ • ω can be excluded (where Θ and ω are the tensor of the moment of inertia density and the angular velocity, respectively) 2 .This is a standard approximation in atmospheric physics, regardless of the averaging scale (Herbert, 1978).
Budgeting l int as a further "vectorial" extensive variable in addition to the momentum density p = ρv would of course be impossible for today's models, due to the computational costs.
In addition, a physically consistent implementation would actually require to consider a contribution from the corresponding frictional heating to Q dyn which is proportional to −F • •∇v 3 .Its magnitude is generally assumed to be small, nevertheless, it is a cumulative process, since the sign of this term is always positive.Due to its computational costs, this term is neglected in ICON.
We refrained from deep-atmosphere modifications of the dissipative terms Q dyn and F dyn for reasons exemplified on F dyn : the expansion of ( 14) in Cartesian coordinates, as currently implemented in ICON, is still managable with regard to its computational costs.However, if ( 14) is expanded in spherical coordinates, we are faced with numerous additional metric terms, due to the non-vanishing spatial derivatives of the unit vectors e t , e n , and e r (compare e.g.Baldauf and Brdar, 2016).Their rigorous implementation into the dynamical core would cause considerable additional computational costs, which are not affordable for us for the time being.Furthermore, terms structurally similar to Cartesian expansions of ( 14) are implemented in the dynamical core of ICON as explicit numerical damping to ensure stability.Since they represent no physical process, at least not explicitly, it is not clear, if one would take any advantage from a deep-atmosphere modification.Not to mention that the numerical damping typically neglects terms from the full Cartesian expansions of ( 14), so that coordinate-independence might be lost and a transformation into spherical coordinates becomes a more or less subjective plausibility consideration.
2 Here, we regard (or define) turbulence as that part of the internal motions in the considered air volume that has neither a net momentum nor a net angular momentum.So the internal angular momentum is not part of the turbulence.
3 Here •• denotes the double scalar product.Given the two dyadic products ab and cd of the arbitrary vectors a, b, c and d, their double scalar product is defined as: Wilson, 1929;Zdunkowski and Bott, 2003).So the double scalar product of two arbitrary tensors A and B, measured in some normal basis e i , i = 1, 2, 3: A = A ij e i e j , B = B kl e k e l reads: Geosci.Model Dev.Discuss., https://doi.org/10.5194/gmd-2018-289Manuscript under review for journal Geosci.Model Dev. Discussion started: 20 December 2018 c Author(s) 2018.CC BY 4.0 License.

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 deepatmosphere 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 above mentioned 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 permilles 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) with the modification factor parenthesized by braces.The formula (15) enters (10) only implicitly via the hydrostatic background state.For its computation (15) is used only once during model initialization.Next, we consider the gradient in edge normal and tangential directions 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).∆χ/∆t, 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 → l{r/a}.Like all other modification factors, it is precomputed during model initialization for the vertical positions, where the gradients have to be evaluated during run time.
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 n is defined.This allows to compute the vorticity from the v n using the Stokes theorem 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 n is parallel or antiparallel to the cyclonic direction of the integration path.Again, the first factor on the right-hand side of Eq. ( 17) 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 ( 16).The Eqs. ( 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 deepatmosphere modifications of the flux divergences result from the quotient of area and volume A/V for the 5 surfaces of a triangular cell.The modification of a cell volume itself reads where V = A c (z top − z bot ) is the cell volume under shallow-atmosphere approximation, with the cell surface area A c , and r bot = a + z bot , r top = a + z 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 where l p is the length of the primal edges, we find for the flux divergences Here, ξ i,top,bot denotes the scalar product of the flux density and e n , or e r at the respective cell face.The first term on the right-hand side of Eq. ( 21) sums the fluxes across the side faces, where the orientation factor f o is either 1 or −1 according to whether e n 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 areaweighted interpolations are used primarily.These interpolations are assumed to take place on coordinate surfaces z = const., i.e.
planes in 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 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 case of the shallow atmosphere.We interpret this as an interpolation along the coordinate lines λ, ϕ = const.in 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 afore-mentioned 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 to modify the corresponding gradient according to Eq. ( 16).
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.

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 4 .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, which is based on a mean vertical temperature profile (neglecting 4 See https://www.ecmwf.int/en/forecasts/documentation-and-support(accessed October 2018).horizontal and temporal variations) obtained from Bates (1959); Hedin (1983); Fleming et al. (1988).Later on, this could be improved by using a climatology generated by UA-ICON itself.In the following, we will shortly describe the method to blend the initial data below the mesopause region and the prescribed atmospheric state above the mesopause region which is derived from the climatological values.
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 The geopotential height below which initial data, e.g., from the IFS model, are available, will be denoted z g in the following.
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 km in 5 km intervals.For our current, simple approach these data sets are averaged temporally and meridionally, to obtain a single vertical temperature profile.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) where T 120km = T F (z g = 120 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. ( 22) by multiplying the right-hand side of the first equation with 1 = (1/z)/(1/z) and applying the limit.The value of T ∞ could be set, for instance, to the mean exospheric temperature of about 1035 K.For a steady transition between T F and T B , the scale height is set to In addition, the temperature blending requires extrapolating the temperature data from below z g to above, which is done by a simple linear extrapolation with T IFS (λ, ϕ) = T IFS (λ, ϕ, z g = z g ), and γ = dT IFS /dz g | zg=zg .To obtain a statically stable stratification, γ is limited by the dry adiabatic lapse rate: γ ≥ −Γ d = −g/c p (e.g., Holton, 2004).The blending reads  24), although somewhat larger values for H blend do likely satisfy this as well.The blending factor α satisfies dα/dz g | zg=zg,zg+Hblend = 0, in order to guarantee a steady transition at z g = z g and z g = z g + H blend .
Given the temperature field (25), the other variables above z g are determined by the hydrostatic and geostrophic equilibrium.
On the one hand this provides a relatively simple way 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 ∂p/∂z g + gp/(RT ) = 0, starting at z g with p IFS (λ, ϕ, z g = z g ) = p IFS , 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 (25) (using the same α).The IFS part for the blending above z g is a simple linear extrapolation of the horizontal velocity, formally identical to Eq. ( 24).The "climatological" part shall satisfy f r e r × v h,clim = −β(ϕ)(∇ h p)/ρ, 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 z g and z 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 with a tuneable decay scale height H v h (we use a value of H v h = 10 km).The value of this factor and its vertical derivative is 1 and 0 at z g = z 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 in order to reduce its magnitude smoothly to zero towards the equator, where the geostrophic balance does not apply.ϕ trop > 0 • is a tuneable tropical latitude (we use ϕ trop = 10 • ).Finally, the vertical wind above z g is computed from a blending, again formally identical to Eq. ( 25), 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.

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.
The processes are categorized into 3 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) Banks and Kockarts (1973) frictional heating Gill (1982) ion drag and Joule heating Hong and Lindzen (1976) gravity wave turbulent mixing Hines (1997a, b) ultraviolet: Schumann-Runge bands and continuum (O2) Strobel (1978) extreme ultraviolet (N2, O, O2) Richards et al. (1994) non-LTE infrared cooling (CO2, NO, O3) Fomichev andBlanchet (1995, 1998); Ogibalov and Fomichev (2003) infrared cooling at 5.3 µm (NO) Kockarts (1980) chemical heating climatology from HAMMONIA the atmosphere up to the thermosphere (1.7 × 10 −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.This increases computational efficiency significantly.

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 are parameterized in the UA package following Huang et al. (1998) and Banks and Kockarts (1973), as in HAMMONIA.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 75 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, b) parameterization included in the ECHAM package, but it is switched off in standard ICON.In UA-ICON we enable the calculation and pass the computed turbulent diffusion coefficient to the turbulent mixing subroutine to account for gravity wave-induced turbulent mixing.

Radiation
In standard ICON, the PSrad radiation package (Pincus and Stevens, 2013) is employed, which is itself an extension to the RRTMG model (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 60 km, ultraviolet solar forcing for the O 2 Schumann-Runge bands (SRB; 175 nm to 205 nm) and continuum (SRC; 125 nm to 175 nm) is calculated based on the model of Strobel (1978).Efficiency factors multiplied to the SRBC 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 nm 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 to 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 to 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 section 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 3 and CO 2 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 2 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.
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.

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 3 Model evaluation

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 section 2.1.2,and at testing the metric terms and the complete Coriolis acceleration in the components of the momentum equation ( 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 ( 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.

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 A1, 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 v = −v F = −Ω × (X − 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) 2 (X − 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 0 and temperature T 0 .As long as the sound waves, propagating with the speed of sound c s = (c p /c v )RT 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 A2, and the solution for the pressure perturbation p = (c p /R)(p 0 /π 0 )π associated with the sound wave reads The solution ( 28) 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 → η 1 a, (η 1 < 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 Ω → η 2 Ω, in order to control the velocity v = −v F with which the sound wave is advected.Further details on the implementation can be found in appendix A2.Apart from that, we followed closely the guidelines given by Baldauf et al. (2014) for the implementation into UA-ICON.
We envisaged two concrete test configurations: The first without rotation (η 2 = 0), to simulate the sound wave propagation in the absolute frame, and the second with rotation (η 2 > 0), to test if the dynamical core is able to maintain the balance of the background velocity in the rotating frame dv F /dt + 2Ω × v F + Ω × [Ω × (X − A)] = 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 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 one 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.8m 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 2 -norm and L ∞ -norm of the pressure difference between the numerical and analytical solutions on the entire circum equatorial 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 andBrdar (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 behaviour, although a relatively small first-order component seems to be present, especially in case of the L ∞ -norm.In the second configuration, with rotation, the convergence rate seems to start with a second-order behaviour for the lower grid resolutions, and changes into a first-order behaviour 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.

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 that 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 η 1 = 1/20 and η 2 = 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 ∆ϕ = 0.95 • .This is close to the value used for the production of the numerical benchmark solution in Ullrich et al. (2014).
The vertical grid is streched, with layer thicknesses increasing from the model bottom to the model top at 30 km.Following behavior, we doubled the horizontal grid resolution twice (see Fig. 6, the middle and lower rows), and did likewise for the vertical resolution in case of the horizontal grid R3B5 (see Fig. 7).
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 5 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 vertical resolution (Fig. 7).This is true for the horizontal resolution as well in the zonal range from 120E to 240E (120W), say.As mentioned before, the tail of the baroclinic wave in the zonal range from about 60E to 120E 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 is not clear yet, but two elements may contribute to it.
First, the non-traditional part of the Coriolis acceleration and the metric terms in the expansion of the velocity advection in a spherical coordinate system (see Eqs. ( 9) and ( 10)) may play a role in the development of the instability, as tests, in which they were switched off, indicate.In the following, we limit the consideration to the Coriolis acceleration as a representative.
We consider Gibbs equation for an infinitesimal volume of dry air, fixed in shape and position in the observer frame (compare Falk, 1968(compare Falk, , 1990;;Straub, 1996) where e, p and s are the energy, momentum and heat densities, respectively, and µ denotes the "chemical" potential.The (energy budget) for the overall system "model atmosphere" where E, P , m and S denote the amounts of energy, momentum, mass and heat contained by the respective subsystem "cell" (index i), or "interface region" (index j).In addition, δ denotes the difference between two successive states, and a quantity δX, with X = E, P , . . .will be called a "flux" in the following.That is, we ragard the P j as complete vectorial degrees of freedom, despite typically referring to them as "components" more or less intentially. 5That the volumes of the cells and the interface regions overlap, is irrelevant here, since the flux of energy δE j is solely associated with the flux of momentum δP j , but not associated with the flux of extensive variables contained by the cells.Likewise the energy flux δE i is only associated with the fluxes of the extensive variables contained by the cells (δm i and δS i ), but not with the flux of momentum, which is contained by the interface regions, so there is no double counting.If the momentum vector P j , contained by the interface region, could point into an arbitrary direction, the formulation of a Coriolis acceleration that changes the direction of P j , but leaves its absolute value unchanged, is possible.However, this does not hold in case of the C-grid staggering of the velocity.In terms of Eq. ( 31), with cells and interface regions as two distinct subsystems, it translates into the definition that the interface regions are capable of containing only momentum that is perpendicular to the interface (maybe in an averaged sense, if the interface is not a plane).
This fixed direction of P j (and v j ) makes it impossible, to satisfy v j • (2Ω × P j ) = 0. Changing from local to global energy conservation, i.e. j v j •(2Ω×P j ) = 0, may be a possible way out.However, we could not follow this alternative route in our implementation of the non-traditional part of the Coriolis acceleration, since it turned out to be very difficult on the triangular grid, if possible at all.So our current implementation conserves energy neither locally nor globally, at least formally.This may affect numerical stability negatively.Especially, since the magnitude of the terms of the non-traditional Coriolis acceleration, such as −wf t on the left-hand side of Eq. ( 9) is increased significantly through the small-Earth approach.Not only due to the factor η 2 = 20 in f t = 2η 2 Ω cos(ϕ)e ϕ • e t , but also since the magnitude of the vertical wind in the baroclinic wave about day 8 is at least one order of magnitude larger than in case of the standard Earth (η 2 = 1, not shown).In addition, the magnitude of the vertical wind increases with increasing horizontal resolution, e.g., roughly by a factor of 2 from R3B4 to R3B6 (not shown).
A second element that may contribute to the observed numerical instability, is the apparent presence of a natural instability in the deep-atmosphere test case that seems to be absent in the shallow-atmosphere version.This instability develops in the vicinity of the equator on a shorter time scale than the baroclinic instability.Typical methods to investigate such instabilities are the wave-dynamical stability analysis of the atmospheric state and the parcel-dynamical stability analysis (e.g., Ertel et al., 1941;Gill, 1982;Holton, 2004;Vallis, 2006).The former may be regarded as a global analysis, where a wave-like perturbation is added to the atmospheric state in the region of interest and the evolution of its amplitude is investigated within the framework 5 An alternative, maybe more common model interpretation is, to regard the v j defined on the interfaces, as the components of an "imaginary" cellbased velocity vector v i , to which P i = v i m i would be the corresponding momentum contained by the cell.So Eq. ( 31) would be replaced by where the v i and δP i depend on the v j and δP j in some way.However, this alternative interpretation entails other problems on a triangular grid, as explained by Gassmann (2011). of the linearized governing equations.In contrast, the parcel-dynamical stability analysis is a local analysis.It investigates, how a displaced fluid parcel behaves due to the state of its closest vicinity.In unstable regions, the displacement would increase exponentially.Due to the differences between the two ansatzes, it is likely possible that an atmospheric region that is unstable for some part of the wave spectrum according to the wave-dynamical analysis, is stable according to the parcel-dynamical analysis and vice versa.Here, we make use of the parcel-dynamical analysis according to Ertel et al. (1941) (compare also Shutts and Cullen, 1987) for basically three reasons.We are only interested in the information, if the atmospheric state is unstable in some regions and how strong this instability may be (which, admittedly, may depend implicitly on the chosen analysis/measure).Our focus is not on the spatial structure of the developing instabilities.Second, the instabilities, as we observe them in the simulations are relatively small-scale and start to develop in the close vicinity of the equator, so that we assume the local nature of the parcel-dynamical analysis to be suitable here.Finally, it is a relatively cheap analysis with regard 10 to computational costs, compared to a wave-dynamical analysis.According to Ertel et al. (1941), the evolution of a parcel displacement δr is governed by the equation where the coefficients a, b and c depend on the background state of the atmosphere, the stability of which shall be investigated.Detailed expressions for them are given in appendix B. Now, we limit the parcel displacement to the form δr(t) = δr exp(−iωt), so that Eq. ( 32) can be written in the form (Ertel et al., 1941) where we multiplied Eq. ( 32) by −1.
. By means of the scaling in the last step, we measure ω in the units of the frequency of a model day η 2 Ω, which may be easier to interpret.
A detailed analysis of the roots of this cubic equation in ω2 is beyond the scope of this work.We limited our analysis to find the smallest real root ω2 min (a cubic equation should have at least one real root (Bronstein et al., 2001)).Regions, in which this root is negative, should be unstable with respect to parcel displacements.We applied the analysis to the zonally symmetric background state of the atmosphere of this test case.An approximation of the root ω2 min is determined by an iterative algorithm.Some caution may be in order with regard to the accuracy of this approximation, since cubic equations can be quite challenging for root finding algorithms, such as the Newton's method.However, sensitivity tests, in which we changed the number of iterations, for instance, indicate that the solution is robust.The result is shown in Fig. 8 for the deep-atmosphere version of this test case (Fig. 8a) and for the shallow-atmosphere version (Fig. 8b).The dipole structures in about the northern and southern mid-latitudes that are visible in Fig. 8b, are present in the deep-atmosphere version as well, but we changed the color scale of Fig. 8a, to highlight the most prominent difference between the two test case versions, which is located in the atmospheric column above the equator region.The dark blue and purple region in Fig. 8a, in a height of about 5 km seems to be the most unstable region in terms of ω2 min .There, parcel displacements accelerate on a typical time scale of a model day, which is about three orders of magnitude larger than in the unstable regions over the equator in case of the shallow-atmosphere version.And it is exactly this dark blue region, where we observe the evolution of small-scale structures, for instance, in the field of the horizontal velocity divergence (Fig. 8c).This fine-grained "belt" spanning the equator is absent in the shallow-atmosphere test case (Fig. 8d).We see only larger-scale structures, varying in meridional direction, whose amplitude is about four orders of magnitude smaller than the amplitude of the afore-mentionde instabilities in the deep-atmosphere version.They are likely the expression of an adjustment process that results from the atmospheric background state, specified in the test case, being an analytical solution to the differential equations of the hydrostatic and geostrophic balances, but only an approximate solution with regard to the respective difference equations of the model.In case of the shallow-atmosphere, we determined the second and third root of Eq. (33) as well.Both are real, greater than zero and greater than ω2 min everywhere, so we can be relatively confident that no instability is missed in that test case version by limiting the analysis to ω2 min .The small-scale instabilities in the deep-atmosphere simulation gradually spread in meridional direction, and may be recognized in the southern half of the surface pressure plots of Fig. 6 as well.It might be the presence of these instabilities on top of the baroclinic instability that triggers or contributes to the development of the numerical instability in our simulation.One may think that the magnitude of the observed natural instabilities could be reduced by scaling the explicit numerical damping of the model, or part of it, with some factor , say, where n may be 1 or 2, and γ is a tunable coefficient.This is, however, unfeasible, since the computation of the coefficients a, b and c in Eq. ( 33), as well as the iterative determination of ω2 min are way too expensive in a three-dimensional, time-dependent simulation.Although a more in-depth analysis of the instability is beyond the scope of this work, we note that it is probably not buoyancy-driven.To show this, we consider the vertical component of the prognostic  12) in Ullrich et al., 2014), we find Ullrich et al. (2014), N is real everywhere for the atmospheric background state (see their Fig.1d).So the instability does not originate from this side, although N has its smallest values in a region that overlies the dark blue region in Fig. 8a.
Since the numerical instability is limited to relatively high horizontal grid resolutions under the extreme conditions of the small-Earth approach, but absent, if we change to the standard Earth, we regard it as acceptable.Apart from that, the comparison of Figs. ( 6) and ( 7) with the benchmark from Ullrich et al. (2014) make us confident that our deep-atmosphere implementation is satisfactory for our purposes.

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, except for the height-dependent gravitational acceleration, due to an unidentified bug in the current version of the interface between the extended dynamical core and the ECHAM physics.We assume that this limitation will not affect the results severely.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 PCMDI AMIP dataset (Taylor et al., 1998) (Kinne et al., 2013) is used; no volcanic or anthropogenic aerosols are used; land-surface parameters for the parameterization of the effects of sub-grid scale orography and for the embedded version of the 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 −2 , and the F10.7 index for the calculation of EUV heating rates is fixed at 150 sfu (1 sfu = 1 × 10 −22 W m −2 Hz −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 −1 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 5 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 parameterization 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 −1 ) 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 atmospheric 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.
For evaluation, 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).

Comparison of simulated and observed climatologies
Figure 9 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.However, UA-ICON's winter mesopause is characterized by a vertically extended region of low temperatures and misses the observed distinct temperature minimum near 100 km.Further sensitivity tests will be needed to identify, if this may be improved by further tuning of the gravity wave parameterizations or may rather be caused by oversimplifications in other physical parameterizations, e.g. the treatment of chemical heating.
UA-ICON and ICON both simulate the stratopause and tropopause fairly accurately, although the winter stratosphere is slightly warmer than suggested by the observations and the tropical tropopause is too cold, which may at least partly be related to the absence of stratospheric aerosol or to the prescribed climatological ozone.Differences between UA-ICON and ICON on the one hand, and between ICON(UA) and ICON on the other hand are presented in Fig. 10, again for the months of January and July.It is clear that the differences below about 60 km are very similar, indicating that in this region they are mostly related  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 results from 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 three 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 URAProject, 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 (except for the vertical variation of the gravitational acceleration).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, finally, 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.By comparing the three simulations, we try to quantify the effect of the upper-atmosphere extension on the middle and lower atmosphere.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.The less pronounced minimum in the winter hemisphere at about 100 km, however, is less well reproduced in its magnitude and vertical extent, which, as we speculate, might be alleviated by a retuning of the gravity wave parameterization.In addition, a slightly too warm winter stratosphere and too cold tropical tropopause 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 satisfactorily in general.One remaining issue is that the magnitude of the jets in the lower thermosphere is too large in UA-ICON, roughly by a factor of two.We Since the wave equation ( A15) is of second order in time, initial conditions for π(x, t = 0) and for its derivative are required.Instead, we require that the incoming wave vanishes (f 1 = 0) and prescribe the initial wind field where δv is a constant wind amplitude, n denotes the number of wave crests, and x = b 1 > 0 and x = b 2 > b 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) or, alternatively, expressed as pressure perturbation 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) follows ) In our implementation we used ρ (x, t = 0) = c v ρ 0 (A21)/(Rπ 0 ) and determined θ (x, t = 0) from Eq. (A23).
Of course, the solution outlined above holds only until the wave front, initially at x = b 2 , 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 and ∇∇ denotes the Hessian, a tensor of second order.In addition, the displacement is assumed to be an adiabatic change of state, so that (δρ/ρ) = (1/κ)(δp/p), with κ = c p /c v .Finally, the pressure of the parcel shall adjust to its new environment, i.e. δp = ∇p • δr.This yields for Eq. ( B2) where we divided Eq. (B2) by ρ.Furthermore, (∇p)(∇p) is the dyadic product of the pressure gradient ∇p with itself, and S is called the stability tensor.It is a symmetric tensor, i.e. S = S.
Eq. ( B3) is of the form A • b = 0, where A denotes a given tensor (operator) and b is a variable vector.In the following, we will use that solutions b = 0 require the determinant of A, denoted by |A|, to vanish where A denotes the adjugate tensor of A, (•) • is the scalar of a tensor 6 , and •• denotes the double scalar product of two tensors 7 .The adjugate is given by (compare Wilson, 1929) where × × denotes the double cross product of two tensors (see footnote 7).Now, our goal is, to find the determinant of the operator on the left-hand side of Eq. (B3), |{•}|, in order to simplify the analysis of that equation.For this purpose we abbreviate 6 Given a dyadic product ab of the two arbitrary vectors a and b, its scalar is defined as (Wilson, 1929;Zdunkowski and Bott, 2003) A ij e i • e j = A ii .If the coefficients of the measurement of the tensor in some coordinate system, A ij , are arranged in an algebraic matrix structure, the value of the scalar of A is equal to the trace of that matrix.For the scalar of the identity tensor, we find: (1)• = (e i e i )• = 1 + 1 + 1 = 3.
7 Given two arbitrary dyadic products ab and cd, typically two kinds of double products can be found in the literature (Wilson, 1929;Pichler, 1984;Sihvola, 1999;Lebedev et al., 2010;Zdunkowski and Bott, 2003;Altenbach, 2012) where we have used that 1 and S are symmetric tensors, and W is an antisymmetric tensor in the second step, and the commutativity of the double cross product in the third step (e.g., S × × W = W × × S).In the following, we use a number of identities, the proof of which is beyond the scope of this work.We have to refer to the literature cited in footnote 7, but typically they can be proven directly using, for instance, the expansion of a tensor in a sum of (three) dyadic products, the laws Following Eq. (B4), the determinant reads then Here, we can use the identity The coefficient c is the determinant of (S) and can be computed from a Laplace expansion or some other rule for the computation of the determinant of a 3 × 3 matrix (e.g.Bronstein et al., 2001).
According to the considerations of Thuburn and White (2013) an application of the above ansatz to the shallow-atmosphere equations might be less straight forward as it may seem.Although the effects may be small, the shallow-atmosphere approximation changes, for instance, the metric of the space into a non-Euclidian one.So the concept of a parcel displacement may be reassessed to some extent.Here, we pass over such considerations and take a simple approach by applying the shallowatmosphere approximations, summarized for instance in White et al. (2005), directly to the components of ({•}) in Eq. (B3).
ζe 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 π = π 0 (r)+π , where −c p θ 0 ∇π 0 − ∇φ g = 0. Thus the right-hand side of Eq. (1) reads − c p θ∇π − ∇φ g = −c p θ∇π − c p θ dπ 0 dr .

Figure 1 .
Figure 1.Schematic illustration of a global horizontal triangular grid as used by ICON (a), and of a grid cell (b).The grid is of R2B0-type (twofold root division, zero bisections).The points in and on the grid cell are shown, where the (prognostic) variables are defined: density ρ, potential temperature θ, and Exner pressure π in the cell center, tangential and normal wind components vt, and vn 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 vt, vn, w > 0 for the cell drawing).
) Geosci.Model Dev.Discuss., https://doi.org/10.5194/gmd-2018-289Manuscript under review for journal Geosci.Model Dev. Discussion started: 20 December 2018 c Author(s) 2018.CC BY 4.0 License.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. ( Geosci.Model Dev.Discuss., https://doi.org/10.5194/gmd-2018-289Manuscript under review for journal Geosci.Model Dev. Discussion started: 20 December 2018 c Author(s) 2018.CC BY 4.0 License.

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

(
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, but the efficiency factor of the shortwave radiation is reduced to 0.23.Between 70 km and 80 km the two heating sources are linearly merged.Geosci.Model Dev.Discuss., https://doi.org/10.5194/gmd-2018-289Manuscript under review for journal Geosci.Model Dev. Discussion started: 20 December 2018 c Author(s) 2018.CC BY 4.0 License.

Figure 3 .
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.
28)where δp = (c p /R)(δT /T 0 )p 0 denotes the pressure amplitude of the wave, determined by a temperature amplitude δT in our implementation, x = |X − B| is the distance from the center of the spherical wave B, x = b 1 > 0, and x = b 2 > b 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 Geosci.Model Dev.Discuss., https://doi.org/10.5194/gmd-2018-289Manuscript under review for journal Geosci.Model Dev. Discussion started: 20 December 2018 c Author(s) 2018.CC BY 4.0 License.etal., 2001) , 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(B) = −u F (B) = −Ωr cos(ϕ)| B = −100 m s −1 in the rotating frame.This is close to the value used by Baldauf et al. (2014) in their test scenario (C).A time step of ∆t = η 1 • 13.2 s fits both, the maximum magnitude of the propagation velocity |v| max = c s = 317 m s −1 in the first configuration, and |v| max = |u F | max + c s = 134 m s −1 + 317 m s −1 = 451 m s −1 in the second configuration.

Figure 4 .
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.Upper row: 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.The parameters of the sound wave are listed in Table 2. Left: configuration without rotation (Ω = 0).Right: configuration with rotation (Ω = η2 • 7.29 × 10 −5 rad s −1 = 6.84 × 10 −4 rad s −1 ).Lower row: pressure difference from subtracting the numerical from the analytical solution for the two respective configurations in the upper row.(The longitude counts positive in eastward direction, in addition the sound wave appearing to have an oval shape is due to the compression of the circular sector into the rectangular plot.

Figure 6 .
Figure 6.Surface pressure (in hPa) at days 8 (left) and 10 (right) 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 (∆ϕ = 0.95 • , ∆t = η1 • 96 s), R3B5 (∆ϕ = 0.48 • , ∆t = η1 • 48 s) and R3B6 (∆ϕ = 0.24 • , ∆t = η1 • 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 ∆zmin = 100 m for the lowermost level up to ∆zmax = 1969 m for the uppermost level.
Coriolis acceleration ∂p/∂t = −2Ω × p + • • • conserves the absolute value of the momentum |p| and is energetically neutral, since it does not contribute to the first term on the right-hand side of Eq. (30), if the equation of state of mechanics v = p/ρ is applied.In ICON, ρ and s are contained by the grid cells.The components of the momentum, or rather its intensive counterpart, the velocity v are defined on the cell interfaces.We may translate this C-grid staggering into the following Gibbs equation Geosci.Model Dev.Discuss., https://doi.org/10.5194/gmd-2018-289Manuscript under review for journal Geosci.Model Dev. Discussion started: 20 December 2018 c Author(s) 2018.CC BY 4.0 License.

Figure 8 .
Figure 8. Upper row: the smallest real root ω2 min of Eq. (33) for the zonally symmetric atmospheric background state of the Jablonowski-Williamson baroclinic instability test for the deep atmosphere (a) and the shallow atmosphere (b).The computation was performed on the grid R3B7L30 (∆ϕ = 0.12 • ).Note the different color scales and the non-linear color scale of the left figure.The one black contour highlighted by thickness is the zero contour.Lower row: the horizontal velocity divergence at height z = 5 km and day 4 for the deep atmosphere case (c) and the shallow atmosphere case (d).A dynamical time step of ∆t = η1 • 12 s was used for the simulations on the R3B7L30-grid.Note the different color scales.
where κ = c p /c v and the • • • stand for further terms that are unimportant for this consideration.If we apply that the prescribed atmospheric background state satisfies the hydrostatic balance ∂p/∂r = −gρ + • • • , where the • • • denote further terms in case of the deep atmosphere (see Eq. ( version 1.1.2over 1979-2014; concentrations of radiatively active gases, namely O 3 , CO 2 , O 2 , O, NO, CH 4 , N 2 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 pptv and 371.1 pptv, respectively; the 1865 condition of the tropospheric background aerosol from the MAC-v1 dataset Geosci.Model Dev.Discuss., https://doi.org/10.5194/gmd-2018-289Manuscript under review for journal Geosci.Model Dev. Discussion started: 20 December 2018 c Author(s) 2018.CC BY 4.0 License.
several tens of Kelvin (max and min values are denoted at the top of each panel).The vertical extension may affect the temperature difference to standard ICON by both the actual vertical extension and by omitting the Rayleigh damping at these altitudes.It is clear from these comparisons that a vertical model extension beyond 80 km as implemented in UA-ICON even influences simulated climatological means down to at least 60 km.Zonal mean zonal wind climatologies are presented in Fig. 11.Patterns simulated by UA-ICON and ICON agree qualitatively 10 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.Experience from test-experiments show that this can very likely be adjusted by tuning the non-orographic gravity Geosci.Model Dev.Discuss., https://doi.org/10.5194/gmd-2018-289Manuscript under review for journal Geosci.Model Dev. Discussion started: 20 December 2018 c Author(s) 2018.CC BY 4.0 License.

Figure 10 .Figure 11 .
Figure 10.Climatological zonal-mean temperature differences for (top) January and (bottom) July averages from the 20-year simulations between (left) ICON(UA) (i.e. with extended dynamics and physics, but lid at 80 km) and ICON, and (right) 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.The max and min differences are indicated at the top of each panel.(The vertical and horizontal axes show the height in km and the latitude in degree, respectively.) Geosci.Model Dev.Discuss., https://doi.org/10.5194/gmd-2018-289Manuscript under review for journal Geosci.Model Dev. Discussion started: 20 December 2018 c Author(s) 2018.CC BY 4.0 License.
are currently investigating, to what extent a retuning of the nonorographic gravity wave parameterization could alleviate the problem without losing accuracy in the middle and lower part of the atmosphere.The comparison of UA-ICON simulations and two different model configurations with a lower model top at 80 Geosci.Model Dev.Discuss., https://doi.org/10.5194/gmd-2018-289Manuscript under review for journal Geosci.Model Dev. Discussion started: 20 December 2018 c Author(s) 2018.CC BY 4.0 License.km has shown that the addition of upper atmosphere physics and dynamics also affects stratospheric temperatures (increasing them by up to 4 K), and that the vertical extension has noticeable effects at least down to about 60 km.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.Geosci.Model Dev.Discuss., https://doi.org/10.5194/gmd-2018-289Manuscript under review for journal Geosci.Model Dev. Discussion started: 20 December 2018 c Author(s) 2018.CC BY 4.0 License.

)
The function (A19) was chosen not least because it satisfies v (x = b 1 , b 2 ) = 0 and dv /dx| x=b1,b2 = 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 δv = (c v /R)(δT /T 0 )c s for this purpose.From the initial wind field (A19) we can calculate the derivative df 2 (x)/dx by Eqs.(A16) and (A18).Integration of df 2 /dx yields the Exner pressure perturbation (compareBronstein et al., 2001) A24) Geosci.Model Dev.Discuss., https://doi.org/10.5194/gmd-2018-289Manuscript under review for journal Geosci.Model Dev. Discussion started: 20 December 2018 c Author(s) 2018.CC BY 4.0 License.beunaffected by the displacement.In the second step, we assumed that Eq. (B1) still holds for the parcel to the extent that we can insert [(B1)−∇p]/ρ into the term with the factor δρ.This is a somewhat vague step, but we use it as a short cut to the same final result as obtained byErtel et al. (1941) in a more elaborate line of argumentation.For the changes in the gradients of the geopotential and the pressure as observed by the displaced parcel, the ansatzes δ(∇χ) = ∇∇χ • δr are used, where χ = φ g , p : (ab)• = a • b = b • a = [(ab) ]•.So the scalar of an arbitrary tensor A, measured in some normal basis e i , i = 1, 2, 3: A = A ij e i e j , reads: (A)• = A ij (e i e j )• = : ab • •cd = (b • c)(a • d) and ab • • cd = (a • c)(b • d), where •, • ∈ {•, ×}.Depending on the combination, the product is either a scalar (for • → • and • → •), a vector (for • → • and • → ×, or vice versa), or a dyadic product (for • → × and • → ×).Here, we make use of the double scalar product of the first kind: ab• •cd = (b • c)(a • d) = (a • d)(b • c) = ba • •dc = (ab) • •(cd) = (d • a)(c • b) = cd • •ab.For the scalar of the product of two dyadic products, we find:(ab • cd)• = [a(b • c)d]• = (b • c)(ad)• = (b • c)(a • d) = ab • •cd.In addition, we need the double cross product of the second kind:ab × × cd = (a × c)(b × d) = (−c × a)(−d × b) = cd × × ab.If we transpose the factors, we find:(ab) × × (cd) = ba × × dc = (b × d)(a × c) = (ab × × cd) .The above applies to tensors as well, since they can be expanded in a sum of dyadic products, e.g., A = a i b i .Geosci.Model Dev.Discuss., https://doi.org/10.5194/gmd-2018-289Manuscript under review for journal Geosci.Model Dev. Discussion started: 20 December 2018 c Author(s) 2018.CC BY 4.0 License. the material derivative by D = d/dt.Since the material derivative in Eq. (B3) acts solely on the displacement δr, but not on the coefficients in {•}, we can treat D as a scalar in the following derivation.We start with the adjugate of { of associativity and commutativity of the respective products, the expansion a× (b × c) = ba • c − ca • b, and a • (b × c) = c•(a×b) = b•(c×a).Since the determinant of the identity tensor is 1, it follows from Eqs. (B4) and (B5) that (1/2)1 × × 1 = 1.In addition,1 × × W = W • 1 − W = W holds,where we have used that W is antisymmetric andW • = (W ) • = −W • = 0. Likewise 1 × × S = S • 1 − S = S • 1 − S, using the symmetry of S. Finally, W × × W = (Ω × 1) × × (Ω × 1) = 2ΩΩ.If we insert them in Eq. (B6), we find

Table 1 .
Physical parameterizations implemented in the UA package of UA-ICON shown with references.

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 Ω.)