## Abstract

The distribution of climatic variance across the frequency spectrum has substantial importance for anticipating how climate will evolve in the future. Here power spectra and power laws (*β*) are estimated from instrumental, proxy, and climate model data to characterize the hydroclimate continuum in western North America (WNA). The significance of the estimates of spectral densities and *β* are tested against the null hypothesis that they reflect solely the effects of local (nonclimate) sources of autocorrelation at the monthly time scale. Although tree-ring-based hydroclimate reconstructions are generally consistent with this null hypothesis, values of *β* calculated from long moisture-sensitive chronologies (as opposed to reconstructions) and other types of hydroclimate proxies exceed null expectations. Therefore it may be argued that there is more low-frequency variability in hydroclimate than monthly autocorrelation alone can generate. Coupled model results archived as part of phase 5 of the Coupled Model Intercomparison Project (CMIP5) are consistent with the null hypothesis and appear unable to generate variance in hydroclimate commensurate with paleoclimate records. Consequently, at decadal-to-multidecadal time scales there is more variability in instrumental and proxy data than in the models, suggesting that the risk of prolonged droughts under climate change may be underestimated by CMIP5 simulations of the future.

## 1. Introduction

In western North America (WNA), the term “megadrought” typically refers to an interval of aridity in the past that is more prolonged than anything seen during the twentieth century (Woodhouse and Overpeck 1998; Cook et al. 1999, 2004; Stahle et al. 2007; Meko et al. 2007; Cook et al. 2010). If such an event were to occur in the future, it would likely put unprecedented strains on existing water resources. Understanding megadrought risk is therefore critical to successfully managing water in the coming decades, which in turn motivates us to characterize the hydroclimatic continuum in WNA on the interannual and longer time scales on which megadroughts occur.

While observational and paleoclimate records reveal underlying statistical qualities of the hydroclimate continuum, climate models are needed to assess how rising greenhouse gases (GHGs) and natural variability will govern risk in the future. It is critical, however, to evaluate the hydroclimate continuum in these climate models to build confidence in their projections of the future. If climate models reproduce a realistic continuum when properly forced and integrated for many centuries, then it is likely that they also capture the full range of future risks originating from external and internal sources of variability.

The paleoclimate record reveals that megadroughts are common in WNA over the last millennium: Woodhouse and Overpeck (1998) showed that droughts like those of the 1950s and 1930s occurred once or twice per century since about 1000 common era (CE) and that far more severe and more persistent dry intervals occurred as recently as the sixteenth century. Tree-ring-based reconstructions of flow on the Colorado River over the last 1200 years (Meko et al. 2007) show about one megadrought per century and at least one severe event that reduced mean flow by 15%. In addition, Woodhouse et al. (2006) showed that the early part of the twentieth century—the climate window upon which Colorado River allocations are based—was one of the wettest in the last 500 years. Reconstructions of the Palmer drought severity index (PDSI; Cook et al. 2004) and a collection of other proxy records (Graham et al. 2007) provide further evidence for a prolonged period of widespread drought between approximately 800 and 1200 years ago, as well as multiple decadal scale intervals with variability matching or exceeding the range of twentieth-century PDSI variability (Cook et al. 2004). Most recently, new research (Routson et al. 2011) has documented another earlier period of exceptional drought, including what was probably the most severe drought in the Colorado headwaters of the last 2000 years.

Several factors complicate outright comparisons of paleoclimate records and model simulations. First, the time evolution of internally driven hydroclimate fluctuations in model simulations will differ from the actual paleoclimate history, even if both types of data exhibit similar underlying behavior. Second, model biases may affect the definition of hydrological extremes: if the model precipitation is greater than observed in a given region, then simulated “droughts” could be wetter than those in nature. Third, proxy records are sensors of climate and filter information in ways that can distort underlying signals. Model output, in contrast, experiences no such filtering.

