This article applies formal detection and attribution techniques to investigate the nature of observed shifts in the timing of streamflow in the western United States. Previous studies have shown that the snow hydrology of the western United States has changed in the second half of the twentieth century. Such changes manifest themselves in the form of more rain and less snow, in reductions in the snow water contents, and in earlier snowmelt and associated advances in streamflow “center” timing (the day in the “water-year” on average when half the water-year flow at a point has passed). However, with one exception over a more limited domain, no other study has attempted to formally attribute these changes to anthropogenic increases of greenhouse gases in the atmosphere. Using the observations together with a set of global climate model simulations and a hydrologic model (applied to three major hydrological regions of the western United States—the California region, the upper Colorado River basin, and the Columbia River basin), it is found that the observed trends toward earlier “center” timing of snowmelt-driven streamflows in the western United States since 1950 are detectably different from natural variability (significant at the p < 0.05 level). Furthermore, the nonnatural parts of these changes can be attributed confidently to climate changes induced by anthropogenic greenhouse gases, aerosols, ozone, and land use. The signal from the Columbia dominates the analysis, and it is the only basin that showed a detectable signal when the analysis was performed on individual basins. It should be noted that although climate change is an important signal, other climatic processes have also contributed to the hydrologic variability of large basins in the western United States.
Previous studies have found hydroclimatological changes in the last 50 yr in the western United States. The changes are evident in the timing of spring runoff (Roos 1987, 1991; Wahl 1992; Aguado et al. 1992; Pupacko 1993; Dettinger and Cayan 1995; Regonda et al. 2005; Stewart et al. 2005), in the fraction of rain versus snow (Knowles et al. 2006), in the amount of water contained in the snow (Mote 2003), and in climate-sensitive biological variables (Cayan et al. 2001). It has been thought that these changes are related mainly to temperature increases as they affect snowmelt-dominated basins in ways predicted in response to warming (Mote 2003; Barnett et al. 2005; Stewart et al. 2005; Maurer et al. 2007), and it has been suspected that the warming trends causing the changes are in part due to anthropogenic effects. Except for a recent study of California rivers (Maurer et al. 2007), though, no other study has attempted formally to detect and attribute those hydrometeorological changes to anthropogenic effects.
The present article is one of a series of papers describing detection and attribution of the causes of hydroclimatological change in the western United States (Barnett et al. 2008; Bonfils et al. 2008; Pierce et al. 2008). In particular, this paper focuses on shifts in the timing of streamflow. We investigate whether the shifts in streamflow over the past 50 yr are unlikely to have come about by natural variability, and, if so, whether these changes can be confidently attributed to human-caused climate change.
The western United States is particularly susceptible to temperature changes as (historically) a large fraction of the precipitation falling in the mountainous regions of the West occurs on days where the temperature is just a few degrees below 0°C (Bales et al. 2006). Presumably, this precipitation would change from snow to rain in warming climatic scenarios that increase the air temperature by a few degrees. As an example, from 1949 to 2004 the warming of less than 3°C in winter wet-day temperatures across the region has resulted in significant negative trends in the snowfall water equivalent (SFE) divided by the fraction of winter precipitation (P) falling on snowy days (SFE/P) (Knowles et al. 2006). Knowles et al. (2006) also found that, although the Pacific decadal oscillation (PDO; Mantua et al. 1997) may have influenced wet-day temperatures and snowfall fractions at interdecadal time scales, longer-term changes also appear to have been occurring. This is also consistent with the findings of Stewart et al. (2004) and Mote et al. (2005, 2008) with respect to streamflow timing and snow water contents, respectively. While the overall volumes of annual streamflow have not changed much over the past 50 yr, the warming-induced changes are manifested in changes in the interseasonal distribution of streamflow. In particular, the March fraction of annual streamflow has increased, while the April–July fraction has decreased in some basins, with the center timing of streamflow (CT) in snow-dominated basins showing significant shifts toward earlier times in the spring (Dettinger and Cayan 1995; Stewart et al. 2005). Moore et al. (2007) mention that the definition of CT used by Stewart et al. (2005) (i.e., centroid) is similar to the date when 50% of the water-year flow has passed (DQF50) but this later index is less sensitive to outliers in the flow. Consistently with Regonda et al. (2005), Maurer et al. (2007), Moore et al. (2007), Rauscher et al. (2008), and Burn (2008), we will use CT as the date of the year when 50% of the water-year flow has passed (DQF50). In Déry et al. (2009), the CT as calculated here is criticized as a method for computing the timing of streamflow. The authors argue that for certain Canadian rivers there is a correlation between the annual flow and CT and that the influence from late season precipitation and glaciers could affect CT. In our case, the correlation between the flow and CT is nonsignificant at the 95% confidence level, there is only a marginal contribution to the flow from glaciers, and the contribution from summer storms to the flow is negligible compared to the winter values, giving us confidence that we can index the timing of streamflow using CT.
In this study, two global climate model control runs are used to characterize CT natural variability. Using these runs, we will determine whether the trends in the observations are to be found in the distribution of trends from natural variability alone. In those cases where the trends in the observations lie (statistically) significantly outside the distribution of the natural variability, “detection” is achieved (Hegerl et al. 2006). Several forced runs will be used to attribute detected signals to anthropogenic or to solar–volcanic forcings on the climate. In all model runs, we will downscale the data from the climate models to a 1/8°-resolution grid and then run the downscaled estimates through the Variable Infiltration Capacity (VIC; Liang et al. 1994) hydrological model. In contrast to detection, attribution of anthropogenic climate change impacts is the process of determining whether the observed impacts are 1) consistent with the type of changes obtained from climate simulations that include external anthropogenic forcings and internal variability and 2) inconsistent with other explanations of climate change (Hegerl et al. 2006). Detection and attribution studies have been conducted for a number of measures of the climate of the atmosphere and ocean (Barnett et al. 2001; 2006; Hegerl et al. 2006; Hoerling et al. 2006; Zhang et al. 2007; Santer et al. 2007; Solomon et al. 2007b). A review of previous detection and attribution studies is available from The International Ad Hoc Detection and Attribution Group (2005). Most of those previous studies examine global- or continental-scale quantities. This study differs by attempting to perform detection and attribution on a regional scale, which is generally more difficult than the larger-scale analyses because the signal-to-noise ratio is proportional to the spatial scale of analysis (Karoly and Wu 2005). This study also is one of the first formal detection and attribution studies of hydrological variables.
2. Methods, models, and data sources
To investigate the detectability and possible attribution of climate change effects on streamflow timing in the western United States, we employed a particular detection method, several climate model simulations (including control runs and also runs that were driven by anthropogenic forcings or by solar and volcanic forcings), two statistical procedures to downscale the global climate model output to a hydrologically suited terrain scale, a macroscale hydrological model with river routing, a set of observed meteorological data over the western United States, and observed streamflow data for key large basins. These elements are each described in the following sections.
a. The optimal detection method
As noted earlier, “detection” of climate change is a procedure to evaluate whether observed changes are likely to have occurred from natural variations of the climate system. The optimal detection method is applied here. Details of the method can be found in Hegerl et al. (1996, 1997), Tett et al. (1999), Allen and Tett (1999), and Barnett et al. (2001). Given a variable for detection, the basic idea is to reduce the problem of multiple dimensions (n) to a univariate or low-dimensional problem (Hegerl et al. 1996). In this low-dimensional space, the detection of signals above the natural variability “noise” can be contrasted. That is, the trends in observed CT will be compared to the distribution of trends from the control run. Also, the detected vector can be compared with the vector obtained from the expected climate change pattern. As in Hegerl et al. (1996), we are using a simplified version of the method in which the signal strength (S) is actually the trend of the climate vector projected into the fingerprint for each of the climate runs. Specifically, S was defined as
where F(x) is the signal fingerprint, (Dx, t) indicates the three regional time series of any ensemble model run or the observations, and ‘trend’ indicates the slope of the least squares best-fit line. Uncertainty in the signal strength is calculated from a Monte Carlo simulation (see following sections). The 95% confidence intervals computed using traditional t test statistics are approximately 0.91 times the confidence intervals obtained from the Monte Carlo method, suggesting that the autocorrelation and ensemble averaging have some effect on the size of the confidence intervals, with the Monte Carlo intervals being the most conservative estimates. For reference, an assessment of S significance using traditional statistics (Student-t probabilities from the slope of the best-fit line and the nonparametric Mann–Kendall test) will be presented in following sections. In the low-dimensional space, the signal strength will be used to 1) determine if the observations contain a significant signal above the natural variability in the system, as determined from two extensive control runs; and 2) to test the hypothesis that the trends of anthropogenically forced runs have the same signs as the trends in the observations and that these signs are different than the solar–volcanic forced runs.
b. Climate models
An 850-yr control simulation using the Community Climate System Model, version 3.0, Finite Volume (CCSM3-FV; Collins et al. 2006; Bala et al. 2008) general circulation model (GCM) will be used here to characterize natural inherent climate variability in the absence of human effects on climate (hereinafter called the CONTROLccsm run). CCSM3-FV is a fully coupled ocean–atmosphere model, run with no flux corrections, and an atmospheric latitude–longitude resolution of 1° × 1.25° and 26 vertical levels. In addition, 750 yr of a control run (run B06.62) from the Parallel Climate Model version 2.1 (PCM; Washington et al. 2000) at a resolution of T42L26 (hereinafter CONTROLpcm run) will be used to verify our results. The anthropogenically forced signal of climate change will be obtained from 4 simulations (runs B06.22, B06.23, B06.27, and B06.28) by PCM under anthropogenic forcings (hereinafter ANTHROpcm runs) as well as from 10 realizations of the Model for Interdisciplinary Research on Climate, medium-resolution version [MIROC(medres) T42L20; Hasumi and Emori 2004; Nozawa et al. 2007; hereafter ANTHROmiroc runs] also with anthropogenic forcings. Two realizations of the PCM forced with solar and volcanic forcings only (runs B06.68 and B06.69) were obtained from the Intergovernmental Panel on Climate Change (IPCC) Web site (http://www.ipcc-data.org/). These latter data are used to test whether observed natural solar and volcanic effects can explain the observed river flow changes. The characteristics of the models and simulations are indicated in Table 1.
The control runs include only internal variability (no forcings). The ANTHROpcm runs include greenhouse gases, ozone, and direct effect of sulfate aerosols. The ANTHROmiroc runs include the previous forcings, plus the indirect effect of sulfate aerosols, direct and indirect effects of carbonaceous aerosols, and land-use change (see supplementary material in Barnett et al. 2008). The PCM and MIROC anthropogenic runs do not include solar and volcanic forcings, so they are not “all forcings” runs. We used these ANTHRO runs because we wanted to separate the effects of solar–volcanic and anthropogenic effects in the detection procedure. All models selected required to have a good representation of typical sea surface temperature patterns of El Niño–Southern Oscillation (ENSO) and the PDO (Bonfils et al. 2008; Pierce et al. 2008). Another requirement was the availability of daily precipitation and temperature data.
c. Downscaling methods
The data from the CONTROLccsm, solar and volcanic runs, and ANTHROmiroc were downscaled to a ⅛° × ⅛° resolution grid over the western U.S. basins using the method of constructed analogs (CA; Hidalgo et al. 2008). Data from the CONTROLpcm were downscaled to the same resolution using the method of bias correction and spatial disaggregation (BCSD; Wood et al. 2004), as were the ANTHROpcm runs. An intercomparison of the methods of downscaling can be found in Maurer and Hidalgo (2008). The most notable difference between the two methods on the decadal time scales of interest here is that trends are modestly weaker in the CA method than in the BCSD (Maurer and Hidalgo 2008).
1) Constructed analogs downscaling
The CA method is an analogous-based statistical downscaling approach described in detail in Hidalgo et al. (2008). In the appendix a description of the mathematical procedures of the method is included. The downscaling is performed on the temperature and precipitation daily anomaly patterns from the GCM. The base period for computing anomalies is 1950–76. A simple bias correction procedure is accomplished by dividing the anomalies at each grid point of the GCM by the standard deviation of the model and multiplying the result by the standard deviation of the observations. The precipitation was transformed by a square root before processing to make its distribution more Gaussian.
The coarse-resolution target pattern to be downscaled from a climate model for a particular day is estimated using a linear combination of previously observed patterns (library) that are similar to the target pattern. The linear estimate at the coarse scale of the target pattern is called the analog. The downscaled estimate is constructed by applying the same linear combination of coefficients obtained at the coarse-scale to the high-resolution patterns corresponding to the same days used to derive the analog. In this application of the CA, the library patterns were composed of the Maurer et al. (2002) daily precipitation and temperature gridded observations, aggregated at the resolution of each climate model, from 1950 to 1976 along with the corresponding ⅛° versions for the same days. As in Hidalgo et al. (2008), the estimation of the target pattern was constructed by using as predictors the best 30 analogs [based on the pattern root-mean-square error (RMSE) distance from the target] selected from all days in the historical record within ±45 days of the day of year of the target. The domain of the downscaled meteorological data contains the four major hydrological regions of the western United States: the Columbia River basin, California, the Colorado River basin, and the Great Basin (Fig. 1). A list of the characteristics of the basins can be found in Table 2.
2) Bias correction and spatial disaggregation downscaling
The BCSD method is described in detail in Wood et al. (2002, 2004) and has been used previously in a number of climate studies for the western United States (Van Rheenen et al. 2004; Christensen et al. 2004; Payne et al. 2004; Maurer et al. 2007). In brief, the two-step procedure first removes bias more rigorously than for the CA method. The bias is removed on a climate model gridcell-specific basis by using a mapping from the probability density functions for the monthly GCM precipitation and temperature to those of observations, spatially aggregated to the GCM scale. The adjusted climate model simulation outputs from this step are then expressed as anomalies from long-term observed means at the climate model scale. Spatial disaggregation is achieved by the second step, in which the monthlong daily sequences of precipitation and temperature minima and maxima at the 1/8° scale are randomly drawn from the historical record, and then scaled (for precipitation) or shifted (for temperatures) so that the monthly averages reproduce the climate model–scale anomalies. Two constraints are applied: that the selected month is the same calendar month as the month being downscaled; and that the same sample year is applied to all grid cells within the basin for each month, which preserves a plausible spatial structure of precipitation and temperature.
d. Hydrological model
The downscaled precipitation (P), maximum temperature (Tmax), and minimum temperature (Tmin) data along with climatological wind speed from all the runs were used as input to the macroscale VIC (Liang et al. 1994) hydrological model. VIC simulates a full complement of hydrological variables making up the land surface water and energy balance such as soil moisture, snow water equivalent (SWE), baseflow, and runoff, using daily meteorological data as time-varying input, based on parameterized soil and vegetation properties. The land surface is modeled using a tiled configuration of vegetation covers, while the subsurface flow is modeled using three soil layers of different thicknesses (Liang et al. 1994; Sheffield et al. 2004). Defining characteristics of VIC are the probabilistic treatment of subgrid soil moisture capacity distribution, the parameterization of baseflow as a nonlinear recession from the lower soil layer, and that the unsaturated hydraulic conductivity at each particular time step is a function of the degree of saturation of the soil (Campbell 1974; Liang et al. 1994; Sheffield et al. 2004). Details on the characteristics of the model can be found elsewhere (Liang et al. 1994; Cherkauer et al. 2002; http://www.hydro.washington.edu/Lettenmaier/Models/VIC/VIChome.html). The model was run using the water-balance mode at ⅛° resolution over the western United States (Fig. 1). VIC has been used extensively in a variety of water resources applications, from studies of climate variability, forecasting, and climate change studies (e.g., Wood et al. 1997, 2002, 2004; Nijssen et al. 1997, 2001; Hamlet and Lettenmaier 1999).
e. Naturalized streamflow data sources
The naturalized streamflow data for California were obtained from the California Data Exchange Center (CDEC; http://cdec.water.ca.gov/). The data for the Colorado River were obtained from James Prairie’s Internet site from the Upper Colorado Regional Office of the Bureau of Reclamation (http://cadswes2.colorado.edu/~prairie/index.html). The data from the Columbia River were obtained from an updated version of A.G. Crook Company (1993). These naturalized streamflow estimates are generated by adding back the consumptive use to the measurement from the streamflow gauges for each month. Although it is difficult to assess the quality of these naturalized flows, we calculated the streamflow climatologies before and after the major dams were built in the rivers. The results showed that the Columbia River has the closest match of these climatologies, followed by the Colorado, and the greatest differences were found for the California rivers (Fig. 2). Overall the differences in the climatologies are not large for the Colorado and Columbia, supporting the use of the naturalized data. The differences are larger for the California rivers, and this may affect the results for individual basins. Although, as it will be seen, because of weighting used the Columbia plays a dominant role in the detection and attribution for the western United States; therefore, the influence of the California rivers in the western U.S.–wide detection and attribution analysis is marginal.
f. River routing
The VIC-simulated runoff and baseflow were routed, as in Lohmann et al. (1996), to obtain daily streamflow data for four rivers: 1) the Sacramento River at Bend Bridge (California), 2) the Colorado River at Lees Ferry (Arizona), 3) the Columbia River at The Dalles (Oregon), and 4) the San Joaquin River (California). The data for the total flow of the San Joaquin River were obtained by adding together the daily streamflow values from the four main tributaries: the Stanislaus, the Merced, the Tuolumne, and the San Joaquin Rivers (Fig. 1). The naturalized monthly observed data from these four rivers showed strong correlations with the values obtained from the VIC model using the meteorological data from Maurer et al. (2002) as input, suggesting that the calibration of the model is good (Fig. 3). These rivers present runoff climatologies consistent with snow-controlled hydrologies with a maximum peak occurring around May or June (Fig. 4), although the Sacramento River is fed in part by runoff from lower elevations as well and peaks in February–March.
Although the variability of streamflow and CT is captured well, the trends in CT by VIC-forced-by-gridded-observations (from Maurer et al. 2002) are all negative but weak (nonsignificant) for all rivers (Fig. 5, left panel). This underestimates the negative trend in the Columbia River, which is highly significant for the CT computed from the naturalized flow (Fig. 5, bottom right). There are three possible explanations for the weaker trends in CT computed from VIC-forced-by-gridded-observations compared to the CT computed from the naturalized streamflow: 1) the naturalized streamflow has errors, 2) the forcing data (precipitation, temperature, and wind speed) for VIC have errors, and 3) a deficiency in VIC to reproduce CT time series. It is important to determine if there are serious deficiencies in VIC to produce CT time series, as this study depends on the correct modeling of CT; for this reason an analysis was developed to look at the sources of errors from CT.
First, in order to discard possible large errors in the naturalized streamflow data for the Columbia, we calculated the 1950–99 trend in CT from a collection of streamflow gauges from the Hydro-Climatic Data Network (HCDN), updated using data from the United States Geological Survey. These gauges are relatively unimpaired by dams. As can be seen in Fig. 6, the trends in CT in the high elevations of the Columbia River have magnitudes on the order of −0.2 days per year, consistent with the CT trend of −0.17 days per year from the naturalized streamflow of the Columbia River at The Dalles (Fig. 5). Note that coastal stations showed a trend toward later streamflow CT (Fig. 6). Stewart et al. (2005) showed this opposite response in CT for coastal stations that are not snow-dominated gauges. Although the consistency of the CT trends at The Dalles compared to Fig. 6 is not a formal indication of the quality of the naturalized data, the previous analysis gives us confidence that there is a significant trend in the CT on the high-elevation parts Columbia basin, a result found also in Stewart et al. (2004).
Second, we looked at possible errors in the trends of the forcing data used to drive VIC. In Fig. 7, the basin average temperature and precipitation from the Maurer et al. (2002) data1 (VIC forcing) and from the area-weighted average of the Climate Divisional data for Washington State are shown. The trends in temperature in the VIC forcing data are slightly weaker than the trends in the Climate Divisional data, while the trends in precipitation are positive in the VIC forcing data and almost zero in the Climate Divisional data. One may ask if these differences are large enough to produce weaker trends in CT in the VIC forced by observations compared to the naturalized flows. We modified the VIC input by increasing the temperature trend in an amount corresponding to the difference in the trends from the Climate Divisional data and the VIC forcing and decreased the precipitation trend in the data using the same criteria. The resulting dataset (modified VIC forcing) has very similar basinwide precipitation and temperature trends compared to the Climate Divisional trends (not shown). When CT was calculated using the modified VIC forcing, the CT trend becomes strongly negative and highly significant (β = −0.45 days yr−1, p < 5.61 × 10−6). Although it is possible that the Climate Divisional data have errors (Keim et al. 2003), this approximate experiment suggests that there may be errors of such magnitude in the trend of the forcing data that are sufficient to diminish the trends in CT as shown in Fig. 5 (left panels) and that the VIC model does not seem to have serious deficiencies in modeling CT trends.
g. Calculation of the signal strength S
For all runs, the Sacramento and San Joaquin Rivers were added to form a single streamflow time series representative of the California region, leaving us with three streamflow time series: the Sacramento/San Joaquin, Colorado, and Columbia Rivers. The CT of streamflow was computed from the simulated daily streamflow time series and was defined according to Maurer et al. (2007) as the day of the water year when 50% of the water-year streamflow has passed through the channel at the calibration points shown in Fig. 1. The CT of the observed naturalized flows was estimated from the monthly data by allocating the monthly values to the middle days of the months and interpolating the CT values between the months that correspond to the point where the fractional flow was below and above 50%. This procedure proved to be accurate using as an example the daily data from the models by comparing the results of computing the CT from the dailies and the monthlies (not shown). We therefore have three time series (one per river), each 50 yr long, for the observations and for every anthropogenically forced model run and 50-yr segment of the control runs.
The fingerprint is defined as the leading empirical orthogonal function (EOF) of the ensemble-averaged CT time series of the PCM and MIROC anthropogenically forced ensemble members. That is, ensemble CTs for each of the three rivers were obtained by averaging 14 members (4 from PCM and 10 from MIROC). We wanted to focus mostly on changes in river flow driven by snowmelt, so, prior to computing the EOF, we weighted each CT series by a factor equaling the basin’s climatological ratio of 1 April SWE divided to water-year precipitation (P) using the data from the control runs [1 April is typically (within 12%) the date of maximum SWE accumulation in the western United States (Bohr and Aguado 2001)]. This choice emphasized rivers driven primarily by snowmelt and deemphasized rivers driven primarily by rainfall. A second set of weights accounted for the area of the basins so time series representing larger areas would have proportionally more influence. Both types of weights, expressed as fractions, are shown in Table 3.
The EOF that comprises the fingerprint explains 78% of the variance. The resulting component is heavily biased toward the variability of the Columbia (because of the weighting), and therefore the majority of the signal comes from this region. The fingerprint pattern is therefore a pattern with large loadings in the Columbia and small loadings for the California and Colorado basins. In a following section the detection for individual rivers will be provided.
The standard deviation from the first 300 yr of the CONTROLccsm run was initially used to optimize the signal-to-noise ratio. Optimization is a process used in certain kinds of detection and attribution studies to accentuate the signal-to-noise ratio, but it requires part of one of the control runs to be used for optimization purposes (and not allow to be part of the detection). In the present study, the optimized results differ little from the nonoptimized version (not shown). We therefore use the nonoptimized solution to allow the entire CONTROLccsm run to be used for detection purposes.
It is important to examine whether the variability of the center timing in the control runs is similar to the variability in the observations. If the variability of the control runs is significantly less than the observations, there is a risk of spurious detection because forced trends in the observations could be significantly higher than the variability from the control run but not the variability in the real world. A comparison of the spectra of a tree-ring reconstruction of the annual streamflow at the upper Colorado River at Lees Ferry (Meko et al. 2007) with the control runs (Fig. 8a) indicated good agreement, suggesting that the low-frequency streamflow variability is well captured by the models (see also supplementary information from Barnett et al. 2008). However a similar comparison of the spectra of the streamflow at the Columbia River at The Dalles from the control runs with a tree-ring reconstruction (Gedalof et al. 2004) show that the control runs generally overpredict the variability of the streamflow at the Columbia for frequencies higher than 0.05 cycles yr−1 or 20-yr periodicities (Fig. 8b). Note that the streamflow low-frequency variability of interest here is captured well by the models, and that, although Fig. 8 is an indicator of the agreement of the annual streamflow spectra between tree-ring data and models, it does not show the agreement between CTs (CT cannot be computed from tree-ring streamflow reconstructions that have annual resolution). For this reason we computed the standard deviations for the 5-yr low-pass-filtered CT for the control runs and naturalized observations and found that the standard deviations are statistically the same. The 5-yr low-pass-filtered CT was used to remove any high-frequency variability, for example, associated with ENSO. The standard deviations for the Columbia River CT are 5.9 days for the CONTROLccsm, 5.5 days for the CONTROLpcm, and 5.5 days for CT from the naturalized flow observations.
The slopes of fitted linear trends of CT from 1950 to 1999 for the naturalized flows are negative in the three rivers, although only in the Columbia are the trends significant (Fig. 5, right panel). In the ANTHRO ensemble all trends are negative with significant trends in the Colorado and Columbia (Fig. 9). By contrast, for the PCM solar–volcanic runs the trends are in all cases positive with significant trends in the Columbia (Fig. 10). Note that the CT trends in the ANTHRO runs are not strictly comparable to the trends in the observations. That is, the ANTHRO runs are not an “all forcings” run (of a similar type to a twentieth-century climate, or 20c3m, run), and solar and volcanic effects for example are not included.
It is interesting to note that California has been one of the first places where the earlier snowmelt has been reported; therefore, the lack of significant CT trends in the observations and in the model is puzzling. However, if we look at the longer CT records from 1907 to 2003, the observed trends in California are highly significant (β = −0.183 days yr−1, p = 0.00046), therefore significant trends can be found when using the long-term data (not shown).
The resulting detection plot is shown in Fig. 11. The fingerprint is shown in Fig. 11a. The detection variable S is shown in Fig. 11b. We used a Monte Carlo test (below) to estimate the likelihood that the model runs and the observations are drawn from the control distribution. The values of S, including the resulting 95% confidence Monte Carlo error bars, are positive for the observed naturalized flow and the ANTHRO runs. By contrast, S is negative for the solar and volcanic runs. This means that the trends toward earlier river flow observed in the model runs and the observations are unlikely to have been obtained from natural internal variability alone. For reference, an assessment of S significance using traditional statistics is shown in Table 4.
It is important to assess the probability that the distribution of S in the ANTHROpcm and ANTHROmiroc runs (or the observations) is significantly different from the distribution of S in the control runs. This was calculated using a Monte Carlo method that estimates the likelihood that a given ensemble mean value of S can be drawn from the control runs, given an ensemble of k members (k = 4 for ANTHROpcm, k = 10 for ANTHROmiroc, k = 2 for the solar–volcanic runs). Groups of k members were randomly selected from among all the 50-yr segments in the control runs, and their ensemble-average S calculated. This was repeated 10 000 times to form a distribution of control S for comparison with the anthropogenic models. The same procedure was used for the observed naturalized flow, although in this case with k = 1, no ensemble averaging is possible. The ANTHROpcm and ANTHROmiroc ensemble means are unlikely to have been drawn from the control distribution (p < 0.05). The observed naturalized flow also differs from the control run at the p < 0.05 level. This implies that the human influences on climate are discernible from the natural variability; that is, detection has been achieved at the 95% confidence level.
Attribution is addressed by determining whether the observed values of S are not inconsistent with the anthropogenic and/or solar–volcanic model results. Assuming a normal distribution, the means and standard deviations of the ANTHROpcm and ANTHROmiroc S values were calculated. The difference between the S from the observations and the ensemble mean of the S values from all the ANTHROpcm and ANTHROmiroc runs is statistically small; hence the trends in the observations are consistent with the trends from the anthropogenic runs. On the other hand, the ensemble-mean S from the solar and volcanic runs is more than four standard deviations away from the observations and more than four standard deviations away from the mean of the anthropogenic runs (not shown), suggesting they are from different statistical distributions than the observations. We conclude that observed changes in river CT are consistent with human forcing of the climate and unlikely to have arisen from natural solar or volcanic variability.
We tested several configurations of the western U.S. model, as shown in Table 5, and find that the model is robust with respect to the choices selected. That is, for all combinations of models used in the fingerprint and control run used for the “noise” detection at the p < 0.05 level was found (Table 5). The only case that we found where formal detection failed occurred when we did not use the SWE/P and area weights (not shown).
Following results from Stewart et al. (2005) and Dettinger and Cayan (1995), which indicate shifts toward greater earlier streamflow fractions in the western U.S. basins, we applied the detection method to the streamflow fractions during the winter–spring, summer, and the March fractions (not shown). For the winter–spring fractions we found detection at the p < 0.05 level for all of the GCMs (not shown). For summer the observations do not show a strong enough trend to result in detectability at the 0.05 level (not shown). For the March fractions, the ANTHROpcm and the observations showed detection at the 0.05 level, but the ANTHROmiroc failed to trend significantly (not shown). The weaker trends in MIROC compared to PCM and the observations indicate that MIROC did not warm realistically (Bonfils et al. 2008).
We also investigated whether the detected changes were associated with temperature or precipitation changes. Pierce et al. (2008) repeated the entire detection and attribution analysis using P instead of SWE/P and found no detection at the p < 0.10 level. That study (Pierce et al. 2008) and previous studies (Mote et al. 2005; Stewart et al. 2005) concluded that the reductions in snowpack are primarily driven by increases in temperature over the western United States (Pierce et al. 2008). Thus we believe that the temperature increases in winter and spring and associated reductions of snow versus rain ratios (Knowles et al. 2006) and in the spring snowpack are responsible for the advanced timing of streamflow.
Maurer et al. (2007) found no detection for California river flow using CT as a detection variable. We repeated our detection and attribution analysis using individual rivers and likewise did not achieve formal detection for the Sacramento and the San Joaquin Rivers (Figs. 12a,b). For the Sacramento River, the CT observations and the ANTHROmiroc fall within the overlap area where they are consistent with both the anthropogenic results and the distribution of natural internal climate variability. However, the PCM anthropogenic model runs (ANTHROpcm) proved to be separate from the distribution of the control runs (CONTROLccsm and CONTROLpcm). In such a case, although there is partial separation, detection of an anthropogenic effect cannot be irrefutably claimed. This example illustrates the present difficulty of regional detection and attribution, when averaging over relatively small regions often produces a distribution of anthropogenically forced responses that is not yet well separated from that of natural internal variability, especially when the character of simulated natural variability varies from model to model. In the case of the San Joaquin River and the upper Colorado River, the ANTHROpcm and ANTHROmiroc are well separated from the distribution of trends from the control run (Figs. 12b,c); however, the observations fall between the overlap area where they are consistent with both the anthropogenic results and the distribution of natural internal variability. In these cases we cannot claim that detection is achieved. It is possible that the lack of trends in the observations is due to the masking of the anthropogenic signal by the opposing effect of solar–volcanic effects (as the observations are an “all forcing” integration). It should be emphasized that, in some cases, the ANTHRO runs predict a significant decline in the CT trends that is not seen in observations (note that the fingerprint has negative loadings so the sign of the trends are inverted in the Fig. 12). It should be kept in mind that the ANTHRO runs are not “all forcings” runs—they lack solar and volcanic effects, for example—and so the PCM model response to greenhouse gases only is not affected by cooling effects from the solar–volcanic forcings. Note also that, for the Sacramento case, the ANTHROmiroc does not separate well from the zero line, while the ANTHROpcm is well separated. The weaker trends in MIROC compared to PCM may be an indication that MIROC did not warm realistically (Bonfils et al. 2008), but this finding awaits a longer record to be conclusive.
In the case of the Columbia (Fig. 12d), the detection and attribution of climate change is evident in the separation of observed and modeled trends from the distribution of trends from the control run. In the observations, the effect of anthropogenic warming is greater than any cooling effect from solar and volcanic sources and therefore results in strong negative trends. One observation regarding the results from the Columbia is the fact that, if we assume that the effects of climate change are linear, one can add the trends from the ANTHRO in Fig. 9 and the solar–volcanic from Fig. 10 and obtain near-zero trends for the observations. At first glance this seems inconsistent with the significant trends found for the Columbia in the naturalized data. It can be argued that either some forcings beside the considered anthropogenic and natural are missing or the model is not able to reliably reproduce the response of the considered forcings. Although these are two valid possibilities, it should be mentioned that we are not claiming that the entire trend in the Columbia is associated with global warming, just a fraction of it (e.g., in Barnett et al. 2008 it was estimated that for multiple variables analyzed that fraction is around 60%). Therefore other effects (e.g., the fact that the PDO switched around the middle of the 1950–99 period, increasing temperatures in the Columbia basin) may have played a role in the observed CT trends in the Columbia. We believe this does not invalidate the attribution part of this analysis as it has been proven in other sources that the shift in the PDO is not enough to explain the change in hydrological measures in the western United States (Stewart et al. 2004; Knowles et al. 2006; Mote et al. 2005, 2008).
A formal attribution and detection procedure can provide insight into the nature of observed changes in the streamflow timing from key snowmelt watersheds over the western United States. If the changes toward earlier timing of streamflow are found to be associated with natural variability, it would be expected that after some time the climate system would rebound toward later streamflow timings, and the hydrology would revert toward more snow (and less rain) in the winter and therefore higher flows in summer. If the changes observed are unequivocally associated with anthropogenic warming (because of changes in the composition of the atmosphere), however, the decreases in winter snow to rain ratios and the timing of snowmelt can only become more pronounced as ongoing changes in atmospheric chemistry become more acute.
Using an optimal detection method, we found that observed streamflow center timing (CT) trends lie beyond a good share of the distribution of trends from simulations of natural variability. We find a detectable signal (at the p < 0.05 level) on the timing of streamflow over the second half of the twentieth century toward earlier streamflow timing.
The changes in streamflow timing are dominated by changes over the Columbia River basin, with lesser signals arising from the California Sierra watersheds and little from the Colorado River basin. This indicates that climate change is an important signal, but also indicates that other climatic processes have also contributed to the hydrologic variability of large basins in the western United States.
The present study employed two control runs (CCSM3-FV and PCM) and two anthropogenically forced models (PCM and MIROC), allowing us to test several configurations of the streamflow timing (CT) detection options by using one or the other control and anthropogenic runs. For CT, the big picture from these experiments indicates a shift toward earlier streamflow, which cannot be explained solely by natural variability. All the options in the selection of the models to use for control and anthropogenic runs resulted in detections at the 95% confidence level, indicating that the results are robust with respect to these choices. For streamflow fractional timing PCM showed positive detection for winter, summer, and March, but MIROC yielded detection only in summer. For all model cases, the Columbia basin was the major contributor to the detection with less influence from California and little from the Colorado. Comparison of two downscaling methods (Barnett et al. 2008) showed little dependence of the detection and attribution results to with respect to the downscaling method. In summary, we can now state with “very high confidence” (Solomon et al. 2007a; box TS.1.1) that recent trends toward earlier streamflows in the Columbia basin are in part due to anthropogenic climate change.
In the cases when detection was positive, we tested the attribution of those changes to two possible explanations: anthropogenic forcings or to natural forcings. In all cases the attribution was consistent with the anthropogenic forcing explanation and inconsistent with the solar–volcanic (natural forcings) explanation. Consistent with the results from Bonfils et al. (2008) for temperature and Pierce et al. (2008) for snow changes, the advance in streamflow timing in the western United States appears to arise, to some measure, from anthropogenic warming. Thus the observed changes appear to be the early phase of changes expected under climate change. This finding presages grave consequences for the water supply, water management, and ecology of the region. In particular, more winter and spring flooding and drier summers are expected as well as less winter snow (more rain) and earlier snowmelt.
This work was supported by the Lawrence Livermore National Laboratory through an LDRD grant to the Scripps Institution of Oceanography (SIO) via the San Diego Supercomputer Center (SDSC) for the LUCSiD project. The MIROC simulations were supported by the Research Revolution 2002 of the Ministry of Education, Culture, Sports, Science and Technology of Japan. The PCM simulation had previously made available to SIO by the National Center for Atmospheric Research for the ACPI project. This work was also partially supported by the Department of Energy and NOAA through the International Detection and Attribution Group (IDAG). The LLNL participants were supported by DOE-W-7405-ENG-48 to the Program of Climate Model Diagnoses and Intercomparison (PCMDI). The USGS and SIO provided partial salary support for DC and MD at SIO; the California Energy Commission provided partial salary support for DC, DP, and HH at SIO. Thanks are also due to the Department of Energy, which supported TPB as part of the IDAG. Thank you to David Meko and Ze’ev Gedalof for providing tree-ring reconstructions of the streamflow of the upper Colorado River and Columbia River basins. We thank two anonymous reviewers, Stephen Déry and the chief editor of the Journal of Climate, Andrew Weaver, for their constructive comments.
Downscaling with Constructed Analogs
Step 1: Fitting a coarse-resolution analog
Once the pool of predictor patterns has been selected for a given coarse-resolution Tmax, Tmin, or P pattern for a certain day and year (𝗭obs), an analog of that pattern (𝗭analog) can be constructed as a linear combination of the (preferred 30-member most suitable subset of) predictor patterns, according to
where 𝗭analog is a matrix of the column vectors comprising the most suitable subset of coarse-resolution patterns identified above specifically for 𝗭obs, and Aanalog is a column vector of fitted least squares estimates of the regression coefficients that are the linear proportions of the contributions of each column of 𝗭analog to the constructed analog. The dimensions of the 𝗭obs matrix are pcoarse × 1, where pcoarse is the number of considered grid points contained in each coarse-resolution weather pattern; that is, 𝗭obs is a column vector. The dimensions of 𝗭analog are pcoarse, x × n where n is the number of patterns in the most suitable predictors subset (i.e., 30), and the dimension of Aanalog is n × 1.
Assuming 𝗭analog has full rank (n) and using the definition of the pseudoinverse (Moore–Penrose inverse), Aanalog is obtained from Eq. (A1) by
where the ′ superscript denotes the transpose of the matrix. The inversion of the matrix was performed using singular value decomposition routine from Press et al. (1992), in which small values of the decomposition were set equal to zero to avoid near-singular matrices.
Step 2: Downscaling a weather pattern
To downscale the 𝗭obs pattern, the coefficients Aanalog from Eq. (A2) are applied to the high-resolution weather patterns corresponding to the same days as the coarse-resolution predictors 𝗭analog, according to
From Eq. (A2),
where is a constructed high-resolution analog (e.g., a P pattern on the VIC ⅛° grid) and 𝗣analog is the set of high-resolution historical patterns corresponding to the same days as the 𝗭analog. The dimension of the downscaled vector is pVIC × 1, and the dimension of the 𝗣analog matrix is pVIC × n, where pVIC is the number of grid points in the high-resolution weather patterns. Note that the matrix, 𝗭′analog𝗭analog, inverted with each application of the procedure is only of dimension n × n, and therefore the numerical computational resources needed to downscale the weather patterns are determined by the number of the patterns included in the most suitable subset and can be quite small.
+ Current affiliation: Universidad de Costa Rica, San José, Costa Rica.
& Current affiliation: Center for Atmospheric and Oceanic Sciences, Bangalore, India.
Corresponding author address: Hugo G. Hidalgo, CASPO Division, Scripps Institution of Oceanography, University of California, San Diego, 9500 Gilman Drive, La Jolla, CA 92093-0224. Email: email@example.com
An analysis using the alternative Hamlet and Lettenmaier (2005) meteorological data was also performed with similar results.