A Comparison of Model Simulations of Asian Mega-droughts during the past Millennium with Proxy Reconstructions

Two PMIP3/CMIP5 climate model ensemble simulations of the past millennium have been analysed to identify the occurrence of Asian mega-droughts. The Palmer drought severity index (PDSI) is used as the key metric for the data comparison of hydro-climatological conditions. The model results are compared with the proxy data of the Monsoon Asia Drought Atlas (MADA). Our study shows that global circulation models (GCMs) are capable of capturing the majority of historically recorded Asian monsoon failures at the right time and with a comparable spatial distribution. The simulations indicate that El Niño-like events lead, in most cases, to these droughts. Both model simulations and proxy reconstructions point to fewer monsoon failures during the Little Ice Age. The results suggest an influential impact of volcanic forcing on the atmosphere–ocean interactions throughout the past millennium. During historic mega-droughts of the past millennium, the monsoon convection tends to assume a preferred regime described as a " break " event in Asian monsoon. This particular regime is coincident with a notable weakening in the Pacific trade winds and So-mali Jet.


Introduction
Mega-droughts are natural phenomena that last for years to decades and emerge at irregular intervals (Cook et al., 2010).They may play a pivotal role in societal changes in Asia: for example, the collapse of the Yuan Dynasty (Zhang et al., 2008), the Ming Dynasty of China (Shen et al., 2007) and the Khmer Empire of Cambodia (Buckley et al., 2010) are attributed to mega-droughts.Mega-droughts have been linked to persistent patterns of sea surface temperature (SST) anomalies in the Indian and Pacific oceans (Ropelewski and Halpert, 1987;Meehl and Hu, 2006;Wahl and Morrill, 2010).Prolonged droughts during the 21st century are expected from global warming (Rind et al., 1990;Wang, 2005;Burke and Brown, 2008;Sheffield and Wood, 2008;Seager et al., 2007;Dai, 2011).In Asia, droughts are always connected to the monsoon circulation.A decline in the rainfall amount or a change in precipitation pattern may lead to droughts in the Asian monsoon region (Shaw and Nguyen, 2011).The detection and analysis of the natural variability of extreme droughts in Asia is complicated by the short period of available observations and the proxy networks with their coarse spatial distribution.
The Monsoon Asia Drought Atlas (MADA) (Cook et al., 2010) provides a proxy-based gridded spatiotemporal reconstruction of the Palmer drought severity index (PDSI) (Palmer, 1965) for the past millennium.The PDSI presents the deviation from a normal hydrological balance level based on a simple water balance model.It quantifies the severity of drought in time and space based on precipitation and temperature values (Wells et al., 2004).Several studies have been devoted to the analysis of the models' capability to simulate the mega-droughts and the processes behind them (Li et al., 2013;Anchukaitis et al., 2010;Dai, 2013).Many simulations capture past mega-droughts at the wrong times (Schiermeier, 2013).Atmospheric simulations using models of intermediate complexity have shown that a warmer tropical Pacific leads to more anomalous south-westerly moisture transport into central Asia (Mariotti, 2007).However, these simulations were driven by the prescribed observed SST.
Global circulation models (GCMs) are able to capture prolonged droughts, El Niño-Southern Oscillation (ENSO) and the global mean aridity trend in recent climate (Dai, 2013).
Published by Copernicus Publications on behalf of the European Geosciences Union.
According to Krishnamurthy and Shukla (2007), inter-annual changes in monsoon precipitation consists of a seasonal mean component and an unpredictable inter-seasonal variability.This leads to the question of to what extent the simulated mega-droughts are linked to external forcing (e.g.greenhouse gases, volcanism) or may be the result of internal variability.
This study evaluates the ability of fully coupled oceanatmosphere climate models in identifying the potential long-term drivers of the moisture anomalies in central and monsoon-dominated Asia.Several coupled atmosphereocean model simulations within the PMIP3/CMIP5 project (Paleoclimate Modelling Intercomparison Project Phase III/Coupled Model Intercomparison Project Phase 5) are used to investigate the atmosphere-ocean interactions for the past millennium.After introducing the model experiments, the proxy data and the statistical methods (Sect.2), the models' capability of simulating mega-droughts is described (Sect.3).Thereafter, we investigate the relationship between the surface temperature and monsoon failures (Sect.4).Additionally, we investigate the long-term evolution of droughts over the past millennium and their relation to global mean temperatures and monsoon convection and circulation, followed by a discussion and conclusions (Sect.5).