We minimize the limitations above by focusing on the power spectra of moisture balance indices derived from observational data, paleoclimate records, and climate model output. In particular, we estimate the power-law-scaling parameter (*β*) that relates hydroclimate variance to time scale such that *S*(*f*) ∝ *f*^{−β} where *f* is frequency and *S*(*f*) is the power spectrum (e.g., Pelletier and Turcotte 1997; Huybers and Curry 2006). This approach allows us to ignore the time evolution of hydroclimate and data model differences in the mean and variance of precipitation. If a model simulates time series with spectra that agree well with observations, this would imply that it is also simulating realistic relative amplitudes of variations across a wide range of time scales. The model, in turn, would therefore faithfully capture the relative magnitude and risk of megadroughts. As it is likely that a power-law distribution is not the *best* fit to the power spectra of all records considered here (e.g., Clauset et al. 2009; Vyushin et al. 2009), for clarity we refer to the values of *β* we estimate as “spectral slopes.” These slopes give us insight into the underlying relationship between time scale and variance across disparate types of data.

Values of *β* near zero describe “white” spectra, which have an approximately uniform distribution of variance across frequencies and no significant autocorrelation. “Red” spectra, and “redness” in the power spectrum, are associated with time series that exhibit more variance at longer time scales. Negative values of *β*, on the other hand, describe “blue” spectra. Such spectra can arise in systems with negative autocorrelation, because from one time step to the next the system is pushed away from its previous value.

Power spectra of temperature fluctuations in the ocean and atmosphere are characterized by spectral slopes (*β*) between 0 and 0.5 on land and near 1.0 for most of the ocean (Fraedrich and Blender 2003; Huybers and Curry 2006). More recently, Vyushin and Kushner (2009) and Vyushin et al. (2009) have used global climate model and reanalysis data to link positive spectral slopes in the troposphere to tropical sea surface temperature (SST) variability and similar behavior in the stratosphere to volcanic activity. Hydroclimate variables have also been shown to exhibit similar values, with *β* estimated near 0.5 for globally averaged spectra of river discharge, tree-ring widths, and precipitation records (Pelletier and Turcotte 1997).

Our study builds on previous attempts to characterize the continuum of climate variability in a number of critical regards. First, whereas previous studies have aimed at estimating spectral slopes of temperature at large spatial scales (e.g., Fraedrich and Blender 2003; Huybers and Curry 2006; Vyushin and Kushner 2009; Vyushin et al. 2009), we are focused on characterizing hydroclimate fluctuations in WNA. We therefore restrict our analysis to reconstructions of hydroclimate variables and paleoclimate records that are indicative of surface moisture balance. Although Pelletier and Turcotte (1997) used long tree-ring chronologies, their approach did not screen these records for moisture sensitivity. Further, in our study, hydroclimate spectra are evaluated against the null hypothesis that they arise solely as a consequence of local autocorrelation (e.g., because of storage in soil or surface water). Finally, we analyze output from state-of-the-art coupled general circulation models (GCMs) that have been forced with the time-evolving boundary conditions of the last millennium. These GCMs are part of phase 5 of the Coupled Model Intercomparison Project (CMIP5) and phase 3 of the Paleoclimate Modeling Intercomparison Project (PMIP3) archives. As such, they are the most scientifically advanced tools for numerically simulating climate variability and change, and they are being used throughout the world to anticipate and prepare for climate change.

The observational, paleoclimatic, and GCM datasets employed in this analysis are described in section 2. In section 2c, we develop analytical and statistical expectations for the power spectrum of hydroclimate variability by considering the hypothetical spectrum of an autoregressive [AR(1)] process forced by white noise and with realistic autocorrelation occurring on the monthly time scale. Our results (section 3) indicate that hydroclimate fluctuations in most tree-ring reconstructions are indeed consistent with the null hypothesis. However, this perspective is incomplete because tree-ring-based reconstructions are usually preprocessed using methods that remove low-frequency nonclimate (and potentially climate) variability. Evidence from a wider range of hydroclimate proxies suggests that variability on decadal and longer time scales is well above what a simple AR(1) process would produce. State-of-the-art GCMs are likewise unable to generate strong variations on these time scales, and the implications of this finding for assessing drought risk are discussed in section 5.

## 2. Data and methods

### a. Data

We estimated the parameter *β* and its significance from the following sources:

