## Abstract

Multidecadal “megadroughts” were a notable feature of the climate of the American Southwest over the Common era, yet we still lack a comprehensive theory for what caused these megadroughts and why they curiously only occurred before about 1600 CE. Here, we use the Paleo Hydrodynamics Data Assimilation product, in conjunction with radiative forcing estimates, to demonstrate that megadroughts in the American Southwest were driven by unusually frequent and cold central tropical Pacific sea surface temperature (SST) excursions in conjunction with anomalously warm Atlantic SSTs and a locally positive radiative forcing. This assessment of past megadroughts provides the first comprehensive theory for the causes of megadroughts and their clustering particularly during the Medieval era. This work also provides the first paleoclimatic support for the prediction that the risk of American Southwest megadroughts will markedly increase with global warming.

## INTRODUCTION

Megadroughts (*1*, *2*) are multidecadal periods of intense aridity that existed over much of the American West from the 9th to 16th centuries [see (*3*) for a review]. The southwestern region of the United States (hereinafter referred to as the Southwest) represents one of the main foci of megadrought activity, the occurrence of which profoundly affected organized societies of the time as well as the landscape and vegetation across the area (*3*). The absence of any megadroughts after the 16th century in the Southwest raises critical questions surrounding their causes, their predominant clustering during the Medieval period, and whether there is a substantial risk of these droughts returning in the near future as a consequence of either natural or anthropogenic causes (*4*, *5*).

The El Niño–Southern Oscillation (ENSO) system and its modulation of storm track variability have been the primary driver of seasonal droughts in much of North America over the Instrumental era (*6*). Therefore, the most prominent Southwest megadrought hypotheses invoke ENSO as the probable cause but in different fashions: Either ENSO was radiatively forced into a persistent La Niña–like mean state [through the dynamical thermostat mechanism (*7*–*9*)] or internal climate variability simply produced multidecadal periods of increased La Niña events (*10*–*15*). Recent work nevertheless provides little evidence for a persistent La Niña–like mean state during the Medieval period (*16*–*21*), and although stochastic ENSO variability can explain isolated megadroughts, their clustering before 1600 CE is extremely difficult to explain with the ENSO variability that is characteristic of the Instrumental era (*13*, *19*, *22*). Additionally, a warm North Atlantic may suppress summer precipitation in the Great Plains and the Southwest [e.g., (*23*–*26*)], but the support for this mechanism causing or increasing the likelihood of megadroughts is relatively tentative (*19*).

## RESULTS

Here, we use a newly available reconstruction product in tandem with a global forcing estimate to directly test the association between Southwest megadroughts and the larger oceanic and radiative forcing conditions. The Paleo Hydrodynamics Data Assimilation product (PHYDA) is a physically and internally consistent reconstruction product that provides gridded hydroclimate and sea surface temperatures (SSTs) globally over the Common era in addition to hydroclimatically relevant dynamical indices. PHYDA is a probabilistic estimate of the time-varying state of the climate over the past 2000 years based on a fusion of both a global climate model and a global collection of proxy time series [see Materials and Methods and (*27*) for more information about the construction of PHYDA]. Using the probabilistic, ensemble reconstructions of the Palmer drought severity index (PDSI) from PHYDA, we identify 14 decadal-scale droughts, called megadroughts herein, in the Southwest during the period 800–1925 (see Materials and Methods); all highly probable megadroughts occur before the year 1600 [>99% ensemble agreement (fig. S1A) and in agreement with (*3*)]. Therefore, we generally split our analyses into pre- and post-1600 periods.

We find that megadroughts are associated with unusually cold values of the NINO3.4 index, with megadrought periods shifted to colder temperatures (Fig. 1D) especially for the lower percentiles of the NINO3.4 probability density functions (Fig. 1G), although there is no discernible shift in the mean state of NINO3.4. The North Atlantic was consistently warmer during megadroughts (Fig. 1, E and H), indicative of a large-scale state shift during the period 800–1600 relative to a reference period of 1601–1925 (Fig. 1, E and H, and fig. S2); this state shift is consistent with the well-known global cooling trend over the preindustrial Common era (*28*). These SST results in Fig. 1 are insensitive to the choice of drought definition [(*29*) versus (*10*)] and the choice of annual or seasonal climate indices [annual versus JJA (June, July, and August) Atlantic multidecadal oscillation (AMO) index and monthly versus DJF (December, January, and February) versus annual NINO3.4 index]. Such a lack of seasonal dependence is likely because the AMO index from 800 to 1600 is elevated across all seasons and the extremes in the NINO3.4 index are dominated by the DJF months, the peak season of ENSO. We also find that the forced climate response (derived from an energy balance model and historical forcing estimates; see Materials and Methods) consists of fluctuations of warmer global mean temperatures during megadroughts relative to the 1601–1925 reference period (Fig. 1, F and I); these warmer temperatures are a result of periods of quiescent volcanoes and elevated solar forcing, filtered by the long-term memory of the deep ocean (fig. S3B).