Data and methods
The "millennium" simulations of two coupled atmosphereocean general circulation models (AOGCMs) from PMIP3 and CMIP5 projects are analysed: ECHAM5/MPIOM (five ensemble members) from Max Planck Institute for Meteorology (MPI-M) (Jungclaus et al., 2010) and GISS-E2-R (eight ensemble members) from NASA Goddard Institute for Space Studies (NASA GISS) (Schmidt et al., 2014).The forcings for these simulations include realistic presentations of solar variability, volcanism and greenhouse gas concentrations (Schmidt et al., 2012).These simulations are used for the Intergovernmental Panel on Climate Change Fifth Assessment Report (IPCC, 2013).They were downloaded from https:// pmip3.lsce.ipsl.fr/.Two of the GISS-E2-R experiments contain no volcanic forcing (r1i1p123 and r1i1p126).In three of the experiments (r1i1p122, r1i1p125 and r1i1p128), volcanic forcing was about a factor of 2 "larger than intended" (http: //data.giss.nasa.gov/modelE/ar5/).Our empirical orthogonal function (EOF) analysis reveals that all these experiments show a completely different response pattern in PDSI compared to MADA (not shown).The forcings in the remaining experiments (r1i1p121, r1i1p124 and r1i1p127) are very similar to the ECHAM5/MPIOM simulations.However, the reconstructed land use/land cover forcing in r1i1p127 experiment is notably larger than the one in r1i1p121 and r1i1p124 (Schmidt et al., 2011).Therefore, we used the r1i1p121 and r1i1p124 experiments in our analysis; these are called GISS-E2-R hereafter.
To evaluate the model results in detecting the historical mega-droughts, we compared the model results with MADA.MADA presents a gridded summer monsoon metric based on the proxy data set of 327 annual tree-ring chronologies (Cook et al., 2010).The number of tree-ring chronology networks in MADA rises strongly in the post-1700 period.To investigate the link between droughts, ocean-atmosphere dynamics and global warming, we used Northern Hemisphere temperature reconstructions (Shi et al., 2013).Three different Niño indices are calculated from simulated area-averaged SST anomalies over the Niño 1+2 region (90-80 • W and 10 • S-0 • ), Niño 3.4 region (170-120 • W and 5 • S-5 • N) and Niño 4 region (160 • E-150 • W and 5 • S-5 • N).To evaluate the ensemble of model experiments, we applied arithmetic averaging to each of the two model experiments (PDSI is calculated for each member separately prior to the averaging).
In order to quantify the severity of droughts in time and space across different climates, we derived the PDSI (Palmer, 1965) from simulated monthly summer (JJA) precipitation and temperature.According to Cook et al. (2010), the first EOF (EOF1) pattern of reconstructed PDSI is associated with a persistent weak monsoon in Asia.Therefore, the first principal component gives us information about "active" and "break" monsoonal phases.In analogy, the first EOF of simulated PDSI from model's ensemble average conveys us similar information to the EOF1 of MADA.
Maximum covariance analysis (MCA) is used to identify the coupled patterns in the climate data.Dai (2013) utilized MCA to explore the relationship between simulated global SSTs and the PDSI.MCA identifies the most important modes of climate data in which the variability of the two fields is coupled (Bretherton et al., 1992).The advantage of MCA compared with coupled EOF analysis is that it focuses on those modes of anomalies which are highly coupled (Bretherton et al., 1992).The first leading MCA modes of simulated SSTA (sea surface temperature anomaly) represent the ENSO-like patterns.Yearly averaged monthly monsoon (JJA) SSTs were used for MCA.A period from 1300 to 1860 was chosen to exclude the impact of anthropogenic climate effects.The GCM data were re-gridded by means of bilinear interpolation to MADA's grid (ca 2.5 • × 2.5 • ).The spatial domain of MADA was considered in our analysis of PDSI (61.25-143.75 • E and 8.75 • S-56.25 • N).The SST data of the GISS-E2-R model are bilinearly remapped to MPIOM's grid (horizontal resolution of GR3.0, ca. 3 • × 3 • ).A 31-year filter is used to smooth the time series.

