Abstract

The western United States was affected by several megadroughts during the last 1200 years, most prominently during the Medieval Climate Anomaly (MCA; 800 to 1300 CE). A null hypothesis is developed to test the possibility that, given a sufficiently long period of time, these events are inevitable and occur purely as a consequence of internal climate variability. The null distribution of this hypothesis is populated by a linear inverse model (LIM) constructed from global sea surface temperature anomalies and self-calibrated Palmer drought severity index data for North America. Despite being trained only on seasonal data from the late twentieth century, the LIM produces megadroughts that are comparable in their duration, spatial scale, and magnitude to the most severe events of the last 12 centuries. The null hypothesis therefore cannot be rejected with much confidence when considering these features of megadrought, meaning that similar events are possible today, even without any changes to boundary conditions. In contrast, the observed clustering of megadroughts in the MCA, as well as the change in mean hydroclimate between the MCA and the 1500–2000 period, are more likely to have been caused by either external forcing or by internal climate variability not well sampled during the latter half of the twentieth century. Finally, the results demonstrate that the LIM is a viable tool for determining whether paleoclimate reconstructions events should be ascribed to external forcings or to “out of sample” climate mechanisms, or if they are consistent with the variability observed during the recent period.

1. Introduction

Megadroughts are prolonged periods of aridity typically defined by their multidecadal duration (Woodhouse and Overpeck 1998; Ault et al. 2014; Cook et al. 2016) and have been linked to the demise of several preindustrial civilizations (Benson et al. 2002, 2003; Hodell et al. 1995; Buckley et al. 2010; Stahle and Dean 2011). Mounting evidence suggests that the risk of these events is increasing throughout the southwestern United States and other subtropical dry zones due to rising temperatures and dynamic circulation changes (Seager et al. 2007b; Woodhouse et al. 2010; Ault et al. 2014, 2016; Cook et al. 2015, 2016; Kelley et al. 2015). If a megadrought were to unfold in the coming decades, it would require water management decisions to be made on scales that have no historical precedent. Understanding the risks of these events is therefore essential for supporting regional adaptation strategies in the face of a changing climate. Furthermore, quantitative assessments of megadrought risk require us to disentangle the influences of external and internal climate variability in the past because future likelihoods will be governed by both of these factors (Ault et al. 2014, 2016; Cook et al. 2015).

The last millennium affords some insights into the range of forced and unforced preindustrial climate variability, yet it is only one realization of many plausible alternative climate trajectories that could have occurred during that same time period. If a given climate forcing (e.g., enhanced solar output) made megadroughts more probable in the past, but they failed to occur due to internal climate variability, then one might erroneously conclude that the forcing was not important to the statistics of megadrought on the basis of the one realization being examined. Alternatively, if megadroughts occur without any clear external forcing, yet are nonetheless more probable during certain climate regimes, one might correctly conclude that megadroughts do not require exogenous changes to climate boundary conditions, but would miss the critical piece of information that such forcings can make megadroughts more likely. These inherent constraints of having only one realization of reality have largely prevented a consensus view from emerging on whether or not megadroughts are forced [see the review by Cook et al. (2016)].

A number of studies have suggested that megadroughts in the western United States, and especially in the southwestern United States, require external forcing mechanisms to alter regional and global circulation patterns in ways that have not been observed over the historical period (e.g., Clement et al. 1996; Graham et al. 2007, 2011; Mann et al. 2009; Seager et al. 2007a). Evidence for this hypothesis stems in part from the observation that the Medieval Climate Anomaly (MCA; 801–1300 CE) occurred during a period of slightly enhanced solar activity (e.g., Lean 2000; Vieira and Solanki 2010) and saw at least three multidecadal intervals of aridity with drought area index (DAI) values above the worst droughts of the twentieth century (Cook et al. 2004, 2016; Coats et al. 2016b). Idealized climate model experiments support this interpretation to varying degrees: when the Cane–Zebiak model (Zebiak and Cane 1987) of tropical Pacific Ocean dynamics is forced with increased solar activity, it simulates a “La Niña–like” mean state (Clement et al. 1996). Prolonged La Niña–like sea surface temperatures (SSTs), in turn, favor relatively arid conditions across the southwestern quadrant of the United States (Seager et al. 2007a; Burgman et al. 2010). In addition to external forcing factors, there is also evidence that regional feedbacks related to dust and land surface dynamics could make megadroughts more probable when the mean state of the western U.S. hydroclimate is drier than today (e.g., Cook et al. 2011, 2016).

Under the interpretations above, prolonged and widespread megadroughts arise from climate regimes that are inherently out-of-sample realizations of forced hydroclimate states not observed during the relatively short instrumental era. Analysis of unforced fully coupled general circulation model (GCM) simulations, in contrast, suggests that megadroughts could occur as a consequence of internal climate variability, if given a sufficiently long period of time. For example, Hunt (2006) argued that megadroughts should be interpreted as normal events in the unforced statistics of a long (10 000 yr) control run. Likewise, Meehl and Hu (2006) linked megadroughts in the U.S. Southwest to low-frequency modes of climate variability using a millennial-scale integration of a moderate-resolution GCM. More recently, Coats et al. (2013b, 2015a,b, 2016a,b) and Stevenson et al. (2015) argued that megadroughts occur as a result of internal variability, and hence they do not require an exogenous mechanism. Furthermore, fully coupled GCMs forced with realistic, time-evolving boundary conditions of the last millennium do not universally simulate clusters of megadroughts during the MCA (Cook et al. 2016; Coats et al. 2016a,b). This result could be interpreted as evidence that models fail to simulate forced megadroughts, or alternatively that these events occurred more frequently during the MCA due to chance (Coats et al. 2016b).

The GCM-based studies above are inherently limited by the rarity of megadroughts even in long control simulations, and by structural biases in GCM moisture fields (Ault et al. 2012; Flato et al. 2013) and teleconnections (Coats et al. 2013a). Even the largest ensembles (e.g., Kay et al. 2014; Otto-Bliesner et al. 2016) only include a few dozen realizations. If megadroughts tend to occur at a rate of 1–2 events per millennium [as originally proposed by Woodhouse and Overpeck (1998) and supported analytically by Ault et al. (2014)], then independent millennial-scale realizations still yield a relatively small sample size to assess megadrought statistics [as does a 10 000-yr control run as in Hunt (2006)]. Likewise, GCMs are known to have systematic biases both in their mean precipitation (Flato et al. 2013), as well as their tropical Pacific variability (Guilyardi et al. 2009; Ault et al. 2013b). Together, these biases imply that megadrought statistics estimated from GCMs are an imperfect approximation of the “true” limiting behavior of the system.

Here we develop an empirical approach for ruling out the possibility that megadroughts could arise from unforced, stationary statistics. Specifically, we use a linear inverse model (LIM) framework to emulate the observable spatiotemporal patterns inherent to the climate system over the instrumental period, when ocean–atmosphere data are widespread and of high quality. In climate studies, LIMs have been used to characterize whether or not the recent clustering of central Pacific El Niño events (also called “Modokis”) is anomalous in the context of El Niño–Southern Oscillation (ENSO) dynamics as they evolve naturally over the course of centuries (Newman et al. 2011). It has also been applied to test the significance of the substantial low-frequency (decadal to centennial) variability in the tropical Pacific (Ault et al. 2013b) and to evaluate significance of “peaks” in the power spectrum of Pacific decadal oscillation (PDO) reconstructions (Newman et al. 2016).

Our null hypothesis is that megadrought characteristics in the western United States are consistent with those of a stochastically forced, linearly damped system with stationary statistics. That is, over any given millennium, the interannual variability in global SSTs and western U.S. (WUS) hydroclimate, combined with spatial and temporal sources of autocorrelation, will produce megadrought statistics that are indistinguishable from reconstructions. We use an LIM to populate the null distribution associated with this hypothesis, which in turn permits us to test the significance of the following characteristics of WUS megadrought (Fig. 1).

  1. Events were widespread, encompassing between 20% and 40% of the domain at their peak (Fig. 1a; Cook et al. 2004).

  2. The worst event in the southwestern United States occurred during the late 12th century and was as severe as the interannual droughts of the twenty-first century, but persisted for decades (Cook et al. 2004, 2016). This event can be clearly identified in the 35-yr running mean of southwestern U.S. Palmer drought severity index (PDSI) in Fig. 1b (red circle).

  3. A “cluster” comprising five megadroughts occurred during a 500-yr period corresponding to the MCA (gray vertical bars in Fig. 1; Coats et al. 2016b; Cook et al. 2016).

  4. And, finally, the mean hydroclimatic state of the Southwest was drier on average during the MCA than any other subsequent period of time (red line in Fig. 1b; Cook et al. 2004, 2016; Graham et al. 2007; Coats et al. 2016b).

Fig. 1.

Key characteristics of WUS megadrought from 801 to 2000 CE. (a) DAI35 computed as the fraction of grid points with 35-yr average PDSI values below −1 (red dot indicates the most widespread event). (b) Smoothed (35-yr running mean) PDSI from the SWUS (32°–41°N, 124°–105°W). The gray vertical bars highlight five MCA droughts (35-yr averages below −0.5 PDSI units), with the lowest 35-yr averaged PDSI value (−1.04) circled in red (1159 CE). The long-term means of PDSI are plotted in red for MCA (−0.19 PDSI units) and 1500-present (0.26 PDSI units). (c) Average PDSI across the WUS domain during the most widespread event (1130–1164 CE). The map of 35-yr averaged PDSI during the most negative (worst) SWUS drought in (b) is very similar due to the map of the most widespread drought because the two periods overlap in time.

Fig. 1.

