The improbable but unexceptional occurrence of megadrought clustering in the American West during the Medieval Climate Anomaly

The five most severe and persistent droughts in the American West (AW) during the Common Era occurred during a 450 year period known as the Medieval Climate Anomaly (MCA—850–1299 C.E.). Herein we use timeseries modeling to estimate the probability of such a period of hydroclimate change occurring. Clustering of severe and persistent drought during an MCA-length period occurs in approximately 10% of surrogate timeseries that were constructed to have the same characteristics as a tree-ring derived estimate of AW hydroclimate variability between 850 and 2005 C.E. Periods of hydroclimate change like the MCA are thus expected to occur in the AW, although not frequently, with a recurrence interval of approximately 11 000 years. Importantly, a shift in mean hydroclimate conditions during the MCA is found to be necessary for drought to reach the severity and persistence of the actual MCA megadroughts. This result has consequences for our understanding of the atmosphere-ocean dynamics underlying the MCA and a persistently warm Atlantic Multidecadal Oscillation is suggested to have played an important role in causing megadrought clustering during this period.


Introduction
Proxy records suggest that the American West (AW) (125°W-105°W, 25°N-42.5°N; hereinafter AW) has experienced multidecadal drought, the severity and persistence of which is beyond the range of hydroclimate variability observed over the instrumental interval (1870-present). While these 'megadroughts' have been a prominent, albeit infrequent, feature of the Common Era (C.E.-the last 2015 years) paleoclimate record in the AW (e.g. Stine 1994, Stahle et al 2000, Cook et al 2004, the occurrence of megadroughts during the Medieval Climate Anomaly (MCA-850-1299 C.E.) is atypical: the five most severe and persistent droughts since 850 C.E. occurred during this short time period (Coats et al 2015a). The atmosphere-ocean dynamics during the MCA therefore have long been a focus of scientific research, with variability in the Pacific and Atlantic Oceans having been invoked to explain the MCA hydroclimate conditions (e.g. Graham et al 2007, Feng et al 2008, Oglesby et al 2011. Nevertheless, much remains to be understood about the MCA, particularly with regard to the relative influence of forced and internal variability in creating the Pacific and Atlantic Ocean conditions underlying MCA hydroclimate, the characteristics of these ocean conditions including their persistence, and the likelihood of such ocean conditions occurring-including in the future. Herein we use a tree-ring derived estimate of AW hydroclimate variability (the North American Drought Atlas (NADA)-e.g. Cook et al 2004Cook et al , 2007 in a timeseries modeling framework to better understand the clustering of megadroughts during the MCA. Two fundamental questions are addressed: (1) Given the characteristics of the last 1156 years of hydroclimate variability in the AW (850-2005 C.E.), what is the likelihood of the five most severe and persistent droughts falling within an MCAlength (450 year) period?
(2) Are certain characteristics of the last 1156 years of hydroclimate variability important for producing clustered megadroughts during the MCA and what does this suggest for the dynamics underlying hydroclimate conditions during this period?
By answering these questions we seek to better understand whether several consecutive centuries of heightened megadrought occurrence during the MCA could arise from our current best estimates of internal variability within the climate system or whether it requires external forcing (e.g. solar and volcanic) and, if not, what modes of climate variability can generate multicentury clustering of megadroughts. Together this information improves understanding of future megadrought risk over the AW.

Data and methods
Reconstructed Palmer drought severity index (PDSI) data are from an updated version of the tree-ring derived NADA version 2a, with improved spatial coverage and resolution, the full details of which can be found in Cook et al (2014a). The data are reconstructed on a 0.5°×0.5°latitude-longitude grid of JJA average PDSI values for the United States, as well as parts of Canada and Northern Mexico.
The AW is defined as the region bounded by 125°W -105°W, 25°N-42.5°N. This definition is consistent with previous research on megadroughts (e.g. Coats et al 2013Coats et al , 2015a and variability within the AW is largely homogenous with an average correlation coefficient of 0.74 between the AW averaged timeseries from the NADA and the individual NADA grid points over the AW for the period 850-2005 C.E. The NADA is thus averaged over the AW to create a single timeseries of hydroclimate variability for the period 850-2005 C.E. that will be analyzed herein. The most severe and persistent droughts within this timeseries are identified and ranked using the two start, two end (2S2E) and cumulative drought severity ranking criteria of Coats et al (2013). These criteria define a drought as commencing after two consecutive years of negative PDSI and continuing until two consecutive years of positive PDSI, with the droughts then ranked by summing the PDSI from the first to the last year of each identified drought feature. For the analyses herein, focus will be restricted to the five highest-ranking (or most severe and persistent) droughts. The cumulative drought severity ranking was chosen over a purely length-based ranking in order to incorporate both the severity and persistence of each drought. The 2S2E and cumulative drought severity method will hereinafter be referred to as the drought identification metric.
All surrogate timeseries are produced using phase randomization (e.g. Schreiber and Schmitz 2000), which takes the continuum of variability described in the frequency domain (power spectral density and phase as a function of frequency) and inverts this information into the time domain after randomizing the phase. In all cases, the power spectral density used for producing surrogate timeseries by phase randomization will be estimated using the Thomson's multitaper method (Thomson 1982). The outputs of phase randomization are timeseries with the time history randomized and the magnitude, variance and autocorrelation properties of the prescribed continuum of variability preserved.
To assess the impact of precipitation and temperature changes on PDSI we also calculate PDSI for 500 year control simulations from the Coupled Model Intercomparison Project Phase 5 (CMIP5-Taylor et al 2012; table 1). Three models are used to test if the characteristics of individual models alter the impact of precipitation and temperature on PDSI; the specific models were chosen to represent a range of characteristics of low-frequency hydroclimate variability (e.g. Coats et al 2015a). Herein we use the Penman-Monteith (PM-Penman 1948) formulation of PDSI (a full treatment of the PDSI formulation can be found in Cook et al 2014b). At each grid point on an even 2.5°×2.5°grid, PDSI was calculated and then standardized against the full 500 year period of the simulations. The PDSI was averaged over June-July-August (JJA) to produce a single average for each year; hereinafter any mention of model PDSI will be with regard to the JJA average values.   (2015a), the drought identification metric was used to identify the five highest-ranking droughts (methods) in the NADA AW timeseries (hereinafter megadroughts) and all five features fall within the MCA. Apart from this, the NADA AW timeseries in figure 1 (top panel) exhibits no obvious differences between the MCA and the 1300-2005 C.E. period (hereinafter the little ice age/modern period-LIA/ Mod). Nevertheless, to explicitly define the characteristics of the MCA relative to the LIA/Mod period, figure 1 also plots the time history of the autocorrelation function (ACF), the scaling exponent of spectral density (beta), variance and mean of AW hydroclimate over the C.E. The latter three characteristics are calculated by sliding windows of varying lengths, from 50 to 500 years, across the NADA AW timeseries. Beta is a measure of the proportion of variance in high and low frequencies, with positive beta values indicating that low frequencies contribute more to the variance of the underlying timeseries than high frequencies; negative beta values indicate the opposite (e.g. Huybers and Curry 2006, Ault et al 2013Ault et al , 2014. There is no consistent autocorrelation structure in the NADA AW timeseries across the range of lags considered in the ACF panel of figure 1. There is, however, a vertical line of enhanced autocorrelation at the termination of the MCA that may be indicative of a shift in the mean hydroclimate conditions between the two periods (see the final paragraph of this section for further discussion).

Assessment of AW hydroclimate variability
For short time windows (50-to-100 years), the beta, mean and variance exhibit a large range in values during the C.E. For beta, these shifts are related to secular trends over the short time windows. The behavior of the variance and mean are less clear, but again appear to be a function of the short time windows for which shifts in the mean state of the NADA AW timeseries are possible.
The longer time windows show no change in beta over the C.E. and, as such, the MCA to LIA/Mod transition is not likely to be explained by a change in . Internannual variability is plotted as the black bar plot and the 20 year low-pass filtered data is plotted as the blue line. The timing of the five most severe and persistent droughts (highest-ranking) as identified using the drought identification metric (methods) is indicated by the gray shaded regions. (second panel from top) The autocorrelation function for the NADA AW timeseries. This is calculated out to 50 years by sliding a 50 year window across the NADA AW timeseries (only the first 20 years, however, are plotted). (bottom three panels) The variance, mean and beta calculated for sliding windows varying in length from 50 to 500 years. Beta is calculated using ordinary least squares in log power-log frequency space on spectra estimated using Thomson's multitaper method (Thomson 1982). In all panels the MCA is defined as the time before the black vertical line (850-1299 C.E.). the spectral properties of the two periods. The variance, likewise, appears to be stationary across the NADA AW timeseries, although a weak shift from lower variance in the MCA to higher variance in the LIA/Mod is apparent. This weak variance shift can potentially be explained by the details of the reconstruction itself, as fewer tree-ring records are available further back in time in the NADA.
The mean hydroclimate conditions are also drier during the MCA as compared to the LIA/Mod, although the magnitude of this shift in PDSI units is only about 0.3 (the variance of the NADA AW timeseries, for reference, is 3.0). To determine if the lower mean during the MCA is indicative of a mean shift rather than a longer-term trend we use the single breakpoint (abrupt change) identification and significance test of Lanzante (1996). A highly significant (at the 99% level) breakpoint is identified at the year 1298 C.E. (just before the MCA termination of 1299 C. E.), suggesting that the MCA had a drier mean climate than the LIA/Mod and that the different mean climates of these two periods were separated by an abrupt change rather than being the consequence of a steady linear trend (e.g. Rodionov 2004).
Likelihood assessment of MCA megadrought clustering Experimental setup Figure 2 plots the raw spectra of the NADA AW timeseries, in addition to the spectra of the NADA AW timeseries with the mean of the MCA removed. An idealized red noise spectrum with the variance and first-order autoregressive coefficient matching the NADA AW timeseries is also plotted in figure 2. The shapes of the raw spectrum of the NADA AW timeseries and the spectrum of the NADA AW timereseries with the mean of the MCA removed are largely characteristic of a red noise process. Nevertheless, removing the mean of the MCA has the effect of significantly decreasing the power spectral density on multi-centennial-to-millennial timescales (at the 90% significance level).
To produce surrogate timeseries, the spectra in figure 2 are inverted into the time domain using standard phase randomization techniques (methods). An ensemble of 10 000 of these 1156 year artificial climate histories (length of the NADA AW timeseries) are produced and subsequently used to estimate the probability of occurrence of MCA megadrought clustering. Differences in the characteristics of the surrogates produced using the raw spectrum of the NADA AW timeseries and the spectrum calculated after removing the MCA mean shift (figure 2), isolate the impact of lowfrequency variability like the MCA mean shift.
The surrogates are used to test the probability that, given the characteristics of the NADA AW timeseries, the five highest-ranking droughts in a 1156 year climate history will fall within a 450 year period (length of the MCA). Hereinafter, MCA-like will be used to describe the clustering of the five highest-ranking droughts during a 450 year subset of a 1156 year surrogate timeseries. We also calculate the average length (drought persistence) and cumulative drought severity (drought severity) of the five highest-ranking droughts for each surrogate timeseries with an MCA-like period. These provide an additional constraint on whether a surrogate timeseries contains a period that is truly representative of the MCA, with droughts that are not only clustered but also characteristic of the severity and persistence of the actual MCA megadroughts.
Experimental results 10% percent of surrogates based on the raw spectrum of the NADA AW timeseries contain an MCA-like period (hereinafter percent occurrence). In the surrogates based on the spectrum with the MCA mean shift removed the percent occurrence is much lower (~three percent) suggesting that low-frequency variability like the MCA mean shift acts to increase the probability of a surrogate timeseries containing an MCA-like period. While this result is unsurprising, figure 3 suggests that without low-frequency variability like the MCA mean shift, the severity and persistence of drought never reaches the level of the actual MCA megadroughts (box plots do not reach the dashed lines in figure 3). In order to get a period that is truly representative of the MCA, therefore, the shift in the mean hydroclimate conditions during the MCA was necessary.

The origin of the MCA mean shift
The results in figure 3 suggest that low-frequency variability manifest in features such as the MCA mean shift in figure 1 is necessary to get drought that is characteristic of the actual MCA megadroughts during an MCA-like period. As was noted in the Introduction, much remains to be understood about the dynamics of the MCA and particularly the dynamics underlying megadroughts during this period. Nevertheless, research has implicated both a centuries-long warm state of the Atlantic Ocean via the Atlantic Multidecadal Oscillation (AMO-e.g. Feng et al 2008, Oglesby et al 2011 and cold state of the tropical Pacific Ocean (e.g. Graham et al 2007 in driving drying over North America during the MCA. In this section we will test whether these Atlantic and/or tropical Pacific Ocean conditions can explain the −0.31 mean shift in PDSI during this period (figure 1).
PDSI is determined by two factors: precipitation and evapotranspiration. In order to have a shift in the mean of PDSI during the MCA, one or both of these factors must have changed relative to the  C.E. interval over which the NADA PDSI was calibrated and standardized. Are the canonical impacts of a warm AMO or cold tropical Pacific sufficient to produce such changes? To address this the canonical impact of a warm AMO and cold tropical Pacific on both precipitation and temperature are defined using the instrumental record. The specific methodology will be outlined below, but we will briefly outline two assumptions of this test. Firstly, PDSI is a function of more than just precipitation and temperature, but only changes in these two variables will be analyzed. Secondly, we will consider any precipitation or temperature changes to be uniform in time and across the AW.
To test the impact of different precipitation and temperature changes on PDSI we superimpose uniform precipitation and temperature changes on model output from three control simulations (i.e. shifting the mean of the precipitation and temperature model output while keeping the variability unaltered). In each case, we calculate the PDSI using a PM framework with all variables besides precipitation and temperature being the unaltered model output. For each combination of precipitation and temperature changes we will calculate the mean PDSI over the AW for the first 450 years (the length of the MCA) of the control simulations. Results are plotted in figure 4, with good agreement across the three models. Unsurprisingly, as precipitation decreases the mean PDSI decreases and a 0.01 mm d -1 change in precipitation is associated with an approximately 0.1 mean shift in the The surrogates for the NADA AW timeseries are produced using spectra estimated with the mean of the MCA removed (MCA Mean Removed) and using the raw data (Raw). All results show the full data range, with the box plots additionally showing the median, interquartile, and 95th and 5th percentile ranges. In both cases the value for the five highest-ranking droughts (the megadroughts) in the actual NADA AW timeseries is plotted as the dashed black line.
PDSI. As temperature decreases the mean PDSI increases, although at a slower rate than that associated with precipitation (a 0.05°C change in temperature is associated with less than a 0.05 mean shift in the PDSI). This slower rate is consistent with PDSI predominantly reflecting precipitation variability over the AW (e.g. St. George et al 2010) and the opposite sign of the PDSI mean shift (relative to precipitation) is consistent with decreased temperature leading to a decreased vapor pressure deficit (and an associated decrease in evapotranspiration), following the Clausius-Clapeyron relationship acting on the saturation vapor pressure (e.g. Cook et al 2014b, Scheff and Frierson 2015).
To specifically test the potential that a persistently warm AMO or cold tropical Pacific produced the MCA mean shift in the NADA AW timeseries (figure 1), we calculate the canonical precipitation and temperature impacts of these ocean states. To do so, we composite detrended (to prevent aliasing a climate change signal) annual (January-December) precipitation from the Global Precipitation Climatology Center  iod 1901-2005 C.E. The temperature and precipitation composites are then averaged over the AW to produce the canonical impacts over the region. For precipitation this impact is −0.036 mm d -1 (a 3% decrease in annual precipitation) for a warm AMO and −0.064 mm d -1 (a 6% decrease) for a cold tropical Pacific and for temperature it is 0.17°C and −0.06°C, respectively. The composite across models in figure 4 suggests that a warm AMO would dry the AW slightly beyond that of the MCA mean shift (of −0.31, figure 1). Given that temperatures over the AW were likely cooler during the MCA than implied by the canonical impacts of a warm AMO-because of cooler than modern Northern Hemisphere temperatures (Moberg et al 2005, Hegerl et al 2007, Mann et al 2009 -the mean PDSI values in figure 4 should be even closer to the MCA mean shift. By contrast, the canonical impact of a cold tropical Pacific would dry the AW far in excess of the MCA mean shift. Importantly, however, weaker cooling in the tropical Pacific-that is to say, much smaller in magnitude than the cold tropical Pacific analyzed herein-would produce drying that is closer to the MCA mean shift. A linear regression of the Niño3.4 index on precipitation over the AW, for instance, suggests that a Niño3.4 index of −0.3°C (32% of the average Niño3.4 index for a cold tropical Pacific) produces drying approximately equal to that of the MCA mean shift.

Discussion of climatic implications and conclusions
A period with clustering of severe and persistent drought like the MCA is not a highly probable event, nevertheless such a feature occurs in approximately 10% of surrogate timeseries that were constructed to have the same spectral characteristics as the NADA AW timeseries. This suggests that while megadrought clustering during the MCA was improbable, with a recurrence interval of approximately 11 000 years, it is a naturally occurring phenomenon that is expected to occur over the AW. Importantly, a shift in the mean hydroclimate conditions during the MCA (850-1299 C.E.) is found to be a critical characteristic of the NADA AW timeseries. Without low frequency variability capable of generating such a mean shift, the probability of drought clustering like the MCA drops to approximately three percent and, more importantly, drought never reaches the severity or persistence of the actual MCA megadroughts.
PDSI over the AW calculated after superimposing the precipitation and temperature impacts of a weak cooling in the tropical Pacific or the canonical impacts of a warm AMO can largely reproduce the shift in mean hydroclimate conditions during the MCA. A centuries-long persistence of one of these ocean states, therefore, may have been a necessary condition for the clustering of severe and persistent drought during the MCA. We hypothesize that there was a persistently warm AMO during the MCA (as proposed by Feng et al 2008, Oglesby et al 2011, in particular, because the AMO is associated with much longer timescales of variability (with a quasi-periodic cycle of 50 to 70 years -e.g. Mann and Park 1994) than the tropical Pacific Ocean, which is dominated by seasonal-to-decadal timescales. For instance, in an additional analysis using an ensemble of 10 000 1156 year surrogate timeseries with the same power law relationship (beta) as the observed AMO, 41% of the surrogate climate histories contain a 450 year period (the length of the MCA) where at least 75% of years are in the top third of observed AMO values. Perhaps unsurprisingly, given the long timescales of variability inherent to the AMO, this suggests that the observed characteristics of this mode of variability-whether forced or internalare sufficient to explain a shift towards a persistently warm state during the MCA.
Importantly, the AMO is only being implicated in the mean hydroclimate shift during the MCA (e.g. figure 1) and the timing of individual MCA megadroughts may still be related to decadal periods of cold conditions in the tropical Pacific Ocean, which occur frequently in observations and can produce severe drying over the AW (e.g. McCabe et al 2004). Because it is often thought that the AMO is an internal mode of variability, and there is little doubt that decadal variability of the tropical Pacific Ocean is, it follows that MCA megadroughts could have arisen via purely internal climate dynamics despite widespread conjecture that they were exogenously forced (e.g. Graham et al 2007, Mann et al 2009. The results presented here can neither prove nor disprove either hypothesis but emphasize the need to better characterize the ocean states during the common era. Such an understanding is particularly important when considering future drought conditions over the AW. If a persistently warm or cold AMO occurs as a response to anthropogenic greenhouse gas forcing, this would have consequences for the probability of MCAlike drought occurring in the future. The models used for future projections, however, are unlikely to simulate the worst-case scenario of a persistently warm AMO and decadal periods of cold conditions in the tropical Pacific because few models have been shown to simulate megadroughts consistently associated with tropical Pacific Ocean (Coats et al 2015a) and most models lack a realistic AMO and associated hydroclimate impacts over the AW (Ting et al 2011, Coats et al 2015b. These model deficiencies can perhaps explain the inability of the CMIP5 models to produce MCA-like drought clustering, despite simulating individual megadroughts that are as severe and persistent as those in the paleoclimate record (Coats et al 2015a). In any case, efforts to characterize the dynamics of past hydroclimate change are critical because they suggest that the CMIP5 ensemble is unlikely to capture the full range of potential future hydroclimate states (as also noted in Ault et al 2013Ault et al , 2014; although the opposite has been implied-Frank et al 2013). Further attempts to characterize model projections in the context of long records of climate variability will therefore continue to improve the basis for accurate predictions of the future climate state.