Global temperature response to the major volcanic eruptions in multiple reanalysis data sets.

. The global temperature responses to the eruptions of Mount Agung in 1963, El Chichón in 1982, and Mount Pinatubo in 1991 are investigated using nine currently available reanalysis data sets (JRA-55, MERRA, ERA-Interim, NCEP-CFSR, JRA-25, ERA-40, NCEP-1, NCEP-2, and 20CR). Multiple linear regression is applied to the zonal and monthly mean time series of temperature for two peri-ods, 1979–2009 (for eight reanalysis data sets) and 1958– 2001 (for four reanalysis data sets), by considering explana-tory factors of seasonal harmonics, linear trends, Quasi-Biennial Oscillation, solar cycle, and El Niño Southern Oscillation. The residuals are used to deﬁne the volcanic signals


Introduction
Explosive volcanic eruptions inject sulphur species to the stratosphere in the form of SO 2 and H 2 S which convert to H 2 SO 4 aerosols.These aerosols are then transported both vertically and horizontally into the stratosphere by the Brewer-Dobson circulation (Butchart, 2014), stay there to perturb the radiative budget on a timescale of a few years, and thus affect global climate (Robock, 2000).The stratospheric volcanic aerosol layer is heated by absorption of near-infrared solar radiation and upward longwave radiation from the troposphere and surface.In the troposphere, the reduced near-infrared solar radiation is compensated by the additional downward longwave radiation from the aerosol layer.At the surface the large reduction in direct shortwave radiation due to the aerosol layer is the main cause of net cooling there.
Stratospheric aerosol optical depth (AOD) is an indicator of volcanic eruptions that affect global climate and has been estimated from various information (e.g.Sato et al., 1993;Robock, 2000;Vernier et al., 2011).Since 1960 astronomical observations such as solar and stellar extinction and lunar eclipses have become available from both hemispheres, and since 1979 extensive satellite measurements have begun with the Stratospheric Aerosol Monitor (SAM) II on the Nimbus-7 satellite.Extending over a longer period, the global radiosonde network that provides global atmospheric (upper-air) temperature data has been operating since the 1940s, with improved spatial resolution since the late 1950s (Gaffen, 1994).Since 1979, global satellite temperature measurements have begun with the Microwave Sounding Unit (MSU) and Stratospheric Sounding Unit (SSU) instruments on the TIROS-N satellite and on the subsequent several National Oceanic and Atmospheric Administration (NOAA) satellites.Since 1998, the Advanced MSU-A (AMSU-A) instruments on several NOAA satellites have provided global temperature measurements.See, e.g.Christy et al. (2003), Wang et al. (2012), Wang and Zou (2014), Zou et al. (2014), and Nash and Saunders (2015) for these satellite temperature measurements.
Since the late 1950s, three major volcanic eruptions occurred that significantly affected global climate, which are Mount Agung (8 • S, 116 • E), Bali, Indonesia in March 1963, El Chichón (17 • N, 93 • W), Chiapas, Mexico in April 1982, and Mount Pinatubo (15 • N, 120 • E), Luzon, Philippines in June 1991.The volcanic explosivity index (VEI) of these eruptions are 6 for Mount Pinatubo, 5 for El Chichón, and 4 for Mount Agung (Robock, 2000).Free and Lanzante (2009) and Randel (2010) used homogenized radiosonde data sets while Santer et al. (2001) and Soden et al. (2002) used MSU satellite data to investigate the tropospheric and stratospheric temperature response to these eruptions.When extracting the volcanic signals, one needs a good evaluation, at the same time, of the components of El Niño Southern Oscillation (ENSO), Quasi-Biennial Oscillation (QBO), and 11-year solar cycle as well as seasonal variations and linear trends.Each of the above four studies used a variety of regression analyses.
An atmospheric reanalysis system provides a best estimate of the past state of the atmosphere using atmospheric observations with a fixed assimilation scheme and a fixed global forecast model (Trenberth and Olson, 1988;Bengtsson and Shukla, 1988).It is an operational analysis system at a particular time (e.g.1995 for the NCEP-1 system and 2009 for the JRA-55 system), which has been continuously improved with the main motivation being to improve the tropospheric weather prediction.Using a fixed assimilationforecast model to produce analyses of observational data that were previously analysed in the context of operational forecasting -hence the "re" in "reanalysis" -prevents artificial changes being produced in the analysed fields due to system changes.But, as described above, the observational data inputs still vary over the period of the reanalysis.Currently, there are about 10 global atmospheric reanalysis data sets available worldwide.Table 1 lists the reanalysis data sets considered in this study.It is known that different reanalysis data sets give different results for the same diagnostic.Depending on the diagnostic, the different results may be due to differences either in the observational data assimilated, the assimilation scheme or forecast model, or any combination of these (see, e.g.Fujiwara et al., 2012, for a list of some examples).It is therefore necessary to compare all (or some of the newer) reanalysis data sets for various key diagnostics for understanding of the data quality and for future reanalysis improvements (Fujiwara and Jackson, 2013).To be more specific to the current study, the major observational sources of atmospheric (upper-air) temperature are basically common for all the reanalysis data sets in Table 1 (except for the 20CR which only assimilated surface pressure reports).They are radiosondes and satellite microwave and infrared sounders (i.e.MSU, SSU, and AMSU-A).There are three components that do differ in different reanalysis systems: (1) detailed bias-correction or quality-control methods for the original observations before the assimilation, (2) the assimilation scheme, and (3) the forecast model.Thus, any differences in the analysis results in this study would be due to the differences in these components (except for the 20CR).
Recently, Mitchell et al. (2015) analysed temperature and zonal wind data from nine reanalysis data sets using a linear multiple regression technique during the period from 1979 to 2009 by considering QBO, ENSO, AOD as a volcanic index, and solar cycle, with a focus on the solar cycle response.However, the volcanic response shown by Mitchell et al. is a combined response due to the major eruptions over the period 1979-2009(i.e. El Chichón in 1982and Mount Pinatubo in 1991).
Investigation of climatic response to individual volcanic eruptions using multiple reanalysis data sets for the purpose of comparison and evaluation of reanalysis data sets is rather limited.For example, Harris and Highwood (2011) showed global mean surface temperature changes following the Pinatubo eruption using NCEP-1 and ERA-40 reanalysis data for comparison with their model experiments.Analysing all available reanalysis data sets for the 20th-century three major eruptions separately and for the region covering both troposphere and stratosphere will provide valuable information for model validation as well as on the current reanalysis data quality for capturing volcanic signals.Such an analysis would also be valuable when assessing one of the proposed geoengineering options, i.e. stratospheric aerosol injection to counteract global surface warming (e.g.Crutzen, 2006;Robock et al., 2013).
In the present study, we analyse zonal and monthly mean temperature data from nine reanalysis data sets to investigate the response to the Mount Agung, El Chichón and Mount a For the version of the operational analysis system that was used for the reanalysis.b The model horizontal resolution has increased in April 2011 in the NCEP-CFSR.
Pinatubo eruptions separately.The temperature response to the Mount Agung eruption is investigated using four reanalysis data sets (JRA-55, ERA-40, NCEP-1, and 20CR) that cover the period back to the 1960s.A multiple regression technique is used to remove the effects of seasonal variations, linear trends, QBO, solar cycle, and ENSO, and the residual time series is assumed to be composed of volcanic effects and random variations.The remainder of this paper is organized as follows.Section 2 describes the data sets and analysis method.Section 3 provides results and discussion.Finally, Sect. 4 lists the main conclusions.