Key characteristics of WUS megadrought from 801 to 2000 CE. (a) DAI35 computed as the fraction of grid points with 35-yr average PDSI values below −1 (red dot indicates the most widespread event). (b) Smoothed (35-yr running mean) PDSI from the SWUS (32°–41°N, 124°–105°W). The gray vertical bars highlight five MCA droughts (35-yr averages below −0.5 PDSI units), with the lowest 35-yr averaged PDSI value (−1.04) circled in red (1159 CE). The long-term means of PDSI are plotted in red for MCA (−0.19 PDSI units) and 1500-present (0.26 PDSI units). (c) Average PDSI across the WUS domain during the most widespread event (1130–1164 CE). The map of 35-yr averaged PDSI during the most negative (worst) SWUS drought in (b) is very similar due to the map of the most widespread drought because the two periods overlap in time.

The list above serves as a basis for developing megadrought test statistics, the significance of which can then be evaluated against the LIM’s distribution of these same quantities. For example, if LIM realizations regularly reproduce the spatial scale, magnitude, and clustering of megadroughts, as well as changes in the long-term mean values of reconstructed PDSI, then we would not be able to reject the null hypothesis. Failing to reject this null hypothesis would imply that megadrought occurrence does not require external forcing mechanisms or climate behaviors that are any different from those observed during the twentieth century (e.g., a multicentury warm phase of the Atlantic multidecadal oscillation; Feng et al. 2008; Oglesby et al. 2012; Coats et al. 2016b; Wang et al. 2017). Rejecting the null hypothesis, on the other hand, would imply that the processes responsible for generating multidecadal periods of aridity in the western U.S. hydroclimate originate from either external forcings (e.g., solar variability, explosive volcanism, or slow-varying orbital changes), from low-frequency climate processes that are not well sampled over the twentieth century, or from nonlinear coupled processes and feedbacks at multidecadal and longer time scales. Accordingly, applying our null hypothesis to the list of western U.S. hydroclimate characteristics from the last 1200 years will help clarify which aspects of megadrought are fundamentally unusual, versus merely unexpected.

Finally, for clarity and consistency with previous studies, we make a distinction between hydroclimate of the entire western United States (WUS; 32°–50°N, 125°–105°W,), and the southwestern (SWUS) portion of that domain (32°–41°N, 124°–105°W). The dual foci of this study therefore permits analysis of both the spatial scale of megadroughts across the WUS and their severity and duration in the SWUS, the latter features having received somewhat more attention in the literature than the former. Moreover, among the leading modes of variability in WUS hydroclimate is a “dipole” pattern with dry and wet anomalies of opposite sign between the Pacific Northwest and the SWUS (e.g., Ropelewski and Halpert 1986; Cook et al. 2017), which further motivates analysis of both domains because megadroughts conditions in the SWUS are not necessarily guaranteed to be representative of WUS-wide conditions.

2. Linear inverse modeling

A linear inverse model presents a simple method to emulate the statistical characteristics of certain dynamical systems (Penland and Matrosova 1994; Penland and Sardeshmukh 1995). In addition to the examples cited in the previous section, the LIM approach has been used to make seasonal predictions of ENSO (Penland and Sardeshmukh 1995; Newman et al. 2011), to quantify predictability of the PDO (Newman 2007; Alexander et al. 2008; Newman et al. 2016), to characterize tropical Pacific ocean–atmosphere coupling across time scales (Newman et al. 2009), and to evaluate decadal forecasts made by fully coupled global climate models (Newman 2013). This suite of studies demonstrates that seasonal to decadal SST variability is well approximated by an LIM (Penland and Sardeshmukh 1995; Newman 2007, 2013; Alexander et al. 2008; Smirnov et al. 2014). Because SSTs, and especially tropical Pacific SSTs, exert a major influence on WUS hydroclimate (Ropelewski and Halpert 1986, 1996; Redmond and Koch 1991; Brown and Comrie 2004), we extend the LIM approach to emulate drought statistics in the WUS and the influence of SSTs on this region across seasonal to centennial time scales.

Using an LIM to generate null distributions for the statistics of widespread megadrought is similar to using an AR(1) model to evaluate the significance of behavior in a scalar time series (e.g., Bartlett 1978). The key difference is that the AR(1) model applies only to a univariate setting , with variations governed by an autocorrelation term (ρ) and stochastic white noise (η), such that ; it therefore cannot incorporate information about both the spatial and temporal scales of autocorrelation. An LIM, in contrast, is appropriate for simulating the variability produced by any multivariate system that evolves according to the following Langevin equation for an Ornstein–Uhlenbeck process (Penland and Sardeshmukh 1995):

 
formula

where X is an m-dimensional state vector and ζ is white-noise forcing term, which itself may have spatial (or temporal) structure, but whose realizations are uncorrelated from one time step to the next. The term is an m × m linear deterministic feedback matrix, which includes information about how stochastic excitations propagate (linearly) through the system. If the system is indeed approximately linear, the eigenvalues of L must have negative real parts; otherwise, it would contain modes that grow exponentially and violate our assumptions about stable linear statistics.

The deterministic linear feedback matrix, , can be estimated empirically from data as

 
formula

Here, the term n is a time lag parameter, is the lag covariance matrix (at lag n), and is the instantaneous covariance matrix. To the extent that the system is approximately linear, the choice of n does not critically alter the structure of because, by definition, the most probable state vector at some lead time is given by Penland and Sardeshmukh (1995) as

 
formula

In a typical application of an LIM, the noise term (ζ) is not known, but instead its covariance matrix is estimated as a residual using the fluctuation–dissipation theorem (Penland and Matrosova 1994):

 
formula

Further information on the assumptions, limitations, and details of implementing an LIM for any generalized system with approximately linear dynamics can be found in Penland and Sardeshmukh (1995) and Penland and Matrosova (1994). Importantly, LIMs partition variability into an approximately linear part, represented by , and a noise forcing part (). The variables contributing to the noise term () might exhibit nonlinear behavior at short time scales, and so an LIM approach is most appropriate for situations where there is temporal averaging inherent to the system of interest (X). Additionally, the time series included in are usually smoothed to emphasize the approximately linear sources of variability in the system (Penland and Sardeshmukh 1995). This averaging (and smoothing), in turn, causes multiple realizations of the nonlinear phenomena to be integrated into each time step, and hence exhibit Gaussian statistics via the central limit theorem (Yuval and Hsieh 2002). This assumption is demonstrably valid for monthly SST anomalies, as shown by many previous studies across weekly to multidecadal time scales (e.g., Penland and Sardeshmukh 1995; Newman et al. 2009, 2011; Newman 2013). As with SST, the drought indices of interest here integrate monthly moisture supply and demand, and hence their variations are well suited to be simulated by an LIM approach.

Given the similarity between Eq. (1) and the equation for univariate red noise, LIM output can be thought of as “multivariate red noise” (Newman et al. 2011; Ault et al. 2013b). Unlike univariate red noise, which smoothly increases in amplitude from short to long time scales (Bartlett 1978), multivariate red noise can be used to generate time series characterized by strong spectral peaks, as is the case in the tropical Pacific. These peaks arise from the linear dynamics of the system, even when it is forced by stochastic white noise, and hence yield an important baseline for variations generated in regions affected by phenomena dominated by quasi-periodic variability (e.g., ENSO).

3. A linear inverse model of western U.S. drought

a. Data

We use global SST anomalies and PDSI values to construct an LIM of the surface ocean and WUS hydroclimate. Monthly SST anomalies from 60°S to 75°N originate from the Kaplan et al. (1998) dataset. This product is generated at a native spatial resolution of 5° × 5° (500 km) lat/lon using a statistical procedure that estimates missing data through time from spatial empirical orthogonal functions (EOFs), statistical interpolation, and Kalman filtering (Kaplan et al. 1998). Although data are available as far back as 1856, we only use data from 1950–2000 for training the LIM. It is also detrended to minimize the contribution of anthropogenic forcing over that period, and smoothed prior to estimating (see section 3b).

Hydroclimate data for the WUS are obtained from the “self-calibrating” PDSI dataset developed and published by Sheffield et al. (2006). This version of the PDSI (termed scPDSIpm) was generated using the Penman–Montieth method for calculating potential evapotranspiration, and the University of East Anglia’s (UEA) Climate Research Unit’s (CRU) TS3.10 precipitation data product as input, as well as near-surface pressure, humidity, and radiation fields derived from reanalysis and remote sensing data (Sheffield et al. 2006). For comparison of the widespread multidecadal megadrought intervals identified by Cook et al. (2004), we restricted our use of this scPDSIpm database to the WUS region (32°–50°N; 125°–105°W). Again, the period from January 1950 through December 2000 was used to train the LIM. Sensitivities of our LIM to the underlying data are presented at the end of the section 4.

b. Data preparation

Following prior LIM studies (Penland and Sardeshmukh 1995; Newman et al. 2011), principal components (PCs) were used to reduce the dimensionality of each field, and a 3-month moving average was applied to smooth each time series prior to computing the PCs. The first 11 PCs were retained from global SST [as in Newman (2013)] and jointly account for 60% of the total variance in that field, although the RMSE between the filtered and raw data, as well as the total fraction of shared variability, can be locally much higher (Fig. 2). While the smoothing has the benefit of stabilizing our estimate of (i.e., making it less sensitive to noise), an obvious limitation is that the LIM cannot reproduce 100% of the variance at all locations (Fig. 2). This limitation is not a serious one for our applications because, locally, at least 90% of the variability in the tropical Pacific and the WUS is captured by the LIM. In other applications, however, such as seasonal prediction, this limitation could present problems if we needed high-accuracy LIM forecasts of future hydroclimate states.

Fig. 2.

(top) Root-mean-square error (RMSE) and (bottom) percent variance accounted for from smoothing and EOF filtering the raw (left) SST and (right) PDSI. White patches over the ocean are areas of missing data.

Fig. 2.

(top) Root-mean-square error (RMSE) and (bottom) percent variance accounted for from smoothing and EOF filtering the raw (left) SST and (right) PDSI. White patches over the ocean are areas of missing data.