We used gridded precipitation and temperature from the University of East Anglia's TS2.1 datasets for WNA (20°–55°N, 125°–100°W). The TS2.1 product is derived from interpolating monthly station data to a 0.5° × 0.5° grid and is continuous from January 1901 to December 2002 (Mitchell and Jones 2005).

We took estimates from a gridded Palmer

*Z*index [a measure of precipitation minus evaporation (*P*−*E*)] before autocorrelation has been introduced) and PDSI (Palmer 1965; Alley 1984) calculated from Climate Research Unit (CRU) precipitation and temperature in WNA.Simulated soil moisture of the twentieth century using the Variable Infiltration Capacity (VIC) model (Liang et al. 1994; Andreadis et al. 2005) provided estimates. The forcing for this data was gridded (0.5° × 0.5° resolution) observations from the U.S. Historical Climate Network (HCN). We report

*β*values here for the second layer in the VIC model, which represents soil depths of 0.2–2.4 m. The spatial domain of data available in WNA from this product covers the coterminous United States.We estimated from gridded (2.5° × 2.5°) tree-ring-based reconstructions of PDSI from the North American Drought Atlas (Cook et al. 2004) spanning 1000–2000 CE. Only grid points in WNA were considered.

We included tree-ring-based reconstructions of hydroclimate variables other than PDSI from the WNA region (Tables 1 and 2). Each of these reconstructions relied on tree-ring data that were preprocessed to remove variability unrelated to climate as described in original references. All data and references are available on the National Oceanic and Atmospheric Administration (NOAA) paleoclimatology web page (http://www.ncdc.noaa.gov/paleo/recons.html).

We took estimates from 10 multicentury moisture-sensitive tree-ring chronologies from a recently developed five needle pine database (Kipfmueller and Salzer 2010). Chronologies were selected if they 1) showed a significant (

*p*value < 0.01) relationship with local precipitation; 2) had at least 900 years of data between 1000 and 2000 CE; and 3) were at least 250 m below tree line, which we required to ensure stable sensitivities through time (e.g., Salzer et al. 2009). We use long chronologies (with high sample depth) to minimize the effects of the “segment length curse,” which limits the amount of variance retained by chronologies at long time scales (i.e., Cook et al. 1995). Thus chronologies constructed from long segments will preserve more low-frequency variance than those constructed from shorter segments. We normalized each chronology to exhibit unit variance over the last 200 years for comparison with modern records. This simple rescaling ensures that variance outside the range of the last 200 years is represented by spectral energies that are correspondingly greater, but it does not alter the shape of the power spectrum.We also analyzed four high-resolution records from lakes and caves (Table 1). These records were included because they were resolved at decadal (or higher) time scales, spanned at least 1000 years, and extended past AD 1800. As with the tree-ring chronologies, each time series was normalized to exhibit unit variance over the last 200 years.

Finally, we estimated from model output from seven “last millennium” and “preindustrial control” simulations performed by different modeling centers throughout the world. The last millennium experiments were conducted using time-evolving boundary conditions described in Schmidt et al. (2012) and include 1) the gradual evolution of seasonal insolation from orbital forcing; 2) volcanic and solar activity; and 3) anthropogenic influences on greenhouse gas concentrations in the atmosphere, as well as surface characteristics on land. Experiments typically covered the period from 850 to 2005, although a few were shorter. The preindustrial control simulations are run without any modifications to their boundary conditions, which isolates the role of internally generated components of variability in the models.

### b. Spectral analysis

We estimated power spectra from annual indicators of instrumental and simulated twentieth-century hydroclimate using the multitaper method (Thomson 1982). Paleoclimate data (excluding tree rings) are unevenly sampled in time and so we used the Lomb–Scargle method (Lomb 1975; Scargle 1982) to estimate their spectra to minimize frequency bias (Schulz and Stattegger 1997). We normalized all paleoclimate data to exhibit unit variance over the last 200 years prior to spectral analysis. As in Huybers and Curry (2006), we estimated the slope (*β*) of individual power spectra by performing a least squares regression of log-transformed spectral density against log-transformed frequency. To avoid biasing the regression toward the higher frequencies (Huybers and Curry 2006), we averaged spectral density estimates to evenly spaced log-frequency bins prior to calculating *β*. We ensured that the values of *β* were not unduly influenced by long-term trends, which would artificially inflate variance at the lowest frequencies by performing the regression over the frequency domain *f _{q}* to 2/