Simulation of Asian mega-droughts
The first leading EOF pattern of simulated PDSI presents a monsoon failure pattern over India and south-eastern Asia.Fig. 1 presents the first EOF (EOF1) pattern from reconstructed PDSI (Cook et al., 2010;Dai, 2011;Li et al., 2013).This pattern is comparable to the principal com-  ponent pattern of modelled PDSI.The pattern correlation between EOF1 of reconstructed and averaged simulated PDSI of ECHAM5/MPIOM (GISS-E2-R) is 0.78 (0.57) (Fig. 1).The explained variances for EOF1 of MADA, ECHAM5/MPIOM and GISS-E2-R are 19.52,24 and 12.28 %, respectively.The EOF1 pattern for each of the model members is very similar to the one from MADA (not shown).However, the time expansions show a large variability and the ensemble averaging improves the agreement between model and proxy (Polanski et al., 2014;Supplement).According to Kalnay et al. (1996), the ensemble average is more accurate than a single deterministic climate simulation.Lamber and Boer ( 2001) concluded that the climatological fields from ensemble average have more agreement with the observations than the fields generated by any single member.
Those periods during which both the model and the reconstruction have the same sign in the selected principal components (PCs) are described as "active" or "break" phases of the monsoon (green/brown bars in Fig. 2).Monsoon "break" phases are more common after the late 17th century in GISS-E2-R simulations.During the Little Ice Age (LIA), active phases of the monsoon were more frequent.This feature is captured by both the model experiments and the proxy reconstruction.Five of the large-scale mega-droughts recorded in MADA (see Cook et al., 2010) are selected for comparisons: (i) mid-14th century drought associated with collapse of the Khmer Empire (1351-1368) in Cambodia (Cook et al., 2010 andBuckley et al., 2010), (ii) the late 16th century severe drought (1560-1587), (iii) the end of the 17th century drought (1682-1699), (iv) the "Strange Parallels" drought (1756-1768) and (v) the East India drought (1790-1796).The simulated PCA1 time series of the GISS-E2-R model detects the last four monsoon failures, but disagrees with MADA during the mid-14th century (Fig. 2a).The ECHAM5/MPIOM model simulation does not simulate the timing of monsoon failure during the Strange Parallels drought.
The EOF analysis is a linear method for reducing the dimensionality of complex data (Hannachi and Turner, 2013).Therefore, it can not capture all possible drought patterns in the data set.By applying EOF analysis, only those droughts which have a similar dipole pattern to that shown in Fig. 1 are detected, and other possible drought patterns can not be identified.
In the following we investigate the historical megadroughts individually (Fig. 3).Based on the reconstructions, the mid-14th century drought (Fig. 3a) exhibits three drought regions over the Asia: (i) most of India; (ii) Vietnam, Thailand, Laos, southern China; and Myanmar and (iii) Mongolia, north-eastern China and Lake Baikal.Figure 3a shows that the GISS-E2-R model captures these three key regions of the mid-14th century drought.The ECHAM5/MPIOM model, however, simulates a wet spell over most of India during this period.
The late 16th century drought (Fig. 3.b) presents two major drought regions in MADA: (i) over India and Bangladesh and (ii) Kazakhstan and northern Mongolia.Both models and MADA indicate the dipole pattern between India and arid central Asia region for this prolonged drought.The ECHAM5/MPIOM simulation, however, disagrees with MADA for the Lake Baikal region.The GISS-E2-R model produces drought patterns similar to the patterns in MADA.
The final stage of the 17th century drought displays a similar pattern to the mid-13th century drought in MADA (Fig. 3c).The GISS-E2-R model captures the two key regions of (i) northern India and Pakistan and (ii) north-eastern China and Mongolia.In the ECHAM5/MPIOM model simulation, the first pattern is shifted to the north.The Strange Parallels drought, the most persistent mega-drought during the last 700 years, had a dominant pattern over India and south-eastern Asia that lasted for more than a decade (1756-1768) (Buckley et al., 2007;Lieberman, 2009;Sano et al., 2009;Cook et al., 2010).
The broad drought patterns captured by reconstructions and simulations are comparable during the Strange Parallels drought (Fig. 3d).In ECHAM5/MPIOM model simulations, all the three key regions of dry conditions (India; southeastern Asia; Siberian plains) are captured.The GISS-E2-R model does not simulate the drought pattern over Kazakhstan and the north of the Lake Balkhash.
The East India drought, an El Niño-induced event, occurred in the late 18th century and led to severe famine in India (Cook et al., 2010).The reconstruction exhibits several dry regions for this period: north-eastern India and the western and northern Himalayas (Fig. 3e).Cook et al. (2010) suggested that the Indian monsoon was not greatly weakened during this period.The model simulations indicate drier conditions over northern and north-eastern India and the northern Himalayas, which are sources of the water supply for the Indian subcontinent.

