## Abstract

The wintertime northern annular mode (NAM) at the surface is known to undergo slow intraseasonal variations in association with stratospheric variability, which leads the surface signal by up to several weeks. The relative contributions, however, of potentially relevant stratosphere–troposphere coupling mechanisms are not yet fully understood.

In this study the relative roles of (i) the downward effect of the zonal-mean secondary circulation induced by quasigeostrophic (QG) adjustment to stratospheric wave drag and radiative damping and (ii) wave drag local to the troposphere are estimated. For this purpose, a spectral tendency equation of the QG zonal-mean zonal wind is derived and used, in a first step, to obtain the external mechanical forcing that, in the QG framework, drives exactly the observed stratospheric and tropospheric daily NAM. In a second step, the equation is then integrated in time to reconstruct the daily NAM, but with the forcing restricted to either stratospheric or tropospheric levels, each case leaving a characteristic NAM surface signal.

The relative roles of the above-mentioned mechanisms are found to be of similar quantitative importance, but to differ in a qualitative sense. The downward effect of stratospheric QG adjustment is responsible for the initiation of the NAM surface signal, whereas subsequently local tropospheric wave drag actively maintains and persists the signal over several weeks. Furthermore, the downward effect of QG adjustment to stratospheric radiative damping is shown to have only a minor impact, compared to that from stratospheric wave drag. The robustness of these conclusions is demonstrated by a sensitivity study with respect to various model parameters.

## 1. Introduction

The wintertime northern annular mode (NAM) (Thompson and Wallace 2000) represents the leading low-frequency variability mode of the zonally averaged circulation in both the troposphere and stratosphere and characterizes deep, hemispheric-scale fluctuations of mid- to high-latitude westerlies. It is closely related to the Arctic Oscillation and North Atlantic Oscillation phenomena and, thus, has a similarly profound impact on Northern Hemisphere surface climate variability (see, e.g. Hurrell 1995; Thompson and Wallace 1998; Wallace 2000).

Intraseasonal predictability, however, of the tropospheric NAM is strongly limited since it is largely driven by relatively fast internal tropospheric processes and exhibits an autocorrelation *e*-folding time scale of only about 10 days, whereas at stratospheric levels much longer time scales of up to 30 days are observed (Mudryk and Kushner 2011). Additionally, a number of studies have shown that stratospheric variability has a noticable impact on the tropospheric circulation, and this coupling is manifested in terms of concurrent NAM anomalies in the stratosphere and troposphere (see, e.g. Baldwin and Dunkerton 1999, 2001). Although the corresponding signal in the troposphere is small, it occurs on time scales of 1–2 months in association with slow variations in the stratosphere, and, hence, the downward effect of stratospheric variability enhances the predictability of the tropospheric NAM on intraseasonal time scales, as demonstrated by, for example, Charlton et al. (2003) and Baldwin et al. (2003).

Different dynamical mechanisms have been discussed and proposed in the literature for the observed stratosphere–troposphere connection, including planetary wave reflection in the stratosphere (Perlwitz and Harnik 2003), internal tropospheric eddy feedbacks (e.g., Wittman et al. 2004; Charlton et al. 2005), and the downward effect of quasigeostrophic (QG) adjustment to stratospheric anomalies (e.g., Ambaum and Hoskins 2002; Black 2002; Charlton et al. 2005; Thompson et al. 2006). However, the relative contributions of the various mechanisms are not yet fully understood and conclusions are still controversial.

The aim of this study is to clarify the relative roles of (i) the downward effect of QG adjustment to stratospheric wave drag and radiative damping and (ii) wave drag local to the troposphere, which may reflect internal tropospheric eddy feedbacks. For this purpose, the stratospheric and tropospheric NAM is computed from a reanalysis dataset, which serves as a proxy for observations, and the contributions of stratospheric and tropospheric forcing to NAM surface variability are then obtained from a simple linear QG model of the zonally averaged circulation.

The outline of this study is as follows. Section 2 presents the data source, the derivation of the model equation, and introduces the method to estimate the relative roles of the above-mentioned stratosphere–troposphere coupling mechanisms. The results obtained from the application of this method are presented in section 3 and further discussed in section 4. A summary of the main conclusions is given in section 5.

## 2. Methodology

This section describes, first, the data and the algorithm used to compute the wintertime NAM. Subsequently, the spectral QG model equation is derived, following the concepts of tidal theory, and, finally, the method to separate the stratospheric from the tropospheric contribution to NAM variability is presented.

### a. Data and computation of the northern annular mode

The wintertime stratospheric and tropospheric northern annular mode is derived from the 40-yr European Centre for Medium-Range Weather Forecasts Re-Analysis (ERA-40) (Uppala et al. 2005), which spans the 45-yr period from 1 September 1957 to 31 August 2002 and is available on 23 pressure levels ranging from 1000 hPa near the surface up to 1 hPa near the stratopause. Daily fields of zonal-mean geopotential at 0000 UTC are used for this study. The climatological seasonal cycle is removed from the time series at each grid point by setting those Fourier modes to zero that correspond to the mean and to frequencies of 1 ×, 2 ×, 3 ×, and 4 × 45 cycles per 45 yr, with 29 February included during the 11 leap years contained in the data series.

The daily NAM index is then computed, separately at each level, by projection of daily anomaly fields north of 20°N onto the first EOF obtained from December to March monthly mean anomaly fields north of 20°N. Subsequently, only the extended winter season from November to April (of length 181 days, omitting 30 Apr during leap years) is retained from the index time series at each level, which is then normalized to have zero mean and unit variance. For both, the projection as well as the EOF computation, fields are weighted with the square root of the cosine of latitude prior to the operation. Finally, the associated global NAM pattern at each level is constructed by linear regression of November–April daily anomalies onto the corresponding index time series. The obtained patterns reveal the well-known meridional structure of atmospheric annular modes, characterized by geopotential anomalies of one sign over the polar cap and of opposite sign at lower latitudes (not shown; see, e.g., Baldwin and Dunkerton 2001). The variance explained by this mode amounts to about 50% in the troposphere and about 80% at stratospheric levels, with a maximum of 84% at 30 hPa, a minimum of 43% at 300 hPa, and 51% at the surface at 1000 hPa, if related to daily zonal-mean geopotential variance north of 20°N. In terms of zonal-mean zonal wind, when regressed onto the same index time series, it explains about 20% in the troposphere and between 50% and 70% at most stratospheric levels.