*N*, where

*f*is the highest resolvable frequency (the Nyquist frequency) and

_{q}*N*is the length of the record in years. This step effectively sets the upper limit of our analyses to half the dataset length (e.g., 50-yr periods for instrumental data). An example of how our methodology is employed is provided in Fig. 1.

### c. A null hypothesis for the hydroclimate continuum

If the underlying continuum of hydroclimate in WNA reflects primarily stochastic climate processes (e.g., “weather noise” or uncorrelated year-to-year variations), then its spectrum would be similar to that of “white noise” and characterized by a spectral slope close to zero. However, positive spectral slopes are expected to arise in some hydroclimate indices as an artifact of how they record information or how they are calculated. For example, the PDSI has a built-in autocorrelation term of 0.897 at the monthly time step (Alley 1984). Although this number is specified a priori in PDSI, it agrees well with autocorrelation in simulated soil moisture during the twentieth century (Fig. 2). Likewise, biological, geochemical, and physical archives of past droughts integrate hydroclimate variability through time: tree-ring widths during a given year may be influenced by previous years (through surface processes that retain moisture or biological processes that store carbohydrates); lakes and caves may integrate information on interannual-to-decadal time scales (e.g., Truebe et al. 2010). Importantly, strong autocorrelation will attenuate high-frequency fluctuations while enhancing low-frequency variance. Here we examine this tendency by testing the possibility that local autocorrelation related to “nonclimate” processes accounts for the low-frequency energy in hydroclimate indicators.

We test the values of *β* by calculating against the null hypothesis that they are not significantly different from those characterizing white noise. To do so requires the null hypothesis to be applied differently depending on the variable being examined and the frequency domain being considered. For example, precipitation and the Palmer *Z* index (a metric of *P* − *E*) do not have a “built-in” source of autocorrelation, hence we evaluate the significance of *β* calculated from these variables as the *p* value of the linear regression of log-transformed spectral density [*S*(*f*)] against log-transformed frequency (*f*).

In contrast to precipitation and the *Z* index, moisture balance indicators such as PDSI, soil moisture, tree-ring records, and other paleoclimate records are subject to autocorrelation that arises locally and does not reflect low-frequency behavior in the climate system. Here, we develop a null hypothesis for the power spectrum of these moisture balance indicators by using the autocorrelation parameter that is built into the PDSI model (Alley 1984). Specifically, we are interested in determining if the PDSI's monthly autocorrelation term is sufficient to explain the power spectrum of simulated soil moisture, instrumental and reconstructed PDSI, reconstructed streamflow and precipitation, normalized paleoclimate records of moisture balance, and PDSI calculated from climate models. We begin by examining the expected power spectrum of an autoregressive process with lag-1 autocorrelation *ρ* (e.g., Bartlett 1978):

where *f* is frequency, *f _{N}* is the Nyquist frequency (in this case frequencies corresponding to 2-month wavelengths), and

*S*is mean spectral density, which is related to the white noise variance (

_{o}*σ*

^{2}) by

Figure 3a shows an example of this analytical expectation with unit variance white noise input. In most observational datasets the white noise variance (*σ*^{2}) is unknown, but the time series have been normalized to unit variance a posteriori. In this case, the value of *S _{o}* will be equal to unity (Fig. 3b). However, here we will be considering time series that have been annualized before being normalized and that are assumed to exhibit monthly autocorrelation close to 0.897 (i.e., the built-in autocorrelation in the PDSI model). In this case, the spectral density on interannual time scales will be given by

where only frequencies corresponding to wavelengths of 2 yr or longer (0 < *f*_{yr} ≤ 0.5) are considered; *S _{o}* is given by Eq. (2) and 〈

*S*(

*f*)〉 is the spectral average of