To assess the relative importance of these climate drivers in explaining hydroclimate in the Southwest, we performed a decadal-scale linear regression analysis over the period 800–1925 using detrended indices (see Materials and Methods). On a decadal time scale, the NINO3.4 index explains the most variance of Southwest PDSI (up to 30%), while the AMO and radiative forcing each explain, on average, less than half the variance explained by the NINO3.4 index (Fig. 2A). When multiple predictors are included, there is a small improvement in explained variance, but the dominant factor is the decadal NINO3.4 index (Fig. 2A). Note that our analysis here focuses on large-scale climate phenomena that act on the annual to multiyear time scales. Other shorter-scale climate phenomena, such as stochastic synoptic variability, likely contribute to drought and megadrought in the Southwest (*6*), but from the perspective of the paleoclimate record, the time history of these phenomena is very likely unrecoverable. Such phenomena likely account for the unexplained variance in these regressions.

We also tested the influence of radiative forcing on the NINO3.4 index in the decadal regression framework as a measure of the potential influence of a dynamical thermostat mechanism (*7*); within this framework, we find no influence of radiative forcing on the NINO3.4 index (Fig. 2B). In contrast, however, we do find that a small, similar amount of variance in both the AMO and local Southwest temperature can be explained by radiative forcing (Fig. 2B). For the Southwest, we interpret the locally forced response as a straightforward mechanism: A global-scale positive radiative forcing leads to an increase in local temperatures, which, in turn, lead to an increase in aridity. In the temperature units of the global forcing estimate, the median local temperature response is 0.12°C; on a decadal time scale, this would represent a decrease in Southwest PDSI of about 0.17 (using the linear relationship between decadal, detrended temperature, and PDSI over the Southwest from 800 to 1925; Fig. 2C). This forced PDSI response is only a fraction of the PDSI excursions seen in PHYDA megadroughts of approximately −1 PDSI (fig. S1D); yet, a forced PDSI value of 0.17 is about half the value required for a proposed mean state shift in Southwest PDSI seen in the North American Drought Atlas during the Medieval megadrought period (*22*).

After exploring the drivers of decadal hydroclimate in the American Southwest generally, we now turn to the causes of megadroughts specifically. For this, we use a logistic regression framework, given that we are treating megadroughts as binary events. We find that megadroughts are most robustly associated with cold NINO3.4 conditions and warm AMO conditions, with the influence of NINO3.4 being again more than twice as important (Fig. 2D; here, the predictors have been standardized, thus allowing for a direct comparison of regression parameter values). The role of a positive radiative forcing is less certain: Approximately three-fifths of the bootstrap regression parameter estimates are positive when four outlier large volcanic eruption data points are removed (β_{1} distribution for “Forc-NLV” in Fig. 2D; see Materials and Methods). Large volcanic eruptions (negative forcing) appear to have no consistent temporal relationship to megadroughts, with, for example, volcanic eruptions bookending a megadrought in the mid-1200s or occurring in the middle of a megadrought in the late 1400s (Fig. 1C). However, 9 of 14 megadroughts are associated with a positive shift in radiative forcing (Fig. 1, F and I), caused in the energy balance model by elevated solar forcing and volcanic quiescence (fig. S3B). Additionally, a positive radiative forcing on drought is supported by the linear regression analysis in Fig. 2C. To test whether the positive forcing response in the regression model could be consistent with a null variable, we reran the logistic regressions including a Gaussian random predictor; we found that the β_{1} distribution for Forc-NLV in Fig. 2D was statistically distinct from that of the random predictor (Kolmogorov-Smirnov test with *n* = 999 and α = 0.01). Given all this evidence and the simplicity of a positive forcing mechanism, we argue that a positive radiative forcing influence on megadroughts is more likely than a negative radiative forcing or no forcing influence at all.

