Impact of Volcanic Aerosols on the Hydrology of the Asian Monsoon and Westerlies‐Dominated Subregions: Comparison of Proxy and Multimodel Ensemble Means

Proxy‐model comparisons show large discrepancies in the impact of volcanic aerosols on the hydrology of the Asian monsoon region (AMR). This was mostly imputed to uncertainties arising from the use of a single model in previous studies. Here we compare two groups of CMIP5 multimodel ensemble mean (MMEM) with the tree‐ring‐based reconstruction Monsoon Asia Drought Atlas (MADA PDSI), to examine their reliability in reproducing the hydrological effects of the volcanic eruptions in 1300–1850 CE. Time series plots indicate that the MADA PDSI and the MMEMs agree on the significant drying effect of volcanic perturbation over the monsoon‐dominated subregion, while disparities exist over the westerlies‐dominated subregion. Comparisons of the spatial patterns suggest that the MADA PDSI and the MMEMs show better agreement 1 year after the volcanic eruption than in the eruption year and in subregions where more tree‐ring chronologies are available. The MADA PDSI and the CMIP5 MMEMs agree on the drying effect of volcanic eruptions in western‐East Asia, South Asian summer monsoon, and northern East Asian summer monsoon (EASM) regions. Model results suggest significant wetting effect in southern EASM and western‐South Asia, which agrees with the observed hydrological response to the 1991 Mount Pinatubo eruption. Analysis on model output from the Last Millennium Ensemble project shows similar hydrological responses. These results suggest that the CMIP5 MMEM is able to reproduce the impact of volcanic eruptions on the hydrology of the southern AMR.


Introduction
Large explosive volcanic eruptions inject a large amount of sulfur into the stratosphere. After being converted to sulfate aerosols, they significantly cool the Earth's surface and warm the stratosphere by reflecting the incoming solar radiation and absorbing both solar and longwave radiation (Robock, 2000(Robock, , 2015. Both observation and model results show that the direct surface cooling effects in summer (Kirchner et al., 1999) lead to a significant reduction in summer precipitation, especially in the African and Asian monsoon regions (Iles et al., 2013;Iles & Hegerl, 2014;Trenberth & Dai, 2007;Zambri & Robock, 2016).
The Asian monsoon region (AMR, 8.75°S-56.25°N, 61.25-143.75°E;Cook et al., 2010) covers China and India, two of the world's most populous countries. The AMR has an uneven precipitation distribution due to different dominant winds, with much larger precipitation in the monsoon-dominated subregion (MDSR, southeast) than in the westerlies-dominated subregion (WDSR, northwest). Understanding the hydrological variation of volcanic perturbation in the AMR is both biophysically and socioeconomically important (Dando, 2005). However, there are only a few studies aimed at this region, for example, Anchukaitis et al. (2010), Zhang et al. (2012), Man et al. (2014), Zhuo et al. (2014), and Stevenson et al. (2016Stevenson et al. ( , 2017, which investigated the hydrological effects of historical volcanic eruptions in past centuries. None of these studies distinguished between the different subregions under consideration. This might average out inverse climate response characteristics in the two dominant subregions. Their results show discrepancy, even inverse spatial distribution of the hydrological effects between proxy reconstruction and single model simulation. The discrepancy was mostly imputed to model uncertainties resulting from a biased trust in proxy data. This can place limitations on studying potential mechanisms of volcanic aerosols' hydrological effects with model simulations. Proxy reconstructions also have uncertainties (PAGES 2k-PMIP3 group, 2015). (2017) suggests an equal consideration of the uncertainties and limitations of proxies and models when comparing them.

