Research ArticleCLIMATOLOGY

Atmospheric dynamics drive most interannual U.S. droughts over the last millennium

See allHide authors and affiliations

Science Advances  07 Aug 2020:
Vol. 6, no. 32, eaay7268
DOI: 10.1126/sciadv.aay7268


The American West exemplifies drought-sensitive regions with growing populations. Paleoclimate investigations have documented severe droughts in this region before European settling, with major implications for water management and planning. Here, we leverage paleoclimate data assimilation to reconstruct past climate states, enabling a large-scale multivariate investigation of U.S. drought dynamics over the last millennium. These results confirm that La Niña conditions significantly influence southwest U.S. drought over the past millennium but only account for, by one metric, ~13% of interannual drought variability in that region. Atlantic sea surface temperatures may also contribute a small influence, but unexplained variability suggests a substantial role for internal atmospheric variability. This conclusion is buttressed by analysis of simulations from the Community Earth System Model Last Millennium Ensemble. While greenhouse gases will increase future drought risk, as shown in other work, interannual U.S. drought variations will also be widely influenced by processes internal to the atmosphere.


Recent droughts across the United States have resulted in loss of life and billions of dollars in damage (13), making drought forecasting and planning a high societal priority. Yet, the paleoclimate record shows that these droughts pale in comparison to the megadroughts of the Common Era (C.E.), many of which appear to be longer lasting than the famed “Dust Bowl” of the 1930s (4, 5) and had similarly marked societal impacts. For instance, the drought of 1276 to 1299 C.E. likely contributed to the migration of the Anasazi people out of the American southwest near the beginning of the 14th century (5, 6). Because continued greenhouse gas emissions are projected to increase drought risk over much of the United States and other regions (79), understanding the dynamical causes of past droughts is of paramount importance.

U.S. droughts have been associated with changes in the El Niño/Southern Oscillation (ENSO) (4, 1012), sea surface temperatures (SSTs) in the Atlantic Ocean (10, 11, 13, 14), internal atmospheric dynamics (1517), and greenhouse forcing (2). Variations in the Atlantic Multidecadal Oscillation (AMO) have also been suggested as an influence on past “megadroughts” (15), especially in conjunction with La Niña–like SSTs (1820). Much of this past work relies either on the instrumental record, which is relatively short (21), or on paleoclimate data processed in different ways (10, 15, 22). Here, we use a data assimilation approach called the Last Millennium Reanalysis (LMR) (23, 24) to combine paleoclimate proxy data with climate model output within a unified framework to explore the climate drivers of U.S. drought.

The LMR methodology (Materials and Methods) combines a network of annually resolved proxy records—including tree rings, corals, and ice cores (Fig. 1)—with output from a general circulation model (GCM), to estimate a gridded multivariate record of climate variability over the past two millennia (23, 24), similar to the Paleo Hydrodynamics Data Assimilation product (PHYDA) reconstruction (25). All temporal variations in the reanalysis are informed by the proxy records weighted against a time-independent climate estimate from the GCM prior. Fundamentally, the network of proxy records sample the climate system at different locations and seasons and for different spans of time, and climate variable relationships from a climate model are used to synthesize these diverse records to reconstruct a dynamically consistent view of past climate. The reanalysis represents a mix of information from the proxies and the GCM, relying more on GCM covariances and remote proxies where local proxies are scarce.

Fig. 1 Proxies used in the LMR reanalysis.

(A) Proxy locations and (B) the number of proxies through time. Proxies come from PAGES2k v2.0.0, Breitenmoser et al. (62), and the National Centers for Environmental Information paleoclimate archives. In total, 2787 proxies are used: 2556 trees, 116 corals and sclerosponges, 105 ice records, 9 lake sediments, and 1 bivalve.


Climate reconstructions via data assimilation

The LMR has been extensively validated against instrumental data and independent proxy data not used in assimilation (24). Compared to instrumental datasets and modern reanalysis products, global-mean temperature correlations range from 0.88 to 0.94 and coefficients of efficiency [(23, 26), here calculated using identical calibration and verification periods] range from 0.77 to 0.87 (fig. S1). Spatial skill is highest over tropical and mid-latitude oceans, similar to previous findings (23). Additional verification in a similar set of experiments, including verification against withheld proxies and tests using separate calibration and verification periods, can be found in recent work (24).

Annual-mean climate indices in the LMR also compare well with observationally based datasets. Three large-scale modes are inspected: the Nino 3.4 index (SST variations in the 5°S to 5°N, 170°W to 120°W region of the Pacific Ocean), the Pacific Decadal Oscillation (PDO; calculated as the first principal component of Pacific SST anomalies north of 20°N in the LMR), and the AMO (calculated as the mean of North Atlantic SSTs detrended over the period 1856 to 2000 in the LMR). A comparison between the Niño 3.4 index derived from the LMR ensemble and from observational data (27) gives a correlation of 0.79 and a coefficient of efficiency of 0.40 over 1873 to 2000 C.E. For PDO, the LMR matches an observationally based time series (28) with a correlation of 0.65. For the AMO, the LMR matches the annualized AMO index from the Earth System Research Laboratory (29) with a correlation of 0.65 (fig. S2). This general agreement provides evidence that the LMR, using data from a network of proxy records, skillfully reconstructs large-scale climate indices given sufficient proxies.

