Recent Arctic sea ice loss in fall has been posited to drive midlatitude circulation changes into winter and even spring. Past work has shown that sea ice loss can indeed trigger a weakening of the stratospheric polar vortex, which can lead to delayed surface weather changes. But the mechanisms of such changes and their relevant time scales have remained unclear. This study uses large ensembles of idealized GCM simulations to identify how and over what time scales the atmospheric circulation responds to short-term surface heat flux changes in high latitudes. The ensemble-mean response of the atmospheric circulation is approximately linear in the amplitude of the surface forcing. It is also insensitive to whether the forcing is zonally asymmetric or symmetric, that is, whether stationary waves are generated or not. The circulation response can be decomposed into a rapid thermal response and a slower dynamic adjustment. The adjustment arises through weakening of vertical wave activity fluxes from the troposphere into the stratosphere in response to polar warming, a mechanism that differs from sudden stratospheric warmings yet still results in a weakened stratospheric circulation. The stratospheric response is delayed and persists for about 2 months because the thermal response of the stratosphere is slow compared with that of the troposphere. The delayed stratospheric response feeds back onto the troposphere, but the tropospheric effects are weak compared with natural variability. The general pathway for the delayed response appears to be relatively independent of the atmospheric background state at the time of the anomalous surface forcing.
The Arctic is warming faster than the rest of the globe (Schneider and Held 2001; Screen and Simmonds 2010; Serreze and Barry 2011). Arctic sea ice cover has undergone unprecedented declines (Comiso et al. 2008; Stroeve et al. 2012b), and the subsurface heat stored in the Arctic Ocean has increased in the past decade (Timmermans et al. 2018). The remote influence of these changes is a pressing topic, in large part because of a possible link to recent episodes of extreme weather in the midlatitudes (e.g., Liu et al. 2012; Tang et al. 2013; Petoukhov et al. 2013). One suggestion for how the Arctic affects the midlatitudes is through the stratosphere: sea ice loss is thought to lead to conditions that favor the upward propagation of planetary waves into the stratosphere, which produce a weakening of the stratospheric polar vortex (Garfinkel et al. 2010; Kolstad et al. 2010; Cohen et al. 2014). The polar vortex signal can then be communicated back down to the troposphere to create a negative North Atlantic Oscillation (NAO) response (Baldwin and Dunkerton 2001; Polvani and Waugh 2004), along with its associated surface climate effects (Kolstad et al. 2010; Scaife et al. 2005; Ambaum and Hoskins 2002; Thompson and Wallace 1998).
A number of studies have investigated the existence of this stratospheric pathway linking Arctic sea ice loss to a midlatitude circulation response, but the results are not clear-cut. There is evidence that, in the present climate, fall sea ice loss is associated with a weaker polar vortex in early winter and a negative northern annular mode (NAM) response in late winter to early spring (Jaiser et al. 2013; Feldstein and Lee 2014; Cohen et al. 2014; García-Serrano et al. 2015; King et al. 2016; Yang et al. 2016; Peings and Magnusdottir 2014; Nakamura et al. 2015). Some modeling studies support these findings (Kim et al. 2014; Nakamura et al. 2016), while others report a strengthening of the polar vortex in response to sea ice loss either earlier or later during the cold season (Cai et al. 2012; Scinocca et al. 2009; Screen and Simmonds 2013; Blackport and Screen 2019). These discrepancies could partly be explained by differences in details of how the sea ice cover is reduced in the model experiments, as it is only ice loss in the Atlantic–Barents Sea sector that seems able to weaken the polar vortex (Screen et al. 2018; Sun et al. 2015; Screen 2017; Zhang et al. 2017). Biases in the background state of models could also play a role. For example, both the strength of the stratospheric response and the sign of the NAO response to sea ice loss can depend on the climatological latitude of the jet stream (Smith et al. 2017).
Furthermore, there are open questions about the exact processes connecting sea ice loss to the polar vortex. It is well known that planetary waves originating in the troposphere and propagating upward into the stratosphere can weaken the polar vortex (Matsuno 1971). In fact, breakdowns of the polar vortex are often preceded by short-term poleward heat flux anomalies, or upward wave activity flux anomalies (Polvani and Waugh 2004; Hitchcock et al. 2013a; Cohen and Jones 2011; Garfinkel et al. 2010). Sea ice loss in the Atlantic–Barents Sea sector is thought to enhance high-latitude upward wave activity fluxes by producing tropospheric circulation anomalies that facilitate upward propagation of waves (Jaiser et al. 2013; Kim et al. 2014; Feldstein and Lee 2014; Sun et al. 2015). Some models reproduce aspects of this stratospheric pathway in response to sea ice loss (Nakamura et al. 2016; Hoshi et al. 2017; Zhang et al. 2017). However, other models find that sea ice loss leads to weaker upward fluxes because the weakened equator-to-pole temperature gradient weakens the baroclinic eddy source (Seierstad and Bader 2009), which is a consistent feature across coupled models (Screen et al. 2018; Smith et al. 2017).
Given these seemingly conflicting lines of evidence, it is legitimate to ask whether fall sea ice loss, or associated surface heat fluxes, can in general modulate the large-scale atmospheric circulation to activate the stratospheric pathway in a consistent and detectable way. And if such a stratospheric bridging mechanism is in play, what is its timing and signature? Why is it that models do not agree on the sign of the response? In reality and in realistic experimental setups, it is difficult to isolate the influence of sea ice on the stratosphere from the influence of other factors (Overland et al. 2016), such as forcing from the tropical Pacific (Feldstein and Lee 2014); substantial noise from internal variability (Sun et al. 2015; Blackport and Kushner 2017); long stratospheric persistence (Shepherd 2000); changes in snow cover and sea surface temperature (Cohen et al. 2014); and interference from winter sea ice loss, which has been shown to dominate the midwinter circulation response in some models (Sun et al. 2015). It is also difficult to separate causes and effects, as, for example, transient sea ice loss may be caused by low-frequency atmospheric circulation anomalies, which then are difficult to disentangle from the atmospheric response to sea ice loss (Blackport et al. 2019). These complications greatly reduce the signal-to-noise ratio of the sea ice influence on the atmosphere and hamper analyses and interpretations (Overland et al. 2016).
Idealized process studies offer a way to clarify the atmospheric response to sea ice loss. The stratosphere responds to Arctic forcing in simplified model setups (Wu and Smith 2016; Zhang et al. 2017), but there are outstanding questions regarding the time scale and dynamics of the adjustment in the context of large internal variability and natural persistence in the stratosphere. Here, we use idealized general circulation model (GCM) simulations to study the atmospheric response to short-term Arctic perturbations designed to represent the transient forcing due to fall sea ice loss. A large ensemble of simulations allows us to estimate mean time scales and amplitudes of the response, and to identify robust mechanisms.
We use the idealized GCM described in O’Gorman and Schneider (2008) and Frierson et al. (2006). The model solves the primitive equations at a spectral resolution of T85 (i.e., with a transformed grid spacing of 1.4° at the equator), with 30 unevenly spaced σ-coordinate levels and the highest level at σ = 4.2 × 10−5. The GCM uses a two-stream gray radiative transfer scheme, in which the radiative flux is only a function of temperature (i.e., there are no cloud radiative or water vapor feedbacks). The radiative forcing at the top of the atmosphere is an approximation of annual-mean insolation, so that variations of the zonal-mean jet position only arise by internally generated variability. The longwave optical thickness τ of the atmosphere has components that depend linearly and quadratically on pressure, and it is proportional to sin(ϕ)2, where ϕ is latitude, with an optical thickness of 7.2 at the equator and 3 at the pole (Frierson et al. 2006; O’Gorman and Schneider 2008). This produces a tropospheric jet in each hemisphere similar to that seen in Northern Hemisphere winter (Fig. 1), despite the annual-mean insolation. The lack of diabatic heating from shortwave absorption by ozone leads to a lack of separation of the tropospheric jet and the stratospheric polar vortex (Fig. 1), but the tropospheric mean state is in good agreement with the observed winter zonal mean. The annular mode is the dominant mode of variability in the stratosphere and troposphere of the model, and its variability is similar to that seen in more complex models (Figs. S1–S3 in the online supplemental material). A better representation of the stratospheric mean state could have been obtained by adjusting the stratospheric radiative forcing (Kushner and Polvani 2004), which may change the dominant adjustment mechanisms to episodic surface forcing.
Earth’s surface is represented by a 1-m-thick slab ocean. This shallow representation of a mixed layer results in a shorter decay time scale of surface energy budget anomalies compared to observations. Turbulent surface fluxes at the lower atmospheric boundary are represented by the standard bulk aerodynamic formulas and Monin–Obukhov similarity theory. A tropical ocean energy flux divergence (Q flux) with an amplitude of 60 W m−2 at the equator is added to mimic the equatorial ocean circulation, which helps to obtain a more realistic Hadley cell and tropospheric jet position (Bordoni and Schneider 2008; Levine and Schneider 2011).
a. Representation of Arctic sea ice anomalies and short-term ocean heat loss
Typically, model studies either prescribe sea ice anomalies [see the overview in Gao et al. (2015)], which artificially fixes part of the surface energy balance, or alter the radiative fluxes to remove sea ice (Screen et al. 2018), which still tampers with the surface energy balance but at least allows for an energetically consistent adjustment of the climate system to the perturbation. A more direct approach is to perturb the surface energy balance itself by imposing Q fluxes that mimic the effect of sea ice variability on the surface heat flux (Fig. 2). There are still some complications here because a substantial portion of surface heat flux variability is not due to sea ice changes (Walsh and Johnson 1979; Fang and Wallace 1994; Deser et al. 2010; Woods and Caballero 2016; Sorokina et al. 2016), but rather drives sea ice changes (Smedsrud et al. 2013; Lee et al. 2017; Timmermans et al. 2018; Blackport et al. 2019). But the approach does allow for more control in the experimental setup, since it circumvents discussions about local ice–atmosphere feedbacks.
Here we prescribe an additional Q flux to mimic temporary sea ice anomalies in the Arctic and the associated perturbation of the surface energy budget. We use a Q-flux convergence/divergence imposed for 20 days (between day −24 and day −5) in an annulus between 70° and 85°N for each ensemble member . The plus sign superscript denotes a heat flux convergence perturbation; the minus sign superscript denotes a heat flux divergence perturbation; the subscript i indexes ensemble members. We use the full annulus for zonally symmetric forcing experiments, while a third of it is used for zonally asymmetric forcing experiments. A third of the annulus is comparable in size to the Barents–Kara Sea area (King et al. 2016); however, the position in longitude (here 0°–120°E) is arbitrary in our setup because there are no other zonally localized features.
The polar boundary layer (BL) stratification affects the upward propagation of perturbations at the surface. A weakening of the BL stratification is observed as part of the long-term trend in the Arctic (Screen and Simmonds 2010). This is expected to continue in the future (Deser et al. 2010). In our GCM, the near-surface static stability in the polar regions is larger than in midlatitudes (Schneider and O’Gorman 2008), but it is still weaker than typical Arctic winter conditions (Dee et al. 2011).
We choose the amplitude of the imposed polar Q fluxes based on an order-of-magnitude analysis of observed surface energy flux anomalies in the Barents–Kara Sea in fall (Fig. 2). The precise value turns out to be relatively unimportant because, as we will discuss, the atmospheric response to surface energy budget perturbations is approximately linear in the perturbation amplitude.
b. Smoothing and output frequency
The model output is sampled as 5-day means and is smoothed between adjacent 5-day blocks. Data at a given time are the centered mean over 10 days. The origin (0) of the time axis is 5 days after the imposed Q flux ends, that is, when no contribution of the Q flux is contained in the centered 10-day mean. The 10-day means centered on day −5 and day −25 contain contributions from the imposed Q fluxes at the beginning and end of the forcing period.
c. Ensemble runs
To isolate the transient response to the polar surface perturbation, a set of N = 50 initial conditions was created from an unperturbed climatological run by saving restart files every 60 days after 1440 days of spinup. From each of the N initial conditions, a pair of ensemble members was started, one (, i = 1, …, N) with a positive imposed Q-flux convergence (warming), and one with a negative imposed Q-flux convergence (cooling). We used Q fluxes on the whole annulus, with amplitude , and Q fluxes on a third of an annulus, with amplitude , to produce three sets of 50 paired ensemble members. The differences of the means , with and , at any time after initialization is the linear ensemble-mean signal. (The linearity assumption is verified in section 3.) The mean of the differences approximates the circulation response irrespective of the background internal variability, because any linear contribution of internal variability vanishes by taking differences within each ensemble pair. Additionally, a nonperturbed ensemble member Ci was started from each initial condition. They were used for one-sided analyses of only positive or negative perturbation responses and with .
Low-frequency variability in the GCM simulation introduces correlations among ensemble members, which increase with decreasing separation of the initial conditions. These correlations reduce the effective ensemble size for estimating a robust response. The decorrelation time scale of most dynamic variables increases with height, so low-frequency variability affects especially variables in the stratosphere. We estimate the ensemble size Nmin needed to distinguish the circulation response to forcing from the background of low-frequency variability at the 5% significance level as follows: First, we take the ensemble mean and standard deviation s of characteristic grid points and variables in a 60-member control ensemble (Ck, k = 1, …, 60). Second, we take a 60-member ensemble of forced runs , with a Q flux of 200 W m−2 in a third of the annulus and with initial conditions spaced 30 days apart, to estimate the amplitude of the forced response . From these, we estimate a minimum ensemble size Nmin needed to distinguish the signal from the control mean at a significance level of 5% assuming t statistics (Screen et al. 2013):
Here, tc is the critical value of the t distribution with the relevant number of degrees of freedom, and s is the pooled standard deviation of and . The results for surface pressure, geopotential near the tropopause, and midlatitude zonal wind are shown in Fig. 3. Significant results in surface pressure (Fig. 3a) and geopotential height at the tropopause level (Fig. 3b) can be expected with about 100 members when their initial conditions are 30 days apart (Figs. 3a,b), except in regions of high wave activity (around 55°N at the surface and 45°N on tropopause level). Figure 3c indicates we need 200 ensemble members to distinguish lower-stratospheric responses from background noise. A similar analysis with ensemble members started at 60-day intervals yields a minimum ensemble size of 100 members. In what follows, we hence use an ensemble of 100 members (50 paired members) with initial conditions spaced 60 days apart.
To test whether the responses or are significantly different from zero, we estimate the PDF of or at each grid point and for each variable by bootstrapping the ensemble, that is, by drawing indices i with replacement to create bootstrap ensembles and estimating for each bootstrap ensemble. The bootstrap PDF was obtained from 1000 bootstrap replications. With this PDF, we test whether the response is significantly different from zero by verifying that the 95% confidence interval does not include zero.
The climatology and standard deviation at each grid point are calculated from a 3600-day-long unperturbed control integration, started from the initial condition of the first member (after 1440 days of spinup). Figure 1 shows the climatological-mean zonal wind and its standard deviation in the latitude band of the lower-stratospheric polar vortex (overbars denote zonal means). As expected, the standard deviation is largest near the core of the jet.
3. Linearity and symmetry
We first analyzed to what extent the response is linear in the forcing by comparing simulations with positive forcing (Q-flux convergence) and negative forcing (Q-flux divergence) of different amplitudes. Figure 4 shows the zonal-mean responses and in north–south geopotential height contrasts in the lower stratosphere for the various sets of simulations (geopotential height averaged from 65° to 90°N minus that averaged from 40° to 55°N). This geopotential height contrast is a measure of the NAM (Thompson and Wallace 1998; Ambaum et al. 2001). The responses for positive (red) and negative (blue) forcings are mirror images of each other (Fig. 4a). Halving the amplitude of the forcing (orange) halves the amplitude of the response (Fig. 4a; see also Hitchcock et al. 2013a). The amplitude of the response is generally, within statistical uncertainties, linear in the amplitude of the forcing (Fig. 4b). The peak amplitude of the response is also approximately linear in the total net energy input to the atmosphere by the forcing; that is, when the length of time over which the forcing is applied is doubled while its amplitude is halved, the NAM response is similar, albeit delayed (see appendix A).
The signal-to-noise ratio of the responses and decreases with time after the forcing is applied (red and blue shaded envelopes). Given the approximate linearity of the response, the signal-to-noise ratio can be improved by using the symmetrized response (black line). As can be seen in Fig. 4a, the symmetrized response remains distinguishable from the background noise of low-frequency variability out to day 90. Therefore, we focus on the symmetrized response in what follows, and unless otherwise noted, we describe the symmetrized response to a positive forcing consistent with reduced sea ice.
We tested to what degree zonal asymmetries in forcing shape the response by comparing simulations with the same zonal-mean forcing (66.6 W m−2), but concentrating it to one-third of a latitude circle, that is, using a Q-flux convergence of 200 W m−2 in a 120°-wide longitude sector. The dashed red and blue lines in Fig. 4a show the NAM response in these experiments with positive and negative forcing, and the dashed black line shows the symmetrized response. Within statistical uncertainties, these responses are indistinguishable from the response to the zonally symmetric forcing (Figs. 4a,b). We found similar results for other variables, such as lower-stratospheric zonal wind. Stationary waves, which have been posited to be important on the basis of other modeling studies (Smith et al. 2010, 2011; Hitchcock et al. 2013b) are absent in our model setup. Thus, our results suggest that they are not essential for the response to high-latitude warming in our simulations. However, they may still play a role in the real atmosphere, where stationary waves are prominent in the Northern Hemisphere.
4. Thermal and eddy-driven response
a. Direct thermal response to zonally symmetric heating
We focus on the zonally symmetric forcing in polar regions, with an imposed Q flux of 66.6 W m−2. The polar heating applied between day −24 and day −5 directly affects the surface and troposphere temperatures in polar regions (Figs. 5a,c). The circulation responds to the forcing with a persistent weakening of the zonal-mean zonal wind (Fig. 5b). Initially, this weakening can be understood from the surface temperature response and thermal-wind balance,
where f is the Coriolis parameter, R the gas constant for dry air, H the scale height, ps the reference pressure, and κ = R/cp the adiabatic exponent with specific heat at constant pressure cp. The overbar denotes the zonal mean. Heating of the surface in polar regions weakens meridional temperature gradients and the vertical wind shear . But the polar surface and near-surface warming decays after around day 30 (Figs. 5a,c), and, after that initial period, meridional temperature gradients and the tropospheric zonal-mean wind relax back to the original state they occupied before the polar heating was applied. Yet the zonal wind in the stratosphere continues to weaken (Fig. 5b). The weakening of the stratospheric wind is largest in absolute amplitude in the upper stratosphere (above σ = 0.01). But relative to the internal wind variations, it has largest amplitude in the lower stratosphere (Fig. 6). The zonal-wind response in the stratosphere peaks between day 0 to day 30 but remains statistically significant in the lower stratosphere over nearly 3 months. This stratospheric response is not simply a result of thermal-wind adjustments to surface temperature gradients, but it is closely related to them.
b. Eddy-driven stratospheric response
While the initial stratospheric response is a result of thermal-wind adjustments to changes in surface temperature gradients, the persistence of the delayed response needs to be understood in the context of eddy–mean flow interactions.
Figure 7 shows the time and height dependencies of the zonal-mean (Eulerian) meridional velocity , zonal-mean vertical velocity (left and center columns), and the zonal-mean zonal wind and potential temperature θ (right column). Color shading indicates the forced responses, and light gray contours indicate the climatology. The climatological stratospheric circulation exhibits northward flow in the subtropics and southward flow in midlatitudes (Fig. 7a). By mass continuity, this is closed by downward flow between 30° and 50°N and upward flow poleward of 50°N (Fig. 7b). The climatological transformed Eulerian mean, or residual, circulation, which takes mass transport by eddies into account, exhibits ascent in low latitudes and descent in middle and higher latitudes, more closely resembling Earth’s Brewer–Dobson circulation (see Fig. B1 in appendix B).
While the residual circulation captures eddy mass transport and hence gives a better description of tracer transport than the Eulerian mean circulation, here we focus on the Eulerian mean because it dominates the response of the circulation to the forcing (cf. the residual circulation response in appendix B). The Eulerian mean stratospheric circulation responds to the surface forcing primarily by weakening, especially after cessation of the surface forcing (e.g., days 31–60; see left two columns of Fig. 7). That is, there is anomalous southward flow in the subtropics and anomalous northward flow in midlatitudes for about 2 months after the surface forcing has ended. By continuity, this is compensated by anomalous upward flow between around 30° and 50°N and downward flow poleward of 50°N.
The weakening of the stratospheric overturning circulation can be understood as a straightforward consequence of the reduced meridional temperature gradient in the troposphere, which weakens meridional eddy heat fluxes and thus weakens vertical wave activity propagation (Edmon et al. 1980). The meridional eddy heat flux is shown in Fig. 8a. It forms the vertical component of the wave activity or Eliassen–Palm (EP) flux, defined, here for simplicity in Cartesian coordinates, as
where primes denote deviations from the zonal mean, . At higher altitudes (above σ ≈ 0.1), the vertical wave activity propagation eventually converts to horizontal propagation (Fig. 8b) and irreversible mixing, which manifests itself as meridional momentum flux convergence (Fig. 8c). At low Rossby number (appropriate in the off-equatorial stratosphere), the dominant balance in the zonal momentum equation is between the Coriolis acceleration of the mean meridional flow and divergence of the meridional eddy momentum flux,
The eddy momentum flux divergence ultimately is a consequence of vertical wave activity propagation, and it weakens because the vertical wave activity propagation weakens (Figs. 8a,c). This causes the weakening of the Eulerian mean meridional circulation (Figs. 7 and 8; see also Fig. B2 in appendix B).
The anomalous subsidence in middle and high latitudes (Fig. 7e) leads to dynamic warming,
because the static stability is positive. What results is dynamic warming of the polar stratosphere, which weakens the zonal mean wind via the thermal wind relation (2) (Fig. 7, right panels). This dynamic warming mechanism is driven by the anomalous Eulerian mean overturning, which warms the polar stratosphere over the course of about a month (green arrow in Fig. 7f). The warm anomaly and associated zonal wind anomaly then decay on stratospheric radiation time scales (about 2 months).
Figure 9 confirms that the polar stratospheric dynamic warming is primarily a result of the dynamic warming associated with anomalous subsidence; dynamic warming associated with eddy fluxes and meridional advection plays a lesser role (see appendix B). The static stability changes only weakly, so the anomalous dynamic warming arises predominantly from the anomalous subsidence (Fig. 9b). Other changes in the dynamic terms entering the thermodynamic budget are small (see dashed line in Fig. 9a and appendix B).
In summary, the delayed response of the zonal wind in the stratosphere arises from eddy-driven changes in the stratospheric circulations. They, in turn, flow from the changes in the tropospheric circulation associated with reduced meridional temperature gradients and reduced vertical wave activity propagation. The crucial reason why the stratospheric response is delayed relative to that of the troposphere is that thermal damping time scales in the stratosphere are longer than in the troposphere.
5. Decay and dependence on the mean flow
a. Decay of stratospheric anomalies and effect on troposphere
Stratospheric zonal wind anomalies respond on two time scales. There is a fast adjustment to surface temperature anomalies through thermal-wind balance, (2), and a slower eddy-mediated response involving an overturning circulation and dynamic heating and cooling, (5). In the ensemble mean, the fast response is evident in the surface warming signal, which decays within a month (Fig. 5c and red contours in Figs. 7c and 7f). The eddy-mediated lower-stratospheric warming lasts at least 50 days, along with the weakening of the stratospheric zonal winds, and these signals decay gradually over radiative time scales (Fig. 7i). The amplitude of the lower-stratospheric warming depends mainly on the integrated amount of heat input from the surface and less on the heating rate (see appendix A).
Stratosphere-to-troposphere linkages beyond 60 days, as seen in Baldwin and Dunkerton (2001), are not clearly recognizable in the ensemble mean in our simulations (Figs. 5b,c). This is because the phases of the downward-propagating signals from the stratosphere have a random component and are weak compared to sudden stratospheric warmings, making it difficult to identify coherent tropospheric signals against the background of internal variability.
b. Dependence on mean flow
To understand the ensemble spread of the circulation response to heating, we investigate how the response depends on the mean flow of each ensemble member.
Figure 10 shows lagged regressions of the potential temperature response in the polar cap (75°–90°N) against the fast stratospheric zonal wind response (days 1–30), both normalized by their climatological standard deviations at each level. Positive values of the regression coefficient indicate a polar cap warming with a weakening of the stratospheric winds, that is, positive potential temperature anomalies associated with negative zonal wind anomalies. We pool ensemble members pairwise to focus on the symmetric response and carry out robust regressions (Sen 1968; Theil 1992) for the ensemble pairs Si. Significance of the regression coefficients is established by estimating a PDF of each local regression coefficient through bootstrapping (section 2c) and establishing significance at the 5% level if zero lies outside the 2.5th and 97.5th percentiles.
The lagged regression analysis shows that more negative anomalies of stratospheric zonal wind lead longer-lasting and more downward-propagating positive polar temperature anomalies (i.e., positive regression coefficients). The regression coefficients reach maximum amplitude in the troposphere and near the surface about 30 days after the early lower-stratospheric zonal wind response (Fig. 10). The early wind response itself depends on the combined effect of the forcing and its (nonlinear) interaction with the background stratospheric zonal wind in any given ensemble member. Members with a weaker stratospheric vortex thus have an increased likelihood of producing significant tropospheric responses (Baldwin and Dunkerton 2001; Hitchcock et al. 2013b). However, the nonlinear interactions also affect whether or not signals propagate downward from the stratosphere, and these are subject to variability generated within individual members.
To examine the sensitivity of the delayed zonal-wind response on the mean-flow sampled by internal variability, we regress the mean zonal wind in the control simulation Ci against the late (negative) zonal wind response Si in the ensemble pair branched off from that control simulation. As a measure of the late response, we use the zonal-wind response in the lower stratosphere averaged over the space–time box spanning 50°–75°N and days 50–75, again normalized by the climatological standard deviation at each level. To establish whether the resulting regression coefficients are significant, we compute a bootstrap PDF of a corresponding regression within the control runs alone, regressing the mean zonal wind in the control simulations Ci against the control zonal wind in the lower-stratospheric space–time box. Significance is then established as in section 5a.
Figure 11 shows the resulting regression coefficients and their significance. The strength of the delayed wind response (green box) is related to the strength of the midlatitude flow during the forcing (days −25 to −5): members that respond to polar surface warming with more weakening of the stratospheric wind between days 50 and 75 tend to be initialized from control states with stronger tropospheric winds. The reason appears to be regression toward the mean, resulting from the model’s natural relaxation toward the mean state. The tropospheric midlatitude wind is related to the meridional temperature gradient in the lower troposphere. Anomalously strong meridional winds are associated with anomalously strong temperature gradients, or, assuming fixed subtropical temperatures, an anomalously cold polar troposphere. The tropospheric temperature at the pole and the meridional gradient naturally relax to the mean state. The natural relaxation toward the mean adds to the response to the forcing.
6. Discussion and summary
We have used a large ensemble of idealized aquaplanet GCM simulations to determine the circulation response to short-term polar surface heating or cooling. We extracted robust and statistically significant aspects of the circulation response from the large background of internal variability using a large ensemble of simulations. The idealized modeling approach circumvents the complexity of the ocean–atmosphere boundary layer, seasonal changes, and stationary waves (Cohen et al. 2014; Overland et al. 2016). Yet many aspects of the circulation response we find in the idealized setting resemble circulation responses found in more complex models (Sun et al. 2015; Kim et al. 2014; Nakamura et al. 2016; Wu and Smith 2016; Zhang et al. 2017).
We chose an idealized GCM setup that produces a tropospheric mean state resembling that observed, while the stratospheric mean state lacks a separation of the zonal-mean jet. The key dynamical mechanisms we identified in this setup can be broken down into a fast thermal-wind response and a slower eddy-mediated response (Fig. 12). Polar surface heating leads to a warmer surface (red stripes on the surface in Fig. 12) and reduces meridional temperature gradients. By thermal-wind balance, the surface heating leads to weakened zonal winds aloft (Fig. 12, top panel). The tropospheric temperature and thermal-wind responses both decay on a radiative time scale (about 30 days; see appendix A) after the forcing ends. But while the tropospheric response persists, the weakened meridional temperature gradients lead to weakened upward wave activity (EP) fluxes in midlatitudes, which eventually lead to a weakened pattern of horizontal wave activity fluxes (momentum fluxes) in the midlatitude stratosphere (Fig. 12, bottom panel). The weakened pattern of eddy momentum flux convergence and divergence is balanced by a weakened Eulerian mean circulation in the stratosphere, with weakened rising motion (anomalous subsidence) in the mid- and high-latitude stratosphere (Fig. 12 lower panel, green arrows). This implies anomalous dynamic warming of the polar stratosphere, which, in turn, modulates the stratospheric zonal flow by thermal-wind balance. The long relaxation times in the stratosphere lead to relatively long (up to 90 days in our GCM) time scales of the stratospheric response, which can in some cases reconnect to the troposphere, when the background conditions set by internal variability are favorable (Fig. 10).
It should be noted that this mechanism differs from that of sudden stratospheric warmings (SSW; or vortex breakdown; Matsuno 1971; Polvani and Waugh 2004; Garfinkel et al. 2010; Hitchcock et al. 2013a). It consequently differs from the usual explanations of Arctic to midlatitude leakage (Cohen et al. 2014; Overland et al. 2016). Our experiment shows a relatively slow (20–30 days) and weak (about 5%–10%) stratospheric warming resulting from a weakening of poleward heat fluxes in midlatitudes and subsequent adjustment of the momentum budget (Hirota and Sato 1969; Kushner and Polvani 2004; Limpasuvan et al. 2005; de la Cámara et al. 2018). In contrast, SSWs are rapid, intense events arising from positive heat flux anomalies by planetary waves that directly perturb the polar stratosphere. Some studies show evidence that polar warming can lead to SSWs (Wu and Smith 2016; Kim et al. 2014; Jaiser et al. 2013), but this pathway is weak compared to the midlatitude pathway in our simulations (Fig. B2b in appendix B; see also Wu and Smith 2016; Seierstad and Bader 2009; Peings and Magnusdottir 2014; Feldstein and Lee 2014; Sun et al. 2015; Nakamura et al. 2015). A better understanding of the conditions under which one pathway is more active than the other could help resolve why models differ in the sign and strength of their stratospheric responses to Arctic sea ice loss (Screen et al. 2018). The outlined mechanism can be also detected in more comprehensive models by analyzing momentum-flux divergences, or residual velocities, rather solely heat fluxes reaching the stratosphere.
A few findings from our simulations are especially noteworthy:
Responses with long time scales can be generated through the long memory of the stratosphere, without memory in surface or near-surface conditions. The cumulative effect of eddies with relatively short time scales on the stratospheric circulation is crucial for generating the late response.
Forcings that are zonally asymmetric and energetically equivalent (in the zonal mean) produce similar responses as forcings that are zonally symmetric in our GCM simulations. That is, stationary waves are not essential for the long-term response. However, stationary waves may still play important roles in the response to high-latitude forcings in Earth’s atmosphere and more complex models (e.g., Smith et al. 2010, 2011; Screen 2017; Smith and Scott 2016).
A robust feature of recent coupled model experiments is a weakening of the poleward flank or shift of the midlatitude jet (Screen et al. 2018; Smith et al. 2017; McCusker et al. 2017; Blackport and Kushner 2016; Deser et al. 2015). We suggest this response is consistent with the thermal adjustment shown in section 4 (Fig. 12, top). Coupled model experiments that show this zonal mean thermal-wind adjustment might also show a common stratospheric response to the reduced surface baroclinicity, as demonstrated here and in idealized models with a more realistic topography (Wu and Smith 2016).
The idealized model responses scale with the amount of heat release and peak around when the surface perturbation terminates, when both the thermal and eddy-driven response are simultaneously significant (section 3 and appendix A). In addition to seasonal changes in the jet position and strength, this may help explain sensitivities of the circulation response to the month and duration of sea ice reduction (King et al. 2016; Blackport et al. 2019).
The stratospheric overturning leads to downward motion in high latitudes, which has to be compensated by southward flow in lower levels, as seen in this model (Fig. 7d). This might be related to cold-air outbreaks in the real atmosphere (Kolstad et al. 2010; Overland et al. 2016).
Model simulations of a future, warmer climate show an ice-free Arctic and increased year-to-year Arctic temperature variability (Wang and Overland 2012; Stroeve et al. 2012a; Holmes et al. 2016; Borodina et al. 2017). In addition, observations in the Arctic Ocean show large amounts of subsurface heat that can be released to the atmosphere (Timmermans et al. 2018). This adds a new perturbation to the circulation that is comparable to the perturbations in this model. The circulation response in our simulations may become more common with an ice-free Arctic.
The idealized GCM simulations were performed on ETH Zurich’s Euler computing cluster. We thank Noel Keenlyside for facilitating this work in an early stage, and Tobias Bischoff, Farid Ait-Chaalal, Robert Wills, and Ori Adam for support with running and analyzing the GCM as well as discussions of the paper. This work was partially supported by the NFR KLIMAFORSK programme (Project 255027).
Change of Perturbation Length
Figure A1 shows the response in four additional ensemble sets in addition to those in Fig. 4. The symmetric response to zonally symmetric forcing is shown for forcings of ±66.6 W m−2 over 20 days (black), ±66.6 W m−2 over 40 days (yellow), ±33.3 W m−2 over 20 days (red), and ±33.3 W m−2 over 40 days (green). That is, the ensembles differ only in the forcing amplitude and the length over which it is applied. The maximum amplitude of the polar surface temperature anomaly in each ensemble closely follows the cumulative forcing integrated over time (Fig. A1a), and the lower-stratospheric NAM anomaly closely follows the surface temperature anomaly (Fig. A1b). This shows that the development of the lower-stratospheric NAM amplitude and its peak are closely related to the integrated surface heating. However, while the surface temperature anomalies always decay on radiative time scales (about 30 days), the stratospheric NAM anomalies decay more slowly (Hitchcock et al. 2013a).
TEM Equations in σ Coordinates on a Sphere
The analyses in section 4b are computed in σ coordinates on a sphere, where σ = p(ϕ, z)/ps(ϕ) is pressure p normalized by surface pressure ps. The transformed Eulerian mean equations in the quasigeostrophic limit in σ coordinates are (Edmon et al. 1980; Andrews and McIntyre 1976, 1978)
where f is the Coriolis parameter; and are the zonal frictional and heating force per unit mass; , with the gas constant for dry air R, reference pressure p0, and adiabatic exponent κ = cp/cυ; y = reϕ is the meridional coordinate with latitude ϕ and Earth’s radius re; and is the surface pressure-weighted zonal mean (with the zonal mean again denoted by the overbar), with primes now denoting deviations therefrom (Held and Schneider 1999).
a. Residual velocities
The residual velocities are defined as
and the residual streamfunction is computed as usual. The vertical component can alternatively be expressed as
and is used in (B3) and (B4). Figure B1 show the residual velocities in the left and center column and the zonal mean wind in the right column. It can be seen that the responses in the residual velocities are mainly driven by the Eulerian component, as discussed in section 4b (see Fig. 7 for direct comparison).
b. Wave activity fluxes
The wave activity or EP flux is defined as (Edmon et al. 1980)
Its divergence is
The wave activity flux and its divergence during the first 30 days (days 0–30) after the end of the forcing are shown in Fig. B2. For plotting purposes, the horizontal wave activity flux component is multiplied by 2πre cosϕ/g, and the vertical component is multiplied by (cf. Edmon et al. 1980). The initial weakening of the near-surface temperature gradient weakens the climatological pattern of eddy heat and momentum fluxes. The tropospheric wave activity flux weakens, as indicated by the downward pointing arrows and the anomalous wave activity flux divergence in midlatitudes (Fig. B2c). As described in section 4b, the warming of the polar stratosphere can be understood by the Eulerian flow that is balancing the anomalous wave activity fluxes in the stratosphere.
c. Thermodynamic equation
The temperature tendency on the left-hand side and the dynamic heating terms on the right-hand side are calculated for the stratospheric polar cap (σ = 5 × 10−3 to σ = 0.1 and 75°–90°N). The subsidence warming term in Fig. 9 (blue line) is the first term on the right-hand side.
Supplemental information related to this paper is available at the Journals Online website: https://doi.org/10.1175/JAS-D-19-0133.s1.