1. Introduction

It is by now widely accepted that short-term variability of the extratropical stratosphere can arise from variability external to the stratosphere, in particular from variability of the tropospheric circulation, a principle source of planetary-scale wave forcing. The most dramatic example of such short-term stratospheric variability is the so-called sudden warming, during which the radiatively driven cyclonic circulation of the winter stratospheric polar vortex undergoes a rapid wave-induced deceleration, in extreme cases even reversing to easterly flow, and polar temperatures increase by several tens of degrees in the space of a few days. Numerous previous numerical and observational studies, dating back to the pioneering work of Matsuno (1971), have linked these events to anomalously large tropospheric wave amplitudes; see, for example, Andrews et al. (1987) and, more recently, the review by Haynes (2005).

Another way in which stratospheric variability can arise is through the action of dynamical mechanisms internal to the stratosphere itself. Studies dating back to Holton and Mass (1976) have shown how vacillating states can arise in zonally and meridionally truncated quasigeostrophic channel models as periodic solutions to the equations of motion. The main requirement of these models is the competing actions of thermal relaxation to a zonally symmetric basic state together with a planetary-scale wave forcing at the lower boundary. Later studies explored the parameter space defined by the radiative and wave forcing amplitudes and described the transition to vacillating states variously in terms of a Hopf bifurcation from a steady solution branch (e.g., Yoden, 1987), or a catastrophe-like transition from one stable solution branch to another (Chao 1985). Christiansen (2000a) demonstrated the existence of further bifurcations leading ultimately to chaotic solutions of the system.

Extending this line of investigation, Scott and Haynes (2000) considered the additional effect of more complete meridional structure in a primitive equation model, yet still within the context of a zonally truncated system, and showed that internal variability persisted over a wide range of forcing amplitudes. Variability has also been found in models with isotropic horizontal resolution. Pawson et al. (1995) found variability in a controlled general circulation model (GCM) simulation, although the character of the variability was not coherent and had a frequency spectrum similar to red noise. The variability obtained by Christiansen (1999, 2000b) in a different GCM, on the other hand, showed more coherence. Even in the latter study however, the variability is not particularly regular in time, and there is a marked decrease in regularity as horizontal resolution is increased from T21 to T42. In a shallow water model of a forced polar vortex, Rong and Waugh (2004) also found a reduction in regularity as horizontal resolution was increased from T42 to T85. Given these ambiguities it remains inconclusive whether coherent variability is an inherent feature of the dynamics of the stratosphere.

Our intent, therefore, is to isolate the dynamical processes involved in internal stratospheric variability in as comprehensive a model as possible that nevertheless allows a simple interpretation of the model evolution. Since we are interested in internal stratospheric variability the natural choice is a GCM dynamical core, run in stratosphere-only mode and with a minimal number of extra physical processes included. Such a model has the additional advantage of permitting a thorough exploration of parameter space. The work presented here extends the preliminary study of Scott and Polvani (2004), which illustrated coherent variability in a GCM dynamical core with a constrained troposphere. We note that, in contrast to recent work by Taguchi et al. (2001), we are not concerned here with the stratosphere–troposphere coupled system, but are rather interested in carefully isolating the dynamical processes within the stratosphere itself. In a subsequent paper, we will address the issue of how these processes are modified under the influence of time-dependent tropospheric forcing.

A major concern and motivation for the present work lies in the fact that variability in dynamical systems often loses coherence as the number of degrees of freedom of the system increases. This was illustrated by Cehelsky and Tung (1987) for the case of multiple flow regimes in the troposphere. The possibility therefore exists that the variability obtained in the Holton–Mass-type models arises purely as an artifact of the severe numerical truncation. Our earlier work (Scott and Polvani 2004) provides a first indication that this is not so in the case of the winter stratosphere, but that coherent (i.e., regular or quasi periodic) variability persists at higher (isotropic) resolution. The present paper extends that work and considers robustness and sensitivity of variability, as well as the dynamical mechanisms involved.

