Pancontinental droughts in North America, or droughts that simultaneously affect a large percentage of the geographically and climatically distinct regions of the continent, present significant on-the-ground management challenges and, as such, are an important target for scientific research. The methodology of paleoclimate-model data comparisons is used herein to provide a more comprehensive understanding of pancontinental drought dynamics. Models are found to simulate pancontinental drought with the frequency and spatial patterns exhibited by the paleoclimate record. They do not, however, agree on the modes of atmosphere–ocean variability that produce pancontinental droughts because simulated El Niño–Southern Oscillation (ENSO), Pacific decadal oscillation (PDO), and Atlantic multidecadal oscillation (AMO) dynamics, and their teleconnections to North America, are different between models and observations. Despite these dynamical differences, models are able to reproduce large-magnitude centennial-scale variability in the frequency of pancontinental drought occurrence—an important feature of the paleoclimate record. These changes do not appear to be tied to exogenous forcing, suggesting that simulated internal hydroclimate variability on these time scales is large in magnitude. Results clarify our understanding of the dynamics that produce real-world pancontinental droughts while assessing the ability of models to accurately characterize future drought risks.
North America spans three countries, nearly 10 million square miles, and multiple geographical regions with distinct climates, seasonalities, and connections to the wider atmosphere–land–ocean system. While each of the distinct geographic and climatic regions in North America (e.g., the southwest, 32°–40°N, 125°–105°W; versus the central plains, 34°–46°N, 102°–92°W) experience recurrent drought (e.g., Cook et al. 2007; McCabe et al. 2004; Nigam et al. 2011; Schubert et al. 2004a,b; Seager et al. 2005; Seager and Hoerling 2014), drought simultaneously affecting multiple regions (hereinafter pancontinental drought) has also been shown to be a consistent, albeit infrequent, occurrence over the last millennium (Cook et al. 2014b). The most recent pancontinental drought event occurred in 2012 (Hoerling et al. 2014) with 62% of the contiguous United States being classified as moderately or extremely dry (NCDC 2013a). Events such as these pose significant management challenges because of the simultaneous drought impacts across regions with distinct water resource constraints (e.g., irrigation from rivers versus groundwater), ecosystems (e.g., forests and grasslands), and crops. Understanding the dynamics that drive pancontinental drought is therefore critical, but there are two fundamental reasons why a comprehensive characterization of pancontinental droughts, and their causes, proves challenging. First, regional hydroclimate is characterized by distinct atmosphere–ocean dynamics, for instance, southwestern hydroclimate is controlled primarily by winter precipitation variability coupled to the tropical Pacific (e.g., Herweijer et al. 2006; Seager et al. 2008; Schubert et al. 2009), while the Great Plains has predominantly summer hydroclimate variability that is driven by the tropical and subtropical Atlantic (e.g., Sutton and Hodson 2005; Kushnir et al. 2010) in addition to the tropical Pacific (Seager et al. 2005). Second, the relative rarity of pancontinental drought and the short (~150 yr) observational record means that there are few events by which to diagnose how relatively distinct regional hydroclimate dynamics can combine to produce pancontinental drought. The limitations of the short observational record were partially addressed by Cook et al. (2014b), who employ a tree-ring-based reconstruction of North American hydroclimate over the last millennium (the North American Drought Atlas, NADA; Cook et al. 2007) to better characterize the statistics of pancontinental drought occurrence. The paleoclimate record, however, does not provide a complete picture of the atmosphere–ocean state during pancontinental droughts. Extending the analyses of Cook et al. (2014b) to the model space using the methodology of paleoclimate-model data comparisons (e.g., Anchukaitis et al. 2012; Ault et al. 2013, 2014; Coats et al. 2013a, 2015; Fernández-Donado et al. 2012; Phipps et al. 2013; Schmidt et al. 2013) can provide additional characterizations of past pancontinental droughts (albeit with the caveat of model bias) and potentially help clarify our understanding of the dynamics that produce these features. Additionally, the frequency of pancontinental drought occurrence is likely to increase in the future as much of North America is expected to dry over the coming century (e.g., Seager et al. 2013; Cook et al. 2014a; Maloney et al. 2014). Because atmosphere–ocean general circulation models (AOGCMs) are used to make these future hydroclimate projections, determining if AOCCMs are capable of reproducing the statistics of past pancontinental drought occurrence, and as a consequence of the correct dynamical drivers, is necessary to assess whether state-of-the-art AOGCMs can accurately constrain future drought risks.
We build off the proxy analyses of Cook et al. (2014b) by employing coupled model simulations to further assess the dynamics of pancontinental drought in North America. We use the same gridded tree-ring-based reconstruction of hydroclimate variability (Cook et al. 2007) over North America from 1000 to 2005 common era (CE) used by Cook et al. (2014b) and compare it to pancontinental drought statistics and dynamics in forced transient simulations of the last millennium (850–1850 CE) and the historical interval (1850–2005 CE), as well as in 500-yr control simulations from the same models—all from the Coupled and Paleo Model Intercomparison Projects Phases 5 and 3 (CMIP5/PMIP3; Taylor et al. 2012; Schmidt et al. 2011). Four fundamental questions are addressed:
Do models simulate the major modes of atmosphere–ocean variability that impact North American hydroclimate (section 3a)?
Are models able to reproduce the statistics of pancontinental drought occurrence (section 3b)?
Do the models have centennial-scale variability in the occurrence of pancontinental droughts, and if so, is this driven by forced or internal variability (section 3c)?
What are the simulated atmosphere–ocean dynamics that drive pancontinental droughts (section 3d)?
2. Data and methods
a. Model and paleoclimate inputs
All employed model output is from the CMIP5/PMIP3 archive (Table 1). We use six “Last Millennium” (LM) simulations spanning the period 850–1850 CE that are forced with reconstructed time-varying exogenous forcings (Schmidt et al. 2011). These simulations have been appended to CMIP5 historical runs that span the period 1850–2005 CE to produce a model record from 850 to 2005 CE Although these simulations are not continuous, both of the paired simulations (historical and LM) have been generated using the same model configurations and resolutions. Consequently, if the simulations have no drift, the discontinuity at 1850 CE should fall within the range of simulated climate variability. A large temperature drift in the MIROC LM simulation (Sueyoshi et al. 2013) likely violates this assumption, while a drift in the early centuries of the GISS LM simulation (Bothe et al. 2013) is likely to have less of an impact. While model drift undoubtedly impacts the hydroclimate variables assessed in this study, the effects are presumed to be moderate given the absence of drift in precipitation (Sen Gupta et al. 2013). The 500-yr control simulations with constant preindustrial forcings (also from CMIP5) were additionally analyzed to aid in the interpretation of the LM model results. All model output has been regridded onto a common 2.5° × 2.5° latitude–longitude grid to allow for homogenous comparisons (this represents a coarsening of the model resolution for four out of six models).
For each model simulation we calculate the Palmer drought severity index (PDSI; Palmer 1965). PDSI is an offline estimate of soil moisture balance, and has been established as a robust estimator of soil moisture variability that compares well with other soil moisture metrics (e.g., the standardized precipitation evapotranspiration index, SPEI; Vicente Serrano et al. 2010; Cook et al. 2014a) and inherent model soil moisture (Cook et al. 2014a; Smerdon et al. 2015). PDSI is calculated from supply via precipitation and losses due to evapotranspiration (ET). In this case, ET is calculated by means of scaling potential evapotranspiration (PET), estimated from surface net radiation (RNET), by a fixed beta function meant to represent vegetative controls on transpiration. There are multiple ways to compute PET, with the ideal method being Penman–Monteith (PM), which includes the effect of the vapor pressure deficit along with the impact of RNET (a more detailed treatment of PM PDSI can be found in Sheffield et al. 2012; Cook et al. 2014a; Smerdon et al. 2015). The necessary model fields needed to compute PM PET were only available for three out of the six analyzed LM simulations and as a consequence RNET instead has been used to calculate PET by assuming that RNET is exactly balanced by latent heat through ET (with sensible heat flux equal to zero). Importantly, PDSI calculated with PET estimated from RNET compares well with PM PDSI on both interannual and decadal time scales (Coats et al. 2015).
For the analyses herein, model PDSI is derived on an even 2.5° × 2.5° latitude–longitude grid. At each grid point, PDSI was calculated and then standardized against an instrumental normalization period (1931–90 CE) for the forced simulations, and the full 500-yr period for the control simulations. The instrumental normalization period is the same time interval used by the National Oceanographic and Atmospheric Administration for normalization of their PDSI calculations, which were subsequently used as the target PDSI for the paleoclimate reconstructions described below. Soil moisture capacities were specified as 25.4 and 127 mm in the top and bottom layers of the PDSI calculation, respectively. The PDSI was averaged over June–August (JJA) to produce a single average for each year; hereinafter any mention of PDSI will be with regard to the JJA average values. PDSI with an absolute value greater than 10 was removed by replacement with the average PDSI of the eight neighboring grid points at that time step as a means of removing unrealistically anomalous PDSI values. This method for removing data errors is consistent with that used by van der Schrier et al. (2011) in the calculation of their observed PDSI dataset. For a more detailed analysis of the PDSI data used herein, see Coats et al. (2015). Additionally, for a full treatment of the impact of inherent model resolution on the simulation of hydroclimate over the study region see Sheffield et al. (2013) and Langford et al. (2014), and for a discussion of how these resolution biases may impact the PDSI data described above, and used herein, see Coats et al. (2015).
Reconstructed PDSI data are from an updated version of NADA (version 2a), with improved spatial coverage and resolution, the full details of which can be found in Cook et al. (2014b). The data are reconstructed on a 0.5° × 0.5° latitude–longitude grid of JJA-average PDSI values for the United States, as well as parts of Canada and northern Mexico.
Observed monthly PDSI data are from a global dataset for the period 1870–2012 CE (Dai et al. 2004). This dataset was derived on an even 2.5° × 2.5° latitude–longitude grid using observed precipitation and temperature data as inputs. For the analyses herein the JJA monthly PDSI values have been averaged to create a single PDSI anomaly for each year. Additionally, only the period after 1950 will be employed herein, as this is the period over which the full North American PDSI grid is available.
JJA PDSI reflects hydroclimate conditions for the past 12–18 months because of persistence built into the PDSI calculation and, as a consequence, will integrate both the winter and summer hydroclimate states in any given year. Nevertheless, tree-ring-reconstructed PDSI has been shown to predominately reflect winter season precipitation variability over the western United States, and both winter and summer precipitation variability over the central and eastern United States (St. George et al. 2010), a characteristic that may be shared by the models. Given uncertainty with respect to the seasonal hydroclimate influences on PDSI, no attempt is made herein to analyze the seasonality of pancontinental drought dynamics.
b. Climate indices
All climate indices are calculated using either the inherent surface temperature output from the control simulations or observations from the National Oceanic and Atmospheric Administration (NOAA) extended reconstructed sea surface temperature (SST) dataset (Smith and Reynolds 2003). The Niño-3.4 index was calculated by averaging December–February (DJF) SST over the region 5°S–5°N, 170°–120°W. The Pacific decadal oscillation (PDO) was evaluated by calculating the EOFs of SST over the extratropical Pacific basin (20°–90°N, 60°W–75°E), and subsequently using the DJF average of the principle component time series corresponding to the EOF that best matches the first EOF of the observed record [which is defined as the observed PDO following Mantua et al. (1997)] when comparing the EOF patterns over the full (Northern and Southern Hemisphere) Pacific basin. This was the first EOF for all models except the BCC model, in which the second EOF was more representative of the observed PDO pattern. In all cases the model patterns exhibit hemispheric symmetry and a tropical expression, as is characteristic of the observed PDO. The Atlantic multidecadal oscillation (AMO) was calculated by averaging JJA Atlantic SSTs over the region 0°–60°N, 80°W–0° and then subtracting the global JJA SST average between 60°S and 60°N (following Enfield et al. 2001). All correlations between simulated PDSI and climate indices will be calculated using the JJA gridpoint PDSI and either the preceding DJF average ENSO or PDO indices or the contemporaneous JJA average AMO index.
Following Cook et al. (2014b), the regional boundaries used in this paper are the southwest (SW), 32°–40°N, 125°–105°W; the central plains (CP), 34°–46°N, 102°–92°W; the northwest (NW), 42°–50°N, 125°–110°W; and the southeast (SE), 30°–39°N, 92°–75°W (see regions in Fig. 1). Similar to results shown for the NADA in Cook et al. (2014b), the designated regions in the models do not have climate variability that is completely independent. Nevertheless, correlation maps between the four regional-average time series and gridpoint PDSI indicate that the hydroclimate within each region is homogenous and that variability between individual regions is largely independent (not shown).
Again following Cook et al. (2014b), droughts are characterized to have occurred in the regional-mean time series when PDSI falls to a value of −0.5 or lower in any individual year. Pancontinental droughts are then defined as occurring when any three [SW, CP, and SE (hereinafter SW+CP+SE); SW, CP, and NW (hereinafter SW+CP+NW); SW, NW, and SE (hereinafter SW+NW+SE); or CP, NW, and SE (hereinafter CP+NW+SE)] or all four [SW, CP, NW, and SE (hereinafter SW+CP+NW+SE)] of the regional mean time series simultaneously have PDSI values of −0.5 or lower in the same year. By this definition, the four-region droughts will overlap with, and also be counted as, three-region droughts. For some of the analyses noted in the results section the three- and four-region droughts were treated as distinct events.
To determine the significance of the pancontinental drought teleconnections (Niño-3.4, PDO, and AMO), climate index composites are computed for all years that exhibit pancontinental drought. A 5000-member ensemble resampling of the climate indices is then performed to generate 90th percentile confidence limits. For example, for the full CCSM control simulation there are 83 yr that qualify as SW+CP+SE droughts. An average of the Niño-3.4 values for these 83 yr gives a composite Niño-3.4 anomaly associated with these events. We then draw 83 random years from the Niño-3.4 time series and average them, repeating this process 5000 times. If the original composited Niño-3.4 anomaly exceeds the 90th percentile thresholds of the ensemble resampling, the association between pancontinental drought and the dynamic teleconnection is characterized as significant.
3. Results and discussion
a. Model dynamics
To diagnose the dynamical drivers of pancontinental drought, the Niño-3.4, PDO, and AMO indices are used to assess the relationships between pancontinental drought and the major modes of atmosphere–ocean variability that impact North American hydroclimate. In this section, we investigate the model expression of the ENSO, PDO, and AMO to determine if the simulated modes of variability are a reasonable representation of real-world dynamics. For ENSO, negative values of the Niño-3.4 index, or La Niña conditions, have been associated with drought in the SW, southern CP, and the SE. The PDO may not be fully separable from ENSO; the two modes of variability are negatively correlated (Newman et al. 2003) and therefore have similar spatial expressions. Nevertheless, the PDO has been shown to have important hydroclimate impacts over North America (e.g., McCabe and Dettinger 2002; McCabe et al. 2004, 2008). A positive AMO, with warm Atlantic SSTs, is associated with drying in the CP and SE (Kushnir et al. 2010; McCabe et al. 2004, 2008; Nigam et al. 2011) and, unlike, the ENSO and PDO has considerable persistence up to multidecadal time scales.
Figure 1 shows the teleconnection patterns calculated as the correlation between the Niño-3.4, PDO, and AMO indices and PDSI over North America. These patterns have been calculated for both the NADA and an observed PDSI dataset (Dai et al. 2004) during the overlapping period with the observed SST dataset (1854–2005 CE and 1950–2005 CE, respectively), and for a sliding 152-yr window (the length of the observed/reconstruction overlap) across the full control model simulations. The model pattern in Fig. 1 plots the 152-yr period in which the simulated teleconnection pattern best represents the pattern between the observed SST dataset and the reconstruction (hereafter observation-to-reconstruction), as determined by the maximum centered pattern correlation statistic (CPCS; Santer et al. 1995) between the two fields. The CPCS is a quantitative measure of the similarity of the simulated and observation-to-reconstruction teleconnection patterns, with the range in the CPCS for all of the 152-yr periods in the models being a measure of the stationarity of that simulated teleconnection (Fig. 1, middle). Additionally, the bottom panel of Fig. 1 shows the strength of the ENSO, PDO, and AMO teleconnections. To do so, the sum of the squared teleconnection correlation coefficients was calculated for each of the model segments and the range in these values was then plotted as a boxplot. For comparison, the sum of the squared teleconnection correlation coefficients was also calculated for the full observation-to-reconstruction and observational records. This analysis was limited to the grid points common to each dataset over the plotted North American domain in the top panels of Fig. 1.
The reconstruction-to-observation ENSO and AMO dynamics are largely characteristic of those in the observed PDSI dataset with CPCS values between the patterns of 0.75 and 0.50 and of nearly equal strengths. The PDO teleconnection, however, is significantly weaker in the reconstruction, despite having a similar spatial pattern (CPCS of 0.52). While this may suggest a deficiency in the reconstruction, it is more likely indicative of an inconsistent impact of the PDO over the much longer reconstructed record (152 versus 56 yr for the observed PDSI).
The models are able to simulate 152-yr periods that have a realistic ENSO teleconnection pattern to North America (with the exception of BCC); however the strength of this teleconnection varies greatly, with CCSM having far too strong of a teleconnection and BCC and GISS having an ENSO teleconnection that is too weak. The stationarity of this teleconnection, likewise, varies significantly between the models; the CCSM ENSO teleconnection, for instance, is highly stationary while the ENSO teleconnections in BCC and MIROC are highly nonstationary. The AMO and PDO teleconnections in the models are much less realistic, with none of the models simulating a 152-yr period with a CPCS value over 0.6 and the teleconnection strength being weaker than the observations for each model and both modes. These teleconnections are also nonstationary, with the largest CPCS range occurring in MIROC for the PDO and BCC for the AMO, but with a CPCS range of at least 0.4 for all of the models.
Figures 2–4 show the SST spatial pattern and autocorrelation of the ENSO, PDO, and AMO for the full model simulations and observations. To assess the model skill in reproducing the observed ENSO, PDO, and AMO spatial patterns, Fig. 5 shows the range in the CPCS between simulated spatial patterns calculated for a sliding 152-yr window (length of the observed SST dataset) across the full control model simulations and the observed spatial patterns. Finally, Fig. 6 shows the variance of, or variance explained by, these dynamic modes. As suggested by the teleconnections, models are generally successful at simulating a reasonable ENSO spatial pattern (with the highest pattern correlation values of the three modes; see Fig. 5), although the SST anomalies extend too far west in all of the models (Fig. 2). Additionally, the simulated ENSO autocorrelation structures are largely characteristic of the observations, with oscillatory behavior that varies between negative and positive. This oscillation looks to be most realistic in CCSM, IPSL, and MPI, with a cycle of variability that is too short and regular in BCC and GISS. The magnitude of the ENSO variability is not consistent across the models, with CCSM having too much variance, while GISS, IPSL, and MIROC have too little, as compared to the observations (Fig. 6). Interestingly, the simulated PDO and AMO behavior is largely consistent between the models with a highly nonstationary, but at times realistic, PDO pattern (every model simulates a PDO pattern CPCS of at least 0.6), and a more stationary, but generally unrealistic, AMO pattern (Fig. 5). In both cases, however, the PDO and AMO patterns in the models are less characteristic of the observed patterns than for ENSO. Additionally, the BCC and GISS models fail at simulating the magnitude of the observed tropical expression of the PDO, with all of the models overestimating the high-latitude North Pacific expression of the PDO relative to the expression in the tropics. This is critical because the PDO forcing of hydroclimate variability has been shown to originate in the tropical Pacific (Seager 2015).
The persistence characteristics of the PDO and AMO are plotted for both the forced and control simulations in the right-hand panels of Figs. 3 and 4, respectively. While the models have a reasonable PDO autocorrelation structure (with the exception of BCC and the forced GISS and IPSL simulations—each having too much persistence), with the exception of CCSM and GISS they struggle at simulating the AMO with enough persistence [this is consistent with the behavior of the CMIP3 model ensemble (Ting et al. 2011)]. This lack of persistence suggests that models will have difficulty in simulating the observed drought recurrence interval (or alternatively drought persistence) in the CP and SE regions, which are tightly coupled to the AMO in the real world (Kushnir et al. 2010; McCabe et al. 2004, 2008; Nigam et al. 2011; Ting et al. 2011). Finally, the magnitude of the PDO and AMO will partially determine the impact of these modes of variability, relative to the impact of ENSO and purely atmospheric variability, on North American hydroclimate. There is a large intermodel spread in the variance or variance explained by both modes, with CCSM and MIROC having too much and BCC, GISS, IPSL, and MPI having too little PDO variability compared to observations, and BCC, CCSM, GISS, and MPI having less AMO variability than observed (Fig. 6).
In aggregate, the models exhibit different teleconnections between the oceanic boundary conditions and North America, with no individual model matching the observed atmosphere–ocean dynamics particularly well. Together, this suggests that models are not likely to agree on the modes of atmosphere–ocean variability that are associated with pancontinental drought. Additionally, the models have a stronger and more realistic ENSO and associated teleconnections, as compared to the PDO and AMO and are therefore expected to be more successful at simulating the pancontinental drought dynamics associated with this mode of variability. Nevertheless, the simulated dynamical relationships are largely nonstationary, and the observed dynamics have been inferred from the short 152-yr instrumental interval. The observed dynamics, thus, may themselves be time variable or inadequately characterized (particularly given the large persistence and consequently the small number of degrees of freedom for the PDO and AMO). It is therefore difficult to attribute differences between the models and observations as solely associated with the model characteristics, as opposed to some combination of model misrepresentations and poorly characterized observational teleconnections due to undersampling of low-frequency modes and nonstationarity.
b. Pancontinental drought occurrence
Figure 7 shows the drought recurrence interval for the individual geographic regions in both the forced and control simulations from the models (dark and light bars, respectively) and the NADA. Models, in general, are able to simulate the correct recurrence interval for drought in each of the regions. The model ensemble, however, slightly overestimates the occurrence of SW and NW drought, and underestimates the occurrence of drought in the CP and SE. This model behavior may be suggestive of more realistic, and in some cases overactive, ENSO variability and teleconnections (e.g., CCSM) relative to other modes of coupled atmosphere–ocean variability, because ENSO-driven hydroclimate variability tends to load heavily on the western United States and thus predominantly affects the SW and NW regions.
The model ensemble is also largely successful at simulating the pancontinental drought recurrence intervals characterized by the NADA (again for both forced and control simulations; see Fig. 8). Taken individually, however, the models appear split into two categories, with CCSM, IPSL, and MPI slightly underestimating the recurrence interval of pancontinental drought of all types and GISS and MIROC overestimating the recurrence interval of these droughts by a much larger margin (BCC has realistic recurrence intervals for three of the five pancontinental drought types). Nevertheless, the spread of the model ensemble encompasses the pancontinental drought recurrence interval of the NADA for each drought type. Furthermore, each model is individually successful at capturing the relative occurrence of the different types of pancontinental drought, for instance, the SW+CP+SE combination being the most common and the SW+CP+NW+SE combination being the least common.
The recurrence intervals for the droughts in the individual regions and for pancontinental droughts are not consistently different for the forced and control simulations (Figs. 7 and 8). This suggests that the simulated pancontinental drought dynamics are not dependent on the exogenous forcing and, as such, provide confidence in the use of control simulations to assess simulated atmosphere–ocean variability and its connection to pancontinental drought (sections 3a and 3d). These findings are consistent with previous work that has specifically focused on the dynamics of persistent droughts in the southwest (Coats et al. 2013, 2015).
The composite PDSI patterns of each pancontinental drought type for the control model simulations and reconstruction are plotted in Fig. 9. The bottom panel of Fig. 9 shows the CPCS calculated between the composite PDSI pattern from the NADA and each individual pattern of that pancontinental drought type. The range in CPCS is thus a measure of the consistency of individual drought patterns (small range being more consistent), with the average magnitude of the CPCS values for each model being indicative of how well the model composite matches the NADA composite. These values have been calculated separately for the three- and four-region droughts (unlike in previous analyses). The composite model patterns compare well with the NADA composite for all but the SW+NW+SE droughts, which are the least common of the pancontinental drought types. The CPCS range in the models and the NADA are likewise consistent for all but the SW+NW+SE combination. Nevertheless, it is important to note that for all drought types there is a large CPCS range in both the models and the NADA. This indicates that individual pancontinental droughts can have different spatial patterns. Pancontinental droughts driven by a consistent dynamical driver might be expected to have a consistent pattern. If this is the case, the large CPCS range may then suggest that multiple dynamical drivers are capable of producing each type of pancontinental drought. Equally likely, however, is that pancontinental droughts are influenced not just by SST variations but also by internal atmospheric variability that can create different spatial patterns, as was argued by Hoerling et al. (2014) for the 2012 drought. A third possibility is that a large CPCS range would also be expected if pancontinental droughts were driven by consistent dynamical drivers but with teleconnection dynamics that are variable through time. The question of pancontinental drought dynamical drivers will be addressed in section 3d.
c. Centennial variability in pancontinental drought occurrence
Perhaps the starkest characteristic of the NADA drought record is the centennial-scale variability in the number of pancontinental droughts, punctuated by an increased rate of occurrence in the Medieval Climate Anomaly (MCA) relative to the Little Ice Age (LIA; Cook et al. 2014b).
The relative timing of hydroclimate change in the models—if such changes are present in the CMIP5 models—and the NADA is of particular interest because exogenous forcing may, or may not, have played a role in driving the MCA-to-LIA transition (e.g., Mann et al. 2009; González-Rouco et al. 2011; Goosse et al. 2012). If, in fact, radiative forcing produced this transition, it would be relevant to our understanding of current and future radiatively forced climate change. Because the models are driven with similar forcing series (see Schmidt et al. 2011), a strong role for exogenous forcing in driving periods with increased pancontinental drought frequency should lead to these periods being contemporaneous in time across the model simulations. The role of radiative forcing in driving the variability of pancontinental drought frequency on centennial time scales can therefore be tested to potentially better understand the origin of this variability in the NADA. It must be noted, however, that the CMIP5 models have different climate sensitivities and, in particular, different parameterizations of land surface and aerosol processes that may drive compensating feedbacks and mask the model response to external forcing. If these differences are large, they would impact our ability to test the hypothesized role of forcing as posed above, a possibility that is outside the scope of this paper.
The number of pancontinental droughts for each century relative to the mean number of droughts per century between 1000 and 2000 CE is plotted for the forced model simulations and the reconstruction in Fig. 10. For the NADA, the increase in the number of pancontinental droughts, relative to mean conditions, during the twelfth century averages to 60%, with a maximum increase of 75% for the SW+CP+NW drought type and a minimum of 40% for the SW+NW+SE type. This period of increased drought frequency does not appear to be captured by the models, nor do the individual models agree on the timing of hydroclimate change, suggesting that these changes are not tied in any coherent way to the exogenous forcing. A possible exception is the CCSM model, which exhibits increased aridity during the eleventh and twelfth centuries (particularly manifest in the number of SW+CP+SE drought occurrences, although the number of droughts in both centuries is within the range from the CCSM control simulation); changes that are contemporaneous with those in the NADA. The models do, however, appear to simulate a large range in the number of pancontinental droughts. To test if this range is of the magnitude observed in the NADA, the number of each pancontinental drought type was calculated for a sliding 100-yr window across the forced model simulations and NADA record and the range is plotted in Fig. 11. Each model is individually capable of simulating centennial-scale variability in the frequency of pancontinental drought occurrence that is characteristic of the NADA. The fact that models simulate large differences in the number of pancontinental drought features for different 100-yr periods, and that these changes are not tied in any coherent way to the exogenous forcing, is suggestive of a large amount of internal variability on centennial time scales. This model behavior is consistent with previous work that has focused on the dynamics of persistent drought in the southwest specifically (Coats et al. 2013, 2015). Additionally, if the model dynamics are in fact representative of the real atmosphere–ocean system, then this result indicates that the observed preponderance of pancontinental droughts in the medieval period could have arisen from internal variability, as opposed to changes in radiative forcing.
d. Simulated pancontinental drought dynamics
Figure 12 shows the composite of the Niño-3.4, PDO, and AMO indices during each type of pancontinental drought for the model control simulations and for the overlapping period between the NADA and the observed SST dataset (1854–2005 CE), with this analysis completed separately for the three- and four-region droughts. Significance at the 90% level using the bootstrapping test described in section 2c is indicated with an asterisk. Importantly, the models do not agree with each other or the observations on the dynamics that drive the different types of pancontinental drought, despite the fact that they produce patterns and recurrence intervals that are consistent with the NADA (Figs. 8–10). The one exception appears to be the importance of ENSO to the SW+CP+SE, SW+CP+NW, SW+NW+SE, and SW+CP+NW+SE pancontinental drought types, which is exhibited to varying degrees by all of the models and for the observation-to-reconstruction. Additionally, with the exception of CCSM and MPI, the models and the observations individually have pancontinental drought types with no significant connection to the major modes of atmosphere–ocean variability. Together these results imply what was suggested by the large CPCS range in the bottom panel of Fig. 9 and by the characteristic model behavior outlined in section 3a, namely that multiple dynamical drivers are capable of producing each type of pancontinental drought.
The left-hand panel in Fig. 13 shows the observed Niño-3.4, PDO, and AMO indices with the timing of each pancontinental drought occurrence marked with a gray bar. While the ENSO behavior appears to be consistent throughout the instrumental interval, the variability in the early part of the PDO and AMO records is subdued. This behavior may explain the weaker PDO and AMO teleconnections in the reconstructed PDSI, relative to observations (Fig. 1, bottom) as the reconstruction spans the length of the 152-yr instrumental SST record. The observed PDSI, by contrast, is analyzed only over the 1950–2005 period, for which both modes exhibit consistently large variability.
The right-hand panels in Figs. 13 and 14 analyze the associations between pancontinental droughts and the dynamical drivers using the observation-to-reconstruction and the models, respectively, but allowing for a more robust statistical assessment of the possibility that it is a combination of dynamical modes that produces pancontinental drought features. To do so, all pancontinental drought types are treated as the same, and considered to be Bernoulli processes (with 1 for drought years and 0 for nondrought years), with the drought frequency then defined as the number of pancontinental drought occurrences over the number of analyzed years (152 yr for the observations–reconstruction and 500 yr for the model control simulations). Within a Bayesian framework, the posterior distribution of the drought frequency can be calculated for subsets of the data that have different phases of the dynamic modes (e.g., a negative or La Niña–like state in the Niño-3.4 index) or some combination of phases of the dynamic modes (e.g., a negative Niño-3.4 index, positive AMO, and negative PDO). If, following Kam et al. (2014), we assume that the prior distribution is a uniform beta distribution or uninformative, then the posterior distribution is easily derived with the alpha and beta parameters being equal to the number of drought occurrences plus one and the number of years minus the number of drought occurrences plus one, respectively. The observation-to-reconstruction posterior distributions in the right-hand panel of Fig. 13 indicate that the frequency of pancontinental drought occurrence is greatest when there are simultaneously negative Niño-3.4 and PDO indices and a positive AMO index (with pancontinental drought occurring nearly 40% of the time when these conditions are met). Interestingly, for the observation-to-reconstruction record the individual impacts of the three modes of variability on the frequency of pancontinental drought are roughly equal. The models, by contrast, tend to overestimate the impact of ENSO on pancontinental drought occurrence (with the main exception being BCC; see Fig. 14). This result is likely indicative of the more realistic, and in some cases overactive, ENSO variability and teleconnections (e.g., CCSM) relative to other modes of coupled atmosphere–ocean variability.
The model split between slightly underestimating the recurrence interval of pancontinental drought of all types (CCSM, IPSL, and MPI) and moderately overestimating the recurrence interval of these droughts (BCC, GISS, and MIROC) in Fig. 8 can also be explained by the results in Figs. 13 and 14. CCSM, IPSL, and MPI all overestimate the impact of ENSO on the frequency of pancontinental drought occurrence and, as a consequence, produce more of these features than is realistic (Fig. 14). In CCSM this behavior appears to result from ENSO variability that is too strong (Figs. 1 and 2), while in IPSL and MPI the ENSO variability is more realistic (though slightly too strong) but the hydroclimate response to ENSO is too homogenous over the North American continent (Fig. 1). BCC, despite underestimating the impact of ENSO on pancontinental drought occurrence, is able to largely reproduce the impact of the AMO, while slightly overestimating the PDO impact (Fig. 1), and has a posterior distribution of ocean-forced pancontinental drought occurrence that is similar to the observation-to-reconstruction (Fig. 14). The same is true of GISS, which exhibits the most realistic impact (as compared to the observations) of the three modes of atmosphere–ocean variability (and the oceanic boundary conditions in general; see Fig. 14) on pancontinental drought occurrence (although it slightly underestimates the AMO impact relative to the PDO). The fact that GISS, and to a lesser degree BCC, underestimate the overall occurrence of pancontinental drought (e.g., Fig. 8), therefore, appears to be related to the frequency with which the simulated ocean produces a simultaneously positive AMO, negative PDO, and negative ENSO. MIROC also underestimates the frequency of pancontinental drought occurrence, and while it simulates a reasonable impact of a negative Niño-3.4 index on the frequency of pancontinental drought occurrence, the PDO and AMO impacts are too weak (Fig. 14).
It must be noted that between the short instrumental record, which limits our knowledge of real-world pancontinental drought dynamics, and the general lack of consistency in the significance of connections between the dynamical modes and pancontinental drought events in both the models and the NADA, it is difficult to make conclusions about the veracity of the simulated dynamics. Nevertheless, attempting to understand the simulated pancontinental drought dynamics within the context of the characteristic model behavior outlined in section 3a may help determine if models will be able to properly constrain the risk of future drought over North America. In particular, the patterns of behavior of CCSM and BCC are an interesting juxtaposition of pancontinental drought dynamics. CCSM has significant connections between four of the pancontinental drought types and the tropical Pacific (Fig. 12) and, in general, greatly overestimates the impact of ENSO on pancontinental drought occurrence (Fig. 14). BCC, on the other hand, has no significant connections between pancontinental drought and the tropical Pacific (Fig. 12) and, generally, underestimates the impact of ENSO on pancontinental drought occurrence (Fig. 14). This behavior can perhaps be understood in terms of the model dynamics outlined in section 3a. While BCC has a somewhat realistic ENSO spatial pattern (Fig. 5), and variability that is moderate (but too regular, e.g., the large negative lag-1 autocorrelation value; see Fig. 2), the ENSO teleconnection to North America is the least realistic, weakest, and most nonstationary of the models analyzed herein (Fig. 1). By contrast, the PDO and AMO teleconnections are more realistic, stronger (Fig. 1), and consequently more strongly connected to pancontinental drought (Fig. 14). CCSM, on the other hand, has a very realistic and stationary ENSO teleconnection (Fig. 1), along with ENSO variability and spatial patterns that are too strong (Figs. 5 and 6). As a consequence, nearly all of the pancontinental drought types exhibit a connection to the tropical Pacific. Interestingly, the PDO and AMO in CCSM are less realistic and exhibit relatively large and small variability (Figs. 5 and 6), respectively; yet both have a significant connection to three of the pancontinental drought types (Fig. 12).
The behavior of the other models is less clear but will be considered separately for models with weak ENSO (GISS and MIROC) and strong ENSO (IPSL and MPI) connections to pancontinental drought in Fig. 14. For weak ENSO, GISS, while only exhibiting a significant connection between ENSO and two of the pancontinental drought types (Fig. 12), largely captures the oceanic impact on the frequency of pancontinental drought occurrence (Fig. 14). This is surprising given that the GISS model underestimates the overall frequency of pancontinental drought occurrence and has generally weak ocean variability (e.g., Coats et al. 2015). MIROC also has a weak connection between the tropical Pacific and pancontinental droughts (with just the connection between SW+CP+SE and ENSO being significant at the 90% level; see Fig. 12), and this is likely related to the highly nonstationary and unrealistic ENSO teleconnection within the model (Fig. 1).
MPI, like CCSM, has a strong connection between ENSO and pancontinental droughts with all of the drought types exhibiting significance (Fig. 12). While the ENSO teleconnection in MPI is less realistic than in CCSM (Fig. 1), the model also exhibits relatively large-magnitude tropical Pacific variability (Fig. 6). IPSL, like CCSM and MPI, overestimates the impact of ENSO on pancontinental drought occurrence (Fig. 14); it also exhibits a moderate connection between the individual pancontinental drought types and ENSO, with the SW+CP+SE and SW+CP+NW+SE time series being significant at the 90% level (Fig. 12). The ENSO teleconnection in IPSL is likewise moderately realistic, as well as highly stationary and driven by moderate variability in the tropical Pacific (Figs. 1 and 6).
Simulated ENSO, PDO, and AMO dynamics, and their teleconnections to North America, differ between models and in their comparisons to observations. As a consequence, models do not agree on the modes of atmosphere–ocean variability that are associated with pancontinental droughts. The models do, however, simulate pancontinental droughts with the frequency and spatial patterns exhibited by the NADA. Additionally, the models display centennial-scale variability in the occurrence of pancontinental droughts that is similar to the magnitude observed in the NADA. These changes do not appear to be tied to the exogenous forcing, suggesting that simulated internal hydroclimate variability on these time scales is large in magnitude.
These results have implications for efforts to project future hydroclimate over North America and to understand past hydroclimate change. First, if centennial-scale internal variability is large in magnitude, the MCA-to-LIA transition and other periods of past hydroclimate change may have resulted from internal variability of the atmosphere–ocean system. If true, this would suggest that the projection of certain characteristics of future hydroclimate is potentially complicated in terms of the impact of trace gas concentrations. Nevertheless, internal variability, if manifest in changes to the oceanic boundary conditions that in turn drive hydroclimate change, may be predictable if interannual-to-multidecadal modes of Pacific and Atlantic Ocean variability can be predicted. While operational prediction of the tropical Pacific, for instance, has yet to exceed the seasonal time scale (e.g., Barnston et al. 2012; Zhu et al. 2012), the slowly varying nature of the PDO and AMO suggests that it is a potential target for predictability on decadal time scales (e.g., Smith et al. 2012). Within the context of the results presented herein, PDO and AMO predictability would help constrain risk assessments of future pancontinental drought occurrence, given that a negative PDO and a positive AMO are associated with an increased frequency of pancontinental drought occurrence (Fig. 13). In particular, decadal AMO and PDO predictability along with seasonal ENSO prediction has the potential to provide robust annual projections of the risk of pancontinental drought occurrence.
Second, models reproduce the statistics of pancontinental drought occurrence, but they do not agree on the dynamics that drive these features. The ENSO (e.g., Clement et al. 1996) and AMO (e.g., Ting et al. 2009) both potentially respond to radiative forcing. If these responses are large in magnitude, the simulated hydroclimate change over North America will be different for each model. In CCSM, for instance, a large forced response in the tropical Pacific would be expected to drive large hydroclimate impacts over North America. The same change in BCC, however, would yield smaller impacts. It is difficult to determine which of these model responses is realistic, and by consequence, which model projection should be considered the most accurate, because even in models that are successful at simulating the observed atmosphere–ocean dynamics, the dynamical relationships are often nonstationary and their connections to radiative forcings are not well constrained. This issue is compounded by the fact that the observed dynamics themselves have been inferred from a 152-yr instrumental interval that cannot provide a full assessment of the stationarity of real-world dynamics. Understanding whether teleconnection nonstationarity is physically reasonable, and whether the instrumental interval is sufficient to characterize the full range of real-world atmosphere–ocean dynamics, is thus essential for constraining the risk of hydroclimate change in North America. This will require, at a minimum, longer records of proxy-estimated SSTs [e.g., Emile-Geay et al. (2013) for the tropical Pacific] and other regional records of hydroclimate variability over the last millennium [e.g., Cook et al. (2010) for monsoonal Asia].
A characterization of the full range of real-world atmosphere–ocean dynamics, and their connections to pancontinental droughts, will help determine which models exhibit realistic dynamics and thus can be interpreted as projecting the future hydroclimate of North America with some accuracy. We have demonstrated that models simulate pancontinental droughts that are characteristic of a paleoclimate estimate (in spatial character and frequency of occurrence), but driven by different atmosphere–ocean dynamics, and that these models simulate a large degree of variability in the occurrence of pancontinental droughts on centennial time scales. These results, and further efforts to characterize relevant model dynamics, will help clarify the interpretation of future hydroclimate projections by providing, in particular, an understanding of the reasons for the differences between model projections of hydroclimate responses to increased greenhouse gas forcing. This understanding, in turn, will help determine what information from the future projections is useful for planning adaptation and management strategies for the impacts of climate change.
This work was supported by NOAA Awards NA10OAR4310137 (Global Decadal Hydroclimate Variability and Change, GloDecH) and NA11OAR4310166, as well as NSF Awards AGS-1243204 and AGS-1401400. We acknowledge the World Climate Research Programme’s Working Group on Coupled Modelling, which is responsible for CMIP, and we thank the climate modeling groups (listed in Table 1) for producing and making available their model output. For CMIP, the Department of Energy’s Program for Climate Model Diagnosis and Intercomparison provides coordinating support and led the development of the software infrastructure in partnership with the Global Organization for Earth System Science Portal. We also thank Haibo Liu and Naomi Henderson for their considerable computational and data management support.
Lamont-Doherty Earth Observatory Contribution Number 7861.