Data and method
Monthly mean pressure-level temperature data from the nine reanalysis data sets listed in  2015).Therefore, 20CR is expected to show volcanic signals even though it did not assimilate upper-air temperature data.
The atmospheric forecast model of the 20CR is nearly the same as used in the NCEP-CFSR but with a lower resolution, and thus the NCEP-CFSR also included the same volcanic aerosols.None of the other reanalysis data sets included radiative forcing due to volcanic aerosols in the forecast model.See Mitchell et al. (2015) for further technical comparisons among different reanalysis data sets.For a complete description of each reanalysis, see the reference papers shown in Table 1.Table 1 also shows the period of data availability for each reanalysis data set.For a direct intercomparison, we define two analysis periods, namely, between 1979 and 2009 (31 years) for eight reanalysis data sets (all except ERA-40) and between 1958 and 2001 (44 years) for four reanalysis data sets (JRA-55, ERA-40, NCEP-1, and 20CR).The former covers the eruptions of El Chichón in 1982 and Mount Pinatubo in 1991, while the latter also covers the eruption of Mount Agung in 1963.Results from JRA-55, NCEP-1, and 20CR for the El Chichón and Mount Pinatubo eruptions for the two different-period analyses also provide an opportunity to investigate sensitivity to the choice of analysis period.
A multiple regression technique is applied to extract volcanic signals (e.g.Randel and Cobb, 1994;Randel, 2010;von Storch and Zwiers, 1999, Chapt. 8.4).First, all major variabilities, except for volcanic effects, were evaluated and subtracted from the original zonal and monthly mean temperature data.The major variabilities include seasonal harmonics of the form, a 1 sin ωt +a 2 cos ωt +a 3 sin 2ωt +a 4 cos 2ωt + a 5 sin 3ωt + a 6 cos 3ωt, with ω = 2π/(12 months), linear trends, two QBO indices, ENSO, and solar cycle.For the latter five climatic indices, the six seasonal harmonics and a constant are further considered to construct seven indices for each of the five indices, as was done by Randel and Cobb (1994).For the two QBO indices, we use 20 and 50 hPa monthly mean zonal wind data taken at equatorial radiosonde stations provided by the Freie Universität Berlin.The crosscorrelation coefficient for these two QBO indices is −0.24 for 1979-2009 and −0.21 for 1958-2001.For the ENSO index, we use the Niño 3.4 index, which is a standardized sea surface temperature anomaly in the Niño 3.4 region (5 • N-5 • S, 170-120 • W), provided by the NOAA Climate Prediction Center.As is often done, a time lag for atmospheric response is considered for the ENSO index.We chose 4 months for the lag, following Free and Lanzante (2009).We confirmed that changing the ENSO lag from 0 to 6 months gives somewhat different ENSO signals particularly in the tropical stratosphere but does not alter other signals, including volcanic signals, significantly.For the solar cycle index, we use solar 10.7 cm flux data provided by the NOAA Earth System Research Laboratory.These climate indices are those considered by Free and Lanzante (2009), Randel (2010), andMitchell et al. (2015), though Free and Lanzante did not consider solar cycle and Mitchell et al. considered the AOD as well.(Note that we do not consider other indices, e.g. the North Atlantic Oscillation index and the Indian Monsoon index because the former is considered to be a response not a forcing and both are considered to be more related to regional response, not zonal mean response.)The multiple regression model that we use in this study is therefore where Y (t) is the zonal and monthly mean temperature time series at a particular latitude and pressure grid point, and a l is the least squares solution of a parameter for climatic index time series x l (t).R(t) is the residual of this model which is assumed to be composed of volcanic signals and random variations (Randel, 2010;Mitchell, 2015).Mitchell (2015) analysed two reconstructions of the SSU data set using model-predicted responses to external forcings as the climatic indices.After regressing the model-predicted response patterns onto observations, it was shown that the noise residual was very small compared with the forcing signal.If the volcanic predictor had been omitted (as in our study), the residual would essentially be the volcanic pattern.Finally, by following Randel (2010), the volcanic signal for each eruption is defined as the difference between the 12month averaged R(t) after each eruption and the 36-month averaged R(t) before each eruption.
There are several other possible minor variations for the methodological details, i.e. for the multiple regression model, the choice of particular index data sets, and the volcanic signal definition.The use of a consistent methodology is important for comparisons of different data sets.Where possible, however, we will discuss the methodological dependence below.