A further motivation for the present work arises from recent studies indicating the possible influence of the stratospheric circulation on both the short and long time scale evolution of the troposphere (e.g., Baldwin and Dunkerton 1999; Polvani and Kushner 2002). For example, if all stratospheric variability is generated by mechanisms whose origin lies entirely within the troposphere, then the stratosphere in itself can never be considered as having any influence on the tropospheric circulation. In this case, it can only act as a mediator of tropospheric variability, though perhaps having an integrating effect in time (e.g., Norton 2003). On the other hand, the existence of internal stratospheric variability raises the possibility of actual downward influence, whereby the stratosphere itself exerts a direct influence on the tropospheric circulation. Finally, an important issue in climate change is the extent to which internal stratospheric variability depends on long-term changes in the radiative forcing of the stratosphere.

This paper presents a detailed account of the persistence of coherent internal variability at a resolution that approaches that which is able to provide a reasonable representation of the dynamical processes involved, and which is comparable to the resolution of current generation comprehensive general circulation models. We present results over a wide range of physical and numerical parameters and find that coherent variability exists under general conditions. Further, we present a detailed description of the vortex structure and dynamical process involved in the variability. Diagnostics based on zonal means, wave fluxes, and three-dimensional potential vorticity fields provide a consistent picture of the variability arising from the competing actions of radiative cooling and wave forcing. In the companion paper we consider how the internal variability is modified when variability is introduced into the external forcing through a simple modulation of the lower boundary wave forcing.

The outline of the paper is as follows. In section 2, we describe our approach and numerical methods. In section 3, we discuss an example of the variability obtained with a particular set of forcing parameters and describe the dynamical processes and the time-dependent evolution of the vortex structure. In section 4, we examine in detail the dependence of the variability on the main parameters defining the thermal and wave forcing. In section 5, we consider the robustness of the variability and its dependence on other physical and numerical parameters, including the radiative relaxation profile and forcing wavenumber, and the horizontal and vertical resolutions, and sponge layer. Section 6 concludes.

2. Method

In this work, our aim is not to reproduce exact features of the stratospheric circulation, and their dependency on details of thermal transfer, chemical process, gravity wave drag, and so on. Rather, we wish to isolate a fundamental dynamical property, namely the ability of the stratosphere to generate variability in itself, to examine the dependence of this property on the main dynamical forcings, and to demonstrate that this property is largely insensitive to other details, both physical and numerical. To achieve such a goal within a numerical modeling study, we need a model configuration that allows us complete control over the relevant dynamical properties.

We use a semispectral, pressure coordinate model to solve the dry primitive equations in a spherical domain. The model is the Built on Beowolf (BOB) dynamical core developed jointly at National Center for Atmospheric Research (NCAR) and Columbia University; numerical details, implementation, and availability are documented fully in Rivier et al. (2002) and Scott et al. (2004). Here it is used as a stratosphere-only model, with no dynamically active troposphere. There are two sources of radiative/dynamical forcing: a relaxation to a prescribed radiative temperature equilibrium, representing radiative cooling over the winter pole; and a wave forcing at the lower boundary representing wave generation in the troposphere and at the ground. Details of each follow.

The radiative forcing is a standard Newtonian cooling to a radiative equilibrium temperature distribution, T_e, that represents a perpetual January state with colder temperatures over the winter pole. The distribution T_e is obtained from a superposition of an isothermal atmosphere at a temperature T₀ = 240 K with a pressure dependent cold anomaly over winter pole T_PV:

View Expanded

where the dependence on latitude ϕ is given by

View Expanded

and where the winter polar vortex equilibrium temperature is defined as

View Expanded

In the above, ϕ₀ = 50°, δϕ = 10°, and p_b = 200 hPa. The parameter γ defines the strength of the associated polar vortex and takes the values 0.5, 1, and 1.5 in the integrations presented below. The time scale, τ at which temperatures are relaxed to T_e is uniform throughout the stratosphere at τ = 10 days. (Sensitivity of the results to τ is discussed in section 5.)