Considering the joint influence of NINO3.4, AMO, and radiative forcing in causing megadroughts, we find that the logistic regression model that has the highest probability of being the most explanatory is the one that includes all three of these climate drivers (Fig. 2E); this is indicated by two different information criteria, both of which include a penalty for adding predictors. These logistic regression results are also consistent with the SST and forcing excursions shown in Fig. 1, D to I. We therefore argue that Southwest megadroughts were primarily driven by anomalously frequent and cold, unforced La Niña conditions, with contributions from a partially forced warm Atlantic and a forced local temperature increase.

The analysis shown in Figs. 1 and 2 assumes from the outset that the specific NINO3.4 and AMO indices are suitable for uncovering the mechanisms of megadroughts. While the choices of these indices are not without warrant based on past work, we nevertheless seek to confirm our results without assuming the form of the variables driving megadroughts. Climate pattern frequency analysis using self-organizing maps (SOMs) (*30*) (see Materials and Methods) reveals spatial SST patterns that are consistent with the Pacific and Atlantic forcing conditions shown previously. The SOM analysis specifically demonstrates a robust increased incidence of cold central and eastern Pacific SSTs during megadroughts relative to the 1601–1925 reference period, with the responses for the stronger nodes 1 and 2 being the most distinct compared to the reference range [Fig. 3, A and B (nodes 1, 2, and 4), and figs. S4, A and B (node 2), and S5, A and B (nodes 1 to 3)]; these specific patterns are concomitant with dry conditions across the western United States, with the strongest SST patterns generally having the most consistently dry PDSI values. Note that these SOM patterns are emergent properties of the underlying dataset, are not imposed a priori, and have no explicit dependence or knowledge of Southwest megadrought conditions. The greatest pattern increases are associated with cold Pacific SSTs co-occurring with areas of warm Atlantic SSTs, although there is some sensitivity to the precise pattern of warm Atlantic SSTs (cf. Fig. 3A, node 2, and fig. S5A, nodes 2 and 3). We note, however, that because of required data preprocessing and detrending (see Materials and Methods), the warming shift in the Atlantic associated with the Medieval Climate Anomaly is partially removed in this SOM analysis. During megadroughts, there was also a robust reduction of some warm Pacific SST patterns (Fig. 3, A and B, nodes 5 and 8), which are concomitant with wetting across the western United States. These results, of an increased frequency in cold Pacific SSTs coupled with warming in the Atlantic and dry conditions in the Southwest, are robust to the choice of the number of SOM nodes (figs. S4 and S5) and drought definition [(*29*) versus (*10*)]. These results are also consistent with a basic compositing analysis across the megadroughts (fig. S6). The SOM analysis also confirms that the dominant climate driver of megadroughts is unusual periods of La Niña conditions.

## DISCUSSION

The new PHYDA dataset, along with modeling of the historical forced climate, has allowed us to construct a more comprehensive theory of Southwest megadroughts than has heretofore been possible. This theory unifies SST and radiative forcing mechanisms that act together to cause megadroughts. We confirm the hypothesis that anomalously frequent periods of cold Pacific SSTs can lead to megadroughts (*10*–*13*), but we also find that these SST patterns in the past were unusually cold and unforced. Over the past 1200 years, this phenomenon is at least twice as important for explaining Southwest megadroughts as an Atlantic forcing or an exogenous radiative forcing (Fig. 2, A and D). Yet, we explain the clustering of megadroughts before circa 1600 CE as a result of changes in each of these climate phenomena: Cold Pacific SST fluctuations are more frequent and more cold (Figs. 1, D and G, and 3), the north Atlantic is warmer (Fig. 1, E and H), and periods of positive radiative forcing are more common (Fig. 1, F and I). Each of these factors contributes to droughts (Fig. 2A), and their changes lead to an increase in the occurrence of megadroughts.

These results imply that ENSO variability, in addition to local positive radiative forcing, has the capacity to induce megadrought conditions at any time in the future. There is evidence that these two factors acting together may explain the drying trend in the Southwest since the early 1980s [e.g., (*31*)]. However, the predictability of a contemporary megadrought may be limited, given that state-of-the-art global climate models have difficulty reproducing important features of ENSO [e.g., (*32*)]. Of perhaps greater concern for the long-term future is the possibility that radiative forcing could gradually come to dominate the hydroclimate of the Southwest, with the recurrence of megadroughts becoming almost assured (*4*, *5*). The results presented here provide paleoclimatic verification of the capacity of radiative forcing to increase the odds of megadrought. Yet, if the past is a useful guide to the next few decades, then megadroughts could also arise through unforced and unpredictable ENSO variability.