The 1979-2009 analysis
Figures 1 and 2 show temperature variations in association with the QBO, solar cycle and ENSO from JRA-55 and MERRA, respectively, for the region from 1000 to 1 hPa.The coloured regions are those evaluated as statistically significant at the 95 % confidence level (von Storch and Zwiers, 1999, Chapt. 8.4.6), with an effective degree of freedom where data are assumed to be independent for every 3 months.Comparing with the results from Mitchell et al. (2015) who used a regression analysis with different details, the setting of this effective degree of freedom may be somewhat too conservative.This is because the regions evaluated as statistically significant are smaller than those in Mitchell et al. (2015) particularly for the solar and ENSO signals in the tropical lower stratosphere, but the general features are quite similar to those shown in Mitchell et al. (2015) although they also considered a volcanic index in the multiple regression analysis.The two QBO variations are displaced vertically by a quarter cycle in the tropics because of their downward phase propagation.The temperature QBO has off-equatorial out-of-phase signals centred around 30 • N and around 30 • S because of the associated secondary meridional circulation (Baldwin et al., 2001).The major response to the solar cycle is the tropical lower stratospheric warming.The ENSO response includes the tropical tropospheric warming and a hint of tropical stratospheric cooling, although the statistical significance of this latter signal is weak.The strength of this cooling signal is sensitive to the choice of the time lag for the ENSO index (4 months in this study and 0 months in Mitchell et al., 2015).There also exists midlatitude lower stratospheric warming in both hemispheres for ENSO.The signals of QBO, solar cycle, and ENSO in the other six reanalysis data sets (ERA-Interim, NCEP-CFSR, JRA-25, NCEP-1, NCEP-2, and 20CR; not shown) are also similar to those in Mitchell et al. (2015).20CR shows no QBO signals (and no zonal-wind QBO; not shown) and no tropical stratospheric solar response.NCEP-CFSR shows weaker tropical lower stratospheric solar cycle warming.The overall agreement with the results in Mitchell et al. (2015), in addition to the finding by Mitchell (2015) as described in Sect.2, supports the assumption that the residual R(t) is composed of volcanic signals and random variations.
Figure 3 shows the residual time series averaged for 30 • N-30 • S at 50 and at 300 hPa together with the lower-tomiddle stratospheric AOD time series averaged for 27.4 • N-27.4 • S provided by the NASA Goddard Institute for Space Studies (Sato et al., 1993).The AOD time series clearly shows the timing of the El Chichón eruption and Mount Pinatubo eruption and the duration of their impact on the stratospheric aerosol loading.At 50 hPa, all reanalysis data sets show 1-2 K peak warming within 1 year after the El Chichón eruption, and most (except 20CR and JRA-25) show 2-2.5 K peak warming within 1 year after the Mount Pinatubo eruption.As described in Sect.2, 20CR does not assimilate upper-air data, but incorporates volcanic aerosols in the forecast model.Thus, 20CR shows a warming signal in association with both eruptions, though the one for Mount Pinatubo is smaller and slower.20CR also shows warming signals in 1989 and in 1990 though none of the other data sets show the corresponding signals.The warming in JRA-25 is ∼ 1 K smaller than other reanalysis data sets except 20CR.This cold bias can be seen at least during the period 1988-1994.This might in part be related to the known stratospheric cold bias in JRA-25 (Onogi et al., 2007).The radiative scheme used in the JRA-25 forecast model has a known cold bias in the stratosphere, and the TOVS SSU/MSU measurements do not have a sufficient number of channels to correct the model's cold bias; after introducing the ATOVS AMSU-A measurements in 1998, such a cold bias disappeared in the JRA-25 data product.It is also possible that the cold bias in JRA-25 during the TOVS era was not constant over time, in particular when unusual, volcanically affected temperature measurements came into the JRA-25 system, which could contribute to the smaller warming signals in our data analysis.As described in Sect.2, except for 20CR, NCEP-CFSR is the only reanalysis that included stratospheric volcanic aerosols in the forecast model, but no clear difference is found in comparison with other recent reanalysis data sets.At 300 hPa, all reanalysis data sets show 0.4-0.8K peak cooling within 1 year after the Mount Pinatubo eruption.No clear signals are found at 300 hPa for the El Chichón eruption.Note that the standard deviation (SD) of the residual time series is ∼ 1 K for tropical 50 hPa and ∼ 0.3 K for tropical 300 hPa for all the data sets; thus, the volcanic signals discussed above are distinguishable from random variations in the sense that these signals are much greater than one SD of the residuals.
Figure 4 shows the temperature signals for the El Chichón eruption from the eight reanalysis data sets.As described in Sect.2, the volcanic signal is defined as the difference between the 12-month averaged R(t) after each eruption and the 36-month averaged R(t) before each eruption.The coloured regions are also defined by following Randel (2010), i.e. as those regions with positive (negative) values more (less) than twice the SD of annual mean residual R(t).The annual mean is taken here because of the use of 12-month average in the volcanic signal definition.For the recent four reanalysis data sets, i.e.JRA-55, MERRA, ERA-Interim, and NCEP-CFSR, the tropical lower stratospheric warming of 1.2-1.6K centred around 50-30 hPa is a common signal.There are also Northern Hemisphere highlatitude middle-upper stratospheric warming and tropical upper stratospheric cooling signals, though the latter is comparable to random variations in some of the four data sets and thus its statistical significance is weak.The tropical and midlatitude troposphere is only weakly cooling, with a maximum cooling (0.4-0.8 K) occurring in the upper troposphere at 20-30 • N.For JRA-25, the tropical lower stratospheric warming is confined around 100-50 hPa with (statistically insignifi- cant) cooling signals around 50-10 hPa.This might be in part related to the cold bias in JRA-25 as described in the previous paragraph.The tropospheric features in JRA-25 are similar to those in the latest four reanalysis data sets.For NCEP-1 and NCEP-2, the tropical stratospheric warming region extends to 10 hPa where it maximizes, and the 20-30 • N upper tropospheric cooling is largely missing.The major differences of the NCEP-1 and NCEP-2 systems from the recent four reanalysis systems include the lower model top height (3 hPa), older forecast model and assimilation scheme (of the 1990s; see Table 1), and the use of retrieved temperature data for the assimilation of SSU, MSU, and AMSU-A data.It is possible that these factors may be responsible for the different signals of the El Chichón eruption in NCEP-1 and NCEP-2.(See also discussion on the results for the Mount Pinatubo eruption below).For 20CR, tropical stratospheric warming is present, but again, this is due to the specified volcanic aerosols in the forecast model.Free and Lanzante (2009) and Randel (2010) analysed the temperature signals for the El Chichón eruption using different homogenized radiosonde data sets globally up to the 30 hPa level.The distribution of the tropical lower stratospheric warming signal is similar, though the peak warm-ing is greater, i.e. 1.6-2 K for Free and Lanzante (2009, their Fig. 3) and 2.5-3 K for Randel (2010, his Fig. 4).(Note that Free and Lanzante defined the volcanic signals as the difference between the 24-month average after the eruption and the 24-month average before the eruption, but we use the same definition of volcanic signals as Randel (2010) and still obtain roughly a factor of two discrepancy in tropical lower stratospheric warming (1.2-1.6K from the reanalyses versus 2.5-3 K from the radiosondes)).Free and Lanzante ( 2009) also show a 20-30 • N upper tropospheric cooling of 0.6-0.9K.
Figure 5 shows the temperature signals for the Mount Pinatubo eruption.For the latest four reanalysis data sets, i.e.JRA-55, MERRA, ERA-Interim, and NCEP-CFSR, the tropical lower stratospheric warming of 2.0-2.8K (depending on data sets) centred around 50-30 hPa is a common signal.In the upper troposphere, a cooling (0.4-0.8 K) at 20-30 • N and at 15-45 • S can be seen, with the latter somewhat greater.JRA-25 shows similar upper tropospheric features and relatively similar lower stratospheric features, though for the latter, the warming magnitude is smaller and the "random" variability becomes large above the 50 hPa level because of the reason described above (i.e. the cold bias and its disappear- ance in 1998).For NCEP-1 and NCEP-2, the tropical tropospheric and stratospheric features are similar to those for the latest four reanalysis data sets, though the lower stratospheric warming magnitude is slightly smaller than in most of the other reanalyses.Comparing with the El Chichón case, the NCEP-1 and NCEP-2 systems worked much better to capture the Mount Pinatubo signals for some reasons.For 20CR, the tropical stratospheric warming is not detected.This is because of the unknown warming signals in 20CR in 1989 and in 1990 (see Fig. 3) that raised the 36-month averaged base in the volcanic signal definition.As in Fig. 3, there are no relevant signals in AOD around 1989-1990.Thus, the unknown warming signals are likely due to unrealistic (unforced) variations in the 20CR system.
The temperature signals for the Mount Pinatubo eruption shown in Randel (2010) are similar to the present results both in the tropical-midlatitude stratosphere and troposphere, though Randel's stratospheric warming peak value is somewhat greater (∼ 3 K) and his upper tropospheric cooling is somewhat greater (0.5-1 K) and more uniform in latitude.On the other hand, Free and Lanzante (2009) show that the lower stratospheric warming signal is split near the equator with two maxima (1.6-2 K at 10 • N and > 2 K at 15 • S, both at 70-50 hPa) and that the upper tropospheric cooling signal has its peak (0.9-1.2 K) around 20 • S. In summary, the recent four reanalysis data sets (i.e.JRA-55, MERRA, ERA-Interim, and NCEP-CFSR) give more consistent signals for both eruptions compared to the two radiosonde data analyses using different homogenized data sets by Free and Lanzante (2009) and Randel (2010).

