Dynamic winter climate response to large tropical volcanic eruptions since 1600

1. Introduction

[2] Studies of the influence of volcanic eruptions on climate have been hampered by the difficulty of determining which eruptions had a global impact and by the limited availability of historical climate data. Owing mainly to the sparse geographical coverage of surface temperature observations, investigations of this topic though the 1980s primarily relied upon analysis of global, hemispheric or zonal means, sometimes annually averaged [Angell and Korshover, 1985; Sear et al., 1987; Bradley, 1988; Mass and Portman, 1989]. The small signals obtained led to a considerable debate as to whether any robust impacts could be discerned at all. More recent studies were able to examine the spatial pattern of the climate response based on composite data sets of late-19th and 20th century meteorological observations. Though coverage was quite sparse for the early part of the record, these analyses indicated that the average winter response to large low latitude volcanic eruptions is a warming over large regions of the Northern Hemisphere (NH) continents (though typically with marginal statistical significance) [Groisman, 1992; Robock and Mao, 1992, 1995; Kelly et al., 1996; Robock, 2000; Jones et al., 2004]. At the same time, sizeable areas experience cool anomalies. Proxy-based studies examining summer temperature anomalies following large eruptions also found large regional differences in the response [Briffa et al., 1994, 1998; Jones et al., 2004]. Thus the global, hemispheric or zonal averages would be expected to show relatively weak responses, and depending upon the areas included, anomalies of either sign could be obtained (accounting for the earlier debates on this topic).

[3] General circulation model (GCM) experiments have successfully reproduced the apparent winter anomaly pattern [Graf et al., 1993, 1994; Mao and Robock, 1998; Kirchner et al., 1999; Shindell et al., 2001; Rozanov et al., 2002; Stenchikov et al., 2002; Collins, 2004; Shindell et al., 2003]. Most of these simulations were done with model versions specifically designed to give a better representation of the stratosphere than in typical climate models. Analyses of the simulations indicate that the anomaly pattern results from the response of atmospheric dynamics to the injection of volcanic aerosols into the stratosphere, which forces a positive shift in the naturally occurring AO/NAM pattern (which in the Atlantic sector is nearly identical to the North Atlantic Oscillation; hereafter we refer to the AO, NAM and NAO collectively as the AO). During the summer, dynamics are less responsive, so that the radiative cooling from volcanic aerosols is the dominant effect [Robock and Mao, 1995; Kirchner et al., 1999; Shindell et al., 2003]. The combination of typically warm wintertime anomalies and cool summertime anomalies leads to an annually averaged response that is much weaker than the wintertime signal over much of the NH land area.

[4] Previous studies of wintertime surface temperature data covered only the past 120 years, including 6 large tropical eruptions [Robock and Mao, 1992, 1995] or 5 large tropical eruptions [Kelly et al., 1996]. Furthermore, those studies relied solely upon thermometer measurements, which become quite sparse towards the early part of the record. It has therefore been difficult to obtain a reliable estimate of the mean winter response to volcanic eruptions over much of the globe, and even more challenging to clearly define the variability of that response. The previous studies [e.g., Robock and Mao, 1992] have suggested that there are large variations between eruptions, an attribute also noted in the longer-term proxy-based studies of past summer temperature responses [Briffa et al., 1998].

[5] We present here an analysis based on eruptions from the 17th century through the present, which provides sufficient statistical power to characterize both the mean and the variability of the climate response to large volcanic eruptions. The analysis is based on newly available cold-season proxy-based reconstructions of surface air temperature anomalies [Rutherford et al., 2004] in addition to instrumental data. The dynamical response to external forcings such as volcanic eruptions is most important during the NH winter, so we focus on the cold-season. Warm-season anomalies are largely radiative.

2. Methodology

