© M.J. Griffith et al., Published by EDP Sciences 2020

1 Introduction

An important focus of many weather forecasting organisations is the development of a complete Sun-to-Earth model in order to enhance forecasting of space weather, and ultimately develop a fully coupled system describing the Earth’s atmosphere (e.g. Tóth et al., 2005). This is particularly since it is now known that the lower atmosphere is an important driving factor in variability and space weather in the upper atmosphere, in addition to solar variability (Akmaev, 2011). It is for this reason that the Met Office wishes to extend its weather and climate model, the Unified Model (UM), into the upper atmosphere (above 85 km).

One of the Met Office’s long term aims to achieve this fully coupled system is the development of the UM into a whole atmosphere model (e.g. Akmaev, 2011; Jackson et al., 2019) – namely one that simulates the Earth’s atmosphere from the surface up to the exobase (≈ 600 km). Such a model is crucial in obtaining accurate prediction of the upper atmosphere as the model resolution increases (e.g. Immel et al., 2006; Yue et al., 2016). This is due to the influence of vertically propagating atmospheric waves, in particular tides, that grow roughly exponentially with height because of the corresponding exponential decrease of the background density. Hence, these waves become particularly large in the Mesosphere and Lower Thermosphere (MLT), with associated influences on the ionosphere and thermosphere at higher altitudes.

Energetic particle precipitation (EPP) and changes in solar radiation associated with the solar cycle can also impact the troposphere. EPP events can lead to changes in MLT and stratospheric NOx and ozone concentrations which can modulate polar surface air temperatures by affecting the radiative budget with consequent effect on atmospheric circulation patterns (e.g. Seppälä et al., 2009; Päivärinta et al., 2013). Ineson et al. (2011) indicates the impact of solar minimum and solar maximum on Northern Hemisphere winter surface temperatures via changes to the North Atlantic Oscillation.

Thus, these waves are influential in atmospheric dynamics at all heights, and so their inclusion is important not only as drivers of upper atmosphere circulation, but also in terms of improving accuracy for current lower atmosphere models.

The forecasting of this region is important for real world applications such as radio communication, satellite drag calculation and Global Navigation Satellite System Positioning, Navigation and Timing (GNSS PNT). These can all be influenced by atmospheric waves in the MLT, and so modelling this region can help mitigate the potential disruption to these systems.

Initial runs of the UM with an increased model lid height performed by Harry (2015) became unstable and crashed unless considerable damping was used. These runs tested the UM with both 100 and 120 km model lids, with experiments focused on tuning parameters to obtain the longest run time possible. In particular, the cause and nature of the instability remained undiagnosed.

These initial results relied on changes in damping using an implicit weighting parameter α (described later) and changes to the time step Δt to produce model runs which were not reliably stable. Both of these modifications are undesirable for an accurate and efficient model. Increasing the implicit weighting parameter α damps out high frequency wave components globally and thus degrades the solution over the entire model domain, whilst decreasing the time step $\Delta t$ leads to a more costly model which is impractical for real world application. The choices made in this research are made with these factors in mind, and are discussed further in Section 2.2.3.

Here, we advance the development of a whole atmosphere model by investigating the instability when the model lid is raised into the MLT, raising the roof of the current 85 km lid. This is achieved by first investigating the output produced using a raised upper boundary of 100 km and considering key prognostic variables such as winds and temperature, as well as diagnostic variables such as short wave radiative heating.

This series of diagnostic tests yielded the important result that the instability is not a result of numerical errors, but is the result of an incorrect application of the parameterization schemes in the MLT, where their physical approximation is no longer valid. By a careful inspection of the parameterization schemes, the cause of the instability was shown to lie in the use of the Local Thermodynamic Equilibrium (LTE) assumption used in the radiation scheme.

New UM chemistry and radiation schemes to supply appropriate radiative forcing in the MLT were not ready at the time of performing this research. Here, we wish to circumvent these delays and press ahead with our aim of producing a stable extended UM. We achieve this by switching off the UM chemistry and radiation schemes (the latter above 70 km in altitude only) and relaxing the model to a prescribed climatological temperature profile. This relaxation approximates the impact of the omitted radiation and chemistry and follows the approach used in other atmospheric models (e.g. Telford et al., 2008).

The non-orographic gravity wave forcing scheme was also suspected to be poorly optimised for the high atmosphere. This scheme is based on the Ultra Simple Spectral Parameterization (USSP) developed by Warner & McIntyre (2001). Several runs comparing the model with this scheme on and off were performed. In agreement with the results of Harry (2015), experimentation revealed that the scheme did not have a large impact on the stability of the UM. Thus, this scheme is left unchanged throughout. The scheme will undoubtedly need modifications to provide accurate non-orographic gravity wave forcing in the MLT, but we do not discuss this in this paper.

The layout of the remainder of this paper is as follows. Section 2 contains a description of the UM and a discussion of parts of the model particularly relevant when considering the instability in the extended model. The instability is diagnosed in Section 3 and is consequently resolved in Section 4 using relaxation to a climatological temperature profile. Further extensions to this stable “nudged” model are made in Section 5 and finally conclusions are drawn along with ideas for future work in Section 6.

2 Background

2.1 The Unified Model

As in any General Circulation Model (GCM), the UM is split into two core sections, one that describes atmospheric dynamics (the dynamical core) and one that describes atmospheric physics (parameterizations). The original UM is documented by Cullen (1993). A new dynamical core was introduced in the early 2000s (Davies et al., 2005) and the UM’s current dynamical core, ENDGame, is described by Wood et al. (2014).

The ENDGame dynamical core solves the non-hydrostatic, fully compressible deep-atmosphere equations of motion on a rotating sphere using a semi-implicit semi-Lagrangian formulation. The primary prognostic variables of the three-dimensional wind components, virtual dry potential temperature, Exner pressure and dry density are used whilst moisture prognostics are advected as free tracers. Importantly this non-hydrostatic framework allows for vertical acceleration; an important consideration in the forcing of vertical winds, which are typically larger in the upper atmosphere than those seen in the lower atmosphere. The discretised equations are solved using an iterative implicit method – more details of which can also be found in Wood et al. (2014).

The model discretisation is split up into horizontal and vertical components. We fix the horizontal resolution at 1.25°N × 1.875°E, or the so called N96 resolution.¹ In the vertical, an 85-level set labelled L85 is used, which has 50 levels below 18 km, 35 levels above this and a model lid 85 km above sea level. This level specification follows a roughly exponential distribution in that the level depth exponentially increases with increase in height. We will add more levels to raise the upper boundary and investigate whether this level specification is appropriate for accurate modelling of the MLT in Section 5.