We discuss the wave forcing in the context of the lower boundary condition, which has the form of a specified, time-independent geopotential, Φ = Φ_s, of the type traditionally used in mechanistic modeling of the stratosphere. First, the zonal mean component, _s fixes the zonal mean velocity, u, at the lower boundary (independently of T_e, which only determines ∂u/∂z) and, in the calculations presented below, has the form

View Expanded

where a = 6371 km is the planetary radius, and Ω = 2π day⁻¹ is the planetary rotation rate. Through thermal wind balance, this form for _s gives a broad midlatitude jet with peak winds of approximately U₀ (neglecting the curvature term in the balance equation). We choose U₀ = 30 m s⁻¹ as representative of the winds in the lower winter stratosphere, but neither this value nor the form of Φ_s significantly affect the results; they are included here simply for reproducibility.

Stationary planetary waves are forced at the lower boundary by specifying the zonally varying component of the geopotential:

View Expanded

with

View Expanded

where g is the gravitational acceleration, θ is longitude, m is the zonal wavenumber of the forcing, and ϕ₁ = 30°. The function G attains its maximum value of one at ϕ = 60°. In the calculations presented below, we consider separately forcing wavenumbers m = 1 and m = 2 and forcing amplitudes h₀ between 0 and 1000 m, varied as an external parameter. Note that the wave forcing amplitude and phase in each calculation is completely time-independent, with the sole exception of an initial 10 day sin² ramp on Φ_s to initialize the integration smoothly.

Having specified the main forcings and boundary conditions, all other numerical and physical modifications to the dynamical equations are kept to a minimum. The model imposes a no-flux condition at the upper boundary p = p_t, which necessitates a sponge layer in the uppermost model levels to prevent spurious wave reflection. The sponge is restricted to above p_sp = 0.5 hPa, above the domain of interest, and damps the waves only at a rate k_sp max{0, [(p_sp − p)/p_sp]²}, with k_sp = 2(day)⁻¹. This is an important departure. Traditionally, sponge layer–type damping has also been imposed on the zonal mean velocity as a form of gravity wave drag parameterization used to close the strong westerly jet in the mesosphere. However, recent studies (Shepherd et al. 1996; Shepherd and Shaw 2004) have demonstrated that such a zonal mean drag term does not respect the angular momentum constraint of the total atmosphere, and as such may lead to unwanted feedback and downward influence effects. To avoid these unwanted side effects, and since the mesospheric winds are outside our immediate domain of interest, we omit the zonal mean drag altogether. As we have verified, however, doing so has almost no effect on the results; similar variability was obtained both with and without a drag on the zonal mean velocity (see section 5b).

To avoid the build up of enstrophy at small scales, a horizontal diffusion is included using a scale-selective ∇⁸ hyperdiffusion (uniform in pressure). A diffusion time scale of τ_h = 0.25 days at the smallest scales is used for most of the calculations presented below. The sensitivity to τ_h as well as to the parameters τ, k_sp, are further discussed in section 5.

The horizontal resolution used for the results presented in the main part of the paper is T42, equivalent to a grid spacing of approximately 2.8°. There are L + 1 half-levels in the vertical between the lower boundary at pressure p = p_b ≡ 200 hPa and the upper boundary at p = p_t ≡ 0.02 hPa. The half-levels are located at pressures p_l+1/2 = p_t(p_b/p_t)^l/L, l = 0, . . . , L, with L = 40. The L full-levels are located at the geometric mean of consecutive half-levels. The vertical domain thus extends from a notional tropopause at 12 km to approximately 75 km. The T42L40 resolution used allows us to perform a comprehensive exploration of the two-dimensional parameter space described below, but is high enough that we are confident that the results obtained are not artifacts of excessive truncation. As confirmation, we separately doubled the horizontal and vertical resolution for a single case, results for which are discussed in section 5b.