PAGES Hydro2k Consortium
Recent studies have reported that an ensemble approach leads to a better estimation of climate change as it averaged out unrelated model errors (Flato et al., 2013;Otto-Bliesner et al., 2016) and El Niño-Southern Oscillation (ENSO) effects (Iles et al., 2013;Stevenson et al., 2016). It also enhanced climate prediction skills (Kadow et al., 2015), which are significantly affected by volcanic aerosols (Timmreck et al., 2016). Multimodel ensemble mean (MMEM) of the fifth phase of Coupled Model Intercomparison Project (CMIP5) showed a large improvement in reflecting global temperature and precipitation variations (Flato et al., 2013;Knutti & Sedláček, 2012) as well as monsoon precipitation variations in East Asia (Kusunoki & Arakawa, 2015;Song & Zhou, 2014). Using the MMEM of CMIP5 model output, responses of reduced temperature and summer monsoon rainfall to volcanic eruptions were also clearly detected in historical simulations (Zambri & Robock, 2016) and even in the "last millennium (LM)" experiment of CMIP5 (Zambri et al., 2017). Zambri et al. (2017) focused on the effect of volcanic eruptions on global climate. However, regional scale studies are crucial for developing future water management and coping strategies following a volcanic perturbation. We compare proxy reconstruction and models in the two dominant subregions of the AMR by also taking into consideration uncertainties in proxy and model simulations, as well as the irregular distribution of available tree-ring chronologies. This comprehensive approach aims to answer the following research questions: What are the similarities and discrepancies between proxy reconstruction and model simulations on reflecting volcanic eruptions' hydrological effects in the different subregions of the AMR? Are the CMIP5 MMEMs able to reproduce volcanic eruptions' hydrology effects in the AMR?
In the following, data and methods are detailed in section 2; a comparison of spatiotemporal hydrological patterns is presented in section 3; the sources of uncertainties are discussed in section 4. Finally, we present our conclusions and address the research questions in section 5.

Proxy Data and Covered Subregions
We adopt the Monsoon Asia Drought Atlas (MADA) proxy reconstruction data . It is a reconstruction of June-August (JJA) Palmer Drought Severity Index (PDSI, Palmer, 1965) based on tree ring chronologies and PDSI reconstruction data (Dai et al., 2004), which has annual recordings from 1300 CE to 2005 CE and 2.5°× 2.5°spatial resolution in the AMR. Hereafter, MADA is referred to as the MADA PDSI. As in the PDSI, positive MADA PDSI values represent wet conditions while negative values indicate dry conditions. Drought (flood) occurs when MADA PDSI falls below −0.5 (goes over 0.5). It has been widely used as a reference data set for proxy-model comparisons on volcanic eruptions' hydrological effects in the AMR Stevenson et al., 2016Stevenson et al., , 2017Wegmann et al., 2014;Zhang et al., 2012).
In previous studies, proxy-model comparisons between MADA PDSI and models were conducted over the AMR Stevenson et al., 2016Stevenson et al., , 2017Wegmann et al., 2014). However, no consideration was given to regional differences in dominant climate and data reliability. The AMR is not simply dominated by the monsoon climate; instead different hydrological conditions prevail on both sides of the modern Asian summer monsoon limit (red dashed line in Figure 1). To the northwest are the westerlies-dominated arid areas, the monsoon-dominated humid areas lie to the southeast (Chen et al., 2008;Dando, 2005;Herzschuh, 2006). It includes two monsoon subsystems-the East Asian summer monsoon (EASM) and the South Asian summer monsoon (SASM), which are usually separated at 100°E longitude (Chiang et al., 2017;Herzschuh, 2006). In view of this, we also perform time series analysis over the separated westerlies and monsoon dominated subregions.
with the most tree ring chronologies especially the ones dating back to 1300 CE . Among the monsoon-dominated subregions, SASM has more tree ring chronologies, followed by EASM and SeA. Among westerlies-dominated subregions, several tree ring chronologies are concentrated in the central part of NA with most of them only dating back to 1700 CE ; western-South Asia (w-SA) and CA have fewer tree ring sites.

