The effect of atmospheric nudging on the stratospheric residual circulation in chemistry–climate models

We perform the first multi-model intercomparison of the impact of nudged meteorology on the stratospheric residual circulation using hindcast simulations from the Chemistry–Climate Model Initiative (CCMI). We examine simulations over the period 1980–2009 from seven models in which the meteorological fields are nudged towards a reanalysis dataset and compare these with their equivalent free-running simulations and the reanalyses themselves. We show that for the current implementations, nudging meteorology does not constrain the mean strength of the stratospheric residual circulation and that the inter-model spread is similar, or even larger, than in the free-running simulations. The nudged models generally show slightly stronger upwelling in the tropical lower stratosphere compared to the free-running versions and exhibit marked differences compared to the directly estimated residual circulation from the reanalysis dataset they are nudged towards. Downward control calculations applied to the nudged simulations reveal substantial differences between the climatological lower-stratospheric tropical upward mass flux (TUMF) computed from the modelled wave forcing and that calculated directly from the residual circulation. This explicitly shows that nudging decouples the wave forcing and the residual circulation so that the divergence of the angular momentum flux due to the mean motion is not balanced by eddy motions, as would typically be expected in the time mean. Overall, nudging meteorological fields leads to increased inter-model spread for most of the measures of the mean climatological stratospheric residual circulation assessed in this study. In contrast, the nudged simulations show a high degree of consistency in the inter-annual variability in the TUMF in the lower stratosphere, which is primarily related to the contribution to variability from the resolved wave forcing. The more consistent inter-annual variability in TUMF in the nudged models also compares more closely with the variability found in the reanalyses, particularly in boreal winter. We apply a multiple linear regression (MLR) model to separate the drivers of inter-annual and long-term variations in the simulated TUMF; this explains up to ∼75 % of the variance in TUMF in the nudged simulations. The MLR model reveals a statistically significant positive trend in TUMF for most models over the period 1980–2009. The TUMF trend magnitude is generally larger in the nudged models compared to their free-running counterparts, but the intermodel range of trends doubles from around a factor of 2 to a factor of 4 due to nudging. Furthermore, the nudged models generally do not match the TUMF trends in the reanalysis they are nudged towards for trends over different periods in the interval 1980–2009. Hence, we conclude that nudging does not strongly constrain long-term trends simulated by the chemistry–climate model (CCM) in the residual circulation. Our findings show that while nudged simulations may, by construction, produce accurate temperatures and realistic representations of fast horizontal transport, this is not typically the case for the slower zonal mean vertical transport in the stratosphere. Consequently, caution is required when using nudged simulations to interpret the behaviour of stratospheric tracers that are affected by the residual circulation.

1 Introduction

The Brewer–Dobson circulation (BDC) is characterized by upwelling of air in the tropics, poleward flow in the stratosphere, and downwelling at mid-latitudes and high latitudes. The circulation can be separated into two branches: the shallow branch in the lower stratosphere and the deep branch in the middle stratosphere and upper stratosphere (Plumb, 2002; Birner and Bönisch, 2011). The BDC affects the distribution of trace species in the stratosphere, such as ozone, and its strength partly determines the lifetimes of long-lived gases such as chlorofluorocarbons (CFCs; Butchart and Scaife, 2001). It also determines stratosphere-to-troposphere exchange of ozone (Hegglin and Shepherd, 2009), which is important for the tropospheric ozone budget (Wild, 2007). In the tropical lower stratosphere, where the photochemical lifetime of ozone is long, variations and trends in the strength of the BDC are the main drivers of ozone within the annual cycle (Weber et al., 2011) for inter-annual and longer-term variability (Randel and Thompson, 2011) and in response to climate change (e.g. Keeble et al., 2017). Here we focus on the advective part of the BDC, or the residual circulation, which is driven by wave breaking in the stratosphere from planetary-scale Rossby waves and gravity waves (Holton et al., 1995). It is important to note that the overall tracer transport in the stratosphere is also affected by turbulent eddy mixing, which has been evaluated separately in previous studies (Garny et al., 2014; Ploeger et al., 2015a, b; Dietmüller et al., 2018; Eichinger et al., 2019; Šácha et al., 2019). The residual circulation is commonly evaluated in model (Butchart et al., 2010, 2011) and reanalysis (Abalos et al., 2015; Kobayashi and Iwasaki, 2016) studies using the transformed Eulerian mean circulation (TEM; Andrews and McIntyre, 2002, 1978; Andrews et al., 1987).

Past studies have shown substantial spread across chemistry–climate models (CCMs) in the mean strength of the residual circulation (e.g. Butchart et al., 2010). Nevertheless, CCMs consistently simulate a long-term strengthening of the residual circulation with an increase of ∼ 2 % per decade (e.g. Butchart et al., 2010; Hardiman et al., 2014), though there are differences across models in the relative contribution to trends from resolved and parameterized wave forcing. Reanalysis datasets also suggest a strengthening of the residual circulation over the past few decades of the order of 2 % per decade to 5 % per decade (Abalos et al., 2015; Miyazaki et al., 2016), apart from one dataset (ERA-Interim – ERA-I) which shows a weakening of the deep branch of the BDC (Seviour et al., 2012; Abalos et al., 2015). However, reanalyses are subject to multiple caveats, particularly in their suitability for trend studies, and there can be substantial differences in residual circulation trends calculated from the same reanalysis using different methods (Abalos et al., 2015).

Given the limitations of reanalyses, evaluating the fidelity of model estimates of residual circulation variability and trends is challenging, since there are no direct measurements of the residual circulation. The only direct estimates of the stratospheric mass circulation come from tracer measurements, which can be used to calculate the stratospheric age of air (AoA; Kida, 1983; Schmidt and Khedim, 1991; Waugh and Hall, 2002). The AoA represents the combined effects of advection and mixing processes and as such cannot be directly related to the residual circulation. While progress has been made in separating the relative effects of advection and mixing for the AoA calculated from models (Garny et al., 2014; Dietmüller et al., 2018; Eichinger et al., 2019; Šácha et al., 2019) from Lagrangian models driven by reanalysis data (Ploeger et al., 2015a, b, 2019; Ploeger and Birner, 2016), and comparing the effects in both CCMs and Lagrangian models (Dietmüller et al., 2017), this is more difficult to achieve in observations. Engel et al. (2009) used balloon-borne measurements of stratospheric trace gases and found a statistically non-significant increase in the AoA in the middle stratosphere at northern mid-latitudes; this has been corroborated in a more recent study using longer measurement records at two mid-latitude sites in the Northern Hemisphere (NH; Engel et al., 2017). It has been hypothesized based on analyses of recent satellite tracer datasets, which have greater spatial and temporal coverage, that subtropical AoA trends can be explained by a weakening of the mixing barriers at the edge of the tropical pipe (Neu and Plumb, 1999) that is masking the effects of an increase in tropical upwelling on the AoA (Stiller et al., 2012; Haenel et al., 2015). In contrast with AoA trends derived from observations, CCMs forced with observed sea-surface temperatures (SSTs), greenhouse gases, and ozone-depleting substances (ODSs) show a decrease in the AoA throughout the stratosphere (Karpechko and Maycock, 2018; Li et al., 2018; Morgenstern et al., 2018; Abalos et al., 2019; Polvani et al., 2019). Theoretical approaches based on the tropical leaky pipe model (Neu and Plumb, 1999) have shown promise for bridging the information on the stratospheric circulation derived from observations with outputs from general circulation models (GCMs) and CCMs (Ray et al., 2016), but differences remain (Karpechko and Maycock, 2018).

More recent theoretical developments offer a means of calculating the diabatic circulation using stratospheric tracers (Linz et al., 2017), which is a promising avenue, as this is more closely related to the residual circulation than the AoA. Linz et al. (2017) showed consistent estimates of the diabatic circulation in the lower stratosphere based on two independent satellite tracer datasets but identified large uncertainties of up to a factor of 2 in the mean circulation strength in the upper stratosphere. Hence, the available tracer datasets are not yet suitable for characterizing trends in the diabatic circulation using these methods. Targeted measurement strategies to better characterize long-term changes in the stratospheric meridional circulation have been proposed (Moore et al., 2014; Ray et al., 2016).

In an attempt to obtain a closer comparison with observed stratospheric trace species, some studies have used model simulations with meteorological fields nudged or relaxed towards analysis or reanalysis datasets (Jeuken et al., 1996). These include studies of stratospheric ozone variability and trends (e.g. van Aalst et al., 2004; Solomon et al., 2016; Hardiman et al., 2017b; Ball et al., 2018) and comparisons between models and satellite-based multi-species observational records (Froidevaux et al., 2019), in particular focusing on specific meteorological events such as the sudden stratospheric warming in the 2009–2010 winter (Akiyoshi et al., 2016) as well as the chemical and climatic effects of volcanic eruptions (Löffler et al., 2016; Solomon et al., 2016; Schmidt et al., 2018). Nudged simulations have also been used to study mechanisms for dynamical coupling between the stratosphere and troposphere (Hitchcock and Simpson, 2014) and to examine the effects of different regions on atmospheric predictability (e.g. Douville, 2009; Jung et al., 2010). Nudging involves adding additional tendencies to the model equations to constrain the modelled variables. Nudged variables can include horizontal winds (or divergence and vorticity), temperature, surface pressure, and latent and sensible heat fluxes. However, vertical winds, which are a small residual from horizontal divergence, are not nudged, and the underlying model physics can yield quite different results from the datasets they are nudged towards (Telford et al., 2008; Hardiman et al., 2017b).