Drought is quantified using the Palmer Drought Severity Index (PDSI), a common hydroclimate metric whose negative values are indicative of moisture deficit (30). First, we examine the full range of PDSI variations to analyze the effects of large-scale climate patterns on U.S. hydroclimate in general. Then, we define a drought threshold (here, all years with PDSI below 1 standard deviation compared to the nearest 51-year period) to evaluate the conditions that are relevant for dry periods alone, for different regions of the United States. For the LMR, PDSI is calculated in the model prior using the Penman-Monteith method for estimating the potential evapotranspiration component (31) and then reconstructed in the past through data assimilation alongside other climate variables (Materials and Methods).

For the purpose of cross-validation, the LMR-derived PDSI is compared to two PDSI datasets based on different information sources. The first dataset, the North American Drought Atlas (NADA), is a summer [June-July-August (JJA)] PDSI reconstruction based on 1845 tree records (32, 33). A comparison between LMR and NADA shows high values of correlation and coefficient of efficiency over most of the United States, especially in densely sampled regions (Fig. 2). This agreement is remarkable despite the differences in seasonality between the two datasets (annual-mean values for LMR versus June to August for NADA) and the very different methodologies. When averaged regionally, the correlation between LMR and NADA from 1001 to 2000 C.E. is 0.84 for the U.S. region as a whole (30°N to 49°N, 130°W to 65°W, land only) and 0.90 for the southwest United States [32°N to 40°N, 125°W to 105°W, land only, as in (10)]. The timing of PDSI variations compares well between the LMR and NADA datasets, especially over the southwest United States (Fig. 2), although the LMR displays reduced variability at basically all time scales. This reduced variability compared to NADA could result from several differences in methodology, including NADA’s use of variance restoration, different approach to tree ring detrending, or proxy prewhitening. Differences between LMR and NADA over Canada and Mexico may be due to the existence of fewer records there; NADA generally shows reduced skill in those regions [see figure 6 in (32)]. A shortcoming of this comparison, however, is an overlap in proxy data sources informing the reconstructions (fig. S3).

Fig. 2 PDSI comparison between LMR and NADA.

Coincident (A) correlation and (B) coefficient of efficiency (CE) calculated between annual-mean LMR and JJA NADA PDSI at every point for the years of common overlap during 1001 to 2000 C.E. Regional-mean time series of LMR (blue) and NADA (black) for (C) an approximate U.S. region (land within 30°N to 49°N, 130°W to 65°W) and (D) the southwest United States (land within 32°N to 40°N, 125°W to 105°W). The box in (B) represents the southwest U.S. region used in (D) and analyses throughout the paper. Blue shading in (C) and (D) represents the 95% uncertainty bands for the LMR PDSI. Uncertainty bands widen back in time, consistent with proxy attrition. A 10-year running mean has been applied to the time series for plotting clarity [(C) and (D)], but all correlation and coefficient of efficiency values are calculated using annual values. Local proxies used in the LMR data assimilation are shown in (A), with count totals representing proxies over the entire globe. Hatched areas in (A) are not significant at the 95% level according to an isospectral test (69).

To compare LMR to a dataset that does not use proxy data, we also validate the LMR-derived PDSI against a PDSI product derived from instrumental observations (Dai PDSI) (fig. S4) (30). While the LMR uses observational temperature and precipitation data [GISS Surface Temperature Analysis (GISTEMP) and Global Precipitation Climatology Centre (GPCC), respectively] to derive relationships between proxy quantities and modern climate quantities, the LMR methodology has no additional knowledge of modern PDSI values, making this an insightful comparison despite only covering 1850 to 2000 C.E. The LMR-derived PDSI shows good agreement with Dai PDSI, with a correlation of 0.65 over the U.S. region and 0.70 for the southwest United States.

The overall agreement between the LMR and these two PDSI datasets, which are formulated using two very different approaches, lends support to the LMR methodology and increases confidence in the multivariate analysis of North American drought. Furthermore, while correlations between LMR and Dai PDSI are slightly worse than between NADA and Dai PDSI, coefficient of efficiency values are better for LMR in the southwest United States, Mexico, and Canada. This comparison shows that, despite methodological differences, the LMR produces PDSI values in line with other datasets. The advantage of LMR is that it provides dynamically consistent reconstructions of other climate fields, promoting a more in-depth multivariable analysis, which is the focus of the following sections.

Drought versus SST

Understanding the causes of drought variability is critical in water-poor regions like nearly all of the United States westward of the 100th meridian (34). Among forcing mechanisms, drought in the western United States is most often related to ENSO, with a La Niña pattern of cooler eastern equatorial Pacific SSTs correlated with drier conditions in the southwest (4, 1012). Heating anomalies associated with equatorial SSTs modify Rossby wave propagation from the tropics into the extratropics; in particular, cool La Niña SSTs reduce tropical convection and upper level divergence, which affects the location of quasi-stationary waves (35, 36). Positive pressure anomalies over the northern Pacific (i.e., a “blocking” high) divert storm tracks northward, reducing precipitation in the southwest United States. Drought development may also depend on transient eddy activity related to the Pacific storm track, as well as land-atmosphere feedbacks such as the soil moisture feedback, with models suggesting that preexisting dry soils help exacerbate subsequent drought (35, 37). A cold phase of the PDO is also associated with drought in the southwest United States and wet conditions in the northwest (10), though the PDO may affect drought primarily in conjunction with ENSO rather than on its own (18).