The model is initialized in a state of rest and run in the above perpetual January configuration for an integration period of 1000 days. Figure 1 shows the steady equilibrium states obtained for values of γ = 0.5, 1.0, 1.5 in the absence of wave forcing (i.e., h₀ = 0), and shows the direct relation of γ to the strength of the polar vortex obtained. Since the only force acting on u is the very weak (at these large scales) hyperdiffusion, the distributions shown are very close to those in exact thermal equilibrium with T_e.

Finally, we wish to emphasize again that all forcing in the model, both at the lower boundary and in the interior of the stratosphere, is time-independent. In addition, since the model lower boundary is at 200 hPa, the troposphere and its accompanying synoptic variability are, by construction, excluded from our calculations. Finally, since the zonal mean profiles shown in Fig. 1 are stable to shear instabilities, this potential source of variability does not exist in our model. Hence any time-dependence in the model must result from the internal dynamics of the stratosphere itself.

3. General characteristics

We begin by presenting and discussing in some detail an illustrative example, chosen to be representative of the stratospheric variability obtained with the model configuration outlined above. Figure 2 shows the zonal mean zonal velocity at 60°, u(z, t) as a function of height and time for wave forcing at wavenumber m = 1 (Fig. 2a) and m = 2 (Fig. 2b). For this particular case the wave forcing amplitude is h₀ = 600 m and the strength of radiative forcing is γ = 1. The character of variability obtained is typical of that found at other parameter values, comprising regular, but not perfectly periodic, vacillation cycles, in which the polar vortex strengthens systematically before reversing throughout the domain in a major sudden warming type event (here defined loosely as the reversal of u throughout a large part of the domain).

Insofar as the variability consists of a series of sudden warming events followed by a more gradual recovery of the vortex, the character of the variability in the two cases m = 1 and m = 2 is similar. In this zonal mean view, the main difference between the two cases is a lower frequency of the sudden warmings for m = 1, approximately one every 100 days compared with one every 60 days for m = 2. Further, for m = 1, the warmings are in general (though not always) less sudden than for m = 2, although in both cases we note that there is substantial variation in the strength and nature of sudden warmings; that is, the vacillations are not perfectly regular. We focus for now on the case m = 2 and note some fundamental properties of the variability; we return to details of the difference in variability between m = 1 and m = 2 in section 5 below.

First, an important observation needs to be made concerning the vertical flux of wave activity into the interior of the model domain: while the lower boundary condition wave forcing is constant in time, the flux of wave activity there is strongly time-dependent. This is illustrated in Fig. 3 for the simulation shown in Fig. 2b. Note the rapid increase in wave flux through the lower boundary that coincides with the sudden deceleration of the upper-level vortex. This effect was noted by Dunkerton et al. (1981), who discussed both the increase in upward wave flux at the lower boundary and the associated increase in westward phase tilt of the waves. We note further that the effect is not necessarily related to the presence of an artificial lower boundary and that similar variability in the wave flux was obtained by the authors in a model with a carefully controlled troposphere (Scott and Polvani 2004). The implication is that a sudden burst of tropospheric wave forcing is not a prerequisite for the sudden increases in wave flux into the stratosphere often accompanying major sudden warmings, but that in fact the upward wave flux into the stratosphere can be controlled by the stratosphere itself. An analytic expression for the vertical wave flux at the lower boundary was given in Esler and Scott [2005, see discussion following Eq. (27) therein] for the case of an axisymmetric vortex subject to a pulse of wave forcing.

Next, we consider in detail the dynamical evolution of the vortex during a complete sudden warming cycle. For definiteness we choose the time period from t = 660 up to the strong sudden warming occurring around t = 700 from the time series shown in Fig. 2b, for which h₀ = 600 and γ = 1. Although we focus on a single vacillation cycle, we emphasize that the dynamical mechanism we describe is quite general. As we show, the vacillation cycle can be considered as resulting from the competing effects of radiative cooling of the vortex and wave forcing, with the effectiveness of the latter being determined by the vertical structure of the vortex.