The approach of nudging a CCM towards reanalysis data follows a similar philosophy to traditional offline chemical transport models (CTMs), though there are fundamental differences between these types of models in terms of their tracer advection. CTMs need to match the mass transport with the evolution of the pressure field. This can be done exactly in isobaric coordinates (often used in the stratosphere) but requires a correction in regions where grid box mass changes (e.g. as surface pressure changes). CCMs are less affected by this mass–wind inconsistency than CTMs (Jöckel et al., 2001), but nudging will add forcings that are inconsistent with the model state. CTMs use the full 3-D circulation from the analyses and reanalyses directly and have been widely developed and used over the past few decades (e.g. Rood et al., 1989; Chipperfield et al., 1994; Lefèvre et al., 1994). They have proven to be very successful at simulating stratospheric tracers on a range of timescales (Chipperfield, 1999), including decadal changes (Mahieu et al., 2014). However, this success has been built on extensive testing of the optimum way to use the reanalysis data to force the CTMs. For example, Chipperfield (2006) showed how different approaches to calculating the vertical velocity in the TOMCAT/SLIMCAT model could lead to very different distributions of stratospheric age of air, while Monge-Sanz et al. (2013a) compared the performance of different European Centre for Medium-Range Weather Forecasts (ECMWF) analyses within the same CTM framework. Krol et al. (2018) recently provided a summary of how current CTMs intercompare for tracer calculations. Monge-Sanz et al. (2013b) compared the approaches of using ECMWF analyses directly in a CTM with the ECMWF CCM nudged using the same analyses. They found that the CTM and nudged CCM were consistent in showing a degraded performance when using older ERA-40 reanalysis compared to the later ERA-Interim. However, they also showed some differences between CTM and nudged-CCM tracers using the same analyses, with the nudged CCM showing stronger upward motion in the tropical stratosphere. Therefore, with regards to the slow residual circulation, one cannot assume that a nudged CCM will behave in a similar way to a CTM even when using the same meteorological analyses. Recently, Ball et al. (2018) showed two nudged CCMs which failed to capture the observed variations in the lower-stratospheric ozone as measured by satellite observations, while Chipperfield et al. (2018), using the TOMCAT CTM, simulated a better agreement of modelled ozone variations with the observations. Overall, the success of some CTM simulations in simulating long-lived stratospheric tracers has been built on many years of model development and testing. In contrast, nudged CCMs are much newer tools and have not yet been evaluated to the same extent. A recent study by Orbe et al. (2018) analysed tropospheric tracers in nudged-CCM simulations and found large differences in the distributions of the tracers, which could be partly traced to differences in the model convection schemes. They urged users to adopt a cautious approach when interpreting tracers in nudged simulations given their dependence not only on large-scale flow but also on sub-grid parameterizations. However, a critical evaluation of the stratospheric residual circulation in nudged-CCM simulations has been lacking to date.

To examine the effect of nudging on the stratospheric residual circulation this study compares hindcast simulations from free-running and nudged versions of the same models that participated in the phase 1 of the Chemistry–Climate Model Initiative (CCMI; Morgenstern et al., 2017). Nudged experiments were not performed in previous chemistry–climate multi-model comparisons (Chemistry–Climate Model Validation Activity 2; CCMVal-2), so CCMI offers a timely opportunity to evaluate the effect of nudging on simulated mean biases, variability, and long-term trends in the residual circulation. For completeness, we also present a comparison between the nudged simulations and the reanalysis datasets the models are nudged towards. The paper is laid out as follows. Section 2 describes the CCMI and the reanalysis data used in the present study along with the diagnostics for the residual circulation; Sect. 3 presents results covering the mean circulation, annual cycle, inter-annual variability, and trends; and Sect. 4 summarizes the results and discusses the implications for using nudged simulations to study aspects of the observational record.

2 Data and methods

2.1 Models and experiments

CCMI is the successor activity to CCMVal-2 and the Atmospheric Chemistry and Climate Model Intercomparison Project (ACCMIP; Lamarque et al., 2013). We use the hindcast free-running simulations, REF-C1, and the nudged specified dynamics simulations, REF-C1SD, which cover the periods 1960–2009 and 1980–2009, respectively. Here we analyse the common 30-year period 1980–2009 that was run by all models for both experiments with prescribed observed SSTs and sea ice concentrations. The CCMI data were downloaded from the British Atmospheric Data Centre (Hegglin and Lamarque, 2015). For an extensive overview of the CCMI models, see Morgenstern et al. (2017). We analyse those CCMI models that output the necessary TEM diagnostics. At a minimum, this requires the residual vertical velocity (${\stackrel{\mathrm{‾}}{w}}^{\ast }$) and the residual meridional velocity (${\stackrel{\mathrm{‾}}{v}}^{\ast }$; Andrews et al., 1987); where available we also use the resolved and parameterized wave forcing fields from the models. This gives results from a total of 10 models, which differ from one another in various aspects, such as their horizontal resolution, ranging from 1.9 to 5.6^∘, their vertical resolution, and their sub-grid parameterizations (see Table 1). The main text concentrates on the 7 out of 10 models that performed both the REF-C1 and REF-C1SD experiments (Table 2). However, the broad conclusions drawn in the main text for the characteristics of the seven-member REF-C1 ensemble are consistent with the behaviour for all 10 models. Hence the three models that only performed the REF-C1 experiment (GEOSCCM, NIWA-UKCA, and ULAQ-CCM) are not discussed further, but for completeness a subset of diagnostics from those models is shown in the Supplement (Figs. S1–S5).

Table 1 CCMI models that provided TEM diagnostic model output used in this study. CP is Charney–Phillips, T21 ≈ 5.6^∘ × 5.6^∘, T42 ≈ 2.8^∘ × 2.8^∘, T47 ≈ 2.5^∘ × 2.5^∘, TL159 ≈ 1.125^∘ × 1.125^∘, TA is hybrid terrain-following altitude, TP is hybrid terrain-following pressure, and NTP is non-terrain-following pressure.

For the REF-C1 simulations we analyse between one and five ensemble members (depending on what was available), and for REF-C1SD the one realization submitted from each model. The REF-C1SD simulations nudge temperature and other meteorological fields such as horizontal winds, vorticity and divergence, and some surface fields (Table 2), while the chemical fields are left to evolve freely. The nudging timescales range from 6 to 50 h, and the height range over which nudging is applied varies (Table 2). The TEM and related diagnostics that were available from each model are shown in Table 3. The models use different reanalysis fields for nudging taken from ERA-Interim (Dee et al., 2011), JRA-55 (Ebita et al., 2011; Kobayashi et al., 2015), or MERRA (Rienecker et al., 2011). The differences in the residual circulation diagnosed from reanalyses have been identified and documented in previous studies (e.g. Abalos et al., 2015).

Table 2 Details of nudging in the CCMI REF-C1SD simulations that provided TEM diagnostics model output used in this study. ERA-I is ERA-Interim, CIRA is Cooperative Institute for Research in the Atmosphere, MERRA is Modern-Era Retrospective reanalysis, and JRA-55 is Japanese 55-year Reanalysis. T (with wave 0) for EMAC refers to the additional nudging of the global mean temperature.

Table 3 Available TEM-related model output for each model from the CCMI-1 archive: ${\stackrel{\mathrm{‾}}{w}}^{\ast }$ (✓), ${\stackrel{\mathrm{‾}}{v}}^{\ast }$ (•), EPFD (+), gravity wave drag (OGWD and NOGWD; ⋆), OGWD (▴), and NOGWD (▪).

2.2 Model diagnostics

2.2.1 TEM residual circulation

The TEM velocities (${\stackrel{\mathrm{‾}}{v}}^{\ast },{\stackrel{\mathrm{‾}}{w}}^{\ast }$) are defined as (Andrews et al., 1987)

$$\begin{array}{}\text{(1)}& {\stackrel{\mathrm{‾}}{v}}^{\ast }={\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}\cdot a\cdot \cos \mathit{\varphi }}}{\displaystyle \frac{\partial {\stackrel{\mathrm{‾}}{\mathrm{\Psi }}}^{\ast }}{\partial z}}{\stackrel{\mathrm{‾}}{w}}^{\ast }={\displaystyle \frac{\mathrm{1}}{{\mathit{\rho }}_{\mathrm{0}}\cdot a\cdot \cos \mathit{\varphi }}}{\displaystyle \frac{\partial {\stackrel{\mathrm{‾}}{\mathrm{\Psi }}}^{\ast }}{\partial \mathit{\varphi }}},\end{array}$$

where ${\stackrel{\mathrm{‾}}{\mathrm{\Psi }}}^{\ast }(\mathit{\phi },z)$ is the residual meridional mass streamfunction, ρ₀ is log-pressure density, α is Earth's radius, and ϕ is latitude. As most of the models analysed here use a hybrid-pressure vertical coordinate, the prognostic variable is the pressure vertical velocity ${\stackrel{\mathrm{‾}}{\mathit{\omega }}}^{\ast }$ (calculated in Pa s⁻¹), which must be converted to metres per second in order to get the residual vertical velocity, ${\stackrel{\mathrm{‾}}{w}}^{\ast }$. The conversion of ω to w in isobaric coordinates is given by the following equation:

$$\begin{array}{}\text{(2)}& \mathit{\omega }={\displaystyle \frac{\mathrm{d}p}{\mathrm{d}t}}={\displaystyle \frac{\partial z}{\partial t}}{\displaystyle \frac{\partial p}{\partial z}}=w{\displaystyle \frac{-pg}{RT}}=w{\displaystyle \frac{-p}{H}},\end{array}$$

