The global precipitation response to volcanic eruptions in the CMIP5 models

We examine the precipitation response to volcanic eruptions in the Coupled Model Intercomparison Project Phase 5 (CMIP5) historical simulations compared to three observational datasets, including one with ocean coverage. Global precipitation decreases significantly following eruptions in CMIP5 models, with the largest decrease in wet tropical regions. This also occurs in observational land data, and ocean data in the boreal cold season. Monsoon rainfall decreases following eruptions in both models and observations. In response to individual eruptions, the ITCZ shifts away from the hemisphere with the greater concentration of aerosols in CMIP5. Models undergo a longer-lasting ocean precipitation response than over land, but the response in the short satellite record is too noisy to confirm this. We detect the influence of volcanism on precipitation in all three datasets in the cold season, although the models underestimate the size of the response. In the warm season the volcanic influence is only marginally detectable.


Introduction
Global precipitation has been found to decrease for a couple of years following large explosive volcanic eruptions (e.g. Robock and Liu 1994, Trenberth and Dai 2007, Schneider et al 2009, Gu et al 2007, Gu and Adler 2011, Broccoli et al 2003, Iles et al 2013. These eruptions inject SO 2 into the stratosphere, where it is converted to sulphate aerosols. The aerosols spread out globally over subsequent months following tropical eruptions, or over the hemisphere of eruption for high latitude eruptions, and reflect incoming solar radiation (e.g. Robock 2000, Timmreck 2012 andreferences therein). This causes widespread surface and tropospheric cooling, lasting a few years, along with changes in atmospheric circulation and precipitation (e.g. Robock and Liu 1994, Robock 2000and references therein, Gillett et al 2004, Trenberth and Dai 2007, Joseph and Zeng 2011, Driscoll et al 2012, Timmreck 2012. This decrease in global mean precipitation is due to a reduction in short-wave radiation reaching the surface, reducing evaporation, stabilizing the atmosphere and reducing the saturation mixing ratio of the air (Bala et al , Cao et al 2012. In addition, a cooler atmosphere undergoes less radiative cooling to space, which allows less condensation and precipitation to occur Ingram 2002, O'Gorman et al 2012). Circulation changes following eruptions modulate this global decrease on a regional scale, e.g. in monsoon regions (e.g. Schneider et al 2009, Joseph and Zeng 2011, Peng et al 2010, Cao et al 2012. There is also a tendency towards a positive phase of the North Atlantic Oscillation (NAO) in the winter following low latitude eruptions with associated winter warming over northern hemisphere continents, although this response is relatively noisy (e.g. Robock and Mao 1992, 1995, Hegerl et al 2011. The NAO response is not well captured by climate models (Stenchikov et al 2006, Driscoll et al 2012, Charlton-Perez et al 2013.

Environmental Research Letters
Environ. Res. Lett. 9 (2014) 104012 (10pp) doi:10.1088/1748-9326/9/10/104012 Content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. Iles et al (2013) (hereafter I13) investigated the precipitation response to volcanic eruptions using last millennium simulations with HadCM3 and compared the response for twentieth century eruptions to observational land precipitation data. They found a significant reduction in global and wet tropical regions precipitation in both the model and observations, whilst dry tropical ocean regions got significantly wetter in the model. This dry get wetter and wet get drier pattern is the opposite of the global warming response and is linked to changes in transport associated with the Hadley circulation (e.g. Soden 2006, Trenberth 2011). Monsoon regions dried (see also Joseph and Zeng 2011, Schneider et al 2009, Trenberth and Dai 2007, with the exception of SE Asia in the observations. The precipitation response lasted longer over ocean than land in the model (see also Joseph and Zeng 2011). The modelled influence of volcanic forcing was detectable in the observations in the boreal cold season, although the magnitude appeared to be underestimated by the model. In contrast, the response to volcanic eruptions was only marginally detectable in the warm season. The model underestimate agreed with previous detection studies that identified a detectable land precipitation response (Lambert et al 2004, 2005, Gillett et al 2004.
Here we investigate whether the main findings of I13 are consistent with results using the CMIP5 models, many of which have higher horizontal and vertical resolutions and extend higher into the stratosphere than HadCM3. Furthermore, whilst I13 only used observational data over land, here we use an additional satellite-gauge dataset to investigate whether the long-lasting ocean response found in HadCM3, along with the wettening response in the dry tropical ocean regions are supported by observations. Finally we examine whether the CMIP5 models underestimate the precipitation response to volcanism, testing sensitivity to using alternative observational datasets.

Model data
We downloaded all twentieth century historical runs from the CMIP5 archive that were available in December 2012 and contain volcanic forcing (see table S1). These are also forced by observed records of other natural and anthropogenic forcings, such as solar variability, greenhouse gases and land use change (see Taylor et al (2012) for details). There is no recommended volcanic forcing dataset, but many modelling groups use Ammann et al (2003), its update (Ammann et al 2007) or an updated version of Sato et al (1993) (available at data.giss.nasa.gov/modelforce/strataer) (see table S1). We did not include HadCM3 in the analysis, for sake of comparison with I13. Multi-model means are constructed by averaging over all 88 available runs, even where there are differing numbers of runs per model. Using only a single simulation for each model yielded qualitatively similar results.

Observational data
As in I13. we used two land precipitation datasets based on gauge data, and additionally a shorter combined satellitegauge record. There are substantial differences in the methods of construction between the datasets, allowing some assessment of the robustness of results to observational uncertainty. The first is the Global Precipitation Climatology Centre's (GPCC) Full Data Reanalysis Version 6 (Becker et al 2013), which is a 2.5 × 2.5°gridded global land precipitation dataset extending from 1901-2010. It is spatially interpolated, resulting in full land coverage and is based on a very large number of station records (67 200 stations with at least ten years of data).
The second is an updated version of that detailed in Zhang et al (2007) (hereafter Z07), used in I13. This is a 5 × 5°gridded dataset covering the period 1900-2009. It is based on a subset of stations from the Global Historical Climatology Network's second dataset (Vose et al 1992) that have at least 25 years of data in the 1961-1990 base period and at least five years of data in every decade from 1950-1999. Z07 does not rely on spatial interpolation, yielding a dataset with less spatial coverage, but that is more homogeneous and constrained by stations than GPCC.
Finally we used the 2.5 × 2.5°version of the Global Precipitation Climatology Project (GPCP) combined satellitegauge dataset (Adler et al 2003). This dataset is spatially complete and, unlike the other two datasets, includes ocean coverage. The dataset begins in 1979, since when there have been only two major eruptions, and only one since the introduction of a microwave based sensor in 1987 which has improved retrieval accuracy relative to the pre-microwave era.

Epoch analysis
As in I13 we used 'epoch analysis', which involves averaging across the precipitation response to several volcanic eruptions in order to reduce internal variability. We used the five largest eruptions since 1900, as defined by global mean aerosol optical depth (AOD): the 1902 Santa Maria eruption: Novarupta in 1912: Agung 1963: El Chichon 1982and Pinatubo 1991 figure 3(b)). For each grid cell and eruption, anomalies for each of the ten years following an eruption, or up until the next eruption if that occurred first, were calculated with respect to a five year pre-eruption mean, to account for multidecadal changes in precipitation due to long-term trends or low frequency variability. We then averaged across all the eruptions available for each grid cell before spatial averaging. Where an eruption occurred too close to the beginning of a dataset (e.g., 1902), a shorter preeruption mean was used (see SI). Following I13, we use half year seasons: the boreal cold season (November to April, NDJFMA), and boreal warm season (May-October, MJJASO). 'Year 1' denotes the season in question commencing after the eruption date, in order to give aerosols time to spread globally (table S2). Where annual data are used for epoch analysis, year 1 is defined as starting three complete months after the eruption date (table S2).
As in I13, significance of results was tested using a Monte-Carlo technique, in which the analysis is repeated 10 000 times using randomly selected years as eruption years (except for CMIP5 masked to Z07, and composite maps, where 1000 cycles were used due to computational constraints). 5-95% Confidence intervals were then calculated from the distribution of these results (see I13 for more detail). Since GPCP has very limited temporal coverage (33 years), we obtained a second set of confidence intervals for it using the CMIP5 runs: for each region we combined the last 33 years of each run into one long time series after converting to anomalies and rescaling each run's standard deviation to the standard deviation of GPCP (which tends to be larger over oceans). We then performed Monte-Carlo analysis on this new time series (see SI for detail).
As in I13, where wet or dry tropical regions are referred to, wet regions are defined as the wettest third of grid cells between 40°N and 40°S, and the dry regions are the remaining two thirds based on climatological precipitation (using the base period 1961-1990 for CMIP5 and Z07, the normals supplied with the dataset for GPCC (see Becker et al 2013) and 1979-2010 for GPCP). These regions are fixed through time (see I13) and are defined separately for each season, model run and observational dataset and after any masking according to land or ocean, or observational data coverage.

Removing the influence of ENSO
The El Nino Southern Oscillation (ENSO) is associated with reduced precipitation over land and increased precipitation over the ocean for its positive phase, and vice versa for la Nina (e.g. Gu et al 2007, Gu and Adler 2011, Liu et al 2012. The 1991Pinatubo, 1982El Chichon, and 1963 Agung eruptions were all followed by El Nino events, albeit a weak one for Agung. Therefore, we repeat the analysis for the observational precipitation datasets after linearly removing the influence of ENSO. As in I13 we use the cold tongue index (CTI) as a measure of ENSO variability and calculate a regression coefficient for each grid cell and season between the detrended time series of CTI and precipitation, avoiding years 0-5 following an eruption. Note that this leaves only a limited number of ENSO events for the regression, particularly for GPCP, hence the removal of ENSO is somewhat noisy in that case (excluding fewer years following eruptions for GPCP did not improve results). We then subtract ENSO related precipitation from the precipitation time series at each gridpoint to arrive at a precipitation dataset with ENSO influence (at least partially) removed. We did not remove ENSO from the models since its signal should average out across the large number of simulations when constructing the multi-model mean. Figure 1 shows the time series of twentieth century precipitation for 50°N-50°S for land (b), ocean (c) and the two combined (a) for the CMIP5 models compared to observations. Latitudes poleward of 50 o are excluded to avoid biases in GPCP over the high latitudes in winter, particularly over oceans (Adler et al 2012), although results using the whole globe are very similar (not shown). A clear decrease in precipitation following volcanic eruptions can be seen in the multi-model mean, particularly over oceans; and over land and ocean combined. This is also reflected in a lowering of the ensemble envelope. Over land the modelled response is noisier but still visible. The observed response is less clear, although a decrease in land precipitation can be seen following the 1991, 1982 and 1912 eruptions in all datasets. The observed datasets appear well correlated over land over the more recent period, but agree less well further back in time (see also Polson et al 2013). Over ocean GPCP shows a noisy  Australian, African and North and South American monsoon regions dry in their respective warm seasons, whilst surrounding areas get wetter and the extratropics get drier on average. There appears to be a southward shift of the ITCZ over the Atlantic and Pacific oceans in MJJASO. Neither the multi-model mean (figure 2(a)) nor individual models (figure S1) show any evidence for a precipitation pattern suggestive of a positive NAO in NDJFMA, consistent with Driscoll et al (2012) and Charlton-Perez et al (2013). The modelled response lasts until year 2 over land and until year 3 (figure S2) or 4 (not shown) over ocean, although the response over the Pacific is less stable over time, as was also the case in HadCM3 in I13.

Results
The observed precipitation response patterns match well between observational datasets, despite GPCP only covering two eruptions (figures 2(c)-(h)). Based on the results presented in I13, observed patterns are expected to be only marginally significant. These observed patterns are of much greater magnitude than the multi-model mean, probably due to the cancellation of noise in the latter. As in the models, the monsoon regions get drier following eruptions in the observations, although the exact location of these drying areas is slightly different and the Asian monsoon regions show a mixed response. Unlike in the models, there is a positive NAO precipitation pattern in NDJFMA in all the observed datasets (see Fischer et al 2007). Removing the ENSO influence generally makes little difference to results ( figure  S3). Comparing observations and models over ocean is more difficult since GPCP results are noisy, but both show a wettening signal in the east equatorial pacific and south pacific convergence zone along with a drying signal in the location of the ITCZ over the Atlantic in both seasons.
Previous studies have identified shifts in the position of the ITCZ in response to volcanic forcing (e.g. Schneider et al 2009, Haywood et al 2013. The ITCZ tends to move away from the cooler hemisphere. In order to keep upper tropospheric temperature gradients small within the tropics, the branch of the Hadley cell in the cooler hemisphere strengthens, transporting heat from the warmer hemisphere to the cooler one. Moisture, which is concentrated in the lower troposphere, is transported in the opposite direction, causing the ITCZ to shift away from the cooler hemisphere (Frierson and Hwang 2012, Hwang et al 2013, Kang et al 2008. For an eruption with a hemispherically symmetric aerosol cloud this causes the ITCZ to move less far into the summer hemisphere, since the summer hemisphere undergoes a greater volcanic cooling relative to the winter one (Yoshimori andBroccoli 2008, Schneider et al 2009). For asymmetric aerosol clouds Haywood et al (2013) found the ITCZ to shift away from the hemisphere with the greatest increase in AOD using HadGEM2-ES. We also find this shift in response to asymmetric forcing in CMIP5, for example a northward shift in both seasons following the southern hemisphere biased Agung eruption (figures S4 (a), (b)), and a southward shift following the northern hemisphere biased 1982 El Chichon (figures S4(c), (d)), 1902 Santa Maria and high latitude 1912 Novarupta eruptions (not shown), although these shifts are only clear in the multi model mean over ocean ( figure S4).
These shifts cause a smaller decrease, or even increase in hemispheric mean precipitation in the hemisphere with fewer aerosols (figure 3). As most of the twentieth century eruptions had stronger aerosol loadings in the NH, the average response resembles such eruptions (figure 2). Figure 4 shows the post-volcanic precipitation response to the two most recent eruptions in GPCP compared to CMIP5 for various regions (the extratropics and dry tropical land regions are shown in figure S5). Results for CMIP5 using all five eruptions are very similar, but more highly significant (not shown). CMIP5 results are very similar to those using HadCM3 in I13; there is a significant decrease in precipitation in the multi-model mean for the global mean and wet tropical regions, over both land and ocean, and a significant increase over dry tropical ocean regions. As was the case for HadCM3, the ensemble mean response over ocean lasts longer than that over land and is smaller in magnitude. Whilst this long ocean precipitation response cannot be seen in all individual ensemble members when averaged over two eruptions (figures 4(e)-(h) grey lines), figure 5 demonstrates that it can be seen in every model when its ensemble mean is taken and the average over five eruptions is calculated. I13 found that the precipitation response over ocean matched the timescale of near-surface air temperature response, suggesting that the ocean precipitation response is driven by changes in sea surface temperature. Over land, precipitation reacted faster than temperature, instead matching the timescale of AOD and a reduction in land-ocean temperature contrast, implying a directly forced component and possible contributions from weakening monsoons respectively.
Over land (figures 4(a)-(d)), there is a significant decrease in observed precipitation in NDJFMA both in the global mean and the wet tropics, independently of the statistical test employed. This result is not sensitive to the removal of ENSO (figure 4), to observational uncertainties (figure S6), or the number of eruptions included in the analysis (2 or 5) (figure S6). While the multi-model mean underestimates the observed response, the latter remains within the ensemble envelop. In MJJASO, the decrease in global and wet tropical regions precipitation is not robust.
Over ocean (figures 4(e)-(j)) GPCP precipitation only really follows model expectations in the cold season, including a significant decrease in the wet tropical regions and a significant increase in the dry tropical regions. Global mean ocean precipitation only decreases once the influence of ENSO is removed. Observed variability is greater in magnitude than that of CMIP5. In the warm season, observed precipitation increases for the global mean and wet tropics contrary to model expectations. An increase in precipitation can be seen in the dry regions as expected, but this is not robust to choice of statistical test or the removal of ENSO. Some of the secondary peaks several years after the eruption in the observations may be due to incomplete removal of ENSO from GPCP (see section 3.2). For example, the peaks in year 7 in NDJFMA coincide with the 1998 El Nino event, and the spatial patterns suggest some El Nino influence even after ENSO removal (not shown). Figure 6 shows the precipitation response in the monsoon regions in GPCC and CMIP5 averaged across five eruptions. Monsoon rainfall is of great importance to people living in these regions. Monsoon regions are defined following Hsu et al (2011) (see methods details in SI) and are shown in figure S7. These regions constitute a substantial part of the wet tropical regions used above, but also include some areas outside (figure S7). A significant decrease in monsoon rainfall can be seen in the southern hemisphere monsoon regions (South American, African and Australian monsoon regions) in austral summer in both the multi-model mean and observations (figure 6(a)). The observed response is not robust to the removal of ENSO. Whilst precipitation also decreases significantly in the multi model mean for the northern hemisphere monsoon regions in boreal summer (Asian, African and North American monsoon regions), the decrease seen in the observations is not significant ( figure 6(b)).
Finally we perform a detection analysis as described in I13 to determine whether the overall response is significant and whether or not the models and data are consistent (figure 7). We first split the global land response into the northern hemisphere extratropics, wet tropical and dry tropical regions to yield a 3 element vector for each year following the eruptions that characterizes the volcanic response. The southern hemisphere extratropics are excluded due to their limited land area. For GPCP a four element vector is used consisting of dry tropical ocean regions, dry tropical land regions, wet tropical ocean regions, and wet tropical land regions. Extratropics are not used in the main analysis for GPCP to avoid less reliable high latitude data (however, results are robust to including the extratropics, see figure S8).
For each year following the eruptions we regress this 3 (land data) or 4 (GPCP) element vector of the observed response averaged across all five eruptions (2 for GPCP) against the equivalent for the multi-model mean fingerprint. This yields a regression coefficient, or scaling factor which indicates whether the modelled response is bigger or smaller than that observed for each year post-eruption. For land data, the model data was first masked according to the data coverage of the observational dataset to which it is being compared. We test whether or not scaling factors are significantly different from those expected by chance through conducting a Monte-Carlo analysis using random years for the observations as discussed above and regressing against the multi-model mean fingerprint. We then calculate 5-95% confidence intervals for the scaling factors. When a scaling factor exceeds the 95th percentile, it is large enough to be unlikely to occur due to climate variability, and the volcanic influence is said to be detected. We repeat the analysis using individual ensemble members instead of the observations, regressing against the mean of the remaining ensemble members to establish whether or not the observed response is consistent with the models. Results are also presented for years 1 and 2 combined, since this is when the clearest response is expected based on the multi-model mean response (figures 4 and S5). Figure 7 shows that the influence of volcanism is detectable (5% significance level) in NDJFMA in year 1 and years 1 and 2 combined in all three datasets. The large regression coefficients in these years suggest that CMIP5 underestimates the magnitude of the response, as HadCM3 also did to a greater extent in I13. In these years less than 5% of the ensemble members exhibit a response as big as that observed in almost all cases, suggesting that the model underestimate is significant. In MJJASO the volcanic influence is only detectable at the 90% level in most datasets in these years. The magnitude of the observed coefficients agrees better with those of the models in MJJASO.
With the influence of ENSO regressed out from the observations, the volcanic influence is still detectable at the 95% level in NDJFMA in year 1 and years 1 and 2 combined in Z07 and GPCP, but not GPCC, which is based on more data but uses interpolation (see discussion in Polson et al 2013). The magnitude of the coefficients decreases when removing ENSO, suggesting that ENSO was partly responsible for the large scaling factors in figure 7(a). In MJJASO removing ENSO decreases detectability making results marginal.

Conclusions
In this study we investigate whether the main features of the precipitation response to large explosive volcanic eruptions found by I13 using HadCM3 are consistent with the response simulated by the CMIP5 models. We also extend their analysis by comparing the modelled response to a satellite-gauge dataset which includes ocean coverage, allowing us to test whether the long-lasting ocean precipitation response found in HadCM3, and the wettening response in the dry tropical ocean regions are supported by observations. Finally, we examine whether the model underestimate of the precipitation response found in I13 is also seen in the CMIP5 models, testing sensitivity to choice of dataset.
The main features of the precipitation response to volcanic eruptions in the CMIP5 models are consistent with those found in HadCM3 by I13. This includes a significant decrease in global, extratropical (see SI) and wet tropical regions precipitation over both land and ocean, and a significant wettening response over dry tropical ocean regions.
The ocean response was longer-lived than that over land in all models with more than one ensemble member. Monsoon regions dried significantly in agreement with other studies (e.g. Joseph and Zeng 2011, Schneider et al 2009, whilst the ITCZ moved away from the hemisphere with the greatest concentration of aerosols in the multi-model mean in agreement with Haywood et al (2013). We use observed gaugebased data over land; and additionally a combined satellitegauge dataset (GPCP) to examine the observed precipitation  d)) Is for land precipitation, ((e)-(j)) is ocean precipitation, ((a), (b), (e), (f)) global mean, ((c), (d), (g), (h)) wet tropical regions, ((i), (j)) dry tropical regions. GPCP is shown in dark blue, light blue once ENSO is removed; CMIP5 multi-model mean is shown in black and individual runs in grey. Vertical black line denotes timing of eruptions. Dashed lines are 5-95% confidence intervals, dark blue for GPCP, light blue for GPCP with ENSO influence removed and green for GPCP confidence intervals calculated from CMIP5. Yellow shading denotes years in which the multi-model mean response is significant (grey for pre-eruption years). From model results, significant responses are only expected in years 0-3.  . GPCC is shown in dark green, light green once ENSO is removed; CMIP5 multi-model mean is shown in black and individual runs in grey. Vertical black line denotes timing of eruptions. Dashed lines are 5-95% confidence intervals, dark green for GPCC and light green for GPCC with ENSO influence removed. Yellow shading denotes years in which the multi-model mean response is significant. If the 10-90% confidence intervals are used the observed response in MJJASO becomes significant in year one (10% significance level based on one sided test), as does the response with ENSO removed in NDJFMA year 1 (not shown). response over oceans. Results using GPCP were noisy due to the short record length, covering only two eruptions. Nevertheless, over land a decrease in global and wet tropical regions precipitation following eruptions could be seen as expected, and this was significant in NDJFMA. Monsoon regions got drier in their respective warm season, with the exception of SE Asia, and this response was significant against internal variability across the southern hemisphere monsoon regions in austral summer, but was not robust for the northern hemisphere monsoon regions in boreal summer. There was a positive NAO response in winter. Regional mean findings were broadly consistent when alternative observational datasets were used, whilst spatial patterns were very similar. Over ocean, observed results were noisy, but the wet tropical ocean regions got significantly drier, and the dry regions significantly wetter as expected in NDJFMA, although with a larger magnitude than any individual ensemble member. It was not possible to determine whether the long ocean response seen in the models was confirmed by observations due to the noisy nature of GPCP results. Bringing all regions with expected changes in precipitation together, a detection analysis found the modelled influence of volcanism on precipitation was detectable in years 1 + 2 combined and year 1 following eruptions in NDJFMA in all observational datasets, and was marginally detectable in MJJASO in these same years. The response in NDJFMA was significantly underestimated by the models, particularly in the wet tropical regions. Removing the influence of ENSO brought the amplitude of the response in models and observations into better agreement, particularly in NDJFMA.
When the next large volcanic eruption occurs, which is likely to be in the next few decades based on the historical record (Crowley et al 2008), the satellite record will be extremely valuable in further constraining the observed response, particularly over oceans. Tomas Figure 7. Detection of the volcanic signal: regression coefficients obtained by regressing the observed average spatial patterns of precipitation response onto the CMIP5 multi-model mean patterns (circles -see text for details). Coloured bars indicate 5-95% range of regression coefficients for internal climate variability. If a coefficient is greater than the 95th percentile the volcanic influence is detected. Grey crosses are coefficients obtained by regressing single ensemble members onto the mean of the remaining members. Numbers indicate the level at which the response is detectable. Red denotes results based on Z07, green GPCC, and blue GPCP. (c) and (d)) Are for results with the influence of ENSO removed from the observations. Asterisks indicate where less than 5% of the ensemble member-based coefficients are larger than the observed coefficients. ((a) And (c)) are for NDJFMA, ((b) and (d)) MJJASO.