15 AUGUST 2013

AULT ET AL.

5863

The Continuum of Hydroclimate Variability in Western North America
during the Last Millennium
TOBY R. AULT,* JULIA E. COLE,* JONATHAN T. OVERPECK,* GREGORY T. PEDERSON,1
SCOTT ST. GEORGE,# BETTE OTTO-BLIESNER,@ CONNIE A. WOODHOUSE,& AND CLARA DESER@* Department of Geosciences, The University of Arizona, Tucson, Arizona
1
U.S. Geological Survey, Bozeman, Montana#
Department of Geography, University of Minnesota, Minneapolis, Minnesota@
National Center for Atmospheric Research,** Boulder, Colorado&
Department of Geography, The University of Arizona, Tucson, Arizona(Manuscript received 11 December 2011, in final form 14 December 2012)
ABSTRACTThe 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 (b) 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 b 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 b 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 ismore low-frequency variability in hydroclimate than monthly autocorrelation alone can generate. Coupledmodel 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 withpaleoclimate 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 climatechange 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 thepast 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 inthe future, it would likely put unprecedented strains on
existing water resources. Understanding megadroughtrisk is therefore critical to successfully managing waterin the coming decades, which in turn motivates us tocharacterize the hydroclimatic continuum in WNA on** The National Center for Atmospheric Research is sponsoredby the National Science Foundation.

Corresponding author address: Toby R. Ault, National Center
for Atmospheric Research, Boulder, CO 80307.
E-mail: ault@ucar.edu
DOI: 10.1175/JCLI-D-11-00732.1
Ó 2013 American Meteorological Society

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 risinggreenhouse 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 fullrange of future risks originating from external and internal sources of variability.
The paleoclimate record reveals that megadroughts arecommon in WNA over the last millennium: Woodhouseand 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

5864

JOURNAL OF CLIMATE

recently as the sixteenth century. Tree-ring-based reconstructions of flow on the Colorado River over the last1200 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 twentiethcentury—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 ofother 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 matchingor 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 headwatersof 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 thedefinition of hydrological extremes: if the model precipitation is greater than observed in a given region, thensimulated ‘‘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 fromobservational data, paleoclimate records, and climate
model output. In particular, we estimate the power-lawscaling parameter (b) that relates hydroclimate varianceto time scale such that S( f ) } f 2b where f is frequency
and S( f ) is the power spectrum (e.g., Pelletier andTurcotte 1997; Huybers and Curry 2006). This approachallows 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 b we estimate as ‘‘spectral
slopes.’’ These slopes give us insight into the underlying

VOLUME 26

relationship between time scale and variance acrossdisparate types of data.
Values of b near zero describe ‘‘white’’ spectra, whichhave an approximately uniform distribution of varianceacross 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 b, on theother hand, describe ‘‘blue’’ spectra. Such spectra canarise in systems with negative autocorrelation, becausefrom one time step to the next the system is pushed awayfrom its previous value.
Power spectra of temperature fluctuations in the ocean
and atmosphere are characterized by spectral slopes (b)between 0 and 0.5 on land and near 1.0 for most of theocean (Fraedrich and Blender 2003; Huybers and Curry
2006). More recently, Vyushin and Kushner (2009) and
Vyushin et al. (2009) have used global climate modeland reanalysis data to link positive spectral slopes in the
troposphere to tropical sea surface temperature (SST)
variability and similar behavior in the stratosphere tovolcanic activity. Hydroclimate variables have also beenshown to exhibit similar values, with b estimated near
0.5 for globally averaged spectra of river discharge, treering widths, and precipitation records (Pelletier andTurcotte 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 andCurry 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 moisturebalance. Although Pelletier and Turcotte (1997) used
long tree-ring chronologies, their approach did not screenthese records for moisture sensitivity. Further, in ourstudy, 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 havebeen forced with the time-evolving boundary conditionsof 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 themost scientifically advanced tools for numerically simulating climate variability and change, and they are being
used throughout the world to anticipate and prepare forclimate change.

15 AUGUST 2013

AULT ET AL.

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 withrealistic autocorrelation occurring on the monthly timescale. 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 thatremove 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 timescales, 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 b and its significancefrom the following sources:(i) We used gridded precipitation and temperaturefrom the University of East Anglia’s TS2.1 datasets for WNA (208–558N, 1258–1008W). The TS2.1
product is derived from interpolating monthly station data to a 0.58 3 0.58 grid and is continuous fromJanuary 1901 to December 2002 (Mitchell andJones 2005).
(ii) We took estimates from a gridded Palmer Z index[a measure of precipitation minus evaporation(P 2 E)] before autocorrelation has been introduced) and PDSI (Palmer 1965; Alley 1984)
calculated from Climate Research Unit (CRU)
precipitation and temperature in WNA.
(iii) 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.58 3 0.58 resolution) observations fromthe U.S. Historical Climate Network (HCN). We
report b 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.
(iv) We estimated from gridded (2.58 3 2.58) tree-ringbased reconstructions of PDSI from the North

5865

American Drought Atlas (Cook et al. 2004) spanning 1000–2000 CE. Only grid points in WNA wereconsidered.
(v) We included tree-ring-based reconstructions of hydroclimate variables other than PDSI from the
WNA region (Tables 1 and 2). Each of thesereconstructions relied on tree-ring data that werepreprocessed to remove variability unrelated toclimate 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).
(vi) We took estimates from 10 multicentury moisturesensitive 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 atleast 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 throughtime (e.g., Salzer et al. 2009). We use long chronologies (with high sample depth) to minimize theeffects 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 overthe last 200 years for comparison with modern
records. This simple rescaling ensures that variance outside the range of the last 200 years isrepresented by spectral energies that are correspondingly greater, but it does not alter the shapeof the power spectrum.
(vii) 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 withthe tree-ring chronologies, each time series was
normalized to exhibit unit variance over the last
200 years.
(viii) Finally, we estimated from model output fromseven ‘‘last millennium’’ and ‘‘preindustrial control’’ simulations performed by different modelingcenters 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