The 1958-2001 Analysis
The multiple regression analysis is applied to the four reanalysis data sets, namely, JRA-55, ERA-40, NCEP-1, and 20CR which cover the period of 1958-2001.Figure 6 shows temperature variations associated with the QBO, solar cycle, and ENSO from JRA-55.Comparing with the 1979-2009 analysis results shown in Fig. 1, all variations are quite similar, with the statistically significant regions for the solar cycle variation being much greater both in the tropical stratosphere and in the tropical troposphere.The same is true for NCEP-1 (not shown).20CR does not have QBO and stratospheric solar-cycle signals, but does show ENSO signals in both 1979ENSO signals in both -2009ENSO signals in both and 1958ENSO signals in both -2001 analyses; the 20CR ENSO signals are similar to those from all other reanalysis data sets.ERA-40 shows similar results to JRA-55 except for the solar cycle variation.In ERA-40, the tropical lower stratospheric warming signal in association with the solar cycle is very weak and not symmetric about the equator, in contrast to the results by Crooks and Gray (2005) and Mitchell et al. (2015) who both applied a regression analysis during the period 1979-2001.
Figure 7 shows the time series of residual R(t) and stratospheric AOD averaged over the tropics for the period between 1958 and 2001.The AOD time series shows the timing of the Mount Agung eruption in March 1963 as well as the El Chichón and Mount Pinatubo eruptions.The features at both 50 and 300 hPa for the El Chichón and Mount Pinatubo eruptions are quite similar to the 1979-2009 analysis results shown in Fig. 3, including the 20CR's smaller and slower Mount Pinatubo signal at 50 hPa.For the Mount Agung eruption, ∼ 2.5 K peak warming is seen within 1 year after the eruption except for 20CR.At 300 hPa, a sudden cooling occurred about 1 year later, i.e. in mid-1964 for all the data sets, which is probably related to the Mount Agung eruption.The cooling might have continued for more than 1 year.ERA-40 shows anomalous ∼ 1 K warming in the mid-1970s at both levels, which are not present in other reanalysis data sets (see also Fig. 14 of Kobayashi et al., 2015).The AOD time series in Fig. 7 shows a small increase in the mid-1970s which is probably due to the eruption of Mount Fuego (14 • N, 91 • W), Guatemala, in October-December 1974 (VEI 4, Smithsonian Institution National Museum of Natural History Global Volcanism Program, http://www.volcano.si.edu/, last accessed August 2015).The magnitude and the sign, however, (i.e.warming) at 300 hPa seem unrealistic.Before the introduction of horizontally dense satellite measurements in 1979, the upper-air temperature is constrained basically only by horizontally inhomogeneous, relatively sparse radiosonde data (see, e.g.Fig. 2 of Uppala et al., 2005).Also, the ERA-40 system is a relatively old system (the 2001 version of the ECMWF analysis system).These two facts are possible reasons for the ERA-40's anomalous warming in the mid-1970s.A stream change of the reanalysis execution could also be a potential reason.For the ERA-40, there were three execu- tion streams, that is, 1989-2002, 1957-1972, and 1972-1988(Uppala et al., 2005)).But the stream change point of 1972 is unlikely to explain the anomalous warming starting around the end of 1974.
Figure 8 shows the temperature signals for the Mount Agung eruption from four different reanalysis data sets.All except 20CR show Southern Hemisphere lower stratospheric warming centred at 40-30 • S and 100-50 hPa, with an extension to equatorial latitudes at 50 hPa.The maximum warming value varies with data set, that is, 1.6-2 K for NCEP-1, 2-2.4 K for JRA-55, and 2.4-2.8K for ERA-40.The reason for the weak signal in 20CR is in that 20CR does not assimilate upper-air temperature observations but does consider volcanic aerosol loading in the forecast model.The modelled aerosol loading was probably too weak to simulate the lower stratospheric warming signals.For all four reanalysis data sets, the 300 hPa cooling shown in Fig. 7 is not captured with the current volcanic-signal definition (i.e.12-month average after the eruption started).
Free and Lanzante (2009) showed a very similar Southern Hemisphere midlatitude lower stratospheric warming signal (> 2 K) in association with the Mount Agung eruption using a homogenized radiosonde data set.Sato et al. (1993) showed that the aerosols emitted from the Mount Agung eruption Figure 9 provides a useful summary plot for the volcanic effects on the temperature at 50 hPa and at 300 hPa using JRA-55 from the 1958-2001 analysis together with the AOD latitudinal time series.The aerosol loading due to the Mount Agung eruption in March 1963 extended primarily to the Southern Hemisphere, that due to the El Chichón eruption in April 1982 was very large in the tropics and extended primarily to the Northern Hemisphere, and that due to the Mount Pinatubo eruption in June 1991 was very large in the tropics and extended to both hemispheres.The tropical lower stratosphere warmed after these three major volcanic eruptions with a timescale of 1-2 years.The warming after the Mount Agung eruption is not equatorially symmetric and is shifted to the Southern Hemisphere and to somewhat lower levels, in association with the distribution of aerosol loading.The tropical troposphere became cooler after the Mount Pinatubo eruption but the tropospheric response is not as clear for the other two eruptions.The high-latitude response is also unclear both in the troposphere and stratosphere due to high random variations that mask any volcanic signals, if they exist.