To examine these links between drought in the contiguous United States and climate features of the surrounding basins, multivariate climate patterns are analyzed over the past millennium. It is important to note that the causes of drought can vary seasonally, with different factors affecting summer versus winter precipitation patterns, but here, we take a broader view and focus on annual-mean anomalies. Because U.S. drought may be influenced by SSTs in the preceding winter, our analysis focuses on connections between annual-mean PDSI and climate patterns in the coincident year as well as in the previous year. In observational datasets, correlation patterns are similar when comparing annual-mean PDSI to Nino 3.4 averaged over the previous December-January-February (DJF) (which is often used as a target for analysis) or over the previous full year, lending support for this approach (note S1 and fig. S5).

To isolate the patterns that constitute the largest amount of covariance between U.S. PDSI and the surrounding climate system, we conduct a maximum covariance analysis (MCA; see Materials and Methods) (38) between annual-mean PDSI over the United States and a joint field consisting of annual-mean SST and 500-hPa heights in the surrounding regions (Fig. 3). MCA isolates orthogonal patterns that explain the maximum amount of covariance between the fields over the analysis period (here, the last millennium), offering a largely impartial assessment of the relationships between chosen fields in the multivariate LMR reconstruction.

Fig. 3 MCA.

MCA between PDSI and SST/500-hPa heights for mode 1. (A and B) Maps showing the spatial patterns of variability and (C) standardized expansion coefficients showing how the magnitudes of the spatial patterns change through time. In (B), positive height anomalies are indicated by solid contours and negative anomalies are indicated with dashed contours, with the thicker line indicating 0. SSTs and 500-hPa heights are standardized before conducting this analysis, so the values shown are not covariances. Calculations are performed on the data in the regions shown. The squared covariance fraction (SCF; measuring the amount of squared covariance for which each mode accounts) as well as the fraction of variance (FOV; measuring the relative amount of variability explained by this mode for the variable under consideration) are listed in the lower right of (A) and (B). The correlation (r) of the expansion coefficients is given in (C).

The first mode of the MCA outlines a clear connection between PDSI and tropical Pacific SSTs, with southwest and western-central U.S. dry conditions over the last millennium corresponding with La Niña and cold PDO SST patterns (Fig. 3). Geopotential heights at 500-hPa increase primarily over the North Pacific, with a band extending across the United States and part of the North Atlantic. Over the equatorial Pacific, 500-hPa heights are slightly reduced. This Pacific response—lower pressure in the tropics and higher pressure over the northern Pacific—fits the canonical view of La Niña producing a blocking high that diverts Pacific storm tracks to the north. In addition, the pattern of height anomalies, with largest increases in the North Pacific and stretching across the contiguous United States, is consistent with (and opposite to) the decreased heights associated with the type of El Niño events that produced the greatest positive precipitation anomalies in California between 1948 and 2016 (39). The first mode of the MCA is fairly robust when analyzing the individual iterations of the LMR (see Materials and Methods).

In the MCA above, the squared covariance fraction quantifies the fraction of squared covariance between two fields represented by a given mode of variability (40), and the values for fraction of variance separately quantify the fraction of total variance represented in each of these fields. These metrics indicate that the first mode of the MCA, discussed above, accounts for 39% of the variance in the U.S. PDSI field, 34% of variance in the joint SST/500-hPa height field, and 83% of the squared covariance between these two fields.

If the MCA is conducted between PDSI and the previous year’s SST and 500-hPa geopotential height anomalies (rather than comparing coincident years), then the patterns are similar to those discussed above (fig. S6). Examining a time-lagged relationship is a good target for analysis, as years are reconstructed individually in the LMR data assimilation and any relationships between different years stem from the proxy records rather than covariances in the model prior. The MCA suggests that the primary link between these fields is that a La Niña/cool PDO pattern is associated with drought in the southwest United States although the correlation of the expansion coefficients is reduced in the lagged case.

As a complement to the previous analysis, which analyzes the full range of hydroclimate variability over the United States, we now explore climate conditions specific to drought states by implementing a drought threshold. “Droughts” are here defined as all years where regional PDSI is more than 1 SD below a 51-year moving window of PDSI, which accounts for any mean state shifts in PDSI values and reductions in variance with the loss of proxies further back in time. To examine the climate patterns associated with drought in different parts of the United States, we calculate drought years for four regions: the northwest, southwest, central, and southeast United States. These regions were chosen in other work to represent regions of greatest statistical drought independence (10), and the mean SST and 500-hPa height anomalies during drought years in each of these regions are shown in Fig. 4 along with the conditional distributions of Nino 3.4, PDO, and AMO values for drought years and nondrought years.

Fig. 4 Climate fields associated with regional drought.