The SST portion of the LIM preserves information about the growth and decay structures of tropical Pacific ENSO variability. Hence, it produces realistic interannual variability and Niño-3.4 power spectra, both of which are critical for emulating hydroclimate statistics in the WUS. The characteristics of the LIM were not particularly sensitive to this choice, and a somewhat larger or smaller number of PCs could have been retained to yield similar results.

In selecting the number of PCs to retain from the eigenvector decomposition of scPDSIpm, we required the total amount of variance in the field be comparable to the full range of the raw dataset (otherwise, any inferences about the statistics of megadrought in the LIM would have limited relevance to paleoclimate records of these events). We therefore retained the leading 22 PCs from scPDSIpm, which accounted for 90.8% of the raw (unsmoothed) variance of the WUS domain.

The final step in constructing X entails normalizing the PCs from SST and scPDSIpm. We required that so that variability in the scPDSIpm would not appreciably affect the structure of . Physically, this would be interpreted to mean that SST variability and regional sources of atmospheric noise can force scPDSIpm in the LIM, but scPDSIpm cannot affect SST. By requiring and , we found that the spatiotemporal structure of the SST patterns in were insensitive to scPDSIpm.

Our state vector is therefore defined as

 
formula

where and are the filtered SST and scPDSIpm sets of PCs, respectively. The state vector defines the instantaneous covariance matrix: , where the angled brackets indicate averaging. The matrix will have diagonal sections that are the covariances of the individual fields (i.e., , and ), and off-diagonal sections that are the covariances between SST and scPDSIpm fields:

 
formula

Finally, there is a well-known trend in global SSTs over the twentieth century, and possible (albeit more uncertain) trends in PDSI over the same period (e.g., Sheffield and Wood 2008; Sheffield et al. 2012; Dai 2011). Trends can influence the structure of because they will appear as weakly damped modes with decorrelation time scales proportional to the length of the training interval (1950 to 2000 in this case). This component, in turn, could impart ultralow-frequency variability on the dynamical system that would be unrepresentative of its unforced statistics. We took several steps to guard against this possibility. First, SST and scPDSIpm fields were both linearly detrended prior to analysis. Second, patterns and time series of the PCs obtained from both fields were examined to ensure they did not include any “residual” trend. Third, the eigenvalues of were inspected to make sure their decay time scales were relatively short (seasonal to interannual). And fourth, by focusing on the interval of 1950–2000, and especially by calibrating our LIM on the instantaneous and lag 3-month covariance matrices, we minimize the effects of secular and long-term variations because seasonal and interannual sources of variability dominate the fluctuations at these time scales.

The LIM was run for 12000 months (1000 yr) a total of 1000 times [using the formula provided by Penland and Sardeshmukh (1995), Eqs. (15a) and (15b) therein]. The PC time series generated from each of these stochastic realizations were projected onto the original subset of spatial loading vectors to generate spatially complete fields of each variable. Drought and tropical Pacific statistics (e.g., the DAI, spatial correlations, and power spectra) were then computed from this output. Although output from the LIM is available at monthly resolution, the PDSI field was averaged over June–August (JJA) for comparison with paleoclimate reconstructions of this same variable over the same season (see next section).

c. Reconstructed PDSI data

We evaluated the synthetic droughts produced by the LIM by comparing them against the North American Drought Atlas (NADA), which is a gridded set of mean JJA PDSI reconstructions derived from tree rings (Cook et al. 2004, 2014). The NADA uses point-by-point principal component regression to estimate JJA PDSI for the hundreds of years before instrumental data are available using nearby tree-ring data, drawing first upon those chronologies within a 450-km search radius of a given grid point, and then expanding the search by 50 km until at least five predictors have entered the model (Cook et al. 1999, 2007). As part of the reconstruction procedure built into the NADA, the tree-ring data are prewhitened to remove temporal autocorrelation believed to be due to biological memory and forest dynamics. After the reconstructions are generated, the estimates are “re-reddened” by adjusting the autocorrelation of the proxy drought series to match the memory observed in instrumental PDSI (Cook et al. 1999, 2007). Here, we used a high-resolution version of the NADA, which was recently updated to 0.5° latitude/longitude resolution (Cook et al. 2014) and incorporates a greater number of predictor chronologies. For comparison with the observational data (Sheffield et al. 2006), we regridded the high-resolution NADA to the same 1° lat/lon grid as the scPDSIpm product derived from observations.

Given our focus on multidecadal time scales, we refer to the 35-yr running average of PDSI as . The specific choice of 35 years for the time averaging is arbitrary; we adopt it here following earlier studies because it is simple and reliably identifies megadroughts as anomalies in preindustrial times in the WUS (Ault et al. 2014, 2016; Cook et al. 2016). It also satisfies earlier, more qualitative definitions of megadroughts as being periods of aridity as severe as the worst droughts of the twentieth century, but much longer lasting (e.g., Woodhouse and Overpeck 1998; Acuna-Soto et al. 2002), and it is amenable to analytical treatment [see the supplementary material of Ault et al. (2016)]. This quantity is computed from reconstructed JJA PDSI in the NADA and compared to the LIM output of JJA PDSI.

d. Megadrought test statistics

We identify megadroughts in the NADA and in the LIM based on the four features of these events described in the introduction (and summarized in Table 1). Megadroughts have been widespread, of comparable magnitude to the twentieth century in their severity, and clustered during the MCA. There has also been a shift in mean hydroclimatic state of 0.47 PDSI units since the MCA (Coats et al. 2016b), from relatively more arid to more humid modern conditions. Each of these characteristics are expressed as the quantitative test statistics described below:

  • Spatial scale. We generate a multidecadal DAI from reconstructions and LIM output by calculating the fraction of WUS grid points with at or below −1 PDSI units. To distinguish this quantity from the DAI employed by Cook et al. (2004), which first computed the fraction of the domain below −1 PDSI units, then smoothed the resulting time series, we term the DAI employed here “DAI35” to emphasize its multidecadal dependence. Standard deviations of the reconstructed PDSI are close to 2 across much of the WUS, so that the threshold used to identify events (−1) is generally consistent with earlier studies that used a threshold to identify events in time series normalized to unit variance over the historical period (e.g., Cook et al. 2004, 2016; Ault et al. 2014, 2016). The DAI35 may therefore be interpreted as the fraction of the WUS experiencing megadrought lasting 35 years or longer at a given time. It could theoretically vary between zero (no grid points in megadrought) and 100% (the entire WUS in megadrought). LIM output is used to compute the DAI35, which is then compared against the NADA’s DAI35 (see below).

  • Magnitude. The significance of multidecadal megadrought magnitude is assessed by comparing LIM output to the worst event in the mean (domainwide) SWUS between 801 and 2000 CE. This event reached a peak amplitude from 1124 to 1159 CE of −1.04 PDSI units on average (red circle in Fig. 1b).

  • Clustering. Five events of relatively low amplitude (−0.5 PDSI units, or approximately −0.25σ) and long duration (35 yr or more) occurred during a 500-yr period of the MCA (Fig. 1b). We tally the total number of events of equivalent amplitude in each possible consecutive 500-yr interval of each LIM realization to evaluate the significance of the clustering observed in the MCA (applying this same criterion to the reconstructed PDSI identifies the five low-amplitude megadroughts of the MCA; Fig. 1).

  • Mean shift. Mean PDSI values of the MCA were 0.47 PDSI units lower than they were during 1500–2000 CE. Differences between the first and last 500 years of each LIM realization are compared against this value.

Table 1.

Summary of key acronyms and quantities used in this study.

Summary of key acronyms and quantities used in this study.
Summary of key acronyms and quantities used in this study.

e. Evaluating significance

We generate empirical p values for each test statistic above by comparing its reconstructed value against the distribution of those same quantities generated by the LIM. The p value is therefore the empirical probability that a given megadrought test statistic was drawn from the LIM’s distribution, and is defined as

 
formula

Empirical p values defined in this way will be close to zero for megadrought characteristics that are unlikely to have been drawn from the null distribution of the LIM, and larger (closer to unity) when associated with events for which the null hypothesis cannot be rejected.

For example, the most widespread event from 801 to 2000 CE reached peak spatial extent (maximum DAI35) over the period from 1130 to 1164 CE (Fig. 1a), with 40.8% of the domain experiencing megadrought conditions ( below −1.0 PDSI units). We use the LIM to test if events of this scale are significant by computing a “fraction of time” test statistic over which megadrought conditions prevailed across some portion of the WUS domain. In the NADA, a total of 3 years (1/400 of all the years from 801 to 2000 CE) saw DAI35 values above 40%. The test statistic characterizing such an event at the DAI35 40% level would therefore be . The values of can likewise be calculated for any other DAI35 threshold.

Output from the 1000 LIM realizations is used to generate the null distribution of . For each run, the fraction of years above various DAI35 values is computed, yielding an empirical probability density function (PDF) of based on the entire LIM ensemble. NADA-based values of are then ranked against this empirical distribution to establish p values, and hence confidence limits for rejecting the null hypothesis. For example, if a reconstructed value of falls outside the range of 95% of the LIM-generated data for that combination of area and fraction of time, we would reject the null hypothesis at the 95% confidence limit.

As with our significance testing of NADA-based values of for the entire WUS domain, we evaluate each of the three other test statistics for the SWUS: the magnitude of SWUS minima , clustering of (low-amplitude) megadroughts during a 500-yr period , and the change in mean PDSI . Because the value of is negative (−0.62σ), we evaluate significance by comparing the distribution of the LIM with minimum values below this threshold. In contrast, the number of low-amplitude events (−0.25σ anomalies in ) that occur during a given 500-yr period is a positive integer. Significance is therefore assigned by tallying the number of LIM realizations with as many or more low-amplitude events over any continuous 500-yr period. Finally, if the shift in the 500-yr mean of the MCA was generated merely by chance (following our null hypothesis), then its sign is necessarily arbitrary. Significance is therefore evaluated against the absolute value of all differences between the first and last 500 years of all realizations.