*S*(

*f*) in Eq. (1). This spectrum is shown in Fig. 3c. Note that autocorrelation, even at the monthly time scale, reddens the power spectrum out to decadal (10 yr) frequencies. We therefore consider the slope of the power spectrum on interannual time scales (

*β*) separately from decadal and longer time scales (

_{I}*β*).

_{D}We used a Monte Carlo (MC) procedure to establish confidence limits for the power spectra of moisture balance indicators and for the overall value of *β*, the interannual spectral slope *β _{I}*, and the decadal spectral slope

*β*. First, we generated 1000 realizations of a random autoregressive process with monthly autocorrelation set to 0.897. Each realization was run for 12 000 months, averaged to produce synthetic annual time series and normalized to exhibit unit variance. We then computed the power spectrum of each of these realizations to generate upper and lower 95% confidence limits of the variance in the power spectrum (Fig. 4a) as well as the spectral slope (

_{D}*β*) on interannual and decadal time scales (Fig. 4b).

## 3. Results

### a. Monte Carlo experiments

Results from our Monte Carlo analysis provide upper and lower confidence limits for the values of *β* calculated from an annually averaged unit variance time series with monthly autocorrelation equal to 0.897. Overall, the upper and lower 95% confidence limits of *β* are 0.46 and 0.08, respectively. On interannual time scales *β _{I}* these limits are 1.61 and 1.0, and on decadal

*β*time scales the limits are −0.32 and 0.34. We deem values outside of these limits significant. Figures 3 and 4 help illustrate the differences in the overall, interannual, and decadal expectations of

_{D}*β*: autocorrelation reddens the spectrum at 1/2–1/10-yr frequencies (2–10-yr periods).

### b. Instrumental data

Frequency-averaged spectral densities from WNA observational data (precipitation, *Z* index, PDSI, and soil moisture) are shown in Fig. 5. As in Pelletier and Turcotte (1997), we find that the spectrum of precipitation does not exhibit much redness (*β _{I}* = 0.03) on interannual time scales. We also observe that the spectra of the

*Z*index and precipitation are virtually identical. PDSI, in contrast, exhibits a much steeper (albeit nonsignificant) slope at interannual time scales (

*β*= 0.74), as would be expected from the analytical AR(1) red noise spectrum (Fig. 5c, dotted line). The change in slope is not as dramatic for soil moisture (

_{I}*β*= 0.27,

_{I}*β*= 0.2). On longer time scales, the slope of precipitation,

_{D}*Z*index, PDSI, and soil moisture (

*β*= 0.27, 0.26, 0.33, and 0.2, respectively) are more similar (Table 3). The higher value of the PDSI slope (0.33) reflects the slight attenuation of high frequencies from autocorrelation.

_{D}Values of *β* are mapped out across WNA for precipitation, *Z* index, PDSI, and soil moisture (Fig. 6). Values of *β* calculated from precipitation and the *Z* index are significantly different from zero (*p* value < 0.1) throughout much of the region (Fig. 6, top left). Values of *β* calculated from PDSI and *Z* index are also significantly different from the range of values expected from our AR(1) null hypothesis (Fig. 6, top right). At interannual time scales, significant values of *β _{I}* are restricted to small regions of the Pacific Northwest and northern Mexico in precipitation and

*Z*index (Fig. 6, middle left). The values of

*β*calculated from PDSI and soil moisture are not significant throughout most of WNA. In contrast, at long time scales there are large regions where the value of

_{I}*β*is significantly different from zero in precipitation and

_{D}*Z*index and significantly different from the AR(1) null hypothesis in PDSI and

*Z*index (Fig. 6, bottom). Importantly, the maps for precipitation,

*Z*index, and PDSI are nearly identical. Positive values (

*β*> 0.5) occur in the Southwest, the Rocky Mountains, and the northern U.S. Great Plains. The pattern of

_{D}*β*is similar for soil moisture (Fig. 6), implying that it is not an artifact of the drought index.