ENSO
The time expansion of the EOF1 mode of PDSI correlates significantly with central Pacific SST anomalies; such a correlation characterizes ENSO (Cook et al., 2010;Li et al., 2013).To investigate the influence of the ocean on the Asian monsoon, an MCA of simulated PDSI and SSTs has been performed.By considering only the leading MCA modes, we excluded several other natural variations.We focus on the first mode of MCA (MCA1) of model simulations.MCA1 of the ECHAM5/MPIOM model accounts for 90 % of explained squared fractional covariance (SFC), and MCA2 of GISS-E2-R for 40 %.This pattern shows a warm anomalous (El Niño-like) SSTA over the central tropical Pacific and Indian oceans (Figs. 4 and 5).Krishnan Kumar et al. (2006) concluded that the warmer central equatorial ("westward-shifted") Pacific Ocean initiates more extreme droughts over India.The time expansions of this mode for PDSI and SST are significantly (p value < 0.01) correlated.The temporal coefficient of the MCA1 pattern of SST from the ECHAM5/MPIOM model (blue line in Fig. 4a) is highly correlated (corr.coeff.= 0.93) with the SST anomalies over the Niño 4 region (Fig. 4a).The red line in Fig. 5a shows the SST anomalies over the Niño 4 region for the GISS-E2-R model.The correlation coefficient between MCA1 of SST and the Niño 4 index for GISS-E2-R is 0.84.
Table 1 shows the correlation coefficients between different Niño indices, MCA1 of PDSI and SSTA, and reconstructed global mean temperatures.The behaviour of individual members is shown in the Supplement.During the late 17th century, MCA time series of simulated SSTs indicate negative patterns in the central Pacific which resemble a La Niña event during this drought.ENSO may have been an important trigger of yet another recorded mega-drought.During the late 18th century, both models simulate an El Niño event.
The great El Niño of this period was the main cause of worldwide societal and economical disturbances (Grove, 2007).
The time series in Figs. 4 and 5 suggest an influence of volcanic external forcing (magenta lines in Figs. 4 and 5), especially in the mid 15th, late 17th and early 19th century.There are clear minima in the time series coincident with major volcanic eruptions.

Monsoon convection and circulation
Turner and Hannachi (2010) defined two preferred regimes in monsoon convection in the ERA-40 (Uppala et al., 2005) reanalysis.They suggested that these regime behaviours may   first EOF of summer (JJA) OLR anomalies.The explained variances for the ECHAM5/MPIOM and GISS-E2-R models are 56 and 13.68 %, respectively.The probability distribution function (PDF) estimate of OLR shows a "shoulder" in the positive PC1 values (Fig. 6b).The skewness of the data for ensemble mean is high (0.42).The distance between the two peaks is larger for any individual members (Fig. 6c-g).The ensemble average produces a better estimate of the mean state; thus its PDF tends towards a normal distribution (Fig. 6b).
Weighting functions as in the study of Turner and Hannachi (2010) are applied for each of these regimes and the composites of 850 hPa wind and OLR are calculated.In contrast with the study of Turner and Hannachi (2010), which used the ERA-40 reanalysis, regime 1 happened more frequently than regime 2 during 1300-1850.
Figure 8 shows the composites of 850 hPa wind and OLR for these two distinct regimes: (8a) regime 1 and (8c) regime 2. The patterns are consistent with negative and positive modes of EOF1 of OLR, respectively.There is a clear regime behaviour in the variability of the Somali Jet with cyclonic/anticyclonic circulations over the northern Arabian Sea during regime 1/regime 2. In regime 1, OLR decreases over peninsular India with an increase in positive rainfall anomalies.The 850 hPa wind field presents a convergence zone over India with strengthening of easterly winds over the Indonesian and central Pacific.This pattern is reversed during regime 2 (Fig. 8c).There is evidence of a negative OLR anomaly over the western equatorial Pacific during regime 2. This pattern is consistent with the observed enhanced rainfall over this region and drought over India during the El Niño events (Krishnan Kumar et al., 2006).These distinctive results agree well with the PC1 time series of PDSI that indicates more "active" phases of monsoon during the LIA.The composites of modelled PDSI for the two regimes also present clear patterns resembling "active" and "break" monsoon phases (Supplement).
To examine the convection regimes of the five megadroughts, we presented the composite (84 years in total) anomalies of 850 hPa wind and OLR for the total period of five mega-droughts (Fig. 9a).During the past millennium, Asian mega-droughts are coincident with increased (decreased) OLR over East India, southern and eastern China

