Projected changes in global rainfall patterns will likely alter water supplies and ecosystems in semiarid regions during the coming century. Instrumental and paleoclimate data indicate that natural hydroclimate fluctuations tend to be more energetic at low (multidecadal to multicentury) than at high (interannual) frequencies. State-of-the-art global climate models do not capture this characteristic of hydroclimate variability, suggesting that the models underestimate the risk of future persistent droughts. Methods are developed here for assessing the risk of such events in the coming century using climate model projections as well as observational (paleoclimate) information. Where instrumental and paleoclimate data are reliable, these methods may provide a more complete view of prolonged drought risk. In the U.S. Southwest, for instance, state-of-the-art climate model projections suggest the risk of a decade-scale megadrought in the coming century is less than 50%; the analysis herein suggests that the risk is at least 80%, and may be higher than 90% in certain areas. The likelihood of longer-lived events (>35 yr) is between 20% and 50%, and the risk of an unprecedented 50-yr megadrought is nonnegligible under the most severe warming scenario (5%–10%). These findings are important to consider as adaptation and mitigation strategies are developed to cope with regional impacts of climate change, where population growth is high and multidecadal megadrought—worse than anything seen during the last 2000 years—would pose unprecedented challenges to water resources in the region.

Implied by Fig. 1 is the possibility that deterministic simulations of climate change using state-of-the-art numerical models may be insufficient for estimating megadrought risk because the ensemble sizes of such experiments are relatively small (tens of realizations per model at most), and the statistics of infrequent events such as megadroughts might not be robust. Using a multimodel ensemble does not completely guard against this limitation because model simulations disagree on the expression of forced changes in hydroclimate at regional scales (e.g., Diffenbaugh and Giorgi 2012 ). Instead, we use large ensembles of stochastic variables to emulate the statistics of interannual to decadal variability, and output from global climate models to estimate how precipitation is expected to change this century. The limitations and possible implications of this assumption are discussed in section 4 .

Schematic illustration of why large ensembles are needed to calculate megadrought risk. (a) The black line shows the original and shifted mean, while (b),(c) the black lines show the original mean for reference. Importantly, the means and variances are the same for the final 50 yr in (a)–(c) but only the realization in (b) experiences a megadrought.

Schematic illustration of why large ensembles are needed to calculate megadrought risk. (a) The black line shows the original and shifted mean, while (b),(c) the black lines show the original mean for reference. Importantly, the means and variances are the same for the final 50 yr in (a)–(c) but only the realization in (b) experiences a megadrought.

The scenario delineated above is shown schematically in Fig. 1 . Here, an idealized time series of some hydrological variable (say P − E) has been generated with unit variance and a mean of zero for the first 100 “years” ( Fig. 1a ). At year 101 the mean is shifted by −0.25σ and an additional 50 years of data are generated while the variance stays the same. Figures 1b and 1c show realizations of 50 yr of data with the same mean and variance as the final 50 yr of the series in Fig. 1a . Although both the time series in Figs. 1b and 1c have the same mean and variance, a prolonged period of time with low values (a “megadrought”) is found in the first realization ( Fig. 1b ), whereas in the second realization ( Fig. 1c ) it is not.

Assumption 3 in the list above deserves further elaboration. We begin by noting that in western North America over the last millennium, stochastic variability and autocorrelation alone may account for the magnitude of hydroclimatic variations on time scales from years to decades ( Hunt 2011 ; Ault et al. 2013 ; Coats et al. 2013b ). Second, one can easily imagine a situation where in a single realization of a given model, climatic forcing enhances overall aridity, but megadroughts do not occur because a few intermittent wet years disrupt their duration. Given the statistics of this model, megadroughts might still be likely, but would not be found in this particular realization.

Justifications for statements 1 and 2 are straightforward: state-of-the-art models agree that semiarid subtropical regions throughout the world will tend to dry under climate change (e.g., Diffenbaugh and Giorgi 2012 ), and paleoclimate records, especially tree rings, are reasonably well validated and widely used to characterize variations of the past for a wide range of water resource management applications ( Meko et al. 2012 ).

Simple models of time series are adequate for simulating the local statistical characteristics of hydroclimate across interannual to multidecadal time scales, regardless of whether these characteristics are externally forced or internally generated.

This paper estimates future prolonged drought risk using information from instrumental records, paleoclimate archives, and climate model simulations in simple Monte Carlo realizations of hydroclimate. The motivation for doing so comes from our notion of risk as a fractional quantity referring to the likelihood of prolonged drought occurrence. We rely on global climate model simulations of change during the twenty-first century as estimates of mean conditions in the future, and we use simple statistical models to build up large ensembles for calculating risk. This technique assumes the following:

Map of projected precipitation over the twenty-first-century (2005–2100) change in the RCP8.5 scenario shown as a percentage of twentieth-century precipitation change [as in the global maps of Diffenbaugh and Giorgi (2012) ]. For each model, the number of available runs from each experiment is shown in parentheses in the following order: historical, piControl, RCP2.6, RCP4.5, and RCP8.5. The red box shows the greater southwestern United States to emphasize the focus of this study (30°–40°N, 120°–103°W).

Map of projected precipitation over the twenty-first-century (2005–2100) change in the RCP8.5 scenario shown as a percentage of twentieth-century precipitation change [as in the global maps of Diffenbaugh and Giorgi (2012) ]. For each model, the number of available runs from each experiment is shown in parentheses in the following order: historical, piControl, RCP2.6, RCP4.5, and RCP8.5. The red box shows the greater southwestern United States to emphasize the focus of this study (30°–40°N, 120°–103°W).

The 11-yr running means of normalized paleoclimate reconstructions for (a) the Colorado streamflow ( Meko et al. 2007 ), (b) reconstructed PDSI from the southwestern United States ( Cook et al. 2004 ), and (c) reconstructed PDSI from Mexico ( Stahle et al. 2011 ). Vertical gray bars indicate decadal-scale drought (a −0.5σ deviation in the 11-yr mean). All time series are standardized to exhibit one unit of standard deviation and a mean of zero over the 1950–2000 reference period.

The 11-yr running means of normalized paleoclimate reconstructions for (a) the Colorado streamflow ( Meko et al. 2007 ), (b) reconstructed PDSI from the southwestern United States ( Cook et al. 2004 ), and (c) reconstructed PDSI from Mexico ( Stahle et al. 2011 ). Vertical gray bars indicate decadal-scale drought (a −0.5σ deviation in the 11-yr mean). All time series are standardized to exhibit one unit of standard deviation and a mean of zero over the 1950–2000 reference period.

The 11-yr running means of normalized paleoclimate reconstructions for twentieth-century precipitation data from (a) the U.S. Southwest and (b) the Great Plains. Precipitation data are from the University of East Anglia’s Climate Research Unit’s time series version 3.1 (TS3.1) dataset ( Mitchell and Jones 2005 ). We identify decadal droughts as −0.5 standard deviations of raw data in the 11-yr mean (vertical gray bars).

The 11-yr running means of normalized paleoclimate reconstructions for twentieth-century precipitation data from (a) the U.S. Southwest and (b) the Great Plains. Precipitation data are from the University of East Anglia’s Climate Research Unit’s time series version 3.1 (TS3.1) dataset ( Mitchell and Jones 2005 ). We identify decadal droughts as −0.5 standard deviations of raw data in the 11-yr mean (vertical gray bars).

Our definition of decadal drought captures major intervals of aridity during the twentieth century as well as others during the last millennium ( Figs. 3 and 2 ). We employ a second and more stringent criterion to identify multidecadal megadrought. In this case, −0.5σ departures in the 35-yr mean are identified. Although this definition is somewhat arbitrary, it is useful because the thresholds employed are both longer in time and greater in magnitude than the descriptions of Meko et al. (2007) and Cook et al. (2007) used to characterize the worst droughts of the past millennium in Colorado streamflow and continental-scale hydroclimate, respectively. By setting the criterion for multidecadal megadrought beyond anything experienced during the last millennium, we suggest that our results will be insightful for developing adaptation and mitigation strategies for addressing worst-case scenarios. We also note that a 35-yr, −0.5σ event would be on par with the consequential late twentieth-century Sahel drought (e.g., Hoerling et al. 2006 ). Other similar thresholds (e.g., 25 or 45 yr) would produce qualitatively similar spatial patterns, but the rates of megadrought occurrence and levels of megadrought risk would of course differ from the ones we report here. In particular, for a given magnitude (say −0.5σ) the risk of shorter events would be higher and the risk of longer events would be lower than that for a 35-yr event.

where F(t) is reconstructed flow and and are the mean and standard deviation, respectively, of the annual data over the reference period of 1950–2000 CE. The time series of is a modified z score of F(t), and its values through time are shown in Fig. 3a . Identifying intervals of −0.5σ departures in the running 11-yr mean highlights the 1150s, as well as several other low-flow decades, which occur about once per century (gray vertical bars). Time series from other recent drought studies ( Cook et al. 2004 ; Stahle et al. 2011 ), normalized in the same way, are also shown in Figs. 3b and 3c . They suggest that the preindustrial rate of comparable decade-long droughts is ~1.5% century −1 , which is quite consistent with the literature-based estimate of 1%–2% century −1 of Woodhouse and Overpeck (1998) .