Model Ensembles and Volcanic Classifications
The "LM" experiment of the CMIP5 was performed only by nine modeling groups (Schmidt et al., 2011). One of the two volcanic forcing data sets-the GRA (Gao et al., 2008) and the CEA (Crowley et al., 2008)-was freely used in the simulations. We separate the models into two groups of MMEMs based on the adopted volcanic forcing data sets. To keep the same number of ensemble members involved in the MMEMs, we adopt six ensemble members of four models in each group, as shown in the green box in Figure 2; more information about the model ensemble members is listed in Table S1 in the supporting information. Only the GISS-E2-R model has three ensemble members, which might predominate the MMEM. Considering this, two sets of MMEMs with four ensemble members (in black in the green box of Figure 2), including only one ensemble member from the GISS-E2-R, are tested. Doubled GRA volcanic forcing was used in the GISS-E2-R model simulations. This exaggerated volcanic forcing might cause excessive climate effects in the GRA-based group of the CMIP5 MMEMs. To verify the model results, we also adopt all the available five ensemble members of the "volcanic only" experiment from the Last Millennium Ensemble (LME; Otto-Bliesner et al., 2016). This project performed large number of LM simulations with the CESM1 (CAM5) model (Hurrell et al., 2013). In the "volcanic only" experiment, the GRA reconstruction (Gao et al., 2008) was adopted as the volcanic forcing data set. Other forcings including solar variability, land use, greenhouse gases, and orbital changes were fixed to the same value as in 850 CE. Following Zhuo et al. (2014), we construct two classifications based on the GRA and the CEA volcanic forcing indices, with the chosen volcanic events that have larger Northern Hemisphere sulfate injection than the 1991 Pinatubo eruption. The two groups of GRA-based and CEA-based Northern Hemisphere volcanic sulfate injection events are named as the GNH and the CNH classifications. As in Zhuo et al. (2014), we adjust the eruption years based on the eruption season. No adjustment is applied for the events with unknown eruption dates as they are assumed to have occurred in spring. We adjust the eruption year to the next year for the events that erupted after August, as their climate impacts are likely to take effect during the next boreal summer. The chosen volcanic events and related aerosol injection magnitude are shown in Figure 2, and the specific values are listed in Table S2. The MADA PDSI covers the period of 1300-2005 CE, while the CMIP5 "LM" experiment covers the period of 850-1849 CE. The overlapping period between 1300 and 1849 CE is chosen as our core study period. Between 1300 and 1849 CE, the GNH classification has 12 volcanic events while the CNH classification has 18 events. Different numbers of classified events may lead to different results between the two classifications. We test this uncertainty using classifications with nine events that are included in both classifications (as shown in red in Figure 2). In order to verify the model results, analyses covering the whole period of 850-1849 CE are also made for the CMIP5 PDSI and the LME PDSI. Twenty volcanic events are used for analysis of the whole period. Except for the same 12 volcanic events in the GNH classification in 1300-1949CE, volcanic eruptions in 901, 933, 1167, 1176, 1195, 1227, 1258, and 1284CE in 850-1299 CE are included in the analysis.