GISS-E2-R
Figure 7 presents the same results as in Sect.4.2.1, but for the GISS-E2-R model simulations.The skewness of PC1 of OLR data in the ensemble average of the GISS-E2-R model is 0.1.Experiment r1i1p121 indicates that the probability of occurrence of regime 2 was greater than for regime 1 (Fig. 7c and d).However, the ensemble average indicates more "active" regimes during the past millennium (α = 0.7, σ 1 = 1, µ 1 = −0.23,σ 2 = 1 and µ 2 = 0.5).The composite anomalies of OLR indicate negative deviations over the Indonesian Pacific, the Arabian Sea and India and positive anomalies over the central equatorial Pacific during regime 1 (Fig 8 .b).The lower tropospheric trade winds during regime 1 are strengthening over the central equatorial Pacific.This pattern is reversed for regime 2 (Fig. 8d). Figure 9b shows similar pattern as in Fig. 9a.A notable weakening of the Somali Jet occurred during the mega-droughts in GISS-E2-R model simulations with positive OLR anomalies over India.The 850 hPa wind pattern of the GISS-E2-R model presents a clear reduction of moisture transport into India, especially for the Somali Jet (Fig. 9b).
The composites of modelled PDSI for the two regimes from GISS-E2-R model present a less clear "active" monsoon phase for regime 1 (Supplement).The ensemble aver- age of the GISS-E2-R model presents a peak near neutral conditions (PC1=∼ 0) in Fig. 7b.