Conclusions
Monthly and zonal mean temperature data from nine reanalysis data sets were analysed to characterize the response to the three major volcanic eruptions during the 1960s to the 1990s.Multiple linear regression analysis was applied to evaluate seasonal variations, trends, QBO, solar cycle and ENSO components, and the residual time series R(t) was assumed to be composed of volcanic signals and random variations.The volcanic signals were defined as the difference between the 12-month averaged R(t) after each eruption and the 36-month averaged R(t) before each eruption.Two separate analyses were performed, that is, one for the period 1979-2009 (31 years) using eight reanalysis data sets and the other for 1958-2001 (44 years) using four reanalysis data sets.The former covered the eruptions of El Chichón (April 1982) and Mount Pinatubo (June 1991), while the latter also covered the eruption of Mount Agung (March 1963).
The general features of the response to QBO, solar cycle, and ENSO were found to be quite similar to those shown in Mitchell et al. (2015) who also used a multiple linear regression with different methodological details, in particular, considering a volcanic index as well.Also, these signals were at least qualitatively similar among reanalysis data sets, with a notable exception that 20CR shows no QBO signals and no tropical stratospheric solar response.
The latitude-pressure distribution of El Chichón and Mount Pinatubo temperature response was quite similar at least among the recent four reanalysis data sets (JRA-55, MERRA, ERA-Interim, and NCEP-CFSR) and between the 1979-2009 and 1958-2001 analyses.For the Mount Pinatubo eruption, tropical lower stratospheric warming and tropical upper tropospheric cooling were observed.For the El Chichón eruption, tropical lower stratospheric warming was observed, but tropospheric cooling was much weaker than the Mount Pinatubo case.For the Mount Agung eruption, JRA-55, ERA-40, and NCEP-1 showed Southern Hemisphere lower stratospheric warming centred at 40-30 • S and 100-50 hPa, with an equatorial extension to 50 hPa.Thus, the Agung signal was asymmetric about the equator and very different from the El Chichón and Pinatubo signals.We suggest that this may be due to differences in the transport of volcanic aerosols (Sato et al., 1993).
Evidently the temperature responses were different for different volcanic eruptions.In particular, wide-spread upper tropospheric cooling was observed only for the Mount Pinatubo case, and the Mount Agung lower stratospheric response was found to be asymmetric about the equator.The characteristics in the temperature response are related to the transport of stratospheric aerosols together with the amount of sulphur species emitted into the stratosphere.Depending on the location, season, and magnitude of the eruption, the climatic response can be very different (e.g.Trepte and Hitchman, 1992).This needs to be taken into account when evaluating the stratospheric sulphur injection as a geoengineering option, and thus accurate estimations of stratospheric circulation and transport are essential for assessing the climate impacts.Also, it should be noted that accurate evaluation of naturally induced variability such as QBO, solar cycle, and ENSO is necessary to detect the effects of artificial injection.
Finally, we conclude that the four most recently developed reanalysis data sets, i.e.JRA-55, MERRA, ERA-Interim, and NCEP-CFSR are equally good for studies on the response to the El Chichón and Mount Pinatubo eruptions.The NCEP-1, NCEP-2, and JRA-25 showed different tropical stratospheric signals particularly for the El Chichón eruption, though the original upper-air temperature observations assimilated are basically common, and this is most probably in association with the use of older analysis systems.The 20CR did not assimilate upper-air observations and gives very different volcanic signals, despite including volcanic aerosols in the forecast model.Of the currently available data sets that extend back far enough (JRA-55, ERA-40, NCEP-1, and 20CR) the JRA-55 data set is probably the most ideally suited for studies of the response to the Mount Agung eruption because it is the only data set that employs the most recent reanalysis system.