Maps show mean SST (°C), 500-hPa heights (hPa), and PDSI for drought years relative to all years in four U.S. regions: (A) northwest, (B) southwest, (C) central, and (D) southeast United States. The contour interval for 500-hPa heights is 2 hPa, with the thicker line indicating 0. Split violin plots show distributions of annual-mean values of Nino 3.4, PDO, and AMO for drought years (brown) versus nondrought years (green) in each region, with lines indicating the medians and interquartile range. Differences in means between drought years and nondrought years that are not significant according to a resampling test at the 95% level (Materials and Methods) are indicated with an asterisk next to the name of the climate index.

Mean climate states for drought years for each of the four U.S. regions are characterized by La Niña and cool PDO SST patterns, although these patterns are weak when considering droughts in the northwest U.S. region. Droughts in each region also generally correspond with warmer temperatures in much of the North Atlantic and increased 500-hPa heights over the North Pacific, continental United States, and North Atlantic. Analysis of climate indices shows that the Nino 3.4 and PDO indices are significantly lower during drought years compared to nondrought years for at least three of the four regions using a resampling test (Materials and Methods), and AMO is significantly higher in drought years for the central and southeast regions, but mean AMO is not significantly different for the two west coast regions. The analysis of climate indices also reveals a considerable range of values in both drought years and nondrought years, with a large degree of overlap between patterns that correspond to drought years and those that do not. This indicates that while certain climate states (i.e., La Niña and cold PDO) are associated with drought states on average, these relationships only emerge when examining mean state differences among considerable amounts of climate variability.

Similar results emerge when computing linear regressions between the full range of PDSI variations and the surrounding climate fields (not shown). Correlations between PDSI and equatorial Pacific SSTs are strongest for the southwest United States (exceeding −0.6 for SSTs just off the equator, larger than for the Nino 3.4 region itself) and weakest for the northwest United States, again indicating the differing effects of these teleconnections on different regions of the United States. To ensure that these results are not overly determined by covariances in the model prior—which may affect results for coincident years but not the lagged analysis, as mentioned above—an alternate experimental design is explored in note S2.

We use self-organizing maps (SOMs) to explore connections between SSTs, 500-hPa heights, and U.S. drought conditions (fig. S7). SOMs isolate characteristic patterns in a given climate field and identify which years are most represented by each pattern (20, 41). Here, eight SOM patterns are computed, as in (20) (see Methods in that paper for details), from the global SST field over years 1001 to 1925, with the post-1925 years removed to eliminate trends in the SOM patterns because of anthropogenic warming. In addition, detrended 500-hPa geopotential height and PDSI anomalies are composited over the years corresponding to each SOM pattern, revealing the geopotential height and PDSI anomalies that correspond to each SST pattern.

The primary SST pattern that emerges through this SOM analysis corresponds to ENSO. In general, drought years in each region have a higher occurrence of La Niña–like patterns and a lower occurrence of El Niño–like patterns, in agreement with the relationships described above (fig. S7I). This connection appears to be strongest for the southwest U.S. region and weakest for the northwest U.S. region. The primary non-ENSO patterns consist of warmer or cooler SST anomalies overall, with warmer SSTs connected with drought in northern North America and a lack of drought in the southern United States, although this connection is relatively weak (fig. S7, D and G).

Together, these results indicate that La Niña is a noisy predictor of reduced precipitation, but much drought variability appears unrelated to simple ENSO metrics such as Nino 3.4 (fig. S8). Considerable variability exists in the analyzed teleconnection patterns, as seen in the large overlap in climate index values for drought versus nondrought years (Fig. 4). Years with drought in the southwest United States have below average Nino 3.4 75% of the time, and, when considering all years, Nino 3.4 only accounts for 13% of the variance in southwest U.S. PDSI [i.e., coefficient of determination (R2) = 0.13; fig. S8]. Some research has suggested that SST anomalies over longer time periods influence longer term drought (10), but only a weak correlation emerges for decadal means of the present reanalysis (R2 = 0.07 for 1001 to 2000 C.E.).

Even observational datasets [Nino 3.4 (27) and Dai PDSI (30)] reveal considerable variability in the relationship between these two quantities; for the years 1874 to 2000 C.E., R2 between southwest U.S. PDSI and Nino 3.4 is 0.10 when Nino 3.4 is calculated during the previous DJF, 0.10 when Nino 3.4 is calculated over the previous year, and only 0.03 for coincident annual means (fig. S5). Despite this, we primarily analyze coincident years in the reanalysis results because they produce a higher correlation on these time scales, potentially a result of the reanalysis methodology. Slightly higher values can be found for PDSI near the coast of Texas, but the general weakness of these relationships suggest that ENSO, measured by standard metrics such as Nino 3.4, is a rather minor influence on U.S. drought. Stronger connections are revealed when considering the northern Pacific and Atlantic basins as a whole, as seen in the MCA analysis, but a considerable portion of drought variability still appears unrelated to the Nino 3.4 metric alone, possibly suggesting the need for a more comprehensive approach when evaluating the ocean’s influence on U.S. drought.

Drought response to external forcings