where p is pressure, R=287 J K⁻¹ kg⁻¹ is the gas constant for dry air, and H is a fixed scale height. Both TEM velocity components were submitted as monthly mean fields to the CCMI data archive. Upon close examination of the CCMI model output, some discrepancies were found in the way that the residual vertical velocity was calculated among the models. Although a fixed scale height of H=6950 m was recommended in the CCMI data request (Eyring et al., 2013; Hegglin and Lamarque, 2015), the TEM output from some models (EMAC and SOCOL3) was calculated incorrectly using a temperature-dependent density, ${\mathit{\rho }}_{\mathrm{0}}=p/RT$, instead of the log-pressure definition of the density, ${\mathit{\rho }}_{\mathrm{0}}={\mathit{\rho }}_s\cdot e^{-z/H}$, such that z has a unique 1:1 correspondence with p. This methodological error leads to artificial spread in the model ${\stackrel{\mathrm{‾}}{w}}^{\ast }$ fields. We note that previous multi-model comparisons of the residual circulation that use ${\stackrel{\mathrm{‾}}{w}}^{\ast }$ taken directly from models may have been subject to the same issue (e.g. Butchart et al., 2010; SPARC, 2010), though we cannot confirm this. To avoid this methodological inconsistency, Dietmüller et al. (2018) recalculated ${\stackrel{\mathrm{‾}}{w}}^{\ast }$ from ${\stackrel{\mathrm{‾}}{v}}^{\ast }$ using the continuity equation, which requires a vertical integration and a derivative along the meridional direction. The recalculation of ${\stackrel{\mathrm{‾}}{w}}^{\ast }$ from ${\stackrel{\mathrm{‾}}{v}}^{\ast }$ was also explored for this study, but it was found to introduce additional errors affecting the latitudinal structure of ${\stackrel{\mathrm{‾}}{w}}^{\ast }$ (not shown) specifically because of the reduced number of CCMI-requested pressure levels compared to the native model levels. We were able to overcome the discrepancy in the submitted ${\stackrel{\mathrm{‾}}{w}}^{\ast }$ fields for the EMAC simulations by reconverting high-frequency ${\stackrel{\mathrm{‾}}{\mathit{\omega }}}^{\ast }$ output to ${\stackrel{\mathrm{‾}}{w}}^{\ast }$ using the log-pressure density as in Eq. (2). However, for SOCOL3 the required output for this was not available, and hence we use the submitted ${\stackrel{\mathrm{‾}}{w}}^{\ast }$ for which the absolute values should be treated with caution. For the other models, the results presented in this study are based on the original diagnostics submitted to the CCMI data archive, which we have verified were calculated in the correct way.

We compute the mass flux across a given pressure surface as (Rosenlof, 1995)

$$\begin{array}{}\text{(3)}& \mathrm{2}\mathit{\pi }\underset{\mathit{\varphi }}{\overset{\mathrm{pole}}{\int }}{\mathit{\rho }}_{\mathrm{0}}{a}^{\mathrm{2}}\cos \mathit{\varphi }{\stackrel{\mathrm{‾}}{w}}^{\ast }\mathrm{d}\mathit{\varphi }=\mathrm{2}\mathit{\pi }a{\mathrm{\Psi }}^{\ast }\left(\mathit{\varphi }\right),\end{array}$$

using the boundary condition in which Ψ=0 is at the poles. By finding, at each pressure level, the latitude at which Ψ_max and Ψ_min occur, which corresponds to the height-dependent turnaround (TA) latitudes, we can calculate the net downward mass flux in each hemisphere. The net tropical upward mass flux, equal to the sum of the downward mass fluxes in each hemisphere, can then be expressed as (Rosenlof, 1995)

$$\begin{array}{}\text{(4)}& \begin{array}{rl}& \mathrm{tropical}\mathrm{upward}\mathrm{mass}\mathrm{flux}\left(\mathrm{TUMF}\right)=\\ & \mathrm{2}\mathit{\pi }a({\mathrm{\Psi }}_{max}^{\ast }-{\mathrm{\Psi }}_{min}^{\ast }).\end{array}\end{array}$$

The tropical upward mass flux (TUMF) has been used widely as a measure of the strength of the BDC (e.g. Rosenlof, 1995; Butchart and Scaife, 2001; Butchart et al., 2006, 2010, 2011; Butchart, 2014, and references therein; Seviour et al., 2012), so its use here enables a direct comparison with earlier studies. Arguably, the strength of the TUMF is a first-order metric for evaluating changes in the stratospheric mass circulation as a consequence of nudging. As mentioned above, by calculating the annual means of TUMF accounting for the seasonal cycle of the TA latitudes, we capture the correct evolution of the intra-seasonal (not shown) and inter-annual variability in the TUMF.

2.2.2 Downward control principle calculations

Under steady-state conditions, ${\stackrel{\mathrm{‾}}{\mathrm{\Psi }}}^{\ast }(\mathit{\phi },z)$ at a specified latitude φ and log-pressure height z is given by the vertically integrated eddy-induced total zonal force, $\stackrel{\mathrm{‾}}{F}$, above that level (Haynes et al., 1991):

$$\begin{array}{}\text{(5)}& {\mathrm{\Psi }}^{\ast }(\mathit{\phi },z)=\underset{z}{\overset{\mathrm{\infty }}{\int }}{\left\{{\displaystyle \frac{{\mathit{\rho }}_{\mathrm{0}}{a}^{\mathrm{2}}\stackrel{\mathrm{‾}}{F}{\cos }^{\mathrm{2}}\mathit{\varphi }}{{\stackrel{\mathrm{‾}}{m}}_{\mathit{\varphi }}}}\right\}}_{\mathit{\varphi }=\mathit{\varphi }\left(z^{\prime }\right)}\mathrm{d}z,\end{array}$$

where in the quasi-geostrophic limit, ${\stackrel{\mathrm{‾}}{m}}_{\mathit{\varphi }}\approx -\mathrm{2}\mathrm{\Omega }{a}^{\mathrm{2}}\sin \mathit{\varphi }\cos \mathit{\varphi }$. The above integration applies along lines of constant mean absolute angular momentum per unit mass, $\stackrel{\mathrm{‾}}{m}=a\cos \mathit{\varphi }(\stackrel{\mathrm{‾}}{u}+a\mathrm{\Omega }\cos \mathit{\varphi })$, where $\stackrel{\mathrm{‾}}{u}$ is the zonal mean zonal wind and Ω is Earth's rotation rate, with boundary conditions of Ψ→0 and ${\mathit{\rho }}_{\mathrm{0}}{\stackrel{\mathrm{‾}}{w}}^{\ast }\to$ as z→∞. These lines of constant angular momentum are approximately vertical except near the Equator (up to $\sim \pm \mathrm{20}$^∘) such that we can calculate the solution of the above integral using the constant φ for the limits of the integral (Haynes et al., 1991). In climate models, $\stackrel{\mathrm{‾}}{F}$ has contributions from resolved waves due to the Eliassen–Palm flux divergence (EPFD) and from parameterized gravity wave drag due to sub-grid-scale waves that originate from orography, convection and frontal instabilities. This enables us to estimate the contribution to the tropical upward mass flux of both resolved planetary wave driving (EPFD) and the orographic and non-orographic parameterized gravity wave drag (OGWD and NOGWD, respectively) from the CCMI model output (Table 3) and compare with the direct estimates derived from ${\stackrel{\mathrm{‾}}{w}}^{\ast }$.

Applying the downward control principle (Haynes et al., 1991) can provide useful insights into the driving mechanisms of the stratospheric residual circulation and therefore explain part of the inter-model spread found in both REF-C1 and REF-C1SD simulations. While the downward control principle enables the contributions of EPFD and OGWD and NOGWD to TUMF to be calculated under various assumptions (Haynes et al., 1991), one has to keep in mind that the different wave forcings can interact and thus are not independent of each other (Cohen et al., 2013).

It is important to note some possible limitations of the diagnostic approaches chosen for this study. Both the direct and downward control principle methods rely on the applicability of quasi-geostrophic theory to interpret the results. In addition to the two approaches used here, the residual circulation can also be estimated using the thermodynamic equation. Studies have shown that the estimates from the different methods for evaluating the residual circulation can differ (Seviour et al., 2012; Abalos et al., 2015; Linz et al., 2019), particularly in reanalyses where standard global conservation laws (e.g. conservation of mass) must generally not be met. Similar issues are likely to beset the nudged model simulations, owing to the additional tendencies included in the model equations. The differences between the calculation methods for the residual circulation can be as large as, or larger than, the differences between reanalysis datasets for the same diagnostic (Abalos et al., 2015; Linz et al., 2019) and may further depend on choices around averaging between fixed latitudes or the TA latitudes (Linz et al., 2019), so it is important to bear this in mind in interpretation of the results presented here. Unfortunately, heating rates were not available from all CCMI model simulations to perform the thermodynamic equation calculation. Nevertheless, we compute the direct and downward control principle diagnostics for the residual circulation in a self-consistent manner in the models and reanalyses to enable comparison with earlier multi-model studies (Butchart et al., 2006, 2010; SPARC, 2010).

2.3 Multiple linear regression model

To investigate the drivers of inter-annual variability in the residual circulation we apply a multiple linear regression (MLR) model (Eq. 6) to the annual mean time series of TUMF. The model includes terms for known drivers of variations in tropical lower-stratospheric upwelling: major volcanic eruptions (Pitari and Rizi, 1993), the El Niño–Southern Oscillation (ENSO; García-Herrera et al., 2006; Marsh and Garcia, 2007; Randel et al., 2009), the quasi-biennial oscillation (QBO; Baldwin et al., 2001), and a linear trend (Calvo et al., 2010):

$$\begin{array}{}\text{(6)}& \begin{array}{rl}\mathrm{TUMF}\left(t\right)& ={\mathit{\beta }}_{\mathrm{0}}+{\mathit{\beta }}_{\mathrm{VOL}}\cdot x_{\mathrm{VOLC}}\left(t\right)+{\mathit{\beta }}_{\mathrm{ENSO}}\cdot x_{\mathrm{ENSO}}\left(t\right)\\ & +{\mathit{\beta }}_{\mathrm{TREND}}\cdot x_{\mathrm{TREND}}\left(t\right)+{\mathit{\beta }}_{\mathrm{QBO}\mathrm{1}}\cdot x_{\mathrm{QBO}\mathrm{1}}\left(t\right)\\ & +{\mathit{\beta }}_{\mathrm{QBO}\mathrm{2}}\cdot x_{\mathrm{QBO}\mathrm{2}}\left(t\right)+\mathit{\epsilon }\left(t\right),\end{array}\end{array}$$