_{D}Trends in precipitation may influence the apparent spectral slope, but we removed frequencies lower than 1/(2*N*) from our calculation of *β _{D}* to avoid this bias. We also computed decadal spectral slopes (

*β*) for precipitation from different twentieth-century time windows, detrended time series, Historical Climate Network station data, standardized precipitation indices (SPI), and wet-day frequency (not shown). Each of these analyses supported a mean value of

_{D}*β*that was positive and significant for WNA precipitation (consistent with Fig. 6).

_{D}### c. Paleoclimate reconstructions and archives

The frequency-averaged power spectrum for tree-ring reconstructed PDSI is shown in Fig. 7a. The spectral slope of the average reconstructed PDSI power spectrum is not significantly different (at the 95% confidence limit) from the AR(1) expectation at interannual (*β* = 0.70) or decadal (*β _{D}* = −0.08) time scales. Overall the value is close to zero (

*β*= 0.09), which is also nonsignificant (at the 95% confidence level). Also shown in Fig. 7 are maps of

*β*,

*β*, and

_{I}*β*calculated from WNA. They indicate that the positive slopes are not significant overall; nor are they significant at interannual and decadal time scales (Figs. 7b–d). Significantly

_{D}*low*slopes occur in the westernmost part of the domain, indicating that interannual variability is energetic in this area (blue and white shading in Fig. 7c).

Figure 8 shows spectral averages from reconstructed normalized precipitation and streamflow, as well as other paleoclimate indicators from across WNA. At high frequencies, average spectra from normalized reconstructions and from the long chronologies are more energetic than the AR(1) expectation. At lower frequencies, the reconstructions are well within the AR(1) expectation, but the mean spectra from the longest 5-needle pine chronologies are not. Spectral variance in other (non-tree-ring) paleoclimate records (Fig. 8, symbols) is within the AR(1) expectation on decadal-to-100-yr time scales. At longer time scales, however, there is more variance in these records than the null hypothesis alone predicts (Tables 4 and 5). The variances and slopes of the paleoclimate records are in good agreement with the variance and slope records of the 5-needle pine chronologies.

### d. CMIP5 last millennium simulations

The PDSI spectra calculated from GCM simulations are within null expectations regardless of whether or not they are forced with time-evolving boundary conditions (Fig. 9). Likewise, values of *β* (Table 6) are nonsignificant overall, at interannual as well as decadal time scales. While unforced simulations are relatively flat at decadal and longer time scales (Fig. 9a), a slight (nonsignificant) reddening can be seen in the forced simulations of the last millennium (Fig. 9b). However, this reddening is evidently due mostly to trends occurring in the twentieth century because spectra calculated from the period 1000–1850 CE are flat at decadal and longer time scales (Fig. 9c) with one exception (MIROC-ESM). Considering twentieth-century simulations in isolation, however, does not provide evidence of redness beyond what is expected from the null hypothesis (Fig. 9d).

## 4. Discussion

We find WNA hydroclimate records exhibit more variance at longer time scales, with *β* between 0.2 and 0.8. This range is broadly consistent with previous studies (e.g., Pelletier and Turcotte 1997), however, *β* varies with hydroclimate indicator, geography, and frequency domain. These results are consistent with the idea that scaling arises from a variety of sources (Klemes 1974) and is not spatially uniform (Kantelhardt et al. 2006). For example, precipitation exhibits almost no persistence on interannual time scales, yet the autocorrelation of the PDSI reddens the spectrum considerably at these frequencies in observations (Fig. 5) as well as in models. Although the built-in PDSI autocorrelation coefficient was derived from local observations of soil moisture at a small spatial scale (Palmer 1965; Alley 1984), it is within the range of simulated soil moisture autocorrelation in WNA (Fig. 2). If this estimate is indeed realistic for WNA, then positive power laws on 2–10-yr time scales in PDSI, soil moisture, and tree-ring records may arise solely from local autocorrelation. This finding is important because it implies that *interannual* persistence in many drought indicators does not exclusively reflect remote sources of climate variability.