According to a recent modeling study using prescribed SSTs (12), SST variations in the global oceans explain 40% of annual-mean precipitation variance in northern Mexico and the southeastern United States but much less in other regions, including the southwest United States. The remaining drought variability is a topic of much interest (12) and may help explain events such as the reduced precipitation in southern California during 2015/2016, which occurred despite a strong El Niño. Along these lines, we now consider the extent to which some drought variability may be driven by external climate forcings, such as greenhouse gases, explosive volcanism, or variations in solar irradiance, which is a question well suited for GCMs. Volcanic eruptions, for instance, reduce global-mean precipitation in the Hadley Centre Coupled Model version 3 (HadCM3) (42) and Coupled Model Intercomparison Project 5 models (43), with the largest modeled precipitation changes taking place in the tropics. Volcanism can also generate abrupt cooling, followed by a recovery of several years (44).

Because climate models can be run with specified forcings, they provide a valuable counterpart to data-driven studies like the LMR. Here, we explore the effects of external forcings on past climate in the Community Earth System Model (CESM)–Last Millennium Ensemble (LME) simulations (45). The CESM-LME consists of a set of transient GCM simulations starting at 850 C.E. and run with changes in one or all of the following: greenhouse gases, volcanic aerosols, solar forcing, orbital forcing, and land use change. In addition to analyzing an ensemble of nine simulations run with all forcings (hereafter called “All” or “fully forced”), we investigate a variety of single- forcing simulations where one forcing varies while all other forcings are set to their 850 C.E. values. These simulations focus on the effects of changes in greenhouse gases (three ensemble members), land use/land cover (three members), Earth’s orbit (three members), solar irradiance (four members), and volcanic forcing (five members) [see table 1 in (45)]. While orbital forcing, greenhouse gases, and land use change have much slower rates of change than the drought variability of interest, they are included here for completeness and to examine whether slow changes in these parameters can affect general drought statistics. These simulations present a useful complement to the LMR reconstruction because they explore climate variations in the presence or absence of certain forcings, which is impossible to fully disentangle in observations alone. In addition, multiple ensemble members allow us to sample many different expressions of internal atmospheric variability, which is important for determining which variations are endogenous to the atmosphere-ocean system and which are exogenous. We examine how different external forcings in these simulations affect PDSI in the southwest United States. Because this analysis spans the years 850 to 1849 C.E., recent anthropogenic changes will not be considered.

The nine fully forced simulations are subjected to identical forcings and differ only in their initial conditions. Because imposed forcings always occur with the same timing and magnitude across these simulations, externally forced responses should exhibit consistent timing across simulations, emerging with averaging, provided that the ensemble size is large enough. Variations that are a function of unforced atmosphere-ocean variability, on the other hand, including variability associated with ENSO or other large-scale teleconnections, need not have consistent timing between ensemble members and tend to cancel out across ensemble members. Put another way, forced responses should emerge as common signals from the otherwise distinct climate variations in each simulation.

We compute PDSI from modeled climate values in the CESM-LME (Materials and Methods). When southwest U.S. PDSI is compared across the nine fully forced simulations, considerable differences are evident, with the mean signal exhibiting relatively small variations (Fig. 5). To quantify similarities between any two simulations, we compute correlations for southwest U.S. PDSI between every pair of fully forced simulations. These correlations have a mean value of 0.02, and no two time series agree with a correlation above 0.09. This lack of consistency indicates that little southwest U.S. drought variability may be explained by external forcing. The largest forced signal relates to explosive volcanism, which tends to produce wetter conditions in the southwest United States after eruptions. This wettening is particularly apparent after the Samalas eruption of 1257 C.E., which is the largest volcanic forcing in these simulations, when all nine ensemble members show positive PDSI regardless of prior conditions (Fig. 5). A more detailed analysis of volcanic responses in these simulations has been presented in past work (46). This volcanically forced signal, however, only accounts for a small part of the overall variability. In other regions of the world, such as northwestern South America and northwestern Africa, volcanic eruptions have a larger relative impact in these simulations (not shown).

Fig. 5 Southwest U.S. PDSI in CESM-LME.

Time series of PDSI averaged over the southwest United States (land within 32°N to 40°N, 125°W to 105°W) in the nine all-forcing simulations, as well as their mean. Years of the 10 largest volcanic forcings are marked; these vertical lines mark the first year of a large volcanic aerosol forcing, so the year listed may not exactly match the year of the actual eruption.

Single-forcing CESM-LME simulations are used to further investigate whether (and to what extent) particular forcings influence southwest U.S. drought. To determine whether much temporal agreement exists between single-forcing experiments and fully forced experiments, we calculate correlations between different sets of simulations (fig. S9). Of the single-forcing experiments, volcanic forcing has the largest correlation with the fully forced simulations, although the median correlation in that case is still only 0.04, suggesting that very little of the total variability is explained by these external forcings.