4. Dependence on forcing parameters

As already mentioned, the dynamical behavior described above is quite generally observed in the simulations carried out in this study. Rather than further document the details of the evolution of different sudden warmings, we proceed instead to examine how the general characteristics of the variability depend on the basic forcing properties of the system. We consider two model parameters as being fundamental to the behavior of the system, namely the radiative cooling, γ, and the wave forcing amplitude, h₀. Since these represent the main forcings of the system, we treat them separately and in more detail than the other physical parameters discussed in section 5 below. This approach is consistent with previous studies considering internal variability in highly truncated models.

In agreement with the results obtained in previous studies with highly truncated models, we find that for a given value of γ the system exhibits a whole range of behavior as the forcing amplitude h₀ is varied, from steady states at small h₀ to strongly time-dependent states at large h₀, such as those illustrated in Fig. 2. In fact, strong variability is found over a large range of wave forcing amplitudes and is generally characterized by recurring major warmings at intervals of approximately 40–100 days. To quantify this variability, we introduce some simple diagnostics based on the zonal mean zonal wind at ϕ = 60° and z = 40 km, u_jet(t), roughly coincident with the polar jet maximum in the upper stratosphere.

A simple, yet informative measure of the variability can be obtained by calculating the frequency spectrum of the time series of u_jet(t). Figure 8 shows the power in each Fourier mode, E(ω) = |û(ω)|² normalized by γ², for various values of γ and h₀ covering the parameter space of interest. The normalization accounts for the fact that strength of the zonal mean flow increases approximately linearly with increasing γ. We do not show E(ω) for h₀ ≤ 200 since all simulations quickly approach a steady state at these low wave forcing amplitudes. Note also that the case discussed in section 3 above (γ = 1, h₀ = 600) appears in the central row. Peaks in E(ω) represent variability dominated by a particular frequency. The height of the peak indicates the strength of the variability, in the sense of the magnitude of the changes to u_jet, for example, from strong westerly to strong easterly in the case described above. The width of the peak indicates the coherence or regularity of the variability: a very narrow peak corresponds to almost periodic variability at the frequency of the peak; a broader peak corresponds to a greater distribution of dominant frequencies and less coherent variability.

Let us start by concentrating on the middle row, which corresponds to γ = 1. We see, first of all, that the case discussed above does not represent the most coherent variability obtained in this framework. At γ = 1, the variability obtained with h₀ = 800 is both more regular and more vigorous, as the time–height plot of u verifies; see Fig. 9b. Note also the increase in frequency with increasing h₀, particularly obvious between h₀ = 600 and h₀ = 800, when there is a single peak (cf. also Fig. 2b and Fig. 9b). At h₀ = 400 (still γ = 1) while the dominant frequency is clearly lower, there are now two strong peaks, indicating intermittency in the flow. In this case, approximately every second sudden warming is reduced in intensity: the spectral peak at T = 80 days corresponds to the type of variability discussed above, consisting of recurring sudden warmings and vortex regrowth, while the larger peak at T ∼ 160 days arises from the lower frequency modulation. At h₀ = 300 (not shown) there is a spectral peak at T = 100 days but also significant power at low frequencies, which corresponds to intermittency as illustrated in the time–height plot shown in Fig. 9a. Here a quiescent regime between t = 400 and t = 650 (and, to a lesser extent, between t = 150 and t = 350) interrupts the otherwise regular sequence of sudden warmings. Despite these variations in dynamical behavior, however, we emphasize that coherent internal variability is present over a very wide range of forcing amplitudes, from h₀ = 300 − h₀ = 800.