where β₀ is a constant, β_i is the regression coefficient for basis function x_i, and ε(t) is the residual. Following Maycock et al. (2018), the volcanic basis function is defined as the tropical lower-stratospheric average volcanic surface-area density (SAD), the ENSO basis function is the time series of eastern–central equatorial Pacific Ocean SST anomalies (Niño 3.4 index; 5^∘ S to 5^∘ N; 170 to 120^∘ W); the two QBO terms are the first two principal-component time series from an empirical orthogonal function (EOF) analysis of the zonal mean zonal winds between 10^∘ S and 10^∘ N and 70 to 5 hPa and a linear trend. The first three regressors, the volcanic, ENSO, and the linear trend, are identical for both REF-C1 and REF-C1SD runs, while the QBO terms are calculated using the model winds for each model and experiment. For the REF-C1 runs, CMAM does not include a QBO; hence when we apply the MLR to the CMAM REF-C1 simulation the QBO terms are omitted. We opted not to include an equivalent effective stratospheric chlorine (EESC) MLR term to account for changes in ozone-depleting substances (Abalos et al., 2019; Morgenstern et al., 2018; Polvani et al., 2018, 2019), as the period considered in the study may not be sufficiently long for the linear trend to be separated properly from EESC. Since we are regressing annual mean TUMF we do not consider a seasonal cycle term or any lag in the terms. The results in Sect. 3.5 focus on the first ensemble member (in the rip-nomenclature, where r stands for realization, i for initialization, and p for physics – r1i1p1), but where applicable the results from the MLR model for the rest of the ensemble members of the REF-C1 runs are presented in the Supplement (Figs. S6–S9).

2.4 Reanalysis data

In order to compare the REF-C1 and REF-C1SD simulations against the reanalysis datasets used for the nudging, we use the SPARC Reanalysis Intercomparison Project (S-RIP) dataset (Martineau, 2017; Martineau et al., 2018). This provides a common gridded version of the reanalysis TEM fields on a 2.5^∘ × 2.5^∘ grid up to 1 hPa. The pressure vertical velocity, ${\stackrel{\mathrm{‾}}{\mathit{\omega }}}^{\ast }$, is converted to the residual vertical velocity, ${\stackrel{\mathrm{‾}}{w}}^{\ast }$, using Eq. (2). A detailed comparison of the stratospheric residual circulation in reanalysis datasets is given by Abalos et al. (2015).

3 Results

3.1 Climatological residual circulation: ${\stackrel{\mathrm{‾}}{w}}^{\ast }$

Figure 1 shows latitude–pressure cross-sections of the climatological (1980–2009) multi-model mean (MMM) annual mean ${\stackrel{\mathrm{‾}}{w}}^{\ast }$ for the REF-C1 (Fig. 1a) and REF-C1SD (Fig. 1b) simulations and their differences (Fig. 1c). In Fig. 1c absolute differences are computed, so positive values indicate where the magnitude of the circulation in REF-C1SD (whether upwelling or downwelling) is larger than in REF-C1. As expected, the climatologies show upwelling in the tropics between around 30^∘ S to 30^∘ N and downwelling at higher latitudes. In the lowermost stratosphere (100–80 hPa), within the region of tropical upwelling, the REF-C1SD MMM generally shows larger ${\stackrel{\mathrm{‾}}{w}}^{\ast }$ values in the subtropics and smaller values at the Equator compared to REF-C1, indicating a tendency for a more double-peaked ${\stackrel{\mathrm{‾}}{w}}^{\ast }$ structure in the tropics in the lowermost stratosphere (Ming et al., 2016a). Above this, between ∼ 70 and 4 hPa, the REF-C1SD MMM shows on average stronger upwelling at the Equator compared to REF-C1, indicating a less pronounced double-peaked ${\stackrel{\mathrm{‾}}{w}}^{\ast }$ structure in the REF-C1SD experiments in the lower stratosphere to middle stratosphere. Between 1 and 2 hPa, the REF-C1SD MMM shows larger inter-hemispheric asymmetry than in REF-C1, with stronger upwelling in the northern tropics (Fig. 1b). At the southern mid-latitudes, between ∼ 30 and 60^∘ S, the REF-C1SD MMM exhibits on average slightly weaker downwelling than in REF-C1, with the largest magnitude differences found in the upper stratosphere. In the Arctic, the REF-C1SD MMM shows significantly stronger downwelling over the poles than in REF-C1. In the Antarctic the picture is more complex, with the REF-C1SD MMM showing weaker downwelling right at the pole in the upper stratosphere (2–10 hPa) but stronger downwelling between around 75 and 88^∘ S. In the middle stratosphere, from 50 to 10 hPa, the REF-C1SD MMM shows stronger downwelling between 60 and 80^∘ S.

Figure 1 Latitude vs. pressure climatology (1980–2009) of MMM annual mean ${\stackrel{\mathrm{‾}}{w}}^{\ast }$ for (a) REF-C1 simulations, (b) REF-C1SD simulations, and (c) the REF-C1SD–REF-C1 absolute differences. Shading denotes statistical significance at the 95 % confidence level, and the red lines in (c) denote the climatological turnaround latitudes in REF-C1SD.

To show the differences in the transition between regions of upwelling and downwelling motion, Fig. 2 shows vertical profiles of the climatological annual mean turnaround (TA) latitudes in each hemisphere for the REF-C1 and REF-C1SD MMM and the three reanalysis datasets used for nudging. Note that since five of the REF-C1SD models were nudged towards ERA-I, the REF-C1SD MMM may be more weighted towards ERA-I than the other reanalyses. In the NH, the REF-C1SD MMM shows a more poleward TA latitude compared to both REF-C1 and the reanalyses throughout almost the entire depth of the stratosphere (Fig. 2b). A more poleward TA latitude for REF-C1SD than in both REF-C1 and the reanalyses is also found in the Southern Hemisphere (SH) at pressures greater than 30 hPa (Fig. 2a). Hence the nudged simulations show, on average, a wider region of tropical upwelling in the lower stratosphere compared to their free-running counterparts by up to around 5^∘ latitude. In the middle stratosphere and upper stratosphere the REF-C1SD MMM shows a narrower upwelling region in the SH. Interestingly, above 10 hPa in the SH (Fig. 2a), the REF-C1SD does not show a progressive widening of the upwelling region with decreasing pressure as seen in the reanalyses. This is reflected in the structural differences in ${\stackrel{\mathrm{‾}}{w}}^{\ast }$ in the SH upper stratosphere found in some models (Fig. S10). It should be noted though that the differences in TA latitudes between the REF-C1 and REF-C1SD MMMs are comparable to the differences found between the three reanalysis datasets.

Figure 2 Vertical profiles of the climatological turnaround latitudes in the stratosphere for the MMM of the REF-C1 runs (black dashes), the MMM of the REF-C1SD runs (grey dashes), and the S-RIP reanalysis datasets (ERA-I, JRA-55, and MERRA) for the (a) Southern Hemisphere and (b) Northern Hemisphere.

Focusing on the lower stratosphere, Fig. 3 shows the climatological annual mean ${\stackrel{\mathrm{‾}}{w}}^{\ast }$ at 70 hPa in the individual models for the (a) REF-C1 and (b) REF-C1SD simulations and (c) their differences. Also plotted in Fig. 3b is ${\stackrel{\mathrm{‾}}{w}}^{\ast }$ from the reanalyses, and Fig. 3d shows the difference between each REF-C1SD simulation and the reanalysis they were nudged towards. Within the upwelling region, all the models show a clear double-peaked ${\stackrel{\mathrm{‾}}{w}}^{\ast }$ structure in the tropics, with the exception of the CCSRNIES-MIROC3.2 and MRI-ESM1r1 models in the REF-C1SD experiment. In those two cases, CCSRNIES-MIROC3.2 simulates a tri-modal ${\stackrel{\mathrm{‾}}{w}}^{\ast }$ structure, while MRI-ESMr1 shows a relatively constant ${\stackrel{\mathrm{‾}}{w}}^{\ast }$ across the tropics. For the REF-C1 experiment, both EMAC simulations, CMAM and SOCOL3, show a narrower double-peaked structure, with EMAC-L47 exhibiting a rather pronounced NH subtropical maximum. Conversely, CESM1-WACCM simulates the broadest region of tropical upwelling in the lower stratosphere, with the SH subtropical maximum occurring at higher latitudes compared with the rest of the models. The other REF-C1 simulations also exhibit a double-peaked ${\stackrel{\mathrm{‾}}{w}}^{\ast }$ structure, which is generally more hemispherically symmetric, but with varying amplitudes.

Figure 3 Mean strength of annual mean ${\stackrel{\mathrm{‾}}{w}}^{\ast }$ (mm s⁻¹) at 70 hPa for (a) REF-C1 free-running models, (b) REF-C1SD nudged models, (c) absolute differences between the REF-C1SD and REF-C1 experiment for each model, and (d) absolute differences between each REF-C1SD simulation and the respective reanalysis used for nudging.

A double-peaked ${\stackrel{\mathrm{‾}}{w}}^{\ast }$ structure in the lower stratosphere has previously been shown in reanalysis datasets (Abalos et al., 2015; Ming et al., 2016a) and some CCMs (Butchart et al., 2006, 2010). This can also be seen in Fig. 3b for the three reanalysis datasets (ERA-I, JRA-55, and MERRA), where ERA-I and JRA-55 show an asymmetric double-peaked structure with stronger upwelling in the NH subtropics compared to the SH. As documented by Abalos et al. (2015), based on the direct calculation of the residual circulation, MERRA exhibits downwelling at the Equator, an issue which was highlighted in Abalos et al. (2015) and manifested as a negative cell in the streamfunction.

Figure 3c shows the absolute differences in ${\stackrel{\mathrm{‾}}{w}}^{\ast }$ at 70 hPa between the REF-C1SD and REF-C1 experiments. Positive values show where the magnitude of the circulation in REF-C1SD is larger than in REF-C1. The largest differences are generally found within the inner tropics, where CCSRNIES-MIROC3.2, CMAM, and MRI-ESM1r1 exhibit significantly stronger upwelling (up to 3 times more for CMAM) near the local ${\stackrel{\mathrm{‾}}{w}}^{\ast }$ minimum at the Equator. There are also larger differences in many models near edges of the upwelling region (30–40^∘ S), which reflect differences in the width of the tropical pipe between the free-running and nudged simulations (Fig. 2 and Sect. 3.3). Around the subpolar and polar latitudes of the SH, the majority of the REF-C1 models simulate stronger downwelling than their nudged counterparts, while in the NH extratropics no consistent picture emerges across the models. EMAC-L47 and EMAC-L90 show markedly different behaviours despite the fact they are nudged towards the same reanalysis (ERA-I) and differ only in their vertical resolution. This indicates that the effect of nudging on the mean residual circulation is likely to be sensitive to a great number of factors that vary from model to model.