## MATERIALS AND METHODS

### Experimental design

*Paleo Hydrodynamics Data Assimilation product*. The PHYDA (*27*, *33*) is based on data assimilation, a method that optimally fuses proxy information with the dynamical constraints of climate models; this method simultaneously estimates both hydroclimate fields and corresponding atmosphere-ocean states (*34*). PHYDA is based on information from a network of 2978 annually resolved proxy-data time series (*35*) together with the Community Earth System Model (CESM) Last Millennium Ensemble of climate model simulations (*36*); as a data assimilation product, PHYDA therefore represents an optimal amalgamation of both model and proxy information [see (*27*) for more details about PHYDA]. The PHYDA reconstructions were performed from the years 1 to 2000 CE targeting three different temporal windows: April to the next calendar year March annual means (indicated as A2M in figure legends), the boreal growing season of JJA, and the austral growing season of DJF. In addition to seasonal and annual variables, PHYDA also includes monthly ENSO indices that are used herein. For computing the probability density functions shown in Fig. 1, we reran the PHYDA reconstruction code (see “Data and materials availability”) to save the full 998-member ensemble for the specific NINO3.4 and AMO index variables.

PHYDA has been extensively validated against observational data in addition to related paleoclimate reconstructions (*27*). Here, we also performed comparisons with the North American Drought Atlas for the specific region of interest (fig. S1). In addition, there may be concerns that the dynamics and temporal variability of PHYDA merely reflect those of the data assimilation prior, namely, the CESM climate model; to explore this possibility, we analyzed the power spectra of the North American Southwest PDSI (NASW PDSI) along with the NINO3.4 and AMO indices and found PHYDA to be much more similar to observations than to the CESM simulation used as the basis of PHYDA (fig. S7). We also found PHYDA to be more similar to observations and less similar to CESM in the cross-correlation of the NASW PDSI with the NINO3.4 and AMO indices (fig. S8). Because the analysis presented here relies on realistic ENSO-hydroclimate teleconnections in PHYDA, we also showed the point correlations between the A2M annual NINO3.4 index and the variables of 2-m surface temperature and PDSI in fig. S9.

*Analysis period*. As in previous work [e.g., (*3*)], we limited the beginning of the analysis to the year 800 due to the relative dearth of proxy information before that time (*27*). For the end date of the analysis, we used the year 1925, thus excluding the period of most pronounced anthropogenic warming (fig. S2).

*NASW region*. For the NASW, we used the average PDSI over land in the regional box bounded by 125°W to 105°W and 31°N to 42°N (fig. S1), which wholly includes the U.S. states of California, Nevada, Utah, Arizona, a southern portion of Wyoming, and both Colorado and New Mexico west of the Rocky Mountains. This box also includes some of the northernmost portions of Mexico where the border deviates from a straight line at 31°N. For nearly all of this region, the majority of the precipitation falls during the winter and spring seasons.

*Megadrought specification*. The specific megadrought periods were determined probabilistically using a 100-member ensemble of NASW PDSI time series based on PHYDA; the PHYDA reconstruction code (see “Data and materials availability”) was reran to output a 100-member ensemble of global PDSI fields (for this analysis, the 100-member ensemble sufficiently characterizes the full ensemble, which is much larger and more cumbersome to analyze), from which the averaged NASW PDSI for each ensemble member was derived. Drought periods in each of these time series were identified as the persistently negative PDSI values relative to an 11-year moving average (*29*). We then found the ensemble agreement about the existence of a drought for each year of the analysis period 800–1925 (fig. S1A). The 14 specific megadroughts were then identified as those periods that had at least 90% ensemble agreement and that were at least 10 years long (fig. S1A). The general results presented here are robust to the particular choice of threshold (for agreement values greater than 50%) or if the megadroughts are determined nonprobabilistically using the ensemble mean NASW PDSI from PHYDA. The results are also robust to the choice of drought definition [(*29*) versus (*10*)].