4. Results

Using an LIM that represents global SST variability and drought in the western United States, we generated 1000 synthetic PDSI and SST realizations that each span one millennium. The LIM reproduces realistic spatial correlation patterns between the tropical Pacific and WUS hydroclimate (Fig. 3). Correlations patterns between the January–March (JFM) Niño-3.4 time series (average SSTs over 5°S–5°N, 170°–120°W) and June–August (JJA) scPDSIpm (Sheffield et al. 2006) are very similar in both the observations and the LIM output (Fig. 3).

Fig. 3.

Correlation maps between boreal winter (JFM) Niño-3.4 and summer (JJA) PDSI calculated from observations (shown at top left) and the first 11 LIM realizations. The iteration number of each LIM realization is shown at the bottom left of each panel.

Fig. 3.

Correlation maps between boreal winter (JFM) Niño-3.4 and summer (JJA) PDSI calculated from observations (shown at top left) and the first 11 LIM realizations. The iteration number of each LIM realization is shown at the bottom left of each panel.

Power spectra of the Niño-3.4 time series generated by the LIM exhibit a broad spectral peak over the interannual time scales associated with ENSO (gray shading, Fig. 4). Indeed, the range of LIM power spectra encompasses three different observational versions of Niño-3.4 (Fig. 4): Kaplan SST (Kaplan et al. 1998), ERSSTv3b (Smith et al. 2008), and HadSST1 (Rayner et al. 2003), which verifies that this version of the LIM generates realistic interannual variability [a result also discussed in Newman (2013) and Penland and Sardeshmukh (1995), among others]. Although one observational product (ERSSTv3b) exhibits spectral densities above the LIM at near-centennial time horizons, this low-frequency variability probably reflects long-term trends in the dataset that are not robust across the tropical Pacific (e.g., Solomon and Newman 2011; Emile-Geay et al. 2013a,b).

Fig. 4.

Power spectra computed from Niño-3.4 time series. The gray shading encompasses the lower 2.5% and 97.5% quantiles of Niño-3.4 time series simulated by the LIM (e.g., 95% of the distribution of all realizations at all frequencies). The three thin lines are spectra computed from Niño-3.4 time series derived from the following gridded observational data products: Kaplan SST (Kaplan et al. 1998), ERSSTv3v (Smith et al. 2008), and HadSST1 (Rayner et al. 2003). Vertical dashed lines bracket the 2–7-yr ENSO band.

Fig. 4.

Power spectra computed from Niño-3.4 time series. The gray shading encompasses the lower 2.5% and 97.5% quantiles of Niño-3.4 time series simulated by the LIM (e.g., 95% of the distribution of all realizations at all frequencies). The three thin lines are spectra computed from Niño-3.4 time series derived from the following gridded observational data products: Kaplan SST (Kaplan et al. 1998), ERSSTv3v (Smith et al. 2008), and HadSST1 (Rayner et al. 2003). Vertical dashed lines bracket the 2–7-yr ENSO band.

To help visualize the output generated by our LIM of global SSTs and western hydroclimate, we illustrate the timing and spatial structure of megadroughts in just one LIM realization (Figs. 5 and 6). In this run, the maximum WUS DAI35 and minimum SWUS values coincide, but this does not occur in all realizations, and high DAI35 values are not always driven by the same SST patterns as low SWUS values. During the first widespread event (culminating at year 54), DAI35 values reached 27.1% (Fig. 5a). Other comparable periods are also present in this run (e.g., years 230–270, 310–390, and 950–970). Here we focus on the year 54 event (vertical dashed line) to distinguish its SST pattern from that of the most severe event (Fig. 6). SST anomalies were cool almost everywhere, and most notably in the northern Atlantic (Fig. 5b). Negative values are apparent in much of the central portion of the domain, with wet anomalies occurring in the Northwest and along much of the Pacific coast (Fig. 5c). The most severe event in the Southwest reached a value of (red circle in Fig. 6; LIM year 946). Average SST anomalies in the LIM during that time were close to neutral in the Indian Ocean, cool in the tropical Pacific, and warm in the western Pacific and North Atlantic (Fig. 6b). Negative PDSI values (between −1 and −2) are found throughout most of the entire Southwest, and most especially in the central part of the domain. Modest positive values are present across the northern regions of the WUS (Fig. 6c).

Fig. 5.

Example output from one LIM realization. (a) DAI35 from this run, with the vertical dashed line marking the first widespread event during the realization (year 54). Also mapped are average (b) SST and (c) PDSI conditions for this event.

Fig. 5.

Example output from one LIM realization. (a) DAI35 from this run, with the vertical dashed line marking the first widespread event during the realization (year 54). Also mapped are average (b) SST and (c) PDSI conditions for this event.

Fig. 6.

Example output from the same LIM realization used in Fig. 5. (a) Southwestern U.S. from this run, with the red circle marking the most severe event (year 929) during the realization. Also mapped are average (b) SST and (c) PDSI anomalies during this drought.

Fig. 6.

Example output from the same LIM realization used in Fig. 5. (a) Southwestern U.S. from this run, with the red circle marking the most severe event (year 929) during the realization. Also mapped are average (b) SST and (c) PDSI anomalies during this drought.

DAI35 time series from individual LIM realizations encompass a similar range as the tree-ring-based reconstructions (Fig. 7). Maps of the most widespread events (maximum DAI35) likewise depict prolonged events commensurate with observations in their spatial scale and severity (Fig. 8). Several of these realizations have at least one event with DAI35 values greater than or equal to the maximum reconstructed value (40.8% of WUS in the NADA). There are also instances of widespread megadroughts that are even more spatially extensive than have been observed (e.g., run 2). In total, 35.4% of the runs included at least one event as widespread as the 1130–1164 megadrought in the NADA.

Fig. 7.

DAI35 calculated from a subset of 19 LIM realizations. The top left panel shows values computed from the NADA (Cook et al. 2004), while the remaining panels are each computed from a single 1000-yr-long realization. The dashed line on each panel marks the highest DAI35 value in the NADA for reference. The realization number is given on the upper left of each panel.

Fig. 7.

DAI35 calculated from a subset of 19 LIM realizations. The top left panel shows values computed from the NADA (Cook et al. 2004), while the remaining panels are each computed from a single 1000-yr-long realization. The dashed line on each panel marks the highest DAI35 value in the NADA for reference. The realization number is given on the upper left of each panel.

Fig. 8.

Reconstructed PDSI during the 1130–1164 CE period (shown at top left) shown with average PDSI from a subset of 19 LIM realizations (the same runs as in Fig. 7) over the 35-yr period leading to maximum DAI35 values during each run.

Fig. 8.

Reconstructed PDSI during the 1130–1164 CE period (shown at top left) shown with average PDSI from a subset of 19 LIM realizations (the same runs as in Fig. 7) over the 35-yr period leading to maximum DAI35 values during each run.

Average SST anomalies from the 35 years associated with each run’s maximum DAI35 are reminiscent of the canonical patterns associated with La Niña and the “cold phase” of the PDO with cool SSTs in the Indian Ocean, equatorial Pacific, and along the western edges of North and South America (Fig. 9). Warm SST anomalies are found in the northwestern Pacific and North Atlantic. The amplitudes of the averaged anomalies, however, are much weaker than during La Niña or PDO cold phases (individual events can have stronger anomalies relative to composites). The westward extent of the cold anomaly along the equator is also less zonal than it is during normal La Niña events.

Fig. 9.

Average SST anomalies from all LIM realizations during maximum DAI35 values from each run.

Fig. 9.

Average SST anomalies from all LIM realizations during maximum DAI35 values from each run.

Global average SST anomalies associated with the lowest values of from each run (Fig. 10) differ from those connected with maximum DAI35 (Fig. 9). In particular, cool SST anomalies along the equatorial Pacific extend less westward in the maximum DAI35 pattern than they do at times of minimum (cf. Figs. 9 and 10). Moreover, SST anomalies associated with SWUS megadrought have a more distinctively La Niña–like pattern, as compared with the anomalies associated with widespread drought.

Fig. 10.

Average SST anomalies from all LIM realizations during the most severe megadroughts (35-yr average PDSI values below ) from each run.

Fig. 10.

Average SST anomalies from all LIM realizations during the most severe megadroughts (35-yr average PDSI values below ) from each run.

Distributions of our test statistic for widespread WUS megadrought from LIM output establish confidence limits for various reconstructed DAI35 thresholds (Fig. 11). In this figure, each panel shows the LIM-based PDF of time spent at or above each DAI35 threshold in 5% increments. The vertical dashed line marks the fraction of time in the NADA spent at or above these thresholds; the gray shaded area of each PDF identifies the portion of LIM realizations that have a lower fraction of years than the NADA at each threshold. For example, in the NADA 1072 out of 1200 years saw DAI35 values of at least 5% ( = 0.893). Of the 1000 LIM realizations, 81.3% simulated a smaller fraction of time below this threshold than in the NADA, meaning that the probability (p value) that the observed value at least as large as was drawn from the null distribution is 0.187 (p = 1 − 813/1000).

Fig. 11.

LIM-based empirical distributions of widespread megadrought for various DAI35 thresholds. In each panel, the ordinate is the fraction of all years spent below a given DAI35 threshold (marked at the top of each panel in bold). The shaded portion of the PDF is the fraction of LIM realizations with a proportionally lower number of years spent below each threshold than the NADA, while the white shading is the remaining portion of the PDF. The vertical dotted line is fraction of time (f) in the NADA spent at or below each DAI35 threshold. The y axis is normalized so that the number of realizations (N) at each value of f is divided by the total number of realizations (r).

Fig. 11.