Another interesting result from Fig. 3 is that the inter-model spread in ${\stackrel{\mathrm{‾}}{w}}^{\ast }$ for both experiments is larger in the NH downwelling region than in the equivalent region of the SH. Specifically, the inter-model spread is 0.14 mm s⁻¹ for the REF-C1 runs for all points between 30 and 80^∘ S and 0.2 mm s⁻¹ for points between 30 and 80^∘ N, while for REF-C1SD the values are 0.12 and 0.19 mm s⁻¹, respectively. This also demonstrates that the inter-model spread in ${\stackrel{\mathrm{‾}}{w}}^{\ast }$ in the REF-C1SD simulations is comparable to that in REF-C1 at extratropical latitudes. In contrast, in the tropics between 30^∘ S and 30^∘ N the REF-C1SD simulations exhibit a slightly larger inter-model spread than the free-running simulations (0.09 mm s⁻¹ vs. 0.07 mm s⁻¹).

Figure 3d shows the absolute differences in ${\stackrel{\mathrm{‾}}{w}}^{\ast }$ between the REF-C1SD simulations and the respective reanalysis dataset used for nudging. In the upwelling region, the REF-C1SD experiments generally show stronger upwelling near the Equator than in the reanalyses. Although CESM1-WACCM is nudged towards MERRA, it does not simulate downwelling at the Equator as seen in the MERRA direct estimate. The relatively larger ${\stackrel{\mathrm{‾}}{w}}^{\ast }$ differences near 10–15^∘ N in CCSRNIES-MIROC3.2, EMAC-L90, and MRI-ESMr1 reflect a lack of inter-hemispheric asymmetry in the double-peaked ${\stackrel{\mathrm{‾}}{w}}^{\ast }$ structure in the REF-C1SD experiment compared to the reanalyses. Outside of the tropics, the REF-C1SD experiments generally show weaker downwelling in the NH mid-latitudes, while at polar latitudes (>65^∘) the REF-C1SD runs consistently show stronger downwelling than in the reanalyses. The difference in ${\stackrel{\mathrm{‾}}{w}}^{\ast }$ at high latitudes between the REF-C1SD and reanalysis datasets extends throughout the depth of the stratosphere (see Fig. S10). More generally, Fig. 3b shows that the different models that all nudge towards ERA-I (CCSRNIES-MIROC3.2, CMAM, EMAC-L47 and EMAC-L90, and SOCOL3) produce very different mean residual circulations.

In summary, we conclude based on the results in Figs. 1 to 3 that nudging meteorology affects the strength and structure of the climatological residual circulation throughout the stratosphere. However, as implemented in these simulations (Table 2), nudging neither strongly constrains the mean amplitude and structure of the residual circulation nor produces circulations that closely resemble the direct estimates from the reanalyses.

3.2 Climatological residual circulation: tropical upward mass flux

Figure 4 shows vertical profiles of the climatological TUMF between 100 and 3 hPa calculated from annual means of ${\stackrel{\mathrm{‾}}{w}}^{\ast }$ for the (a) REF-C1 and (b) REF-C1SD experiments and (c) their difference. Note the logarithmic x-axis scale and that the CCMI and S-RIP fields have been interpolated from their native model levels to a set of predefined common pressure levels, which are rather sparse in the upper stratosphere; hence the TUMF calculation could be different if it were performed on the native model grid of both CCMI models and the reanalyses.

Figure 4 Vertical profiles of climatological (1980–2009) tropical upward mass flux (10⁹ kg s⁻¹) averaged between the turnaround latitudes for (a) REF-C1 and (b) REF-C1SD, (c) differences (%) between REF-C1SD and REF-C1, and (d) % differences between REF-C1SD and the respective reanalysis used for nudging. Note the logarithmic x axis in panels (a) and (b).

In terms of the differences between the REF-C1SD and REF-C1 simulations (Fig. 4c), there is no consistent picture of the effect of nudging on the TUMF at different stratospheric levels. In the lowermost stratosphere between 70 and 100 hPa, most models (apart from EMAC-L90) simulate stronger TUMF in the REF-C1SD runs than in REF-C1. The largest TUMF differences in the lower stratosphere due to nudging occur in EMAC-L90 and SOCOL3, which show differences at 90 hPa of around −20 % and +25 %, respectively. In the middle stratosphere, between 10 and 70 hPa, some models show almost no difference in TUMF due to nudging (MRI-ESMr1), some show a stronger mass flux (CCSRNIES-MIROC3.2, CESM1-WACCM, and CMAM) and others show a weaker mass flux (EMAC-L47 and SOCOL3). In the upper stratosphere (above 10 hPa) the picture is also mixed, as half of the models show higher TUMF in the nudged experiments (CESM1-WACCM and EMAC-L47 and EMAC-L90) and the others show weaker TUMF (CCSRNIES-MIROC3.2, MRI-ESM1r1, and SOCOL3). CMAM shows the smallest change in TUMF in the upper stratosphere due to nudging. CESM1-WACCM is the only model to show a consistent sign of the TUMF differences between REF-C1SD and REF-C1 at all levels, with higher TUMF found throughout the stratosphere. There is no apparently simple relationship between the free-running model TUMF climatologies (Fig. 4a) and the effect of nudging (Fig. 4c).

We now compare the TUMF in each REF-C1SD experiment with the reanalysis it was nudged towards (Fig. 4d). Taking at first a broad view of the entire profiles, there is a resemblance between the profiles of TUMF differences in EMAC-L47 and SOCOL3 as compared to ERA-I, which may be related to the similarities in the implementation of nudging in these models; for example, vorticity and divergence were nudged with the same relaxation parameters (see Table 2). The CESM1-WACCM REF-C1SD simulation generally shows larger TUMF values than MERRA by up to 10 %–15 % apart from in the upper stratosphere, where they start to converge. MRI-ESM1r1 exhibits relatively better agreement of TUMF with JRA-55 throughout the stratosphere. Looking across the models, most of the REF-C1SD simulations simulate stronger upwelling than their respective reanalysis in the upper stratosphere, with differences reaching up to 30 %–35 % in the two EMAC models. In fact, EMAC-L47 and EMAC-L90 show a high degree of similarity in the vertical structure of the TUMF differences between REF-C1SD and ERA-I at pressures less than 30 hPa, despite showing substantial differences in the lower stratosphere. This could be because in EMAC nudging is only imposed strongly up to 10 hPa, while higher model layers have weakening nudging coefficients, as they serve as transition layers. In the middle stratosphere (50–20 hPa), most of the REF-C1SD models simulate a lower TUMF compared to the reanalysis. Again, a key message is that the nudged REF-C1SD simulations show a comparable, if not a slightly larger, spread in the climatological TUMF compared to the free-running REF-C1 simulations throughout almost the whole depth of the stratosphere.

To understand the dynamical factors that contribute to the modelled climatological residual circulation and its spread, Fig. 5 shows the annual mean TUMF at 70 hPa along with the downward control calculations (Sect. 2.2.2) to quantify the contribution of resolved and parameterized wave forcing to the TUMF. The black bars on the left show the TUMF diagnosed from ${\stackrel{\mathrm{‾}}{w}}^{\ast }$, and the grey bars on the right show the estimated contribution to TUMF from the EPFD (dark grey) and the orographic (medium grey) and non-orographic (light grey) gravity wave drag. Note that SOCOL3 did not provide wave forcing fields (Table 3), so we cannot perform the downward control calculations for that model.

Figure 5 Tropical upward mass flux at 70 hPa (left bars) along with downward control calculations (right bars) showing contributions from EPFD (dark grey), OGWD (medium grey), and NOGWD (light grey) for (a) REF-C1 and (b) REF-C1SD and the reanalyses. For CMAM, the NOGWD contributes negatively to TUMF and is indicated with two red horizontal lines inside the lighter grey bar.

In the free-running REF-C1 simulations (Fig. 5a), the estimated TUMF from the total wave forcing for the majority of the models (apart from CESM1-WACCM and EMAC-L90) slightly exceeds the TUMF calculated directly from ${\stackrel{\mathrm{‾}}{w}}^{\ast }$. Since these simulations are internally consistent, the imperfect match indicates that the downward control principle as applied here relies on the close but inexact applicability of certain assumptions, such as the system being in a steady state in response to a steady mechanical forcing (Haynes et al., 1991). The REF-C1 inter-model range in TUMF at 70 hPa is 5.74×10⁹ to 6.62×10⁹ kg s⁻¹ (inter-model standard deviation of 0.29×10⁹ kg s⁻¹). Comparing the CCMI results in Fig. 5a with the results from CCMVal-2 models (see Fig. 4.10; SPARC, 2010), the MMM TUMF at 70 hPa for the seven REF-C1 model simulations analysed here (6.05×10⁹ kg s⁻¹) is within the inter-model range of the 14 CCMVal-2 models, which show a MMM TUMF around 4 % weaker (5.8×10⁹ kg s⁻¹; SPARC, 2010). In terms of the contribution of the resolved wave forcing to the TUMF in the free-running simulations, there appears to be a decreased inter-model range (3.26×10⁹ to 5.33×10⁹ kg s⁻¹) in the present study compared with the CCMVal-2 models, albeit that study included more models (1.5×10⁹ to 5.5×10⁹ kg s⁻¹; SPARC, 2010). Some CCMI models have increased their horizontal resolution by up to a factor of 2 (CMAM, MRI-ESM1r1, SOCOL3) and also their vertical resolution by up to 80 vertical levels (MRI-ESM1r1) compared with CCMVal-2 models (Dietmüller et al., 2018), which could improve their ability to simulate resolved wave forcing. There is a notable feature of CMAM which shows that the NOGWD contributes negatively to TUMF (indicated with two red horizontal lines on Figs. 5 and S11); this was also found for CMAM in CCMVal-2 (Fig. 4.10; SPARC, 2010).