Identifying −0.5σ events in the 11-yr means of paleoclimate records requires us to normalize these time series to exhibit unit variance over the twentieth century, so that fluctuations in the past are scaled relative to this baseline period. To that end we represent the entire Colorado streamflow record as normalized departures [ ] from the late twentieth-century mean:

Here we develop a systematic approach to normalizing hydroclimate fluctuations so that they retain their essential meaning whether they originate from climate model simulations, observational datasets, or paleoclimate reconstructions. We further seek to distinguish between decadal droughts, which have been experienced during the instrumental era (e.g., the 1930s Dust Bowl), and multidecadal megadrought events that are outside the range of variability experienced during the twentieth century. To begin, we consider two of the worst decade-scale droughts during the twentieth century: the 1930s Dust Bowl and the 1950s Southwest drought. Both of these intervals can be identified as −0.5σ departures in the decadal (11-yr) running mean of precipitation ( Fig. 2 ).

With the definitions of “decadal drought” (an 11-yr, −0.5σ event) and “multidecadal megadrought” (a 35-yr, −0.5σ event) that we have outlined above, we now develop “null” expectations for the rate at which these events would occur from random chance under a stationary climate, but with three different assumptions about the underlying frequency characteristics of hydroclimate variability on interannual to centennial time scales.

We begin by examining the statistics of prolonged drought when interannual hydroclimate fluctuations are simulated as normally distributed white noise with unit variance and standard deviation. An example of one such time series, X w (t) is shown in Fig. 5. The decadal drought statistics of this type of noise, obtained from 1000 white noise realizations (each of length 100 years), are summarized in Fig. 6. If the distribution of variance across the hydroclimatic continuum were indeed white, then decadal droughts would be expected to occur at a rate of slightly <1 (100 yr)−1 (Fig. 6a), and the risk of such an event occurring during any given 50-yr period would be around 45% (Fig. 6b). The likelihood of a multidecadal megadrought during any given 50-yr period would be only about 0.45% (Fig. 6b). This preliminary Monte Carlo result establishes a benchmark for the minimum rate at which decadal droughts and megadroughts would occur in a climate with only stochastic interannual variability and no sources of long-term persistence.

Fig . 5. View largeDownload slide Examples of Monte Carlo (MC) time series used to simulate decadal drought with three different underlying frequency characteristics [white noise, AR(1), and power law]: (a) 50 yr of a single MC realization of each type of noise, (b) the 1000-yr full realization of each noise type, and (c) the 11-yr running mean of each type of noise, shown with the dashed line denoting a decadal drought (e.g., a −0.5σ in the 11-yr running mean). Importantly, each MC time series has a mean of zero and unit standard deviation, and differ only in the distribution of their variances across the power spectrum. Moreover, the AR(1) and power-law time series are generated by rescaling the white noise data, which is why the different realizations of noise appear so strongly correlated with each other. Fig . 5. View largeDownload slide Examples of Monte Carlo (MC) time series used to simulate decadal drought with three different underlying frequency characteristics [white noise, AR(1), and power law]: (a) 50 yr of a single MC realization of each type of noise, (b) the 1000-yr full realization of each noise type, and (c) the 11-yr running mean of each type of noise, shown with the dashed line denoting a decadal drought (e.g., a −0.5σ in the 11-yr running mean). Importantly, each MC time series has a mean of zero and unit standard deviation, and differ only in the distribution of their variances across the power spectrum. Moreover, the AR(1) and power-law time series are generated by rescaling the white noise data, which is why the different realizations of noise appear so strongly correlated with each other.

Fig . 6. View largeDownload slide Statistics summarizing Monte Carlo simulations: (a) distributions of years spent in decadal (≥11 yr) drought conditions for each type of noise (as a percentage of all realizations), (b) risk of at least one decadal (≥11 yr) drought during any given 50-yr window in any realizations, and (c) risk of a multidecadal (≥35 yr) drought during any 50-yr window. Risk in (b) and (c) is expressed as a percentage of the total number of simulations. Fig . 6. View largeDownload slide Statistics summarizing Monte Carlo simulations: (a) distributions of years spent in decadal (≥11 yr) drought conditions for each type of noise (as a percentage of all realizations), (b) risk of at least one decadal (≥11 yr) drought during any given 50-yr window in any realizations, and (c) risk of a multidecadal (≥35 yr) drought during any 50-yr window. Risk in (b) and (c) is expressed as a percentage of the total number of simulations.