Figure 1 .
Figure 1.Latitude-pressure distribution of the temperature variations in association with (top left) QBO 20 hPa zonal wind index, (top right) QBO 50 hPa zonal wind index, (bottom left) solar cycle index, and (bottom right) ENSO index from JRA-55 reanalysis data for the period 1979-2009.The units are in Kelvin per standard deviation (SD) of each index (note that each index time series was standardized before the regression analysis).Solid and dashed lines denote positive and negative values, respectively.The contour interval is 0.2 K for QBO, and 0.1 K for solar cycle and ENSO.Coloured regions denote those greater (orange) and smaller (blue) than random variations with the 95 % confidence interval at each location.

Figure 3 .
Figure 3.Time series of temperature residual R(t) (including volcanic signals and random variations) averaged for 30 • N-30 • S for the 1979-2009 regression analysis from eight reanalysis data sets at (a) 50 hPa and (b) 300 hPa.(c) Time series of aerosol optical depth at 550 nm averaged for 27.4 • N-27.4 • S and integrated for the region 15-35 km.Vertical dotted lines indicate the starting date of the two volcanic eruptions.

Figure 4 .
Figure 4. Latitude-pressure distribution of the temperature response to the El Chichón eruption in April 1982 for the 1979-2009 analysis from eight reanalysis data sets.Solid and dashed lines denote positive and negative values, respectively.The contour interval is 0.4 K. Coloured regions denote those with positive and greater (orange) and negative and smaller (blue) than twice the SD of annual mean residual R(t) at each location.