The MMM TUMF at 70 hPa in the REF-C1SD simulations (Fig. 5b) is 6.32×10⁹ kg s⁻¹, or around 5 % higher than in REF-C1. The REF-C1SD model range is larger than in REF-C1, being 5.39×10⁹ to 7.08×10⁹ kg s⁻¹ (inter-model standard deviation of 0.51×10⁹ kg s⁻¹). A notable feature is that the contribution from the individual and total wave forcing contributions shows reduced inter-model spread in the REF-C1SD simulations (Fig. 5b; darker grey bars). For example, the inter-model standard deviation of the EPFD contribution to TUMF at 70 hPa is around 40 % smaller than in REF-C1 (0.44×10⁹ and 0.72×10⁹ kg s⁻¹, respectively). Nonetheless, the residuals (i.e. the difference between the directly calculated TUMF and the total downward control estimated contribution from the wave forcing) are substantially larger and more positive (except for EMAC-L90) in the REF-C1SD experiment than in REF-C1. This shows that nudging adds an additional non-physical tendency in the model equations which acts to decouple the wave forcing from the residual circulation; this means that the physical constraint that the divergence of the angular momentum flux due to the mean motion is balanced over some sufficient time average by that of all eddy motions does not apply in the nudged models (Haynes et al., 1991). The details of how this decoupling is manifested are likely to vary from one model to another, depending on multiple factors such as nudging timescales, nudging parameters, nudging height range, and model resolution. Comparison of the TUMF at 10 hPa for the REF-C1SD experiment (see Fig. S11b) also reveals substantial differences in some models between the direct and downward control TUMF estimates in the middle stratosphere. Variations in the residuals as a function of height may indicate differences in the effect of nudging on the connection between the climatological wave forcing and the shallow and deep branches of the circulation (Birner and Bönisch, 2011). However, the inter-model ranges in the directly calculated TUMF at 10 hPa are more comparable in the two experiments than those found at 70 hPa (1.45×10⁹ to 1.70×10⁹ and 1.51×10⁹ to 1.72×10⁹ kg s⁻¹ for REF-C1 and REF-C1SD, respectively; Fig. S11b).

Interestingly, for the single simulations that were nudged towards MERRA and JRA-55 (CESM1-WACCM and MRI-ESM1r1, respectively), the TUMF at 70 hPa in the REF-C1SD runs appears to be close to the estimates from the reanalyses they are nudged towards (compare black bars in Fig. 5b). This may simply be a coincidence given that there are substantial differences in the structure of ${\stackrel{\mathrm{‾}}{w}}^{\ast }$ between the REF-C1SD simulations for those models and the reanalyses (Fig. 3b and d), and this is not found for all five models that were nudged towards ERA-I. Indeed, given that there is substantial spread in TUMF amongst the five REF-C1SD models nudged to ERA-I, it is likely that the differences between the REF-C1SD and reanalysis datasets are related to how nudging was implemented in each model; a wide variety of relaxation timescales and vertical nudging ranges were used by the models (Table 2). Despite this, the lower TUMF calculated directly from ${\stackrel{\mathrm{‾}}{w}}^{\ast }$ in EMAC-L90 compared to EMAC-L47, seen in both the REF-C1 and REF-C1SD experiments, is consistent with the results of Revell et al. (2015b), who also find that an increase in the model vertical resolution for SOCOL3 results in a slowdown of the BDC.

In summary, the results from Figs. 4 and 5 further demonstrate that nudging imparts an external and non-physical tendency in the model equations, which in turn might cause violations of the normal constraints on the global circulation, such as conservation of momentum and energy. This is found to alter the residual circulation but in a manner that cannot be understood from a closure of the circulation through the integrated wave forcing, as would ordinarily apply in the downward control principle (Haynes et al., 1991).

3.3 Annual cycle

We now evaluate the representation of the annual cycle in the residual circulation. Figure 6 shows the MMM climatological annual cycle of ${\stackrel{\mathrm{‾}}{w}}^{\ast }$ at 70 hPa for the REF-C1 and REF-C1SD simulations and their difference. Both experiments show similar broad features in the annual cycle, with stronger tropical upwelling in boreal winter, a latitudinal asymmetry in the region of upwelling, with the TA latitude being further poleward in the summer hemisphere, and stronger downwelling over the winter pole. These features resemble the annual cycle found in other multi-model studies (e.g. Hardiman et al., 2014). Figure 6c shows that on average the nudged models simulate stronger upwelling in the subtropics, particularly in the NH in boreal winter, with a few exceptions, the most prominent one being the narrow band between the Equator and 10^∘ N, where the REF-C1 simulations exhibit stronger upwelling in austral winter. Consequently, the nudged models simulate substantially stronger downwelling in the mid-latitudes in winter. In the NH mid-latitudes in the summer months, nudged runs show weaker downwelling, which reverses for the SH mid-latitudes in the austral winter. At polar latitudes there is a distinct seasonality to the differences between the REF-C1SD and REF-C1 simulations, with the nudged models simulating stronger downwelling in boreal winter and weaker downwelling in the Arctic during the rest of the year, corresponding to an amplified annual cycle. Conversely in the Antarctic, the REF-C1SD simulations generally simulate weaker downwelling, particularly during austral summer and spring.

Figure 6 Climatological MMM annual cycle in ${\stackrel{\mathrm{‾}}{w}}^{\ast }$ (mm s⁻¹) at 70 hPa for (a) REF-C1, (b) REF-C1SD, and (c) the REF-C1SD minus REF-C1 absolute differences. The shading in (c) denotes regions where the differences are statistically significant above 95 % using a two-tailed Student's t test. The turnaround latitudes (${\stackrel{\mathrm{‾}}{w}}^{\ast }=\mathrm{0}$) are shown by the thick black lines in (a) and (b) and by the thick red lines for the REF-C1SD MMM in (c).

To compare the annual cycle in residual circulation in the individual models, Fig. 7a and b show the mean tropical (30^∘ S–30^∘ N) ${\stackrel{\mathrm{‾}}{w}}^{\ast }$ at 70 hPa for the REF-C1 and REF-C1SD simulations, respectively. Comparing the MMM annual cycle of the REF-C1 runs (Fig. 7a) with the MMM REF-C1SD (Fig. 7b) reveals that on average the nudged models show a slightly larger peak-to-peak annual cycle amplitude (0.16 mm s⁻¹ vs. 0.13 mm s⁻¹). In general, the amplitude of the annual cycle in tropical mean ${\stackrel{\mathrm{‾}}{w}}^{\ast }$ is slightly more constrained across the REF-C1SD simulations with the spread in peak-to-peak amplitude, as measured by the inter-model standard deviation, being around 25 % smaller than in REF-C1 (σ=0.015 mm s⁻¹ vs. 0.020 mm s⁻¹, respectively). In terms of seasonal mean behaviour, the nudging appears to constrain the tropical mean ${\stackrel{\mathrm{‾}}{w}}^{\ast }$ in boreal summer (June–July–August – JJA), which exhibits ∼20 % less spread than in the REF-C1 experiments, but it does not constrain the tropical mean ${\stackrel{\mathrm{‾}}{w}}^{\ast }$ in boreal winter (December–January–February – DJF), which shows a larger spread than the free-running models by a factor of 2. Furthermore, the differences in tropical mean ${\stackrel{\mathrm{‾}}{w}}^{\ast }$ between the REF-C1SD runs and the respective reanalysis they are nudged towards are generally larger in boreal winter than in boreal summer for most models. In terms of spatially resolved differences in ${\stackrel{\mathrm{‾}}{w}}^{\ast }$ between REF-C1SD and the reanalyses (Fig. S12), some consistent features include the REF-C1SD simulations showing stronger downwelling in the Arctic in boreal winter compared to the reanalyses and showing weaker upwelling in the northern subtropics in boreal summer and autumn. Overall, the REF-C1SD minus reanalysis differences for the individual models highlight a wide variety in both the magnitude and the spatial patterns of their absolute differences, with no consistent picture emerging even for the models nudged towards the same reanalysis dataset.

Figure 7 (a, b) Climatological annual cycle in ${\stackrel{\mathrm{‾}}{w}}^{\ast }$ (mm s⁻¹) at 70 hPa between 30^∘ S and 30^∘ N in (a) REF-C1 and (b) REF-C1SD. (c, d) Climatological annual cycle in turnaround latitudes at 70 hPa for each model in (c) REF-C1 and (d) REF-C1SD.

Figure 7c and d show the climatological annual cycle in the TA latitudes at 70 hPa for the REF-C1 and REF-C1SD runs, respectively. This further breaks down the MMM annual mean perspective shown in Fig. 2 by model and by season. In the SH, the spread in seasonal mean TA latitude across models, as measured by the intermodel standard deviation, is increased in the REF-C1SD experiment in all seasons by up to 30 % compared to REF-C1. Conversely in the NH, the spread in seasonal mean TA latitude is decreased for REF-C1SD in all seasons except boreal spring (MAM), where it is increased. There are also substantial differences between the TA latitudes in the REF-C1SD experiment and the reanalyses in all months, which shows that nudging does not produce consistent structures of regions of upwelling and downwelling to those in the reanalysis. To summarize the results of Fig. 7, there is substantial inter-model spread in the TA latitudes and in the amplitude of the annual cycle in ${\stackrel{\mathrm{‾}}{w}}^{\ast }$, highlighting significant inter-hemispheric differences in the upwelling region between both sets of simulations and between the nudged experiment and the reanalyses.

3.4 Interannual variability in the tropical upward mass flux