Methods
For a better comparison between proxy reconstruction and models, the CMIP5 "LM" experiment outputs are regridded to the same spatial resolution as the MADA PDSI. Subsequently, using the MATLAB program produced by Jacobi et al. (2013), model precipitation and temperature data, together with latitude and water-holding capacities (Webb et al., 2000), are transferred into the PDSI. Finally, the MMEM of PDSI is calculated. Hereafter, it's referred to as the CMIP5 PDSI. Model ensemble members from LME have the same resolution. These model outputs are directly transferred into the PDSI, and the multimember mean is referred to as the LME PDSI in this study.
Considering that the PDSI combines both temperature and precipitation, we also adopt another widely used hydrological drought index: 12 months of Standardized Precipitation Index (SPI12; Mckee et al., 1993), which is calculated only with model precipitation data. As an index of meteorological drought, SPI12 can reflect water supply in streams and reservoirs. It indicates mild drought once SPI12 falls below zero, while mild wet occurs when SPI12 goes over zero. As in the CMIP5 PDSI, MMEM of SPI12 is calculated and referred to as the CMIP5 SPI12. Since the MADA PDSI only reflects the hydrological condition of the boreal summer season, for consistency, summer JJA means of the CMIP5 PDSI, the CMIP5 SPI12 and the LME PDSI are analyzed in this study.
After pretreatment of the classifications and hydrological data, we conduct Superposed Epoch Analysis (SEA; Haurwitz & Brier, 1981) on hydrological indices (MADA PDSI, CMIP5 PDSI, CMIP5 SPI12, and LME PDSI) for 11 years (−5 to 5) around the eruption year (Year 0) in each classification. To study the significance of the hydrological effects, we conduct Monte Carlo model tests (Adams et al., 2003) based on the null hypothesis that there is no relationship between the volcanoes and the hydrological conditions. Each volcanic event is randomly reassigned a new eruption year in the study period, and then the average values of the hydrological indices are calculated for 11 years. For significance tests of time series analysis, 10,000 times resampling are made on regionally averaged hydrological indices. For spatial analyses, 1,000 times resampling are made on each grid. This builds a random distribution, against which our SEA results are considered to be statistically significant at the 95% (99%) confidence level when they exceed the 95% (99%) range of the Monte Carlo results. To quantify the same drought and wet areas between proxy and the MMEMs, we count the number of grid cells that have the same sign between the MADA PDSI and the CMIP5 PDSI/SPI12 and calculate their percentage in each subregion. Figure 3 shows the SEA results of the MADA PDSI, the CMIP5 PDSI, and the CMIP5 SPI12 over the Asian monsoon region for the GNH and the CNH volcanic classifications. As shown in Figure 3a, the MADA PDSI decreases 1 year after the eruption (Year 1), and a significant drying effect emerges 2 and 3 years after the eruption (Years 2 and 3). The CMIP5 PDSI decreases promptly and sharply in the eruption year (Year 0). The significant drying effects last for 3 years and gradually return to normal conditions in Year 4. Similarly, the CMIP5 SPI12 decreases rapidly in Years 0 and 1; after the strongest drying effect in Year 1, it gradually recovers in Year 2 and returns to normal condition in Year 3. This indicates an agreement between the MADA PDSI and the CMIP5 PDSI/SPI12 on the drying effects of the volcanic eruptions, although with 1 year time lag in the MADA PDSI compared to the CMIP5 PDSI/SPI12; also, the magnitude shown in the CMIP5 PDSI/SPI12 is much larger than that in the MADA PDSI. This is probably due to exaggerated doubled GRA forcing used in the GISS-E2-R model simulations. Figure 3b shows hydrological responses to the CNH volcanic eruptions. The trend of response is similar to that in the GNH classification, that is, the MADA PDSI increases before the eruption and decreases in Years 1 and 2; the CMIP5 PDSI and the CMIP5 SPI12 decrease promptly in Year 0 and reach the lowest value in Year 1 and gradually recover from Year 2. Compared to the significant results (even at the 99% confidence level) in the GNH classification, the results are less significant in the CNH classification, but the magnitude of change between the MADA PDSI and the CMIP5 PDSI/SPI12 is closer. The different scales between MMEMs of two classifications might also result from different number of superposed volcanic events. But even with the same number of superposed volcanic events (the same nine events in both volcanic forcing indices, shown in red in Figure 2), the differences between the two classifications persist. Crowley and Unterman (2013) suggested that the volcanic forcing in the GRA index is overestimated. When reconstructing the CEA index, they used a scaling of two third to calculate the forcing of the explosive eruptions, which are larger than 1991 Pinatubo eruption. Volcanic events included in the CNH classification are affected by this scaling process, which result in minor hydrological responses.

Temporal Hydrological Responses to Volcanic Classifications
The exaggerated doubled GRA forcing used in the GISS-E2-R model simulations also caused excessive climate response in the GNH classification. As illustrated in Figure 4, temporal SEA results of the LME PDSI over two periods show significant drying effects in Year 0, which last up to Year 1. CMIP5 PDSI shows excessive drying effects over the whole period, which is similar to that shown in Figure 3a over the core study period. This confirms the findings shown in Zhuo et al. (2014) that larger magnitude volcanic aerosols lead to larger drying effect. The significant drying effects 2 to 3 years after the volcanic eruptions are consistent with 10.1029/2020JD032831

Journal of Geophysical Research: Atmospheres
previous studies Liu et al., 2016;Man et al., 2014;Zhuo et al., 2014) and is prominent even against the general background of a significant reduction in global precipitation (Iles & Hegerl, 2014).
Temporal SEA analysis over the whole region confirms different climate conditions in the westerlies-and monsoon-dominated subregions. Additionally, temporal SEA results over the separated westerlies and monsoon dominated subregions are presented. Figure 5a shows different hydrological responses in the WDSR, as  This might indicate that the hydrological response to volcanic perturbations in the WDSR is insensitive to volcanic forcing. However, the CMIP5 PDSI and the CMIP5 SPI12 in the GNH classification show highly significant drying effects in Years 0 to 3 at the 99% confidence level. In the CNH classification, the CMIP5 SPI12 shows significant drying effects in Years 1 and 2, but the CMIP5 PDSI does not indicate drying effect, instead only significant wetting variations in Year −1. This indicates a large difference between the MADA PDSI and the CMIP5 PDSI/SPI12 in the WDSR. Considering the exaggerated volcanic forcing used in the GNH classification, this might suggest that the wetting or drying effect in this insensitive area depends largely on the magnitude of the injected volcanic aerosols. In the MDSR (Figure 5b), the MADA PDSI and the CMIP5 PDSI/SPI12 in both classifications agree on the drying effects in Years 0 and 1 and recovery from Year 2 onward. The MADA PDSI decreases most in Year 1, whereas the CMIP5 PDSI/SPI12 decreases most in Year 0. This probably indicates a time lag effect in proxy data. Figure 6 shows the temporal SEA results of the LME PDSI over both periods and the CMIP5 PDSI for the whole period over the separated westerlies and monsoon dominated subregions. In the WDSR, the