At still larger forcing amplitudes, for h₀ = 1000, the coherent variability is strongly attenuated. Inspection of time–height plots of u (not shown) indicate that the reason for this is that the vortex remains persistently weak in the middle stratosphere, and does not recover sufficiently to allow the burst of upward wave propagation and wave-induced deceleration of the upper vortex observed at lower forcing amplitudes. A manifestation of this can be seen in the time averaged u in the meridional plane, shown in Fig. 10c, which reveals predominantly easterly winds in the middle stratosphere. A systematic reduction in the time averaged winds with increasing h₀ can be seen from the other panels (Figs. 10a,b), although part of the reason the winds are so much stronger for h₀ = 300 is simply the prolonged period of vortex stability seen in Fig. 9a.

A similar pattern is found for other values of γ, shown in the top and bottom rows of Fig. 8. In particular, note that the dominant period of the response appears independent of γ, decreasing steadily with h₀ in each case from around T = 100 days below h₀ = 400 to around T = 45 days for h₀ = 800. All γ cases also exhibit intermittence at h₀ = 400 (and below), with a broad or double peak structure in E(ω). Note that the apparent loss of coherence for γ = 1.5 and h₀ = 800 is again due to intermittence, evidence for which can be seen in the broader and double peak structure, as well as the power at very low frequencies. In fact, plots of u(z, t) for γ = 1.5 and h₀ = 1000 are similar to Fig. 9a, but squashed in time corresponding to the high frequency dynamics.

Finally, we summarize all simulations considered here by considering the following basic quantities, again based on the time series u_jet(t) at ϕ = 60° and z = 40 km: the time mean, the standard deviation, the average extrema or envelope,² and the dominant frequency. We present these in Fig. 11 as a function of wave forcing amplitude, (with simulations at h₀ = 0, 100, 200, . . . , 1000) for each of γ = 0.5, 1.0, and 1.5.

First consider the dependence of the diagnostics on h₀. There is a certain similarity across all values of γ. All γ cases have steady states for h₀ = 200 and there is a sudden transition to vacillating states at h₀ = 300. Further, all γ show a systematic decrease in time mean u with increasing h₀, although the standard deviation and envelope remain roughly constant around mean. This means that the intensity of variability is roughly independent of h₀, while the average vortex strength decreases with h₀. Finally, all γ cases show a systematic increase in dominant frequency with increasing h₀ (the period decreases from about 100 days at h₀ = 300 to about 40 days at h₀ = 1000); that the maximum frequency is not exactly monotonic appears to be a consequence of the fact that it is not always well defined (e.g., at γ = 1.5, h₀ ≥ 900, or γ = 1, h₀ = 400).

Concerning the dependence on γ, we note the following. First, there is a general increase in time mean u with increasing γ, a simple consequence of the colder equilibrium temperature profile. Second, there is an increase in the envelope or extreme departures from the mean with increasing γ, particularly the large negative values seen for γ = 1.5 at higher forcing amplitudes. There is also an increase in the standard deviation. Both these quantities indicate an increase in vigor of variability as the radiative cooling becomes stronger. Finally, we note that the weaker standard deviation for γ = 0.5 is likely related to the fact that the variability is less regular in these cases: this is not readily apparent from Fig. 8 but visible in height–time plots of u (not shown). In these cases, the vortex is generally weaker, and there appears less scope for sudden changes in the wave transmission properties of the flow of the form described above.

5. Sensitivity and robustness

Since this paper is concerned with the robustness of internal variability, we must now investigate the effect of the other main physical and numerical parameters on the variability. After the parameters γ and h₀, the radiative relaxation rate, τ, and the forcing wavenumber, m, are arguably the next most influential physical parameters, and we examine the dependence on these next. Again we are interested in whether the variability persists over a broad range of these forcing parameter values, and, if it does, which parameters determine qualities of the variability such as the peak frequency and magnitude of sudden warmings. Similarly we wish to establish that the results presented here are not sensitive to the numerical procedure chosen to carry out the computation of the physical model, including the particular resolution used and any unphysical parameterizations such as the sponge layer; this is considered in section 5b.

a. Physical parameters

1) Radiative relaxation rate