Discussions and conclusions
Our analysis of the climate during the past millennium is based on ensemble simulations of two state-of-the-art coupled atmosphere-ocean models and published proxy reconstructions of the climate (i.e.tree-ring reconstructions from Cook et al., 2010).
We applied EOF analysis of PDSI to identify the monsoon failure periods within the model-proxy space.Then MCA of PDSI and SSTA was used to capture the coupled spatiotemporal modes of atmosphere-ocean interactions during the past millennium.Finally, by using a statistical mixture model method on an index of convection, we investigated the regime structures in monsoon convection during the past millennium.
When the PDSI is used as metric, the two model and the proxy data sets agree in terms of timing and duration of the Asian mega-droughts.Time expansions of ENSO-like patterns agree well with Northern Hemisphere temperature anomalies based on palaeoclimate data (Shi et al., 2013), both in magnitude and timing.Our analysis showed that the mega-droughts are mostly linked to El Niño-like patterns in the central equatorial Pacific.The linkage between droughts and temperature may originate from PDSI, which contains information about the simulated temperature.
The IPCC Fifth Assessment Report (IPCC, 2013) shows that volcanic and solar forcing have a weak influence on monsoon.Adams et al. (2003) and Mann et al. (2005) concluded that ENSO may have been influenced by solar and volcanic radiative forcing.The correlation between the time expansion of MCA1 mode for SSTA of the ECHAM5/MPIOM (GISS-E2-R) model the Northern Hemisphere temperature reconstruction from Shi et al. ( 2013) is 0.57 (0.51).These time expansions show distinctive negative anomalies after explosive tropical volcanic eruptions for the mid-15th century (Pinatubo eruption) and early 19th century (Tambora, Galunggung, Babuyan Claro and Cosiguina volcanic eruptions).These anomalous cooling SSTAs are consistent with the reconstructed temperature changes in the Northern Hemisphere (Shi et al., 2013).The early part of the 1800s was the coldest period during the past four centuries (D'Arrigo and Wilson, 2009).The simulations show that the volcanic eruptions in the early 19th century were the major cause of reductions in ocean heat content (similar to LIA).
During the mega-droughts, there is evidence of decreased monsoon convection over India and the northern Arabian Sea and increased convection in the western and central equa-torial Pacific.The MCA1 patterns for SSTAs are in good agreement with the correlation patterns of MADA and observed SSTAs for the period 1856-2004 (Cook et al., 2010).In agreement with England et al. (2014), we conclude that the weakening in the Pacific trade winds and Somali Jet explains the decrease in monsoon convection.Weakening of the Pacific trade winds leads to El Niño-like patterns of SSTAs over the central and eastern equatorial Pacific by decreasing the ocean upwelling in these regions ("wind-driven circulation").This increases the surface temperatures over the central and eastern Pacific, which further drives the Asian monsoon failures.
An increase in the likelihood of very wet conditions is projected over the Indian monsoon region for the next century under the scenarios of rising greenhouse gas concentration (Dai, 2013).Further studies are necessary to investigate the ENSO-monsoon relationship under increased greenhouse gas concentrations.It needs to be stressed that the methods applied in this study can only capture the linear hydro-climatic interactions of the complex Asian summer monsoon system.The nonlinear mechanisms require further investigations (Dai, 2013;Hannachi and Turner, 2013).
The Supplement related to this article is available online at doi:10.5194/cp-11-253-2015-supplement.

Figure 1 .
Figure 1.Asian monsoon failure patterns: (a) first EOF of PDSI for MADA and (b) for the ensemble mean of GISS-E2-R and (c) ECHAM5/MPIOM.Pattern correlation coefficients between (a) and (b) and between (a) and (c) are 0.57 and 0.78, respectively.

Figure 2 .
Figure 2. Time expansions of monsoon failure patterns for (a) MADA (red) and the ensemble mean of GISS-E2-R (blue), and (b) MADA (red) and the ensemble mean of ECHAM5/MPIOM (blue).Time series have unit variance and zero mean and are smoothed using a 31-year moving average filter.Green (brown) shadings indicate the times when both time series are greater (smaller) than zero.Five historically recorded mega-drought periods are shaded in yellow.

Figure
Figure 6a presents the first EOF pattern of summer OLR for ECHAM5/MPIOM ensemble mean and Fig. 6b-g the histogram of the first principle component (PC1) of OLR for ensemble mean and each of the experiments (mil0010 to mil0015).The EOF1 pattern of OLR (Fig 6a) shows a positive anomaly pattern over the western and central North

Figure 6 .
Figure 6.Regime behaviour in ECHAM5/MPIOM model simulations.(a) Dominant OLR pattern (EOF1).Positive (negative) contours are shown by solid (dotted) lines.(b) PDF estimate of PC1 of OLR (red solid line) and each regime (dashed lines) for ensemble average and (c-g) for each of the experiments separately.

Figure 7 .
Figure 7. Regime behaviour in GISS-E2-R model simulations.(a) Dominant OLR pattern (EOF1).Positive (negative) contours are shown by solid (dotted) lines.(b) PDF estimate of PC1 of OLR (red solid line) and each regime (dashed lines) for ensemble average and (c-d) for each of the experiments separately.

Figure 8 .
Figure 8. Model simulations' composite anomalies of 850 hPa wind (m s −1 ) and OLR (W m −2 ) for the first (a and b) and second (c and d) regime from ECHAM5/MPIOM (left panels) and GISS-E2-R (right panels).The largest wind vector in the left panels is 0.6 m s −1 ; in the right panels it is 0.1 m s −1 .

Figure 9 .
Figure 9. Model simulations' composite anomalies of 850 hPa wind (m s −1 ) and OLR (W m −2 ) for the five mega-droughts (84 years in total) of the past millennium from (a) ECHAM5/MPIOM and (b) GISS-E2-R.The largest wind vector in (a) is 0.33 m s −1 ; in (b) it is 1.02 m s −1 .

Table 1 .
Correlation coefficients between the MCA1 time series, Niño indices and reconstructed global temperatures from Shi et al. (2013).Bold numbers indicate the largest value in each row.Numbers in parentheses indicate the statistical significance levels.NS stands for not significant.