There are many atmospheric processes which operate on smaller length scales than those used in the dynamical scheme, and therefore are not resolved on the numerical grid. Furthermore, some physical processes such as radiation are not represented in the equations of motion used in the UM. Such processes are therefore approximated within the atmospheric physics section of the model by parameterization. These parameterizations are key to obtaining realistic and useful forecasts for the operational model, and their addition to the dynamics provides a reasonably complete description of the atmosphere.

Amongst these parameterizations, it is the radiation and chemistry scheme which are most relevant for this paper, and thus we outline their parameterization here.

2.1.1 The radiation scheme

The purpose of the radiation scheme, which is based on the work of Edwards & Slingo (1996), is to simulate the effect of the Sun on the Earth’s atmosphere. It parameterizes the effect of incoming short wavelength (0.2–0.5 μm) radiation and the subsequent emission of long wavelength radiation (>3 μm) into space, along with the redistribution of heat within the atmosphere.

In order to do this, radiative fluxes (i.e. the amount of energy transferred by incoming solar radiation per unit area per unit time) are calculated and are modelled as upward and downward fluxes of energy relative to a model grid point. From these fluxes, quantities such as atmospheric heating rates are derived, which can then be used to influence the model.

Also captured are the effects of various radiative processes, such as the absorption and emission of radiation by various airborne molecules. In particular, absorption of short wave (SW) and long wave (LW) radiation by various atmospheric constituents must be considered, as this is important in achieving the correct atmospheric temperature structure.

However, the current implementation of the radiation scheme makes the assumption of LTE – the condition under which matter emits radiation based on its intrinsic properties and its temperature, uninfluenced by the magnitude of any incident radiation. This assumption breaks down above ≈65 km (Fomichev et al., 2004). In addition, the radiation scheme does not consider wavelengths less than 0.2 μm, which are significant for radiative heating (see e.g. Chamberlain & Hunten, 1987) and to drive exothermic heating from chemical reactions in the MLT. Both of these factors make the current radiation scheme unsuited for use in the MLT. We diagnose this problem in Section 3 and address the consequent stability issues with the introduction of a relaxation scheme to a climatological temperature profile in Section 4.

2.1.2 The chemistry scheme

The UKCA (United Kingdom Chemistry and Aerosols; Morgenstern et al., 2009 and O’Connor et al., 2014) is a framework that encompasses several atmospheric chemistry and aerosol schemes, which is broadly referred to as the chemistry scheme. These schemes are responsible for modelling the effects of airborne molecules and aerosols, and their transport, interactively. They also account for the interaction of these particles with other parameterization schemes such as the radiation scheme described above.

However, the UKCA does not currently² include the relevant reactions above the 85 km model lid height. Chemical heating from exothermic reactions is key for the large increase in temperature with altitude seen in the lower thermosphere (see e.g. Marsh et al., 2007). However, these reactions are not currently represented in the UKCA, nor are their corresponding photolysis rates included in the radiation scheme. Thus, the chemistry scheme is switched off and in its place, monthly zonal mean profiles are provided to represent important missing atmospheric constituents. Due to its importance as a significant heat source in the stratosphere, profiles for atmospheric ozone are provided to retain this effect in the model. The various atmospheric aerosols in the model primarily influence only tropospheric dynamics, and so these are switched off in the extended upper boundary runs.

2.2 Stability mechanisms in the Unified Model

2.2.1 Model damping/sponge layer

In atmospheric models, it is often necessary to add a damping or sponge layer to the upper boundary of the model in order to prevent spurious reflection of vertically propagating waves from the model lid, and resulting numerical instabilities. These occur due to the no-flow (or zero vertical velocity) upper boundary condition which results in a rigid, reflecting lid.

It is also necessary to add damping at the poles to deal with numerical instabilities caused by the clustering of points due to the latitude-longitude grid, see e.g. Thuburn & Staniforth (2004).

In the UM, these instabilities are dealt with by adding a simple damping term −μw to the forcing term of the vertical wind (w) equation – see Wood et al. (2014) for an in depth explanation of how it is incorporated into the model setup.

As an example, μ can be defined as a height dependent function³ μ = μ(z) so that the sponge is implemented gradually from the base of the sponge layer z _B until the top of the model z _TOP, namely