We first examine the sensitivity to the radiative relaxation time scale, τ. To do so, we consider three variations of the foregoing simulations, identical in all respects apart from the value of τ, which is (i) halved to τ = 5 days, (ii) doubled to τ = 20 days, and (iii) set to a height-dependent value τ = τ_H(z) representative of the height dependence of radiative time scales in the stratosphere (Holton 1976), where τ_H(z) is given by

View Expanded

with z in kilometers; τ_H varies from around 20 days below 35 km to around 4 days above 35 km.

The results for the series of simulations with γ = 1, h₀ = 0, 100, 200, . . . , 1000 are summarized in Fig. 12. We see, first of all, that in terms of time mean u and standard deviation the general variability obtained in all cases is generally similar to that obtained above (cf. Fig. 11b). In particular, strong variability is once again found over broad range of h₀.

Comparing more closely with the τ = 10 case discussed above, we find that the dominant frequency of the variability is higher for τ = 5 and lower for τ = 20; that is, the frequency of the variability decreases with increasing radiative time scale. In addition, the dominant frequency for the case τ = τ_H is closer to that for the case τ = 20 than for τ = 5. Since τ_H is approximately 20 days in the lower stratosphere and 4 days in the upper stratosphere, this suggests that it is the lower stratospheric radiative time scale that is important in determining the frequency of the variability. Conversely, the time mean u for τ = τ_H is closer to that for τ = 5; hence we may conclude, perhaps unsurprisingly, that it is the upper stratospheric radiative time scale that is important in determining the upper-stratospheric jet strength.

Finally, we find intermittency in certain forcing ranges, which for τ = 20 and τ = τ_H is reflected in the nonmonotonic dependence of the dominant frequency on the wave forcing amplitude. In particular, the dominant frequency as determined by the spectral analysis is not always the frequency of a sudden warming–vortex growth cycle, but sometimes represents lower frequency modulations. Having said this, however, we emphasize that even in these cases coherent variability is still evident. The main point of this exercise is to demonstrate that the variability is robust and persists across a wide range of radiative time scales.

Since the radiative profile τ_H(z) is closer to the actual vertical profile of radiative time scales in the stratosphere, we show as a final example, the variability of the zonal mean winds for this τ, with γ = 1 and the two cases h₀ = 300, 800; see Fig. 13. Note the remarkable regularity of the variability. Comparing this with Fig. 9 we see that the variability is in fact more regular with τ = τ_H than with the constant τ = 10. Also, we see directly that the frequency of sudden warmings is lower in the case τ = τ_H, especially for h₀ = 800. We may even speculate that it is a consequence of the vertical variation of τ_H that leads to the stronger regularity of the variability: on the one hand, the longer relaxation time scale in the lower stratosphere allows a lower frequency variability; on the other hand the shorter radiative time scale in the upper stratosphere means that the vortex more nearly approaches its equilibrium state, presumably resulting in less variation between warming cycles. In any case, the approximation of constant τ used in section 4 does not appear to be a significant shortcoming, and we can say with confidence that coherent internal variability is also a feature of more realistic relaxation rates.

2) Forcing wavenumber

As was seen in Fig. 2, the character of the variability obtained is roughly similar between forcing wavenumber m = 1 and m = 2. In this section, we wish to establish whether the variability obtained with m = 1 is also robust over a large range of forcing parameters. To this end, we have repeated the main series of simulations presented in section 4, namely with γ = 0.5, 1, 1.5, h₀ = 0, 100, . . . , 1000, with identical parameters as before but forcing instead with zonal wavenumber m = 1. To summarize the results we show in Fig. 14 the mean, standard deviation, envelope of extrema, and dominant frequency as a function of h₀ for each of the three values of γ.

The similarities with the m = 2 series in all quantities are clear. As with the m = 2 series, as h₀ increases there is a change from a steady to time-dependent vortex around h₀ = 300, with a reduction in the zonal mean jet maximum and an increase in the variance and envelope. Throughout the range of forcing amplitudes between h₀ = 500 and h₀ = 800 the dominant frequency of the variability generally increases with h₀, but otherwise the behavior is roughly independent of h₀. As h₀ is increased still further, the variability as measured by the envelope of extrema and standard deviation begins to weaken as we approach a more continually disturbed vortex state. These features are all in common with the m = 2 case discussed above.