LIM-based empirical distributions of widespread megadrought for various DAI35 thresholds. In each panel, the ordinate is the fraction of all years spent below a given DAI35 threshold (marked at the top of each panel in bold). The shaded portion of the PDF is the fraction of LIM realizations with a proportionally lower number of years spent below each threshold than the NADA, while the white shading is the remaining portion of the PDF. The vertical dotted line is fraction of time (f) in the NADA spent at or below each DAI35 threshold. The y axis is normalized so that the number of realizations (N) at each value of f is divided by the total number of realizations (r).

The shape of the distributions for is clearly more skewed for higher DAI35 thresholds because most of the realizations only simulate a small fraction of years with DAI35 values greater than 35% (Fig. 11). The LIM also produces a small number of events (not shown) with even greater values (i.e., the DAI35 thresholds of 45% and 50% both include a few nonzero values when considering the full 1000-member ensemble).

None of the reconstructed DAI35 values are significant at the 95% confidence limit (Fig. 12). In this case, we applied a one-sided significance test because values of can only be positive, and we are interested in the fraction of LIM realizations with more time spent below various DAI35 thresholds. A significance level value of 0.3 (y axis of Fig. 12) means that 3/10 of all LIM realizations simulated less time at a given DAI35 level than in the NADA, while 7/10 spent more time at that threshold. Even the most widespread reconstructed event (1130–1164 CE), when approximately 40.8% of the WUS was affected by megadrought, is not significant at the 90% confidence limit (p value = 0.354).

Fig. 12.

Empirical confidence limits for rejecting the null hypothesis that widespread megadroughts of the last millennium are caused by SST variability, combined with spatial and temporal autocorrelation. For each DAI35 threshold (abscissa), the fraction of years in the NADA found at that value is compared against the total fraction of years in all LIM realizations. The line with open circles traces the significance limits for the amount of time spent at each DAI35 threshold. The 50% (black), 66% (dashed), 90% (dash-dotted), and 95% (dotted) confidence limits are also shown for reference.

Fig. 12.

Empirical confidence limits for rejecting the null hypothesis that widespread megadroughts of the last millennium are caused by SST variability, combined with spatial and temporal autocorrelation. For each DAI35 threshold (abscissa), the fraction of years in the NADA found at that value is compared against the total fraction of years in all LIM realizations. The line with open circles traces the significance limits for the amount of time spent at each DAI35 threshold. The 50% (black), 66% (dashed), 90% (dash-dotted), and 95% (dotted) confidence limits are also shown for reference.

We assign significance to the other test statistics , , and as follows:

  • Magnitude. During the worst SWUS drought, tree-ring reconstructed reached −0.62σ ( = −0.62σ). Out of the 1000 LIM realizations, 34.8% of the runs included events as severe or worse (more negative) than this period (Fig. 13), while the other runs did not include such high magnitude events. The null hypothesis can only be rejected at about the 65% confidence limit (p value = 0.348) for the magnitude of the worst megadroughts in the SWUS over the last 1200 years.

  • Clustering. Five low-amplitude megadroughts occurred during the MCA (Fig. 1 in Coats et al. 2016b), which is more than occurred during most 500-yr periods of the LIM (Fig. 14). In total, 2.1% of LIM runs had 500-yr periods with as many or more megadroughts, suggesting that the null hypothesis for megadrought clusters can be rejected at the 90% confidence level (p value = 0.021).

  • Mean shift. The long-term difference in the MCA mean versus the 1500–2000 CE mean (0.47 PDSI units) is highly significant against the null hypothesis (Fig. 15). Out of the 1000 LIM realizations, none had higher magnitude 500-yr mean shifts, hence the reconstructed values are significant at the 99.9% confidence limit or higher (p value < 0.001).

Fig. 13.

PDF of normalized minimum , averaged over the SWUS. Black vertical bars indicate realizations with runs that produced minimum values below the NADA (−0.63σ).

Fig. 13.

PDF of normalized minimum , averaged over the SWUS. Black vertical bars indicate realizations with runs that produced minimum values below the NADA (−0.63σ).

Fig. 14.

Distribution of maximum number of low-amplitude megadroughts to have occurred over a 500-yr period in each LIM realization. Black vertical bars identify the fraction of runs with 500 year periods that include an equal or greater number of events than occurred during the MCA (801–1300 CE), as inferred from tree-ring reconstructions.

Fig. 14.

Distribution of maximum number of low-amplitude megadroughts to have occurred over a 500-yr period in each LIM realization. Black vertical bars identify the fraction of runs with 500 year periods that include an equal or greater number of events than occurred during the MCA (801–1300 CE), as inferred from tree-ring reconstructions.

Fig. 15.

PDF of SWUS shifts in 500-yr means, expressed in PDSI unit differences . The bars are white because all runs had smaller amplitude differences than inferred from the NADA for the SWUS between 801–1300 CE and 1500–2000 CE (this difference is marked by an asterisk on the x axis).

Fig. 15.

PDF of SWUS shifts in 500-yr means, expressed in PDSI unit differences . The bars are white because all runs had smaller amplitude differences than inferred from the NADA for the SWUS between 801–1300 CE and 1500–2000 CE (this difference is marked by an asterisk on the x axis).

Our results were robust with respect to a number of other alternative assumptions made as part of constructing of the LIM, identifying megadrought, and evaluating significance. For example, other calibration periods (e.g., 1980–2010) yielded structurally similar linear deterministic matrices, meaning that estimating is not particularly sensitive to the 1950–2000 interval used here. Likewise, the scPDSIpm product we included in the LIM uses precipitation from the CRU TS3.10 data product [for details, see Sheffield et al. (2006)]. There are 15 alternative PDSI products, however, including ones with different underlying precipitation inputs and approaches to estimating evapotranspiration (Sheffield et al. 2006). Use of the other 15 versions of PDSI did not appreciably change the structure of the LIM or the confidence limits we established. Furthermore, the scPDSIpmdiffers from its predecessor, the PDSI, in that the local surface storage coefficient is estimated empirically from data, rather than prescribed as a spatially invariant constant. Use of the original PDSI did not critically alter the nature of variability in the LIM, nor did use of alternative SST datasets (e.g., ERSSTv3b; Smith et al. 2008).

There was not any strong dependence of the LIM on the details of the detrending; use of a spline fit, high-pass filter, or linear fit all yielded structurally similar LIMs. We did find, however, that failing to detrend the data had a modest impact on the significance levels, in that it made more prolonged and more widespread events somewhat more common in the LIM. The likely explanation for this particular sensitivity is that the global warming pattern becomes included in the LIM as a “least damped” mode that grows and decays on multidecadal time scales, imparting additional low-frequency variability. The limitation of detrending is, of course, that there is some possibility that we have inadvertently removed unforced low-frequency climate variability. If this is the case, rejecting the null hypothesis would be even more difficult, as it would raise the confidence limits even further. Our results should therefore be interpreted as being slightly “conservative” because of our use of detrended data in constructing the LIM. Note, however, that the tree-ring data are also detrended, so there is also a possibility that the true amplitude of climate fluctuations is greater in reality than inferred from the reconstruction (see discussion). In either case, the results presented here are at least consistent, in that both the LIM and the tree-ring data have had long-term trends removed.

5. Discussion

Several key features of Pacific SST variability and its impact on WUS hydroclimate can be reproduced by the relatively simple, low-order linear dynamical system represented by the LIM. For example, the correlation patterns between the Niño-3.4 time series and JJA PDSI are nearly identical between observations and the LIM (Fig. 3). Likewise, the LIM simulates spectral densities in the central Pacific Ocean that encompass the range of observational estimates (Fig. 4). Both of these features are important for verifying that the LIM’s underlying statistics in space and time are representative of observations over the short instrumental period, and hence relevant to the longer time scales of megadrought.

The LIM generates synthetic megadroughts, but this outcome was not guaranteed from its formulation because we used only seasonal data from the late twentieth century (an interval when such events have not been observed). We also emphasize the LIM did not include any information from tree-ring reconstructions of past droughts in the western United States. Megadroughts in the LIM are therefore to be regarded as an emergent feature of the low-order dynamical system we have investigated. This finding could be further used as a basis for predicting an optimal forcing pattern (or combinations of patterns) that would elevate the risk of widespread megadrought if it were to appear across ocean basins this century.

We used the LIM to evaluate our null hypothesis, namely that widespread WUS megadroughts occur fundamentally from stochastic forcing. To test the significance of reconstructed events against this null hypothesis, we compared the fraction of years spent below each DAI35 threshold (; see section 3e) in the NADA to the LIM’s distribution at those same levels. This comparison, in turn, permits us to establish confidence limits for rejecting the null hypothesis across reconstructed DAI35 thresholds. For instance, the maximum DAI35 value of the NADA is 40.8%, and just 3 years over the last 1200 years were spent at or above the 40% threshold. The fraction of time megadrought conditions prevailed in the NADA at this level is therefore 3/1200 (0.0025). While 64.6% of the LIM realizations did not yield DAI35 values this high, 35.4% of the runs include megadroughts that are equally (or more) widespread (Fig. 12). We may therefore only reject our null hypothesis at the nominal 64% confidence level. Similar considerations apply to other lower thresholds, none of which may be said to be significant at the 95% confidence limit.

As with the spatial scale of reconstructed WUS droughts, the null hypothesis cannot be rejected with very much confidence for the magnitude of the worst droughts in the SWUS. That is, about a third (34.8%) of the LIM realizations generated events as severe as, or worse than, the worst drought in the SWUS. Hence, such events should be considered relatively normal occurrences in the hydroclimatic statistics of the region on millennial time scales.

Two features of SWUS megadrought are more significant than the scale and severity of past events. First, the clustering of low-amplitude megadroughts during the MCA is unlikely to have occurred from chance alone (p value = 0.021), and hence we reject the null hypothesis for this characteristic of SWUS hydroclimate with fairly high (>95%) confidence. We can be even more confident in rejecting the null hypothesis for the shift in 500-yr mean reconstructed PDSI values (p value < 0.001), implying that such long-term changes are extremely unlikely to be driven by chance alone.