Journal of Geophysical Research: Atmospheres
CMIP5 PDSI indicates significant drying effects for the whole period (850-1849 CE). The LME PDSI increases in Year 1 and decreases from Years 2 to 5 over both periods. This is similar to the response of the MADA PDSI in the GNH classification (Figure 5a), but both results did not pass the significance tests even at the 95% confidence level. In the MDSR, all the model results suggest consistent drying effects in the Years 0 and 1 and then gradual recovery in Year 2.

Spatial Patterns of Hydrological Response
Considering uncertainties in spatial responses arising from the estimated aerosol magnitude in volcanic forcing reconstructions, following discussions focus on the horizontal distribution of the hydrological tendencies. To quantify the similarity in drought and wet areas between proxy and model in Figure 7, we show percentages of grid cells that have same sign between the MADA PDSI and the CMIP5 PDSI/SPI12 in different subregions in Years 0 (in magenta) and 1 (in red). When separated into two dominant subregions, it is hard to find out consistent tendency in these variations, except that the percentage increases from Years 0 to 1 in the WDSR, decreases in the MDSR, and different ensemble members show larger difference in the WDSR than in the MDSR. When the dominant regions are separated into seven subregions, taking the spatial coverage of tree-ring chronologies into consideration, the percentages show large differences in different subregions. Four subfigures indicate largest similarity between the MADA Figure 7. Histogram of percentages of grid cells that have same sign between the MADA PDSI and the CMIP5 PDSI/SPI12 in the GNH classification (a and b) and the CNH classification (c and d). Columns indicate the percentages in the westerlies-dominated subregion (WDSR) and monsoon-dominated subregion (MDSR), as well as in the seven subregions in Years 0 (in magenta) and 1 (in red). Different markings indicate the percentages between the MADA PDSI and the PDSI/SPI12 of different single ensemble members.

Journal of Geophysical Research: Atmospheres
PDSI and the CMIP5 PDSI/SPI12 in the w-EA, where the most tree-ring chronologies are available, in both Years 0 and 1. They also show fewest differences between different ensemble members. However, single models and the MADA PDSI have large uncertainties. The consistency between different groups and ensemble members improves the reliability of both proxy and models in reproducing the hydrological effects of volcanic eruptions. Better agreements are then shown in SASM and EASM, with fewer differences between different ensemble members than in w-SA and SeA. These results suggest that proxy and models agree better in the subregions with more tree-ring chronologies confirming the importance of tree-ring chronologies in providing accurate proxy reconstruction data. The percentages are larger in Year 1 than in Year 0, except for NA and CA with fewest tree-ring chronologies. This is consistent with the temporal SEA results and spatially shows that the MADA PDSI and the CMIP5 MMEMs agree better in monsoon-dominated subregions in Year 1.
To investigate the spatial hydrological variation, we show the spatial patterns of the superposed hydrological responses to the GNH classification in Figure 8.  (Wu et al., 2015) and the dating uncertainty of volcanic eruptions might also contribute to the difference. General agreements are shown in the w-EA, which has a dense coverage of tree-ring chronologies. In the westerlies-dominated subregions with scarce tree-ring chronologies, the MADA PDSI shows wetting to drying transition, the CMIP5 PDSI/SPI12 show continuous drying effects in the CA and the NA, but wetting effects in the w-SA. This is consistent with the temporal SEA results shown in Figure 5, indicating that the MADA PDSI and the CMIP5 PDSI/SPI12 agree on the tendency of the hydrological response to volcanic perturbation in the MDSR, while discrepancies exist in the WDSR.
In the CNH classification (Figure 9), the MADA PDSI ( Figure 9a) and the CMIP5 PDSI (Figure 9b) indicate weaker effect of volcanic perturbations. It shows similar hydrological patterns as that in the GNH classification (Figures 8a and 8b), except that the CMIP5 PDSI shows limited response in the CA and the NA (Figure 9b). The drying effects shown in the GNH classification (Figure 8b) might be a response to the exaggerated volcanic forcing used in GISS-E2-R model simulations. The CMIP5 SPI12 (Figure 9c) indicates an even weaker effect, but the obvious drought areas agree well with those in the CMIP5 PDSI patterns. Better agreement between the MADA PDSI and the CMIP5 PDSI/SPI12 occurs in the subregions with more tree-ring chronologies. Highly significant results of the CMIP5 PDSI/SPI12 in the GNH and the CNH classifications indicate the consistency of MMEMs in reproducing the effect of volcanic aerosols on hydrology in the southern AMR. In the CA and the NA, discrepancies between the MADA PDSI and the CMIP5 PDSI/SPI12 do not allow drawing of definite conclusions.
To verify model results, we show spatial patterns of the LME PDSI over both periods (1300-1849 CE and 850-1849 CE) and the CMIP5 PDSI over the whole period (850-1849 CE) in Figure 10. The CMIP5 PDSI shows similar patterns even when the period is extended to the whole 1,000 years ( Figure 10c). Similar patterns are also shown in the LME PDSI over both periods, especially in the southern Asian monsoon region. The drought and wet areas in different model results are not same. However, it is difficult to get a complete 10.1029/2020JD032831