In a simplistic sense, year-to-year persistence can be described as a first-order autoregressive [AR(1)] process [X AR(1) (t)]:

where α is the lag-1 (i.e., 1 yr) autocorrelation coefficient and is derived empirically from data, X AR(1) (t) is autoregressive red noise, and X w (t) is the white noise input. In WNA, the value of α is about 0.3 on interannual time scales for the three (tree ring based) paleoclimate reconstructions shown in Fig. 3, as well as for many other hydroclimate indicators (Ault et al. 2013). A single realization of this type of noise (normalized to exhibit unit variance overall) is shown in Fig. 5, and the statistical characteristics of megadroughts in this type of noise are shown in Fig. 6.

Despite the intuitive and simple representation of hydroclimate as an AR(1) process—moisture deficits tend to persist through time—there is some evidence that such an approximation misses key characteristics of variability on longer time scales (Pelletier and Turcotte 1997; Kantelhardt et al. 2006; Koscielny-Bunde et al. 2006; Ault et al. 2013). As a complementary approach, we also simulate hydroclimate as a process with underlying frequency characteristics that are described by a weak power-law relationship between frequency f and variance S(f), such that S(f) ∝ f−β. Power spectra with higher values of β correspond to time series that exhibit more variance at lower frequencies. To generate time series with this type of frequency behavior, we employ a method similar to the one described by Pelletier and Turcotte (1997) and explained thoroughly in Pelletier (2008). First, we calculate the discrete Fourier transform of a white noise time series X w (t), and filter it to conform to a predefined value of β:

where k are the standard Fourier frequencies and N is the length of the time series. The term ψ k rescales the Fourier coefficients so that they are approximately power-law distributed:

Here the value of β is divided by because it is being applied to the raw Fourier coefficients, which have amplitudes proportional to the square root of the power spectrum.

The rescaled Fourier series is then used to generate power-law time series X p (t) by taking the real part of the inverse Fourier transform of ⁠:

Finally, the mean and variance are restored to the values of the original white noise data (zero and unity, in this case).

We used a value of 0.5 for β to rescale each realization of X w (t), which was suggested as an appropriate estimate by Ault et al. (2013) from synthesis of tree-ring reconstructions of precipitation, PDSI, and streamflow as well as non-tree-ring estimates of hydroclimate. As a check, we calculated the power laws of the noises after they had been rescaled. We found that the actual values of β varied from one realization to the next, but were generally between 0.4 and 0.6. This range agrees well with instrumental and paleoclimate estimates of this parameter for the region, and is certainly within the observational uncertainty (Ault et al. 2013). Importantly, time series with spectra scaled by power laws of ~0.5 will also appear to exhibit autocorrelation of about 0.3, which in turn implies that the AR(1) and power-law realizations will behave very similarly on short time scales, but not necessarily on longer ones (e.g., Pelletier and Turcotte 1997; Ault et al. 2013). Finally, our use of power-law noises does not make any assumptions about the underlying climate dynamics governing the shape of the power spectrum of hydroclimate: linear and nonlinear processes alike may produce such spectral distributions (Milotti 1995; Penland and Sardeshmukh 2012).

Table 1 highlights a few key features of the two models employed here. In particular, the noise models used to estimate drought risk use parameters that do not vary across space, and all are scaled to the twentieth-century mean and variance. The autocorrelation parameter of 0.3 is a middle-of-the-road value from the time series shown in Fig. 3, and is well within the range of estimates for autocorrelation in the region from other paleoclimate and observational datasets (Ault et al. 2013). The value used for β (0.5) is from the analysis of proxy records in Ault et al. (2013) as well, and is supported by Fig. 7.

Fig . 7. View largeDownload slide Power-law estimates (β) from (a) twentieth-century instrumental data and (b) and CMIP5 historical (1850–2005 CE) simulations. Instrumental data originate from the University of East Anglia’s TS3.1 data product and, like the CMIP5 records, were annualized prior to calculating β. Fig . 7. View largeDownload slide Power-law estimates (β) from (a) twentieth-century instrumental data and (b) and CMIP5 historical (1850–2005 CE) simulations. Instrumental data originate from the University of East Anglia’s TS3.1 data product and, like the CMIP5 records, were annualized prior to calculating β.