*Forced climate response*. We estimated the historical forced climate response with a recently developed energy balance model (*37*, *38*) that has been shown to capture both the fast and slow components of atmosphere-ocean global climate models as they respond to time varying forcings (*37*). This coupled linear model is given by*T*_{s} and *T*_{d} are the temperature anomalies of the surface components of the climate system and the deep ocean below the mixed layer; *F* is the anomalous radiative forcing (W m^{−2}); λ is the global climate feedback factor (W m^{−2} K^{−1}); γ is an exchange coefficient determining the strength of the coupling and the rate of temperature exchange between the surface and deep ocean (W m^{−2} K^{−1}); *c*_{s} is the effective heat capacity of the surface components of the climate system, which are dominated by the ocean mixed layer (W yr m^{−2} K^{−1}); *c*_{d} is the effective heat capacity of the deep ocean, approximately set by the depth of heat uptake or release during transient forcing events (*39*) (W yr m^{−2} K^{−1}); and ε represents the ocean heat uptake efficacy (unitless).

We used this model to generate an ensemble of probable forced temperature histories over the Common era. The process of finding realistic model parameter combinations first began by randomly drawing 10,000 sets of parameters (γ, *c*_{s}, *c*_{d}, and ε), sampling with uniform likelihood from the range of parameter fits found in (*38*) (γ = [0.49 to 1.06], *c*_{s} = [6.2 to 9.7], *c*_{d} = [56 to 271], and ε = [0.83 to 1.82]; these fits were computed using 16 atmosphere-ocean global climate models that simulated an abrupt 4 × CO_{2} experiment within the CMIP5 framework). In addition, we used 10,000 randomly drawn values of λ that correspond to an equilibrium climate sensitivity (ECS) range of 1° to 10°C, consistent with the potential range from (*40*), and used the relation that λ_{eq} is equal to the forcing from CO_{2} doubling (*F*_{2×}) divided by the ECS, λ_{eq} = *F*_{2×}/ECS = 3.44 W m^{−2}/ECS [e.g., (*41*)], so that the range of λ = [0.344 to 3.44]. We then ran numerical simulations of the energy balance model using the total instantaneous radiative forcings from (*42*), *F* in Eq. 1, for the years 1850–2005. We kept the parameter combinations from the best 10% (or *n* = 1000) of these simulations, as measured by the mean absolute error relative to the historical global mean temperature over the period 1850–2005 (fig. S3A) (*43*). These parameter combinations were then used in simulating the years 1–2000 using historical estimates of greenhouse gases (CO_{2}, N_{2}O, and CH_{4}) (*44*) converted to radiative forcing via equations from (*45*), solar irradiance (*46*), and global volcanic forcing (*47*). The volcanic forcing from (*47*) is given as year-specific forcing events, which we converted to a continuous time series by assuming a volcanic forcing *e*-folding time scale of 13 months, consistent with an observed range of about 12 to 14 months (*48*, *49*); we also incorporated volcanic forcing uncertainty information by generating an ensemble of 1000 forcing time series with Gaussian random perturbations of the forcing estimate based on the magnitude of event-specific uncertainties from (*47*). Using both the parameter combinations and total forcing estimates, we generated 1000 forced climate response time series for the years 1–2000; fig. S3B shows the total forced climate response in this model and the response to individual forcing components.

While this model is not spatially resolved, we used it here for two reasons: (i) There are no local, Southwest forcing estimates available before the Instrumental era. All Common era forcing estimates are available on the hemispheric or global scale (with any reconstructed spatial pattern of volcanic forcing based primarily on ice core records in only Greenland and Antarctica). A global-scale model is therefore a conservative use of the available forcing information. (ii) The use of a low-order model, as opposed to a global climate model, allows us to use an updated volcanic forcing estimate (*47*) and to include uncertainties associated with the volcanic forcing. Including uncertainties (here within a Monte Carlo framework) is particularly important, given that volcanoes are the dominant forcing over the Common era and that their magnitude and spatial extent can be quite uncertain.

Despite the global nature of the energy balance model, it can still be representative of regional information. For example, the decadal mean correlation between the energy balance model global mean temperature and the mean temperature of the Southwest from PHYDA is 0.4 over the analysis interval of 800–2000.

### Statistical analysis