PDSI spectra calculated from tree-ring-based reconstructions and GCM simulations of the last millennium are fundamentally in agreement with each other and consistent with the AR(1) null hypothesis for hydroclimate variability. However, this agreement is misleading and the underlying spectrum of hydroclimate variability may indeed be much redder on decadal-to-centennial time scales than indicated by either reconstructed or simulated PDSI. Evidence for this interpretation comes from the spectrum of the longest moisture-sensitive chronologies in the 5-needle pine database and from the spectra of several other high-resolution paleodrought indicators, as well as instrumental data (Tables 3–5; Fig. 8). Moreover, an AR(1) process cannot account for the redness observed in twentieth-century instrumentally based PDSI and soil moisture at time scales longer than 10 years. The spectrum at these frequencies is relatively free of local autocorrelation effects, and the spectral slope reflects the behavior of the underlying climate system. Finally, the null hypothesis cannot account for the low-frequency scaling of precipitation (Fig. 5a) because precipitation fluctuations are not autocorrelated at time scales longer than a few days (Kantelhardt et al. 2006). The importance of redness at decadal-to-centennial time scales in observations and proxies supports our argument that this behavior is a robust feature of precipitation in western North America.

State-of-the-art GCMs appear unable to simulate decadal-to-centennial variability beyond what would be expected from local autocorrelation. This is apparent when the models are run with or without dynamical boundary conditions, implying that redness on these time scales is either an internally generated feature of WNA hydroclimate that the GCMs do not currently simulate or a forced response to external boundary conditions that is not captured by these models. On the other hand, limitations of the existing paleoclimate archive cannot be ignored: reconstructed forcing fields may be unrealistic, causing models to simulate forced responses that are inconsistent with the true evolution of climate during the last millennium or nonclimate processes that we do not yet understand may be introducing variability into proxy records at decadal and longer time scales. If the disagreement between proxies and GCMs is indeed driven by climatic variability that models do not simulate, then this finding has critical implications for assessing future risks because these models are actively used to understand climate change and its impacts.

## 5. Implications for drought risk

The differences between modeled hydroclimate and observational hydroclimate spectra are not very noticeable when considering only interannual-to-decadal time scales: the redness imposed by the PDSI in the model mimics redness in the observed system. However, at longer time scales, the discrepancies between modeled hydroclimate variability and proxy-inferred hydroclimate become much more important. To illustrate this point, we generated two different sets of 1000-yr Monte Carlo “red noise” time series with variance equal to the Southwest PDSI average. The first set was generated by calculating AR(1) time series with a monthly autocorrelation of 0.897. The second set was generated from 1000-yr realizations of random annually resolved white noise time series that were then filtered to exhibit spectral slopes that were at the lower (i.e., conservative) end of the instrumental and proxy range (*β* = 0.4–0.6). We generated these red noise realizations using the same technique employed by Pelletier and Turcotte (1997) and described more thoroughly in Pelletier (2008). We examine drought regimes across different time scales in these two sets of red noise time series by tallying the number of years in each realization when ⅗ (60%) of the antecedent years were below a given threshold. We then calculated the percentage of time each Monte Carlo realization spent at a given duration–magnitude combination and averaged these percentages across all realizations (Fig. 10).

At interannual-to-decadal time scales, the differences in drought risk distributions between AR(1) time series and those scaled by realistic spectral slopes are not very noticeable. However, at longer time scales, the two types of red noise data present very different views of drought risk. For instance, according to the AR(1) average distribution there is no risk of a 30- or 50-yr drought regime when 60% of the years experience PDSI values at or below −1. In contrast, the Monte Carlo realizations that have been rescaled by a realistic spectral slope tend to spend about 15% of the time in such a regime at the 30-yr time scale, and they spend about 10% of their time in this regime at the 50-yr time scale. We can also interpret these findings in terms of other hydroclimate variables, such as streamflow on the Colorado River, by rescaling the Monte Carlo realizations to exhibit the same variance as the variable of interest. The −1 threshold (1 standard deviation) would correspond to 4.37 million acre-feet, which, for reference, is roughly the state of California's entire lower basin allocation (Bureau of Reclamation 1948). Clearly, any decadal drought regime at this magnitude would pose unprecedented challenges to water resource managers. Our results suggest that this type of behavior is realistic for the region. Furthermore, inferences about future megadrought risk based on climate model simulations with an AR(1) distribution of spectral variance would therefore likely underestimate the risk of these events.