5866

JOURNAL OF CLIMATEVOLUME 26

TABLE 1. Estimates of the parameter b (see section 2c) from a variety of precipitation reconstructions in WNA. Bolded values are
significant against the AR(1) null hypothesis at the 95% confidence limits. In the third column, N refers to the length of the time series.
RegionReference

N/2

b

bI

bD

Durango, CO
Chihuahua, Mexico
D7, Oregon
D8, Oregon
Bighorn, WY
Uinta, UT
Malpais, NM
D6, Nevada
D1, Arizona
D2, Arizona
D3, Arizona
D4, Arizona
D5, Arizona
D6, Arizona
D7, Arizona
D1, New Mexico
D2, New Mexico
D3, New Mexico
D4, New Mexico
D5, New Mexico
D6, New Mexico
D7, New Mexico
D8, New Mexico
Saltillo, Mexico
Colorado Plateau

Cleaveland et al. (2003)
Diaz et al. (2002)
Garfin and Hughes (1996)
Garfin and Hughes (1996)
Gray et al. (2004a)
Gray et al. (2004b)
Grissino-Mayer (1995)
Hughes and Gramlich (1996)
Ni et al. (2002)
Ni et al. (2002)
Ni et al. (2002)
Ni et al. (2002)
Ni et al. (2002)
Ni et al. (2002)
Ni et al. (2002)
Ni et al. (2002)
Ni et al. (2002)
Ni et al. (2002)
Ni et al. (2002)
Ni et al. (2002)
Ni et al. (2002)
Ni et al. (2002)
Ni et al. (2002)
Pohl et al. (2003)
Salzer and Kipfmueller (2005)

304
173
138
138
370
388
1065
1498
495
495
495
495
495
495
495
495
495
495
495
495
495
495
495
110
709

0.15
20.15
0.25
0.28
20.32
0.42
0.53
0.12
0.35
0.21
0.42
0.47
0.29
0.3
0.6
0.32
0.52
0.49
0.53
0.62
0.56
0.53
0.54
0.07
20.06

0.08
0.23
0.27
0.47
0.08
0.83
0.1
0.05
0.23
0.69
0.57
0.79
0.46
0.45
0.53
0.59
0.62
0.4
0.59
0.49
0.42
0.41
0.48
0.13
0.55

0.28
20.19
0.45
0.4
20.8
0.21
0.55
0.22
0.37
0
0.3
0.31
0.16
0.21
0.52
0.11
0.48
0.39
0.48
0.48
0.41
0.41
0.46
0.47
20.51

and solar activity; and 3) anthropogenic influenceson 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, whichisolates 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 theirspectra 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 estimatedthe slope (b) of individual power spectra by performinga 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 b. We ensured that the values of b were notunduly influenced by long-term trends, which would

TABLE 2. Estimates of b, bI, and bD from several tree-ring reconstructions of streamflow for watersheds in western North America.
Bolded values indicate significance at the 95% confidence level.
Site

Reference

N/2

b

bI

bD

Yellowstone River, MT
Clear Creek, CO
Salinas, CA
Boulder Creek, CO
Feather River, CA
Upper Colorado River
Lees Ferry, AZ
Sacramento River, CA

Graumlich et al. (2003)
Woodhouse (2000)
Griffin (2005)
Woodhouse (2001)
Meko (2001)
Meko et al. (2007)
Woodhouse et al. (2006)
Meko et al. (2001)

147
152
298
143
539
622
255
555

20.76
20.15
0.1
0.28
20.55
0.37
0.52
0.05

0.36
0.23
0.27
0.47
0.08
0.83
0.1
0.05

21.76
20.19
0.05
0.4
21.09
0.16
0.56
0.11

15 AUGUST 2013

AULT ET AL.

5867

FIG. 1. (left) Time series and (right) spectra for hydroclimate indicators averaged for ‘‘Montana’’ (458–508N, 1148–1058W) as an example of how b is calculated from (top–bottom) annually
averaged precipitation, Z index, PDSI, and soil moisture. Gray lines show the original monthly
time series and spectra and black lines show the time series and spectra of the annual average.
Time series were normalized to unit standard deviation before computing the spectrum. Annual
spectra are offset from monthly spectra because they exhibit less variance. For reference, theanalytical AR(1) spectrum is shown on the bottom right (thin, dotted line). Thick dashed lines onthe (right) show the overall spectral slope (b) calculated from 1/ 2- to 1/ 50-yr frequencies.

artificially inflate variance at the lowest frequencies byperforming the regression over the frequency domain
fq to 2/N, where fq is the highest resolvable frequency(the Nyquist frequency) and N is the length of the record in years. This step effectively sets the upper limitof our analyses to half the dataset length (e.g., 50-yr
periods for instrumental data). An example of how ourmethodology 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 arecalculated. 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 archivesof 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 thistendency by testing the possibility that local autocorrelation related to ‘‘nonclimate’’ processes accounts forthe low-frequency energy in hydroclimate indicators.
We test the values of b by calculating against the null
hypothesis that they are not significantly different from

5868

JOURNAL OF CLIMATEVOLUME 26

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 issufficient 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 processwith lag-1 autocorrelation r (e.g., Bartlett 1978):1 2 r2