*Linear regression analysis*. We performed a decadal-scale linear regression analysis to assess the relative importance of climate drivers in explaining decadal hydroclimate over the period 800–1925 (Fig. 2A); for this panel, decadal-scale PDSI was the predictand and decadal-scale NINO3.4, AMO, and total forcing were the predictors. We performed the univariate and multivariate regressions within a bootstrap sampling framework that sampled both the choice of averaging and the values in the time series: The ensemble mean time series shown in Fig. 1 were averaged at 9-, 10-, and 11-year block averages and randomly sampled with replacement such that 999 (3 × 333) total regressions were performed for each predictor/predictand(s) option. The bootstrap sampling used here (and for all the panels in Fig. 2) allowed us to incorporate the uncertainty associated with the choice of averaging time and to more accurately estimate regression parameters with a relatively small sample size (*n* ≈ 110 at decadal time scales). Before computing the regressions, the averaged time series were detrended to remove autocorrelation (the removal of which was confirmed). The adjusted *r*^{2} values (adjusted for the number of predictors) for all 999 bootstrap regressions are summarized with box plots in Fig. 2A.

For Fig. 2 (B and C), we performed a similar decadal-scale, bootstrap regression analysis but with NINO3.4, AMO, and NASW averaged temperature as predictands and the total forcing as the predictor. We also showed the distribution of regression parameter values β_{1} (Fig. 2C), as in the regression equation *y* = β_{0} + β_{1}*x* + ε. These bootstrap-generated distributions provide an estimate of the confidence in β_{1}, although they are not identical to formal linear regression confidence intervals, which rely on assumptions that we did not make with the bootstrap approach. Note that the predictors here have not been standardized because we used the magnitude of the regression parameters in the main text; thus, the relative importance of the predictors cannot be deduced from Fig. 2C. Note also that because this analysis only captures linear relationships, there may be further nonlinear causal factors that are not accounted for.

*Logistic regression analysis*. We performed a logistic regression analysis to directly assess the role of climate drivers in explaining megadroughts in the Southwest (Fig. 2, D and E). For this analysis, the predictand is the binary variable of being in a megadrought or not being in a megadrought and the predictors are the ensemble mean NINO3.4, AMO, and total forcing time series as in Fig. 1C (local temperature is not included in the regression model because PDSI is explicitly dependent on temperature). Because of the strong autocorrelation in the annual megadrought time series, we binned all the relevant annual time series according to the boundaries of the start and end years of the megadroughts. Each megadrought period is thus a megadrought bin (*n* = 14), and the surrounding periods are nonmegadrought bins that are divided according to the following rule: If a given nondrought period is larger than the average megadrought bin size of 22 years, then that period is divided into the largest possible equal (or equal ±1)–sized bins such that no bin is larger than 22 years; the continuous nondrought periods at the beginning and end of the time series also follow this rule. The predictor variables of NINO3.4, AMO, and total forcing are averaged over each bin segment. The logistic regression for the probability of megadrought during a given bin, π = Pr (*Y* = 1 ∣ *X* = *x*), thus takes the form logit(π) = log [π/(1 − π)] = β_{0} + β_{1}*x*_{1} + … + β* _{k}x_{k}* for

*k*predictors. Because the bins are of varying sizes, the regression is weighted by the size of the bin. In addition, each variable is standardized so that the relative influence of each variable can be gauged by comparing the β

*values. To address the issue of small sample size (*

_{i}*n*= 54 total bins) and to be commensurate with Fig. 2 (A to C), we performed the regression using a bootstrap sampling with replacement 999 times for each predictor/predictand(s) option. From this, we showed the β

_{1}distributions (the β

*distributions for multivariate logistic regression were very similar), and on the basis of the bootstrap sampling, we showed two different information criteria for determining the best regression model: the Akaike information criterion (AIC) corrected for the sample size and the Bayesian information criterion (BIC); both information criteria include a penalty term for the number of model parameters. For each bootstrap regression, the same randomly drawn data points were used for all the predictor/predictand(s) options, and the regression model with the minimum AIC and BIC was recorded; the bars in Fig. 2E show the percentage of the bootstrap iterations where a given model had the lowest AIC and BIC and was thus the “best” model for those randomly drawn data points.*

_{i}A further complication is that a few very large volcanic eruptions in the forcing time series introduce outliers that can influence the regression. Thus, in Fig. 2 (D and E), we showed two logistic regression analyses with one using the full time series and another where four outlying bins were removed for all regression variables (these corresponded to two megadrought bins and two nondrought bins); the outliers in the total forcing time series were considered as those values that were more than three scaled median absolute deviations away from the median (the default outlier detection method of MATLAB).