What does accepting the null hypothesis for two of the four megadrought test statistics, and rejecting it for the other two, mean for the nature of hydroclimate in the WUS and SWUS? First and foremost, our results help clarify what aspects of megadrought occurrence in the past are truly exceptional given what we can observe in the modern system. Clusters of megadroughts, and especially changes in long-term means, are unlikely to be generated by stationary statistics and stochastic forcing. Their existence in the NADA therefore requires further examination. One possibility is that their odds of occurrence were externally forced during the MCA [e.g., as discussed extensively in Cook et al. (2016)], either by slightly enhanced summer insolation (Lean 2000) or by coupled feedbacks between the land surface, vegetation, and atmosphere (Cook et al. 2011). Alternatively, a multicentury warm phase of the Atlantic multidecadal oscillation (AMO) (Coats et al. 2016b; Wang et al. 2017) could have shifted the mean state of hydroclimate in the SWUS and increased megadrought risk (Feng et al. 2008; Oglesby et al. 2012). Such persistence itself would be unusual in the context of seasonal dynamics of global SST over the late twentieth century, and could require either ocean–atmosphere coupling as originally proposed by Bjerknes (1969) or the influences of external forcings [including “quiet” periods of relatively low variability in forcings as described by Bradley et al. (2016)]. It is clear that a smaller fraction of total SST variance is captured by the LIM in the Atlantic than in the Pacific (Fig. 2). Consequently we could be undersampling AMO variability in construction of the LIM due to the short time scale of the observational record. If the AMO is responsible for megadrought clustering during the MCA, then we cannot completely rule out the possibility of megadroughts being internally generated because our LIM does not fully characterize variability in the North Atlantic.

The change in the long-term mean of reconstructed PDSI is the most statistically significant feature of SWUS hydroclimate. If it is climatic in origin (Cook et al. 2004; Coats et al. 2016b), it too likely arises from either nonlinear or externally forced sources of variability. If it is not driven by climate, then it could arise from phenomena related to the proxy itself (e.g., biological growth trends not completely removed through detrending; Fritts 1976) or the methodology used to develop the reconstruction. For example, reconstructed PDSI values of the NADA are associated with calibration and measurement uncertainties, and the effect of these uncertainties on our rejection or acceptance of the null hypothesis cannot be directly addressed using the NADA. In contrast, more recent drought atlases (Cook et al. 2010, 2015) have been constructed using methods that explicitly generate an ensemble of probable reconstructions over the respective domains, and a similar approach could be incorporated to derive an ensemble of WUS reconstructions to compare against the LIM. Further work should therefore test ensembles of drought atlas reconstructions against our LIM-based significance levels, although we have focused here exclusively on the available expected value reconstructions over the NADA domain.

Notwithstanding the inherent methodological limitations of reconstructing past climates, the differences in the multicentury mean PDSI values between the modern period and the MCA are extremely unlikely to have arisen from stationary statistics. The LIM therefore provides important insight into the details of the data that require further investigation even if they are related to the proxy or reconstruction methodology.

Under the changing climates of the present day, it is often stated that individual extreme weather events cannot be unambiguously attributed to climate change, but that the occurrence of more frequent extremes reflects the net effects of anthropogenic activity on the statistics of the climate system (e.g., Trenberth et al. 2015). A similar principle applies to megadroughts during the MCA: such events should occur naturally even without any changes to the boundary conditions (e.g., solar, orbital, or land cover) of WUS hydroclimate. However, their frequency from 801 to 1300 CE suggests that exogenous factors may have played a role in making megadroughts more probable during that time.

Our results suggest that the LIM is a viable tool to establish baselines (or expectations) when evaluating long model simulations or suites of paleoclimate reconstructions. For example, it is not entirely clear if GCMs simulate megadroughts at realistic rates with or without externally varying boundary conditions (e.g., Ault et al. 2014; Coats et al. 2013b). A LIM could be trained on model output, and the PDF of megadrought occurrence derived from it could be compared to the observation-based LIM we used. Such an analysis would therefore not be limited by the size of millennial-scale GCM ensembles because the LIM would empirically fill out the PDF of megadrought from a limited set of runs. Likewise, an LIM could be adapted to simulate the risk of widespread megadrought, building on the work of Ault et al. (2014) and Cook et al. (2015), both of which focused on single time series methods of estimating risk. Finally, and perhaps most importantly, this effort should motivate wider adoption of the LIM framework in paleoclimate significance testing. In high-resolution paleoclimatology, researchers frequently observe low-frequency variability in time series or networks of paleoclimate archives, and it is often unclear whether such fluctuations require external forcings (e.g., volcanic, solar, or land surface processes), “out of sample” climate mechanisms, or if they are simply consistent with the variability observed during the late twentieth century. Adapting an LIM to these settings would shed new light on the importance of internal variability and forced climate change across a broad range of reconstructed phenomena.

6. Conclusions and implications

We constructed a simple, low-order empirical model (an LIM) of SST and WUS hydroclimate to evaluate the odds that megadroughts occur by chance—that is, arising from observable spatiotemporal modes of late-twentieth-century climate variability, stochastic forcing, and sufficient time. Despite being trained only on seasonal data from the late twentieth century, the LIM successfully reproduces realistic power spectra in the tropical Pacific and correlation patterns between the Niño-3.4 region and WUS hydroclimate. It also produces megadroughts that are comparable in their duration, spatial scale, and magnitude as the most severe events of the last 12 centuries (801–2000 CE). For example, the most widespread megadrought in the NADA lasted from between 1130 and 1164 CE and reached DAI35 values of approximately 40%. About 35% of LIM realizations included events of this scale or worse. Because the LIM is an empirical model of dynamics with fundamentally stationary statistics, we cannot reject the null hypothesis that even the most widespread megadroughts are internally generated. Similar considerations apply to other smaller DAI35 thresholds, and to the magnitude of the most severe SWUS megadroughts.

In contrast to the spatial scale and magnitude of megadroughts, the clustering of these anomalies and the change in the 500-yr mean since the MCA are nominally, and highly, significant against the null hypothesis, respectively. Accordingly, these characteristics of the NADA may require nonlinear processes, external forcing mechanisms, or exceptionally persistent climate states that have not been observed during the twentieth century [e.g., a multicentury warm phase of the AMO, as discussed by Coats et al. (2016b)]. These results help clarify which aspects of megadrought are unusual but unremarkable in the context of the observable system (widespread intervals and high magnitude events), and what characteristics are truly unusual and require further investigation (clustering and the lower mean of the MCA).

In summary:

  • The most widespread megadrought in the last millennium occurred between 1130 and 1164 CE and encompassed 40% of the WUS domain. About 35% of the LIM output simulates megadroughts that are equally (or more) widespread (Fig. 12), meaning we cannot reject the null hypothesis of internal variability as causative at the 90% or 95% confidence limit (p value = 0.646).

  • The worst drought in the reconstructed SWUS PDSI average was −0.62σ. Because 34.8% of LIM runs had values below this level, we cannot reject the null hypothesis for the severity of the worst megadrought in the NADA (p value = 0.652).

  • During the MCA five low-amplitude megadroughts ( −0.25σ) occurred over a 500-yr period; only 2.1% of LIM runs had 500-yr periods with as many or more megadroughts, suggesting that we can reject the null hypothesis at the 95% confidence limit (p value = 0.021).

  • Mean PDSI values from 801 to 1300 CE were 0.47 PDSI units lower than from 1500 to 2000 CE. All of the LIM runs had smaller mean shifts between their first and last 500 years, meaning we can reject the null hypothesis with a high degree of confidence (p value < 0.001).

Finally, our results hold a number of key implications for stakeholder and research communities alike. First, events as arid as and even more arid than the most widespread drought of the last 1200 years are possible, even without any changes to the boundary conditions. Given their spatial extent, these events would require large-scale efforts to manage water across much of the western United States, as well as to protect ecosystems and agriculture. A LIM could be used to support water resource management strategies because it permits us to generate a large number of statistically plausible outcomes even without climate change. These scenarios could, in turn, be used by policy makers and water resource managers to prepare for possible future water stress. Second, given that the LIM simulates events that are slightly worse than those in the NADA, and given that the NADA is based on tree-ring data that are themselves detrended, prewhitened, and then “re-reddened” [see the supplementary material in Cook et al. (2004)], it is possible that there is more low-frequency variability inherent to the climate system than would be inferred from the NADA alone [see also Ault et al. (2013a)]. Even if this is the case, our results clearly show that we cannot rule out the possibility that such variability is internally generated by the climate system itself. Future widespread megadrought risks are therefore likely to emerge as a combination of this internal low-frequency variability and long-term forced changes.

Acknowledgments

This material is based on work partially supported by NSF EaSM2 (Earth System Model 2) Grants AGS1243125 and AGS1243204 and NSF Grants AGS1602564, AGS1401400, and AGS1602512. BI Cook is supported by the NASA Modeling, Analysis, and Prediction program. LDEO Contribution No. 8165.

REFERENCES

