This study examines the geographic structure of observed trends in key hydrologically relevant variables across the western United States at ⅛° spatial resolution during the period 1950–99. Geographical regions, latitude bands, and elevation classes where these trends are statistically significantly different from trends associated with natural climate variations are identified. Variables analyzed include late-winter and spring temperature, winter-total snowy days as a fraction of winter-total wet days, 1 April snow water equivalent (SWE) as a fraction of October–March (ONDJFM) precipitation total [precip(ONDJFM)], and seasonal [JFM] accumulated runoff as a fraction of water-year accumulated runoff. Observed changes were compared to natural internal climate variability simulated by an 850-yr control run of the finite volume version of the Community Climate System Model, version 3 (CCSM3-FV), statistically downscaled to a ⅛° grid using the method of constructed analogs. Both observed and downscaled model temperature and precipitation data were then used to drive the Variable Infiltration Capacity (VIC) hydrological model to obtain the hydrological variables analyzed in this study. Large trends (magnitudes found less than 5% of the time in the long control run) are common in the observations and occupy a substantial part (37%–42%) of the mountainous western United States. These trends are strongly related to the large-scale warming that appears over 89% of the domain. The strongest changes in the hydrologic variables, unlikely to be associated with natural variability alone, have occurred at medium elevations [750–2500 m for JFM runoff fractions and 500–3000 m for SWE/Precip(ONDJFM)] where warming has pushed temperatures from slightly below to slightly above freezing. Further analysis using the data on selected catchments indicates that hydroclimatic variables must have changed significantly (at 95% confidence level) over at least 45% of the total catchment area to achieve a detectable trend in measures accumulated to the catchment scale.
A growing number of studies have investigated recent trends in the observed (and simulated) hydrometeorological variables across the western United States. The main changes observed in this region include a large increase of winter and spring temperatures (Dettinger and Cayan 1995; Karoly et al. 2003; Bonfils et al. 2008a,b), a substantial decline in the volume of snowpack in low and middle altitudes (Lettenmaier and Gan 1990; Dettinger et al. 2004; Knowles and Cayan 2004; Hamlet et al. 2005), a significant decline in 1 April snow water equivalent (SWE; Mote 2003; Mote et al. 2005; Mote 2006; Mote et al. 2008; Pierce et al. 2008), and a reduction in March snow cover extent (Groisman et al. 2004). A reduction of the proportion of precipitation falling as snow instead of rain has also been observed (Knowles et al. 2006), as well as an earlier streamflow from snow-dominated basins (Dettinger and Cayan 1995; Cayan et al. 2001; Stewart et al. 2005; Regonda et al. 2005), and a sizeable increase of winter streamflow fraction (Dettinger and Cayan 1995; Stewart et al. 2005). These changes are likely to have important effects on western U.S. water resources management and distribution if they continue into future decades, as is projected for greenhouse-forced warming trends (Barnett et al. 2004; Christensen et al. 2004; Christensen and Lettenmaier 2007; Cayan et al. 2008a,b). This is because much of the water in the western United States is stored as snow in winter, which starts to melt during late spring and early summer. An earlier snowmelt and more precipitation falling as liquid instead of stored as snow could provide new stresses on the existing water resources management structures in the western United States in coming decades.
Some of these studies have indicated that such changes are partially linked with rising greenhouse gas (GHG) concentrations, which alter temperature and thus affect the snowpack distribution in the western United States, and partly from natural climatic decadal fluctuations over the North Pacific Ocean (Dettinger and Cayan 1995). Pacific decadal oscillation (PDO; Mantua et al. 1997) fluctuations, the dominant decadal natural variability in this region, however, can only partially explain the magnitude of the recent changes in snowfall fractions (Knowles et al. 2006), spring snowpack (Mote et al. 2005), and center timing from snow-dominated basins (Stewart et al. 2005). Knowles et al. (2006), Mote et al. (2005), and Stewart et al. (2005) argued that the remaining parts of the variability might be due to large-scale anthropogenic warming.
Only recently have formal efforts been undertaken (Knutson et al. 1999; Karoly et al. 2003; Maurer et al. 2007; Bonfils et al. 2008a) to distinguish whether the recent changes occurred because of internal natural variations of the climate system or because of human influence, using rigorous detection-and-attribution procedures (Hegerl et al. 1996, 1997; Barnett et al. 2001; Zwiers and Zhang 2003; The International Ad Hoc Detection and Attribution Group 2005; Zhang et al. 2007; Santer et al. 2007). Detection is the determination that a particular climate change or sequence is unlikely to have occurred solely because of natural causes. In the present study, climate from a long model control run is used to characterize the long-term variations that can arise solely from the internal fluctuations of the global climate system. Other external but natural forcings of the climate system, such as solar-irradiance changes and volcanic emissions, were not tested here. [Barnett et al. (2008) tested hydroclimatic trends from a simulation with climate forced only by historical solar and volcanic influences and found that observed trends could not be attributed to those influences.] Attribution (not undertaken here) is a later step in which the particular causes of the “unnatural” parts of observed trends are rigorously identified. Detection studies are important because if the recent changes are found to be due to internal natural variations alone, one can reasonably anticipate that the climate system will return to its past states after some time has passed.
Karoly et al. (2003) carried out a comparison of temperature trends in observations and three model simulations at the scale of North America. They found that the temperature changes from 1950 to 1999 were unlikely to be due to natural climate variation alone, whereas most of the observed warming from 1900 to 1949 was naturally driven. Accounting for uncertainties in the observational datasets, Bonfils et al. (2008a) observed increases in California-averaged annual mean temperature for the time periods 1915–2000 and 1950–99. They found these warmings are too large and too prolonged to have likely been caused by natural variations alone. In their study, natural variations were characterized using multiple control simulations (no change in greenhouse-gas concentrations) by multiple global climate models to develop multimodel 86- and 50-yr trend distributions. The authors also indicated that the recent warming in California is particularly fast in winter and spring and is likely associated with human-induced changes in the large-scale atmospheric circulation pattern occurring over the North Pacific Ocean. The hypothesis that human activities have influenced the circulation over the North Pacific Ocean is strengthened by a recent study (Meehl et al. 2009) that has identified an anthropogenic component in the phase shifts of the PDO mode.
More recently, a series of three formal fingerprint-based detection and attribution studies have been performed for the western U.S. region. The first study focused on hydrologically relevant temperature variables from late winter to early spring (Bonfils et al. 2008b). The second study examined SWE as a fraction of precipitation (SWE/P) over nine mountainous regions in the western United States (Pierce et al. 2008). The third study analyzed center timing of streamflow (CT; defined as the day when half of the water-year flow has passed a given point) in three major tributaries areas of the western United States (the California region was represented by the Sacramento and San Joaquin Rivers, Colorado River at the Lees Ferry, and Columbia River at The Dalles; Hidalgo et al. 2009). Bonfils et al. (2008b) showed that the changes in the observed temperature-based indices across the mountainous regions are unlikely, at a high statistical confidence, to have occurred because of natural variations. They concluded that changes in the climate as the result of anthropogenic GHG, ozone, and aerosols are causing part of the recent changes. Similarly, Pierce et al. (2008) and Hidalgo et al. (2009) showed that the observed changes in SWE/P and in CT are unlikely to have arisen exclusively from natural internal climate variability. Barnett et al. (2008) performed a multiple variable detection and attribution study and showed how the changes in minimum temperature (Tmin), SWE/P, and CT for the period 1950–99 covary. They concluded, with a high statistical significance, that up to 60% of the climatic trends in those variables are human related.
In regions with complex topography such as the western United States, there are strong gradients in temperature and associated hydrologic structure. These gradients motivate investigating responses to climate variability and climate change at high resolution (e.g., ∼12 km), scales that are much finer than are provided by global climate models. However, the detection of climate change at fine scales is challenging because less spatial averaging means “weather noise” increases with deceasing scale (Karoly and Wu 2005). On the other hand, when a variety of elevational settings are lumped together, the response to warming may be diluted because of the strong variations that are mixed together. For example, although Hidalgo et al. (2009) were able to detect changes in CT that were different from background natural variability at a high level of confidence in the Columbia basin, changes aggregated over the California Sierra Nevada and in the Colorado basins were only marginally significant or not at all. Maurer et al. (2007) examined whether the decreases in CT at four river points in the Sierra Nevada are statistically significantly different from changes associated with internal natural variability, and they concluded that the recent observed trends are still within simulated natural variations. This suggests that, in settings that contain strong topographic variation, climate responses can usefully be evaluated at finer, rather than coarser, spatial units despite the increase in weather noise and the amount of uncertainty in the forcing data and modeling.
The present study investigates the hypothesis that there are detectable climate changes that can be delineated over a complex topographic setting using a high-resolution ⅛° (∼12 km) spatial network over the western United States. (Fig. 1a). Because of the increased noise-to-signal issues that plague evaluations at this scale, we do not attempt to formally attribute the causes of the unnatural trends at every grid cell. Rather, we use fine-resolution simulations to investigate the spatial structure of detectable trends across the snow-dominated western United States. Our objective is to find the fraction of the regions of the western United States where we should expect to see detectably unnatural trends. We focus on several indices that are hydrologically relevant in the area of interest, including late-winter and spring temperature, total number of wintertime snowy days as a fraction of all wet winter days, 1 April SWE as a fraction of October–March (ONDJFM) precipitation total [SWE/Precip(ONDJFM)], and January–March runoff as fraction of water-year total. We also extend the analysis to consider how the detectability of trends over a whole catchment depends upon the fraction of individual grid cells within the catchment that exhibit detectable changes. This dependence provides useful rules of thumb for use in designing monitoring networks or helping to decide whether detectable trends in a catchment of interest should even be expected.
Following the Intergovernmental Panel on Climate Change (IPCC) Fourth Assessment Report, there is a nuanced definition of a two step “joint attribution” that begins by distinguishing between detecting a change relative to the variability in an observed record and a change that exceeds natural variability (Rosenzweig et al. 2007). The latter conclusion implies that some external forcing must be at work, although not attributing it to specific causes, and is the first step in joint attribution. The next step is to demonstrate that the changes are best (and only) explained by a specific forcing. In this study we carried out the first step of the joint attribution by showing that hydrological measures in the period 1950–99 trended within the historical variability and then determining whether those trends were consistent with natural variability as represented by a long control run of the combination of a climate model and a hydrologic model.
Section 2 presents the datasets and models used in our study. A description of the methodology and definitions of various climate indices analyzed in this study are given in section 3. Section 4 presents results we obtained for the different indices analyzed. The relationship between total significant area and detectability at the catchment scale is also presented in section 4. A summary and conclusions are given in section 5.
2. Datasets and models
a. Observed data and global climate model results
Gridded meteorological observations were used to characterize observed climate changes across the western United States during the period 1950–99. Daily P, Tmin, maximum temperature (Tmax), and wind speed at ⅛° spatial resolution were obtained from the Surface Water Modeling Group at the University of Washington (available online at http://www.hydro.washington.edu; Hamlet and Lettenmaier 2005). We also repeated the analyses using a different meteorology, the Maurer et al. (2002) dataset, to determine the extent to which our results depended on the inputs. The Maurer et al. (2002) and Hamlet and Lettenmaier (2005) datasets are closely related in that both used the same station data, although the Hamlet and Lettenmaier (2005) dataset focused much more on the Historical Climatology Network (Easterling et al. 1996) subset of stations, which are chosen and corrected to eliminate most temporal inhomogeneities, whereas the Maurer et al. (2002) dataset treated all stations more or less equally. The datasets also share a reliance on monthly Precipitation-elevation Regressions on Independent Slopes Model (PRISM) data fields (Daly et al. 1994) to adjust for elevation effects on precipitation and temperature. Such corrections are necessary because topography in the study region strongly determines not only spatial patterns of precipitation (and temperature) but also basin scale–to–regional scale totals of precipitation. In particular, Pan et al. (2003) found, in comparisons of North American Land Data Assimilation System (NLDAS; Mitchell et al. 1999) fields with snowpack telemetry (SNOTEL) data, that snow accumulation was underestimated by as much as half when no PRISM-like adjustments were included. We have used wind speed data from the Surface Water Modeling Group at the University of Washington [which Maurer et al. (2002) obtained from the National Centers for Environmental Prediction/National Center for Atmospheric Research (NCEP–NCAR) reanalysis (Kalnay et al. 1996)], which are obviously very low resolution (at T62 Gaussian grid, approximately 1.9°) and are likely erroneous for complex topography and might have an effect on the results (e.g., sublimation). The results that follow did not depend sensitively on the choice between these two admittedly closely related meteorological datasets. In the following sections, only the results using the Hamlet and Lettenmaier (2005) dataset are presented, because this dataset was produced with attention to accounting for station and instrument changes that would otherwise add nonclimatic noise to the long-term trend signals (Hamlet and Lettenmaier 2005).
Internal climate variability in the western United States in the absence of any anthropogenic effects is characterized using precipitation and temperature data from an 850-yr preindustrial control simulation of the NCAR/Department of Energy (DOE) Community Climate System Model, version 3 (CCSM3; Collins et al. 2006). The simulation was performed at Lawrence Livermore National Laboratory and used the finite volume (FV) dynamical methods for the atmospheric transport (CCSM3-FV; Bala et al. 2008a,b). The horizontal spatial resolution of the atmospheric model was 1° × 1.25° with 26 vertical levels. This preindustrial control simulation used constant 1870-level atmospheric composition to force the model. Bala et al. (2008a) have evaluated the fidelity of a 400-yr present-day control climate simulation that used this FV configuration for CCSM3. They found significant improvement in the simulation of surface wind stress, sea surface temperature, and sea ice when compared to a spectral version of CCSM3.
b. Downscaling of the control run
Daily P total, daily Tmax, and daily Tmin from the CCSM3-FV model were downscaled to ⅛° resolution using the constructed analogs (CANA; Hidalgo et al. 2008) statistical downscaling method. The CANA procedure starts with a simple variance correction to ensure the same variability of the GCM data as observations. Then, the bias-corrected global model fields are downscaled using a linear combination of previously observed patterns1 (Maurer and Hidalgo 2008; Hidalgo et al. 2008). The 30 most similar previously observed patterns are used in a linear regression to obtain an estimate that best matches, on the coarse grid, the GCM pattern to be downscaled. The downscaled values of precipitation and temperatures are obtained by applying the linear regression coefficients to the finescale versions of the previously observed patterns. Results using CANA and those obtained with another statistical downscaling methodology (bias correction and spatial downscaling; Wood et al. 2004) are qualitatively similar (Maurer and Hidalgo 2008). An advantage of the CANA method over the bias correction and spatial downscaling method is that CANA can capture changes in the diurnal cycle of temperatures; the disadvantage is that to do this it requires daily rather than monthly data. Details of the CANA method can be found in Hidalgo et al. (2008).
c. Hydrological model
Runoff and SWE, major variables of interest to hydrological studies, have not been readily observed at the temporal and spatial scales required for this study. Likewise, they cannot be obtained by downscaling global model results because no library of observed fine-resolution daily fields exists to be used in the downscaling scheme. Accordingly, to produce both the “observed” and climate model–driven SWE and runoff fields on the fine spatial scale, we use the Variable Infiltration Capacity (VIC; Liang et al. 1994, 1996) model (version 4.0.5 beta release 1). To estimate the observed trends, we drove VIC with observed daily P, Tmin, Tmax, and wind speed fields on the ⅛° grid; to estimate the downscaled climate model trends, we drove VIC with the downscaled model daily P, Tmin, and Tmax, along with climatological wind speed fields, on the ⅛° grid. Derived variables such as radiation, humidity and pressure are estimated within the model based on the input P, Tmax, and Tmin values using the algorithms of Kimball et al. (1997) and Thornton and Running (1999). How well the algorithms used to estimate these variables will apply in the future is uncertain and could not be addressed here. VIC uses a tiled representation of the land surface within each model grid cell and allows subgrid variability in topography, infiltration, and land surface vegetation classes (Maurer et al. 2002). The subsurfaces are modeled using three soil layers with different thicknesses. Surface runoff uses an infiltration formulation based on the Xinanjiang model (Wood et al. 1992), whereas baseflow follows the Arno model (Liang et al. 1994). Subgrid variability in soil moisture storage capacity is represented through the use of a spatial probability distribution function, and a nonlinear function is used to model the baseflow component from the lowest soil layer (Liang et al. 1994; Sheffield et al. 2004). VIC has been successfully applied at spatial scales ranging from regional to global (Hamlet and Lettenmaier 1999; Nijssen et al. 2001; Maurer et al. 2002; Christensen et al. 2004; Wood et al. 2004; Christensen and Lettenmaier 2007; Hamlet et al. 2007; Maurer et al. 2007; Sheffield and Wood 2007; Barnett et al. 2008; Pierce et al. 2008; Hidalgo et al. 2009).
The calibrated soil parameters for VIC were obtained from Andrew W. Wood at the University of Washington, presently at 3TIER in Seattle. The vegetation cover was obtained from the NLDAS and was held static through all the simulations. In particular, leaf area index (LAI) values were specified from average values during the period 1981–94 of Myneni et al.’s (1997) monthly global LAI database. Realistically, the land cover probably changed in many areas, but these changes were not explored here. In this study we did not include a frozen soils component because of the large computational costs that would have been required in the very long control simulations made here. The VIC model was run at a daily time step, with a 1-h snow model time step and five snow elevation bands. The first nine months of the simulations were used for model initializations and not considered for further analysis, as suggested by Hamlet et al. (2007). The VIC model uses the gridded observed and model control run meteorologies, along with the physiographic characteristics of the catchment (e.g., soil and vegetation), to calculate runoff, baseflow, soil moisture at three soil layers, and SWE. The ability of the model to simulate monthly streamflow at some of the calibration points across the study domain is satisfactory when compared with the naturalized streamflow (Maurer et al. 2002; Hamlet et al. 2007; see Fig. 3 of Hidalgo et al. 2009). Additionally, Mote et al. (2005) found reasonable agreement between the spatial pattern of observed SWE and the VIC simulated values.
d. Definition of climate variables
Our study focused on five hydrologically relevant detection variables:
Monthly and seasonal precipitation as a fraction of total precipitation over the water year (October through September)
Monthly and seasonally averaged temperatures
Seasonal (January–March) accumulated runoff (as simulated by VIC), calculated as the fraction of accumulated runoff over the water year
1 April SWE as a fraction of October through March precipitation total (SWE/Precip(ONDJFM)), chosen to reduce the influence of precipitation on snowpack and produce a snow-based climate index that is more directly sensitive to temperature changes (Pierce et al. 2008)
- The number of winter days with precipitation occurring as snow divided by the total number of winter days with precipitation. A given wet day (precipitation >0.1 mm), in the period November through March, was classified as a snowy day if the amount of snowfall (S) was greater than 0.1-mm water equivalent. Here, S was calculated using the same equation as VIC:
e. Natural variability in the control run
The strength of the conclusions of any detection analysis rely on the ability of the control model to represent the strength and key features of the natural internal climate variability in the absence of anthropogenic effects. In particular, the ability to simulate decadal variability is crucial for the identification of slow-evolving climate responses to slow-evolving external forcings. To compare the low-frequency variability in the model control run simulation to observations, we computed standard deviations in each grid cell for each index after the application of a 5-yr low-pass filter. The observations were linearly detrended before the calculation in an attempt to remove the linear part of possible anthropogenic influence. The low-frequency variability in the control simulation is reasonably well represented with no evidence that the model systematically underestimates or overestimates the observed variability for all climate indices (Fig. 2). Thus, we conclude that the CCSM3-FV model used here provides an adequate representation of natural internal climate variability for our detection work. Barnett et al. (2008), Pierce et al. (2008), and Bonfils et al. (2008b) have also addressed this issue using the CCSM3-FV data (i.e., Barnett et al. 2008, their Fig. S3) and reached similar conclusions.
At each grid cell and for each variable, the linear trend over 50-yr segments (with the start of each segment offset by 10 yr from the previous segment’s start) was calculated from the 850-yr control run. This produced 80 partially overlapping estimates of what the 50-yr trend could be in the absence of anthropogenic forcing. An Anderson–Darling test2 (Anderson and Darling 1952) showed that the distribution of control run trends was Gaussian in the great majority of the grid cells, except for only some grid cells of the JFM runoff fractions. Specifically, the percentage of grid cells in which the distributions were non-Gaussian was about 7% for JFM average temperature, snowy days as a fraction of wet days, and SWE/Precip(ONDJFM). For JFM runoff fractions, the percentage of grid cells that were non-Gaussian was as large as about 20% of the region of interest. Nonetheless, we included the non-Gaussian cells in the subsequent analysis because we wanted to maintain the same spatial domain for all of the variables analyzed. In addition, to maintain consistency among analyses, we used the mean and standard deviation from the control run to fit a Gaussian distribution at all grid cells.
We evaluated the observed trends mainly over water years 1950–99 and later over different starting and ending years within this period. The probability of finding the observed trend in the estimated Gaussian distribution of unforced trends is computed using a two-tailed test. We used a two-tailed test because we did not make any a priori assumption on the direction of the trends of the indices analyzed because we wanted to evaluate, for example, a significant lack of negative temperature trends as well as a significant surplus of positive temperature trends. Figure 3 shows the schematic diagram of the methodology we employed to compute the probability. The bars represent the distribution of the 50-yr unforced trends in the model control run. If an observed trend (short vertical line) falls within the shaded region (showing the two-tailed p = 0.05 level), which indicates the amplitude of naturally driven trends that occur only 5% of the time, then we can conclude that this trend is unlikely to be the result of internal natural variations. Probability maps for each variable were obtained by applying this procedure to all grid cells across the western United States.
We also examined the effect of spatial coherence on our results using a Monte Carlo simulation as in Livezey and Chen (1983) and Karoly and Wu (2005). Because there is a high spatial coherence of the hydrometeorological variables, this can lead to spurious detection, as described in those references. The Monte Carlo approach we used accounts for the effects of this spatial coherence: We analyzed all 800 possible 50-yr segments (i.e., moving 50-yr windows with 1-yr shifts) from the 850-yr control run to compute probabilities (based on 50-yr means and standard deviations) at each grid cell for each of the hydrometeorological variables. This resulted in 800 probability maps for each variable. The fraction of grid cells exhibiting apparent trends that rose above the higher-frequency natural variations in each 50-yr segment (at 95% confidence level) was counted in each probability map, giving us 800 values with which to estimate the distribution of the fractions of grid cells that might, by chance or natural variability, appear to yield seemingly detectable trends in a 50-yr segment. Although this number would be 5% on average over the 850-yr control run if all grid cells varied independently of each other, the lack of independence between nearby grid cells means that, in any particular 50-yr segment, considerably more or less grid cells can show seemingly significant trends. Consequently, the 95th percentile of the distribution of possible “trends” arising from natural variations is considerably broader than 5% of the grid cells, as will be shown in section 4b. The number of grid cells in the observational data with trends was then compared to the distribution of trending grid cells from the control run to decide whether observed trends can be explained by natural variability.
Because our main focus is to investigate the changes in hydrology, we begin by focusing our analysis on the mountainous western United States, where warming-related effects are particularly important (Mote et al. 2005) and for which hydrological changes may have large implications for the water supply, ecology, or likelihood of wildfire in the region. As in Hamlet et al. (2007), we include locations where long term (1950–99) mean 1 April SWE simulated by VIC is greater than 50 mm.
In the last section, we extended the analysis using the data on the selected catchments across the western United States to identify relationships, for each of the climate variables, between the fraction of catchment area within which significant changes have occurred and the significance of detectability at the whole-catchment scale. Such information can be of practical use to resource managers trying to understand local climate changes. Trends in 66 catchments across the western United States were analyzed (Fig. 1a). The areas of the catchments range between 720 and 679 250 km2, with a median value of 19 000 km2. The average elevations of the catchments range between 360 and 2900 m, with a median value of 1765 m. The catchment-average spring [March–May (MAM)] temperatures range between −2° and 14°C, with a median value of 3°C.
4. Results and discussion
a. Spatial pattern of observed trends
We analyzed observed monthly precipitation (for January through March) as a fraction of water-year total precipitation, and monthly average temperatures, for the period 1950–99. We found trends in monthly precipitation fraction that were well within the distribution of natural variability as estimated from the control model run (not shown). This agrees with the results of Barnett et al. (2008), who also found that natural variability could account for changes in water-year total precipitation for the mountainous western United States during this period.
Observations show warming temperatures since 1950 over the western United States during the months of January–March (Fig. 4a). Among these months, March average temperature shows the strongest and most widespread upward trends, with larger warming in the interior west than along the coast. Notable warming in January is concentrated along the coast of California region and the Columbia River basin, and February average temperature shows widespread but only mild warming trends; see Knowles et al. (2006) for more detail on these patterns.
In view of the considerable warming trends for the study domain during January and March, we investigated changes in observed JFM average temperature. A linear trend calculation using the JFM average temperature shows a considerable upward trend across most parts of the snow-dominated western United States, with notably larger warming trends across the high mountains of the Columbia River basin (Fig. 5a).
A chain of hydrologic responses to warming is evident in the trends. Reductions in observed winter-total snowy days as a fraction of winter-total days with precipitation (indicating a decrease in days with snowfall) are also common across many parts of the snow-dominated region in the observations, except in regions in the northern Rockies that show no trend (Fig. 5b). There are widespread downward trends in observed SWE/Precip(ONDJFM) across most parts of the snow-dominated western United States, with stronger downward trends in the northern Rockies of the Columbia River basin along with some upward trends in the southern Sierra Nevada and part of the northern Rockies (Fig. 5c). These findings are in agreement with those of Pierce et al. (2008), who described declining fractional SWE/P from snow course data across the nine mountainous regions of the western United States. These trend patterns are also consistent with the results in Mote et al. (2005), who analyzed 1 April SWE from 824 snow stations for the period 1950–97, and in Hamlet et al. (2005), who analyzed VIC-simulated 1 April SWE. Using regression analyses, those two studies attributed the widespread downward trend in SWE to a warming trend and a more regional upward trend in SWE in the southern Sierra Nevada (in the California region) to an increase of precipitation over the period. Changes in snowmelt initiation and changes in the snow-to-rain ratio should concur with large changes in runoff. Indeed, upward trends in JFM runoff fractions predominate across the snow-dominated western United States, except for some weaker trends in the Canadian part of the Columbia River basin and Colorado Rockies and some weaker downward trends in the southern Sierra Nevada (Fig. 5d).
b. Comparison of observed trends with model control run trends distribution at the grid scale
Figures 4b and 6 illustrate the probability of the observed trends in Figs. 4a and 5 arising in absence of any external forcings. There are considerable regions over which the observed trends in January and March average temperature are unlikely to have arisen from internal natural variability alone, at 95% significance level (Fig. 4b). By contrast, the mild warming trends in February are not detectably different from internal natural variability (Fig. 4b).
The observed trends toward warmer JFM average temperature across nearly all (89%) of the snow-dominated regions of the western United States cannot be explained (at 95% confidence level) by internal natural variability alone, except for relatively small areas of the southern Sierra (California region) and southern Rockies (lower Colorado River basin; Fig. 6a). The downward trends of the snow day fraction of wet days (Fig. 6b) also exhibit detectable signals for 42% of the grid cells over the mountainous western United States. The decline in SWE/Precip(ONDJFM) found in the observations, at 40% of the snow-dominated grid points, is also unlikely to be associated with natural variations alone in many regions (Fig. 6c). However, opposite changes in regions containing upward trends in SWE/Precip(ONDJFM) (e.g., southern Sierra Nevada and Utah) cannot be confidently distinguished from natural internal variability. Consistent with the warming and reduction in fraction of snowy days and SWE/Precip(ONDJFM), increases in JFM runoff fraction exceed those expected from natural variations alone over broad mountainous regions, at some 37% of the snow-dominated grid points, especially in the Columbia River basin (Fig. 6d). Changes in regions such as the southern Sierra Nevada (California region) and southern Rockies (Colorado River basin) cannot be distinguished confidently from natural variability. Notably, when we analyzed water-year runoff totals, historical trends were well within the distribution of natural variability as estimated from the control model run (not shown), so that the changes in JFM runoff fractions reflect seasonal timing shifts rather than overall increases in runoff. We have presented the results for JFM accumulated runoff as a fraction of water-year runoff here because JFM is the season in which there is the strongest temperature trend (Dettinger and Cayan 1995; Bonfils et al. 2008b). We also analyzed April–June (AMJ) runoff as a fraction of water-year runoff because AMJ runoff is important to water resources in the western United States. Widespread reductions of AMJ runoff fractions were found (not shown). However, there is more fraction-runoff variability in the spring season, and the changes found are not yet distinguishable from natural variability for the most part, except in scattered locales.
There is high spatial coherence in the meteorological and hydrological variables, which may overstate how widespread the statistically significant trends are (Livezey and Chen 1983) in Fig. 6. To estimate the sampling distribution of the percentage of the grid cells that could simultaneously show a statistically significant trend in the model control run, taking the observed spatial coherence into account, we have performed a Monte Carlo experiment based on resampling from the model control run, as described previously. The Monte Carlo–derived value is noted for each hydroclimatic variable in parenthesis in the panel titles of Fig. 6. The 95th percentile limits are still much less than the observed fractions of grid cells exhibiting significant trends for each variable, indicating that more grid cells contain significant trends than would be expected by chance, even taking the spatial coherence into account (Fig. 6).
An important property of the changes depicted in Fig. 6 is that they depend on elevation. To illustrate the dependence of the changes on elevation, we computed the total number of observed grid cells showing significant trends for each elevation class. Results are shown in Fig. 7, where the gray regions indicate results not significantly different from the control run at the 95% level, based the Monte Carlo resampling. The gray regions include zero; the wideness of the sampling distribution, noted above, means that even finding no grid cells with a significant trend does not indicate a statistically significant lack of trends. For example, finding no grid points at all with a statistically significant decrease in temperature is still consistent with the control run. Consequently, all significant results presented here arise from a surfeit of trends, not a deficit of trends.
In Fig. 7, the solid triangles in the left-hand panel show the numbers of positive trends, and the solid squares in the right-hand panel show the number of negative trends. The JFM warming (Fig. 7b, left) is detectable at all elevations, but the very small number of downward trends is not inconsistent with natural variability (Fig. 7b, right). The fraction of cells exhibiting significant upward trends decreases monotonically with elevation. The decline of the snowy days as a fraction of wet days from elevations near sea level up to 3000 m also exhibits a high tendency of being statistically significantly different from the distribution of trends from natural variations alone (Fig. 7c, right). Conversely, the grid cells with increasing trends—which show up mostly in small patches in the Rocky Mountains (e.g., see also Knowles et al. 2006)—are not inconsistent with natural variability (Fig. 7c, left). The reduction in SWE/Precip(ONDJFM) is particularly detectable at the lower elevations, but it is also detectable at medium altitudes (below 3000 m; Fig. 7d, right). The grid cells with positive trends (Fig. 7d, left) for all elevation classes, and the highest grid cells with negative trends (higher than 3000 m), exhibit trends in numbers that could be expected because of natural variability. The upward trends in the JFM runoff fractions in the regions with elevation ranging between approximately 750 and 2500 m tend to be statistically significantly more common than the model-estimated natural trends (Fig. 7e, right); however, the downward trends for all elevation classes and the upward trends at lower altitudes (below 750 m), and higher altitudes (above 2750 m) are not statistically significant in numbers than those that would occur because of natural variability (Fig. 7e).
To parse the geographical distribution of trends still further, we divided the study domain into latitudinal bands 6.25° wide. We performed the same analysis as shown in Fig. 7, but for each latitudinal band. In general, the majority of significant trends are north of 36°, except for JFM average temperature and snowy days as a fraction of wet days (Fig. 8), which are present across the entire latitudinal range. Significant trends in SWE/Precip(ONDJFM) are only found north of 36°, and in that range, are present mostly at low and medium altitudes. Most of the significant JFM runoff fraction trends were found north of 42° and at altitudes between 500 and 2750 m (Fig. 8).
Figure 9 demonstrates another aspect of the hydrological changes—the number of grid cells that show significant trends, stratified by 1950–99 climatological spring average temperature classes (instead of elevation classes). In general, the results in terms of temperature should be opposite of the results in terms of elevation (as shown in Fig. 7) because temperatures decrease with altitude at the resolutions considered here. Nonetheless, we evaluated the trends in terms of temperature because temperatures also generally decrease with increasing latitude, so that neither altitudes nor latitudes alone could describe the complete relationships. Significant trends were nearly all found in locations having mean temperatures above −4°C. Interestingly, the changes for snowy days, SWE/Precip(ONDJFM), and runoff fractions are consistent with natural variability for cells where spring temperatures are below −4°C. The results support the findings of Knowles et al. (2006) that showed that regions at low-to-medium elevations with temperature near freezing are more likely to have a decrease in the fraction of precipitation falling as snow, and are also consistent with Mote et al. (2005), who found these elevations to have incurred unusual reductions in spring snowpack. JFM runoff fraction has trended most significantly at middle elevations—high enough to have significant snowmelt contributions but low enough so that temperatures are close to freezing during critical times. As noted above, decreasing trends in temperature and runoff, as well as the rare increasing trends in snowy days and SWE/Precip (ONDJFM), cannot be shown to be different from natural variability with this dataset. Also, we did not find precipitation trends to be distinguishable from natural variability, except around 1500-m elevation and winter temperatures around −4°C (Figs. 7a and 9a).
Thus hydrological trends driven by temperatures are the trends most distinguishable from natural variability. Figures 7 and 9 also show that changes in the sense a priori expected from warming conditions (e.g., a decrease of days with snowfall) are more prevalent than those in the opposite sense. Again, the changes in the JFM precipitation fraction at different temperature ranges are not beyond what could be expected because of internal natural variability, except in temperature class −4°C. Previous detection and attribution studies of regionally averaged variables (Barnett et al. 2008; Bonfils et al. 2008b; Pierce et al. 2008; Hidalgo et al. 2009) have successfully attributed the temperature trends and associated hydrologic responses, that we detect here at finescales, to forcing from greenhouse gases.
We also investigated the sensitivity of these results to the period analyzed. For example, the results from the JFM average temperature are shown in Fig. 10a. In this experiment we used three different analysis periods, all starting in 1950, to compute the observed trends: 30 yr (1950–79), 40 yr (1950–89) and 50 yr (1950–99). The results show that the longer periods contain more grid cells exhibiting a detectable warming trend (Fig. 10a, left). This is different from what is expected for natural variability in an equilibrated climate system, where the period of averaging will make no systematic difference to the fraction of grid cells deemed to have significant trends. Interestingly, the grid cells located at higher elevations (above approximately 1500 m) exhibit more detectable trends as the period increases in length. Also, the changes at the grid cells located at high elevations are not inconsistent with natural variability for the shorter period (1959–79; Fig. 10a, left). Two potential reasons can explain these results: increases in noise when trends are calculated over shorter periods, or the strength of the trend becomes stronger at the end of the period (which can occur if the climate responds to the slow-evolving anthropogenic forcing).
To investigate these possibilities, we reanalyzed the trends using a fixed period length of 30 yr but with three different starting years: 1950, 1960, and 1970 (Fig. 10, right). Starting in 1950, cells with warming that is greater than would be expected locally from the natural variability are all below 750-m elevation. Starting in 1960, grid cells with locally detectable warming are above 2250-m elevation, but the Monte Carlo resampling suggests that the numbers of trends seemingly distinguishable from natural variability are not, yet, any larger than might be expected from the spatially coherent natural-variability fields. Starting in 1970, though, cells above 2250-m elevation experienced a detectable warming (Fig. 10a, right). Thus the warming trends appear to have begun at lower elevations earlier than at higher elevations. Longer observational records also contributed to our growing ability to detect the long-term trends. Similar patterns were also found in the hydrological variables analyzed in this paper [SWE/Precip(ONDJFM) and JFM runoff fractions; Figs. 10b and 10c], indicating the crucial role of the length of the time series in analyses such as this.
c. Detection at catchment scale
In the real world, the hydroclimatic trends evaluated here are more likely to be addressed or observed at watershed-to-catchment scales than on the 12-km grid used here. For example, runoff is measured and managed primarily as streamflow accumulated to the catchment scale rather than as a distributed runoff pattern. In light of the strong elevation dependence of the detectability of trends discussed above, it is natural to ask, How much of a catchment must lie within the critical elevation bands and yield runoff with detectable trends before the observations from the catchment as a whole are likely to show detectable trends? To address this question and in an attempt to develop some rules of thumb for where to expect detectability of unnatural trends thus far, in this section we analyze the relations between fractions of catchment areas with detectable trends and corresponding detectability of trends at the whole-catchment scale.
Trends in 66 catchments across the simulation domain of the western United States were analyzed (Fig. 1a). Hydroclimatic variables from all grid cells within a given catchment were averaged for the observed (or simulated using the observed meteorology) and control run data. The probabilities of any resulting trends of the catchment-averaged observed time series were then computed using the same procedure previously applied at the gridcell scale (described in section 3). The detectability of unnatural trends within each catchment-averaged series was then compared to the fractions of grid cells within that catchment that were locally detectably distinguishable from the control-run natural variability.
This analysis indicates that approximately 25% of the catchment area must have trended significantly (at 95% confidence level) before there are detectable changes (at 95% confidence level) in the catchment level for snowy days as a fraction of wet days and SWE/Precip (ONDJFM). Approximately 45% of the catchment area must have trended significantly before there are detectable trends in JFM runoff fractions at the catchment scale (Fig. 11). We believe that the reason for the higher threshold for catchment-scale detectability for runoff is that runoff fractions are noisier and that the most significant runoff trends are more restricted in space (at least as reflected by elevation bands, Fig. 7), so that catchment-scale significance may be challenged from both the higher and lower parts of catchment (unlike the other hydrological variables considered here).
Because we have found that certain elevation zones or average spring temperature bands are most likely to yield detectable trends (thus far), it would be useful to know whether the (known) fraction of a catchment area within these ranges dictates detectability at the catchment scale better than the area with locally “detectable” trends, which generally is not known a priori. Unfortunately, no clearly preferred mean spring temperature ranges or elevation ranges that characterize the significant catchment were found, except with respect to JFM runoff fractions. Catchments with significant trends in JFM runoff fractions all have catchment-average spring temperatures between −2° and 6°C, and those catchments are located at the medium-elevation range (ranging approximately between 1400 and 2500 m). Fractions of catchment areas within such ranges, rather than catchment-average values, did not relate usefully to whole-catchment detectability.
5. Summary and conclusions
This study has used a finescale (⅛° latitude × ⅛° longitude) analysis of meteorological and hydrological variables to investigate the structure of observed trends from 1950–99 in some key hydrologically relevant measures across the western United States. Combined with estimates of natural variability from an 850-yr GCM control simulation, observations were evaluated to determine which elevations and locations have experienced trends that are unlikely to be derived entirely from internal natural climatic variations. The VIC hydrologic model was used to simulate the surface hydrological variables, both during the observational period (when driven by observed meteorology) and from the global climate models (when driven by downscaled model fields). Using key hydrologic measures—including JFM temperature, fraction of days with snow, SWE/Precip(ONDJFM), and JFM runoff fractions—we find that the observed winter temperature and each of the hydrologic measures have undergone significant trends over considerable parts (37%–89%) of the snow-dominated western United States (Fig. 6). These trends are not likely to have resulted from natural variability alone, as gauged from the distribution of trends produced from the long control simulation. In a relatively large portion of the Columbia River and to a lesser extend in the California Sierra Nevada and in the Colorado River basin, trends in snow accumulation and runoff timing across many middle altitudes are unlikely to have been caused by natural variations alone (Fig. 7). These trends are caused by the warming of regions with mean spring temperature close to freezing. The majority of the significant trends for SWE/Precip(ONDJFM) and JFM runoff as a fraction of water-year runoff occur north of 36°.
In all cases, the significant changes occurred in a direction consistent with the sign of the changes associated with warming, for example, JFM average temperature increases, days with snowfall decreases, snowpack decreases, and JFM runoff increases. Reinforcing this result is that trends that occurred in the opposite direction are no more frequent than would be expected from spatially coherent natural variability.
For SWE and JFM runoff fractions that we have evaluated here, good observational datasets do not exist for the spatial scales we considered. We have used the VIC hydrologic model forced by observed meteorological conditions to simulate these variables, a limitation of this study that should be kept in mind. Though the VIC model performance has been evaluated for the domain of interest for a number of variables (Maurer et al. 2002; Mote et al. 2005), there could be uncertainties arising from several factors, including the lack of ability to simulate accurate observed trend or uncertainties in the preparation of the gridded forcing dataset (particularly in the mountains as a result of fewer stations available for the interpolation); that is, there may be some biases due to specific stations used to construct the gridded meteorological dataset.
Experiments that considered different start and end points of the 1950–99 interval suggest that significant warming and associated hydrological trends, not explained by natural variations, have begun earlier at lower elevations than at higher elevations. Longer observational records contribute a growing ability to detect the trends.
We also analyzed the finescale data in snow-influenced catchments across the western United States. To find a detectable trend (at 95% confidence level) at the catchment scale, at least 25% of the total catchment area must have trended significantly for snowy days as a fraction of wet days and SWE/Precip(ONDJFM) but at least 45% of the area for JFM runoff fractions (Fig. 11). These thresholds provide a context to understand the behavior observed in the major tributaries of the western United States, Columbia at The Dalles, Colorado at the Lees Ferry and the California Sierra Nevada (used in Barnett et al. 2008; Hidalgo et al. 2009), as well as many smaller river basins. Among the three major tributaries analyzed there, the Columbia contains the largest percentage area with significant decreasing trends for 1 April SWE/Precip(ONDJFM) and for the increasing fraction of annual runoff in JFM, as shown in Table 1. Although the portion of the Sierra Nevada and Colorado with significant trends in these measures is 15%, or less, those in the Columbia exceed 25%. Stronger signatures observed in the Columbia basin are quite clearly a reflection of the greater proportion of low-to-middle elevations and, in association, a preponderance of late-winter and early-spring temperatures in the sensitive −2°C to +4°C category. Lower-to-middle altitudes (near sea level to nearly 3000 m) of California showed the second highest percentage area exhibiting significant trends, but these signals are diluted by the much larger number of grid cells that are located in an elevational environment where warming has not been great enough to produce a significant effect. Warming of a few degrees in the higher altitudes, above 3000 m, where the temperature is currently much below the freezing point in winter is not sufficient yet to produce detectable changes.
In addition to conducting climate detection on a very fine scale the present study differs from most previous trend significance studies, in which a more traditional significance test (parametric or nonparametric) is performed to assess whether or not an observed trend is significantly different from zero. Naturally occurring climate phenomena, such as the Pacific decadal oscillation, can give statistically significant trends over long periods, so the presence of nonzero trends is not necessarily inconsistent with the hypothesis that the trends are caused by natural variability. Instead we used long model control simulations to quantify the trends in our variables likely to arise from natural internal climate variability and compared the observed trends to those.
The present study describes hydrologic changes over the last five decades (1950–99), at fine scales (∼12 km) over the western United States, that rise above the level expected if these changes resulted solely from natural variability. Although this study establishes the detectability of these changes, we did not conduct experiments to attribute these changes to particular external forcings; that is, we have performed the first step of the “joint attribution” as outlined in the IPCC Fourth Assessment Report (Rosenzweig et al. 2007). However, given the conclusions of Barnett et al. (2008), Bonfils et al. (2008b), Pierce et al. (2008), and Hidalgo et al. (2009) using the same domain but at a much larger spatial scale (nine regions over the western United States), we can reasonably predict that the origin of a substantial portion of the trends is anthropogenic warming. If this warming continues into future decades as projected by climate models, there will be serious implications for the hydrological cycle and water supplies of the western United States. The present results usefully bring the results of regional-scale detection and attribution down to scales needed for water management, studies of ecosystem diversity, and anticipation of wildfires.
We thank two anonymous reviewers and Guido Salvucci, editor for the Journal of Hydrometeorology, for their constructive suggestions. This work was supported by the Lawrence Livermore National Laboratory through a LDRD grant to the Scripps Institution of Oceanography (SIO) via the San Diego Supercomputer Center for the LUCiD project. The USGS provided salary support for MD and DC. The Scripps Institution of Oceanography and the California Energy Com-mission PIER Program, through the California Climate Change Center, provided partial salary support for DC, DP, and HH. The NOAA RISA Program provided partial salary support for DC, Mary Tyree, and Jennifer Johns at SIO. The Program for Climate Model Diagnosis and Intercomparison (PCMDI) supported CB and GB (the former via a DOE Distinguished Scientist Fellowship awarded to B. Santer). TD was partially supported by a CALFED Bay-Delta Program–funded postdoctoral fellowship grant. We thank Andrew W. Wood, Alan Hamlet, Dennis Lettenmaier at the University of Washington, and Edwin P. Maurer at Santa Clara University for sharing VIC and VIC input data. The authors thank Mary Tyree for her support in processing some of the data. We thank Jennifer Johns for her comments on the initial version of the manuscript.
+ Current affiliation: Universidad de Costa Rica, San José, Costa Rica.
& Current affiliation: Centre for Atmospheric and Oceanic Sciences, Indian Institute of Science, Bangalore, India.
Corresponding author address: Tapash Das, Scripps Institution of Oceanography, UCSD, 9500 Gilman Drive–0224, La Jolla, CA 92093. Email: firstname.lastname@example.org
The coarsened gridded meteorological observations of Maurer et al. (2002) from the period 1950–76 and their corresponding high-resolution patterns were used as the library.
The Anderson–Darling test is a modification of the Kolmogorov–Smirnov test in which a test statistic (p) was calculated to assess whether the distribution of the trends in the climate indexes computed using the control run data were drawn from a population with a normal distribution. The null hypothesis that the data (trends in the climate indexes computed using the control run) came from a normal distribution was rejected when the calculated p value was less than a chosen alpha (0.05).