While external forcing appears to have little influence on the timing of droughts, it is worth investigating whether imposed forcings affect PDSI characteristics in other ways, such as the length or severity of a drought. To prevent long-term trends from overly affecting the estimated variability, we remove a linear trend from each southwest U.S. PDSI time series, and droughts are then calculated as periods beginning with two consecutive years with PDSI below 1 SD relative to the nearest 51-year periods and ending with two consecutive years with PDSI above that threshold, similar to the definition used by Coats and coauthors (47, 48). Using this definition, the average frequency, length, and magnitude of droughts are calculated over the 1000-year interval for each simulation. This analysis is similar to the one performed in other work (16), which showed that a considerable portion of drought variability may be unrelated to SSTs. Comparison of these drought statistics reveals a high degree of similarity across CESM-LME experiments, indicating that the applied external forcings do not have a large impact on the longevity or magnitude of droughts in this region (Fig. 6). This is consistent with past work, which has shown that external forcings are not required to explain the magnitude, spatial, and temporal extent of severe droughts such as those seen in the proxy record (i.e., megadroughts), though these forcings may be necessary to explain the clustering of these droughts during the medieval era (49).

Fig. 6 PDSI statistics in different experiments.

Statistics of annual-mean PDSI in the southwest United States in different CESM-LME simulations, LMR, and NADA. Bar plots show the (A) number of droughts, (B) average drought length, (C) and average drought strength for years 850 to 1849 C.E. Colored bars show the ensemble-mean values for each experiment type, with black dash marks showing the values for each ensemble member. The LMR and NADA results are shown in different colors to call attention to the different methodologies used. GHG, greenhouse gases; LULC, land use/land cover; Orbit, Earth’s orbit; Solar, solar irradiance; Volc, volcanic forcing.

For comparison, statistics of the LMR show southwest U.S. droughts that are generally longer and weaker, a characteristic that can be seen in the PDSI time series (fig. S10). The statistics of southwest U.S. PDSI in NADA, on the other hand, show more intense droughts. These differences are difficult to explain, but some of the disparity may stem from the aforementioned methodological differences between LMR and NADA regarding variance restoration, tree ring detrending, and proxy prewhitening.

In general, external forcings appear to have only minor effects on southwest U.S. drought in the CESM-LME simulations. Volcanic eruptions encourage wetter conditions in the short term, but these forced variations make up only a small portion of the total PDSI variability. Our focus on the years 850 to 1849 C.E. makes it difficult to comment on the effect of anthropogenic forcing during the industrial period, but this analysis suggests that natural forcings have only exerted a minor influence on southwest U.S. drought in the centuries before 1850, when greenhouse gases were more constant.


The LMR constitutes a powerful methodology for creating a physically consistent multivariate climate reconstruction from a diverse array of proxy records. The proxies provide data about specific regions, climate fields, spans of time, and seasons, and we use proxy system models (PSMs) and the covariance structure from a GCM to synthesize these diverse perspectives into a cohesive view of past climate. While more work is needed, analyses show that the LMR climate reconstruction compares well with established datasets for temperature, PDSI, and large-scale climate indices, providing evidence of reconstruction skill. In accord with a wide body of published work, the LMR reconstruction also finds a clear connection between southwest drought and a La Niña SST pattern over the past millennium. This connection emerges as the primary mode of covariability between PDSI, SST, and 500-hPa height fields, even though the analysis methodology (MCA) is not directed to focus specifically on the equatorial Pacific. While this pattern is robust, teleconnections to SST variations appear to explain only a part of U.S. drought variability, leaving the larger portion of drought variability unexplained. Fundamentally, this data assimilation approach presents data-based evidence for the importance of internal atmospheric variability in determining past hydroclimate variability, in agreement with other work (1517, 39, 50).

The present work contains important caveats, however. In particular, the data assimilation methodology relies on GCM output to partly quantify relationships between different variables and locations, as well as provide a first estimate of past climate. Using a model prior like this is necessary, as it provides the framework for synthesizing information from diverse proxies—which differ in their climate sensitivities, locations, seasonal biases, and temporal coverage—into a physically consistent multivariate climate reconstruction. However, model bias in spatial climate covariance patterns does affect the reconstructions. While this work has used alternate analyses (e.g., lagged correlations and alternate experimental designs) to minimize the impact of these potential biases in results, future work should explore this topic in more depth. The use of model priors from several different GCMs, for instance, may help mitigate the effect of biases in any one particular model.

Another area for future improvement is the incorporation of additional proxies into the data assimilation product, particularly from poorly sampled regions and using additional archive types. One data assimilation advance, which is being explored in current and past work (51), is the incorporation of lower-resolution proxies into the data assimilation methodology. While proxies that lack annual temporal resolution will require additional considerations within the data assimilation framework, these proxies can provide information about sparsely sampled regions (such as continental margins, in the case of marine cores) and should more accurately capture low-frequency variations compared to tree rings (52, 53), which are heavily represented in the current data assimilation approach. This could provide additional information about slower climate variations and trends, refining our understanding of climate variations such as drought and potentially making this approach more relevant for studies of past megadroughts [e.g., (54)]. Considering the potential for anthropogenic changes to worsen future droughts in many regions (79), better understanding of the climate dynamics behind drought variations is critical for future planning.