## 6. Conclusions

Our analysis uses a wide range of observed and simulated datasets to examine the distribution of variance in the hydroclimatology of WNA. By comparing spectral slopes against an appropriately defined null hypothesis of variance distribution, we are able to identify where hydroclimatic variance at low frequencies is significantly stronger than that at high frequencies. One of our most compelling findings is that nearly all observational datasets indicate substantial low-frequency variability at decadal and longer periods, with the exception of tree-ring reconstructions (in which such variance has been limited or reduced by preprocessing).

Several different processes likely explain the significant positive spectral slopes of North American hydroclimate. These processes vary depending on frequency interval and geography. At interannual time scales, precipitation and the Palmer *Z* index exhibit spectral slopes that are near zero (*β* = 0.03) and nonsignificant, whereas soil moisture and PDSI are considerably redder (*β* = 0.52 and *β* = 0.42, respectively). At decadal and lower frequencies, the spectrum of precipitation converges with the PDSI and soil moisture spectrum, implying that at low frequencies, terrestrial hydroclimate variability, and precipitation input are fundamentally similar and governed by climate. Moisture-sensitive tree-ring chronologies and lake and cave records all support significant, positive values of *β*. Such characteristics are generally not present in GCM simulations considered here. This finding is consistent with Ault et al. (2012), which found that the strength of decadal-to-multidecadal precipitation variability was underestimated in most CMIP5 simulations of the twentieth century. Considering that future megadroughts will likely emerge as consequences of both forced changes and natural variability, the relatively low amplitude of hydroclimate fluctuations in CMIP5 models suggests that risk of such events may be underestimated in the Intergovernmental Panel on Climate Change Fifth Assessment Report (IPCC AR5).

Future modeling efforts may help clarify the origins of hydroclimate redness in WNA on decadal-to-centennial time scales. One possibility is that differences in spectral slopes reflect the role of decadal-to-centennial variability in the tropical Pacific, which has been implicated as a major driver of WNA hydroclimate (e.g., Ropelewski and Halpert 1987; Seager et al. 2005a,b; Herweijer and Seager 2008) and yet is not well resolved in instrumental records (Cole et al. 1993; Urban et al. 2000; Cobb et al. 2003; Deser et al. 2004; Conroy et al. 2008; Ault et al. 2009). While efforts to reconcile the origins of WNA redness on decadal and longer time scales are under way, assessments of future drought risk may need to incorporate paleoclimatic insights if they are to encompass the full range of variability apparent in hydroclimatic observations.

## Acknowledgments

We thank Steve Gray, Adam Phillips, and Justin Sheffield for assistance and for data, and we are especially grateful to Dave Meko for sharing his Matlab code to calculate PDSI and the Palmer *Z* index. We appreciate the helpful reviews we received from Mike Evans, Zack Holden, and Andrew Ray. This project was supported in part by an NSF doctoral fellowship to T. Ault, a NCAR-ASP fellowship to T. Ault, NOAA CCDD (NA07OAR4310054) and NSF P2C2 (0903093) support for J. Cole, and NOAA Climate Program Office support for CLIMAS (J. Overpeck). Any use of trade, product, or firm names is for descriptive purposes only and does not imply endorsement by the U.S. Government.

## REFERENCES

*An Introduction to Stochastic Processes: With Special Reference to Methods and Applications.*CUP Archive, 388 pp.

*Geological Society of America Annual Meeting,*Denver, CO, Geological Society of America, 46-3. [Available online at https://gsa.confex.com/gsa/2005AM/finalprogram/abstract_97073.htm.]

*Climatic Variations and Forcing Mechanisms of the Last 2000 Years,*P. D. Jones et al., Eds., NATO ASI Series, Vol. 41, 109–124.

*Quantitative Modeling of Earth System Processes.*Cambridge University Press, 304 pp.

^{18}O and application to speleothem records

## Footnotes

The National Center for Atmospheric Research is sponsored by the National Science Foundation.