Figure 8 shows time series over 1980–2009 for the annual, DJF, and JJA mean TUMF at 70 hPa for the REF-C1 (left column) and REF-C1SD (right column) simulations. As expected, the TUMF is larger in DJF compared to the annual and JJA means in both the REF-C1 and REF-C1SD runs because the average tropical upwelling is stronger in boreal winter. The individual REF-C1SD simulations show remarkably similar temporal variability in contrast to REF-C1, where the modelled inter-annual variability is very diverse despite the models all being forced with observed SSTs. Hence, although nudging does not constrain the mean TUMF in the lower stratosphere, it constrains the inter-annual variability; this is even more apparent for the DJF and JJA seasonal means (Fig. 8d, f). Additionally, the REF-C1SD simulations show a relatively high agreement in their temporal variability to the reanalysis datasets they were nudged towards, albeit with differences in magnitude and trend at the beginning of the 21st century, where ERA-I and MERRA show a decrease in TUMF.

Figure 8 Top to bottom: time series of annual, DJF, and JJA means of tropical upward mass flux (×10⁹ kg s⁻¹) at 70 hPa for (a, c, e) REF-C1 and (b, d, f) REF-C1SD.

To investigate the cause of the high temporal coherence of the REF-C1SD TUMF time series, Fig. 9 presents the annual mean TUMF anomalies at 70 hPa along with the relative contributions from EPFD, OGWD, NOGWD, and the total parameterized wave forcing (from top to bottom panels) for REF-C1 (left column) and REF-C1SD (right column), respectively. Figure 9b shows again the remarkably similar temporal variability in TUMF across the REF-C1SD runs, which can be contrasted against the weak inter-annual coherence in the REF-C1 runs (Fig. 9a). Figure 9d and j show that both the EPFD and the total parameterized wave forcing contributions to the TUMF show a high degree of temporal coherence in the REF-C1SD simulations. The fact that the individual OGWD and NOGWD terms do not show such a strong inter-model agreement, while the total parameterized wave forcing does, could suggest that there is some compensation occurring between the different parameterized wave forcing components (e.g. Cohen et al., 2013). It should be noted that the reanalyses have been shown to exhibit strong similarities in their resolved EP fluxes as shown by the linear correlation in the time series of tropical upwelling at the 70 hPa level when considering the momentum balance estimates of ${\stackrel{\mathrm{‾}}{w}}^{\ast }$ (Abalos et al., 2015). This result indicates that although nudging does not constrain the mean residual circulation, it constrains the inter-annual variability and produces similar contributions to variability across models from both resolved and parameterized wave forcing. In contrast, the REF-C1 simulations show a highly variable pattern of the estimated TUMF anomalies from EPFD and parameterized wave forcing (Fig. 9c and i), despite the fact that they use the same observed SSTs and that some nudge the phase of the QBO (CCSRNIES-MIRCO3.2, CESM1-WACCM, EMAC-L47 and EMAC-L90, and SOCOL3). In summary, the remarkably coherent inter-annual variability in the annual TUMF time series in the REF-C1SD simulations is due to both the resolved and parameterized wave forcing being constrained by nudging; this is in strong contrast to the climatological strength of the TUMF, where there were large differences between the directly calculated TUMF and that due to wave forces (Fig. 5b). The reasons for the difference in the effect of nudging on the behaviour of the residual circulation between the long-term mean and inter-annual variability are unclear.

Figure 9 Time series of the annual tropical upward mass flux anomalies (×10⁹ kg s⁻¹) calculated from (top to bottom) ${\stackrel{\mathrm{‾}}{w}}^{\ast }$ (a, b), and the downward control principle inferred contributions from resolved (EPFD) wave driving (c, d), orographic gravity wave drag (OGWD) (e, f), non-orographic gravity wave drag (NOGWD) (g, h), and from the total parameterized (OGWD and NOGWD) gravity wave drag (i, j) for REF-C1 (left panels) and REF-C1SD (right panels).

3.5 Multiple linear regression analysis

Figures 10 and 11 show time series of annual TUMF anomalies at 70 hPa attributed to each of the basis functions in the MLR model described in Sect. 2.3 and the regression residuals for the REF-C1 and REF-C1SD runs, respectively. Also shown in Figs. S13 and S14 are the regression coefficients for each term and for each model along with their uncertainties. Figure 10a shows a large spread in the diagnosed signal of volcanic eruptions in the TUMF time series. The majority of the REF-C1 simulations analysed here show a negative TUMF anomaly around the time of the El Chichón (1982) and Mount Pinatubo (1991) eruptions; however, the magnitude is within the estimated uncertainty range for all models except SOCOL3 (Fig. S13). In contrast to the REF-C1 results, most REF-C1SD simulations (except EMAC-L47 and EMAC-L90 – see above discussion) show a positive anomaly in TUMF attributed to volcanic eruptions (Fig. 11), consistent with earlier studies (Garcia et al., 2010; Diallo et al., 2017). However, there is still a considerable range of amplitudes, and only the CESM1-WACCM and MRI-ESM1r1 regression coefficients are highly significant (Fig. S14). The issue of establishing a robust response of the TUMF to volcanic forcing over a short period is demonstrated by the range in amplitudes of the volcanic regressors for different REF-C1 ensemble members from the same model (see Figs. S6–S9). This highlights that in a free-running climate simulation, internal variability can overwhelm the response to forcing over short timescales. The “true” volcanic signal in TUMF will also depend on the representation of stratospheric heating due to aerosol in the various models. We note that the EMAC-L47 and EMAC-L90 models contained a unit conversion error where the extinction of stratospheric aerosols was set too low by a factor of ∼500 (see Appendix B4 of Morgenstern et al., 2017); hence the stratospheric dynamical effects of the eruptions were not properly represented in the EMAC simulations (Jöckel et al., 2016).

Figure 10 Time series for REF-C1 simulations of the components of the annual mean tropical upward mass flux (×10⁹ kg s⁻¹) attributed to (a) volcanic aerosol, (b) ENSO, (c) linear trend, (d, e) QBO, and (f) regression residuals.

Figure 11 Time series for REF-C1SD simulations of the components of the annual mean tropical upward mass flux (×10⁹ kg s⁻¹) attributed to (a) volcanic aerosol, (b) ENSO, (c) linear trend, (d, e) QBO, and (f) regression residuals.

The REF-C1 models all show a positive best-estimate regression coefficient for the TUMF response to ENSO (Fig. 10), which is quite consistent in amplitude, but it is only strongly statistically significant in CCSRNIES-MIROC3.2 and SOCOL3 (Fig. S5). This is in contrast to the REF-C1SD models, which all show a larger and more significant positive ENSO regression coefficient. The linear trend regression coefficient over 1980–2009 is positive in all REF-C1 models and is statistically significantly different from zero at the 95 % confidence level in five out of the seven models. The magnitude of the linear trend term varies by around a factor of 2 for REF-C1. In REF-C1SD, the amplitude of the linear trend regression coefficient increases in all models, but the intermodel spread increases to around a factor of 4. Hence, in these simulations nudging increases the disparity across models in the magnitude of the long-term TUMF trend.

Figure 12 Tropical upward mass flux trends at 70 hPa (×10⁹ kg s⁻¹ decade⁻¹) for different start (abscissa) and end (ordinate) dates over the period 1980–2009 for the REF-C1 (r1i1p1) simulations. Trends are not shown for periods of less than 10 years. Values with statistical significance greater than the 95 % level are shaded.

Figure 13 As in Fig. 12 but for the REF-C1SD simulations.

As expected, the variations in TUMF attributed to the QBO is quite different in the REF-C1 and REF-C1SD runs for those models that do not nudge the QBO in REF-C1, as shown in Figs. 10 and 11. The nudging of zonal winds in REF-C1SD constrains the phase of the QBO, and hence there is strikingly similar variability in the TUMF anomalies attributed to the QBO in the REF-C1SD runs.

The overall R² values from the MLR model for the REF-C1 simulations vary between 0.16 (CMAM) and 0.67 (CESM1-WACCM). REF-C1SD runs generally give more consistent R² values across the models, ranging from 0.62 (CCSRNIES-MIROC3.2) to 0.77 (EMAC-L47). This means that there is still a substantial fraction (>23 %) of unexplained variance in the annual TUMF time series in the REF-C1SD simulations after applying the MLR model, and the residuals exhibit a remarkable degree of temporal correlation. In contrast, the MLR residuals in the REF-C1 runs (Fig. 10f) show much less temporal coherence apart from a drop around 1989. The residuals in the REF-C1SD simulations (Fig. 11f) show a high degree of coherent inter-annual variability, another manifestation of the fact that the nudged runs do reproduce a much more consistent inter-annual variability. This makes a substantial contribution to the coherence of the TUMF time series in Fig. 9b, but it cannot be attributed to any of the terms included in the MLR model.

For completeness, the MLR model was also applied to the reanalysis TUMF at 70 hPa (Figs. S15 and S16). This highlights significant discrepancies in attributing the variance in TUMF in the different reanalysis datasets to the various basis functions in the MLR model. Both the volcanic activity and ENSO contributions to the variance in the TUMF are rather weak compared to the REF-C1SD runs. The negative linear trend in ERA-I is in strong contrast to the positive trends found in the other reanalyses and the REF-C1SD models. The negative trend in ERA-I found in the TUMF in the lower stratosphere over 1980–2009 corroborates the findings of Abalos et al. (2015), who showed a negative trend in the direct ${\stackrel{\mathrm{‾}}{w}}^{\ast }$ estimate in ERA-I over 1979–2012. Despite this difference in the representation of the long-term TUMF trend, ERA-I shows the highest percentage of TUMF variance explained by the MLR model (66 %), with MERRA showing a substantially lower R² (0.3) compared to the other reanalyses and the REF-C1SD models. The residuals are generally less correlated between the reanalyses on inter-annual timescales than those found in the REF-C1SD simulations but are broadly similar on inter-decadal timescales. However, the regression residuals in the reanalyses show a different temporal behaviour from those in the REF-C1SD simulations (Fig. 11; note that the y-axis scale for the residuals in Fig. S15 is double that for the CCMI models in Figs. 10 and 11). In summary, although nudging constrains the inter-annual variability in the TUMF at 70 hPa, the attribution to some specific drivers differs across the models and, in comparison, to the reanalyses they were nudged towards.