If, as we argue, internal atmospheric variability has been a leading cause of multiyear drought over the Common Era, then this bears unfavorably on the prospect for forecasting these droughts. This may have been at play in the 2011–2017 California drought, which was, by some measures, the most severe Californian drought of the past 1200 years (2). Much drought relief was expected from the 2015/2016 El Niño, which rivaled in magnitude the extreme El Niño events of 1982/83 and 1997/98 (55). However, while those previous events brought abundant rainfall to California (39), the 2015/2016 event produced average rainfall throughout most of California (39, 55), defying the expected teleconnection pattern in this region (39). Consequently, this event failed to end the prolonged drought, which had persisted since 2011/12. Southern California had to wait until the following year, characterized by mild La Niña conditions, to receive enough rainfall to end the drought (56). Our results suggest that this situation may have been a common occurrence throughout the past millennium, making current limitations in interannual drought predictability especially important in a warming climate with greater evaporative demand for moisture (2, 57).


Experimental design: Paleoclimate data assimilation

Paleoclimate data assimilation offers a powerful approach for synthesizing a vast array of proxy observations (here, thousands of records) with the aid of a model’s climate covariance structure. In particular, the LMR is a data assimilation approach that uses information from proxy records and output from a GCM to estimate climate variability over the past two millennia. An earlier version of this method is described in past work (23), and updates to the methodology are described in a recent paper (24). The data assimilation methodology is composed of four primary components: (i) GCM output, which serves as a “first guess” at the range of possible climate states and quantifies covariances within the climate system; (ii) PSMs, which relate the model quantities to proxy quantities; (iii) proxy records, which provide the temporal information for the reconstruction; and (iv) a Kalman filter, which is used to perform the data assimilation.

The methodology works by first randomly selecting a collection of annual-mean climate states from the output of an existing GCM simulation. Here, we use 100 years from the Community Climate System Model 4 (CCSM4) last millennium simulation (58). These randomly selected model states (i.e., the prior) are initially identical for every year of the assimilation, serving as the first guess of the real climate state for any given year. In other words, before any assimilation takes place, the real climate is suspected to be somewhere in this range of modeled climate states. The LMR is then run for 20 iterations; within a given iteration, the same prior is used for every year of the reconstruction, so the model provides no temporal information to the reconstruction. In addition to providing an initial range of plausible climate states, the 100-member prior is also used to quantify the covariances within the climate system, which forms the mathematical scaffolding that relates climate variations at different locations and between different fields.

To perform the data assimilation, proxy and model quantities must be compared in the same units, so PSMs are needed. Here, relationships between proxy quantities and climate variables are computed by regressing proxy records onto instrumental fields over the period 1880 to 2015 C.E. For all proxies except tree rings, a linear regression is computed against temperature [GISTEMP version 4 (59)]. For tree rings, a bivariate regression is computed against both temperature (GISTEMP) and precipitation [GPCC (60)]. Regarding seasonality, proxy records are not assumed to record annual-mean quantities. Instead, proxy records are regressed onto climate quantities averaged over the entire year as well as multiple subsets of the year: summer, winter, and four different two-season half-year periods (as well as the season specified in the proxy metadata if different from the previously mentioned seasons). The averaging windows that produce the best regression between instrumental data and proxy records are used in the data assimilation, and this is determined separately for each proxy (24). Tree ring width proxies, which are regressed onto both temperature and moisture, are allowed to have different seasonal sensitivities for temperature and moisture (24). Linear regressions are simpler than a process-based model, but they provide good results that can provide a baseline for more complex PSMs in the future (24).

The proxy network used in the present work includes records from three sources: the PAGES2k v2.0.0 database (61), several thousand tree chronologies compiled by Breitenmoser et al. (62), and a selection of other proxy records from the National Centers for Environmental Information (formerly National Climatic Data Center) paleoclimate archives. The proxies are available at Zenodo (63), and the non-PAGES2k records are described in a recent data report (64). To be assimilated, records must be at least annual in resolution and must have at least 25 years of overlap with the instrumental records. A total of 2787 records are assimilated in the present reanalysis, the spatial coverage of which is shown in Fig. 1.

For each year of the reanalysis, the prior is used as a starting point, and the climate state is updated through assimilation of annually resolved proxy records one by one via a Kalman filter. The Kalman filter compares each proxy value against an estimate of the proxy value computed from the model prior and then adjusts the climate state to produce a better fit for the given year. Because the model prior quantifies the climate covariance structure (between locations as well as between climate variables), it provides the mathematical framework for updating more distant locations as well as a variety of climate variables in a uniform framework. In general, climate is reconstructed through comparison with both local and remote proxies. Locations closest to proxies, as well as variables that are most closely related to proxy measurements (23), are expected to be better informed by the proxy network, while other locations and climate fields rely more on model covariances. For example, previous research has shown that the LMR has higher skill in reconstructing surface temperature than 500-hPa height (23), although the skill of both has been improved with recent methodological innovations (24). These qualifications should be kept in mind when interpreting results.

The ability to reconstruct multiple variables has clear benefits and facilitates the analysis of climate teleconnections over an extended period, with some qualification (note S2). In the LMR, the proxy records provide temporal information and some spatial information (by making use of multiple records in space), while the covariance structure of the GCM prior is used to propagate information between locations and between climate fields. A localization radius is used to ensure that proxies cannot influence the climate farther than 25,000 km from their location, a value that was chosen to produce the expected variance characteristics in the reconstructed temperature [see table 1 in (24)]. Further details of this methodology are explained in other work (23, 24).