S( f ) 5 So

1 2 2r cosp
FIG. 2. Monthly autocorrelation parameter estimated froma VIC soil moisture model simulation of the twentieth century
using instrumental data (Andreadis et al. 2005).

those characterizing white noise. To do so requires thenull hypothesis to be applied differently depending on
the variable being examined and the frequency domain
being considered. For example, precipitation and thePalmer Z index (a metric of P 2 E) do not have a ‘‘builtin’’ source of autocorrelation, hence we evaluate thesignificance of b calculated from these variables as thep 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, moisturebalance indicators such as PDSI, soil moisture, tree-ringrecords, and other paleoclimate records are subject toautocorrelation that arises locally and does not reflect
low-frequency behavior in the climate system. Here, wedevelop a null hypothesis for the power spectrum of

f
1 r2
fN

,

(1)

where f is frequency, fN is the Nyquist frequency (in thiscase frequencies corresponding to 2-month wavelengths),and So is mean spectral density, which is related to thewhite noise variance (s2) by
So 5

s2
.
1 2 r2(2)

Figure 3a shows an example of this analytical expectation with unit variance white noise input. In mostobservational datasets the white noise variance (s2) is
unknown, but the time series have been normalized tounit variance a posteriori. In this case, the value of So
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., thebuilt-in autocorrelation in the PDSI model). In thiscase, the spectral density on interannual time scales
will be given by

FIG. 3. Examples of the analytical AR(1) spectrum to illustrate the effects of different
normalizations. (a) An AR(1) process with unit white noise variance and autocorrelation of
0.897 at the monthly time step [Eq. (1)]. (b) An AR(1) process that has been normalized to
exhibit unit variance overall [Eq. (2)]. (c) An AR(1) process with monthly autocorrelation of
0.897 that has been averaged to annual resolution and normalized to unit variance [Eq. (3)].

15 AUGUST 2013

AULT ET AL.

S( fyr ) 5SohS( f )i

1 2 r2,fyr
2
1 2 2r cosp 1 r
fN

5869

(3)

where only frequencies corresponding to wavelengths of
2 yr or longer (0 , fyr # 0.5) are considered; So is given
by Eq. (2) and hS( f )i 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 (bI) separately from decadal
and longer time scales (bD).
We used a Monte Carlo (MC) procedure to establish
confidence limits for the power spectra of moisturebalance indicators and for the overall value of b, theinterannual spectral slope bI, and the decadal spectral
slope bD. 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 thencomputed 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 (b) on interannual and decadal
time scales (Fig. 4b).

3. Resultsa. Monte Carlo experiments
Results from our Monte Carlo analysis provide upperand lower confidence limits for the values of b calculated
from an annually averaged unit variance time series withmonthly autocorrelation equal to 0.897. Overall, theupper and lower 95% confidence limits of b are 0.46 and
0.08, respectively. On interannual time scales bI theselimits are 1.61 and 1.0, and on decadal bD time scales thelimits are 20.32 and 0.34. We deem values outside ofthese limits significant. Figures 3 and 4 help illustrate thedifferences in the overall, interannual, and decadal expectations of b: 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, andsoil moisture) are shown in Fig. 5. As in Pelletier andTurcotte (1997), we find that the spectrum of precipitation does not exhibit much redness (bI 5 0.03) oninterannual time scales. We also observe that the spectraof the Z index and precipitation are virtually identical.

FIG. 4. Results of Monte Carlo simulations of an AR(1) process
to establish confidence limits for the spectral slope overall (b) andat interannual (bI) and decadal (bD) time scales. (a) The upper andlower 95% confidence limits from all Monte Carlo spectra (gray)shown with the analytical expectation (black line). (b) The distribution of bI values at interannual time scales (2–10 yr). (c) The
distribution of bD on .10-yr time scales.

PDSI, in contrast, exhibits a much steeper (albeit nonsignificant) slope at interannual time scales (bI 5 0.74),
as would be expected from the analytical AR(1) red
noise spectrum (Fig. 5c, dotted line). The change inslope is not as dramatic for soil moisture (bI 5 0.27,
bD 5 0.2). On longer time scales, the slope of precipitation, Z index, PDSI, and soil moisture (bD 5 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 theslight attenuation of high frequencies from autocorrelation.
Values of b are mapped out across WNA for precipitation, Z index, PDSI, and soil moisture (Fig. 6).

5870

JOURNAL OF CLIMATEVOLUME 26

TABLE 3. Summary of spectral slopes calculated from mean
spectra for WNA (208–558N, 1258–1008W). Slopes were calculated
as the linear least squares regression between frequency andspectral density. Estimates of b were made from 2- to 50-yr timescales (i.e., the shortest resolvable frequencies to frequencies corresponding to one-half of the lowest resolvable wavelength), as
well as from 2- to 10-yr time scales (bI) and 10- to 50-yr time scales(bD). In the case of precipitation and the Z index, significance ofthe slope was determined by the p value of the regression line.
Significance of the slope of PDSI and soil moisture was determinedusing a Monte Carlo method. Bolded values indicate they weresignificant at the 95% confidence limit or higher.
Variable

SourcebbI

bD

Precipitation

CRU TS2.1 (Mitchelland Jones 2005)
CRU TS2.1 (Mitchelland Jones 2005)
CRU TS2.1 (Mitchelland Jones 2005)
VIC (Andreadis et al. 2005)

0.2

0.03

0.27

0.19

0.03

0.22

0.52

0.74

0.33

0.42

0.27

0.2

Z index
PDSI
Soil moisture

FIG. 5. MNA spectra (gray lines) and corresponding values of b
computed over interannual and decadal frequencies (black lines)for (a) CRU TS2.1 precipitation, (b) the Palmer Z index (a measure of the P 2 E balance) calculated from CRU TS2.1 precipitation and temperature, (c) PDSI also calculated from CRU
TS2.1 data, and (d) soil moisture simulated by the VIC model(Andreadis et al. 2005). Diamonds indicate the mean value of thespectrum at the points used to perform regression calculations, andthey are evenly spaced in log-transformed unit frequency. The
dotted line on panel (c) shows the theoretical AR(1) spectrum forred noise with autocorrelation equal to the built-in monthly autocorrelation of the PDSI.

Values of b 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). Valuesof b 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 bI are restricted to small regions of the Pacific Northwest andnorthern Mexico in precipitation and Z index (Fig. 6,
middle left). The values of bI calculated from PDSI andsoil moisture are not significant throughout most of
WNA. In contrast, at long time scales there are large
regions where the value of bD is significantly different
from zero in precipitation and 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 (bD . 0.5) occur in the Southwest, theRocky Mountains, and the northern U.S. Great Plains.
The pattern of bD is similar for soil moisture (Fig. 6),
implying that it is not an artifact of the drought index.
Trends in precipitation may influence the apparent
spectral slope, but we removed frequencies lower than
1/(2N) from our calculation of bD to avoid this bias. We
also computed decadal spectral slopes (bD) for precipitation from different twentieth-century time windows,
detrended time series, Historical Climate Network station data, standardized precipitation indices (SPI), andwet-day frequency (not shown). Each of these analyses
supported a mean value of bD that was positive and significant for WNA precipitation (consistent with Fig. 6).