$$\mu \left(z\right)=\{\begin{array}{ll}0,& \mathrm{if} z<z_{\mathrm{B}}\\ \overline{\mu }{\sin }^2\left[\frac{\pi }{2}\left(\frac{z-z_{\mathrm{B}}}{z_{\mathrm{TOP}}-z_{\mathrm{B}}}\right)\right],& \mathrm{if} z_{\mathrm{B}}\le z\le z_{\mathrm{TOP}},\end{array}$$where $\overline{\mu }\in \left[\mathrm{0,1}\right]$ is the vertical damping coefficient which can be tuned according to the application to increase/decrease the strength of the sponge layer – the default value is $\overline{\mu }=0.05$ s⁻¹.

This example however does not incorporate the latitudinal variation that can also be added to the sponge layer to introduce damping around the problematic polar regions. The form of the sponge layer implemented in the default configuration of the UM provides damping of strength $\overline{\mu }$ towards the upper boundary of the model domain, as well as providing damping of strength $2\overline{\mu }$ at the poles for all altitudes. This is depicted for a 100 km model lid in Figure 1.

Fig. 1 The default sponge layer with vertical damping coefficient $\overline{\mu }=0.05$ s−1 and 100 km model lid. The sponge layer can be seen to be at $\overline{\mu }=0.05$ s−1 for most latitudes and for high altitudes but is doubled at the polar latitudes for all altitudes to help numerical stability in this area.

This parameter $\overline{\mu }$ can be tuned appropriately to deal with larger vertical velocities at the upper boundary as the upper boundary height and number of vertical levels increase. The downside to increasing $\mu$ in this way is that it introduces artificial damping. However, this damping is highly localised to the model top and poles, and so is favourable to applying the numerical damping discussed in Section 2.2.3.

As well as this, it is already used in the standard UM, and so it is natural to include it in our experiments with an extended upper boundary. We will see that the damping must be increased in order to stabilise the model with the finer vertical resolutions and greater upper boundary heights used in Section 5.

2.2.2 Model halos

The UM is parallelised for use on a large high performance computer. The model domain is consequently split up across separate processors; each processor has access to a limited latitude-longitude square on its main CPU grid. The halo or ring of neighbouring grid points is then put in place to account for advection outside this region.

However, if the horizontal wind is large enough, the advection scheme will try to access data that is very far away from the arrival point and outside this additional halo. If this occurs, then the model fails with a halo error. For global model runs, this problem only occurs in the North-South wind direction – the periodic structure in the East-West direction allows data to be communicated from one processor to another whenever necessary for East-West winds. A schematic of this problem can be seen in Figure 2.

Fig. 2 Schematic illustrating the halo error. On the right hand side of the figure, the advection scheme can be seen trying to access data that is outside of the halo region which causes an error. The halo size can be increased to accommodate the larger winds present in the upper atmosphere.

This halo can be extended in order to deal with larger wind speeds, with some impact on computational efficiency. However, this importantly does not affect the model solution, and is a necessary additional expense to deal with the larger wind speeds encountered on raising the model lid (see for example those observed in Fig. 1 of Hedin et al., 1991). We will see that the extension of this halo is also necessary with the finer vertical resolutions and greater upper boundary heights used in Section 5.

2.2.3 Numerical stabilisation

The model can also be stabilised more directly by adjusting the implicit weighting parameter α as well as the time step $\Delta t$. This implicit weighting parameter (or off-centring parameter) is used in the approximation made in the semi-implicit semi-Lagrangian discretization – discussed in Wood et al. (2014) and in further detail in e.g. Rivest et al. (1994). However, both of these options are detrimental to the model as a whole.

Making the model more implicit by increasing the implicit weighting parameter α improves the stability of the model by introducing more implicit damping into the model. However, it therefore also decreases the accuracy of the entire solution because waves (in particular their high-frequency components) are artificially damped throughout the whole atmosphere, without physical motivation.

Decreasing the time step $\Delta t$ allows the model to cope better with larger wind speeds, as the advection scheme then advects less per time step. However, this is only really a viable solution for research purposes. For use as an operational forecasting model, the run time must be kept short to achieve timely forecasts. This is an important motivation for retaining a larger time step.

We favour the more tailored approaches discussed above and avoid changing α and Δt for the extended UM experiments.

3 Diagnosing the instability

We start with the standard General Atmosphere (GA) 7.0 configuration of the UM as described in Walters et al. (2017), with the standard 85 km lid. We begin by raising the model lid to 100 km in order to determine the nature of the instability causing the model to crash. To do this, we undertake an empirical analysis of model runs of the UM using different configurations.

3.1 Extended model setup

The changes that are made to the aforementioned standard GA7.0 configuration of the UM are summarised below

1. The model chemistry scheme is entirely switched off – this is not yet configured to work above 85 km.

2. Atmospheric aerosols are switched off and ozone background files are switched on, accounting for the lack of the chemistry scheme.

3. The model upper boundary is raised from 85 to 100 km.

Otherwise, we conserve as many model features as possible in order to appropriately compare the output to that of the 85 km model. For the latter of these points, a new vertical level set is provided that extends the 85 vertical levels in the 85 km model up to 88 vertical levels in the 100 km model. We adhere to the default level spacing to add an additional two levels, and then add the final level at 100 km. A comparison of the two sets of vertical levels can be seen in Figure 3.

Fig. 3 Vertical level sets for the standard model height (black) and those added for the raised model height (blue). Two levels are added conforming to the 85 km level distribution, whilst the final is placed at 100 km.

This level spacing – where the levels increase in depth according to an approximately exponential distribution – is chosen with the lower atmosphere in mind; the majority of levels are placed in this region and the vertical resolution at the top of the model is much larger – nearing 6 km. To diagnose the cause of the aforementioned instability, we keep this model level spacing so as to keep the first 85 km of the model domain the same and provide a good comparison between the configurations. However, to obtain a model that resolves appropriate wave scales in the MLT, we must have a finer vertical resolution. We will investigate different vertical level sets that are better suited to achieving this objective in Section 5.1.

As discussed previously, monthly zonal mean profiles are prescribed for atmospheric ozone, which are extended to have an upper boundary at 100 km. An example profile for June is depicted in Figure 4.

Fig. 4 A zonal mean profile produced to prescribe the initial configuration of atmospheric ozone for the 100 km model lid. Its inclusion is important in order to represent the correct vertical temperature structure (via radiative heating).

Finally, it is important to note that in all cases, the experimental runs are performed multiple times in order to validate the consistency of the output.

3.2 Extended model results

We summarise the runs performed using the extended UM in Table 1.

In both model runs, an error related to the North-South advection scheme is produced and causes the model to crash. This means that the prognostic variables related to a grid point have been advected further than the data halo for the given grid point. This is indicative of excessive wind speeds in the model. This instability is localised in space and time. We proceed to investigate these abnormal wind speeds in order to determine their cause.

To do this, we consider the September model run, which crashed in the following December, therefore running with no issues and normal model fields for several thousand model time steps. We subsequently plot instantaneous westerly (u) and southerly (v) winds at the last time step before the model crash.

The problem arises over the South Pole of the model. Thus, we take a stereographic plot over this pole, as well as a slice through the atmosphere at this most southerly latitude circle (i.e at 89.375°S for u wind and 90°S for v wind). This can be seen in Figure 5. Here, the greatly excessive wind speeds are evident, on the order of 600 m/s for the $u$ wind and 400 m/s for the $v$ wind. This is in comparison to $u$ winds of at most 90 m/s in e.g. Swinbank & Ortland (2003) and $v$ winds of at most 100 m/s in e.g. Hedin et al. (1991). The problematic winds can also be seen to be constrained to the uppermost model level.

Fig. 5 Modelled wind in the 100 km September model run at the time step before the model crash. Stereographic plots at the uppermost model level (a) and (b). Plots of the wind at the southernmost latitude circle (c) and (d). This makes evident the abnormally high wind speeds causing the model to crash.

The same phenomenon is observed for the March model run which crashes in the following June, but with the problem occurring around the North Pole.

It is evident therefore that the problematic model winds occur in the summer season for a particular hemisphere. Namely, in the southern hemisphere summer (December) for a September start date and the northern hemisphere summer (June) for a March start date. We deduce that the problem is therefore a seasonal one.

At this stage, it is important to consider that the instability could be numerical in nature. On increasing the model upper boundary, larger density perturbations associated with exponential wave growth are introduced which could cause the underlying numerical scheme to become unstable. However, the long run times of the model and the seasonal pattern of the crash dates indicate that the observed instability is dominated by problems in the physical schemes used in the model.

It is also pertinent to observe that this instability is local in nature. A global instability is more likely to be caused by problems in the dynamical core of the model, whereas a local instability is more illustrative of a problem with localised forcing provided by the parameterization schemes in the model. We indeed see localised instabilities, which consistently occur in the polar regions of the model.

Both of these factors mean we are confident in excluding numerical issues as the cause of the instability in the first instance. We therefore focus our investigation on the seasonal forcing provided by the radiation scheme – which we already suspect will have problems when extended to the MLT.

The radiation scheme provides the solar energy input responsible for seasonal weather variations. Hence, we look at the SW radiation provided by the radiation scheme which governs this input of energy through heating.

We plot the zonal mean of instantaneous SW radiative heating at the last time step before the model crashes for both start dates in Figure 6. The heating from the ozone layer is evident in the summer hemisphere of both figures, but more pertinent are the abnormal values at the top of the model. The values are focused in the summer pole of the uppermost layer; exactly the region where the problematic winds were previously observed. Thus, we deduce that it is this incorrect SW radiative forcing that is driving the high polar wind speeds and ultimately therefore causing the observed instability in the UM with a raised upper boundary.

Fig. 6 Plots illustrating the issue evident with the SW radiation produced by the radiation scheme. We see abnormal values in exactly the regions where the high wind speeds occur in the previously discussed model runs. (a) Zonal mean of the SW radiative heating for the model run starting in September and crashing in December. (b) Zonal mean of the SW radiative heating for the model run starting in March and crashing in June.

The observed SW radiative heating is in fact in keeping with LTE. However, this assumption is no longer valid in the MLT, as is demonstrated in Figure 7.

Fig. 7 Total heating due to solar energy absorption in the near-IR CO2 (NIR) bands. This illustrates the difference in heating in K day−1 when the Local Thermodynamic Equilibrium (LTE) approximation is made. This approximation breaks down at around 0.1 hPa (approx. 65 km) as can be seen. After Fomichev et al. (2004).

We strongly suspect that it is indeed the assumption of LTE that is leading to the observed instability in the UM, and that non-LTE (NLTE) effects must be considered when extending the upper boundary to the MLT.

4 Resolving the instability

The introduction of a NLTE radiation scheme is an important step towards a whole atmosphere model. However, such a scheme was not ready at the time of performing this research. In the meantime, it is important to verify that the radiative forcing is the sole cause of the observed instability, and not one of several factors. We confirm this by verifying that with appropriate radiative forcing in the upper atmosphere, the model is stable under the current 100 km model configuration.

To supply appropriate radiative forcing for the upper atmosphere without access to the NLTE radiation scheme, it is necessary to develop the model to account for NLTE effects which become prevalent in the MLT. Thus, an interim solution of relaxation – or “nudging” – to an analytic temperature profile is added to the model to confirm that the extended UM is stable under corrected forcing. We shall refer to the model with this nudging included as the nudged model.

4.1 Nudged model setup

The erroneous radiation scheme is replaced at high altitudes by relaxation to a temperature profile based on climatology. The scheme begins to drift away from LTE at around 65–70 km, and so it is from 70 km and above that the scheme is replaced. For simplicity, this profile is globally uniform and varies only with height so that T _PROFILE = T _PROFILE(z). This is sufficient as an interim solution before the introduction of the full NLTE radiation scheme.

This relaxation acts to push – or “nudge” – the temperature field from the model dynamics toward the chosen temperature profile, over a given timescale. This timescale τ is chosen as to allow some variation about this profile, so that the dynamics of the model can also influence the temperature field.

More concretely, the nudged temperature field is computed in an additional step after each dynamics step as

$$T_{\mathrm{NUDGED}}={\alpha }_{\mathrm{NP}}{T}_{\mathrm{PROFILE}}+(1-{\mathrm{\alpha }}_{\mathrm{NP}})T_{\mathrm{DYN}},$$where $T_{\mathrm{NUDGED}}$ is the nudged temperature field, $T_{\mathrm{PROFILE}}$ is the chosen analytic profile, $T_{\mathrm{DYN}}$ is the temperature field computed by the dynamics and ${\alpha }_{\mathrm{NP}}=k\Delta t/\tau$ is the nudging parameter. The timescale $\tau$ is chosen to be 24 hours, in keeping with work by Song et al. (2018). Here, using a hydrostatic model, a physical reasoning was used to determine appropriate relaxation timescales for a high atmosphere model being nudged to data. They found that between 8, 24 and 40 h timescales, the 24 h gave the most reasonable model fields and conservation properties out of the three options. Thus, this informed our choice of using τ as 24 h. Finally, $k\in \left[\mathrm{0,1}\right]$ is the nudging ramp parameter – this gives a smoothed transition to the nudging as it is introduced at 70 km.

The temperature profile $T_{\mathrm{PROFILE}}$ used is based on climatological and satellite data so that

• Between 70 and 86 km the profile is based on the US Standard Atmosphere (USSA) (COESA, 1976).

• Between 86 and 119.7 km the profile is based on the Committee on Space Research (COSPAR) International Reference Atmosphere (CIRA) (Fleming et al., 1990).

• Above 119.7 km the temperature asymptotes to a selected exobase temperature (here 1000 K).

This leads to the expression given in Equation (1) for this profile

$$T_{\mathrm{PROFILE}}\left(z\right) =\{\begin{array}{ll}{T}_{70\mathrm{ }\mathrm{km}}-{\mathrm{\Gamma }}_{\mathrm{MESO}}\left(z-z_{\mathrm{USSA}-\mathrm{bottom}}\right)& , \mathrm{if} 70\le z\le 86 \mathrm{km}\\ T_{70 \mathrm{km}}-{\mathrm{\Gamma }}_{\mathrm{MESO}}(z_{\mathrm{USSA}-\mathrm{top}}-z_{\mathrm{USSA}-\mathrm{bottom}})& \\ -{\mathrm{\Gamma }}_{\mathrm{CIRA}}\left(z-z_{\mathrm{USSA}-\mathrm{top}}\right)& , \mathrm{if} 86<z\le 93.3 \mathrm{km}\\ T_{70 \mathrm{km}}-{\mathrm{\Gamma }}_{\mathrm{MESO}}(z_{\mathrm{USSA}-\mathrm{top}}-z_{\mathrm{USSA}-\mathrm{bottom}})& \\ -{\mathrm{\Gamma }}_{\mathrm{CIRA}}\left(z-z_{\mathrm{USSA}-\mathrm{top}}\right)& \\ +(3\times 10^{-7})(z-z_{\mathrm{Tmin}-\mathrm{CIRA}}{)}^2& , \mathrm{if} 93.3<z\le 119.7 \mathrm{km}\\ T_{\mathrm{exobase}}-\left(T_{\mathrm{exobase}}-T_{\mathrm{CIRA}-\mathrm{top}}\right)\exp \left(-{\mathrm{\Gamma }}_{\mathrm{THERMO}}\sigma \left(z\right)\right)& , \mathrm{if} z>119.7 \mathrm{km}\\ & \end{array}$$(1)where

$$\sigma \left(z\right)=\frac{\left(z-z_{\mathrm{CIRA}-\mathrm{top}}\right)\left(R_{⨁}+z_{\mathrm{CIRA}-\mathrm{top}}\right)}{\left(z+R_{⨁}\right)}$$with the radius of the Earth denoted $R_{⨁}$ and with the following parameter values:

The $T$ values represent temperatures, the $\mathrm{\Gamma }$ values represent lapse rates and the $z$ values represent heights, all at locations indicated by the subscripts. Note in particular the exobase temperature at 1000 K. In reality this can vary from 700 to 1600 K dependent on the solar cycle. This parameter can be easily tuned but is just set to a typical value here, given that the upper boundaries considered are much lower than the exobase.

The resulting profile is plotted in Figure 8. The main features of the profile are that it lapses linearly up to $\approx$ 93 km, then is quadratic to $\approx$ 120 km (including the mesopause temperature minimum) and finally asymptotes towards the exobase temperature.

Fig. 8 Temperature profile used for the relaxation of the temperature field.

Importantly, the chosen temperature profile shows good agreement with the standard 85 km model temperatures at 70 km, where the nudging is introduced, as can be seen in Figure 9. This therefore minimises the possibility of discontinuities in the temperature field in the area of transition.

Fig. 9 Comparison of the 85 km temperature field for various latitudes and the forcing profile (a), with pertinent area zoomed in (b). Good agreement between the two is seen.

Finally the erroneous values from the radiation scheme are zeroed above 70 km. This removes the erroneous forcing that was evident in Figure 6.

4.2 Nudged model results

We now perform the same model runs as in Section 3.2 using the nudged model. Now the model differs from the standard GA7.0 configuration of the UM in the following ways:

1. The model chemistry scheme is entirely switched off – this is not yet configured to work above 85 km.

2. Atmospheric aerosols are switched off and ozone background files are switched on, accounting for the lack of the chemistry scheme.

3. The model upper boundary is raised from 85 to 100 km.

4. The forcing from the radiation scheme is zeroed above 70 km.

5. The temperature field above 70 km is nudged towards a prescribed climatological temperature profile.

Importantly, with this in place, the model is consequently stable for both start dates and can be run for several years without the appearance of the instability observed in Section 3. In Figure 10, we plot instantaneous $u$ and $v$ winds produced using the September nudged model run for comparison with the winds from the initial extended model run in Figure 5. These are plotted at the same model time as previously, namely the December crash date of the initial extended model. We also plot the winds in the same manner – stereographic plots over the South Pole, as well as a slice through the atmosphere at the most southerly latitude circle (at around 90°S).

Fig. 10 Modelled wind in the 100 km September model run with the nudging scheme in place. Stereographic plots at the uppermost model level (a) and (b). Plots of the wind at the southernmost latitude circle (c) and (d). We see the removal of the anomalous wind speeds present previously by comparison with Figure 5.

In comparing the two figures, it is evident that the instability in the polar region has been removed, with the resultant winds of a much more reasonable order of magnitude; the 400–600 m/s winds are no longer evident and the unphysical winds of opposite signs in neighbouring grid cells are also removed. The winds are also much more on par with the magnitudes seen in Swinbank & Ortland (2003) and Hedin et al. (1991).

We also see that the excessive SW radiative heating has been removed in Figure 11.

Fig. 11 SW radiation showing zeroed field above 70 km.

These results indicate that the LTE assumption in the radiation scheme is indeed responsible for the model instability with a lid at 100 km and importantly, we are able to conclude that there are no additional factors causing instability with the current configuration – such as problems resulting from the polar singularity of the latitude-longitude numerical discretisation.

We would now like to build on this stable configuration and extend the model further to different vertical resolutions and upper boundary heights. However, it is first diligent to compare important model fields with the 85 km model to illustrate that the nudging is not having a detrimental effect on the modelled fields in the lower atmosphere.

To do this, we compare climatological fields rather than the instantaneous fields presented previously. Thus, we present monthly mean zonal mean plots for westerly ($u$) wind, vertical ($w$) wind and temperature from the nudged model. These are plotted for December in the second column of Figure 12 and for June in the second column of Figure 13. They are compared to fields from the 85 km model in the first column of the respective figures.

Fig. 12 Zonal mean monthly mean climatologies in December of westerly ($u$) wind (top), vertical ($w$) wind (middle) and temperature (bottom) comparing the original 85 km model with the 88, 94 and 113 level configurations.

Fig. 13 Zonal mean monthly mean climatologies in June of westerly ($u$) wind (top), vertical ($w$) wind (middle) and temperature (bottom) comparing the original 85 km model with the 88, 94 and 113 level configurations.

We observe good consistency of the model fields for both December in Figure 12 and June in Figure 13 up to around 70 km where the nudging is implemented. At this height, there is evidence of the globally uniform approximation of the temperature field in comparing plot (i) with plot (j) in both figures at around 70 km. Here, we observe temperature changes on the order of 10–20 K at the summer poles in both figures. However, this is to be expected given the simplified physical approximation made in the temperature profile. We reiterate that the nudging is put in place to observe if the model is consequently stable rather than as a permanent solution – the more physically accurate NLTE radiation scheme and the updated UKCA chemistry scheme will address this issue.

We comment on the other plots in this figure (with different vertical resolutions) in Section 5.1.

At this stage, we also make a preliminary comparison with data in order to test the realism of the nudged model – in reality given the simplified physical approximations made we expect the model to drift away from realistic values.

We compare climatologies for monthly mean zonal mean westerly ($u$) wind and temperature for the 100 km lid model run against data. We use the 3 km vertical resolution 94 level configuration which will be developed in Section 5 for comparison, given the finer structure observed in the model fields. The wind data used for comparison is the Upper Atmosphere Research Satellite (UARS) Reference Atmosphere Project (URAP) (Swinbank & Ortland, 2003). The temperature data used for comparison is from the Earth Observing System Microwave Limb Sounder (EOS MLS) on the Aura satellite (Waters et al., 2006). The data is recorded on pressure levels and we make an approximate conversion to height using $z=-H\ln (p/p_0)$, with $p_0$ = 1000 hPa and H = 7200 m.

We plot the climatologies and data for December in Figure 14 and for June in Figure 15. A more in depth comparison to data will come as the model is improved further in the MLT.

Fig. 14 Zonal mean monthly mean climatologies in December of westerly ($u$) wind (top) and temperature (bottom). These are compared to URAP data for zonal wind and MLS data for temperature. The approximate height is calculated using a scale height of H = 7200 m.

Fig. 15 Zonal mean monthly mean climatologies in June of westerly ($u$) wind (top) and temperature (bottom). These are compared to URAP data for zonal wind and MLS data for temperature. The approximate height is calculated using a scale height of H = 7200 m.

A slight temperature gradient towards the summer pole is seen at around 90 km in the modelled temperature, indicative of a cold summer mesopause. This is however less pronounced than the cold summer mesopause in the MLS data, but is encouraging given the simplified physics. It is also worth noting that there are also doubts in the robustness of the MLS data at low pressures (e.g. Schwartz et al., 2008), meaning that the cold summer mesopause is not necessarily as cold as indicated by the MLS data, but is certainly colder than the model currently predicts. In the model winds, the summer-time westerly mesospheric jet is not present at 90 km when compared to the URAP data, however this is to be expected given the lack of a strong temperature gradient provided by a cold summer mesopause.

It is evident therefore that the model still requires several additions before being able to produce realistic dynamics in the MLT. However, as a developmental tool the nudging has served its purpose in that it has stabilised the model for this configuration, and has in fact given evidence of producing more complex dynamics, even with its simplified setup. Furthermore, implementation of the correct radiative and chemical forcing will hence provide a complete and stable model of the mesosphere and lower thermosphere up to 100 km.

5 Extending the nudged model

With a stable model produced up to 100 km, we can now extend this new configuration in order to more appropriately model atmospheric waves in the MLT as well model greater atmospheric heights. To this end, we investigate changes to both the vertical resolution and the height of the model upper boundary.

We discuss the implementation of these extensions and analyse the model fields produced with them in place.

5.1 Experimentation with vertical resolution

The standard model vertical level set is set up so that model levels increase in height with each new level introduced, so that at the upper boundary of the 88 level model configuration, the vertical level height is approximately 6 km (cf. Fig. 3). However, this resolution is not appropriate to accurately represent vertically propagating gravity waves in the MLT; these waves are significant carriers of energy and momentum in the upper atmosphere and are important contributors to general circulation in the MLT – see for example Becker & Vadas (2018).

To see why this resolution is not sufficient in the MLT, we must consider the model scale height H. This can be related to the atmospheric temperature so that $H=R_{\mathrm{dry}}T\left(z\right)/g\left(z\right)$, where $R_{\mathrm{dry}}$ is the gas constant for dry air⁴ and $g\left(z\right)$ is the acceleration due to gravity as a function of height above the surface $z$. This is given by $g\left(z\right)=g_0{R}^2/(R_{\oplus }+z{)}^2$, where $g_0$ is the acceleration due to gravity at sea level and $R_{\oplus }$ is the radius of the Earth. This quantity dictates the height for which the atmospheric pressure drops by a factor $e$ at a given temperature.

In order to accurately model vertical wavenumbers in atmospheric dispersion relations – see Griffin & Thuburn (2018) for a more in depth explanation – the vertical resolution must be finer than this scale height. Thus we must make our choice of vertical resolution based on the scale height H so as to be able to resolve these waves. To do this, we plot this scale height as a function of altitude in Figure 16. Since we would like the grid to be finer than the scale height $H$, we also plot $H/2$ and $H/4$ for comparison – these would give two and four vertical levels per scale height respectively. The temperature field used in the calculation of the scale height is the same as that used in the nudging – namely as shown in Figure 8.

Fig. 16 Plot illustrating the scale height of the upper atmosphere H. We require multiple vertical levels per scale height and thus H/2 and H/4 are also plotted to illustrate the vertical level height necessary to have 2 and 4 vertical levels per scale height respectively.

This informs our choice of vertical level heights. With a vertical level height of 1.5 km, we place approximately four vertical levels per scale height. This resolution will give an accurate representation of waves throughout the MLT.

However, this fine resolution leads to concerns regarding consequent numerical instability in the vertical motion as well as computational cost. We therefore also consider a vertical level height of 3 km. This gives approximately two vertical levels per scale height in the MLT and is a good middle ground; some finer vertical wave scales are resolved and the computational cost incurred of resolving at a much finer scale is reduced.

In theory the vertical resolution could increase after around 100 km in keeping with the corresponding increase in scale height. This would make the model less costly and less prone to numerical instability as the model is extended up into the thermosphere. However this is not investigated in this paper and will be the subject of further research on the extended UM.

These two new vertical level sets are shown in Figure 17, and we proceed to test the model using these vertical level sets.

Fig. 17 Vertical level sets for the 3 km maximal vertical spacing (red) and the 1.5 km maximal vertical spacing (blue). These give approximately 2 and 4 vertical levels per scale height in the MLT respectively.

5.1.1 3 km vertical resolution

We first experiment with the vertical level set that has a maximal vertical resolution of 3 km. To do this, we set up the vertical levels so that the model levels increase in depth as before, up until they near 3 km. We then fix the vertical level height at this value. This yields a 94 level configuration going up to a height of 102 km – this can be seen in Figure 17.

With this vertical level set and the nudging scheme employed as in Section 4.1, we observe no effect on the stability of the model and once more have a stable model with the lid at approximately 102 km.⁵ The model is able to run for several years with no issues.

In this paper, we focus on climatological mean fields rather than instantaneous or short term fields used to examine finer atmospheric features such as atmospheric tides. This is because our priority here is model stability, and improving the realism of the UM in the MLT will be the subject of future research. The mean fields provide a good starting point to examine the fields produced by the nudged UM.

Thus, we compare climatologies for monthly mean zonal mean westerly ($u$) wind, vertical ($w$) wind and temperature for the 88 and 94 level configuration. We plot this for December in Figure 12 and for June in Figure 13, with column 2 and column 3 in both figures showing the 88 and 94 level configurations respectively.

In particular, we observe the same general structure in the model fields but with some finer scale vertical features evident towards the upper boundary with the 94 level configuration.

The primary goal here is stability and it is encouraging that no additional model changes such as modification of the sponge layer are necessary for the model to run with the new 94 level configuration. We will see that for a finer vertical resolution, additional modifications are necessary to stabilise the model.

5.1.2 1.5 km vertical resolution

We now experiment with a vertical level set that has a maximal vertical resolution of 1.5 km. The construction is as in Section 5.1.1 but using 1.5 km as the maximal depth – this can be seen in Figure 17. This yields a 113 level configuration going up to a height of around 101 km.

With this configuration, and the nudging implemented as in Section 4.1, the model does not initially run successfully, and crashes due to excessive wind speeds in the North-South and vertical directions.

To combat each of these issues, we make two modifications.

Firstly, we adjust the halo region (cf. Sect. 2.2.2) to account for the larger wind speeds in the lower thermosphere. In order to maximise the wind speeds that the model can cope with, we set the halo region to its maximum permitted size. This will allow for significantly larger wind speeds in the North-South direction. Thus, any further issues with N-S halos must be put down to unphysical winds developing elsewhere in the model.

Secondly we adjust the sponge layer of the model, as described in Section 2.2. To do this we shall increase the value of the vertical damping coefficient $\overline{\mu }$ to damp out large vertical velocities near the model upper boundary until the model is able to run. However, despite this targeted damping, we wish to avoid over damping modelled vertical velocities. Thus, we use as low a value as possible so that the model is consequently stable.

Adding an extended halo gives a stable model run with a vertical damping coefficient of $\overline{\mu }=0.3$ s⁻¹. As in Section 5.1.1, we compare the 88 level climatologies to the corresponding climatologies for the new 113 level configuration. Climatologies for December are shown in Figure 12, and for June in Figure 13, with column 2 and column 4 in both figures showing the 88 and 113 level configurations respectively.

Here, many finer scale vertical and horizontal features become evident from about 40 km upwards where the 1.5 km resolution begins. For example, looking at the westerly ($u$) wind there appears to be more structure in the equatorial winds in the upper stratosphere and mesosphere. This suggests possible impacts on the Semiannual Oscillation (SAO) (e.g. Shepherd et al., 2006), but more detailed study is required to confirm this. In the vertical ($w$) wind, the appearance of an equatorial upwelling is apparent, but this is likely an effect of the upper boundary of the model. It is evident therefore that this finer resolution does indeed give more features in the model fields. By adjusting the halo region and increasing the vertical damping coefficient, we are able to achieve a stable model with this finer vertical resolution.

5.2 Experimentation with the height of the model upper boundary

Now that experiments with resolution have been performed, we look to extend the upper boundary of the model to greater heights to evaluate the performance of the model in preparation for its use as a whole atmosphere model. In the medium term, this will be done by coupling the extended atmospheric model to an appropriate thermospheric model, such as the Drag Temperature Model (DTM) developed by Bruinsma (2015), which spans the 120–1500 km range. This coupling is a goal of the Space Weather Atmosphere Models and Indices (SWAMI) project – see http://swami-h2020.eu/project-swami/. At an absolute minimum, the extended UM will therefore need to run up to 120 km. However, a merging region between the two models would be advantageous. Ideally, we would like to raise the upper boundary as high as possible to give a large merging region between the two models, however the globally uniform temperature structure in the current nudged model means it has a limited ability to represent latitudinal and longitudinal variability, which becomes more significant with increasing altitude in the lower thermosphere. Therefore we present model runs with 120 and 135 km upper boundary heights.

The 1.5 km resolution would be preferable due to its ability to resolve finer scale features in the model fields. However, on extension of the model upper boundary with this resolution, the nudged model could not be run for even short periods without becoming unstable even whilst using the stability mechanisms discussed in Section 2.2. Further additions to the model will be required to allow stable runs with this resolution at higher altitudes, such as more complete radiation and chemistry schemes.

We therefore use the 3 km resolution to extend the model upper boundary. We saw that this resolution still resolves some finer scale features as well as needing no additional stability mechanisms for the 100 km run. In fact, the 3 km resolution will give four vertical levels per scale height above around 120 km, and so this resolution is sufficient to accurately capture waves as the upper boundary is extended higher into the thermosphere.

5.2.1 120 and 135 km model lid

We now extend the 102 km, 94 level configuration to accommodate the new model lid heights.

For the 120 km lid this results in adding six more 3 km levels to give a 100 level configuration going up to exactly 120 km. For the 135 km lid this results in adding eleven more 3 km levels to give a 105 level configuration going up to exactly 135 km.

With these configurations, and the nudging implemented as in Section 4.1, the model does not initially run successfully. This is to be expected as greater wind speeds are observed at greater heights. Thus, as in the 1.5 km resolution setup, we adjust the halo region and sponge layer to combat these issues, and observe the following results:

• 120 km

  • Without the extended halo region, a stable run is achieved with $\overline{\mu }=0.6$ s⁻¹;

  • With the extended halo region, a stable run is achieved with $\overline{\mu }=0.2$ s⁻¹.

• 135 km

  • Without the extended halo region, a stable run cannot be achieved even when $\overline{\mu }$ is adjusted;

  • With the extended halo region, a stable run is achieved with $\overline{\mu }=0.4$ s⁻¹.

It is evident therefore that the extended halos help considerably in preventing the model from crashing due to large instantaneous meridional wind speeds, without degrading the solution. We therefore use the extended halo runs for comparison. We compare the 3 km resolution 100 km lid climatologies presented in Section 5.1.1 (column 1) with the climatologies for both the 120 km (column 2) and 135 km (column 3) lids. Climatologies for December are shown in Figure 18 and for June in Figure 19.

Fig. 18 Zonal mean monthly mean climatologies in December of westerly ($u$) wind (top), vertical ($w$) wind (middle) and temperature (bottom) comparing the nudged 100 km model run with the successful 120 and 135 km lid model runs.

Fig. 19 Zonal mean monthly mean climatologies in June of westerly ($u$) wind (top), vertical ($w$) wind (middle) and temperature (bottom) comparing the nudged 100 km model run with the successful 120 and 135 km lid model runs.

With the modifications to the halo region and sponge layer we are able to achieve stable model runs with an extended upper boundary. In these figures we observe increasing variability from expected fields as the upper boundary increases; however we stress that the focus here is on attaining a stable model. This is an encouraging first step towards a whole atmosphere model, and it remains to add the physical parameterizations necessary to improve the realism of the fields at these heights.

In particular, strong unrealistic easterlies (cf. Figs. 5a and 6a of Liu et al., 2010) and vertical upwelling are evident at the equator in the upper levels of the model. These move upwards and strengthen as the upper boundary is raised. Thus, these features appear to be related to the placement of the upper boundary and as a result of the model specification rather than through realistic mechanisms. We do see the appearance of the westerly summer-time mesospheric jet at around 90 km. However, it is unclear as to whether this is realistic or also an artefact of the placement of the upper boundary.

More realistic damping mechanisms such as molecular viscosity and diffusion have been developed by Griffin & Thuburn (2018). This offers a route to model stability at higher levels which is preferable to ever stronger applications of the model sponge layer. This type of damping becomes important above around 130 km as the molecular diffusion timescale becomes shorter than the time scale for wave growth. As the UM upper boundary is raised to greater heights, it is envisaged that this will replace the current UM sponge layer detailed in Section 2.2.

It is evident therefore that the model still requires several additions before being able to produce realistic dynamics in the MLT. Consequent development of the extended UM will be able to build on the success of the nudging scheme and add the necessary physics to approach realistic simulation of the mesosphere and lower thermosphere.

6 Conclusions and future work

The Met Office’s UM with an extended upper boundary of 100 km was tested and it was discovered that high wind speeds near the poles were causing the model to crash. On investigation of the radiative forcing, it was found that the LTE approximation made in the radiation scheme was leading to unphysical forcing in the MLT. This was the leading factor in causing model instabilities with an extended upper boundary.

A NLTE radiation scheme was not available at the time of performing this research. Thus, a nudging scheme was successfully implemented to replace the radiation scheme above 70 km. This scheme uses relaxation to push the temperature field from the dynamics of the model towards the prescribed temperature profile given in Figure 8. With this in place, the large winds evident from initial model runs disappeared and the extended UM was stabilised with a 100 km upper boundary. This key result provides a valuable developmental tool to further extend the model’s capabilities in the MLT.

To this end, the model vertical resolution and upper boundary height were then scrutinised. On changing the vertical resolution with a fixed 100 km lid, the use of the 3 km resolution provided a stable run with no additional changes necessary, whereas use of the 1.5 km resolution required an increase to the vertical damping coefficient and halo parameter. The model fields produced for both resolutions were in keeping with the original nudged model and additional detail was evident, particularly in the case of the 1.5 km resolution.

On extending the upper boundary however, the 1.5 km resolution made the model unstable. As commented, the 3 km resolution is in fact sufficient to provide four vertical levels per scale height above 120 km in any case, which allows vertically propagating waves to be resolved in the MLT. Thus the 3 km resolution was favoured for experiments to raise the upper boundary height. Successful, stable model runs were performed when the lid was raised to both 120 and 135 km, but an increase in the vertical damping coefficient and halo parameter was necessary. Figure 16 shows that scale height increases with altitude throughout the thermosphere, making the requirement for a 1.5 km or a 3 km vertical resolution too stringent. Future work will include experiments which use coarser vertical resolution in the thermosphere that match the increasing scale height.

Importantly, the various model configurations were all stable with the use of the nudging scheme and adjustment of the various pre-existing stability mechanisms. This is an encouraging first step towards the development of a whole atmosphere UM.

While the model showed some signs of realistic physics, the difference when compared to data was still significant. However, this is to be expected given the simplification of physics used here. In fact, the purpose of the implementation of the nudging scheme was as a developmental tool; to prove that a stable model run was possible with an extended model lid and to provide a base for further model extensions.

It is now the task to further develop the physics, chemistry and dynamics within the model to gradually replace the nudging to climatology used here. With this a more thorough comparison to data can be carried out, looking at more complex atmospheric features such as tides.

In particular, work will be carried out to introduce a NLTE radiation scheme to provide realistic heating in the model (see Jackson et al., 2020, for initial results), as well as introduce molecular viscosity and diffusion to provide a more realistic sponge layer at the upper boundary. This will make good progress into appropriately capturing realistic dynamics in the mesosphere and lower thermosphere in the Met Office’s UM.

Acknowledgments

MJG and CJB are supported by a NERC GW4+ Doctoral Training Partnership studentship from the Natural Environment Research Council [NE/L002434/1] and are thankful for the collaborative support of the Met Office, UK. DJG and DRJ received funding for this work from the European Unions Horizon 2020 research and innovation programme under grant agreement No. 7776287. The editor thanks two anonymous reviewers for their assistance in evaluating this paper.

References

All Tables

All Figures

	Fig. 1 The default sponge layer with vertical damping coefficient $\overline{\mu }=0.05$ s−1 and 100 km model lid. The sponge layer can be seen to be at $\overline{\mu }=0.05$ s−1 for most latitudes and for high altitudes but is doubled at the polar latitudes for all altitudes to help numerical stability in this area.
In the text

	Fig. 2 Schematic illustrating the halo error. On the right hand side of the figure, the advection scheme can be seen trying to access data that is outside of the halo region which causes an error. The halo size can be increased to accommodate the larger winds present in the upper atmosphere.
In the text

	Fig. 3 Vertical level sets for the standard model height (black) and those added for the raised model height (blue). Two levels are added conforming to the 85 km level distribution, whilst the final is placed at 100 km.
In the text

	Fig. 4 A zonal mean profile produced to prescribe the initial configuration of atmospheric ozone for the 100 km model lid. Its inclusion is important in order to represent the correct vertical temperature structure (via radiative heating).
In the text

	Fig. 5 Modelled wind in the 100 km September model run at the time step before the model crash. Stereographic plots at the uppermost model level (a) and (b). Plots of the wind at the southernmost latitude circle (c) and (d). This makes evident the abnormally high wind speeds causing the model to crash.
In the text

Fig. 6 Plots illustrating the issue evident with the SW radiation produced by the radiation scheme. We see abnormal values in exactly the regions where the high wind speeds occur in the previously discussed model runs. (a) Zonal mean of the SW radiative heating for the model run starting in September and crashing in December. (b) Zonal mean of the SW radiative heating for the model run starting in March and crashing in June.

In the text

	Fig. 7 Total heating due to solar energy absorption in the near-IR CO2 (NIR) bands. This illustrates the difference in heating in K day−1 when the Local Thermodynamic Equilibrium (LTE) approximation is made. This approximation breaks down at around 0.1 hPa (approx. 65 km) as can be seen. After Fomichev et al. (2004).
In the text

	Fig. 8 Temperature profile used for the relaxation of the temperature field.
In the text

	Fig. 9 Comparison of the 85 km temperature field for various latitudes and the forcing profile (a), with pertinent area zoomed in (b). Good agreement between the two is seen.
In the text

	Fig. 10 Modelled wind in the 100 km September model run with the nudging scheme in place. Stereographic plots at the uppermost model level (a) and (b). Plots of the wind at the southernmost latitude circle (c) and (d). We see the removal of the anomalous wind speeds present previously by comparison with Figure 5.
In the text

	Fig. 11 SW radiation showing zeroed field above 70 km.
In the text

	Fig. 12 Zonal mean monthly mean climatologies in December of westerly ($u$) wind (top), vertical ($w$) wind (middle) and temperature (bottom) comparing the original 85 km model with the 88, 94 and 113 level configurations.
In the text

	Fig. 13 Zonal mean monthly mean climatologies in June of westerly ($u$) wind (top), vertical ($w$) wind (middle) and temperature (bottom) comparing the original 85 km model with the 88, 94 and 113 level configurations.
In the text

	Fig. 14 Zonal mean monthly mean climatologies in December of westerly ($u$) wind (top) and temperature (bottom). These are compared to URAP data for zonal wind and MLS data for temperature. The approximate height is calculated using a scale height of H = 7200 m.
In the text

	Fig. 15 Zonal mean monthly mean climatologies in June of westerly ($u$) wind (top) and temperature (bottom). These are compared to URAP data for zonal wind and MLS data for temperature. The approximate height is calculated using a scale height of H = 7200 m.
In the text

	Fig. 16 Plot illustrating the scale height of the upper atmosphere H. We require multiple vertical levels per scale height and thus H/2 and H/4 are also plotted to illustrate the vertical level height necessary to have 2 and 4 vertical levels per scale height respectively.
In the text

	Fig. 17 Vertical level sets for the 3 km maximal vertical spacing (red) and the 1.5 km maximal vertical spacing (blue). These give approximately 2 and 4 vertical levels per scale height in the MLT respectively.
In the text

	Fig. 18 Zonal mean monthly mean climatologies in December of westerly ($u$) wind (top), vertical ($w$) wind (middle) and temperature (bottom) comparing the nudged 100 km model run with the successful 120 and 135 km lid model runs.
In the text

	Fig. 19 Zonal mean monthly mean climatologies in June of westerly ($u$) wind (top), vertical ($w$) wind (middle) and temperature (bottom) comparing the nudged 100 km model run with the successful 120 and 135 km lid model runs.
In the text