Extending the Modular Earth Submodel System (MESSy v2.54) model hierarchy: the ECHAM/MESSy IdeaLized (EMIL) model setup

As models of the Earth system grow in complexity, a need emerges to connect them with simplified systems through model hierarchies in order to improve process understanding. The Modular Earth Submodel System (MESSy) was developed to incorporate chemical processes into an Earth System model. It provides an environment to allow for model configurations and setups of varying complexity, and as of now the hierarchy ranges from a chemical box model to a fully coupled chemistry–climate model. Here, we present a newly implemented dry dynamical core model setup within the MESSy framework, denoted as ECHAM/MESSy IdeaLized (EMIL) model setup. EMIL is developed with the aim to provide an easily accessible idealized model setup that is consistently integrated in the MESSy model hierarchy. The implementation in MESSy further enables the utilization of diagnostic chemical tracers. The setup is achieved by the implementation of a new submodel for relaxation of temperature and horizontal winds to given background values, which replaces all other “physics” submodels in the EMIL setup. The submodel incorporates options to set the needed parameters (e.g., equilibrium temperature, relaxation time and damping coefficient) to functions used frequently in the past. This study consists of three parts. In the first part, test simulations with the EMIL model setup are shown to reproduce benchmarks provided by earlier dry dynamical core studies. In the second part, the sensitivity of the coupled troposphere–stratosphere dynamics to various modifications of the setup is studied. We find a non-linear response of the polar vortex strength to the prescribed meridional temperature gradient in the extratropical stratosphere that is indicative of a regime transition. In agreement with earlier studies, we find that the tropospheric jet moves poleward in response to the increase in the polar vortex strength but at a rate that strongly depends on the specifics of the setup. When replacing the idealized topography to generate planetary waves by mid-tropospheric wave-like heating, the response of the tropospheric jet to changes in the polar vortex is strongly damped in the free troposphere. However, near the surface, the jet shifts poleward at a higher rate than in the topographically forced simulations. Those results indicate that the wave-like heating might have to be used with care when studying troposphere–stratosphere coupling. In the third part, examples for possible applications of the model system are presented. The first example involves simulations with simplified chemistry to study the impact of dynamical variability and idealized changes on tracer transport, and the second example involves simulations of idealized monsoon circulations forced by localized heating. The ability to incorporate passive and chemically active tracers in the EMIL setup demonstrates the potential for future studies of tracer transport in the idealized dynamical model.

1 Introduction

Earth system models continue to incorporate more processes to enable a more complete simulation of the climate system and thus produce the best possible climate projections. In practice, this increases the complexity of model codes as new compartments are added to represent new processes and interactions. However, with models gaining more and more complexity, it becomes difficult to isolate and understand the role of individual processes. This “gap between simulation and understanding in climate modeling” was pointed out by Held (2005), and it was suggested that the way forward is to work with a hierarchy of models with reduced to full complexity. Two recent overview papers (Jeevanjee et al., 2017; Maher et al., 2019) give surveys of current concepts and activities in building hierarchical model systems.

The basic concept in constructing a simplified model is to include only those processes that are (absolutely) relevant for the question to be addressed. Thereby, the behavior of those processes can be isolated in an idealized environment, and the interaction of the limited number of processes chosen can be investigated.

A frequently used idealized model setup for studying global large-scale dynamics is the dry dynamical core model proposed by Held and Suarez (1994, HS94 hereafter). While originally developed and used for testing dynamical cores of atmospheric models, the elegance of the model makes it an ideal tool for dynamical process studies, and it is widely used for this purpose (see Maher et al., 2019, for a review of applications). This Held–Suarez-type model uses the full dynamical core of a general circulation model (GCM) but replaces all thermodynamical processes (e.g., radiation, convection) by relaxation towards a prescribed equilibrium temperature and the surface boundary layer by relaxation of near-surface winds towards zero (as described in detail in Sect. 2). Thus, with this model setup, the thermodynamic forcing of the atmosphere can be easily modified and the response of the large-scale circulation to those isolated modifications can be studied. Examples are changes in equilibrium meridional temperature gradient or thermal damping timescale (Gerber and Vallis, 2007) or changes in surface friction (Chen et al., 2007).

The functions for the equilibrium temperature and relaxation coefficients suggested in HS94 are widely used, and the HS94 model setup was extended to study the dynamics of the stratosphere–troposphere system by modifying the equilibrium temperature of the stratosphere (Polvani and Kushner, 2002, PK02 hereafter) and later by adding topography to include planetary wave generation that is essential for the stratospheric circulation (Gerber and Polvani, 2009). This model setup was used among others to study stratosphere–troposphere coupling (Gerber and Polvani, 2009), the structure of the Brewer–Dobson circulation (Gerber, 2012) and the circulation's response to idealized heating resembling the thermal response to greenhouse forcing (e.g., Butler et al., 2010; Wang et al., 2012). Recently, it was suggested that the forcing of the planetary waves relevant for stratospheric dynamics can also be achieved by inserting diabatic heating in the middle to upper troposphere (Lindgren et al., 2018), which leads to a climatology similar to that of the topographically forced simulations but to changes in the sudden stratospheric warming properties.

While the dry dynamical core model has proven useful in advancing our understanding of the dynamical response to given thermodynamic forcing, the application of the model hinges on a realistic representation of the Earth atmosphere's behavior of the modeled dynamics. Gerber and Polvani (2009) and Chan and Plumb (2009) showed that the strong response of the surface jet location to stratospheric polar vortex changes found in the original study by PK02 resulted from a regime shift of the tropospheric jet. With a changed setup, e.g., by including topography (Gerber and Polvani, 2009), or with enhanced meridional temperature gradients in the winter hemisphere (Chan and Plumb, 2009), the regime-like behavior of the jet location is suppressed, and thus the response of the jet location to stratospheric polar vortex changes is damped strongly. However, the regime shift can re-emerge for experiments with strong additional forcing (e.g., tropical heating, as shown by Wang et al., 2012). Overall, those results indicate that the dynamical response to a given forcing is highly (non-linearly) dependent on the basic state of the model. Whether this sensitivity to the basic state due to dynamical regimes is relevant for the real atmosphere will have to be evaluated with care. If the regime behavior proves to be an artifact of the idealized models, this would impede its application to advance the understanding of dynamical processes of the real atmosphere.

Beyond the purely dry dynamical core models, which are useful to study aspects of the global circulation, a question that motivates the expansion to another level of complexity is the interaction of moisture with large-scale dynamics, either by latent heat release or by its role as a greenhouse gas. Frierson et al. (2006) expanded the dry dynamical core (Held–Suarez) model by adding moisture and convection with latent heat release to the model, including simplified (gray) radiation that is insensitive to water vapor, thus tackling the question of the role of latent heat release for large-scale dynamics. This model setup has been extended by including the radiative effects of ozone in an idealized manner, resulting in a more realistic simulation of stratospheric dynamics (Davis and Birner, 2019). In a step further, the role of water vapor as radiatively active gas is included by using more comprehensive radiation schemes, as done by, e.g., Merlis et al. (2013), Jucker and Gerber (2017) and Tan et al. (2019). In those setups, treatment of radiatively relevant fields as clouds, ozone and aerosol forcing is mostly based on simple assumptions such as constant values.

As stated above, the nature of the hierarchy that is to be constructed depends on the scientific question at hand. Our aim is to study the large-scale dynamical variability of the stratosphere–troposphere system and its response to idealized forcings, and in particular the impact of dynamical variability and forced changes on the transport of passive and chemically active trace gases. The latter is motivated by a variety of research questions on the distribution of trace species in the atmosphere, for example, on how changes in the circulation in a changing climate will affect stratospheric ozone. This question has received a lot of attention recently in the light of observed lower stratospheric ozone trends that are not fully understood (Ball et al., 2018). Another question we aim to tackle with the idealized model is the efficiency of troposphere–stratosphere transport in monsoonal circulation systems via different pathways. The idealized setup allows us to study the role of different transport pathways depending on the details of the forcing of the circulation system. To enable those studies, a well-suited model setup is a dry dynamical core model with the utilities for tracer transport and the possibility to include chosen chemical reactions (simplified to the needs of the user). Therefore, we implement such a model setup within the Modular Earth Submodel System (MESSy; Jöckel et al., 2005) framework, which provides the needed utilities in a modular manner.