Figure 5 .
Figure 5.As in Fig. 4 but for the Mount Pinatubo eruption in June 1991.

Figure 7 .
Figure 7.As in Fig. 3 but for the 1958-2001 regression analysis from four reanalysis data sets.Vertical dotted lines indicate the starting date of the three volcanic eruptions.

Figure 8 .
Figure 8.As in Fig. 4 but for the Mount Agung eruption in March 1963 for the 1958-2001 analysis from four reanalysis data sets.

Figure 9 .
Figure 9. Time-latitude distribution of temperature residual R(t) (including volcanic signals and random variations) for the 1958-2001 regression analysis from JRA-55 reanalysis data at (a) 50 hPa and (b) 300 hPa.Thirteen-month running average has been taken for R(t).The contour interval is 1.0 K for (a) and 0.25 K for (b).The regions with 0-1 K (> 1 K) are coloured in orange (red) in (a).The regions with 0 to −0.25 K (< −0.25 K) are coloured in light (dark) blue.(c) Time-latitude distribution of aerosol optical depth at 550 nm integrated for the region 15-35 km.The contour interval is 0.04.The regions with 0.04-0.12(> 0.12) are coloured in orange (red) in (c).

Table 1 .
List of global atmospheric reanalysis data sets considered in this study.

Table 1
(1993)were translated to the optical depth values for ultraviolet, visible, near infrared, and infrared spectral bands (Y.-T.Hou, personal communication,