6. Conclusions

The main purpose of this paper was to determine the extent to which coherent modes of internal stratosphere variability exist and can be found over a wide range of forcing parameters. We further wished to characterize the dominant dynamical processes responsible for the variability and determine which physical parameters influenced features of the variability such as its dominant frequency and intensity. In so doing, we have adopted the expedient of using a stratosphere-only model, with a simple lower boundary condition to represent wave forcing from the troposphere. This is of course an artificial construction that does not provide a complete representation of the troposphere–stratosphere coupled system, but which does have the advantage of allowing the dynamics of the stratosphere to be considered to some extent in isolation.

Coherent variability exists over the full range of forcing parameters that gives rise to a realistic polar vortex, including the strength of radiative cooling and amplitude of wave forcing. The ranges over which variability was found more than span physically reasonable values. Further, the existence of coherent variability does not depend on the radiative time scale or wavenumber of the lower boundary forcing, nor on details of the numerical scheme used to solve the physical system. In particular, coherent variability appears to persist as the number of degrees of freedom increases, although current computational resources are insufficient for a complete analysis.

The variability arises because of the competing effects of radiative cooling, which acts to strengthen the vortex, and wave forcing, which acts to destroy it. It is thus characterized by a cycle of gradual intensification of the vortex followed by a rapid reduction in circulation resembling a stratospheric sudden warming. Vertical structure is important for this cycle: when the vortex is destroyed at low levels, further upward wave propagation is prevented and the upper vortex strengthens radiatively. When the lower vortex also eventually recovers upward wave propagation resumes. Importantly, when this happens, the total wave flux into the domain through the lower boundary increases dramatically, despite the fact that the lower boundary geopotential wave amplitude is fixed. Thus, the lower stratosphere appears to act as a valve, preventing or allowing upward wave flux into the domain according to the state of the vortex, and, in our case, independently of the actual lower level geopotential wave amplitudes. Note that a similar dependence of the upward wave flux at the tropopause on the state of the stratosphere was observed in simulations which included a constrained (time-independent) troposphere (Scott and Polvani 2004). The variability in wave flux does not, therefore, appear to be simply an artifact of the choice of lower boundary condition.

Christiansen (1999) described similar vacillations in a GCM in terms of a descending critical level. Although the wave breaking in Fig. 7 begins at upper levels and progresses downward, it is not clear here to what extent the concept of a critical level applies, given the deep vertical structure of the wave. On the other hand, waves do appear to be prevented from propagating upward when the vortex is destroyed at low levels. Moreover, the dependence of the vacillation frequency on the radiative and wave forcing parameters that is obtained here agrees qualitatively with Christiansen’s predictions based on scaling arguments, which gives a vacillation time scale proportional to the radiative time scale and decreasing with increasing wave forcing amplitude. The dependence of time scale on forcing amplitude found here is also consistent with earlier studies using the quasigeostrophic channel model (Holton and Mass 1976; Yoden 1987) and zonally truncated primitive equation model (Scott and Haynes 2000).

The above results have implications for the potential role for stratospheric influence on the troposphere. If internal stratospheric variability did not exist, then the only way in which the stratosphere could be considered as having an influence on the troposphere would be in its capacity as a mediator of tropospheric variability. We have shown here that this is not the case. Further, the fact that the wave flux entering the model domain is controlled entirely by the state of the lower vortex supports the possibility, previously put forward by Hartmann et al. (2000), that the lower stratosphere may also control the amount of wave flux leaving the troposphere, and hence may influence other aspects of the tropospheric circulation. We will return to this question in Part II, where we examine the relative importance of internal variability and variability arising from time-dependence in the wave forcing.