3.6 Trend sensitivity analysis

Following the results of the MLR analysis described in Sect. 3.5, which showed a statistically significant positive linear trend in most REF-C1 and REF-C1SD models for the 30-year period 1980–2009, we now explore the sensitivity of the linear trend to the time period considered. We apply the same MLR model as discussed in Sect. 3.5 to the annual mean 70 hPa TUMF time series of the first ensemble member for both REF-C1 and REF-C1SD runs as well as the reanalyses but systematically vary the start and end dates to cover all time periods in the window 1980–2009 that are at least 10 years in length. We then extract the linear trend coefficient from the MLR model and its associated p value. Figures 12 and 13 present the linear trend calculations for the REF-C1 and REF-C1SD runs, respectively, as a function of trend start and end date. The same trend sensitivity analysis for the reanalyses is presented in the Supplement (Fig. S17). Statistically significant trends at the 95 % confidence level are marked with black shading.

None of the periods considered in either the REF-C1 or REF-C1SD experiments show a significant negative TUMF trend. A statistically significant positive trend emerges in almost all of the REF-C1SD models for trends beginning in the mid-1980s to early 1990s extending to the mid-2000s. The trends are mainly significant for periods of 20 years or more and no less than around 12 years. This result broadly corroborates the findings of Hardiman et al. (2017a), who used a control run to estimate the period required to detect a BDC trend with an amplitude of 2 % per decade against the background internal variability. There is range of different structures in the diagnosed trends among models, particularly for the REF-C1 simulations, where a consistent pattern of positive trends only emerges across most models for the entire time period. This is because internal variability can mask BDC trends over short periods (Hardiman et al., 2017a). However, the REF-C1SD runs simulate more consistent variations in TUMF trends as a function of time period but generally show stronger positive trends than their free-running counterparts. Interestingly, the reanalysis trend sensitivity analysis highlights that nudging does not constrain the underlying trends of the REF-C1SD models in the TUMF at 70 hPa, as the reanalysis datasets exhibit a wide range of different trends from one another (Fig. S17) and differences compared to the trends in the REF-C1SD simulations (Fig. 13). For example, none of the REF-C1SD models simulate a statistically non-significant negative trend in TUMF, starting around mid-1990s up to 2009, as seen in all the reanalyses. However, it should also be noted that any trend combination starting around the end of 1990s in almost all cases of both REF-C1 and REF-C1SD runs exhibit no statistical significance possibly pointing towards the role of declining ODSs due to the implementation of the Montreal Protocol (Polvani et al., 2018).

4 Conclusions

This study has performed the first multi-model intercomparison of the impact of nudged meteorology on the representation of the stratospheric residual circulation. We use hindcast simulations over 1980–2009 from CCMI, with identical prescribed external forcings in two configurations: REF-C1SD with meteorological fields nudged towards reanalysis data (specified dynamics – SD) and REF-C1 that is free-running. The nudged simulations use one of three different reanalysis datasets (ERA-Interim, JRA-55, and MERRA), nudge different variables (u, v, T, vorticity, divergence, and surface pressure), and use different time constants to impose the additional nudging tendencies in the model equations.

The key findings of this study are as follows:

1. Nudging meteorology does not constrain the mean strength of the residual circulation compared to free-running simulations. In fact, for most of the metrics of the climatological residual circulation examined, including residual vertical velocities and mass fluxes, the inter-model spread is comparable or in some cases larger in the REF-C1SD simulations than in REF-C1.

2. Nudging leads to the models simulating on average stronger upwelling at the Equator in the lower stratosphere to middle stratosphere and a wider tropical pipe in the lower stratosphere. In most cases, the magnitude and structure of the climatological residual circulation in the REF-C1SD experiments differ markedly from those estimated for the reanalysis they are nudged towards.

3. In most of the nudged models there are large differences of up to 25 % between the directly calculated tropical upward mass flux in the lower stratosphere and that calculated from the diagnosed total wave forcing using the downward control principle (Haynes et al., 1991). However, the spread in the contributions from the resolved and parameterized wave forcing to the tropical mass flux is slightly reduced in the REF-C1SD simulations compared to REF-C1.

4. Despite the lack of agreement in the mean circulation, nudging tightly constrains the inter-annual variability in the tropical upward mass flux (TUMF) in the lower stratosphere. This is associated with constraints to the contributions from both the resolved and parameterized wave forcing despite the fact that the models use different reanalysis datasets for nudging. The reanalysis datasets themselves exhibit broadly similar inter-annual variability in TUMF in the lower stratosphere, albeit with different long-term trends.

5. A multiple linear regression (MLR) analysis shows that up to 77 % (67 %) of the inter-annual variance of the lower-stratospheric TUMF in the REF-C1SD (REF-C1) experiments can be explained by volcanic eruptions, ENSO, the QBO, and a linear trend. The remaining unexplained TUMF variance in the nudged models shows a high degree of a temporal coherence, but this is not the case for the free-running simulations.

6. The results of the MLR analysis applied to the TUMF in the reanalyses show differences in the total variance explained and the attribution of variance to the different physical proxies. There are also marked differences between the individual regression coefficients derived for the REF-C1SD models and the reanalysis dataset used for nudging.

7. Most nudged simulations show a statistically significant positive trend in TUMF in the lower stratosphere over 1980–2009, which is on average larger than the trends simulated in the free-running models. This is despite the fact that five out of the seven models analysed were nudged towards ERA-Interim, which shows a negative long-term trend in TUMF (see also Abalos et al., 2015), while JRA-55 and MERRA show a positive trend. However, the magnitude of the TUMF trend varies by up to a factor of 4 across the nudged models, which is larger than the spread in the free-running simulations. This is an important limitation for using nudged-CCM simulations to interpret long-term changes in stratospheric tracers.

8. A sensitivity analysis of the time period for calculating lower-stratospheric TUMF trends shows that a statistically significant (at the 95 % confidence level) positive trend in TUMF takes at least 12 years and in most cases around 20 years to emerge in the REF-C1SD runs. Despite the three reanalysis datasets showing different 30-year trends (1980–2009), they show a striking agreement in the statistically non-significant negative trends starting from the late 1990s up to 2009.

Our findings highlight that nudging strongly affects the representation of the stratospheric residual circulation in chemistry–climate model simulations, but it does not necessarily lead to improvements in the circulation. Similar disagreement in the characteristics of tropospheric transport in the CCMI nudged simulations has also been reported (Orbe et al., 2018). The differences found in the nudged runs compared with the free-running simulations suggest that although nudging horizontal fields can remove model biases in, for example, temperature and horizontal wind fields (Hardiman et al., 2017b), the simulated vertical wind field will not necessarily be similar to the reanalysis. A particularly interesting finding of our study is that while nudging does not constrain the mean strength of the residual circulation, it constrains the inter-annual variability. The reason for the distinct effects of nudging on the residual circulation across these different timescales is currently unknown.

Multiple factors are likely to determine the effect of nudging on the residual circulation in a given model, including model biases, nudging timescales, nudging parameters, nudging height range, and model resolution. The differences in the stratospheric residual circulation between the REF-C1SD and the REF-C1 runs may not arise solely from the dynamics but can also be partly influenced by the indirect effects of nudging the temperatures, which in turn affect the diabatic heating (Ming et al., 2016a, b). In addition to nudging the horizontal winds (mechanical nudging), nudging the temperature (thermal nudging) might be systematically creating a spurious heat source in the model, which leads to a stronger BDC in the lower stratosphere, as suggested by Miyazaki et al. (2005) with the MRI GCM. Our results highlight that in the method by which the large-scale flow is specified and more specifically the choice of the reanalysis fields, the relaxation timescale and the vertical grid (pressure level versus model level) in which the nudging is applied need to be better understood and evaluated for their influence on the stratospheric circulation. Discrepancies between the vertical grid of the models and the reanalysis pressure levels they are interpolated onto or unbalanced dynamics are possible explanations for the differences found between the directly inferred circulation and that diagnosed from the wave forcing in the nudged simulations. Nudging would either violate continuity, or if continuity is maintained, it will come at the expense of the vertical fluxes, which are not nudged. The interesting aspect here seems to be that this results in substantial change to the net fluxes across a range of timescales; i.e. it does not only increase numerical noise in the ${\stackrel{\mathrm{‾}}{w}}^{\ast }$ component. In order to reduce discrepancies between nudged and free-running simulations, various nudging techniques have been investigated. The role of gravity waves in the error growth that the nudging introduces over time has been highlighted for a single model (Smith et al., 2017). Constraining just the horizontal winds without the temperature was found to be a good strategy when investigating the aerosol indirect effects without affecting significantly the mean state (Zhang et al., 2014). The relaxation timescale when applying the nudging has been found to play an important role in single-model studies (Merryfield et al., 2013), but there is no general consensus for the value of the relaxation constant, which is specific to the model for the simulations considered here (Morgenstern et al., 2017). Given the varying implementations of nudging in the models analysed here, our study is ill-suited to investigate in detail the mechanisms for how nudging affects the residual circulation. A dedicated study of the sensitivities within one model to relaxation timescales, nudging parameters, nudging height range, and vertical resolution would help with offering a detailed explanation for these differences.

The large spread in climatological residual circulation in nudged-CCM simulations is an important limitation for those wishing to use them to examine tracer transport, for example in stratospheric ozone trends (Solomon et al., 2016), volcanic aerosols (Schmidt et al., 2018), and diagnostics for the age of air (Dietmüller et al., 2018). Despite the limitations for transport within the stratosphere described here, some success has been reported in studies that used nudged simulations to investigate specific meteorological events such as sudden stratospheric warmings and, in particular, for exploring processes beyond the top of the nudging region in the mesosphere–lower thermosphere (e.g. Tweedy et al., 2013; Chandran and Collins, 2014; Pedatella et al., 2014). In conclusion, owing to the limitations of the current techniques for nudging models highlighted here, we urge caution in drawing quantitative comparisons of stratospheric tracers affected by the residual circulation in nudged simulations against stratospheric observational data.