Several initiatives are aiming to build modeling frameworks with setups of varying complexity within the same model system (Vallis et al., 2018; Polvani et al., 2017), an approach that will advance both the usability of idealized models as well as the connectedness of the simple and the more complex model setups. In the same spirit, the MESSy framework was developed explicitly with the goal to provide “a framework for a standardized, bottom-up implementation of Earth System Models (or parts of those) with flexible complexity” (see https://www.messy-interface.org/, last access: 29 October 2020). The motivation of the MESSy framework was originally to incorporate chemical processes of varying complexity into an Earth system model. The MESSy framework couples a base model (dynamical core) to submodels that contain the physical parameterizations as well as diagnostics. Among other base models, the ECHAM dynamical core is available in MESSy. The MESSy framework includes model configurations ranging from a zero-dimensional box model of atmospheric chemistry (Sander et al., 2019) to the complex chemistry–climate ECHAM/MESSy Atmospheric Chemistry (EMAC) model, coupled to a deep ocean model (Jöckel et al., 2016). An illustration of a selection of available model complexities is shown in Fig. 1, as a function of the complexity in physical processes/compartments included (horizontal axis) and of the complexity of atmospheric chemical processes included (vertical axis). The lowest complexity on the chemical axis is prescribed concentrations for radiatively active species (e.g., ozone), followed by a simplified parameterization to include effects of methane oxidation on stratospheric water vapor. The chemistry module MECCA (Module Efficiently Calculating the Chemistry of the Atmosphere; Sander et al., 2019) contains a large set of reactions relevant in the troposphere and stratosphere, but it can be configured to the user's needs by choosing any subset of reactions, thus allowing for simplified to very comprehensive chemical setups. The chemical calculations can be performed as a box model (denoted CAABA; see Sander et al., 2019) or within a full general circulation model either without feedback between dynamics and chemistry (the so-called “quasi-chemistry–transport model” (QCTM); see Deckert et al., 2011) or with feedback, i.e., as a full chemistry–climate model (Jöckel et al., 2006, 2010, 2016). Besides the prescribed sea surface temperature setup, a mixed-layer ocean (Dietmüller et al., 2014) or a full ocean model (Pozzer et al., 2011) can be used.

One advantage of the MESSy framework is its modular nature; i.e., individual processes are implemented as independent submodels that can be easily exchanged or complemented by new processes, and each submodel can be easily switched on or off (by namelist choice). Therewith, the hurdle of code modifications to build a model tailored to the necessary complexity is rather low. Moreover, the design of the model system allows the creation of model hierarchies in which the same code can be used in a simple model setup as well as in the full Earth system model. Any developments in model components can be transferred easily up- and downward in the model hierarchy.

Figure 1 Schematic of the MESSy model hierarchy with existing (blue) model setups and the model setup described in this paper (red). Model setups are displayed as function of their complexity in dynamics/physics/compartments (horizontal axis) versus complexity in chemical mechanism (vertical axis). The horizontal axis ranges from (left) a zero-dimensional box model to (middle) models with an atmospheric dynamical core (ECHAM or other implemented dynamical cores in MESSy), but with varying physical complexity, and to (right) models with an additional ocean model (full 3-D or mixed-layer ocean). The vertical axis displays the chemical complexity that can gradually be increased from prescribed tracer concentrations for the radiation scheme to a more and more comprehensive set of chemical reactions. The chemistry can be used diagnostically only or in a coupled manner (yellow box).

Aims and structure of paper

The aim of this paper is three-fold. Firstly, it serves the documentation of the model and its performance (Sects. 2 and 3). Secondly, we study the sensitivity of the simulated troposphere–stratosphere dynamics to the model setup (Sect. 4), and thirdly we present application examples that serve to highlight the capabilities with the EMIL implementation (Sect. 5). These three parts are stand-alone sections, and the reader may choose to focus on the section of her/his interest.

In the first part (Sects. 2 and 3), we document the implementation of the dynamical core setup within the MESSy framework and its performance. While the Held–Suarez test simulations with the same dynamical core (ECHAM) were previously performed to study the resolution sensitivity of the model core (Wan et al., 2008), the here-presented implementation is new in that it is part of the MESSy framework. The implementation within MESSy ensures an easily accessible idealized model setup that is consistently integrated in the MESSy model hierarchy, and that enables the use of all tracer utilities, including the utilization of diagnostic chemical tracers. The implementation is achieved by adding a simple submodel for Newtonian cooling and Rayleigh friction that replaces the complex physics (see Sect. 2 and the Supplement for technical details including a user manual). We present standard test cases with forcings given by Held and Suarez (1994) and its stratospheric extension (Polvani and Kushner, 2002) in Sect. 3.

In the second part (Sect. 4), we study the sensitivity of the simulated troposphere–stratosphere dynamics with respect to several modifications to previously used setups. Most importantly, this includes modifications of the equilibrium temperatures in the tropical troposphere and in the winter high latitudes (Sect. 4.2). The latter lead to more realistic temperature profiles in the lower stratosphere. We further test the sensitivity of the simulated dynamics to the generation of large-scale waves by zonally asymmetric heating instead of idealized topography, as suggested recently by Lindgren et al. (2018). We end the section with a discussion on the different states of the tropospheric jet and of the stratospheric polar vortex, as well as their relation, in the suite of different sensitivity experiments (Sect. 4.4).

Finally, in the third part (Sect. 5), we present two application examples of the model: first, we present simulations including a small set of chemical reactions (namely photolysis of chlorofluorocarbons) and demonstrate the potential of the model to study the role of dynamical variability and idealized changes on tracer transport (Sect. 5.1). Secondly, the simulation of an upper tropospheric anticyclone forced by simple, localized constrained heating that resembles the Asian monsoon anticyclone is presented in Sect. 5.2.

2 Model description

The ECHAM/MESSy IdeaLized (EMIL) model setup is based on MESSy version 2.54 (Jöckel et al., 2006, 2010, 2016) and will be available for users in the next release, i.e., version 2.55. In the idealized Held–Suarez-type model setup, all physics (radiation, clouds, convection and surface processes) are switched off and are replaced by the newly implemented submodel (“RELAX”) that relaxes the variables temperature and horizontal winds to given background values. The submodel RELAX is described in the next subsection, and we provide technical details of the model setup (namelist choices, etc.) and implementation in the Supplement.

2.1 The submodel RELAX

The submodel RELAX calculates

1. Newtonian cooling, i.e., temperature relaxation towards a given equilibrium temperature with a given relaxation timescale;

2. Rayleigh friction, i.e., horizontal wind relaxation towards zero with a given damping coefficient; and

3. additional diabatic heating over selected regions.

The three processes are switched on/off via namelist parameters, as described in the Supplement. The submodel is called from the physics routine physc through messy_physc. The full call tree including all subroutines is provided in the Supplement. In the following, the implemented options in the routines are described, with the full equations given in Appendix A.

2.1.1 Newtonian cooling

The temperature tendency calculated by Newtonian cooling is given by $\mathit{\delta }T/\mathit{\delta }t=-\mathit{\kappa }(T-T_{\mathrm{eq}})$, where κ is the inverse relaxation timescale, T the actual temperature calculated by the model, and T_eq the prescribed equilibrium temperature. The inverse relaxation timescale and the equilibrium temperature have to be specified in the model setup, either by setting them to fields imported from an external file or by setting them to values given by pre-implemented functions. Currently, the implemented functions for the inverse relaxation timescale and the equilibrium temperature are firstly those given by HS94 (option “HS”; see Eq. A1), but with the possibility to include hemispheric asymmetry, and secondly those given by PK02 (option “PK”; see Eq. A4) but with the following extension: we include the possibility to vary the transition pressure between tropospheric and stratospheric temperature from summer to winter hemisphere. This latitudinal variation is implemented by using the same weighting function as is used for the transition to the polar vortex equilibrium temperatures (see Eqs. A5 and A6). The transition pressure in the remaining area is held constant at p_Ts=100 hPa, as in the original PK setup. Figure 2 shows an example of the equilibrium temperature with modified winter transition pressure, here for p_Tw=400 hPa and γ=2 K km⁻¹. In Sect. 4.2.2, sensitivity simulations with respect to variations in the transition pressure over winter high latitudes (p_Tw) are presented.

Additionally, we performed simulations with a modified prescribed tropospheric equilibrium temperature, differing in the strength of the tropical tropospheric vertical temperature gradient. In this setup of T_eq, in the formulation for the tropospheric equilibrium temperature, the term that reduces the vertical temperature gradient in the tropics (see Eq. A1, fourth term) was inadvertently implemented with a logarithm with base 10 instead of the natural logarithm. As a result, the control parameter of the tropical vertical temperature gradient δ_z=10 K is reduced by a factor of $\mathrm{1}/\log \left(\mathrm{10}\right)\approx \mathrm{0.43}$, so that in the sensitivity experiments ${\mathit{\delta }}_z=\mathrm{0.43}\times \mathrm{10}\mathrm{K}=\mathrm{4.3}\mathrm{K}$. The resulting difference in the equilibrium temperature, as displayed in Fig. 2b, maximizes at around 5.5 K in the tropical upper troposphere. The simulated temperatures in the same region are about 3.5 K lower. While the set of simulations with this modified setup of T_eq was produced inadvertently by an implementation oversight, they provide an interesting sensitivity to the original PK setup and will thus be used in this study to test the sensitivity of the tropospheric jet response to forced polar vortex changes (see Sect. 4.2).

Figure 2 (a) Equilibrium temperature (in K) for p_Tw=400 hPa and γ=2 K km⁻¹ together with the transition pressure p_T(ϕ) (dashed blue line) and the pressure above which damping sets in (dashed–dotted gray line). (b) T_eq as given by the original PK02 implementation (black contours) and difference in T_eq between simulations with reduced δ_z=4.3 K and with standard δ_z=10 K (colors).

2.1.2 Rayleigh friction

Horizontal winds are relaxed to zero (i.e., damped) with a given damping coefficient k_damp by $\mathit{\delta }\mathit{v}/\mathit{\delta }t=-k_{\mathrm{damp}}\mathit{v}$. As for the Newtonian cooling, the damping coefficient can be selected via the namelist with the same options. The implemented functions that can be chosen are

1. the surface layer damping as specified by HS94 (option HS; see Eq. A7),

2. the damping of a layer at the model top as specified by PK02 (option PK; see Eq. A8), and

3. a newly introduced option for damping of a layer at the model top that follows the function as implemented in the original ECHAM code (option EH; see Eq. A9).

For the EH option, the drag is enhanced by a given factor for each level going upward. Sensitivity simulations with respect to the newly implemented form of the upper level damping are presented in Sect. 4.1. Note that as damping at the model surface (option 1) and at the upper layers (options 2 or 3) are complementary; more than one option can be chosen, in which case the profiles of the damping coefficients are added.

2.1.3 Diabatic heating routines

Next to the zonally symmetric temperature tendency calculated by Newtonian cooling, additional temperature tendencies (diabatic heating and cooling) can be added to mimic particular thermodynamic forcings in the atmosphere. Currently, three options are implemented:

1. The first is a function for zonal-mean heating (tteh_cc_tropics) that allows the user to apply a temperature tendency with a Gaussian shape in latitude and pressure, as detailed in Eq. (A10). This zonal-mean heating can be applied to mimic greenhouse gas induced temperature changes, as has been done by, e.g., Butler et al. (2010).

2. The second option is wave-like heating varying with longitude (tteh_waves) that can be used for diabatic generation of planetary waves, as introduced by Lindgren et al. (2018) and detailed in Eq. (A11).

3. The third option is a function for localized heating (tteh_mons) that can be applied to simulate monsoon-like circulation systems, as previously done by, e.g., Siu and Bowman (2019). The formulation for the localized heating is given by Eqs. (A12)–(A16).

3 Model benchmark tests

In this section, results obtained with the EMIL setup are compared to results of earlier studies with identical setups (both with the Held–Suarez setup, Sect. 3.1, and the Polvani–Kushner setup, Sect. 3.2) to test whether the EMIL implementation is able to reproduce the results of those earlier studies.

The benchmark simulations presented in this section are performed for 1825 d for the Held–Suarez forcing, and for 10 950 d for the Polvani–Kushner forcing (see figure captions and Table B1 for details).

3.1 Held–Suarez forcing

Held–Suarez test simulations with the same dynamical core (ECHAM) as the one used in EMIL were previously performed (Wan et al., 2008). We ran a simulation with identical setup and resolution to test whether our implementation of the Held–Suarez setup with the same base model can reproduce the results of Wan et al. (2008). The resolution consists of a spectral horizontal resolution of T63 and 19 vertical levels extending up to 10 hPa, denoted as T63L19. As shown in Fig. 3, the climatologies of zonal wind, temperature and eddy fluxes are closely reproduced when compared to Fig. 1 of Wan et al. (2008). In both model setups, the wind jet maxima are around 30 m s⁻¹, the eddy heat flux maxima around 20 K m s⁻¹ and the eddy momentum flux maxima around 70 m² s⁻². Furthermore, the EMIL temperature and wind climatologies also compare well to the simulations shown in the original HS94 study (see their Figs. 1 and 2).

Figure 3 (a, b, c) Results from a HS simulation at T63L19 resolution, showing mean temperature (K) and zonal-mean zonal wind (m s⁻¹; contour interval is 5 m s⁻¹) (a, d), mean eddy momentum fluxes (m² s⁻²) (b, e) and mean eddy heat fluxes (K m s⁻¹) (c, f) averaged over 1500 d (after spin-up of 325 d). (d, e, f) As above but difference of a simulation at T42L90MA resolution minus the T63L19 simulation (with a wind contour interval of 2 m s⁻¹).

In the remainder of the paper, a vertical resolution with higher top (0.01 hPa) and with 90 levels (L90MA, where MA stands for “middle atmosphere”) will be used together with T42 as spectral resolution (one of the standard resolutions of EMAC; see Jöckel et al., 2016). The differences in the climatologies between the T42L90MA and the T63L19 simulation (for the HS setup) are shown in Fig. 3d–f. The jets are shifted equatorward with T42L90MA resolution, and the eddy variance is generally reduced. This is likely a combined effect of lower horizontal and higher vertical resolution, consistent with the results by Wan et al. (2008): they reported an equatorward shift of the jet and reduced eddy variance both with decreasing horizontal resolution, as well as with increasing vertical resolution. However, in our experiments, not only does the vertical resolution increase but also the model top in the L90MA experiment, and the latter might play a role here. The issue of resolution sensitivity will not be touched further as it is not the subject of this paper, but it should be kept in mind that the results might in general be dependent on the chosen resolution.

3.2 Polvani–Kushner setup

In the study by PK02, an equilibrium temperature is introduced that enables the simulation of an active stratosphere with a polar vortex in the winter hemisphere.

As a test case, EMIL simulations are performed with forcing identical to that in PK02, i.e., with the same choice of the prescribed equilibrium temperature and the damping layer at the top of the model. The results for simulations with the polar temperature lapse rate γ set to 4 K km⁻¹ are shown in Fig. 4a. The polar vortex strength maximizes at around 90 m s⁻¹ for γ=4 K km⁻¹, and at 30 m s⁻¹ for γ=1 K km⁻¹ (not shown), similar to the wind maxima shown in PK02. Also the structure of the polar vortex, and the subtropical jets agree well between our simulation and the ones presented by PK02. Based on the same model as the one used by PK02 (namely GFDL's spectral dynamical core), Jucker et al. (2013) show climatologies of wind and temperature for the PK02 setup with γ=4 K km⁻¹. The temperature climatology of the EMIL simulation with γ=4 K km⁻¹ agrees well with the one shown by Jucker et al. (2013), with both models simulating a tropical lower stratospheric temperature minimum of 210 K and a pronounced minimum in temperature (T<180 K) at the winter pole around 10 hPa.

For a second test case, we include the generation of planetary waves by an idealized topography, as proposed by Gerber and Polvani (2009). Figure 4b shows the simulated climatologies with a wavenumber 2 (WN2) mountain with amplitude h=3 km and γ=4 K km⁻¹. Following Gerber and Polvani (2009), the mountain is centered at 45^∘ N and falls off to zero at 25 and 65^∘ N (see their Eq. 1). This setup of the mountain was found to lead to most realistic simulation of the mean state of the polar vortex and its variability by Gerber and Polvani (2009). The resulting climatologies of zonal wind, with a polar vortex strength of about 50 m s⁻¹, and of temperature, with a minimum temperature over the winter pole at 10 hPa of around 180 K, again closely reproduce the results by Gerber and Polvani (2009) and the equivalent simulation shown by Jucker et al. (2013).

The variability of the polar vortex is diagnosed by the time series of the zonal-mean zonal wind at 10 hPa and 60^∘ N for the PK simulation with γ=4 K km⁻¹ and a WN2 mountain amplitude of h=3 km in Fig. 5 (black line). The polar vortex is highly variable with winds between −10 and 60 m s⁻¹, with sudden decreases in the wind speeds, known as sudden stratospheric warmings. The time series of the EMIL simulation presented here closely resembles that shown by Gerber and Polvani (2009) in terms of variability.

In the study by PK02, it was shown that an increased polar vortex strength, forced by an enhanced stratospheric meridional temperature gradient (i.e., via parameter γ), induces a poleward shift of the tropospheric jet. The tropospheric jet response to stratospheric polar vortex changes will further be discussed in Sect. 4, but it is noted here that EMIL model simulations with the same setup as in the PK02 study reproduce the behavior of the poleward shift of the tropospheric jet with increasing polar vortex strength (see Fig. 8, left, solid yellow line).

Overall, the results of this section show that the EMIL setup is able to reproduce earlier results of simulations performed with dynamical core models under an identical setup of equilibrium temperature, relaxation time, damping layer and topographically generated planetary waves.

Figure 4 Climatologies of wind (black contours, contour interval is 10 m s⁻¹; solid indicates positive and dashed negative) and temperature (K) (colored contours) of (a) an EMIL simulation with the PK02 setup with γ=4 K km⁻¹ and (b) an EMIL simulation with the PK02 setup with γ=4 K km⁻¹ and with WN2 topography with h=3 km. The gray dashed horizontal lines in the EMIL climatologies mark the lower boundary of the damping layer. Averages are performed over 10 000 d.

Figure 5 Time series of zonal-mean zonal wind (m s⁻¹) at 10 hPa and 60^∘ N for different configurations of the PK setup with a WN2 topography with h=3 km. The black line displays the reference simulation with γ=4 K km⁻¹ as in Fig. 4b. The colored lines display the sensitivity simulations discussed in Sect. 4.2, with lowered winter transition pressure p_Tw and for γ=1, 2 and 3 K km⁻¹ (see legend). In the legend, the average wind speed at this location over the whole simulation is given (denoted μ).

4 Sensitivity of coupled troposphere–stratosphere dynamics to modified setups

In this section, the response of the simulated troposphere–stratosphere dynamics to three different types of modifications are studied: (1) modifications of the shape of the upper atmospheric sponge layer (see Sect. 4.1), (2) modifications of the equilibrium temperature, namely of the tropical tropospheric vertical temperature gradient (Sect. 4.2.1) and of the winter high-latitude equilibrium temperature profile (Sect. 4.2.2), and (3) planetary wave generation by wave-like heating instead of topography (Sect. 4.3). The section concludes with a discussion of the sensitivities of the stratospheric polar vortex and the tropospheric jet to the different kinds of modifications.

The simulations presented here are performed for at least 1825 d, and a number of simulations are extended up to 10 950 d. The simulation length is specified for each simulation in Table B1 and in the figure captions. To reduce the uncertainty of the results, it would be favorable to extend each simulation until convergence of the climatologies is reached. In particular, for climate states with multiple dynamical regimes, this would, however, require very long integration times. To reduce computational and data storage costs, we used the strategy of variable simulation length; i.e., we extended only a chosen set of simulations to test for the robustness of the results. In the shorter simulations, considerable uncertainty in the climatologies due to variability can be present. However, as shown in the following, the results are qualitatively robust when comparing the short and long simulations. Details on the simulation setup and integration length can be found in Table B1.

4.1 Sensitivity to the shape of the upper atmospheric damping layer

The damping layer at levels above 0.5 hPa is included to account for the strong damping of winds that in the real atmosphere (or the full model) is due to drag by breaking gravity waves (GW). The simplified manner of damping the entire horizontal wind fields introduces a non-physical sink of momentum, as not only the zonal-mean wind but also all waves are damped. When analyzing results obtained with the model, this has to be kept in mind.

The damping layer as introduced by PK02 uses a damping coefficient that increases quadratically with decreasing pressure. The profile of the PK02 damping coefficient is shown in Fig. 6 together with the profile of zonal-mean zonal wind tendencies due to parameterized gravity waves divided by the zonal-mean wind (averaged over 40–60^∘ N/S) from a model simulation with the full atmospheric EMAC setup, i.e., an equivalent damping coefficient of the zonal-mean wind by the parameterized GW drag. The “damping” by GW drag varies between years and hemispheres but generally increases exponentially with decreasing pressure, not quadratically. Therefore, we argue that a damping coefficient with exponential increase mimics the net effects of parameterized GW drag better.

A sensitivity simulation is performed in which the damping coefficient in the upper model domain follows the exponential function given by Eq. (A9) (option EH; this is the shape of the “sponge” layer originally implemented in the ECHAM model). The damping coefficient of this sensitivity simulation is shown in Fig. 6 as a red line.

Figure 6 Damping coefficient (s⁻¹) of the sponge layer in the EH (red) and in the PK02 setup (blue) together with the effective damping timescale of the zonal winds by gravity wave drag (i.e., −GWD/u) from an ECHAM simulation with a non-orographic and orographic GW scheme, averaged over 40–60^∘ S from July 1960 to 2010 (dashed gray line; average over all years shown as a dashed black line) and averaged over 40–60^∘ N from January 1960 to 2010 (solid gray line; average over all years shown as a solid black line).

Figure 7 (a) Differences in zonal-mean temperature (contour interval is 2 K) and zonal-mean zonal wind (contour interval is 2 m s⁻¹) between a model simulation with exponentially increasing damping coefficient (EH) and a model simulation with quadratically increasing damping coefficient (PK) for p_Tw=100 hPa and γ=4 K km⁻¹ without topography. Averages are performed over 1500 d of simulation. Stippling indicates non-significant wind differences (on a 95 % level, based on a t test performed on 30 d means). (b) Polar winter temperature profiles of same model simulations averaged from 70^∘ N to 90^∘ N, with temperature profiles from the SPARC climatology in northern hemispheric winter conditions (solid black line) and southern hemispheric winter conditions (dashed black line) as well as the equilibrium temperature profile (dashed gray line). The dashed–dotted line marks the lower boundary of the sponge layer.

The simulated climate states with the two different setups naturally differ within the sponge layer, with maximum differences in zonal winds of 30 m s⁻¹ around 0.5 hPa (see Fig. 7). Considerable differences in wind and temperature also extend below the damping layer, in particular at high latitudes but also in the tropics. In the tropics, alternating jets form, reminiscent of the quasi-biennial oscillation but with very long periods (about 5 years). As the tropical winds have thus not converged for the given simulation length, and differences are not significant, we will not regard them here.

Differences are mostly insignificant below 10 hPa; however, small (significant) differences in zonal winds of 2 m s⁻¹ extend down into the troposphere. The increase in high-latitude zonal winds, which maximizes at the lower bound of the damping layer, is accompanied by an upward shift of the temperature maximum at the stratopause. This brings the temperature maximum closer to realistic values, as evident from the comparison to the “SPARC” climatology (Randel et al., 2004; SPARC, 2002) in Fig. 7b.

Since the EH sponge is weaker, the increase in zonal-mean winds within the damping layer can be expected. The weaker sponge and changed zonal wind structure modifies planetary wave propagation (stronger upward propagation between from about 3 hPa upward; not shown), thus influencing the mean climate also below the damping layer (decreased wave driving, leading to stronger zonal winds and lower polar temperatures). The effect of the modified damping coefficients is similar in simulations with WN2 topography (albeit with weaker absolute differences; not shown).

As the exponentially increasing damping coefficient (EH) resembles the vertical structure of GW drag, and since for both a flat surface and idealized topography, the height at which the polar winter temperature profile reaches its maximum is more realistic in the case of the EH damping layer (see Fig. 7b), we chose to use the exponentially increasing damping coefficient (EH) in the following as our reference setup.

4.2 Sensitivity to modification of the equilibrium temperature

In addition to the benchmark simulations with identical setup as in PK02, we performed sensitivity simulations with modified prescribed equilibrium temperature, with an altered tropospheric vertical temperature gradient (Sect. 4.2.1) and with reduced winter high-latitude lower stratospheric temperatures (Sect. 4.2.2).

4.2.1 Sensitivity to reduced tropical tropospheric vertical temperature gradient

In the standard HS94 setup, the vertical temperature gradient in the (sub)tropical troposphere is reduced to mimic the effects of latent heat release. This reduction is controlled by the parameter δ_z in the fourth term of Eq. (A1), and in the sensitivity simulations presented here, δ_z was reduced from 10 to 4.3 K, as described in Sect. 2 (see also Fig. 2b for the resulting difference in T_eq). These simulations with reduced static stability in the (sub)tropics resulted from an implementation oversight but proved to provide an interesting sensitivity to the original PK02 setup.

EMIL model simulations with an identical setup to that in PK02 reproduce the strong poleward shift of the tropospheric jet in response to a forced increase in the polar vortex strength via the prescribed stratospheric meridional temperature gradient (i.e., via parameter γ). The increase in polar vortex strength (measured by the zonal-mean wind speed at 60^∘ N and 10 hPa) with increasing values of the polar vertical lapse rate γ is shown in Fig. 8 (bottom left), together with the location of the tropospheric jet (Fig. 8, second row from the bottom), both for the simulations with original PK02 setup and the setup with the altered tropical tropospheric temperature gradient (dashed lines). The tropospheric jet location is measured here as the latitude of the maximum of climatological mean zonal-mean zonal winds at 500 hPa. Note that this value can differ from the time mean of the jet latitude determined at individual days (compare to mean values given in the Appendix, Figs. C1 to C3).

For increasing polar lapse rates γ, the polar vortex strengthens at a similar rate in the original setup and the simulations with reduced δ_z (this also applies to variability, as evident from probability distributions; see Fig. C1). Consistently, also the residual circulation changes in a similar manner in the two sets of simulations. This is diagnosed here as deviation of temperature from the equilibrium temperature, a valid measure of the strength of the residual meridional circulation in the idealized model (see, e.g., Jucker et al., 2013). We choose to average this temperature difference from 40 to 90^∘ N, as this is the region of diabatic heating associated with downwelling. Larger values of this temperature difference therefore imply a stronger circulation. In Fig. 8 (top two rows), these temperature differences are displayed for 10 and 100 hPa to represent the strength of the circulation in the middle and lower stratosphere, respectively.

Despite the similar response of stratospheric dynamics and polar vortex strength in the simulations with reduced δ_z compared to the simulations with original PK setup, the response of the tropospheric jet is strongly damped compared to the original PK02 setup. At 500 hPa, the maximum wind location shifts only by a few degrees from the simulation without a polar vortex to the one with γ=5 K km⁻¹, while in the original PK02 setup this shift amounts to more than 10^∘ latitude (see Fig. 8). Note that in the reduced δ_z simulation, the tropospheric jet is shifted slightly equatorward compared to the reference simulations already for the basic state without polar vortex (see Fig. 8), which might be the reason for the different response, as will be discussed in Sect. 4.4.

As has been shown by Gerber and Polvani (2009), the response of the tropospheric jet location to stratospheric forcing is strongly damped in simulations with idealized topography compared to those without topography. The EMIL model simulations presented here reproduce the damped response to stratospheric forcing (changes in γ) under same setup as in Gerber and Polvani (2009), as shown in Fig. 8 (middle, yellow line). In the simulations including topography, the tropospheric jet is likewise shifted equatorward with reduced δ_z. However, we detect little difference in the tropospheric jet response to stratospheric polar vortex changes between the two sets of simulations with topography. Thus, in the presence of planetary wave forcing, the different tropospheric equilibrium temperatures appear to play a smaller role for the stratosphere–troposphere coupling. This will further be discussed in Sect. 4.4.

Figure 8 Top row: difference between temperature and equilibrium temperature T−T_eq averaged from 40 to 90^∘ N at 10 hPa. Second row: same but at 100 hPa. Third row: latitude ϕ_max of the zonal-mean zonal wind speed maximum of the tropospheric jet. Fourth row: zonal-mean zonal wind u at 60^∘ N and 10 hPa. The left column shows results from model simulations without planetary wave forcing, the middle column with WN2 topography of height h=3 km and the right column with WN2 tropospheric heating of amplitude q₀=6 K d⁻¹. Solid lines represent simulations with standard δ_z=10 K, dashed lines with reduced δ_z=4.3 K.

4.2.2 Sensitivity to modification of the equilibrium temperature in the winter high-latitude lower stratosphere

The simulated winter high-latitude temperature profiles for EMIL simulations with PK setup and WN2 topography are shown in Fig. 9a for varying γ, compared to temperature profiles from the observationally based SPARC climatology for northern winter. The comparison of the simulations to the observational climatology reveals a positive temperature bias in the upper troposphere–lower stratosphere (UTLS) region of the winter high latitudes (70 to 90^∘ N), when using the standard PK setup with a constant transition pressure of p_T(ϕ)≡100 hPa. The positive temperature bias remains unchanged for different polar vortex lapse rates γ. Even for strong decreases of the equilibrium temperature above the 100 hPa level, the positive temperature bias in the UTLS region cannot be compensated. This is essentially because the equilibrium temperature already exceeds the observational temperatures in that region. Due to the general-circulation transport of heat from the tropics to polar regions throughout the troposphere and stratosphere, the temperature bias even increases. Therefore, every simulation with p_Tw=100 hPa necessarily has a too-warm UTLS region in the winter high latitudes compared to observations. The warm bias is associated with an unrealistic “step” in the temperature profile, forced by the constant equilibrium temperature profile in the UTLS up to 100 hPa.

Figure 9 Polar winter temperature profiles of model simulations with WN2 topography of height h=3 km and different polar vortex lapse rates γ for p_Tw=100 hPa (a) and p_Tw=400 hPa (b), together with the temperature profiles obtained from the SPARC climatology (black line) as well as the equilibrium temperature profiles (dashed colored lines). Averages are based on about 10 000 d.

In order to approach a more realistic temperature profile in the UTLS region of the winter high latitudes, the transition pressure p_Tw is increased. A similar approach was used by Sheshadri et al. (2015), who lowered the transition pressure globally to 200 hPa and showed that this led to an improvement in lower stratospheric zonal winds. Here, we systematically vary the the transition pressure in polar winter high latitudes across p_Tw=100 to 450 hPa as well as γ across γ=1 to 4 K km⁻¹ (see Table B1).

As discussed in the previous section, a modified version of the tropospheric equilibrium temperature function with a changed vertical temperature gradient was implemented in a previous model version. The whole parameter sweep was performed in this modified model setup with reduced δ_z, and we repeated simulations for p_Tw=100 and 400 hPa with the standard setup to test the combined sensitivity of the results on modifications in p_Tw and on the tropical tropospheric equilibrium temperature. In Fig. 8, the results of both setups are shown (with the modified simulations as dashed lines). The simulations with reduced δ_z are performed for 1825 d (with the first 300 d considered as spin-up), while the simulations under the standard setup were extended to 10 950 d (with 1000 d of spin-up). While the 1825 d simulations are mostly too short to establish convergence of the climatologies, the qualitative behavior in those simulations is in general in agreement with the results from the extended simulations (as presented in the following).

The polar temperature profiles shown in Fig. 9 for the simulations with standard setup are very similar to those for the simulations with the modified setup (not shown). In both setups, at p_Tw=400 hPa, all equilibrium temperatures in the polar winter UTLS region fall below the corresponding temperatures obtained from observations except the one for γ=1 K km⁻¹ (see Fig. 9b). For the simulations with γ=3 K km⁻¹ and γ=4 K km⁻¹, the winter high-latitude temperatures are lower than observational temperatures throughout the UTLS region and follow the equilibrium temperature up to about 30 hPa. Above, the temperature increases strongly, reaching a maximum at around 0.7 hPa. In contrast, the UTLS temperature in the simulation with γ=1 K km⁻¹ is well above the corresponding equilibrium temperature in the UTLS, and the temperature maximum at around 0.5 hPa is weaker. The simulation with γ=2 K km⁻¹ lies in between: its temperature in the UTLS is higher than the equilibrium temperature but less so than for γ=1 K km⁻¹.

The non-linear behavior of the deviation from the equilibrium temperature is illustrated for a variety of values of γ and p_Tw in Fig. 8. For low values of γ, we find an increase in T−T_eq with γ in the mid-stratosphere but a decrease in the lower stratosphere, in agreement with the result of Gerber (2012) that a stronger polar vortex leads to a strengthened circulation in the upper stratosphere and to a weakened circulation in the lower stratosphere (see Fig. 3 of Gerber, 2012, for comparison). However, above a certain threshold of γ, the circulation strength decreases with γ and then stagnates also in the mid-stratosphere (10 hPa). This critical value of γ depends on p_Tw, in line with lower polar T_eq values for higher values of p_Tw. In the upper stratosphere (1 hPa), a monotonic increase in T−T_eq both with larger γ and p_Tw is found (not shown), but we exclude this analysis because of the likely influence of the upper damping layer.

The strength of the polar vortex increases for larger γ and p_Tw values, as expected from the stronger meridional temperature gradient (see Fig. 8). However, the polar vortex increases non-linearly with increasing γ, with stronger acceleration above a critical value. This is in line with the change in behavior of T−T_eq at 10 hPa when reaching this critical value (e.g., for p_Tw=400 hPa, between γ=2 K km⁻¹ to 2.5 K km⁻¹ in the modified setup). Thus, for increases in the prescribed meridional temperature gradient in the polar stratosphere (i.e., via γ) below a certain threshold, the polar vortex strength increases only very weakly. At the same time, midlatitude to high-latitude temperatures increase above the corresponding equilibrium temperature (i.e., T−T_eq increases with γ). Thus, the residual circulation is strengthened, and the associated high-latitude warming counteracts the increase in the prescribed meridional temperature gradient, explaining the weak changes of the polar vortex. Once a certain threshold in the prescribed meridional equilibrium temperature is reached, the polar vortex increases strongly, and at the same time T−T_eq decreases, indicating a reduction in wave driving and thus additional dynamical strengthening of the meridional temperature gradient.

In response to the polar vortex increase for larger γ and p_Tw values, the tropospheric jet shifts poleward (see Fig. 8). This poleward shift is found to be similarly strong in the standard setup and the modified (reduced δ_z) setup. Thus, the strong dependence of the strength of the tropospheric jet response on the setup found for the simulations without topography (see Sect. 4.2.1 and Fig. 8, left) is not present in the simulations with WN2 forcing, even under the stronger stratospheric forcing in the simulations with increased p_Tw. The sensitivity of the rate of the tropospheric jet response will further be discussed in Sect. 4.4.

Overall, the stratospheric circulation responds non-linearly to modifications of the winter equilibrium temperature profile. Lowering the height at which the equilibrium temperature starts to decrease can diminish the high-latitude lower stratospheric temperature bias. To more or less completely remove the warm bias and the associated unrealistic “step”, p_Tw has to be increased to 400 hPa. In the simulation setup with p_Tw=400 hPa and γ=2 K km⁻¹, the winter high-latitude temperature profile is close to reanalysis data (SPARC climatology and ERA-Interim; the latter not shown) in the UTLS region and a moderate oscillation of the temperature in the upper atmosphere is simulated.

4.3 Planetary wave generation with topography versus heating

In the experiments presented in the preceding Sect. 4.2, an idealized topography was used to generate planetary waves. Recently, Lindgren et al. (2018) suggested an alternative method to generate planetary waves in a setup with an active stratosphere: they introduced a tropospheric wave-like thermal forcing of the form of Eq. (A11), which is added to the temperature tendency of Newtonian cooling.

For the equilibrium temperature, Lindgren et al. (2018) employ a constant transition pressure of p_T(ϕ)=200 hPa, i.e., $p_{\mathrm{Ts}}=p_{\mathrm{Tw}}=\mathrm{200}$ hPa, and ϵ=0, i.e., a hemispherically symmetric temperature distribution in the troposphere. Figure 10 shows the temperature profiles in the winter high latitudes for different simulations that were thermally forced by Eq. (A11). The model simulation with the original Lindgren setup exhibits a too-high temperature in the winter high-latitude UTLS region compared to the SPARC climatology for the same reason as was explained for the topographically forced simulations with p_T(ϕ)=100 hPa (PK setup) in the previous subsection: the decrease of the equilibrium temperature due to γ starts too high to be able to compensate the warm bias. This motivated the investigation of model simulations with a larger transition pressure p_Tw in the winter high latitudes for the thermally forced simulations as well.

In our model simulations with WN2 tropospheric heating, we similarly use ϵ=0 but return to p_Ts=100 hPa and vary p_Tw.¹ In addition to the profile obtained from the Lindgren setup, Fig. 10 contains the winter high-latitude temperature profiles for different polar vortex lapse rates, γ and for p_Tw=400 hPa.

Figure 8 (right) shows results from the simulations with thermally forced planetary waves, and again both the modified (reduced δ_z) simulations as well as simulations with the standard tropospheric equilibrium temperatures are included. As discussed for the topographically forced simulations, we also find a non-linear behavior of the stratospheric circulation in the thermally forced simulations: a state with a weak polar vortex (for γ smaller than about 3 K km⁻¹; see Fig. 8, bottom right) manifests in a positive temperature bias in the UTLS region of the winter high latitudes (Fig. 10), and increasing temperature deviation from the equilibrium temperature with increasing γ (Fig. 8, top right). A state with a strong polar vortex arises for γ≥3 to 3.5 K km⁻¹ (Fig. 8, bottom right), with temperature following the equilibrium temperature very closely in the UTLS region (see Fig. 8).

Figure 10 Same as Fig. 9 but for model simulations with WN2 tropospheric heating of amplitude q₀=6 K d⁻¹ and different polar vortex lapse rates γ for p_Tw=400 hPa. The temperature profiles obtained from a simulation with the original Lindgren setup (dotted black line) are added for comparison.

The response of the polar vortex to changes in the equilibrium temperature is similar between the topographically versus thermally forced model simulations in that a transition from a weak to a strong polar vortex is found for both cases. Thermally forced model simulations also show an increase of the strength of the meridional circulation at 10 hPa up to a certain threshold of γ, similar to the topographically forced simulations (see Fig. 8, top). Note, however, that the threshold is higher for the thermally forced simulations for identical equilibrium temperature. The change in the behavior of the meridional circulation in the model simulations with p_Tw=450 hPa and p_Tw=400 hPa (both for the modified and standard setup) appears at the same polar vortex lapse rates, at which the polar vortex starts to strengthen. At 100 hPa, the topographically forced simulations show a (non-linear) decrease of the circulation strength with increasing γ for all values of p_Tw, while in the thermally forced simulations the circulation in the lower stratosphere responds in a similar non-linear way as at 10 hPa.

Further, we compare the response of the tropospheric jet to changed equilibrium temperatures in topographically forced simulations to the response in the thermally forced simulations. As discussed in the last section, in the case of the topographically forced model simulations, the location of the free tropospheric jet shifts poleward in the simulations with a stronger stratospheric polar vortex. However, when the planetary waves are thermally forced, the free tropospheric jet maximum remains at an almost constant latitudinal location or even moves equatorward (see lower panels of Fig. 8). Even strong increases in the stratospheric polar vortex for γ>3 K km⁻¹ at p_Tw=400 hPa and for γ>2.5 K km⁻¹ at p_Tw=450 hPa, respectively, are not accompanied by a northward shift of the free tropospheric jet maximum.

Thus, while the non-linear increase in the stratospheric polar vortex strength is overall similar in the topographically and diabatically forced simulations, other aspects of the circulation response show distinct differences. Overall, the different behavior of model simulations with topographically and thermally forced circulations outlined here indicates that the thermally forced simulations might have to be used with caution, in particular for studying troposphere–stratosphere coupling.

4.4 Discussion on dynamical states of the stratospheric polar vortex and the tropospheric jet

The sensitivity simulations with respect to modifications of the equilibrium temperature (Sect. 4.2) and to planetary wave generation (Sect. 4.3) revealed the following. Firstly, the stratospheric polar vortex responds non-linearly to the enhancement of the meridional temperature gradient in the simulations with planetary wave forcing. Secondly, the strength of the tropospheric jet response to this stratospheric vortex strengthening depends on the model's basic state. These two results are discussed in the following.

4.4.1 Stratospheric polar vortex regimes

We found strong non-linear strengthening of the polar vortex with enhanced meridional temperature gradients (via increasing γ and/or p_Tw), i.e., the polar vortex transitions from a weak state to a strong state. In between, climate states are found in which the vortex appears to alternate between those two states, indicative of a transition between different vortex regimes. This is reflected in changes in the polar vortex variability, as shown in Fig. 5 for the topographically forced simulations with p_Tw=400 hPa and γ=1, 2, and 3 K km⁻¹. While the simulation with the weakest polar vortex (γ=1 K km⁻¹) exhibits large variability with frequent crossings of the zero-wind line (indicative of sudden stratospheric warmings), in the simulation with γ=3 K km⁻¹, the wind oscillates around its large mean value without crossing the zero-wind line. With an intermediate polar lapse rate (γ=2 K km⁻¹), episodes with strong stable winds are disrupted by sudden decelerations, and the polar vortex remains in a weak state for an extended period (up to a several hundred days) thereafter; i.e., the vortex alternates between a strong and a weak regime. The regime behavior is further supported by the shape of the probability distribution of the maximum wind at 10 hPa for those three simulations (see Fig. C2): while in the simulation with γ=1 K km⁻¹, the polar vortex strength is bound below 50 m s⁻¹, and in the simulation with γ=3 K km⁻¹, the vortex strength is always above 75 m s⁻¹, for γ=2 K km⁻¹, a broad, nearly bimodal distribution is found. The distribution functions for the modified (reduced δ_z) simulations (not shown) are noisier due to the shorter simulation length but show a similar behavior to those shown in Fig. C2. Further, also the diabatically forced simulations indicate a regime shift of the polar vortex: the probability distribution functions of the polar vortex strength change strongly between the weak and strong vortex state (e.g., change in sign of skewness; see Fig. C3), supporting further that we see a regime transition. In terms of polar vortex changes with increasing prescribed polar stratospheric temperature gradient, the simulations with a different setup show only minor differences (both for modifications in δ_z and in planetary wave generation).

Analysis of the residual circulation strength (as measured by deviation from equilibrium temperature; see Sect. 4.2) indicates that the sudden strengthening of the polar vortex is associated with a strong decrease of the residual circulation in the mid-stratosphere. In the weak vortex regime, the weak westerly winds allow for vertical planetary wave propagation, as expected from linear wave theory (Charney and Drazin, 1961). The waves dissipate in the vicinity of the polar vortex and drive the residual circulation that in turn reduces the meridional temperature gradient. When increasing the prescribed meridional equilibrium temperature gradient up to a certain threshold, the wave forcing appears to be increasing. Thus, the stronger residual circulation counteracts the decreasing equilibrium temperatures in the polar stratosphere, and the polar vortex changes little. However, when lowering polar equilibrium temperatures to a certain threshold, the wave driving seems not to be able to counteract the strengthening meridional temperature gradient any longer: the polar vortex strengthens, and thereby vertical wave propagation is inhibited by the strong winds (again by simple arguments following Charney and Drazin, 1961). The positive feedback between strong winds and suppressed wave forcing might explain the suddenness of the polar vortex strengthening, or in other words, it explains why the vortex behaves regime-like. The non-linear interaction of planetary waves and the mean flow is less pronounced in simulations without planetary wave forcing, explaining that the polar vortex responds more linearly to the prescribed temperature changes. However, the above line of argument will have to be tested by more thorough analysis of wave fluxes in the simulations.

4.4.2 Tropospheric jet location and its response to stratospheric forcing

The strength of the tropospheric jet response to the stratospheric vortex strengthening is summarized in Fig. 11, where the tropospheric jet location is shown as a function of the polar vortex strength for increasing polar lapse rate γ for a number of different experiment setups. The slope of the lines in Fig. 11 thus indicates the sensitivity of the tropospheric jet location to stratospheric polar vortex changes. The strongest response is found in the experiments with the setup as in the original PK02 study (black line), while the response of the tropospheric jet is considerably weaker in the setup with reduced (sub)tropical static stability (reduced δ_z, dashed black line). This holds for locating the tropospheric jet in the mid-troposphere (Fig. 11a) as well as near the surface (Fig. 11b).

The inclusion of topographically generated planetary waves dampens the response of the tropospheric jet (as previously shown by Gerber and Polvani, 2009), both in the original and the reduced δ_z setup (see blue lines). With reduced polar lower stratosphere temperatures (i.e., p_Tw=400 hPa), the stronger polar lapse rate leads to enhanced polar vortex strengths, and the tropospheric jet is shifted northward with an intermediate, almost linear rate (red lines). Again, the jets diagnosed near the surface and in the mid-troposphere reveal similar rates of change (compare Fig. 11a and b). Under same setup of T_eq but with diabatically forced planetary waves, the tropospheric jet diagnosed in the mid-troposphere is almost insensitive to the stratospheric polar vortex increase. However, when diagnosed near the surface, the jet shifts poleward at an even higher rate than in the topographically forced simulations (Fig. 11, purple lines).

As mentioned in the introduction, the strong poleward jet displacement in the original PK02 experiments has been shown to be associated with a shift between subtropical and midlatitude jet regimes (e.g., Chan and Plumb, 2009). While we do not find the bimodal distribution of the near-surface jet location as shown by Chan and Plumb (2009), a broadening of the probability distribution of the tropospheric jet location in the simulation with γ=3 K km⁻¹ and a change in skewness of the distributions (from positive to negative between γ=3 K km⁻¹ and γ=4 K km⁻¹) indicate a regime shift in the jet location also in our simulations. The probability distributions are appended for reference; see Fig. C1.

In the simulations without topography and with reduced δ_z, the tropospheric jet shifts slightly poleward with increasing γ but remains in the subtropical regime (see Fig. 11 and Fig. C1). We presume that the more equatorward location of the tropospheric jet in the basic state inhibits the regime transition to a poleward-located tropospheric jet in the reduced δ_z simulations in response to the stratospheric forcing. The equatorward shift of the jet in the basic state is consistent with the previously reported jet response to the insertion of diabatic heating in the tropical upper troposphere (e.g., Butler et al., 2010). The anomalies in the equilibrium temperature extend into the subtropics (see Fig. 2), so that the static stability is reduced in the tropical and subtropical troposphere. The reduced static stability effectively enhances the tendency for baroclinic instability in the subtropics (e.g., Lu et al., 2008), which could favor a subtropical eddy-driven jet location. This is consistent with the persistent subtropical jet regime in the simulations with reduced δ_z. Whether the jet would move to the higher-latitude regime in our reduced δ_z simulations under stronger stratospheric forcing remains to be investigated.

In the simulations with topography, the jet is likewise shifted equatorward with reduced δ_z. However, the jet response to the strengthening stratospheric polar vortex is not different in the simulations with reduced δ_z (compare solid and dashed red lines in Fig. 11a). This could be due to the additional effects of planetary waves, again consistent with the result by Lu et al. (2008) that changes in the jet location are more tightly related to the subtropical stability in the Southern Hemisphere (SH) (with little planetary wave activity) than in the Northern Hemisphere (NH).

When the planetary waves are forced diabatically, the basic state mid-tropospheric jet is likewise shifted equatorward in the simulations with reduced δ_z. The response to stratospheric forcing is again unaltered, as in both sets of simulations the location of maximum winds in the mid-troposphere appears to be fixed in the midlatitude regime throughout the range of polar vortex strengths. Consistent with the static stability argument, the diabatically forced simulations exhibit the highest gross static stability in the subtropical troposphere compared to the other sets of simulations (not shown). Thus, the reduced tendency towards baroclinic instability in the subtropics favors the midlatitude jet regime. However, as the equilibrium temperature is identical in the topographically and diabatically forced simulations, there is no obvious reason for the enhanced stability and the different basic states of the tropospheric jet.

Overall, the set of simulations presented here confirms that the tropospheric jet in the dry dynamical core model tends to fall into either a subtropical or a midlatitude regime. This extends the result of Chan and Plumb (2009), in that different kinds of modifications of the setup can lead to states with a strongly damped response of the tropospheric jet location to the stratospheric forcing. In Chan and Plumb (2009), enhancement of the surface equilibrium temperature Equator-to-pole gradient led to states with a weak jet response, because the tropospheric jet is located in the midlatitude regime already for a weak stratospheric polar vortex, similar to our diabatically forced simulations. However, in the diabatically forced simulations, the wind maximum in the mid-troposphere and near the surface seem to be decoupled: while the former remains at a rather constant latitude, the latter strongly moves poleward under strong stratospheric forcing, with signs of a regime transitions to an even higher latitude regime (indicated by bimodality; see Fig. C3). Next to this state with the jet remaining in the midlatitude regime, we also found a state in which the tropospheric jet remains in the subtropical regime even under strong forcing (namely in the setup with reduced δ_z and no planetary wave forcing). However, the jet might likely move to the midlatitude regime if the stratospheric forcing was increased further, but that remains to be investigated.

Figure 11 Latitude of the zonal-mean zonal wind maximum at 500 hPa (a) and 850 hPa (b) displayed against the maximal zonal-mean zonal wind u_max at 10 hPa for various simulation setups: simulations under the PK setup without topography (labeled WN0), and with WN2 planetary waves forced topographically (labeled WN2T) and diabatically (labeled WN2H) for various values of the winter transition pressure p_Tw and tropospheric tropical vertical temperature gradient δ_z (see legend; if no value is given, the parameter is set to default). Each symbol displays the value diagnosed from the climatology of simulations with varying polar vortex lapse rate γ. The values for the WN2T simulations with δ_z=4.3 K are not shown on the right, because data at 850 hPa were not saved appropriately.

5 Application examples

In the previous sections, the implementation of the EMIL model was documented and tested, and modified setups were introduced, showing that the model is well suited for further applications. In the following, two examples of research applications with the dynamical core model are shown. First, variability and changes in tracer transport in response to changes in the polar vortex are analyzed, using the simulation setup with the modified equilibrium temperature (see Sect. 4.2). Secondly, the localized heating routine (see Sect. 2.1.3) is used to force an idealized monsoon circulation system.

5.1 Chemistry and tracer transport

With the implementation of the idealized model setup in the MESSy framework, all tracer utility and chemistry submodels can be easily used to study the tracer transport in the idealized model. Within the chemistry submodel MECCA (Sander et al., 2019), tailor-made chemical mechanisms can be selected to the users' needs, allowing for selection of simplified chemistry setups. As a proof of concept, we present results from simulations where the only selected chemical reactions are the photolysis of chlorofluorocarbons (CFCs, namely CFC-11 and CFC-12).

Technically, this simulation setup requires, in addition to the “standard” EMIL setup, to switch on submodels for solving chemical kinetics (MECCA; Sander et al., 2019), calculating photolysis rates (JVAL; Sander et al., 2014) and determining orbital parameters (ORBIT; Dietmüller et al., 2016), as well as submodels for tracer definition (TRACER and PTRAC; Jöckel et al., 2008) and tracer boundary condition nudging (TNUDGE; Kerkweg et al., 2006). CFC mixing ratios were set to values representative of the year 2000 at the surface, and tracers were initialized with a mean distribution from an earlier EMAC simulation. To obtain constant January conditions of solar irradiance (compatible with the idealized thermodynamical forcing of the dynamics), in the TIMER namelist, a perpetual month simulation can be selected.

Figure 12 (a) Time series of zonal-mean zonal wind $\stackrel{\mathrm{‾}}{u}$ at 60^∘ N and 10 hPa (black) and mean $\stackrel{\mathrm{‾}}{w^{*}}$ at 50 hPa and 60–90^∘ N (red). (b) Zonal-mean CFC-11 mixing ratios (color, in ppbv) at 50 hPa as a function of simulated day and latitude, and zonal-mean zonal wind at 50 hPa (white contours, interval 15 m s⁻¹). Vertical gray lines mark dates of SSWs (centered here at the dates of the strongest wind decelerations). Right: time-mean CFC-11 mixing ratios as function of latitude over days with a strong vortex (days 400–600; 1200–1350; 2000–3000, solid black line), over days following SSWs with strong downwelling (days 600–780; 1380–1500; 3100–3280, dashed) and over days with eroded polar vortex (days 800–1000 and 1580–1700, dotted).

With the given model setup including chemical tracers, the influence of idealized dynamical variability on chemically active species can be studied. Shown in Fig. 12 are zonal-mean CFC-11 mixing ratios at 50 hPa as function of latitude and time in a simulation with the PK setup, with the reduced value of δ_z=4.3 K and with p_Tw=400 hPa and γ=2 K km⁻¹. The polar vortex variability leads to variability in CFC-11 mixing ratios in particular at high latitudes. As diagnosed from the time series of zonal-mean zonal wind at 60^∘ N and 10 hPa (Fig. 12a, black line) sudden stratospheric warming (SSW) events occur at around day 600 and day 1350, as indicated by the dashed gray lines. While the most common definition of the zero crossing of the 10 hPa zonal wind is met a few days later, the lines are inserted at the dates of the strongest deceleration of the wind. Both events are followed by an extended period with a weak polar vortex.

For both SSW events, the CFC-11 mixing ratios drop at high latitudes simultaneously with the drop of zonal winds at 10 hPa. However, around 200 d after the SSW events, high-latitude mixing ratios increase again. This behavior can be explained as follows. Simultaneously with the SSW, strong downwelling occurs at high latitudes (north of 60^∘ N), driven by the strong wave dissipation that affected the SSW (see red line in Fig. 12a). The enhanced downwelling transports CFC-depleted air from higher altitudes downward. However, due to the diminished vortex in the period after the SSW, air from midlatitudes with higher CFC mixing ratios can be mixed towards the pole, thus leading to an enhancement of CFC mixing ratios at high latitudes. This is evident in Fig. 12 around days 800–1000 and days 1500–1700, when zonal winds are below 15 m s⁻¹.

The transport anomalies are evident in the latitudinal profiles of CFC-11 mixing ratios, as shown on the right side of Fig. 12: during episodes with a strong polar vortex (solid line), there is a local minimum in mixing ratios close to the polar vortex edge (in agreement with strongest downwelling at the vortex edge; see Fig. 13c), denoting the separation between midlatitude and high-latitude air by the vortex. Just during and after the SSW events, CFC mixing ratios drop at high latitudes (dashed line), while in the episodes with eroded vortex, CFC-11 mixing ratios are enhanced at midlatitudes to high latitudes and no mixing barrier can be identified (dotted line).

Two additional simulations were performed with idealized changes in the polar vortex (intermediate vortex: γ=2 K km⁻¹, weak vortex: γ=1 K km⁻¹, strong vortex: γ=3 K km⁻¹). The resulting climatological mean CFC-11 mixing ratios at 50 hPa are shown in Fig. 13a. The differing dynamical states of the three simulations are clearly reflected in the tracer mixing ratios. The simulation with the weak vortex (γ=1 K km⁻¹, red) shows highest CFC-11 mixing ratios in the tropics to midlatitudes, with a smooth transition from tropics to high latitudes, in line with strongest upwelling (see Fig. 13c; see also results in Sect. 4.2) and strong midlatitude wave driving that results in mixing (see Fig. 13d). In the simulation with a strong vortex (γ=3 K km⁻¹, blue), mixing ratios in the tropics are lower, due to weaker upwelling in the lower stratosphere, and the gradient to midlatitudes is steep. This can be explained due to weaker mixing (see Fig. 13d) and also because the region of downwelling is shifted towards lower latitudes (see Fig. 13c). In the simulations with stronger polar vortex (γ=2 K km⁻¹ and 3 K km⁻¹), downwelling is maximized at the equatorward flank of the polar vortex and is weak within the vortex. In contrast, in the simulation with a weak polar vortex (γ=1 K km⁻¹), downwelling is maximized more poleward and is stronger also at high latitudes. The maximum of downwelling in the midlatitudes as well as the high isolation of vortex air in the strong vortex case likely explains why CFC-11 mixing ratios are elevated within the vortex. The intermediate simulation with γ=2 K km⁻¹ lies in between the other two simulations but shows highest variability (largest standard deviation) in most quantities, as expected, since this simulation oscillates between states with a weak and strong vortex (see Fig. 12 and Sect. 4.2).

As demonstrated here, the idealized setup of the simulation allows us to study the role of vortex variability or specifically forced polar vortex strength changes on tracer mixing ratios in an isolated manner, i.e., the absence of other chemical processes or variability like the annual cycle.

Figure 13 (a) Zonal-mean CFC-11 mixing ratios, (b) zonal-mean zonal wind, (c) mean vertical mass flux and (d) Eliassen–Palm (EP) flux divergence, all at 50 hPa as function of latitude for EMIL simulations with the PK setup with p_Tw=400 hPa and γ=1 (red), 2 (black) and 3 K km⁻¹ (blue).

5.2 Monsoon anticyclone forced by localized idealized heating

Idealized models have been widely used to understand the basic processes occurring in the monsoon regions (e.g., Gill, 1980; Yano and McBride, 1998; Bordoni and Schneider, 2008). In particular, the development and dynamics of the monsoon anticyclones in the UTLS over Asia (e.g., Gill, 1980; Hoskins and Rodwell, 1995; Liu et al., 2007; Wei et al., 2014, 2015; Hsu and Plumb, 2000; Amemiya and Sato, 2018; Siu and Bowman, 2019) and North America (Siu and Bowman, 2019) have been investigated using simplified modeling approaches. Here, we impose an idealized heating field to force monsoonal anticyclones. The analyses presented here document the capability of the model system to apply such a forcing and to simulate anticyclones with realistic properties. This will enable more rigorous, in-depth analyses of the dynamics and of transport processes in such idealized monsoon simulations. These future analyses will exploit further capabilities of the MESSy infrastructure, in particular the inclusion of idealized and realistic tracers to study transport processes.

In the following, we show results from a T42L90MA simulation with the standard HS setup of equilibrium temperature and NH summer constellation; i.e., hemispheric asymmetry is caused by setting the asymmetry factor ϵ to −10 K. The first 2 years of this simulation have been neglected and here results from the third simulation year are presented. On top of the basic state, a regionally confined heating source is imposed in the NH tropics to subtropics (following Eqs. A12–A16 with ${\mathit{\varphi }}_{\mathrm{0}}^{\mathrm{m}}=\mathrm{20}$^∘ N, ${\mathit{\lambda }}_{\mathrm{0}}^{\mathrm{m}}=\mathrm{90}$^∘ E, $\mathit{\delta }{\mathit{\varphi }}^{\mathrm{m}}=\mathrm{10}{}^{\circ }$, $\mathit{\delta }{\mathit{\lambda }}^{\mathrm{m}}=\mathrm{30}{}^{\circ }$). In the vertical, the heating extends from $p_{\mathrm{bot}}^{\mathrm{m}}=\mathrm{800}$ hPa to $p_{\mathrm{top}}^{\mathrm{m}}=\mathrm{100}$ hPa. The heating is turned on at day 0 of the simulation with a spin-up of $t_{\mathrm{s}}^{\mathrm{m}}=\mathrm{20}$ d. Other temporal variations are not considered as $q_{\mathrm{temp}}^{\mathrm{m}}=\mathrm{0}$ K d⁻¹. The temporally constant (neglecting the spin-up period) heating is imposed with $q_{\mathrm{0}}^{\mathrm{m}}=\mathrm{8}$ K d⁻¹. After the spin-up period, the average total energy per day that is added into the model due to this additional heat source (deduced from 6 h model output) is slightly below 21×10¹⁹ J d⁻¹. This heating is of the same order of magnitude as the idealized heat source of 6×10¹⁹ J d⁻¹ prescribed in Siu and Bowman (2019) to model the North American monsoon anticyclone (see their experiments 5a–5e).

The mean geopotential height field at 100 hPa for this T42L90MA simulation with the described idealized heating is shown in Fig. 14a. A clear anticyclone is produced in response to the additional heating. This anticyclone is similar to the Asian monsoon anticyclone (e.g., Hoskins and Rodwell, 1995; Zhang et al., 2002; Randel and Park, 2006; Nützel et al., 2016). Figure 14b shows a latitude versus pressure cross section of zonal winds averaged over all longitudes. The zonal winds averaged over the anticyclone region are overlaid in black contours. The positive wind speed in the north and the negative wind speeds towards the Equator, marking the edges of the anticyclone, are clearly visible (cf. Figs. 2 and 1 in Randel and Park, 2006; Garny and Randel, 2016, respectively).

Figure 15 shows the temporal variation of the monsoon anticyclone during a 5 d period of the simulation. The daily geopotential height fields in Fig. 15 show an example of a splitting event of the anticyclone. On the first day of the depicted period, the anticyclone is elongated (Fig. 15a). After that, the anticyclone splits, and 2 d later two anticyclone centers can be identified (red dots in Fig. 15b). Then, 4 d after the elongated phase, the western center decays and the eastern center moves slightly westwards to roughly 90^∘ E (Fig. 15c). Such splitting events (sometimes also denoted as westward eddy shedding; Figs. 15 and 16 in Hsu and Plumb, 2000), as shown in Fig. 15, are (typical) features during the monsoon period (e.g., Fig. 13 in Garny and Randel, 2013; Vogel et al., 2015; Nützel et al., 2016; Pan et al., 2016).

Coincidentally, an example of eastward eddy shedding was found during the same 5 d period and is indicated via arrows in Fig. 15. This phenomenon has been previously investigated in several publications (e.g., Dethof et al., 1999; Garny and Randel, 2013; Vogel et al., 2014; Nützel et al., 2016) and constitutes a major mode of variability observed in the monsoon anticyclone. During the depicted period on the eastern edge of the anticyclone, a filament is torn off. On the first day, the anticyclone stretches to the east (Fig. 15a). Then, 2 d later, this development is even more pronounced (Fig. 15b), and again 2 d later a filament is separated from the main anticyclone (Fig. 15c).

Those examples show that the EMIL model implementation is suited to simulate a monsoon-like anticyclone reminiscent of the realistic mean state by imposing an idealized localized heating. Also, typical features of the variability of the monsoon anticyclone are reproduced in the presented EMIL simulation. Details of the variability of the idealized anticyclone, also under constant versus time-varying forcing, and its impact on troposphere–stratosphere tracer transport will be the subject of future studies.

Figure 14 (a) Mean anticyclone structure via geopotential height (km) at ∼100 hPa from 366 d of the simulation (map included for orientation and scale purposes only; i.e., the simulation features no orography). (b) Vertical cross section of zonal-mean wind (color coded) and wind in the anticyclone region averaged over 60–120^∘ E (black contours; in steps of 8 m s⁻¹; negative values are denoted by dashed lines). White contours show the maximum along the longitudes of the implied heating function (in K d⁻¹).

Figure 15 Evolution of geopotential height at ∼100 hPa showing an example of a splitting event and an eastward shedding event over a 5 d period. Red dots indicate the approximate positions of the anticyclone centers, while arrows highlight the eddy shedding event.

6 Summary and outlook

The implementation of a dry dynamical core model setup within the MESSy framework is documented. This setup, denoted EMIL (ECHAM/MESSy IdeaLized), is shown to perform consistently with established dry dynamical core benchmarks, both earlier configurations of the ECHAM core and those developed by other modeling centers.

The implementation of the submodel for temperature and wind relaxation (submodel RELAX), necessary for the dynamical core setup, includes pre-implemented functions for the parameters for Newtonian cooling and Rayleigh friction as described in HS94 and PK02. Extensions to those functions are added, namely the option to change the transition pressure between tropospheric to stratospheric equilibrium temperatures in the winter hemisphere. Further, the submodel includes the possibility to include additional diabatic heating either by pre-implemented functions for zonal mean, localized or wave-like heating, or by externally read files. Thus, the implementation provides a toolkit for the users to choose model setups to their needs.

Modifications to the setup by PK02 and Gerber and Polvani (2009), which were used frequently in the past, are presented with respect to the shape of the upper sponge layer and with respect to the equilibrium temperature profile in the winter stratosphere. The damping coefficient of the upper sponge layer is set to increase exponentially with height instead of quadratically, resembling more closely parameterized drag by GW in the full model and leading to more realistic temperature profiles in the stratopause region. However, the impact outside the sponge layer is considerable only above 10 hPa at high latitudes. Modifications of the equilibrium temperature in the high-latitude UTLS region are performed by increasing the transition pressure (p_Tw) between troposphere and stratosphere at high latitudes (thus, the decrease of temperatures forming the polar vortex starts at lower altitudes). We find that increasing the transition pressure from 100 to 400 hPa results in a realistic mean temperature profile (with the polar lapse rate γ=2 K km⁻¹), thus correcting for the UTLS warm bias of the PK02 or Gerber and Polvani (2009) simulations. With the increased transition pressure, we find a regime-like behavior of the polar vortex. The polar vortex appears to transition from a weak to a strong regime with increasing stratospheric polar lapse rate (i.e., increasing γ and p_Tw). The simulation with the most realistic mean temperature profile is at the transition point between those regimes. While we presented evidence here that the polar vortex changes reflect a regime transition, we will address the polar vortex regimes in more detail in a follow-up study.

In the past, there has been discussion (Chan and Plumb, 2009; Gerber and Polvani, 2009; Wang et al., 2012) about the regime behavior of the tropospheric jet that led to a very strong response of the tropospheric jet to changes in the polar vortex in the original work by PK02. We find a similarly strong response of the tropospheric jet to stratospheric forcing in our model under same setup as in PK02. In line with previous results that have shown sensitivity of the tropospheric jet response to the tropospheric equilibrium temperatures (Chan and Plumb, 2009), we found a strongly damped poleward shift of the tropospheric jets in response to stratospheric forcing under a setup with lower tropical upper tropospheric temperatures. We hypothesize that the more equatorward position of the tropospheric jet in those simulations, possibly related to decreased static stability in the subtropics, leads to the damped response, but this remains to be analyzed in more detail. In simulations with topographically forced planetary waves, we find a similar tropospheric jet response to stratospheric forcing to the one in Gerber and Polvani (2009), as can be expected from the previously shown robustness of the tropospheric jet behavior in different model configurations (resolution and different dynamical cores; see Wang et al., 2012). The sensitivity of the jet response to the tropical tropospheric equilibrium temperatures is lower in the topographically forced simulations, possibly because the presence of planetary waves limits the control of the subtropical static stability on the jet location. In general, it remains to be understood to what extent the regime changes of the polar vortex are connected to regime changes of the tropospheric jet, and whether those regime changes found in the idealized models are also relevant for the real atmosphere. If this behavior occurs only in idealized models, this would put their application to advance the understanding of the jet's locations in our atmosphere into question.

Simulations with planetary wave generation by topography and by wave-like heating (as suggested by Lindgren et al., 2018) are contrasted. Generally, similar basic states can be simulated with the two different setups, and in both cases increases in γ lead to increases in the polar vortex strength. However, the heating-forced simulations respond more non-linearly to increases in γ both in terms of polar vortex strength and lower stratospheric downwelling. Furthermore, while in the topographically forced simulations, both the jet maximum in the free troposphere and near the surface move poleward with a stronger polar vortex (in agreement with Gerber and Polvani, 2009), the wind maximum in the free troposphere remains at an almost constant latitude in the simulations with wave-like heating. Despite this almost constant location of the jet in the free troposphere, the jet is strongly displaced poleward near the surface in response to stratospheric forcing, transitioning from a regime of jet locations around 35–40^∘ N to a regime with jet locations north of 50^∘ N for strong polar vortex increases (similar to previously reported regime transitions, e.g., by Chan and Plumb, 2009). The prescribed wave-like heating extends throughout the free troposphere and into the lower stratosphere, leading to the hypothesis that the prescribed heating damps the ability of the wind maximum in the free troposphere to shift. Another possible explanation of the constant location of the free tropospheric jet is that it is located already at far higher latitudes than in the topographically forced simulations for weak polar vortices (see Fig. 11). Available observational evidence on troposphere–stratosphere coupling indicates that zonal wind anomalies usually occur in a vertically coherent manner, for example, due to thermal forcing by stratospheric ozone depletion (e.g., Son et al., 2010) or in connection with SSW events (e.g., Baldwin and Dunkerton, 1999). This puts the behavior of the wave-like heating experiments into question, and further work will be necessary to understand the different behavior of troposphere–stratosphere coupling in the different model versions and setups. Overall, we recommend to use the thermally forced wave generation with caution.

As a first application example of the dry dynamical core model we present, as proof of concept, a simulation with very basic chemistry (here only photolysis of CFCs), and the potential of such simulation setups to study the impact of dynamical variability and changes on tracer transport in an idealized fashion is shown. The set of chemical reactions can be expanded to the user's needs to study transport of more complex chemical tracers, such as ozone.

Secondly, we present a simulation of a monsoon-like upper tropospheric anticyclonic circulation with realistic variability forced by imposed localized heating. Such a setup can be used to study the dynamics of diabatically forced circulation systems such as monsoon anticyclones under different forcings and background states.

With the dry dynamical core model setup, the model hierarchy within the MESSy framework is extended by a commonly used model setup for studying dynamical processes. With the implementation in MESSy, the tracer utilities including the possibility to consider diagnostic chemically active tracers are available in the dry dynamical core model. As a next step, we envision an expansion to account for chemistry–dynamics interaction in a simplified manner as an intermediate step between the dry dynamical core model and the full chemistry–climate model. This next step in constructing a consistent model hierarchy of chemistry–dynamical coupling is motivated by the research question on how circulation-induced anomalies in radiative trace gases (e.g., ozone) feed back on the dynamics. This question is relevant both on climate timescales as well as on intra-seasonal timescales (e.g., during sudden stratospheric warmings). This extended setup would require radiative calculations depending on the actual tracer concentrations. While this expansion of the coupled idealized setup will be the subject of future work, we note here that all necessary components are available already in the MESSy framework. The radiation scheme from the full EMAC model (Dietmüller et al., 2016) can be used with setting the input to either the online simulated values of the trace gas of interest (i.e., ozone), while the other relevant species can be set to climatological values (e.g., water vapor) or zero (e.g., clouds and aerosols). The envisioned model setup, basically an idealized “chemistry–dynamical model”, would thus consist of a dry dynamical core with thermodynamic forcing by an idealized prescribed latent heating and radiative calculations that depend on the chemical species of interest (e.g., ozone).

Appendix A: Implemented functions in the RELAX submodel

A1 Newtonian cooling

The inverse relaxation timescale κ and the equilibrium temperature T_eq have to be specified in the model setup via the RELAX namelist file (see the Supplement). The following pre-implemented functions are available.

The functions for the HS setup, as defined by Held and Suarez (1994) but including the option of hemispheric asymmetry (as introduced by PK02), are

$$\begin{array}{}\text{(A1)}& {\displaystyle }\begin{array}{rl}& T_{\mathrm{eq}}^{\mathrm{HS}}=max\left\{T_{\mathrm{0}},\left[T_{\mathrm{1}}-{\mathit{\delta }}_y{\sin }^{\mathrm{2}}\mathit{\varphi }-\mathit{ϵ}\sin \mathit{\varphi }\right.\right.\\ & \left.\left.-{\mathit{\delta }}_z\log \left({\displaystyle \frac{p}{p_{\mathrm{0}}}}\right){\cos }^{\mathrm{2}}\mathit{\varphi }\right]{\left({\displaystyle \frac{p}{p_{\mathrm{0}}}}\right)}^k\right\},\end{array}\text{(A2)}& {\displaystyle }\mathit{\kappa }={\mathit{\kappa }}_{\mathrm{a}}+({\mathit{\kappa }}_{\mathrm{s}}-{\mathit{\kappa }}_{\mathrm{a}})max\left(\mathrm{0},{\displaystyle \frac{p/p_{\mathrm{s}}-{\mathit{\sigma }}_{\mathrm{b}}}{\mathrm{1}-{\mathit{\sigma }}_{\mathrm{b}}}}\right){\cos }^{\mathrm{4}}\mathit{\varphi },\end{array}$$

where ϕ is the geographical latitude, p is the actual pressure, p_s is the current surface pressure, and $k=R/c_p=\mathrm{2}/\mathrm{7}$. All constants can be set via namelist entries, with defaults set to the values given in HS94 (see Table S1 in the Supplement for a description of parameters and default values). The parameter ϵ sets the hemispheric asymmetry, and its sign is controlled by the namelist parameter h_fac. If h_fac is zero, the equilibrium temperature is symmetric between the hemispheres (i.e., ϵ=0). If h_fac≠0, then

$$\begin{array}{}\text{(A3)}& \mathit{ϵ}=\mathrm{sign}\left(h_{\mathrm{fac}}\right)\times \left|\mathit{ϵ}\right|;\end{array}$$

i.e., the sign of h_fac determines which hemisphere is the winter hemisphere (positive h_fac: northern hemispheric winter, negative h_fac: southern hemispheric winter).

The inverse relaxation timescale in the PK setup is identical to that in the HS setup. The equilibrium temperature in the PK setup is similar to the one of HS in the troposphere but uses the following function in the stratosphere above a given transition pressure p_T(ϕ):

$$\begin{array}{}\text{(A4)}& \begin{array}{rl}& T_{\mathrm{equ}}^{\mathrm{PK}}(\mathit{\varphi },p)=\\ & \left\{\begin{array}{l}\max \left\{T_{\mathrm{US}}\left(p_{\mathrm{Ts}}\right),\left[T_{\mathrm{1}}-{\mathit{\delta }}_y{\sin }^{\mathrm{2}}\left(\mathit{\varphi }\right)\right.\right.\\ \left.\left.-\mathit{ϵ}\sin \left(\mathit{\varphi }\right)-{\mathit{\delta }}_z\log \left(\frac{p}{p_{\mathrm{0}}}\right){\cos }^{\mathrm{2}}\left(\mathit{\varphi }\right)\right]{\left(\frac{p}{p_{\mathrm{0}}}\right)}^k\right\}& \text{for }p\ge p_{\mathrm{T}}\left(\mathit{\varphi }\right)\\ [\mathrm{1}-W(\mathit{\varphi }\left)\right]T_{\mathrm{US}}\left(p\right)+W\left(\mathit{\varphi }\right)T_{\mathrm{US}}\left(p_{\mathrm{Ts}}\right){\left(\frac{p}{p_{\mathrm{T}}\left(\mathit{\varphi }\right)}\right)}^{\frac{R\mathit{\gamma }}{g}}& \text{for }p<p_{\mathrm{T}}\left(\mathit{\varphi }\right).\end{array}\right.\end{array}\end{array}$$

The stratospheric temperature profile is based on the US standard atmosphere temperature profile T_US(p) in the summer hemisphere (U.S. Standard Atmosphere, 1976) and exhibits a temperature decrease with lapse rate γ in the winter hemisphere representing the region of the polar vortex. This transition is performed by the weighting function

$$\begin{array}{}\text{(A5)}& W\left(\mathit{\varphi }\right)={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\left[\mathrm{1}+\mathrm{sign}\left(h_{\mathrm{fac}}\right)\tanh \left({\displaystyle \frac{\mathit{\varphi }-{\mathit{\varphi }}_{\mathrm{0}}}{\mathit{\delta }\mathit{\varphi }}}\right)\right].\end{array}$$

The transition latitude ϕ₀ is set, similar to ϵ, to ${\mathit{\varphi }}_{\mathrm{0}}=\mathrm{sign}\left(h_{\mathrm{fac}}\right)\times \left|{\mathit{\varphi }}_{\mathrm{0}}\right|$. The smooth transition between tropospheric and stratospheric temperatures is ensured by bounding the tropospheric temperature to the temperature in the transition region T_US(p_Ts).

As an extension to the original PK setup, we include the possibility to vary the transition pressure from summer to winter hemisphere, using the weighting function W(ϕ):

$$\begin{array}{}\text{(A6)}& p_{\mathrm{T}}\left(\mathit{\varphi }\right)=(p_{\mathrm{Tw}}-p_{\mathrm{Ts}})W\left(\mathit{\varphi }\right)+p_{\mathrm{Ts}},\end{array}$$

where p_Ts and p_Tw are the transition pressures over the summer and winter hemisphere, respectively. Again, all constants can be set in the namelist with default values that correspond to the original PK02 setup (i.e., with constant transition pressure p_T(ϕ)≡100 hPa), as detailed in the Supplement (Table S1).

A2 Rayleigh friction

The following implemented functions are available for setting the horizontal wind damping coefficient k_damp:

1. Damping of the surface layer as specified by HS94 (option HS):

   $$\begin{array}{}\text{(A7)}& k_{\mathrm{damp}}=k_{\max }^{\mathrm{HS}}max\left(\mathrm{0},{\displaystyle \frac{\frac{p}{p_{\mathrm{s}}}-{\mathit{\sigma }}_{\mathrm{0}}}{\mathrm{1}-{\mathit{\sigma }}_{\mathrm{0}}}}\right),\end{array}$$

   with default values $k_{\max }^{\mathrm{HS}}=\mathrm{1.16}\times {\mathrm{10}}^{-\mathrm{5}}$ s⁻¹, σ₀=0.7 and p_s the current surface pressure.

2. Damping of a layer at the model top as specified by PK02 (option PK):

   $$\begin{array}{}\text{(A8)}& k_{\mathrm{damp}}=\left\{\begin{array}{ll}\mathrm{0}& \text{for }p>p_{\mathrm{sp}}\\ k_{\max }^{\mathrm{PK}}{\left(\mathrm{1.0}-\frac{p}{p_{\mathrm{sp}}}\right)}^{\mathrm{2}}& \text{for }p\le p_{\mathrm{sp}},\end{array}\right.\end{array}$$

   with default values $k_{\max }^{\mathrm{PK}}=\mathrm{2.3148}\times {\mathrm{10}}^{-\mathrm{5}}{\mathrm{s}}^{-\mathrm{1}}$ and p_sp=0.5 hPa.

3. Damping of a layer at the model top with the function as implemented in the original ECHAM code (option EH):

   $$\begin{array}{}\text{(A9)}& k_{\mathrm{damp}}=\left\{\begin{array}{ll}\mathrm{0}& \text{for }{i}_{\mathrm{lev}}>i_{\mathrm{lev}}^{\mathrm{sp}}\\ k_{\mathrm{drag}}{c}^{i_{\mathrm{lev}}^{\mathrm{sp}}-i_{\mathrm{lev}}}& \text{for }{i}_{\mathrm{lev}}\le i_{\mathrm{lev}}^{\mathrm{sp}},\end{array}\right.\end{array}$$

   where i_lev is the number of the hybrid level counted from the top of the model for a vertical resolution of L90MA. Thus, the drag k_drag is enhanced by a factor of c for each level going upward. Default values are c=1.5238, $k_{\mathrm{drag}}=\mathrm{5.02}\times {\mathrm{10}}^{-\mathrm{7}}{\mathrm{s}}^{-\mathrm{1}}$ and $i_{\mathrm{lev}}^{\mathrm{sp}}=\mathrm{10}$, corresponding to a pressure of 0.43 hPa for the L90MA vertical resolution. If the model is run at a different vertical resolution, the damping coefficients are first calculated according to Eq. (A9) for L90MA within the routine and then interpolated to the current vertical levels. The equivalent pressure levels can be found in the data set that accompanies this publication (see data availability section).

A3 Diabatic heating

The implemented heating function for the zonal-mean heating (tteh_cc_tropics), as given by Butler et al. (2010), reads

$$\begin{array}{}\text{(A10)}& \begin{array}{rl}{Q}_{\mathrm{0}}(\mathit{\lambda },\mathit{\varphi },p)=& q_{\mathrm{0}}^{cc}\exp \left[-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\left({\displaystyle \frac{\mathit{\varphi }-{\mathit{\varphi }}_{\mathrm{0}}^{cc}}{{\mathit{\delta }}_{\mathit{\varphi }}^{cc}}}\right)}^{\mathrm{2}}\right.\\ & \left.-{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}{\left({\displaystyle \frac{p/p_{\mathrm{s}}-{\mathit{\sigma }}_z^{cc}}{{\mathit{\delta }}_z^{cc}}}\right)}^{\mathrm{2}}\right],\end{array}\end{array}$$

with p_s being the surface pressure, and default values are set to those by Butler et al. (2010) (see Table S3 in the Supplement).

The temperature tendency tteh_waves, used here for the generation of planetary waves introduced by Lindgren et al. (2018), reads

$$\begin{array}{}\text{(A11)}& \begin{array}{rl}& Q_{\mathrm{0}}(\mathit{\lambda },\mathit{\varphi },p)=\\ & \left\{\begin{array}{l}{q}_{\mathrm{0}}^w\sin \left(m\mathit{\lambda }\right)\exp \left[-\frac{\mathrm{1}}{\mathrm{2}}{\left(\frac{\mathit{\varphi }-{\mathit{\varphi }}_{\mathrm{0}}^w}{{\mathit{\delta }}_{\mathit{\varphi }}^w}\right)}^{\mathrm{2}}\right]\\ \sin \left(\mathit{\pi }\frac{\log (p/p_{\mathrm{bot}})}{\log (p_{\mathrm{top}}/p_{\mathrm{bot}})}\right)\text{for }{p}_{\mathrm{top}}\le p\le p_{\mathrm{bot}},\\ \mathrm{0}\mathrm{\text{otherwise}},\end{array}\right.\end{array}\end{array}$$

where λ is the geographical longitude, and all parameters are set to default values as used by Lindgren et al. (2018); see Table S3 in the Supplement.

The function describing the localized heating field, tteh_mons, is given as

$$\begin{array}{}\text{(A12)}& Q_{\mathrm{loc}}(\mathit{\lambda },\mathit{\varphi },p,t)=Q_{\mathrm{temp}}\left(t\right)Q_{\mathrm{pres}}\left(p\right)Q_{\mathrm{lat}}\left(\mathit{\varphi }\right)Q_{\mathrm{lon}}\left(\mathit{\lambda }\right).\end{array}$$

Here, the individual factors are used to describe the temporal and spatial (horizontal and vertical) dependence of the heating function. The temporal evolution of the heating is given by

$$\begin{array}{}\text{(A13)}& \begin{array}{rl}& Q_{\mathrm{temp}}\left(t\right)=\\ & \left\{\begin{array}{ll}\frac{t}{t_{\mathrm{s}}^{\mathrm{m}}}\times (q_{\mathrm{0}}^{\mathrm{m}}+q_{\mathrm{temp}}^{\mathrm{m}}\sin (\mathrm{2}\mathit{\pi }\frac{t}{\mathit{\delta }{t}^{\mathrm{m}}}\left)\right)& \text{ for }\mathrm{0}\le t\le t_s^{\mathrm{m}},\\ \mathrm{1}\times (q_{\mathrm{0}}^{\mathrm{m}}+q_{\mathrm{temp}}^{\mathrm{m}}\sin (\mathrm{2}\mathit{\pi }\frac{t}{\mathit{\delta }{t}^{\mathrm{m}}}\left)\right)& \text{ otherwise. }\end{array}\right.\end{array}\end{array}$$

To slowly increase the heating after the start of the simulation, a spin-up factor of $\frac{t}{t_{\mathrm{s}}^{\mathrm{m}}}$ is included until the end of the spin-up time ($t_{\mathrm{s}}^{\mathrm{m}}$). After the spin-up time ($t_{\mathrm{s}}^{\mathrm{m}}$), the temporal variation of the heating is only given by a periodic oscillation (period δt^m) with amplitude ($q_{\mathrm{temp}}^{\mathrm{m}}$) around a constant base heating ($q_{\mathrm{0}}^{\mathrm{m}}$).

In the vertical, the heating is assumed to be of the same form as in Eq. (A11); i.e.,

$$\begin{array}{}\text{(A14)}& \begin{array}{rl}& Q_{\mathrm{pres}}\left(p\right)=\\ & \left\{\begin{array}{ll}\sin \left(\mathit{\pi }\frac{\log \left(p/p_{\mathrm{bot}}^{\mathrm{m}}\right)}{\log (p_{\mathrm{top}}^{\mathrm{m}}/p_{\mathrm{bot}}^{\mathrm{m}})}\right)& \text{ for }{p}_{\mathrm{top}}^{\mathrm{m}}\le p\le p_{\mathrm{bot}}^{\mathrm{m}},\\ \mathrm{0}& \text{ otherwise.}\end{array}\right.\end{array}\end{array}$$

Here, $p_{\mathrm{bot}}^{\mathrm{m}}$ and $p_{\mathrm{top}}^{\mathrm{m}}$ denote the maximum and minimum pressure to which the heating is confined in the vertical. The latitudinal dependence for $\mathit{\varphi }\in [-\mathrm{90},\mathrm{90}]$ follows the function suggested by Schubert and Masarik (2006, their Eq. 4.1) and is given as

$$\begin{array}{}\text{(A15)}& Q_{\mathrm{lat}}\left(\mathit{\varphi }\right)=\exp \left(-{\left({\displaystyle \frac{\mathit{\varphi }-{\mathit{\varphi }}_{\mathrm{0}}^{\mathrm{m}}}{\mathit{\delta }{\mathit{\varphi }}^{\mathrm{m}}}}\right)}^{\mathrm{2}}\right).\end{array}$$

Finally, the longitudinal dependence for $\mathit{\lambda }\in [\mathrm{0},\mathrm{360})$ is given by

$$\begin{array}{}\text{(A16)}& \begin{array}{rl}& Q_{\mathrm{lon}}\left(\mathit{\lambda }\right)=\\ & \left\{\begin{array}{ll}\mathrm{0.5}(\mathrm{1}+\cos (\mathit{\pi }\frac{g(\mathit{\lambda },{\mathit{\lambda }}_{\mathrm{0}}^{\mathrm{m}})}{\mathit{\delta }{\mathit{\lambda }}^{\mathrm{m}}}\left)\right)& \text{ if }g(\mathit{\lambda },{\mathit{\lambda }}_{\mathrm{0}}^{\mathrm{m}})\le \mathit{\delta }\mathit{\lambda }\\ \mathrm{0}& \text{ otherwise, }\end{array}\right.\end{array}\end{array}$$

where $g(\mathit{\lambda },{\mathit{\lambda }}_{\mathrm{0}}^{\mathrm{m}})=min\left(\right(\mathit{\lambda }-{\mathit{\lambda }}_{\mathrm{0}}^{\mathrm{m}}\left)\mathrm{mod}\right(\mathrm{360}),({\mathit{\lambda }}_{\mathrm{0}}^{\mathrm{m}}-\mathit{\lambda }\left)\mathrm{mod}\right(\mathrm{360}\left)\right)$ and the modulo function mod(360) maps ℝ to [0,360); i.e., the function returns the smallest angle between the longitude λ and the central longitude ${\mathit{\lambda }}_{\mathrm{0}}^{\mathrm{m}}$ with accounting for the crossing of the 0^∘ line. Again, the longitudinal function is based on the heating described by Schubert and Masarik (2006, their Eq. 4.1). However, as Schubert and Masarik (2006) were aiming to investigate the Madden–Julian oscillation, they included a movement of the localized heat source which we do not include here (i.e., we use their equation with propagation speed 0). Overall, this heating structure is similar to other idealized heating structures used for studying monsoon anticyclones (e.g., Siu and Bowman, 2019).

Appendix B: List of simulations

In Table B1, a list of all simulations presented in this study is given with details on their setup, resolution and simulation length. The simulations without planetary wave forcing are labeled “WN0”, the ones with topographic wavenumber 2 wave forcing with “WN2T” and the ones with diabatic wave forcing with “WN2H”. The values of the winter transition pressure p_Tw are given in hPa, and the values of the polar vortex lapse rate γ in K km⁻¹. “npv” refers to simulations with “no polar vortex”; i.e., the summer stratospheric equilibrium temperature is extended to the winter pole. The values of the parameter in control of the tropical vertical temperature gradient in the troposphere (δ_z) are given in K. As stated in Sect. 2, the reduced value of δ_z=10 K$/\log \left(\mathrm{10}\right)=\mathrm{4.34}$ K in a number of simulations is the result of a previous implementation of the term ${\mathit{\delta }}_z\log \left(\frac{p}{p_{\mathrm{0}}}\right){\cos }^{\mathrm{2}}\mathit{\varphi }$ in Eq. (A1) as logarithm to base 10 instead of the natural logarithm. Note that in the accompanied data set (see data availability), the simulations are labeled accordingly as “ln” versus “log10”. The upper atmospheric damping coefficients are set to the formulation by PK (see Eq. A8) or to the formulation EH (see Eq. A9). The total length of the simulations is given in the table in days, with the number of analyzed days given in the figure captions.

Table B1 List of all simulations presented in this study, with information on all relevant parameter settings and variations, and simulation length in days (see text for details). Simulations with multiple values of γ are listed in one row for brevity. ^a Also p_Ts=200 hPa. ^b Including additional localized heating with parameter settings given in Sect. 5.2.

Appendix C: Probability distributions of polar vortex strength and tropospheric jet location

The probability distributions of the polar vortex strength (maximum zonal-mean zonal wind at 10 hPa) and the position of the tropospheric jet (latitude of zonal-mean zonal wind maximum at 850 and 500 hPa) are shown for a number of selected experiments: for the original PK02 setup, and the equivalent reduced δ_z setup (Fig. C1), for the simulations with a WN2 topography with differing p_Tw (Fig. C2), and for simulations with a WN2 diabatic heating for the standard and reduced δ_z setups (Fig. C3).

Figure C1 Probability distributions of (a, d) the maximal zonal-mean zonal wind at 10 hPa, (b, e) the latitude of the maximum of zonal-mean zonal winds at 850 hPa, and (c, f) at 500 hPa for simulations without planetary wave forcing, with (a, b, c) for the original PK02 setup and (d, e, f) for the modified setup with reduced δ_z. The mean value of the distribution is given in the legend (denoted as μ).

Figure C2 As Fig. C1 but for simulations with the standard PK setup and WN2 topography (a, b, c) for p_Tw=100 hPa and (d, e, f) for p_Tw=400 hPa.

Figure C3 As Fig. C1 but for simulations with WN2 heating and p_Tw=400 hPa for (a, b, c) the standard setup and (d, e, f) the modified setup with reduced δ_z.