REFERENCES
Acuna-Soto
,
R.
,
D. W.
Stahle
,
M. K.
Cleaveland
, and
M. D.
Therrell
,
2002
:
Megadrought and megadeath in 16th century Mexico
.
Emerging Infect. Dis.
,
8
,
360
362
, doi:.
Alexander
,
M. A.
,
L.
Matrosova
,
C.
Penland
,
J. D.
Scott
, and
P.
Chang
,
2008
:
Forecasting Pacific SSTs: Linear inverse model predictions of the PDO
.
J. Climate
,
21
,
385
402
, doi:.
Ault
,
T. R.
,
J. E.
Cole
, and
S.
St. George
,
2012
:
The amplitude of decadal to multidecadal variability in precipitation simulated by state-of-the-art climate models
.
Geophys. Res. Lett.
,
39
,
L21705
, doi:.
Ault
,
T. R.
,
J. E.
Cole
,
J. T.
Overpeck
,
G. T.
Pederson
,
S. S.
George
,
B.
Otto-Bliesner
,
C. A.
Woodhouse
, and
C.
Deser
,
2013a
:
The continuum of hydroclimate variability in western North America during the last millennium
.
J. Climate
,
26
,
5863
5878
, doi:.
Ault
,
T. R.
,
C.
Deser
,
M.
Newman
, and
J.
Emile-Geay
,
2013b
:
Characterizing decadal to centennial variability in the equatorial Pacific during the last millennium
.
Geophys. Res. Lett.
,
40
,
3450
3456
, doi:.
Ault
,
T. R.
,
J. E.
Cole
,
J. T.
Overpeck
,
G. T.
Pederson
, and
D. M.
Meko
,
2014
:
Assessing the risk of persistent drought using climate model simulations and paleoclimate data
.
J. Climate
,
27
,
7529
7549
, doi:.
Ault
,
T. R.
,
J.
Mankin
,
B. I.
Cook
, and
J. E.
Smerdon
,
2016
:
Relative impacts of mitigation, temperature, and precipitation on 21st-century megadrought risk in the American Southwest
.
Sci. Adv.
,
2
,
e1600873
, doi:.
Bartlett
,
M. S.
,
1978
: An Introduction to Stochastic Processes: With Special Reference to Methods and Applications. Cambridge University Press, 388 pp.
Benson
,
L.
, and Coauthors
,
2002
:
Holocene multidecadal and multicentennial droughts affecting Northern California and Nevada
.
Quat. Sci. Rev.
,
21
,
659
682
, doi:.
Benson
,
L.
,
B.
Linsley
,
J.
Smoot
,
S.
Mensing
,
S.
Lund
,
S.
Stine
, and
A.
Sarna-Wojcicki
,
2003
:
Influence of the Pacific decadal oscillation on the climate of the Sierra Nevada, California and Nevada
.
Quat. Res.
,
59
,
151
159
, doi:.
Bjerknes
,
J.
,
1969
:
Atmospheric teleconnections from the equatorial Pacific
.
Mon. Wea. Rev.
,
97
,
163
172
, doi:.
Bradley
,
R. S.
,
H.
Wanner
, and
H. F.
Diaz
,
2016
:
The medieval quiet period
.
Holocene
,
26
,
990
993
, doi:.
Brown
,
D. P.
, and
A. C.
Comrie
,
2004
:
A winter precipitation ‘dipole’ in the western United States associated with multidecadal ENSO variability
.
Geophys. Res. Lett.
,
31
,
L09203
, doi:.
Buckley
,
B. M.
, and Coauthors
,
2010
:
Climate as a contributing factor in the demise of Angkor, Cambodia
.
Proc. Natl. Acad. Sci. USA
,
107
,
6748
6752
, doi:.
Burgman
,
R.
,
R.
Seager
,
A.
Clement
, and
C.
Herweijer
,
2010
:
Role of tropical Pacific SSTs in global medieval hydroclimate: A modeling study
.
Geophys. Res. Lett.
,
37
,
L06705
, doi:.
Clement
,
A. C.
,
R.
Seager
,
M. A.
Cane
, and
S. E.
Zebiak
,
1996
:
An ocean dynamical thermostat
.
J. Climate
,
9
,
2190
2196
, doi:.
Coats
,
S.
,
J. E.
Smerdon
,
B. I.
Cook
, and
R.
Seager
,
2013a
:
Stationarity of the tropical Pacific teleconnection to North America in CMIP5/PMIP3 model simulations
.
Geophys. Res. Lett.
,
40
,
4927
4932
, doi:.
Coats
,
S.
,
J. E.
Smerdon
,
R.
Seager
,
B. I.
Cook
, and
J. F.
González-Rouco
,
2013b
:
Megadroughts in southwestern North America in ECHO-G millennial simulations and their comparison to proxy drought reconstructions
.
J. Climate
,
26
,
7635
7649
, doi:.
Coats
,
S.
,
B. I.
Cook
,
J. E.
Smerdon
, and
R.
Seager
,
2015a
:
North American pancontinental droughts in model simulations of the last millennium
.
J. Climate
,
28
,
2025
2043
, doi:.
Coats
,
S.
,
J. E.
Smerdon
,
B. I.
Cook
, and
R.
Seager
,
2015b
:
Are simulated megadroughts in the North American Southwest forced?
J. Climate
,
28
,
124
142
, doi:.
Coats
,
S.
,
J. E.
Smerdon
,
B. I.
Cook
,
R.
Seager
,
E. R.
Cook
, and
K. J.
Anchukaitis
,
2016a
:
Internal ocean–atmosphere variability drives megadroughts in western North America
.
Geophys. Res. Lett.
,
43
,
9886
9894
, doi:.
Coats
,
S.
,
J. E.
Smerdon
,
K. B.
Karnauskas
, and
R.
Seager
,
2016b
:
The improbable but unexceptional occurrence of megadrought clustering in the American West during the Medieval Climate Anomaly
.
Environ. Res. Lett.
,
11
,
074025
, doi:.
Cook
,
B. I.
,
R.
Seager
, and
R. L.
Miller
,
2011
:
The impact of devegetated dune fields on North American climate during the late Medieval Climate Anomaly
.
Geophys. Res. Lett.
,
38
,
L14704
, doi:.
Cook
,
B. I.
,
J. E.
Smerdon
,
R.
Seager
, and
E. R.
Cook
,
2014
:
Pan-continental droughts in North America over the last millennium
.
J. Climate
,
27
,
383
397
, doi:.
Cook
,
B. I.
,
T. R.
Ault
, and
J. E.
Smerdon
,
2015
:
Unprecedented 21st century drought risk in the American Southwest and central plains
.
Sci. Adv.
,
1
,
e1400082
, doi:.
Cook
,
B. I.
,
E. R.
Cook
,
J. E.
Smerdon
,
R.
Seager
,
A. P.
Williams
,
S.
Coats
,
D. W.
Stahle
, and
J. V.
Díaz
,
2016
:
North American megadroughts in the Common Era: Reconstructions and simulations
.
Wiley Interdiscip. Rev.: Climate Change
,
7
,
411
432
, doi:.
Cook
,
B. I.
,
A. P.
Williams
,
J. S.
Mankin
,
R.
Seager
,
J. E.
Smerdon
, and
D.
Singh
,
2017
:
Revisiting the leading drivers of Pacific coastal drought variability in the contiguous United States
.
J. Climate
, https://doi.org/10.1175/JCLI-D-17-0172.1, in press.
Cook
,
E. R.
,
D. M.
Meko
,
D. W.
Stahle
, and
M. K.
Cleaveland
,
1999
:
Drought reconstructions for the continental United States
.
J. Climate
,
12
,
1145
1162
, doi:.
Cook
,
E. R.
,
C. A.
Woodhouse
,
C. M.
Eakin
,
D. M.
Meko
, and
D. W.
Stahle
,
2004
:
Long-term aridity changes in the western United States
.
Science
,
306
,
1015
1018
, doi:.
Cook
,
E. R.
,
R.
Seager
,
M. A.
Cane
, and
D. W.
Stahle
,
2007
:
North American drought: Reconstructions, causes, and consequences
.
Earth-Sci. Rev.
,
81
,
93
134
, doi:.
Cook
,
E. R.
,
K. J.
Anchukaitis
,
B. M.
Buckley
,
R. D.
D’Arrigo
,
G. C.
Jacoby
, and
W. E.
Wright
,
2010
:
Asian monsoon failure and megadrought during the last millennium
.
Science
,
328
,
486
489
, https://doi.org/10.1126/science.1185188.
Dai
,
A.
,
2011
:
Drought under global warming: A review
.
Wiley Interdiscip. Rev.: Climate Change
,
2
,
45
65
, doi:.
Emile-Geay
,
J.
,
K. M.
Cobb
,
M. E.
Mann
, and
A. T.
Wittenberg
,
2013a
:
Estimating central equatorial Pacific SST variability over the past millennium. Part I: Methodology and validation
.
J. Climate
,
26
,
2302
2328
, doi:.
Emile-Geay
,
J.
,
K. M.
Cobb
,
M. E.
Mann
, and
A. T.
Wittenberg
,
2013b
:
Estimating central equatorial Pacific SST variability over the past millennium. Part II: Reconstructions and implications
.
J. Climate
,
26
,
2329
2352
, doi:.
Feng
,
S.
,
R. J.
Oglesby
,
C. M.
Rowe
,
D. B.
Loope
, and
Q.
Hu
,
2008
:
Atlantic and Pacific SST influences on medieval drought in North America simulated by the Community Atmospheric Model
.
J. Geophys. Res.
,
113
,
D11101
, doi:.
Flato
,
G.
, and Coauthors
,
2013
: Evaluation of climate models. Climate Change 2013: The Physical Science Basis. T. F. Stocker et al., Eds., Cambridge University Press, 741–866.
Fritts
,
H. C.
,
1976
: Tree Rings and Climate. Academic Press, 567 pp.
Graham
,
N. E.
, and Coauthors
,
2007
:
Tropical Pacific–mid-latitude teleconnections in medieval times
.
Climatic Change
,
83
,
241
285
, doi:.
Graham
,
N. E.
,
C. M.
Ammann
,
D.
Fleitmann
,
K. M.
Cobb
, and
J.
Luterbacher
,
2011
:
Support for global climate reorganization during the “Medieval Climate Anomaly.”
Climate Dyn.
,
37
,
1217
1245
, doi:.
Guilyardi
,
E.
,
A.
Wittenberg
,
A.
Fedorov
,
M.
Collins
,
C.
Wang
,
A.
Capotondi
,
G. J.
van Oldenborgh
, and
T.
Stockdale
,
2009
:
Understanding El Niño in ocean–atmosphere general circulation models: Progress and challenges
.
Bull. Amer. Meteor. Soc.
,
90
,
325
340
, doi:.
Hodell
,
D.
,
J.
Curtis
, and
M.
Brenner
,
1995
:
Possible role of climate in the collapse of classic Maya civilization
.
Nature
,
375
,
391
394
, doi:.
Hunt
,
B. G.
,
2006
:
The Medieval Warm Period, the Little Ice Age and simulated climatic variability
.
Climate Dyn.
,
27
,
677
694
, doi:.
Kaplan
,
A.
,
M.
Cane
,
Y.
Kushnir
,
A.
Clement
,
M.
Blumenthal
, and
B.
Rajagopalan
,
1998
:
Analyses of global sea surface temperature, 1856–1991
.
J. Geophys. Res.
,
103
,
18 567
18 589
, doi:.
Kay
,
J. E.
, and Coauthors
,
2014
:
The Community Earth System Model (CESM) Large Ensemble project: A community resource for studying climate change in the presence of internal climate variability
.
Bull. Amer. Meteor. Soc.
,
96
,
1333
1349
, doi:.
Kelley
,
C. P.
,
S.
Mohtadi
,
M. A.
Cane
,
R.
Seager
, and
Y.
Kushnir
,
2015
:
Climate change in the Fertile Crescent and implications of the recent Syrian drought
.
Proc. Natl. Acad. Sci. USA
,
112
,
3241
3246
, doi:.
Lean
,
J.
,
2000
:
Evolution of the sun’s spectral irradiance since the Maunder Minimum
.
Geophys. Res. Lett.
,
27
,
2425
2428
, doi:.
Mann
,
M. E.
, and Coauthors
,
2009
:
Global signatures and dynamical origins of the Little Ice Age and Medieval Climate Anomaly
.
Science
,
326
,
1256
1260
, doi:.
Meehl
,
G. A.
, and
A. X.
Hu
,
2006
:
Megadroughts in the Indian monsoon region and southwest North America and a mechanism for associated multidecadal Pacific sea surface temperature anomalies
.
J. Climate
,
19
,
1605
1623
, doi:.
Newman
,
M.
,
2007
:
Interannual to decadal predictability of tropical and North Pacific sea surface temperatures
.
J. Climate
,
20
,
2333
2356
, doi:.
Newman
,
M.
,
2013
:
An empirical benchmark for decadal forecasts of global surface temperature anomalies
.
J. Climate
,
26
,
5260
5269
, doi:.
Newman
,
M.
,
P. D.
Sardeshmukh
, and
C.
Penland
,
2009
:
How important is air–sea coupling in ENSO and MJO evolution?
J. Climate
,
22
,
2958
2977
, doi:.
Newman
,
M.
,
M. A.
Alexander
, and
J. D.
Scott
,
2011
:
An empirical model of tropical ocean dynamics
.
Climate Dyn.
,
37
,
1823
1841
, doi:.
Newman
,
M.
, and Coauthors
,
2016
:
The Pacific decadal oscillation, revisited
.
J. Climate
,
29
,
4399
4427
, doi:.
Oglesby
,
R.
,
S.
Feng
,
Q.
Hu
, and
C.
Rowe
,
2012
:
The role of the Atlantic multidecadal oscillation on medieval drought in North America: Synthesizing results from proxy data and climate models
.
Global Planet. Change
,
84–85
,
56
65
, doi:.
Otto-Bliesner
,
B. L.
, and Coauthors
,
2016
:
Climate variability and change since 850 CE: An ensemble approach with the Community Earth System Model
.
Bull. Amer. Meteor. Soc.
,
97
,
735
754
, doi:.
Penland
,
C.
, and
L.
Matrosova
,
1994
:
A balance condition for stochastic numerical models with application to the El Niño–Southern Oscillation
.
J. Climate
,
7
,
1352
1372
, doi:.
Penland
,
C.
, and
P. D.
Sardeshmukh
,
1995
:
The optimal growth of tropical sea surface temperature anomalies
.
J. Climate
,
8
,
1999
2024
, doi:.
Rayner
,
N. A.
,
D. E.
Parker
,
E. B.
Horton
,
C. K.
Folland
,
L. V.
Alexander
,
D. P.
Rowell
,
E. C.
Kent
, and
A.
Kaplan
,
2003
:
Global analyses of sea surface temperature, sea ice, and night marine air temperature since the late nineteenth century
.
J. Geophys. Res.
,
108
,
4407
, doi:.
Redmond
,
K. T.
, and
R. W.
Koch
,
1991
:
Surface climate and streamflow variability in the western United States and their relationship to large-scale circulation indices
.
Water Resour. Res.
,
27
,
2381
2399
, doi:.
Ropelewski
,
C. F.
, and
M. S.
Halpert
,
1986
:
North American precipitation and temperature patterns associated with the El Niño/Southern Oscillation (ENSO)
.
Mon. Wea. Rev.
,
114
,
2352
2362
, doi:.
Ropelewski
,
C. F.
, and
M. S.
Halpert
,
1996
:
Quantifying Southern Oscillation–precipitation relationships
.
J. Climate
,
9
,
1043
1059
, doi:.
Seager
,
R.
,
N.
Graham
,
C.
Herweijer
,
A. L.
Gordon
,
Y.
Kushnir
, and
E.
Cook
,
2007a
:
Blueprints for medieval hydroclimate
.
Quat. Sci. Rev.
,
26
,
2322
2336
, doi:.
Seager
,
R.
, and Coauthors
,
2007b
:
Model projections of an imminent transition to a more arid climate in southwestern North America
.
Science
,
316
,
1181
1184
, doi:.
Sheffield
,
J.
, and
E. F.
Wood
,
2008
:
Global trends and variability in soil moisture and drought characteristics, 1950–2000, from observation-driven simulations of the terrestrial hydrologic cycle
.
J. Climate
,
21
,
432
458
, doi:.
Sheffield
,
J.
,
G.
Goteti
, and
E. F.
Wood
,
2006
:
Development of a 50-year high-resolution global dataset of meteorological forcings for land surface modeling
.
J. Climate
,
19
,
3088
3111
, doi:.
Sheffield
,
J.
,
E. F.
Wood
, and
M. L.
Roderick
,
2012
:
Little change in global drought over the past 60 years
.
Nature
,
491
,
435
438
, doi:.
Smirnov
,
D.
,
M.
Newman
, and
M. A.
Alexander
,
2014
:
Investigating the role of ocean–atmosphere coupling in the North Pacific Ocean
.
J. Climate
,
27
,
592
606
, doi:.
Smith
,
T. M.
,
R. W.
Reynolds
,
T. C.
Peterson
, and
J.
Lawrimore
,
2008
:
Improvements to NOAA’s historical merged land–ocean surface temperature analysis (1880–2006)
.
J. Climate
,
21
,
2283
2296
, doi:.
Solomon
,
A.
, and
M.
Newman
,
2011
:
Decadal predictability of tropical Indo-Pacific Ocean temperature trends due to anthropogenic forcing in a coupled climate model
.
Geophys. Res. Lett.
,
38
,
L02703
, doi:.
Stahle
,
D. W.
, and
J. S.
Dean
,
2011
: North American tree rings, climatic extremes, and social disasters. Dendroclimatology, M. Hughes, T. Swetnam, and H. Diaz, Eds., Developments in Paleoenvironmental Research Series, Vol. 11, Springer, 297–327.
Stevenson
,
S.
,
A.
Timmermann
,
Y.
Chikamoto
,
S.
Langford
, and
P.
DiNezio
,
2015
:
Stochastically generated North American megadroughts
.
J. Climate
,
28
,
1865
1880
, doi:.
Trenberth
,
K. E.
,
J. T.
Fasullo
, and
T. G.
Shepherd
,
2015
:
Attribution of climate extreme events
.
Nat. Climate Change
,
5
,
725
730
, doi:.
Vieira
,
L. E. A.
, and
S. K.
Solanki
,
2010
:
Evolution of the solar magnetic flux on time scales of years to millennia
.
Astron. Astrophys.
,
509
,
A100
, doi:.
Wang
,
J.
,
B.
Yang
,
F. C.
Ljungqvist
,
J.
Luterbacher
,
T. J.
Osborn
,
K. R.
Briffa
, and
E.
Zorita
,
2017
:
Internal and external forcing of multidecadal Atlantic climate variability over the past 1,200 years
.
Nat. Geosci
,
10
,
512
517
, doi:.
Woodhouse
,
C. A.
, and
J. T.
Overpeck
,
1998
:
2000 years of drought variability in the central United States
.
Bull. Amer. Meteor. Soc.
,
79
,
2693
2714
, doi:.
Woodhouse
,
C. A.
,
D. M.
Meko
,
G. M.
MacDonald
,
D. W.
Stahle
, and
E. R.
Cooke
,
2010
:
A 1,200-year perspective of 21st century drought in southwestern North America
.
Proc. Natl. Acad. Sci. USA
,
107
,
21 283
21 288
, doi:.
Yuval
, and
W. W.
Hsieh
,
2002
:
The impact of time-averaging on the detectability of nonlinear empirical relations
.
Quart. J. Roy. Meteor. Soc.
,
128
,
1609
1622
, doi:.
Zebiak
,
S. E.
, and
M. A.
Cane
,
1987
:
A model El Niño–Southern Oscillation
.
Mon. Wea. Rev.
,
115
,
2262
2278
, doi:.

Footnotes

© 2018 American Meteorological Society. For information regarding reuse of this content and general copyright information, consult the AMS Copyright Policy (www.ametsoc.org/PUBSReuseLicenses).