*SOM analysis*. SOMs are a neural network-based cluster analysis (*30*) that has been used recently in the climate sciences to describe the continuum of climate patterns within a dataset [e.g., (*50*, *51*)]. An analysis of climate modes using SOMs is conceptually similar to traditional principal component analysis: The “nodes” of a SOM represent the primary modes of the underlying climate field. A critical advantage of a SOM analysis compared to a principal component analysis is that it does not impose linearity and orthogonality on states that may be related in nonlinear, nonorthogonal ways; practically, this feature makes SOMs more likely to correspond to physically plausible patterns [e.g., (*52*)].

Here, we used the SOM algorithm to assign annual SST anomaly fields to spatial patterns of a preset number. The SOM algorithm creates spatial patterns that maximize their similarity to the underlying SST fields (by minimizing their Euclidian distance) and then assigns each annual SST field to the best matching pattern. The SOM patterns are approximately the mean of the assigned SST fields and are thus approximately a composite of relatively similar SST fields (*50*). In addition, the SOM analysis organizes the patterns such that similar patterns are assigned to nearby locations within a regular two-dimensional grid. Thus, this full process allows one to visualize a reduced space continuum of patterns in the dataset (Fig. 3 and figs. S4 and S5).

For the SOM analysis, the SST fields are preprocessed by detrending the SST fields (removing the linear trend at each grid point over the full analysis period, 800–1925; c.f. fig. S2); this detrending prevents spurious trends in the pattern occurrence (*51*), although it does remove somewhat the influence of the warm AMO in the early centuries of the analysis. As in (*50*), we applied an area weighting of the SST fields according to the cosine of latitude because PHYDA is on a regular latitude-longitude grid. The SOM algorithm keeps track of which year’s SST field has been assigned to which best matching pattern (or “best matching unit”), and we made composites of JJA PDSI over the years assigned to each pattern, as shown in Fig. 3 and figs. S4 and S5.

For the bootstrap resampling for Fig. 3B, we performed a 500 iteration bootstrap resampling of the same number of years as in the megadroughts, drawn from the reference period; corresponding ranges for the nondrought years are very similar. Thus, they are not shown for clarity.

## SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/5/7/eaax0087/DC1

Fig. S1. Ensemble drought agreement and comparison of NASW PDSI indices.

Fig. S2. Global mean temperature from PHYDA.

Fig. S3. Energy balance model estimates of the historical, global mean forced climate response.

Fig. S4. SOM analysis using six SOM nodes.

Fig. S5. SOM analysis using 12 SOM nodes.

Fig. S6. Temperature and PDSI composite over the megadroughts.

Fig. S7. Power spectra comparisons.

Fig. S8. Cross-correlation comparisons.

Fig. S9. NINO3.4 correlation maps.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is **not** for commercial advantage and provided the original work is properly cited.

## REFERENCES AND NOTES

**Acknowledgments:**

**Funding:**This work was supported in part by the NOAA Climate and Global Change Postdoctoral Fellowship Program administered by UCAR’s Visiting Scientist Programs. This work was also supported in part by the NSF under grants AGS-1243204, AGS-1401400, AGS-1602581, AGS-1602920, OISE-1743738, AGS-1805490, and AGS-1703029, and in part by the NASA Modeling, Analysis, and Prediction Program, NASA-80NSSC17K0265. LDEO contribution number 8325.

**Author contributions:**N.J.S. developed the code, designed experiments with J.E.S., ran experiments, performed statistical analyses of experiments, and produced figures. N.J.S. and J.E.S. prepared the manuscript with assistance from B.I.C., R.S., A.P.W., and E.R.C. In addition, B.I.C. and R.S. reviewed and provided the feedback on experimental results, while E.R.C. aided in the interpretation of the North American Drought Atlas.

**Competing interests:**The authors declare that they have no competing interests.

**Data and materials availability:**All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. PHYDA is publicly available at the Zenodo data repository as NetCDF4 files: https://doi.org/10.5281/zenodo.1154913. The code for computing the SOMs is publicly available from https://github.com/ilarinieminen/SOM-Toolbox. The code for rerunning PHYDA to output ensemble estimates is available at https://github.com/njsteiger/PHYDA, and the 100-member ensembles of PHYDA used here are available at http://clifford.ldeo.columbia.edu/nsteiger/recon_output/phyda_ens/. Additional data related to this paper may be requested from the authors.

- Copyright © 2019 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).