[6] We use instrumental surface temperature data from 1856 to 2001. Up until 1980, these data are based upon the Climate Research Unit (CRU) of the University of East Anglia data set [Jones et al., 1999; Jones and Moberg, 2003], including instrumental land air and sea-surface temperatures (SSTs), as described in Rutherford et al. [2003]. From 1980 to 2001, we use the GISS surface temperature data set which contains land temperatures from meteorological observations [Hansen et al., 1999, 2001] and SSTs from the Hadley Centre [Rayner et al., 2003] (as in the CRU data set). To extend our analysis back before the advent of modern measurements with broad spatial coverage, we use reconstructions of large-scale surface temperature patterns which extend back several centuries. These reconstructions are based on networks of diverse proxy information (e.g., tree rings, ice cores, corals, thermometers, and historical records). This data set uses an updated methodology [Rutherford et al., 2003] applied to an earlier data set [Mann et al., 1998], and includes independent cold (October–March) and warm-season mean temperature patterns. These seasonal reconstructions have been successfully cross-validated against other proxy data and modern observations [Rutherford et al., 2004]. The data maintain relatively complete areal coverage through the centuries. The amount of input proxy-data decreases towards the early part of this record, however. Additionally, the dating and identification of volcanic eruptions becomes less reliable further back in time (e.g., the unknown source of the 1809 eruption). We therefore have been extremely conservative in use of this data for years prior to 1800, and include only two well-dated and positively identified eruptions (see below). Cold-season averages are used for all data sets. While analysis of December–February would have facilitated comparison with earlier work, the proxy-data are not available at that resolution. Model results and the recent meteorological observations indicate that the two seasonal periods are broadly comparable, though they suggest that the extended cold-season anomalies are somewhat muted compared with the December–February average. All years refer to the year of the October through December, i.e., the eruption year (note that this convention varies, so care should be taken in comparisons between various studies).

[7] Years containing large volcanic eruptions were selected using a reconstruction of radiative forcing from volcanic eruptions [Crowley, 2000]. To create that reconstruction, ice core aerosol information from multiple sites in Greenland, augmented with Antarctic data, was calibrated against optical depth information for modern eruptions. This provides estimates of the aerosol loading of past globally influential eruptions, which were then converted to equivalent radiative forcing. Beginning with the 1960s, modern observations [Sato et al., 1993] were used instead of ice core derived data. We then used objective criteria to select eruption years from that data set. We chose all years for which the negative radiative forcing in year n + 1 minus that in year n was greater than a threshold value. The results using a threshold of 1 or 3 W/m² are denoted by RF1 and RF3, respectively. Since the aerosol loading decays very rapidly, usually in one or two years, there is no bias against eruptions during periods when volcanic activity was greater than normal. The use of an abrupt jump works extremely well based on a close examination of the ice core derived record. The only time when this gave unusual results was in selecting both 1830 and 1831. As several volcanoes erupted during these years, we used only 1831 when the forcing was several times larger than 1830. Selection based on an abrupt jump in forcing does not successfully isolate eruption years using modern optical measurements, however, as these show more gradual changes, presumably due to their higher temporal resolution. We therefore used surface temperatures from the known years of the large eruptions of Mt. Agung, El Chichón and Mt. Pinatubo, but took the radiative forcing (for selection purposes) as the maximum value following the eruption, which in all three cases was the subsequent year.

[8] We further winnowed the data set by selecting only those years that could be definitively matched to large (volcanic explosivity index (VEI) >=3), known eruptions at tropical latitudes based on a historical catalogue [Simkin and Siebert, 1994]. We note that eruptions with a VEI of 3 often do not inject much material into the stratosphere, as a large VEI is a necessary but not sufficient condition for stratospheric injection and global radiative impact. Our use of estimated radiative forcing should select only globally important eruptions, however. Years with large radiative forcing and either a known large eruption at middle or high latitudes, or no known eruption (implicitly assumed to have been high-latitude) were thus excluded. Being extremely conservative with the early portion of the record, only two years were chosen prior to 1800, as noted earlier. Lastly, we evaluated a set of 12 eruptions since 1883 including high-latitude and smaller eruptions that were examined previously [Robock and Mao, 1992] for comparison, which we call RM92. Table 1 lists the years included in each analysis.

[9] Surface temperature anomalies were calculated by taking the temperature during the cold season following the eruption (i.e., October to December of the eruption year and January to March of the subsequent year) and subtracting the cold season average over twenty neighboring years without large aerosol forcing. We took the nearest ten years in each direction whose negative forcing was less than 0.5 W/m² compared with the average of the preceding 50 years (except near the end of the record, when 10 following years were not available). For example, the background for Tambora (1815 eruption) was calculated by averaging the years 1818–1827 and 1802–1814, without 1809–1811 (following the aforementioned unknown eruption) which had too much volcanic forcing (as did 1816 and 1817). The anomaly was then created by subtracting this background from the 1815 temperatures. The removal of the 50 year average forcing was designed to account for the greater background levels present in the more sensitive modern optical measurements, but in practice only changes a very few years in the 20-year averages. The comparison against 20-year averages centered approximately around the eruption year should ensure that the signal does not contain any long-term trend components.