Journal of Geophysical Research: Atmospheres
match with different model resolutions. This indicates that in the study periods, the number of the superposed events and the aerosol magnitude does not much affect the spatial pattern of hydrological changes. These similar patterns support the reliability of models in reproducing the hydrological effects of volcanic eruptions in the southern Asian monsoon region. The LME PDSI also suggests a slight drying effect in the NA from Years 0 to 2 over both periods. These PDSI patterns might suggest that drying effects can emerge in the NA and the CA with a strong enough volcanic forcing.

10.1029/2020JD032831
Journal of Geophysical Research: Atmospheres eruptions. In the southern Asian monsoon region, spatial patterns of the MMEMs in Years 0 and 1 agree well with precipitation anomaly pattern after Krakatau and Pinatubo eruptions shown in Zambri and Robock (2016). The identified wet areas in the EASM are close to that in Gao and Gao (2018), which showed an increased precipitation over the Yangtze-Huaihe River valley using Feng et al. (2013) precipitation reconstruction. The patterns are also consistent with the observed precipitation and the PDSI variations shown in Trenberth and Dai (2007), with a drying effect in the SASM, the SeA and the northern EASM, and wetting effect in the w-SA and the southern EASM after the Mount Pinatubo eruption. In the northern Asian monsoon region, except for the PDSI, which suggests drying effects (Trenberth & Dai, 2007), limited changes in precipitation (Trenberth & Dai, 2007;Zambri & Robock, 2016) and runoff (Trenberth & Dai, 2007) are seen.
These results indicate the reliability of the MMEMs in reflecting the spatiotemporal patterns of hydrological response to volcanic perturbations in the Asian monsoon region, except for Central Asia and North Asia; definite conclusions in these two subregions cannot be drawn as the CMIP5 PDSI and the CMIP5 SPI12 in the CNH classification display no impact, and there are limited available modern observations in these subregions to validate the results.