c. Paleoclimate reconstructions and archives
The frequency-averaged power spectrum for tree-ringreconstructed 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 (b 50.70) or decadal (bD 5 20.08) time scales. Overall thevalue is close to zero (b 5 0.09), which is also nonsignificant (at the 95% confidence level). Also shown in15 AUGUST 2013

AULT ET AL.

5871

FIG. 6. (top–bottom) Spectral slopes (b) calculated for (left–right) various instrumental datasets fromWNA (208–558N, 1258–1008W) at 2–50 yr, interannual (2–10 yr), and decadal (10–50 yr) time scales.

Fig. 7 are maps of b, bI, and bD 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 low slopes
occur in the westernmost part of the domain, indicatingthat interannual variability is energetic in this area (blueand white shading in Fig. 7c).
Figure 8 shows spectral averages from reconstructednormalized precipitation and streamflow, as well asother paleoclimate indicators from across WNA. At
high frequencies, average spectra from normalized reconstructions and from the long chronologies are moreenergetic 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 (nontree-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 inthese 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 varianceand 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 b (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 evidenceof redness beyond what is expected from the null hypothesis (Fig. 9d).

4. Discussion
We find WNA hydroclimate records exhibit morevariance at longer time scales, with b between 0.2 and0.8. This range is broadly consistent with previous
studies (e.g., Pelletier and Turcotte 1997), however, bvaries with hydroclimate indicator, geography, andfrequency 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5872

JOURNAL OF CLIMATE

(a)

(b)

(c)

(d)

VOLUME 26

FIG. 7. (a) Mean spectral density (i.e., the WNA average at each frequency) calculated fromtree-ring reconstructed PDSI (Cook et al. 2004). As in Fig. 5, the gray line shows the mean power
spectrum, the straight black lines show the spectral slopes at interannual and decadal frequencies(bI and bD), and the black diamonds show the spectral bins that were used to calculate thoseslopes. The thin dotted line shows the analytical AR(1) spectrum. (b) The overall value of b(there are no significant values at the 95% confidence limit). (c) Estimates of bI. Significant
spectral slopes are outlined in black (values to the west of this line are significant at the 95%
confidence limit). (d) The decadal spectral slope bD, which is not significant anywhere in WNA.

persistence on interannual time scales, yet the autocorrelation of the PDSI reddens the spectrum considerably atthese frequencies in observations (Fig. 5) as well as inmodels. 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 scalesin PDSI, soil moisture, and tree-ring records may arise
solely from local autocorrelation. This finding is importantbecause it implies that interannual persistence in many
drought indicators does not exclusively reflect remotesources of climate variability.
PDSI spectra calculated from tree-ring-based reconstructions and GCM simulations of the last millennium are fundamentally in agreement with each otherand consistent with the AR(1) null hypothesis forhydroclimate variability. However, this agreement ismisleading 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 thisinterpretation comes from the spectrum of the longest
moisture-sensitive chronologies in the 5-needle pine
database and from the spectra of several other highresolution paleodrought indicators, as well as instrumental data (Tables 3–5; Fig. 8). Moreover, an AR(1)process cannot account for the redness observed intwentieth-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, thenull 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 importanceof 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15 AUGUST 2013

5873

AULT ET AL.

FIG. 8. Analytical AR(1) spectrum (black) and upper and lower 95% Monte Carlo confidence limits shown with spectral averages from a variety of paleoclimate archives. The red line
shows the spectral average of reconstructed streamflow from 8 sites in WNA, and the purple
line shows the spectral average of 25 precipitation reconstructions. Mean spectral densityfrom 10 moisture-sensitive long 5-needle pine chronologies are shown in gold. Symbols show
the spectra of four paleoclimate time series that have been averaged into evenly spaced logtransformed bins. These records include: a diatom-inferred salinity record from Moon Lake,ND (Laird et al. 1998; circles); an oxygen isotope record from Pyramid Lake, CA (Benson et al.
2003; triangles); a speleothem band thickness record from Carlsbad Caverns, NM (Polyak andAsmerom 2001; stars); and a cave d18O record from Hidden Cave, NM (Rasmussen et al. 2006;
squares). The thick gray line is the average spectral density from the four paleoclimate records,
and the dashed black line shows the regression line for that average.

expected from local autocorrelation. This is apparent
when the models are run with or without dynamicalboundary conditions, implying that redness on thesetime 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 theother 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 ofclimate during the last millennium or nonclimate processes that we do not yet understand may be introducing
variability into proxy records at decadal and longer timescales. 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 andobservational hydroclimate spectra are not very noticeable when considering only interannual-to-decadal
time scales: the redness imposed by the PDSI in themodel mimics redness in the observed system. However, at longer time scales, the discrepancies between
TABLE 4. Estimates of bD from four high-resolution moisturesensitive paleoclimate records from WNA. Since the shortest resolvable frequencies in these records are decadal, bD and b would
be the same. For clarity, we report the values in terms of bD.
Bolded values indicate significance at the 95% confidence level.
Variable

ReferencebD

Moon Lake salinity
Pyramid Lake d18O
Carlsbad Cave band thickness
Hidden Cave band thickness

Laird et al. (1996)
Benson et al. (2002)
Polyak and Asmerom (2001)
Rasmussen et al. (2006)

0.84
1.16
0.63
0.82

5874

JOURNAL OF CLIMATETABLE 5. As in Table 3, but for paleoclimate records. Since thehighest resolvable frequencies from the lake and cave records usedhere correspond to 10-yr wavelengths, estimates of bD and b would
be the same. For clarity, we report these values in terms of bD.
Bolded values indicate significance at the 95% confidence level.
Variable

Reference

b

bI

bD

North American
Drought Atlas
Precipitation
reconstructions
Streamflow
reconstructions
5-needle pine
chronologies
WNA lakes andcaves

Cook et al. (2004)

0.09

0.7

20.08

See Table 1

0.1

0.24

0.01

See Table 2

20.02

0.19

20.14

0.59

0.3

0.57

—

—

0.86

Kipfmueller andSalzer (2010)
See Table 4

modeled hydroclimate variability and proxy-inferredhydroclimate 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

VOLUME 26

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 (b 5 0.4–0.6).
We generated these red noise realizations using the same
technique employed by Pelletier and Turcotte (1997) anddescribed more thoroughly in Pelletier (2008). We examine drought regimes across different time scales inthese two sets of red noise time series by tallying thenumber of years in each realization when 3/ 5 (60%) of
the antecedent years were below a given threshold. Wethen 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 differencesin drought risk distributions between AR(1) time series

FIG. 9. North American spectral means for PDSI computed from seven CMIP5 models thathave run last millennium experiments. In all panels, the smooth black line shows the theoreticalmean spectrum of PDSI, and the dashed straight line shows the regression line that was fit to themultiproxy spectrum in Fig. 8, for reference. Spectral means are shown for (a) control integrations for each model, (b) complete last millennium simulations, (c) the forced preindustrial portion (1000–1850 CE) of each simulation, and (d) the twentieth-century portion ofeach simulation.

15 AUGUST 2013

5875

AULT ET AL.

TABLE 6. As in Tables 3 and 5, but for PDSI calculated from climate model data. Values that are significant at the 95% confidence limit are
shown in bold.

Model
Community Climate System Model,
version 4.0 (CCSM4)
CCSM4
Flexible Global Ocean–Atmosphere–Land
System Model (FGOALS-gl)
Goddard Institute for Space Studies
Model E, coupled with Russell
ocean model (GISS-E2-R)
GISS-E2-R
L’Institut Pierre-Simon Laplace Coupled
Model, version 5, coupled with NEMO,
low resolution (IPSL-CM5A-LR)
IPSL-CM5A-LR
Model for Interdisciplinary Research on
Climate, Earth System Model(MIROC-ESM)
MIROC-ESM
Max Planck Institute Earth System Model,
paleo (MPI-ESM-P)
MPI-ESM-P
Beijing Climate Center, Climate System
Model, version 1.1 (BCC-CSM1.1)
BCC-CSM1.1

Approximate resolution(lat 3 lon)

Experiment

Nyr

b

bI

bD

18 3 1.258

past1000

1156

0.28

1.01

0.18

18 3 1.258
48 3 58

piControl
past1000

501
1000

0.17
0.12

0.83
0.9

20.02
20.04

28 3 2.58

past1000

1156

0.24

0.87

0.09

28 3 2.58
28 3 3.758

piControl
past1000

850
1156

0.19
0.26

0.97
0.9

0.05
0.13

28 3 3.758
2.88 3 2.58

piControl
past1000

1000
1156

0.18
0.31

0.85
0.94

0.05
0.19

2.58 3 2.58
1.98 3 1.98

piControl
past1000

630
1156

0.29
0.32

0.85
0.98

0.11
0.22

1.98 3 1.98
2.88 3 2.88

piControl
past1000

1156
1151

0.22
0.21

1.04
0.84

0.03
0.07

2.88 3 2.88

piControl

500

0.24

0.79

0

and those scaled by realistic spectral slopes are notvery noticeable. However, at longer time scales, thetwo types of red noise data present very different viewsof 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 21. In contrast, the Monte
Carlo realizations that have been rescaled by a realistic
spectral slope tend to spend about 15% of the time insuch a regime at the 30-yr time scale, and they spend
about 10% of their time in this regime at the 50-yr timescale. We can also interpret these findings in terms ofother hydroclimate variables, such as streamflow on
the Colorado River, by rescaling the Monte Carlo realizations to exhibit the same variance as the variableof interest. The 21 threshold (1 standard deviation)
would correspond to 4.37 million acre-feet, which, forreference, 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 ofspectral 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 variancein 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 ofour most compelling findings is that nearly all observational datasets indicate substantial low-frequency variability at decadal and longer periods, with the exceptionof 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 intervaland geography. At interannual time scales, precipitationand the Palmer Z index exhibit spectral slopes that arenear zero (b 5 0.03) and nonsignificant, whereas soil
moisture and PDSI are considerably redder (b 5 0.52
and b 5 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 governedby climate. Moisture-sensitive tree-ring chronologies

5876

JOURNAL OF CLIMATEVOLUME 26

of both forced changes and natural variability, the relatively low amplitude of hydroclimate fluctuations in
CMIP5 models suggests that risk of such events may beunderestimated in the Intergovernmental Panel on Climate Change Fifth Assessment Report (IPCC AR5).
Future modeling efforts may help clarify the origins ofhydroclimate 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 amajor driver of WNA hydroclimate (e.g., Ropelewski andHalpert 1987; Seager et al. 2005a,b; Herweijer and Seager2008) and yet is not well resolved in instrumental records(Cole et al. 1993; Urban et al. 2000; Cobb et al. 2003; Deseret al. 2004; Conroy et al. 2008; Ault et al. 2009). Whileefforts to reconcile the origins of WNA redness on decadal and longer time scales are under way, assessments offuture 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 toT. Ault, a NCAR-ASP fellowship to T. Ault, NOAA
CCDD (NA07OAR4310054) and NSF P2C2 (0903093)
support for J. Cole, and NOAA Climate Program Officesupport for CLIMAS (J. Overpeck). Any use of trade,
product, or firm names is for descriptive purposes only anddoes not imply endorsement by the U.S. Government.
FIG. 10. Duration vs magnitude diagram of Monte Carlo PDSI
simulations of droughts in the Southwest. In this case, a drought isdefined as a year in which 3/ 5 of the antecedent years exhibited
PDSI values below a given threshold. The thresholds are shown on
the x axis and the length of the time window considered is shown on
the y axis. The shading indicates the percentage of time (averaged
across MC realizations) spent at different duration–magnitude
combinations. (a) Average distribution for AR(1) time series.
(b) The average distribution for time series that have been rescaled
to exhibit a realistic spectral slope.

and lake and cave records all support significant, positive values of b. Such characteristics are generally
not present in GCM simulations considered here. Thisfinding is consistent with Ault et al. (2012), which foundthat 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

REFERENCES
Alley, W. M., 1984: The Palmer Drought Severity Index: Limitations
and assumptions. J. Climate Appl. Meteor., 23, 1100–1109.
Andreadis, K. M., E. A. Clark, A. W. Wood, A. F. Hamlet, and
D. P. Lettenmaier, 2005: Twentieth-century drought in theconterminous United States. J. Hydrometeor., 6, 985–1001.
Ault, T. R., J. E. Cole, M. N. Evans, H. Barnett, N. J. Abram, A. W.
Tudhope, and B. K. Linsley, 2009: Intensified decadal variability
in tropical climate during the late 19th century. Geophys. Res.
Lett., 36, L08602, doi:10.1029/2008GL036924.
——, ——, and S. St. George, 2012: The amplitude of decadal
to multidecadal variability in precipitation simulated by stateof-the-art climate models. Geophys. Res. Lett., 39, L21705,
doi:10.1029/2012GL053424.
Bartlett, M. S., 1978: An Introduction to Stochastic Processes: With
Special Reference to Methods and Applications. CUP Archive,
388 pp.
Benson, L., and Coauthors, 2002: Holocene multidecadal andmulticentennial droughts affecting northern California and
Nevada. Quat. Sci. Rev., 21, 659–682.

15 AUGUST 2013

AULT ET AL.

——, B. Linsley, J. Smoot, S. Mensing, S. Lund, S. Stine, and
A. Sarna-Wojcicki, 2003: Influence of the Pacific decadal oscillation on the climate of the Sierra Nevada, California and
Nevada. Quat. Res., 59, 151–159.
Bureau of Reclamation, 1948: Upper Colorado River Basin
Compact of 1948, 12 pp.
Clauset, A., C. R. Shalizi, and M. E. J. Newman, 2009: Power-law
distributions in empirical data. SIAM Rev., 51, 661–703,
doi:10.1137/070710111.
Cleaveland, M. K., D. W. Stahle, M. D. Therrell, J. VillanuevaDiaz, and B. T. Burns, 2003: Tree-ring reconstructed winterprecipitation and tropical teleconnections in Durango, Mexico.
Climatic Change, 59, 369–388.
Cobb, K. M., C. D. Charles, H. Cheng, and R. L. Edwards, 2003: El
Ni~no/Southern Oscillation and tropical Pacific climate during
the last millennium. Nature, 424, 271–276.
Cole, J. E., R. G. Fairbanks, and G. T. Shen, 1993: Recent variability
in the Southern Oscillation—Isotopic results from a Tarawa
atoll coral. Science, 260, 1790–1793.
Conroy, J. L., J. T. Overpeck, J. E. Cole, T. M. Shanahan, and
M. Steinitz-Kannan, 2008: Holocene changes in eastern tropical Pacific climate inferred from a Galapagos lake sedimentrecord. Quat. Sci. Rev., 27, 1166–1180.
Cook, E. R., K. R. Briffa, D. M. Meko, D. A. Graybill, andG. Funkhouser, 1995: The segment length curse in long treering chronology development for paleoclimatic studies. Holocene, 5, 229–237.
——, D. M. Meko, D. W. Stahle, and M. K. Cleaveland, 1999:
Drought reconstructions for the continental United States.
J. Climate, 12, 1145–1162.
——, C. A. Woodhouse, C. M. Eakin, D. M. Meko, and D. W.
Stahle, 2004: Long-term aridity changes in the western UnitedStates. Science, 306, 1015–1018.
——, R. Seager, R. R. Heim, R. S. Vose, C. Herweijer, andC. Woodhouse, 2010: Megadroughts in North America: Placing IPCC projections of hydroclimatic change in a long-term
palaeoclimate context. J. Quat. Sci., 25, 48–61.
Deser, C., A. S. Phillips, and J. W. Hurrell, 2004: Pacific interdecadal climate variability: Linkages between the tropics andthe North Pacific during boreal winter since 1900. J. Climate,17, 3109–3124.
Diaz, S. C., M. D. Therrell, D. W. Stahle, and M. K. Cleaveland,
2002: Chihuahua (Mexico) winter-spring precipitation reconstructed from tree-rings, 1647–1992. Climate Res., 22, 237–244.
Fraedrich, K., and R. Blender, 2003: Scaling of atmosphereand ocean temperature correlations in observations andclimate models. Phys. Rev. Lett., 90, 108501, doi:10.1103/
PhysRevLett.90.108501.
Garfin, G. M., and M. K. Hughes, 1996: Eastern Oregon divisional
precipitation and Palmer Drought Severity Index from treerings. U.S. Forest Service Intermountain Research Station
Rep., USDA Forest Service Cooperative Agreement PNW
90-174.
Graham, N. E., and Coauthors, 2007: Tropical Pacific—Mid-latitude
teleconnections in medieval times. Climatic Change, 83, 241–
285.
Graumlich, L. J., M. F. J. Pisaric, L. A. Waggoner, J. S. Littell, andJ. C. King, 2003: Upper Yellowstone River flow and teleconnections with Pacific basin climate variability during thepast three centuries. Climatic Change, 59, 245–262.
Gray, S. T., C. L. Fastie, S. T. Jackson, and J. L. Betancourt, 2004a:
Tree-ring-based reconstruction of precipitation in the Bighorn
Basin, Wyoming, since 1260 AD. J. Climate, 17, 3855–3865.

5877——, S. T. Jackson, and J. L. Betancourt, 2004b: Tree-ring based
reconstructions of interannual to decadal scale precipitation
variability for northeastern Utah since 1226 AD. J. Amer.
Water Resour. Assoc., 40, 947–960.
Griffin, R., 2005: Paleohydrology of the Salinas River, California,
reconstructed from blue oak tree rings. Geological Society ofAmerica Annual Meeting, Denver, CO, Geological Society ofAmerica, 46-3. [Available online at https://gsa.confex.com/
gsa/2005AM/finalprogram/abstract_97073.htm.]
Grissino-Mayer, H. D., 1995: The climate and fire history of El
Malpais national monument, New Mexico. Ph.D. dissertation
The University of Arizona, 407 pp.
Herweijer, C., and R. Seager, 2008: The global footprint of persistent extra-tropical drought in the instrumental era. Int.
J. Climatol., 28, 1761–1774.
Hughes, M. K., and L. J. Gramlich, 1996: Multi-millennial dendroclimatic studies from the western United States. Climatic
Variations and Forcing Mechanisms of the Last 2000 Years,
P. D. Jones et al., Eds., NATO ASI Series, Vol. 41, 109–124.
Huybers, P., and W. Curry, 2006: Links between annual, Milankovitch and continuum temperature variability. Nature,
441, 329–332.
Kantelhardt, J., E. Koscielny-Bunde, D. Rybski, P. Braun,
A. Bunde, and S. Havlin, 2006: Long-term persistence andmultifractality of precipitation and river runoff records.
J. Geophys. Res., 111, D01106, doi:10.1029/2005JD005881.
Kipfmueller, K. F., and M. W. Salzer, 2010: Linear trend and climate
response of five-needle pines in the western United States related to treeline proximity. Can. J. For. Res., 40, 134–142.
Klemes, V., 1974: Hurst phenomenon: A puzzle? Water Resour.
Res., 10, 675–688.
Laird, K. R., S. C. Fritz, E. C. Grimm, and P. G. Mueller, 1996:
Century-scale paleoclimatic reconstruction from Moon Lake:
A closed-basin lake in the northern Great Plains. Limnol.
Oceanogr., 41, 890–902.
——, ——, and B. F. Cumming, 1998: A diatom-based reconstruction of drought intensity, duration, and frequency
from Moon Lake, North Dakota: A sub-decadal record of thelast 2300 years. J. Paleolimnol., 19, 161–179.
Liang, X., D. P. Lettenmaier, E. F. Wood, and S. J. Burges, 1994: A
simple hydrologically based model of land-surface water andenergy fluxes for general-circulation models. J. Geophys. Res.,99 (D7), 14 415–14 428.
Lomb, N. R., 1975: Spectrographic study of beta-centauri. Mon.
Not. Roy. Astron. Soc., 172, 639–647.
Meko, D. M., 2001: Reconstructed Sacramento River system runoff from tree rings. California Department of Water Resources Rep., 282 pp.
——, M. D. Therrell, C. H. Baisan, and M. K. Hughes, 2001: Sacramento River flow reconstructed to AD 869 from tree rings.
J. Amer. Water Resour. Assoc., 37, 1029–1039.
——, C. A. Woodhouse, C. A. Baisan, T. Knight, J. J. Lukas, M. K.
Hughes, and M. W. Salzer, 2007: Medieval drought in theupper Colorado River basin. Geophys. Res. Lett., 34, L10705,
doi:10.1029/2007GL029988.
Mitchell, T. D., and P. D. Jones, 2005: An improved method of
constructing a database of monthly climate observations
and associated high-resolution grids. Int. J. Climatol., 25,
693–712.
Ni, F., T. Cavazos, M. Hughes, A. Comrie, and G. Funkhouser, 2002:
Cool-season precipitation in the southwestern USA since AD
1000: Comparison of linear and nonlinear techniques for reconstruction. Int. J. Climatol., 22, 1645–1662, doi:10.1002/joc.804.

5878

JOURNAL OF CLIMATE

Palmer, W. C., 1965: Meteorological drought. U.S. Department ofCommerce Research Paper 45, 65 pp.
Pelletier, J. D., 2008: Quantitative Modeling of Earth System Processes. Cambridge University Press, 304 pp.
——, and D. L. Turcotte, 1997: Long-range persistence in climatological and hydrological time series: Analysis, modeling and application to drought hazard assessment. J. Hydrol., 203, 198–208.
Pohl, K., and Coauthors, 2003: A cool season precipitation reconstruction for Saltillo, Mexico: 11th North American dendroecological fieldweek: Climatic reconstruction researchgroup: Saltillo, Mexico, August 2001. Tree-Ring Res., 59, 11–19.
Polyak, V. J., and Y. Asmerom, 2001: Late Holocene climate and culturalchanges in the southwestern United States. Science, 294, 148–151.
Rasmussen, J., V. Polyak, and Y. Asmerom, 2006: Evidence for
Pacific-modulated precipitation variability during the late
Holocene from the southwestern USA. Geophys. Res. Lett.,
33, L08701, doi:10.1029/2006GL025714.
Ropelewski, C. F., and M. S. Halpert, 1987: Global and regional
scale precipitation patterns associated with the El Ni~no/Southern Oscillation. Mon. Wea. Rev., 115, 1606–1626.
Routson, C. C., C. A. Woodhouse, and J. T. Overpeck, 2011: Second century megadrought in the Rio Grande headwaters,
Colorado: How unusual was medieval drought? Geophys. Res.
Lett., 38, L22703, doi:10.1029/2011GL050015.
Salzer, M. W., and K. F. Kipfmueller, 2005: Reconstructed temperature and precipitation on a millennial timescale from treerings in the Southern Colorado Plateau, USA. Climatic
Change, 70, 465–487.
——, M. K. Hughes, A. G. Bunn, and K. F. Kipfmueller, 2009:
Recent unprecedented tree-ring growth in bristlecone pine at
the highest elevations and possible causes. Proc. Natl. Acad.
Sci. USA, 106 (48), 20 348–20 353.
Scargle, J. D., 1982: Studies in astronomical time-series analysis.
II—Statistical aspects of spectral-analysis of unevenly spaced
data. Astrophys. J., 263, 835–853.
Schmidt, G. A., and Coauthors, 2012: Climate forcing reconstructions for use in PMIP simulations of the last millennium (v1.1).
Geosci. Model Dev., 5, 185–191, doi:10.5194/gmd-5-185-2012.

VOLUME 26

Schulz, M., and K. Stattegger, 1997: SPECTRUM: Spectral analysis
of unevenly spaced paleoclimatic time series. Comput. Geosci., 23, 929–945.
Seager, R., N. Harnik, W. Robinson, Y. Kushnir, M. Ting,
H. Huang, and J. Velez, 2005a: Mechanisms of ENSO-forcingof hemispherically symmetric precipitation variability. Quart.
J. Roy. Meteor. Soc., 131, 1501–1527, doi:10.1256/qj.04.96.
——, Y. Kushnir, C. Herweijer, N. Naik, and J. Velez, 2005b:
Modeling of tropical forcing of persistent droughts and pluvials over western North America: 1856–2000. J. Climate, 18,4065–4088.
Stahle, D. W., F. K. Fye, E. R. Cook, and R. D. Griffin, 2007: Treering reconstructed megadroughts over North America since
AD 1300. Climatic Change, 83, 133–149.
Thomson, D. J., 1982: Spectrum estimation and harmonic-analysis.
Proc. IEEE, 70, 1055–1096.
Truebe, S., T. Ault, and J. Cole, 2010: A forward model of cave
dripwater d18O and application to speleothem records. IOP
Conf. Ser., 9, 012022, doi:10.1088/1755-1315/9/1/012022.
Urban, F. E., J. E. Cole, and J. T. Overpeck, 2000: Influence ofmean climate change on climate variability from a 155-year
tropical Pacific coral record. Nature, 407, 989–993.
Vyushin, D. I., and P. J. Kushner, 2009: Power-law and longmemory characteristics of the atmospheric general circulation.
J. Climate, 22, 2890–2904.
——, ——, and J. Mayer, 2009: On the origins of temporal powerlaw behavior in the global atmospheric circulation. Geophys.
Res. Lett., 36, L14706, doi:10.1029/2009GL038771.
Woodhouse, C. A., 2000: Extending hydrologic records with treerings. Water Resour. Impact, 2, 25–27.
——, 2001: A tree-ring reconstruction of streamflow for the Colorado Front Range. J. Amer. Water Resour. Assoc., 37, 561–569.
——, and J. T. Overpeck, 1998: 2000 years of drought variability inthe central United States. Bull. Amer. Meteor. Soc., 79, 2693–
2714.
——, S. T. Gray, and D. M. Meko, 2006: Updated streamflow reconstructions for the upper Colorado River basin. Water Resour. Res., 42, W05415, doi:10.1029/2005WR004455.