The sign of the NAM index is chosen such that positive values correspond to stronger than normal mid- to high-latitude westerlies, implying, in the troposphere, a strengthened and poleward displaced eddy-driven jet and, at stratospheric levels, an intensified polar night jet on the edge of an anomalously cold polar vortex and vice versa for negative values. A sample time series of the daily NAM index from winter 1998/99 is shown later (see Fig. 2b). The same winter has been highlighted by Baldwin and Dunkerton (2001) and was characterized by two major sudden stratospheric warmings (see Charlton and Polvani 2007), associated with an anomalously warm polar vortex and a negative NAM index (as indicated by red colors in the figure).

### b. Model equations

The approach to estimate the relative contributions of stratospheric versus tropospheric forcing to NAM surface variability is based on assuming QG transformed Eulerian mean (TEM) flow, driven by zonally uniform thermal and mechanical forcing, in spherical and log–pressure coordinates, given by (e.g., Garcia 1987)

where zonal averaging is denoted by overbars, is the zonal wind, and are the northward and upward TEM residual velocities, is temperature, and is geopotential. The nonconservative terms on the rhs of Eqs. (1) and (4) represent the mechanical and thermal forcing, respectively, where and are external forcings, *κ*(*z*) is a Rayleigh drag coefficient, and *α*(*z*) a Newtonian cooling coefficient. Hereafter, the above dependent variables shall represent anomalies from the climatological mean seasonal cycle. The term *f* = 2Ω sin*φ* is the Coriolis parameter, Ω is the angular frequency, and *a* is the radius of the earth, and *N*^{2}(*z*) = *R*/*H*(*T*_{0z} + *κT*_{0}/*H*) is the square of the buoyancy frequency, with a reference temperature profile *T*_{0}(*z*), and a reference density profile is given by *ρ*_{0} ∝ exp(−*z*/*H*). The meridional and vertical coordinates are latitude *φ* and log–pressure height *z* = −*H* ln(*p*/*p _{s}*), respectively, with the scale height

*H*= 7 km and a reference pressure

*p*= 1000 hPa. Also,

_{s}*R*is the gas constant for dry air. Partial derivatives with respect to

*φ*,

*z*, and time

*t*are indicated by subscripts. For zero mechanical and thermal forcing, Eqs. (1)–(5) represent the conservation laws for eastward, northward, and upward momentum, for thermodynamic energy, and for mass, respectively.

These equations can be combined into an elliptic partial differential equation for the zonal-mean zonal wind tendency [e.g., Haynes and Shepherd 1989, their Eq. (3.5)],^{1}

with the linear differential operators

and where the modified variables

have been used, with and . Then, following the concepts of tidal theory (Chapman and Lindzen 1970; Andrews et al. 1987), by separation of the horizontal from the vertical dependence according to

one may from Eq. (6), with *F* = *Q*_{*} = 0 and with uniform damping coefficients *κ* = const and *α* = const, obtain two eigenvalue equations, namely, the horizontal structure equation (HSE; a special case of Laplace's tidal equation)

and the vertical structure equation (VSE)

with the eigenvalues *ε* = *λ*(*α* − *σ*)(*κ* − *σ*)^{−1}, where *λ* = 4Ω^{2}*a*^{2}(*gh*)^{−1} is the Lamb parameter, and *γ* = −(*gh*)^{−1}, which are related through the separation constant *h*, known from tidal theory as the equivalent depth, and *g* is the gravitational acceleration.

The eigenmodes of the HSE, when subject to proper horizontal boundary conditions, form an infinite set of horizontal structure functions *U*^{(n)} [the so-called Hough functions; see Longuet-Higgins (1968), for reference] with associated eigenvalues *ε*^{(n)}. For the VSE, however, the form taken by the lower boundary condition depends on the type of forcing applied. Under purely mechanical forcing the eigenmodes of the VSE form a single infinite set of vertical structure functions with associated eigenvalues , whereas under purely thermal forcing they form, separately for each *n*, one infinite set of vertical structure functions with associated eigenvalues , since in the latter case the lower boundary condition also depends on the horizontal eigenvalue *ε*^{(n)}. The vertical modal indices of the *Z _{F}* and the

*Z*modes are denoted by

_{Q}*m*and , respectively. The choice of the upper boundary condition and, in particular, the lower boundary condition is not trivial, however, as it can have a significant impact on the surface response to a given forcing above [as discussed in detail by Haynes and Shepherd (1989)]. Moreover, in the context of the present study, additional assumptions have to be made so as to obtain boundary conditions (and, thus, vertical structure functions) that are independent of surface eddy heat fluxes, which, in turn, have implications for the separation between stratospheric and tropospheric wave drag (for the separation method see section 3c below). Further details regarding the boundary conditions and related issues can be found in the appendix.

The horizontal and vertical structure functions, when normalized accordingly, fulfill the normalization and orthogonality conditions

with *δ _{ij}* = 0, if

*i*≠

*j*, and

*δ*= 1, if

_{ij}*i*=

*j*, and the top of the model domain

*z*. Orthogonality follows from using the modified variables (11) and can be demonstrated by two integrations by parts of the respective eigenvalue equation for one mode, multiplied by another mode, and application of the corresponding boundary conditions. Subsets of the various structure functions are illustrated in Fig. 1, and eigenvalues are presented in Tables 1 and 2. Note, that the first vertical mode under mechanical damping, , is external and represents the Lamb mode, whereas under thermal damping all modes are internal, as is also the case for the horizontal modes. These horizontal and vertical structure functions can then be combined into two-dimensional spectral modes

_{t}*UZ*(

*φ*,

*z*), each of which, when subject to uniform damping, has its own fixed decay time scale

*σ*

^{−1}, given by the corresponding pair of eigenvalues.

The spectral tendency equation for the zonal wind can now be derived from the elliptic equation (6) by projection onto the two-dimensional spectral modes *UZ*. For this purpose, quantities need to be expanded in terms of those modes as

in case of purely mechanical forcing (*Q*_{*} = *α* = 0), or as

in case of purely thermal forcing (*F* = *κ* = 0), such that the spectral coefficients are given by the projections

and

respectively [using conditions (15)], and where *X* stands for any of the terms *u _{t}*,

*F*,

*Q*

_{*},

*κu*, or

*αu*in Eq. (6). Then, by multiplication of Eq. (6) by 4Ω

^{2}

*a*

^{2}, expansion of the aforementioned terms into their spectral components, use of the HSE (13) and VSE (14), projection onto either the (

*m*,

*n*)th or onto the th spectral mode, and by use of conditions (15), for each type of forcing a spectral zonal wind tendency equation is obtained, that is, an equation for either or for , valid under arbitrary damping profiles

*κ*(

*z*) or

*α*(

*z*). These equations can be combined into a single spectral tendency equation,

with the operator

and where the terms

represent the effect of the secondary circulation (,) induced by QG adjustment of the zonal-mean atmosphere to either mechanical or thermal forcing, respectively; and *Q*_{*} has been set to zero. This latter assumption means, in the context of the present study, that we suppose the external heating term *Q*_{*} to be of only minor importance to intraseasonal NAM variability, and it also enables us to invert the spectral tendency equation (20) to obtain the external mechanical forcing *F* that drives, in the QG framework, a given time series of the zonal wind *u*,

### c. Method to estimate the contribution of stratospheric versus tropospheric forcing

The above spectral tendency equation (20) and its inverse counterpart equation (23) are now applied, as follows, to analyze NAM variability driven by either stratospheric or tropospheric external mechanical forcing *F*. First, the discretized versions of the HSE (13) and VSE (14) are solved numerically, using a horizontal discretization interval of Δ*φ* = *π*/102, an interval of Δ*z* = 1 km in the vertical, and with the model top at *z _{t}* = 48 km. This yields a set of 102 horizontal structure functions and sets of 48 vertical structure functions. Subsequently, the method proceeds in two steps.

In step 1, the inverse equation (23) is used to obtain the forcing *F* that drives precisely the zonal wind time series associated with the ERA-40 wintertime NAM. The latter time series is constructed from two components as *u*(*φ*, *z*, *t*) = *u*_{NAM}(*φ*, *z*) × NAMI(*z*, *t*). The first component is given by the geostrophic wind relation , with , and where is obtained by interpolation of the ERA-40 global NAM patterns onto the model grid. For simplicity, we set the Southern Hemisphere equal to minus the Northern Hemisphere to make *u*_{NAM} antisymmetric about the equator. Additionally, to avoid nonzero geostrophic zonal winds near the equator, *u*_{NAM} is set equal to zero between 20°S and 20°N. The resulting discontinuities at those bounding latitudes are removed by retaining only the first 15 odd horizontal modes (*n* = 1, 3, 5, … , 29), which leaves the extratropical pattern virtually unchanged. The second component is obtained by interpolation of the ERA-40 daily NAM index onto the vertical model grid and onto time steps with an interval of 0.1 day. The resulting zonal wind time series, in spectral representation, is then prescribed to the rhs of Eq. (23), separately for each November–April extended winter season.

In step 2, the spectral tendency equation (20) is integrated in time, separately for each winter season, but with the forcing time series *F*, previously obtained from step 1, restricted to either stratospheric or tropospheric levels, by setting *F* = 0 at *z* < *z _{c}* or

*z*>

*z*, respectively, with a cutoff height at

_{c}*z*= 16 km (corresponding to the 100-hPa pressure level). Initial conditions for each season are set to the 0000 UTC 1 November value of the zonal wind time series used in step 1. From the resulting geopotential time series , implied from

_{c}*u*(

*φ*,

*z*,

*t*) by the geostrophic relation, the NAM index is then computed for each case, but always normalized in relation to the ERA-40 NAM index time series, which, by construction, is identical to the NAM index obtained from the case forced at all levels.

Finally, the static stability and damping parameters need to be specified. The static stability profile, given by *N*^{2}(*z*), is chosen such as to correspond to the *U.S. Standard Atmosphere 1976* (COESA 1976).^{2} Mechanical damping is confined to near the surface to mimic surface drag within the planetary boundary layer. Specifically, the vertical profile of the damping coefficient decreases linearly from *κ* = *κ _{s}* at

*z*= 0 toward

*κ*= 0 at and above

*z*= 3 km. The effects of radiative damping and of small-scale turbulent heat fluxes are parameterized by the Newtonian cooling term −

*αT*and the damping coefficient is set to

*α*=

*α*at

_{s}*z*= 0, to

*α*= (40 days)

^{−1}at the tropopause, and to

*α*= (5 days)

^{−1}at the stratopause and varies between these levels according to a half-cosine mode. The values for stratospheric thermal damping time scales are based on, and are in close agreement with, estimates of radiative damping time scales in the stratosphere, derived from both observational data (Newman and Rosenfield 1997) as well as a chemistry climate model (Hitchcock et al. 2010). For the key scenario, described in section 3, we choose as default surface values

*κ*= (1 day)

_{s}^{−1}and

*α*= (40 days)

_{s}^{−1}.

Sample time series of the forcing *F*, obtained from step 1, and of the NAM index, obtained from step 2, under the above static stability and damping profiles, are illustrated in Fig. 2, for winter 1998/99. The equator-to-pole average of the forcing *F* is shown in Fig. 2a. The bulk of this forcing is due to planetary wave drag and, at upper stratospheric levels, partly due to gravity wave drag. However, since here we use the QG set of equations whereas the ERA-40 reanalysis is the output of a primitive equation model, *F* includes other terms like nonlinear advection of relative angular momentum and of temperature anomalies, diffusive processes due to small-scale mixing, and the effects of any other processes that are not accounted for by our QG approximation. However, in the context of the present study, the effects of these further terms are assumed to play only a minor role owing to the large spatial and long time scales of NAM variability compared to, for example, individual baroclinic systems where deviations between QG and primitive equation dynamics become much larger. Therefore, we refer to *F* simply as wave drag hereafter, although this interpretation is not exact. Evident from Fig. 2a are two episodes with strong easterly wave drag in the middle to upper stratosphere, which triggered the associated major sudden stratospheric warmings during that winter, as seen by the negative NAM anomalies in Fig. 2b; also, the shorter time scale of the tropospheric forcing is obvious. Note that, in contrast to the NAM index, the forcing is not normalized but is shown in units of meters per second per day and, thus, appears less emphasized in the figure at lower stratospheric levels. The fact that the upper boundary cuts through the levels of maximum forcing does indeed not seriously affect our results since wave drag at these altitudes near the stratopause has only a negligible impact on the surface. Finally, Figs. 2c and 2d exemplify those components of the NAM index during the same winter, which are due to either stratospheric or tropospheric wave drag, respectively. This demonstrates the successful separation of these two components of wintertime NAM variability as obtained by our method.

## 3. Results

In the following, the NAM surface signal associated with stratospheric variability, as derived from ERA-40, is first illustrated and briefly discussed. Subsequently, the method, introduced in the previous section, is applied to deduce the respective contributions of stratospheric versus tropospheric forcing, and, finally, the robustness of the obtained results is demonstrated through a sensitivity study with respect to our assumptions on damping time scales and static stability.

### a. Northern annular mode surface signal associated with stratospheric variability

The wintertime NAM is investigated by means of a lagged linear regression analysis. Specifically, the NAM index at each individual level, from the surface to the stratopause, is regressed onto the 10-hPa midstratospheric NAM index time series. Here, we choose to regress onto the negative 10-hPa time series to ease comparison with related studies that investigate stratosphere–troposphere interactions during negative stratospheric NAM episodes, associated with sudden stratospheric warmings (Charlton and Polvani 2007) or with weak vortex events (Baldwin and Dunkerton 2001), although due to linearity the time–height structure of the results do not depend on the sign. This analysis is similar to that of Thompson et al. (2006), with the exception that here we focus on the NAM itself rather than on the zonal-mean zonal wind within a given latitude band.

The results of our analysis are presented in Fig. 3a. (By construction, this corresponds to both, the ERA-40 NAM index as well as the NAM index from the case forced at all levels.) First, the signature of a stratospheric negative NAM anomaly that progresses downward from the upper into the lower stratosphere on a time scale of several weeks is evident. This downward progression is a feature that is well known from composites of sudden stratospheric warmings and of weak and strong vortex events, although the time scales of such events are observed to differ from what we find here. As shown by the analysis of Limpasuvan et al. (2005), sudden stratospheric warmings are associated with downward progression time scales of only 1–2 weeks, whereas for so-called polar vortex intensification events it takes up to 8 weeks to progress from the upper into the lowermost stratosphere. Since our linear regression analysis includes both positive and negative NAM anomalies, our results must represent an average signature of strong and weak vortex events. Furthermore, in contrast to composites based on only a few and large events, the regression analysis is based on all times in the series, and it is expected that small anomalies of either sign descend on some intermediate time scale. Thus, the analysis presented here combines different types of stratospheric events, and this should be kept in mind for interpretation of the results. Second, in addition to this stratospheric signature, there is a signal of the same sign at tropospheric levels that maximizes near the surface at positive lags. Note, at this point, that due to the separate normalization of the NAM index at each level, effects of the density stratification of the atmosphere are eliminated from this analysis. Thus, the smallness of the surface signal compared to that at stratospheric levels is only a consequence of the fact that most tropospheric NAM variability—driven by internal dynamical processes—is unrelated to the stratosphere.

The magnitude and time evolution of the surface signal are emphasized in Fig. 3b (thick line). The signal is found to be statistically significant at the 99% level over a period of about 8 weeks, lasting from lag −8 until lag 47 (with lags in units of days), and attains its peak amplitude at lag 19, that is, with a delay of about 3 weeks relative to the midstratopheric NAM.^{3} The peak amplitude of 0.22 standard deviations implies an explained variance of 5% if the full surface NAM index time series is considered, which also contains variability at much higher frequencies than that of the surface signal. If only low-frequency NAM surface variability (with periods > 30 days) is considered, the fraction of surface variance explained by the midstratosphere increases to 10% [for related results, see also Baldwin et al. (2003) and Charlton et al. (2003)]. For the remaining analysis we, nevertheless, proceed with the unfiltered NAM index since we are interested in the dynamical evolution and some dynamical features may be lost, or smoothed out, when using the low-frequency time series.

### b. Contribution of stratospheric versus tropospheric forcing

We now consider the cases with the forcing restricted to either stratospheric or tropospheric levels, each case of which is leaving a characteristic NAM surface signal, also shown in Fig. 3b (thin gray and thin black line, respectively). The surface signal in response to stratospheric wave drag exhibits only one phase, with significantly negative values from lag −18 to lag 25, whereas the signal in response to tropospheric wave drag has a two-phase character, with significantly positive values from lag −14 to lag −3 and negative values from lag 18 to lag 46. Hence, the role of the first phase of the response to tropospheric wave drag is to largely cancel out, until lag −8, the surface signal due to stratospheric wave drag. The role of the second phase, by contrast, is to prolong the signal by 3 weeks beyond the response to stratospheric wave drag.^{4} This indicates that wave drag local to the troposphere is indeed necessary to accomplish the delay of several weeks between the NAM surface signal and associated stratospheric variability since, otherwise, step 1 of our method would not have produced the respective tropospheric forcing, and it is precisely this delay by which the stratosphere may add any intraseasonal predictability potential to the NAM at the surface. Nevertheless, aside from the wrong timing, the amplitude of the surface signal in response to stratospheric wave drag compares well with the full NAM surface signal.

To further investigate the dynamics behind the NAM surface signal, we compute the tendency of the NAM index and of its components due either to wave drag, surface drag, or thermal damping. These tendency components correspond respectively to the first, second, and third additive term on the rhs of the spectral tendency equation (20). The same lagged linear regression analysis onto the 10-hPa midstratospheric NAM index is then performed for these quantities, as described above for the NAM index itself. The result is illustrated in Fig. 4. The NAM surface tendency due to tropospheric wave drag (Fig. 4e and red line in Fig. 4n) exhibits the same two-phase character as discussed above for the corresponding NAM surface signal and, thus, first acts to delay the buildup of this signal at negative lags, whereas at positive lags it helps to maintain and persist the signal over several weaks. The tendency due to wave drag is largely counterbalanced by the effect of surface drag (Fig. 4h and green line in Fig. 4n), leading to a positive NAM tendency until lag −7 and to negative tendencies thereafter until about lag 30 (Fig. 4b and black line in Fig. 4n).

By contrast, the NAM surface tendency due to stratospheric wave drag (Fig. 4f and red line in Fig. 4o) is exclusively negative and shows a gradual increase in magnitude from lag −40 toward lag −5, followed by a sudden drop across lag 0 to about 20% of its peak value, and then further reduces and is near zero after lag 20. This sudden drop in stratospheric wave driving can be explained by nonlinear wave–mean flow interactions in the stratosphere and is likely to reflect the occurrence of sudden stratospheric warming events, which are associated with a complete breakdown of the polar vortex, indicating highly nonlinear wave dynamics. This sudden drop, together with the effect of surface drag (Fig. 4i and green line in Fig. 4o) acting on the previously generated tropospheric flow anomaly, leads to a switch from negative to positive NAM surface tendencies (black line in Fig. 4o) across lag zero, and only positive tendencies thereafter. Hence, the downward effect of the secondary circulation induced by QG adjustment to stratospheric forcing cannot actively contribute to the late NAM surface signal. Note, however, that there is a significant contribution due to QG adjustment to stratospheric radiative damping (blue line in Fig. 4o), which produces negative surface tendencies at both negative and positive lags. At this point, it is noteworthy that, as seen in Figs. 4f and 4l, respectively, the effect of the secondary circulation induced by stratospheric wave drag is to spread the corresponding negative NAM anomalies downward to the surface, leading to anomalies of the same sign there, whereas the secondary circulation due to the concurrent thermal relaxation of the anomalously warm polar vortex induces positive tendencies at stratospheric levels and tendencies of opposite sign at the surface. This mechanism by which stratospheric radiative damping contributes to NAM surface tendencies has been nicely illustrated by Thompson et al. (2006, see their Fig. 3). However, this contribution is found to be insufficient to create negative NAM surface tendencies against surface drag. Since, for the key scenario presented here, the thermal damping time scale at the surface is set to 40 days and thus thermal relaxation in the troposphere is rather weak, the according overall surface tendencies due to thermal damping (blue line in Fig. 4m) are largely dominated by the contribution of radiative damping in the stratosphere.

Regarding the effect of wave drag local to the troposphere it is important to note that this contribution cannot be viewed as due to a completely independent process from stratospheric wave drag. Since virtually all wave drag at stratospheric levels is generated by upward-propagating planetary (and partly gravity) waves of tropospheric origin, there must, for reasons of momentum conservation, always be an oppositely signed wave drag signature in the troposphere, although this signature may partly project onto other modes than the NAM [because waves may also propagate horizontally within the troposphere before entering the stratosphere, as shown by Thompson et al. (2006, see their Fig. 9)]. The positive surface tendencies due to tropospheric wave drag at negative lags (Fig. 4e and red line in Fig. 4n), therefore, may be expected, to some extent, to represent the tropospheric tail end of the same wave propagation process that leads to oppositely signed tendencies at stratospheric levels (Fig. 4f) and, by QG adjustment, also at the surface (red line in Fig. 4o). On the other hand, the negative surface tendencies due to tropospheric wave drag at positive lags must be largely due to internal tropospheric wave dynamics since there is no similar anticorrelation with the corresponding surface tendencies in response to stratospheric wave drag. An additional integration of the spectral tendency equation (20) with the forcing further restricted to only the lower half of the stratosphere, between the 100- and the 10-hPa level (not shown), indicates that those near-zero surface tendencies at positive lags are, indeed, not the result of a cancelation between the effect from the upper versus the lower stratosphere, as one may want to anticipate from the stratospheric tendencies shown in Fig. 4f. Thus, the downward effect of QG adjustment to stratospheric wave drag at positive lags appears, in fact, to be weak.

### c. Sensitivity to damping time scales and static stability

The robustness of the above findings is now investigated by variation of the model parameters that control damping time scales and static stability. Each parameter configuration is consistently used in both steps 1 and 2 of our method, described in section 2c, since otherwise the second step would not reconstruct the true ERA-40 NAM index. Thus, this sensitivity study does not investigate the response of the annular mode to changes in those parameters under the same external forcing but, rather, investigates how the stratospheric and tropospheric forcings, implied by step 1 of the method, do actually change if we modify our assumptions regarding the values of damping time scales or static stability.

When the surface thermal damping coefficient is increased from *α _{s}* = (40 days)

^{−1}to

*α*= (5 days)

_{s}^{−1}, then tropospheric flow anomalies driven by the downward effect of QG adjustment to stratospheric forcing are damped out more efficiently. Thus, the NAM surface signal in response to stratospheric wave drag is reduced in this scenario and is significant only until lag 7 (see Fig. 5a; cf. Fig. 3b), beyond which the effect of tropospheric wave drag is now already needed to persist the surface signal. Hence, our above finding that wave drag local to the troposphere is necessary for the persistent surface signal is even further corroborated the shorter we assume the surface thermal damping time scale. The surface tendency components shown in Fig. 6 illustrate how the balance changes in this scenario. The negative peak value of the NAM surface tendency due to stratospheric forcing (black line in Fig. 6c) reduces by about 50% as a consequence of the additional positive tendency due to tropospheric thermal damping, compared to the key scenario (see Fig. 4o). Accordingly, the surface tendency due to tropospheric wave drag (red line in Fig. 6b) increases in magnitude to balance the stronger thermal damping as well as to account for the reduced NAM surface tendencies due to stratospheric forcing, such that the full NAM surface tendency (black line in Fig. 6a, which, by construction, is identical to the black line in Fig. 4m) is recovered.

As an alternative to changing the thermal damping coefficient only at the surface, we may change the overall structure of its vertical profile. As an extreme test we simply remove all vertical structure by choosing a height independent thermal damping coefficient according to *α* = const = (20 days)^{−1}. The results (not shown) are very similar to those from the key scenario. The only noteworthy exception is that the fraction of the NAM surface signal in response to stratospheric wave drag that comes from the lower stratosphere (between 100 and 10 hPa), as compared to that from the upper stratosphere, increases since thermal damping is strengthened at lower- and weakened at upper-stratospheric levels.

When the surface drag coefficient is increased from *κ _{s}* = (1 day)

^{−1}to

*κ*= (0.5 day)

_{s}^{−1}, then the results (not shown) are similar to those from the scenario with increased

*α*, although the changes are weaker since

_{s}*κ*is increased by only a factor of 2. When, on the other hand, surface drag is reduced, then wave drag local to the troposphere becomes less important in maintaining and persisting the NAM surface signal, and with

_{s}*κ*= (3 days)

_{s}^{−1}its contribution vanishes completely in the sense that the corresponding surface signal (thin black line in Fig. 5b) is only positive at positive lags. Note, however, that this is an unrealistic value for a surface drag coefficient. Atmospheric circulation models that employ a Rayleigh drag term to mimic the effect of small-scale momentum fluxes within the planetary boundary layer generally produce realistic climatologies for surface drag time scales not longer than 1 day [see also the benchmark calculation by Held and Suarez (1994)]. Furthermore, here we use

*α*= (40 days)

_{s}^{−1}, which represents an upper bound for lower tropospheric thermal damping time scales since this value corresponds to lower stratospheric radiative damping times (see Newman and Rosenfield 1997; Hitchcock et al. 2010) and, therefore, does not account for the additional effects of small-scale turbulent heat fluxes and of higher temperatures at lower tropospheric levels. Hence, the role of tropospheric wave drag for maintaining the NAM surface signal vanishes only if this upper bound is combined with unrealistically weak surface drag.

Finally, the impact of our assumption for the static stability profile is investigated. As an extreme test we set *N*^{2} = const such as to correspond to an isothermal atmosphere at *T*_{0} = 240 K. Since this scenario (not shown) implies an increased static stability throughout the troposphere, the downward penetration into the troposphere of the NAM signal induced by QG adjustment to stratospheric forcing is less efficient and, thus, the role of tropospheric wave drag for persisting the NAM surface signal is again further emphasized, as was the case in the above scenarios with increased surface damping time scales.

## 4. Discussion

For interpretation of the results presented in the previous section, it is important to note that the ERA-40 in the upper stratosphere is based on only few observations, in particular before the end of the 1970s when satellite measurements were not included. Since our analysis indicates that part of the NAM surface signal in response to stratospheric forcing comes from the upper half of the stratosphere, one may speculate to what extent our results might be influenced by any misrepresentation of the upper-stratospheric circulation in the reanalysis. On the other hand, since we are concerned only with the zonally averaged circulation, such errors are possibly relatively small compared to a situation when the detailed two-dimensional flow during individual sudden stratospheric warming events were considered. Moreover, even if the ERA-40 NAM index at upper-stratospheric levels exhibits significant errors during some episodes, the linear relationship between the stratosphere and the troposphere, as revealed by our lagged regression analysis, may still be little affected in a statistical sense as long as the underlying dynamical processes are well represented by the circulation model that has been used for the reanalysis.

Regarding the dynamical mechanisms behind the NAM surface signal associated with stratospheric variability, it is interesting to note, from Fig. 4a, that the short episodes in the troposphere with negative NAM tendencies during the onset (lag −10 to lag 0) and positive tendencies during the decline (lag 40 to lag 55) of the surface signal occur both at the end of longer lasting episodes with tendencies of the same sign in the stratosphere. Despite the apparent symmetry of these two episodes, however, the results presented in the previous section suggest that the onset is associated with planetary wave propagation processes between the troposphere and the stratosphere, whereas the decline is likely to be caused by the end of a period with anomalous wave drag due to internal tropospheric dynamics, although possibly in response to lowermost stratospheric NAM anomalies. At this point, however, dynamical interpretation is limited since our approach does not allow us to distinguish between different classes of wave processes, but only accounts for their overall effect on the zonally averaged circulation, indicating eddy feedbacks on the zonal flow.

For further interpretation of our results it is of interest to relate the approach used here to those employed in other studies, which also investigate stratospheric impacts on the troposphere induced by QG adjustment. Black (2002), for example, who applies the QG potential vorticity inversion approach of Hartley et al. (1998), diagnoses the tropospheric wind field that is associated with stratospheric potential vorticity anomalies. This is, in fact, closely related to but not the same as our approach, which implicitly also performs an inversion of QG potential vorticity anomalies, although not entirely restricted to stratospheric potential vorticity. To make this explicit, it is noteworthy that the tendency equation of the zonal mean QG potential vorticity , defined by (see, e.g., Andrews et al. 1987)

can be easily obtained by deriving an equation analogous to Eq. (6) in terms of Φ, rather than *u*, and multiplication by the Coriolis parameter, as

with the modified variables , , and . From Eq. (25) it is obvious that the production of potential vorticity by the external mechanical forcing *F* is local in *z* since this term involves only a horizontal derivative. Thus, in the case forced at stratospheric levels, potential vorticity anomalies produced directly by the external mechanical forcing are also restricted to the stratosphere. Likewise, potential vorticity anomalies produced by surface drag are restricted to the boundary layer and, therefore, cannot contribute directly to stratospheric potential vorticity tendencies. Since, on the other hand, the thermal damping term in Eq. (25) includes a vertical derivative, thermal damping just below the cutoff height *z _{c}* could in principle lead to a direct contribution to stratospheric potential vorticity tendencies just above, through a change in static stability and, thus, the production of stretching vorticity—although this effect might be weak owing to the long thermal damping times in the lowermost stratosphere. However, both surface drag and tropospheric thermal damping can have an indirect upward effect on stratospheric potential vorticity tendencies. This indirect effect arises from lower-stratospheric thermal damping that acts on previously generated temperature anomalies driven by the upward secondary circulations induced by the tropospheric damping.

Hence, the stratospheric potential vorticity anomalies inverted by Black (2002) do not include the direct effect of surface drag and probably not much of a direct effect of tropospheric thermal damping, but may include their weak indirect effects. The same is true for the results of Hartley et al. (1998) and Ambaum and Hoskins (2002),^{5} which are also based on the inversion of purely stratospheric potential vorticity anomalies; in general, similar arguments apply to the study of Charlton et al. (2005), although they use a more accurate Ertel potential vorticity inverter which is not restricted to QG scaling. The approach of Thompson et al. (2006), on the other hand, is conceptually similar to our case forced at stratospheric levels. However, their surface drag and thermal damping are specified from observed anomalies and, thus, include the damping that acts on anomalies driven by tropospheric wave drag although their model is forced only at stratospheric levels (this would correspond to combining the red line in our Fig. 4o with the green and blue lines from Fig. 4m). We assume that this may partly explain their different results.

In summary, the above-mentioned approaches based on the inversion of stratospheric potential vorticity anomalies must, by construction, largely ignore the effects of surface drag and thermal damping acting on the associated tropospheric wind field. As a consequence, the associated surface signal always tends to be proportional to (and, thus, to be in phase with) the stratospheric potential vorticity anomalies above and, accordingly, does not exhibit separate phases of initiation and maintenance. When tropospheric dissipative processes are included, however, as in Thompson et al. (2006) and in the present study, it turns out that the surface signal needs to be actively maintained against this dissipation in order to obtain a lagged surface signal relative to stratospheric variability. Or, in terms of potential vorticity, the tropospheric damping produces tropospheric potential vorticity anomalies that oppose the downward effect of the stratospheric potential vorticity anomalies above; thus, in order to maintain the surface signal against the effects of dissipation, a mechanism is needed that removes, or at least reduces, those anomalies at tropospheric levels.

Finally, the dynamical justification for neglecting the external heating term *Q*_{*} in Eq. (20) is discussed. Since the squared horizontal and vertical wavenumbers −*ε* and −*γ* are all positive (and, thus, *λ* is positive) and also increase monotonically with their modal index (see Tables 1 and 2 for the first few eigenvalues), it turns out that the terms and , as defined by Eq. (22), are always positive and smaller than one—which precisely represents the counteracting effect of the secondary circulation against any forcing applied to the zonal mean atmosphere. Owing to the opposite scale dependence, however, of this counteracting effect for mechanical versus thermal forcing, it can be deduced from Eqs. (20) and (22) that the external mechanical forcing is most efficient in driving deep modes with small horizontal scales, which are subject to inefficient thermal damping, whereas the external heating , when retained in Eq. (20) by the additional term , is most efficient in driving shallow modes with large horizontal scales, which are subject to efficient thermal damping. Thus, as long as thermal damping dominates over mechanical damping, which is true at least for the stratosphere, wave drag will have a tendency to be more efficient than external heating in perturbing zonal wind anomalies, which generally reduces the importance of external heating anomalies for intraseasonal NAM variability. Nevertheless, external heating anomalies may become important, for example, in late winter and early spring when positive feedbacks between the strength of the polar vortex and ozone destruction in the presence of polar stratospheric clouds play a role (e.g., Austin et al. 1992), and the effects of such feedbacks are, therefore, not accounted for by our approach.

## 5. Conclusions

In this study, the wintertime NAM surface signal associated with stratospheric variability, as derived from ERA-40, is investigated by means of a lagged linear regression analysis, with reference to the midstratospheric daily NAM. The main findings are as follows.

In the troposphere a statistically significant NAM signal occurs, which maximizes at the surface and lasts for about 8 weeks, consistent with findings of previous studies (e.g., Baldwin and Dunkerton 2001).

The peak time of the NAM surface signal lags the midstratosphere by about 3 weeks, similar to the results of previous studies (e.g., Thompson et al. 2006).

The explained variance of the surface signal amounts to 5% of daily NAM surface variability, and to 10% of low-frequency (periods > 30 days) NAM surface variability [related results were obtained in previous studies by, e.g., Baldwin et al. (2003) and Charlton et al. (2003)].

To estimate the relative contributions to this surface signal, of (i) the downward effect of QG adjustment to stratospheric wave drag and radiative damping and (ii) of wave drag local to the troposphere, a method is introduced to obtain, in a first step, the external mechanical forcing that, in the QG framework, drives precisely the ERA-40 stratospheric and tropospheric daily NAM. This is accomplished by inversion of an elliptic partial differential equation for the zonal-mean zonal wind tendency, where the latter is given by the daily NAM time series, and where the inversion is achieved via a spectral tendency equation that is derived following the concepts of tidal theory. In a second step, this spectral tendency equation is then integrated in time to reconstruct the daily NAM, but with the forcing, interpreted as wave drag, restricted to either stratospheric or tropospheric levels, each case leaving a characteristic NAM surface signal. From the application of this method we can conclude the following:

The NAM surface signal that is due to QG adjustment to stratospheric wave drag and radiative damping has an amplitude that compares well with the full surface signal, but peaks 3 weeks earlier with no lag relative to the midstratosphere.

The NAM surface signal that is due to wave drag local to the troposphere has a two-phase character. At negative lags, relative to the midstratosphere, its role is to counteract the effect of stratospheric forcing and, thus, to delay the buildup of the surface signal, whereas at positive lags it actively maintains and persists the surface signal over several weeks and is responsible for the lag of 3 weeks relative to stratospheric variability.

The tropospheric wave drag during the first phase that delays the buildup at negative lags, may be expected, to some extent, to be associated with planetary wave propagation between the troposphere and stratosphere and, thus, to be intimately linked to the concurrent and oppositely signed stratospheric wave drag anomalies and their corresponding effect on the NAM surface signal. Further research, however, is needed to confirm this conclusion.

By contrast, the tropospheric wave drag during the second phase at positive lags that helps to maintain the surface signal results from wave dynamics internal to the troposphere since it is found to be unrelated to stratospheric wave drag.

This indicates that tropospheric eddy feedbacks are needed to achieve the prolonged NAM surface signal and its delay relative to stratospheric variability.

The downward effect of QG adjustment to stratospheric radiative damping alone is found to be of only minor importance in the sense that it is not sufficient to even compensate the oppositely signed tendencies due to surface drag that acts on the previously generated NAM surface anomaly in response to stratospheric wave drag.

These conclusions are found to be robust, in a qualitative sense, against changes of model parameters that control damping time scales and static stability, as long as values are varied within realistic bounds. Thus, our conclusions do not depend on the details of the specific model configuration. Exact quantitative statements, however, are difficult to achieve with our approach owing to the involved simplifications, which are reflected in the assumptions regarding model parameters, vertical boundary conditions, and the conceptual model setup.

The above results help to reconcile the conclusions from previous studies. While some studies (e.g., Charlton et al. 2005) claim that the QG adjustment of the troposphere to stratospheric anomalies is inadequate to explain the tropospheric response to stratospheric variability and that internal tropospheric eddy feedbacks are needed, other studies, most explicitly Thompson et al. (2006), conclude that such eddy feedbacks are not needed and the surface signal associated with stratospheric variability can be fully explained by the QG adjustment mechanism. Our results, however, indicate that both types of mechanisms are equally relevant to the observed connection between the stratosphere and the troposphere. Whereas the QG adjustment mechanism is responsible for the initiation of the NAM surface signal, internal tropospheric eddy feedbacks are needed to explain the time lag of the surface signal relative to the stratosphere. Since it is this time lag that is relevant in the intraseasonal prediction context, further research should focus on such eddy feedbacks, whether they are due to high-frequency synoptic-scale processes or due to slow variations of quasi-stationary Rossby waves, or both, and whether they occur in response to the initial tropospheric signal or directly to lower-stratospheric anomalies above the tropopause.

## Acknowledgments

TK and RJG are grateful to support from GEOMAR during the time this work was carried out. The authors also thank two anonymous reviewers for their helpful comments and Mojib Latif for discussions of the results.

### APPENDIX

#### Boundary Conditions

##### a. Horizontal boundary conditions

The boundary conditions to the horizontal structure equation (13) require at the poles, which implies

##### b. Vertical boundary conditions

###### 1) Mechanical forcing

The upper-boundary condition to the vertical structure equation (14) requires that the log–pressure coordinate vertical velocity vanishes at *p* = 0. Here, we apply this condition at the top of the model domain *z _{t}*, which can be shown to be equivalent to posing

based on the additional assumption that = at *z* = *z _{t}*. Since the model top is located at

*z*= 48 km, this is only a crude approximation to the true upper boundary condition. As a test, we have tried more advanced configurations of the upper part of the model domain, with the model top located several tens of scale heights above the stratopause, tapering

_{t}*N*

^{2}smoothly to zero across this interval [such as to meet the criteria for obtaining a discrete vertical eigenvalue spectrum, according to Cohn and Dee (1989)] and using constant extrapolation of the NAM pattern beyond the stratopause. However, the impact of such changes, although important for the upper stratosphere, was found to be negligible for the NAM surface response to tropospheric as well as stratospheric forcing. Therefore, we simply use Eq. (A2) as the upper boundary condition at

*z*= 48 km.

_{t}The lower boundary condition states that there is no mass flux across the surface [see Andrews et al. (1987), assuming a flat bottom],

with *T*_{00} = *T*_{0}(*z* = 0). In case of purely mechanical forcing (*Q*_{*} = *α* = 0), using Eqs. (3) and (4), the lower boundary condition is then obtained as

Here, again, we assume that = at *z* = 0, although eddy heat fluxes are generally nonzero at the surface and, thus,

where the primes indicate deviations from the zonal mean. The above assumption would imply wrong values of at *z* = 0 and, thus, of at and above the surface, with the consequence of incorrect at and above the surface. The corresponding additional Coriolis torque above the surface will be taken up by the forcing *F*, obtained from step 1 of our method described in section 2c. Since, however, in step 2 the same lower boundary condition is applied, the correct NAM is, nevertheless, obtained as long as *F* is prescribed to Eq. (20) at all levels. When *F* is restricted, however, to either stratospheric or tropospheric levels, the result might indeed be influenced if the additional forcing *F* that takes up the additional Coriolis torque is significant both below and above the cutoff height at *z _{c}* = 16 km. Then, the separation between stratospheric and tropospheric forcing would indeed be influenced by the assumption that = at

*z*= 0. The following argument, however, indicates that this impact is small above

*z*= 16 km, and therefore, the lower boundary condition (A4) is applicable in the context of the present study.

_{c}The time series of the eddy flux term (A5) is computed from ERA-40 and regressed onto the surface NAM index, and the resulting regression pattern is then interpolated onto the horizontal model grid and transformed into spectral representation. The majority of the variance of the term (A5) is found to project onto horizontal mode *n* = 7 with some smaller fraction still contained on horizontal mode *n* = 3, and the undamped atmospheric responses in these horizontal modes, if forced at a given level, are shown [by Plumb (1982), see his Eq. (4.4) and Table 1] to have upward vertical decay scales of only 4 and 12 km, respectively. Thus, the upward vertical extent of the secondary circulation induced by an eddy heat flux forcing at the surface according to Eq. (A5) is sufficiently less than *z _{c}* = 16 km.

Finally, as an alternative to Eq. (A4), we also tested the more simplified boundary condition at *z* = 0 [for a detailed discussion of the effects of this simplification, see Haynes and Shepherd (1989)], but this was found to have a significant impact on the results, even on the NAM surface response to forcing applied only at stratospheric levels. The surface response, in some cases, amplifies by more than a factor of 2, and also the relative roles of the various terms analyzed in section 3 are changed. Therefore, we apply Eq. (A4) at the lower boundary.

###### 2) Thermal forcing

In analogy to the HSE (13) and VSE (14) for the zonal wind, an equivalent pair of eigenvalue equations can be obtained for the vertical velocity. With the separation according to *w**(*φ*, *z*, *t*) ∝ exp(−*σt*)Θ(*φ*)*W*(*z*), where , the corresponding HSE is obtained as

with the horizontal structure functions (which satisfy at the poles), and the VSE may be written as

with the vertical structure functions . The additional dependence on the horizontal modal index *n* results from the lower boundary condition, obtained in case of purely thermal forcing (*F* = *κ* = 0) from Eq. (A3), also using Eqs. (1), (2), (5), and (A6), as

which depends on the horizontal eigenvalue *ε*^{(n)}.

The upper boundary condition is identical to the case of purely mechanical forcing, given by Eq. (A2), and in terms of *W* takes the form

The vertical structure functions in terms of zonal wind *u* are then given by the relation

## REFERENCES

*Middle Atmosphere Dynamics.*Academic Press, 489 pp.

_{2}climate

*Atmospheric Tides*. D. Reidel, 200 pp.

*U.S. Standard Atmosphere, 1976*. NOAA, 227 pp.

**A**

## Footnotes

^{1}

It might be instructive to note that an analogous but simplified version of Eq. (6), with *F* = *Q*_{*} = 0, appears in the appendix of Scott and Haynes (1998), see their Eq. (A.5), for a Boussinesq atmosphere in *f*-plane geometry and with uniform static stability *N*^{2}.

^{2}

We use a slightly smoothed version of the U.S. Standard Atmosphere to avoid any discontinuities in static stability that, in a time and zonal mean sense, are rather unrealistic. Technically, the method works equally well without smoothing, although this has the effect to focus the forcing, obtained from step 1, to those discontinuities located at the interface levels between the various layers of the standard atmosphere.

^{3}

Statistical significance is estimated by a bootstrap approach in Fourier space (see Ebisuzaki 1997). Specifically, the surface time series during each individual season is resampled by randomizing the phases of the corresponding Fourier coefficients, and is then regressed onto the time series at 10 hPa. This procedure is repeated for 1000 times and the upper and lower 0.005% quantiles are then used as significance thresholds.

^{4}

Note that a hint at this two-phase character was implicitly contained already in the results of Black (2002); compare his Figs. 7a and 7b. In that analysis, upper-tropospheric potential vorticity anomalies also had the effect to first reduce and later to amplify the direct tropospheric response to stratospheric variability.

^{5}

These authors speculate that the omission of surface drag in their simple model may partly explain the quantitative mismatch between their model results and observations.