Discussion on Sources of Uncertainty
Our results indicate large discrepancies between the MADA PDSI and models in the westerlies-dominated subregions with fewer available tree-ring chronologies and a better agreement 1 year after the eruption instead of in the eruption year. These discrepancies indicate the uncertainty in the results derived from the data sources and analysis processes. As suggested by the PAGES Hydro2k Consortium (2017), we treat proxy reconstruction and model data in the same manner and discuss uncertainties and limitations of both the MADA PDSI and the CMIP5 PDSI/SPI12. From the temporal SEA results of the Asian monsoon region (Figure 3), we can see that the CMIP5 PDSI/SPI12 agrees with the MADA PDSI on the drying effects of explosive volcanic eruptions. The CMIP5 PDSI shows stronger effect than the MADA PDSI in the GNH classification. Stronger effects are also shown in the GNH classification as compared to the CNH classification. This is caused by both the exaggerated volcanic forcing used in the GISS-E2-R model ensemble members and the reduced amplitude of the forcing in the CEA reconstruction (Crowley & Unterman, 2013). This can be verified by the results of the LME PDSI. Besides, faster responses are shown in the CMIP5 PDSI/SPI12 than in the MADA PDSI in both classifications. This reflects the time lag effect in the tree-ring-based ecological response compared to the meteorological response in the model simulation (Wu et al., 2015).
Volcanic years identified in the volcanic forcing indices deviate from reality. Superposed volcanic classification averages out the effect of single event, but the dating uncertainty in volcanic events can cause large uncertainty on the hydrological effects reproduced by the MADA PDSI. Uncertainty in the eruption month of the volcanic forcing indices also causes uncertainty in defining the eruption year. This might explain the abnormal wetting effect in Year 0 shown by the MADA PDSI (Figure 3), which was also identified in Anchukaitis et al. (2010) after different superposed volcanic events. To investigate these uncertainties, we test several different classifications. Volcanic years and the number of included events in different classifications are listed in Table S3. We show response of the MADA PDSI to these different classifications in Figure 11. As in the GNH and the CNH classifications, the MADA PDSI also suggests wetting effect in Year 0 in the GCNH classification. However, the MADA PDSI starts to decrease in Year 0 and drops to the lowest value in Years 2 and 3 in the SNH and the A07 classifications, respectively. The SNH classification is based on the latest start-of-the-art volcanic forcing reconstruction with largely improved dating accuracy (Sigl et al., 2015), while the A07 classification includes only five well-known explosive eruptions in the past centuries. These two classifications have minimum dating uncertainty among the volcanic classifications used in this study. This indicates that the dating uncertainty largely affects the climate response especially in Year 0. The wetting effects shown by the MADA PDSI probably result from dating uncertainty of the volcanic events.
The temporal SEA results of two separated subregions ( Figure 5) suggest an agreement between the MADA PDSI and the CMIP5 PDSI/SPI12 in the MDSR while large discrepancies exist in the WDSR. Quantification of the grid cells with same sign between the MADA PDSI and the CMIP5 PDSI/SPI12 indicates better agreement in subregions with more tree-ring chronologies and in the second summer after the volcanic perturbations. This might partly explain the spatial proxy-model discrepancies suggested by previous studies Stevenson et al., 2016Stevenson et al., , 2017Zhang et al., 2012), because comparisons were only made in the first summer after the eruptions. Comparison of the spatial patterns of the hydrological effects shows large discrepancies in the westerlies-dominated subregions with limited tree-ring chronologies revealing the limitation of the MADA PDSI. The drying tendencies in the NA and the CA reflected by the CMIP5 PDSI/SPI12 in the GNH classification might be realistic, but it could be that the misleading patterns result from the exaggerated volcanic forcing used in the GISS-E2-R model ensemble members. This exaggerated forcing also contributes to the faster and longer drying effects of volcanic perturbation in the monsoon-dominated subregions.
Spatial SEA results ( Figure 8) indicate significant wetting effect in the NA by the MADA PDSI, which is opposite to the drying effects shown by the CMIP5 PDSI and the LME PDSI in the GNH classification. Similar discrepancies were also presented in Liu et al. (2016). This may indicate data uncertainties in the MADA PDSI, especially in the westerlies-dominated subregions where few tree-ring chronologies extend up to 1300 CE . However, we would like to point out that the models also suggest wetting effects in the western areas. They show variations in different groups of model ensemble means with different volcanic forcing magnitudes. Thus, the discrepancies can be also caused by uncertain aerosol magnitudes and the consequent uncertain effects shown in the models. The differences in resolution between proxy reconstructions and models also introduce uncertainties.
A limited number of ensemble members might also introduce uncertainty in the model results.
Especially, three ensemble members of the GISS-E2-R model might have a predominant effect on the MMEMs. However, when testing the MMEMs with only four members (members in black in Figure 2), which include only one member of the GISS-E2-R model, temporal and spatial patterns remain largely unchanged. Deviation of the model-based analysis between two classifications can come from the number of classified events based on the volcanic forcing indices. When testing the classifications with the same nine events in both indices (marked in red in Figure 2), temporal and spatial patterns remain largely constant.
The internal variability of the climate system often brings uncertainty in detecting the hydrological effects of volcanic eruptions, especially the barely constrained effect of the concurrent ENSO events (Adams et al., 2003;Khodri et al., 2017;Li et al., 2013;Stevenson et al., 2016Stevenson et al., , 2017. The effect of the eruption season on the circulation and the ENSO can introduce additional uncertainties (Stevenson . All these factors might contribute to the proxy-model discrepancies, especially in the initial phase and the phase-out period of the hydrological effects. Following the method in Iles et al. (2013), we test this uncertainty by repeating the SEA analysis after regressing out the effect of ENSO. Consistent with Iles et al. (2013) and Iles and Hegerl (2014), it only results in a lower amplitude response, but the temporal and spatial patterns remain largely unchanged. In addition, previous researches show that volcanic eruptions can affect the hydrological condition by affecting the temporal evolution of the ENSO but with widely contradictory findings (Adams et al., 2003;Li et al., 2013;Liu et al., 2018;Stevenson et al., 2016;Sun et al., 2018;Wang et al., 2017). This is an additional source contributing to proxy-model discrepancies. Future improvement of volcanic forcing reconstructions, model simulations, proxy reconstructions, and observations will lead to a better understanding and reconciling of the proxy-model discrepancies.

Conclusions
Previous studies show large discrepancies between proxy and model data on hydrological response patterns to volcanic aerosols in the Asian monsoon region. In this study, we use tree-ring based proxy data, that is, MADA PDSI and a number of model ensemble members from the CMIP5 and the LME, to compare their spatiotemporal hydrological response to two classified volcanic events in 1300-1850 CE in the subregions of monsoonal Asia.
Our temporal SEA results show that the MADA PDSI and models agree on the significant drying effects of volcanic aerosols in the MDSR, while disagreement exists in the WDSR. Spatial comparisons indicate better agreement in subregions with more available tree-ring chronologies. Especially in the w-EA, where most available tree-ring chronologies dating back to 1300 or even earlier are available, the MADA PDSI agrees with models on the significant drying effects of volcanic aerosols. In monsoon-dominated subregions, the MADA PDSI and models show similar drying to wetting variations after the volcanic perturbations but with faster and prolonged drying effects in the latter, which might result from an overestimated aerosol magnitude in the volcanic forcing index and the time lag effect of tree-ring-based proxy reconstruction data. The effect of uncertain eruption season on the atmospheric circulation and the definition of the eruption year might also contribute to their differences. Because of these uncertainties, MADA PDSI and models show better consistency in Year 1, with significant drying effects in the northern EASM, the SASM, and the SeA, and opposite wetting effects in the southern EASM. In westerlies-dominated subregions, with fewer tree-ring chronologies MADA PDSI and models show larger discrepancies. Since two groups of the CMIP5 MMEMs and the LME PDSI all show similar patterns, and with verification from previous studies, we propose the reliability of the CMIP5 MMEMs in reproducing the wetting effects in the w-SA. In the CA and the NA subregions, MADA PDSI shows significant wetting effects. The CMIP5 MMEMs in the GNH show significant drying effect; the LME PDSI shows some drying effect only in the NA, whereas the CMIP5 MMEMs in the CNH show limited response. Considering the lack of cross verification, we do not draw definite conclusion in these two subregions.
Through spatiotemporal comparison, we examine the reliability of the MADA PDSI and the CMIP5 MMEMs in reproducing the pattern of hydrological responses to volcanic perturbations. Our results suggest higher reliability of the MADA PDSI in subregions with more tree-ring chronologies. Comparisons between proxy, observation and models indicate that the CMIP5 MMEMs are reliable in reproducing the hydrological effect of volcanic aerosols in the southern Asian monsoon region. Further analysis on the CMIP6 and improved proxy reconstruction data will contribute to verification of these results.
This study discusses the long-standing problem of discrepancy between proxy-model data. We treat the uncertainties and limitations of both proxy and models equally, thus contributing to better interpretation of results. This approach sheds new light on the reliability of both proxy data and the CMIP5 model simulations in reproducing the hydrological effects of historical volcanic eruptions in the Asian subregions. Based on our results, suggesting the reliability of the MMEM, additional studies aimed at exploring the mechanisms of hydrological changes can be pursued in future. This is important for better evaluating the effects of both historical and future volcanic eruptions and assessing the feasibility of future stratospheric aerosol injection engineering (Simpson et al., 2019).