3. Data Analysis

[10] The mean surface temperature anomalies for each individual eruption in the RF1 and RF3 analyses are shown in Figure 1. It is clear that the variability between eruptions is large, demonstrating the need for statistical analyses of multiple eruptions. We note that using the GISS temperature data set during 1880 to 1980 gives extremely similar results to this instrumental analysis. The mean surface temperature anomalies in the composite RF1 and RF3 analyses (Figure 2) show very similar spatial patterns and magnitudes to one another (note that they include 8 eruptions in common out of a total of 11 and 8, respectively). The so-called ‘winter warming’ pattern evident in Figure 2 is similar to that obtained by other analyses of tropical eruptions [Robock and Mao, 1992, 1995; Kelly et al., 1996] that basically used a subset of the eruptions included here (5 of the 6 eruptions they examined are included here, with the somewhat smaller 1974 eruption of Mt. Fuego being the only one excluded from our analyses). The areas for which the anomalies are statistically significant at the 95% confidence level are substantially larger than those obtained from the earlier analyses, however (all statistical significances are reported at the 95% confidence level base on a student's t-test). Previous studies used both more limited data sets and a shorter three-month seasonal average. Figure 2 shows that in both our analyses, the broad warming over northern Eurasia is significant over much of the region, while the cool anomaly located to the south is significant over northern Africa, the Middle East, and even as far eastward as China. A significant warm anomaly is found over nearly the entire eastern United States (especially in the RF3 analysis), while a significant cool anomaly occurs over northern Canada, the Labrador Sea area, and southern Greenland. In the RF1 analysis, this cool anomaly extends westward all the way to the eastern tip of Siberia, and is significant over much of Alaska. Over the oceans, cool anomalies over the tropical and subtropical western Pacific and the tropical Atlantic are significant. The similarity between the two analyses shows that nearly all the statistically significant areas are robust features even with the addition of 38% more years in the RF1 case.

[11] The surface temperature anomalies at midlatitudes are consistent with increased westerly circulation in the NH associated with a forced positive shift in the mean state of the AO. This favors a more northerly storm-track across the Atlantic, bringing warm, moist air to Northern Europe and Russia, and cool, dry air to the Mediterranean basin (as in the NAO paradigm [Hurrell, 1995]). Colder continental air is carried to the Bering Strait region. Over North America, a combination of increased westerly advection and the resulting reduced meridional flow yields a warmer eastern United States, while northern Canada, the Labrador Sea area and southern Greenland are likewise more isolated from lower latitudes and receive greater outflows of cold air from the continental interior. Indeed the spatial pattern obtained by a regression of surface temperature onto the strength of the AO [Thompson and Wallace, 2000] is nearly identical to the anomaly patterns shown in Figure 2.

[12] Both analyses show statistically significant warm anomalies in the eastern equatorial Pacific consistent with a positive El Niño Southern Oscillation (ENSO) signal (Figure 2). The coincidence of ENSO events is a perennial problem in determining the impacts of volcanic eruptions, particularly those that occurred during recent decades (similarly, an ENSO component could perhaps also be present in the AO pattern). To minimize these effects, we examined an alternate RF1 data set. By happenstance we discovered that an analysis of the volcanic years 1809, 1815, 1831, 1835, 1840, 1843, 1883, 1902, 1912, 1924. 1982 and 1991 showed a very similar pattern to the RF1 pattern presented in Figure 2, but that when 1982 and 1991 were excluded, this showed less than 0.1 C anomaly in the eastern equatorial Pacific. The ENSO signal is thus minimized, though due to nonlinearities, the cancellation of El Niño and La Niña in the tropical eastern Pacific does not necessarily fully remove their influence elsewhere. While these years did not all meet the selection criteria (some eruptions were high-latitude), the comparison nevertheless provides a rough approximation of the response pattern without ENSO. The resulting surface temperature anomaly spatial patterns were very similar to those seen in the RF1 and RF3 analyses over the Eastern Hemisphere, though with a slightly reduced magnitude. The patterns themselves were somewhat different over western North America, however. In the regions where the anomalies are statistically significant in the RF1 and RF3 analyses (the eastern United States, northern Canada and the Labrador Sea area), the analysis without ENSO gave a similar response. Over the western United States and Canada, however, that analysis showed statistically significant large cool anomalies. Over the western United States, Alaska and northwestern Canada, there are hints of similar anomalies in the RF1 and RF3 analysis, but over western Canada anomalies of the opposite sign are seen in the RF1 and RF3 analyses. Since these patterns are not robust, we place less faith in their veracity.

[13] Very similar ENSO effects were found in other observational studies [Robock and Mao, 1995; Mao and Robock, 1998; Yang and Schlesinger, 2001], which also showed the largest impacts on their volcanic analyses in western North America, and little impact on the Eurasian signal. Kelly et al. [1996] argued that removal of the ENSO signal at the grid-point, monthly level was likely to be highly uncertain due to the difficulty in unequivocally defining the ENSO signal at that resolution. We believe that ENSO remains a confounding factor in analyses of the effects of volcanic eruptions on climate over North America outside of the northeast quadrant due to the difficulty of reliably separating these two patterns. Though irregular, large El Niños have occurred close to several of the large tropical eruptions examined here, making the two signals tightly interwoven. It is intriguing to consider whether this is more than coincidental [see, e.g., Hirono, 1988; Robock et al., 1995; Self et al., 1997; Adams et al., 2003].

[14] The second winter following large tropical eruptions is thought to show a similar surface temperature response, though somewhat less AO-like than during the first winter [Robock and Mao, 1995; Kelly et al., 1996]. Figure 3 shows the RF3 analysis for the second winter. Clearly the Eurasian pattern and the cooling of the northern portion of North America and the Labrador Sea area are quite similar to those seen during the first winter. Independent analyses of European temperatures appear to show similar results (E. Fischer, personal communication, 2004). Cooling over the ocean basins is more pronounced, as would be expected given the longer adjustment time to the negative radiative forcing. The persistent El Niño signal has diminished, and most strikingly, the warming over the eastern United States and central Canada has vanished, as in previous analyses [Robock and Mao, 1995; Kelly et al., 1996]. This provides further evidence that anomalies in this region may be more reflective of ENSO than of the volcanic signal (although the canonical AO pattern does show a warm anomaly there).

[15] Restricting the analysis to the 1600–1856 period, during which the proxy reconstructions are entirely independent of the later instrumental data, yields similar anomaly patterns, but with magnitudes of only 25–60% that seen in the full analysis. Anomalies are statistically significant for at least part of the large cool anomalies over the northern part of North America and the NE Africa-Middle East cooling, but not for the warming of northern Eurasia or the eastern United States. The reduced magnitude of the response in the reconstruction is likely a result of the decreased resolution of variance at regional spatial scales prior to the introduction of instrumental data (see variability discussion below). Additionally, two eruptions with noticeably reduced continental warming response are included in that time period (1641 and 1835, see Figure 1). However, the qualitative agreement between the earlier and later centuries demonstrates that the overall response pattern is robust throughout the record. We note also that significant stratospheric ozone loss in response to volcanic eruptions has occurred since the late 1970s due to the presence of anthropogenic halogens. This may have amplified the effects of the most recent eruptions [Robock, 2000; Stenchikov et al., 2002; Shindell et al., 2003], which show larger magnitude responses than any of the earlier eruptions in our analysis.

[16] We have compared our RM92 analysis (Figure 4) with the analyses of Robock and Mao [1992, 1995]. Several differences exist between the analysis techniques: (1) we used October to March averages, while they used December to February, (2) we calculated anomalies relative to 20-year averages over neighboring years, while their anomalies are with respect to the 1886 to 1992 average, and (3) we did not remove the ENSO signal, while they did. Despite these differences, we obtain fairly similar results to theirs, and these results also resemble our RF1 and RF3 analyses. Thus the response pattern appears to be fairly robust to changes in surface temperature baselines, seasonal averaging and the surface temperature data set (which was slightly different in our analysis, especially over the oceans).

[17] Using the RM92 set of volcanoes gives statistically significant anomalies only over a few small areas, however, as in the original RM92 analysis. The total significant area is much less than that seen in our other analyses. These particular eruptions were chosen based on dust veil index and VEI values, which are less directly related to radiative impact, and indeed they show a much lower mean forcing in the reconstruction used here (Table 1). Additionally, they include 6 middle and high-latitude volcanoes along with 6 tropical ones. The surface temperature response during the second year following high-latitude eruptions is somewhat similar to the pattern during the first winter following tropical eruptions, while for midlatitude eruptions a similar response seems to occur in either the first or second year, making the selection rather arbitrary. For both middle and high latitude eruptions, however, the response in some areas, such as Europe, southern Asia and the Labrador Sea area, is quite different from the response to tropical eruptions [Robock and Mao, 1995]. Indeed, the responses to low and high latitude eruptions may arise in different ways [Graf and Timmreck, 2001]. Thus the averaging of these two patterns together, along with the weaker average forcing, gives a lower signal-to-noise ratio than in our RF1 or RF3 analyses. The analysis of Kelly et al. [1996] shows similar surface temperature anomaly patterns as well. Again, however, these are statistically significant over a much smaller area, having been derived from only 4 tropical eruptions (results without Pinatubo are shown) and having used a three-month seasonal average. One other analysis extending back to the 17th century has been published [Lough and Fritts, 1987]. They show only the winter response averaged over the three years following eruptions, however, and only over the United States. Their result indicates a warm anomaly over the central and western part of the United States. Since their analysis covers several subsequent years, and a large fraction of the eruptions they include were middle and high latitude volcanoes, we would not expect the results to agree with ours.

4. GCM Simulations

[18] The response seen in the RF1 analysis, with a mean radiative forcing of −3.77 W/m² can be reasonably compared with results from ensemble model simulations of the eruption of Mt. Pinatubo, with a forcing −3.70 W/m² (all forcing values are from the data of Crowley [2000], for consistency, though the forcing in a given GCM simulation may be different). Simulations by several groups have demonstrated that GCMs are capable of capturing the dynamical response to volcanic eruptions and reproducing the ‘winter warming’ pattern of surface temperature response [Graf et al., 1993, 1994; Mao and Robock, 1998; Kirchner et al., 1999; Shindell et al., 2001; Rozanov et al., 2002; Stenchikov et al., 2002; Collins, 2004; Shindell et al., 2003]. As an example, Figure 5 shows the RF1 pattern compared with the cold season surface temperature response seen in two GISS GCMs. We show results from the state-of-the-art GISS ModelE run at 4 by 5 degree horizontal resolution with 20 vertical layers extending to ∼62 km. In this case, results are taken from an ensemble of 5 simulations, and the average response following the eruptions of Pinatubo, Santa Maria, and Krakatau (mean forcing −3.68 W/m²) is presented (thus 5 × 3 = 15 eruptions are included, giving us good statistics). That model was run with prescribed SSTs, so that the surface air temperature response over the oceans is minimal. We also show results from an ensemble of 10 Pinatubo simulations with the coarse resolution (8 by 10 degrees) older GISS model II GCM in a version with 23 vertical layers extending up to ∼85 km, well into the mesosphere [Shindell et al., 2003]. The atmospheric model was coupled to a mixed-layer ocean in these runs. Both models were driven by time-varying stratospheric aerosol optical properties (spatial and temporal distribution, effective radius, and optical thickness) taken from the GISS data set [Sato et al., 1993; Hansen et al., 1996]. This data set is primarily based upon SAGE II measurements during the early 1990s, with available optical and historical information for earlier periods. The vertical profile of the aerosols covers 15–35 km in 5 km steps and the optical properties are described at several key wavelengths (further information available at http://www.giss.nasa.gov/data/strataer/).

[19] Both models show a positive shift in the mean state of the AO, and are able to capture much of the warming over northern Europe and Russia and the cooling over the Middle East seen in the observations. The magnitude of the anomalies in both regions agrees fairly well with the data analyses. Though the response patterns in the two models are broadly similar, the changeover from cooling to warming is too far north over Eurasia compared to the observed pattern in the older GCM. Its coarse resolution limits its ability to accurately simulate the location of North Atlantic storm tracks. Indeed other low resolution simulations show a similar bias [Kirchner and Graf, 1995]. Most higher- resolution GCMs obtain ensemble-mean spatial patterns that more closely resemble the observations over Eurasia [Kirchner et al., 1999; Stenchikov et al., 2002], as do the results from modelE. Nevertheless, the coarse resolution model results do not differ dramatically from the higher resolution model. The coarse resolution GISS model with a mixed-layer ocean also reproduces the radiatively driven coolings over the equatorial Atlantic and western Pacific (since the model does not simulate ENSO, it is not practical to compare the eastern Pacific).

[20] Over North America, the models capture the cool anomaly beginning in easternmost Russia and extending across North America to the Labrador Sea region. Neither model reproduces the warm anomaly over the eastern half of the United States seen in the historical analyses. Those observational analyses, however, show an overall positive ENSO state, while the mixed-layer GCM shows a persistent La Niña-like condition in these simulations (and in fact the model results look fairly similar to the alternate RF1 analysis without ENSO). Most other GCM studies [Kirchner et al., 1999; Stenchikov et al., 2002; Collins, 2004] also do not show a warm anomaly over the eastern United States, despite the fact they also report a positive shift in the mean state of the AO. Furthermore, this feature is no longer apparent in the second cold-season (Figure 3). Thus the anomaly in this region appears to be dominated by ENSO, with the smaller volcanic contribution difficult to isolate reliably.

[21] GISS model results for the second cold-season (not shown) display the quadrupole pattern of cooling over northeastern North America and Greenland, warming over some parts of the United States and over northern Eurasia, and cooling over the Mediterranean, north Africa and south Asia, as do the simulations of Stenchikov et al. [2002]. In the GISS model II case, these do not project strongly onto the AO for the second cold-season following Pinatubo (+0.1 ± 0.3 mb, where the AO is defined as the opposite of the mean sea-level pressure change poleward of 60°N), though they do somewhat more for simulations with double and triple the Pinatubo aerosol (+1.1 ± 0.6 and +0.5 ± 0.6 mb, respectively) [Shindell et al., 2003] and in another GCM [Stenchikov et al., 2002].

[22] The ensemble-mean response to volcanic eruptions in other GCMs is also an AO-like pattern during the NH winter immediately following the Pinatubo eruption [Kirchner et al., 1999; Rozanov et al., 2002; Stenchikov et al., 2002]. While results for December or December–February are shown, the general pattern is fairly consistent through the entire cold-season used here. All three of these previous ensemble studies were done with models containing at least a fairly well-resolved stratosphere (though less so in the Kirchner et al. simulations than in the other two), as in the GISS model, but with prescribed SSTs (as in the ModelE simulations). Since both prescribed and calculated SSTs reproduce the AO-like climate response, it appears that how the ocean surface conditions are specified is not a critical factor in the generation of an AO-like response, though the SSTs play a significant role in how that response is manifested [Kirchner et al., 1999]. This is as one would expect given the ocean's long timescale for heat adjustment. Since the models were all forced with stratospheric aerosols, it is clear that the stratospheric composition changes are responsible for the altered surface climate. The primary mechanism by which this is thought to take place is: (1) aerosol heating of the sunlit portion of the lower stratosphere enhances the meridional temperature gradient, (2) this strengthens the westerly zonal winds near the tropopause, (3) planetary waves propagating upwards in the troposphere are refracted away from the pole due to the altered wind shears, further allowing the westerlies to strengthen, (4) the enhanced westerlies propagate down to the surface via wave-mean flow interaction reinforced by a positive feedback between the zonal wind anomalies and tropospheric eddies, (5) strengthened westerly flow near the ground creates the surface temperature response pattern typical of the AO. This mechanism is consistent with the behavior of GCMs [Shindell et al., 2001; Stenchikov et al., 2002], mechanistic models [Eichelberger and Holton, 2002; Polvani and Kushner, 2002], and observations [Kodera, 1994; Perlwitz and Graf, 1995; Robinson, 2000; Lorenz and Hartmann, 2003; Perlwitz and Harnik, 2003]. This mechanism is also consistent with the observation of a strong dynamical response only during the NH cold-season, when dynamical connections between the troposphere and stratosphere are strongest. The robust ability of GCMs to capture the surface climate response to the injection of volcanic aerosols into the stratosphere is thus a clear demonstration of stratospheric forcing of surface-level climate dynamics.

[23] Stenchikov et al. [2002] proposed that in addition to the above mechanism, cooling of the tropical troposphere by the overlying aerosols leads to a reduced meridional temperature gradient within the troposphere. This decreases the generation of planetary waves which then allow the zonal winds in the lower stratosphere/upper troposphere region to strengthen and affect surface climate as in the mechanism outlined above. Though this tropospheric mechanism alone reproduced much of the Eurasian climate response in their simulations, and therefore it seems plausible that it contributes to the overall climate response, it yielded a poorer match to the historical analyses over North America than their full stratospheric aerosol simulations. While the additional mechanism does not arise from stratospheric heating, it still takes place via wave-mean flow interactions involving the stratosphere. Additionally, the simulations of Kirchner et al. [1999] showed that the imposed SSTs can modulate the stratosphere-troposphere interaction. Thus there is ample evidence that the climate response to volcanic forcing is not solely dependent upon the stratosphere, but rather is a function of both the local and the tropospheric/surface radiative effects of the stratospheric aerosols and the interaction of those effects with the existing climate state (including SSTs and perhaps the quasi-biennial oscillation). Nevertheless, the dynamic response to large eruptions appears to be among the clearest indications thus far of a significant role for stratospheric perturbations in surface-level climate change.

[24] While three of the four GCMs discussed above included a well-resolved stratosphere, the ensemble simulations of Collins [2004], in contrast, used a model with only 6 levels in the stratosphere, and produced the weakest winter warming response of any of the GCMs discussed here. Additionally, while the simulations of Rozanov et al. [2002] captured the AO-like response to volcanoes very well, prior simulations with the same model failed to do so [Yang and Schlesinger, 2002]. The main difference between the model versions was a modification to the gravity-wave drag scheme which improved the simulation of the NH polar night jet [Rozanov et al., 2002]. Similarly, simulations with different versions of the GISS GCM showed that the representation of the stratosphere was crucial to simulation of anything more than quite weak changes in the AO in response to greenhouse gases [Shindell et al., 1999]. These results provide a further strong indication of the critical role of the stratosphere in the dynamic response to external forcing.

5. Mean Response Versus Internal Variability

[25] The variability within the GISS model ensemble and within the observational analyses is extremely large. By examining a great many responses the standard deviation can be determined, giving the likelihood for a value to be within a particular distance from the mean if one had only a single response pattern. The physicist August Beer pointed out that this is a ‘peculiar bit of hindsight.’ In this case, however, the likelihood of the response to a single eruption lying near the mean is of great interest in interpreting individual past and future eruptions. We therefore now examine the standard deviation for an individual eruption, though the statistical significance of the mean response shown earlier is of course based upon the standard deviation of the mean.

[26] Examining surface temperatures using the RF1 case, the anomaly-to-standard deviation ratio shows that the standard deviation exceeds the signal over nearly all locations (Figure 6). However, even a shift that is within the internal variability can produce a sizeable increase in the probability of a particular anomaly. Picturing a probability distribution function shifted by the mean anomaly at each point (i.e., almost always less than its half width) gives a good idea of how the distribution of probable temperatures can be biased even though the interannual variability is quite large. For example, large regions of the NH continents show an anomaly-to-standard deviation ratio greater than 0.5, meaning those areas have more than a 69% probability of showing a warm anomaly during the winter following a large tropical eruption. Similarly, negative anomaly-to-noise ratios exceeding −0.75 occur over much of NE Africa and the Middle East, indicating more than a 77% chance of a negative anomaly.

[27] In the coarse resolution GCM simulations, the standard deviation exceeds the mean response at all locations (not shown). The mean signal is within 25% of the standard deviation over northernmost Siberia, near the southern tip of Greenland, over part of the Arabian peninsula, and over some tropical ocean regions (the broad pattern resembles that shown in Figure 5). Those areas thus exhibit more than a 77% probability for an anomaly of the same sign as in the mean response pattern. These values correspond reasonably well with the results obtained by Collins [2004]. He analyzed 20-member ensemble simulations with the Hadley Centre GCM of the response to the El Chichón and Pinatubo eruptions and found that over broad areas of northern Eurasia, the chance for warmer than normal temperatures increased from the default 0.5 to 0.7–0.8. Thus the observations and the GCMs demonstrate that the forced dynamical response to volcanic eruptions is weak compared with natural variability, though it can be isolated with a sufficient number of eruptions.

[28] The model's variability in surface temperatures is roughly 40% greater than that in the observations over NH continents. The historical data sets are able to resolve approximately 60% of the cold-season variance at the hemispheric mean level, but this fraction decreases at smaller spatial scales. At the broad regional scales of interest here, the reconstructions retain about one-third of the variance prior to the introduction of instrumental data, a value that changes little through the past several centuries. Thus we would expect that the variability in the reconstructions is a lower limit. As more proxy records become available, it may become possible to improve reconstructions and allow better estimates of historical variability at small spatial scales. At present, the model's variations appear to be consistent with the observational constraints, though the latter are relatively weak. Since the model does not include a dynamic ocean with an ENSO, however, it's variability is also negatively biased. Furthermore, the model has repeatedly simulated a single eruption with identical forcing each time, while the actual eruptions would have had variable forcing.

[29] We note that the observed surface temperature anomalies following the eruption of Pinatubo were in fact fairly similar to the mean response from the historical eruptions. This is quite fortunate, because this single case is often used in analyses, since is was very well-observed, and is implicitly assumed to be representative of eruptions in general. Even for Pinatubo, however, the match with the mean response is poor in some areas. The southern portions of Europe show more cooling than in the historical analyses, while western Canada shows a large warming (probably ENSO related) which is not present in the mean response.

[30] The shift in the AO is a useful indicator of the overall NH extratropical response as it is the dominant variability mode there. In the GISS 10-member ensemble Pinatubo simulations the response was +1.8 ± 1.5 mb. Thus the standard deviation is comparable to the typical response, and there is even a 15% chance of getting a response that is of the opposite sign to the mean. A similar conclusion was reached by Collins [2004], who showed that in their ensemble simulations there was a positive bias in the NAO following large eruptions, but it was small compared to natural variability. These results further demonstrate that analyses based upon many eruptions are required to get an adequate picture of the response even at large spatial scales.

6. Discussion and Conclusions

[31] Through an analysis of the largest tropical eruptions of the past four centuries, we have confirmed that the mean cold-season climate response to volcanic eruptions in the NH is a distinct pattern of surface temperature anomalies consistent with a dynamical shift in the AO, the leading variability mode of the NH. The anomalies are apparent during the first and second cold-seasons following these large eruptions. Using such a long data set has allowed us to obtain sufficient sample size that the results are statistically significant at the 95% confidence level over very large areas of the NH continents. Our results also demonstrate that the interannual variability is larger than the mean response to volcanic eruptions nearly everywhere. Nevertheless, taking into account the typical AO shift following eruptions should yield sizeable gains in predictive capability.

[32] Comparison with simulations of volcanic eruptions by several GCMs shows that climate models are capable of reproducing the main features of the observed mean surface temperature response pattern and that they appear to capture the variability as well. In such models, the volcanic aerosols heat the sunlit portions of the stratosphere, enhancing the strength of the wintertime Arctic vortex through the thermal wind relationship. The strengthened westerly winds in the lowermost stratosphere then propagate down in to the troposphere via interactions with planetary waves, and the enhanced surface westerlies create the typical AO spatial pattern of temperature anomalies. Since stratospheric variability is strongly correlated with tropospheric variability only during the cold-season, when planetary waves are strongest, it is not surprising that these dynamical effects do not appear to be important during the warm-season. Instead, the volcanic aerosols have a more purely radiative impact during summer, causing a general cooling response. Similarly, the effect of eruptions on cold-season temperatures in the Southern Hemisphere (SH) is likely to be much smaller than for the NH, as planetary wave activity is generally much less owing to the arrangement of the continents. This leads to the presence of a strong coupling between the stratosphere and troposphere occurring only briefly during the southern spring and fall as opposed to the entire cold-season in the NH. A small signal, though, could be present in the SH via enhancement of the Southern Annular Mode. Analysis of sparse instrumental records for recent eruptions does not clearly show this pattern [Robock and Mao, 1995], however (and there is not yet enough proxy-data available for reliable reconstructions of the SH). Temperatures following El Chichón seem to display an annular mode-type response, but this is not present following Pinatubo [Collins, 2004]. This suggests that any signal is indeed too weak to appear in the limited data available.

[33] Unlike more controversial cases in which the above mechanism may also be important, such as the impacts of increased greenhouse gases or solar forcing, the role of the stratosphere is much clearer in the case of volcanoes. Unfortunately, the data record cannot provide a direct confirmation of stratospheric influence on tropospheric climate as not enough historical circulation and pressure data are available. However, the impacts of volcanic aerosols on other potentially important factors such as SSTs simply do not have sufficient time to be dominant by the winter immediately following eruptions (though the tropical SST state seems to influence the response). This is evidenced by the ability of models with fixed SSTs to simulate the dynamic response to eruptions quite well. Thus the consistency of the response to volcanic eruptions between models and the observational analyses, and between different GCMs, is a powerful validation of the ability of GCMs to simulate forced dynamical changes and of the role of stratospheric temperature and wind anomalies in affecting surface climate via modulation of the AO.