In the present analysis, 20 iterations of the LMR were run. Each iteration uses a different random selection of 100 model years for the prior and a different random selection of 75% of the proxies for assimilation. Variety in the priors and assimilated proxies helps sample uncertainty in the results. All climate fields are output as annual quantities averaged from January to December. Exact methodological choices are explored in past work (24), and the data have been made available as the v2.0 release of the LMR dataset (see “Data and materials availability”).

The number of annually resolved proxy records used in the reanalysis decreases back in time (Fig. 1B), and larger differences emerge between the LMR results and the NADA for the first millennium compared to the second, so the analysis in this paper focuses on years 1001 to 2000 C.E. The tree proxies used in the LMR and the NADA have considerable overlap (fig. S3), although the two approaches have numerous differences: The methodologies are distinct, the LMR uses additional proxy types, and the methods used to remove tree ring growth curves are likely different as well, among other differences. Still, because trees are the most numerous proxy used in this study, both over the United States and globally, this overlap in data sources should be considered when comparing drought in the LMR and the NADA.

Using this multivariate reconstruction of past climate, relationships in the climate system can be explored over an extended period of time, with the qualifications mentioned in note S2. Pure modeling studies, which are used for exploring possible future drought changes (7, 8), are deficient in modeling some aspects of drought variability; models may have inadequate low-frequency hydroclimate variability (65, 66), although at least one study finds a similar number of long droughts in the southwest United States in models as compared to NADA (48) and the present analysis finds similar values between NADA and the CESM-LME simulations. In addition, because reanalysis uses real proxy data, this method can provide insight on actual past droughts.

PDSI calculations

To reconstruct the PDSI over the past two millennia, PDSI values are first calculated from quantities in the CCSM4 last millennium simulation, which is used as the model prior. This was done using the Penman-Monteith equation for potential evapotranspiration and monthly climate model output of 2-m air temperature, precipitation, vapor pressure, surface pressure, net surface radiation, and surface wind (estimated from 10 m down to 2 m using the wind profile power law). The computations of PDSI were carried out using the MATLAB code from Jacobi et al. (67), which produces the standard formulation of PDSI as opposed to self-calibrating versions [e.g., (68)]. Once PDSI is calculated in the model prior, it can be included in the LMR data assimilation to calculate proxy-informed PDSI values over the past two millennia. The same method was used to calculate PDSI in the CESM-LME simulations.

Maximum covariance analysis

To examine how U.S. drought covaries with large-scale patterns of the surrounding climate system, an MCA [also called singular value decomposition (40)] is used to isolate the mode that explains the largest amount of covariance between two fields. Here, one field is PDSI over the United States and the second field is a concatenation of SSTs and 500-hPa heights over a larger region (see the regions displayed in Fig. 3). To ensure that neither the temperature nor the 500-hPa height anomalies dominate the second term of the MCA, all climate anomalies have been standardized by the mean and SD over their entire regions, an alternate approach mentioned in past work (40). This MCA analysis is conducted on variables on their reconstructed 2° latitude-longitude grid, with spatial weighting applied. To illustrate common patterns in SST/500-hPa height and their impact on drought, figures show the homogeneous map of SST/500-hPa height and the heterogeneous map of PDSI. Maps are scaled to have standardized expansion coefficients before plotting.

Resampling test for significance

The split violin plots in Fig. 4 show the range of annual-mean climate indices (Nino 3.4, PDO, and AMO) for years when different regions of the United States were experiencing drought or nondrought. To see whether the mean of each climate index was significantly different during drought and nondrought years for each of these cases, a resampling test was done. In this test, the continuous spans of years spent in drought or nondrought were identified for each case, and then analogous sets of values were randomly selected from the whole time series of the climate index, using sampling with replacement. This was repeated 1000 times, and the original drought versus nondrought climate index anomaly was compared against the anomalies in these 1000 randomly sampled cases. Cases in which the original difference was outside of the 2.5th or 97.5th percentile were deemed significant.


Supplementary material for this article is available at

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.


Acknowledgments: We would like to extend thanks to E. Cook for invaluable information about tree ring proxies and R. Tardif for help with the LMR code. Funding: This work was supported by the National Oceanic and Atmospheric Administration (grants NA14OAR4310175 and NA14OAR4310176) and the National Science Foundation (grants NSF AGS-1702423 and NSF AGS-1805490). M.P.E. was supported by the University of Southern California and Northern Arizona University, and computational resources were also provided by both universities. LDEO contribution number 8427. Author contributions: M.P.E. provided most of the analysis and writing of the paper. J.E.-G. provided additional analysis and experimental design. N.S. provided PDSI calculations and discussion of drought considerations. G.J.H. and E.J.S. provided much of the methodological design of the LMR, as well as feedback on the paper. All authors contributed to the writing. Competing interests: The authors declare that they have no competing interests. Data and materials availability: The LMR v2.0 dataset, as well as documentation, proxies, and other sample data to run additional reanalyses, has been made publicly available at and The reanalysis data consist of spatially reconstructed climate fields for a variety of variables over the years 1 BCE to 2000 C.E. The source code for the LMR is available at

Stay Connected to Science Advances

Navigate This Article