Results from a regional climate model simulation show substantial increases in future flood risk (2040–69) in many Pacific Northwest river basins in the early fall. Two primary causes are identified: 1) more extreme and earlier storms and 2) warming temperatures that shift precipitation from snow to rain dominance over regional terrain. The simulations also show a wide range of uncertainty among different basins stemming from localized storm characteristics. While previous research using statistical downscaling suggests that many areas in the Pacific Northwest are likely to experience substantial increases in flooding in response to global climate change, these initial estimates do not adequately represent the effects of changes in heavy precipitation. Unlike statistical downscaling techniques applied to global climate model scenarios, the regional model provides an explicit, physically based simulation of the seasonality, size, location, and intensity of historical and future extreme storms, including atmospheric rivers. This paper presents climate projections from the ECHAM5/Max Planck Institute Ocean Model (MPI-OM) global climate model dynamically downscaled using the Weather Research and Forecasting (WRF) Model implemented at 12-km resolution for the period 1970–2069. The resulting daily precipitation and temperature data are bias corrected and used as input to a physically based Variable Infiltration Capacity (VIC) hydrologic model. From the daily time step simulations of streamflow produced by the hydrologic model, probability distributions are fit to the extreme events extracted from each water year and flood statistics for various return intervals are estimated.
1. Introduction and background
The combined effects of land use changes and projected climate change have raised substantial concerns about the future vulnerability of ecosystems and infrastructure worldwide to extreme hydrologic events including drought and flooding. Over the past several decades, there has been a clear trend toward increased heavy precipitation in many regions (Easterling et al. 2000; Groisman et al. 2005; Min et al. 2011), including parts of the northwestern United States (Dulière et al. 2013; Mass et al. 2011). In contrast, negative trends are observed in some locations, and the observational evidence for changes in the frequency, duration, and intensity of extreme precipitation events is only now emerging.
There is considerable uncertainty whether local changes in extreme precipitation and flooding are related to global climate change and will persist over the next few decades. Natural variability is clearly important and is likely to continue to be a primary driver of local changes in the near term. Furthermore, there is a complex array of interacting processes than can influence changes in hydrologic extremes. For example, the direct effect of warmer temperatures, elevated freezing levels, and reduced mountain snowpack would have a substantial effect on the hydrologic response to heavy precipitation, irrespective of trends in precipitation itself (Tohver et al. 2014). Thus, over the next few decades, several climatic processes could change the frequency or intensity of extreme hydrologic events.
Heavy precipitation in the Pacific Northwest (PNW) is most frequently associated with atmospheric river events (Colle and Mass 1996; Garvert et al. 2007; Neiman et al. 2011; Warner et al. 2012). Changes in the jet stream and storm-track dynamics (Chang 2007; Salathé 2006; Ulbrich et al. 2009) or atmospheric rivers (Leung and Qian 2009; Neiman et al. 2008a,b) associated with climate change could have substantial implications for future extreme hydrological events. Global climate models project an increased risk of more frequent extreme precipitation in the PNW by the second half of the twenty-first century (Tebaldi et al. 2006) with intensifying atmospheric rivers along the West Coast (Dettinger 2011) and enhanced precipitation in region complex terrain (Salathé et al. 2010). Studies of regional climate model simulations of future climate (Salathé et al. 2010; Dominguez et al. 2012) have found consistent increases in extreme precipitation over the northwestern United States despite modest and sometimes negative changes in total precipitation. Alterations in heavy precipitation events could affect 1) the temporal and spatial extent of heavy precipitation, 2) the temperature anomaly and freezing level associated with storms, 3) the likelihood of storms occurring with specific antecedent conditions (e.g., above-normal soil moisture or snow), and 4) the orientation of storms relative to local terrain. These effects could have profound implications for the hydrologic response to extreme events that go beyond simple changes in rainfall amount, requiring an integrated approach to modeling precipitation and surface hydrology.
Tohver et al. (2014) projected future changes in flood risk in the Pacific Northwest for climate change scenarios using statistically downscaled global climate model output as driving data for a physically based hydrologic model (Hamlet et al. 2013). Flood risk was assessed from daily streamflow simulations at 297 river locations under natural conditions by fitting generalized extreme value (GEV) probability distributions to simulated annual peak flows and estimating the 20-, 50-, and 100-yr flood following methods developed by Hamlet and Lettenmaier (2007). River locations for Tohver et al. (2014) were selected over the Columbia River basin (CRB) and the smaller basins of Washington and Oregon west of the Cascade Range (see Fig. 1). These simulations project widespread increases in flooding for the twenty-first century because of the combined effects of increasing cool season precipitation and rising snow levels during storms associated with warmer temperature (Tohver et al. 2014). The largest increases in flooding in the PNW were generally simulated for river basins with temperature sensitive snowpack, typically west of the Cascade crest.
Earlier research provides a good sense of future PNW flood magnitudes, including the sign of changing flood risks, and identified important mechanisms contributing to the impacts of climate change on flood risk. The results from these studies, however, are very conservative because the statistical downscaling methods were based on monthly climate model output, and thus, no explicit changes in daily precipitation probability distributions were incorporated. Daily precipitation was assumed to simply scale with projected changes in monthly statistics. In contrast, using high-resolution, regional-scale climate models to project future climate change in the PNW, Salathé et al. (2010) showed robust increases in precipitation intensity on the windward slopes, and decreases on the leeward slopes, of the Cascades and Rockies. These effects cannot be captured by global models, which do not resolve realistic terrain features and mesoscale precipitation processes. Statistical downscaling methods only weakly represent these effects, and only when large-scale winds are taken into account (e.g., Widmann et al. 2003).
To improve the simulation of future extreme storms and their effects on precipitation and runoff production over complex topography in the PNW and Intermountain West, we have applied regional climate model (RCM) simulations for streamflow projections over the PNW. In this paper, we present methods for correcting bias in an RCM simulation and linking the results to a hydrologic model. The regional climate simulations are then used to assess future flood risk over the entire PNW for the same river basins assessed by Tohver et al. (2014).
2. Model configuration
a. Regional climate model
The Weather Research and Forecasting (WRF) Model is a state-of-the-art mesoscale numerical weather prediction system designed to serve both operational forecasting and atmospheric research needs (www.wrf-model.org). This model has been developed and used extensively in recent years for regional climate simulation (Leung et al. 2006). WRF is a nonhydrostatic model with multiple choices for physical parameterizations suitable for applications across scales ranging from meters to thousands of kilometers. The WRF physics package includes microphysics, convective parameterization, planetary boundary layer (PBL), land surface models (LSMs), and short- and longwave radiation. The WRF Model has been implemented as a regional climate model over the northwestern United States at 12-km grid spacing as described in Zhang et al. (2009). In particular, the simulations use the following parameterization choices: the WRF single-moment 5-class (WSM5) microphysics scheme (Hong et al. 2004), the Kain–Fritsch convection (Kain and Fritsch 1993), the Noah land surface model (Chen and Dudhia 2001), the Yonsei University (Hong and Pan 1996) boundary layer, and the National Center for Atmospheric Research (NCAR) Community Atmosphere Model (CAM) short- and longwave radiation (Collins et al. 2004).
The simulations in this paper used nested grids at 36- and 12-km spacing. The 12-km nest covers the Pacific Northwest region shown in Fig. 1 and includes the locations of the 297 river locations used for assessing flood risk (see Fig. 4, described in greater detail below). Climate simulations are performed using 6-hourly forcing fields from global reanalysis or global climate models, typically with grid spacing of 150–200 km. The outer 36-km nest receives boundary conditions and interior nudging from the global fields; the inner 12-km nest is forced only at its boundary by the outer nest (i.e., without nudging). Gridded analysis nudging was applied only to the outer nest, only on upper air wind, temperature, and moisture fields, and only on wind in the PBL. To maximize the model’s ability to simulate mesoscale features and to obtain stability for long simulations, nudging coefficients were set to ⅓ the default values for temperature and wind (10−4) and to near zero for moisture (10−6). This study uses a 100-yr (1970–2070) simulation with the WRF Model using boundary conditions from the ECHAM5 global climate model (see section 2b).
The WRF Model configuration used here has been extensively evaluated using simulations forced by reanalysis fields (Zhang et al. 2009) and global climate models (Dulière et al. 2013). In addition, a similar configuration is used for numerical weather prediction over the same domain (Mass et al. 2003). This model configuration is capable of resolving the finescale structure of storms and their effects on precipitation in complex terrain (Zhang et al. 2009; Dulière et al. 2011). In particular, the model successfully simulates important large-scale features of PNW winter storms, such as atmospheric rivers, that have been shown to be the major cause of the largest floods in rivers that drain the windward (western) slopes of the Cascades (Neiman et al. 2011; Warner et al. 2012).
For this study, we use simulated WRF daily output of total precipitation, maximum and minimum temperature, and mean wind speed. As described in section 3, precipitation and temperature variables are bias corrected using quantile mapping approaches to remove deficiencies in the simulated climate relative to the observed climate.
b. Large-scale forcing data
Large-scale forcing for the WRF Model was taken from a simulation of the ECHAM5/Max Planck Institute Ocean Model (MPI-OM) global climate model completed as part of phase 3 of the Coupled Model Intercomparison Project (CMIP3; Covey et al. 2003). The atmospheric component is the fifth-generation general circulation model developed at the European Centre for Medium-Range Weather Forecasts and the Max Planck Institute for Meteorology (Roeckner et al. 1999, 2003), and the ocean component is the MPI-OM. Here, we will refer to the coupled model simply as ECHAM5. For the period 1970–99, we used an ECHAM5 simulation of the twentieth century with historical greenhouse gas, aerosol, and solar forcing; for the twenty-first century, we used an ECHAM5 simulation with the Special Report on Emissions Scenarios (SRES) A1B emissions scenario (Nakićenović et al. 2000). ECHAM5 was run at T63 spectral resolution, which corresponds to a horizontal grid spacing of approximately 140 km × 210 km at midlatitudes, and 32 levels in the vertical. Model output at 6-hourly intervals was obtained from the Climate and Environmental Retrieval and Archive (CERA) Gateway (http://cera-www.dkrz.de/CERA/index.html); the data are managed by the World Data Center for Climate (www.dkrz.de/daten/wdcc/). Simulations with WRF forced by ECHAM5 were performed for the period 1970–2069 and will be referred to as ECHAM5-WRF.
ECHAM5 was the highest-ranked global model based on criteria established by Mote and Salathé (2010) for reproducing the observed PNW climate of the twentieth century, and the ECHAM5 A1B results fall close to the average temperature and precipitation changes for an ensemble of 20 climate models from the CMIP3. The Intergovernmental Panel on Climate Change (IPCC) A1B emission scenario is a medium-high scenario that reflects “business as usual” in the first half of the twenty-first century (i.e., little greenhouse gas mitigation) followed by greater mitigation in the second half of the century as impacts intensify.
c. Regional climate simulations
The WRF Model simulates mesoscale dynamics and orographic effects consistent with the large-scale forcing. Simulated changes in heavy precipitation between the current and future climate depend both on changes in the large-scale forcing and in the consequent mesoscale effects, which control the local intensity of precipitation and its distribution relative to the terrain. The storm simulated for 27 November 2030 (Fig. 2) is characteristic of an atmospheric river event. Such atmospheric rivers are the predominant cause of heavy precipitation along the West Coast and a common cause of flooding in PNW river basins, like the Skagit (see Fig. 1), that drain the western slopes of the Cascades (Neiman et al. 2011; Warner et al. 2012). Global models project that more intense atmospheric rivers will occur earlier in autumn under future climate scenarios, including the ECHAM5 simulation used in this study (Dettinger 2011). Thus, the realistic simulation of the mesoscale properties of such systems is critical in assessing flood risk under climate change scenarios.
In addition to changes in the intensity of atmospheric rivers, changes in the nature of future storms may lead to regional redistributions of flood risk. In particular, the wind direction has a substantial effect on the distribution and intensity of precipitation in atmospheric river events because of their interactions with orography (Neiman et al. 2011). Thus, zonally oriented atmospheric rivers result in a very different distribution of streamflow in river basins across the region than a more southwesterly atmospheric river. For these reasons, the ability to simulate the mesoscale precipitation features using the WRF Model provides the ability to incorporate multiple pathways for climate simulations to influence projected flood risk in a region of complex terrain.
Figure 3 shows the simulated change in the 30-yr average annual maximum daily precipitation from the current (1970–99) to future (2040–69) climate periods from ECHAM5-WRF. Substantial increases are simulated over Oregon and British Columbia. Consistent with previous modeling studies (e.g., Salathé et al. 2010; Dominguez et al. 2012), the simulation shows a substantial increase in heavy precipitation on the windward slopes of coastal terrain and the Oregon Cascades. These simulated changes in heavy precipitation are a primary driver of projected increases in flood risk discussed below.
d. Hydrologic model
The Variable Infiltration Capacity (VIC) model (Liang et al. 1994, 1996; Nijssen et al. 1997) is a macroscale, fully distributed hydrologic model that solves the water and energy balance at each model grid cell, producing (among other water balance variables) daily time step runoff and base flow at each model grid cell. The VIC model implemented at 0.125° latitude–longitude resolution (~130 km2 per grid cell) has been used extensively to assess the impacts of climate change on hydrology of the CRB (Hamlet and Lettenmaier 1999; Payne et al. 2004; Lee et al. 2009). The spatial resolution of the model over the PNW has recently been increased to 0.0625° latitude–longitude grid spacing (°, approximately 5 km × 7 km or 35 km2 per grid cell; Elsner et al. 2010; Hamlet et al. 2013).
The VIC model requires daily inputs of total precipitation, maximum and minimum air temperature, and mean wind speed. The model estimates other needed driving variables such as solar (shortwave) and longwave radiation, relative humidity, vapor pressure, and vapor pressure deficit from the temperature and precipitation data (Hamlet et al. 2013). A dataset of historical meteorological parameters from January 1915 to December 2006 was developed at 0.0625° spatial resolution based on methods developed by Hamlet and Lettenmaier (2005) as described by Elsner et al. (2010). This observed dataset is used as input to the VIC model to produce simulations of historical streamflow, which shall be referred to as the OBS-VIC simulation.
In addition to these meteorological forcing variables, the VIC model uses inputs of four parameter files that are static and predefined: the soil parameter file, vegetation parameter file, vegetation library, and snowband parameter file. Each parameter file contains information specific to each model grid cell, while the vegetation library contains general information for vegetation types that may be referenced by the vegetation parameter file. These parameter files were developed for the 0.0625° VIC implementation (Elsner et al. 2010). The soil parameter file, which contains information such as soil layer depth and infiltration and drainage parameters, was used for model calibration (Hamlet et al. 2013).
The model incorporates topographically varying soil depths. The top and middle soil depths were fixed to 0.1 and 0.3 m, respectively, and the soil depth of the bottom VIC model soil layer was calculated using a digital elevation model (DEM) within a predefined range of 0.5–2.5 m (Elsner et al. 2010). The VIC model produces runoff and base flow at each model grid cell for each model time step. A separate streamflow routing model, developed by Lohmann et al. (1996), is used to generate daily time step streamflows at various river locations from the simulated runoff and base flow from each cell.
Calibration of the VIC hydrologic model at 0.0625° spatial resolution was carried out for monthly streamflow in 11 major subbasins of the CRB (Hamlet et al. 2013). The performance of the calibrated model was evaluated at smaller watershed scales where naturalized flow data were available. Model calibration and validation used a split sample approach in which calibration was performed over a 15-yr period (typically water years 1975–89) and model validation was performed over a separate 15-yr period (typically water years 1960–74). Calibration and validation periods were chosen to include a range of wet, dry, and average years to test and/or evaluate model performance. Figure 4 shows Nash–Sutcliffe efficiency (NSE) and correlation R2 statistics for each site for which observed natural streamflow was available. NSE indicates the ability of the model to predict observations, with values above 0 indicating skill superior to climatology and 1 indicating a perfect match with observations. Generally, a well-calibrated model produces NSE and R2 values that exceed 0.7 (Liang et al. 1996; Nijssen et al. 1997). A few locations showed negative NSE coefficients because of substantial bias in the simulations (to which the NSE metric is very sensitive). Most of these sites showed acceptable R2 statistics, however, which supports the assumption that relative changes (e.g., the ratio of the future to historical 100-yr flood reported below) are useful metrics even for those sites with relatively poor NSE scores. Additional evaluation of simulated flood statistics using observed streamflow data is reported below.
A detailed analysis of the uncertainties resulting from hydrologic modeling is beyond the scope of this initial investigation; however, a brief discussion of these issues is in order since this study uses a single hydrologic model implementation to assess changes in flood risk. To begin with, previous studies (e.g., Wilby and Harris 2006) have demonstrated that in rain-dominant basins, uncertainties in simulated hydrologic extremes are dominated by uncertainties in precipitation inputs from driving climate model simulations. Thus, for lower elevation sites included in this study, uncertainties in simulated hydrologic extremes are mostly related to uncertainties in simulated storm statistics from WRF and the large-scale forcing scenario. In basins with snow as an important part of the hydrologic cycle, however, important uncertainties can be introduced by the simulated temperature response of snowpack in the model. In particular, differences in snow modeling approach (e.g., energy balance versus temperature based) can result in different sensitivity to increasing temperature, and the resulting uncertainty in antecedent snowpack and soil moisture ultimately influence extreme high flows in the simulations. Thus, for colder basins, hydrologic modeling uncertainty may be important in assessing changing hydrologic extremes, and studies like this one that use a single hydrologic model may ultimately underrepresent the uncertainty in quantitative outcomes. These issues can be important in understanding the potential range of outcomes for a particular river location. In this paper, however, we focus on the broad response of flood risk to changing conditions over a large number of river sites with different historical temperature and precipitation regimes, and these issues are less important to the conclusions we make. Furthermore, the metrics used to assess change in peak flows (ratios of the future to historical 100-yr flood) are, by design, relatively insensitive to model calibration choices and other sources of hydrologic modeling uncertainty.
3. Bias correction and downscaling of 12-km WRF simulations
Although regional-scale climate models represent the important topographic features of the PNW and the mesoscale structure of storms that control flooding in PNW rivers better than global models, RCM simulations are still subject to biases resulting from deficiencies in both the global forcing fields and the regional model (Wood et al. 2004; Christensen et al. 2008). To obtain acceptable hydrologic simulations, these biases must be removed when using RCM results in impacts studies. In addition, to link the WRF results to the VIC hydrologic model, the simulations require additional downscaling from the 12-km WRF grid to the 0.0625° VIC grid. These downscaling and bias correction steps are performed as follows. The daily WRF data were mapped to the 0.0625° grid using the Synagraphic Mapping System (SYMAP) algorithm [Shepard (1984), as applied by Maurer et al. (2002)], which uses a four-nearest-neighbor, inverse-distance weighted average approach. Regridded precipitation and temperature data were then bias corrected using a quantile mapping approach (Wood et al. 2002; Oettli et al. 2011) applied at daily time steps to the WRF output data. This approach maps between simulated and observed cumulative distribution functions (CDFs) to remove bias, while largely preserving local climate signals and time series behavior in the simulation. These steps remove systematic bias in simulated meteorological variables resulting from combined bias in the large-scale forcing and WRF simulations, as well as the bias introduced by obtaining daily maximum and minimum temperatures from simulated 6-hourly data. Simulated wind speed data were not bias corrected; previous 0.0625° wind speed datasets used for VIC simulations were derived from regridded reanalysis data (Elsner et al. 2010), and the 12-km WRF data used here are believed to be superior because the effects of local topography are explicitly represented.
For the bias correction, daily CDFs were generated for observed and simulated temperature and precipitation at each grid point and each calendar month. The observed CDFs were developed using a training period of 1970–99 from the gridded historical temperature and precipitation data (described above), and model CDFs were obtained from the regridded ECHAM5-WRF simulation for the same period. For a given time step in the WRF time series, the quantile mapping algorithm looks up the simulated value in the model CDF from the training period. The corresponding quantile value from the observed CDF then replaces the simulated value for that time step.
For the historical period (1970–99), the simulated values fall within the range of the model CDFs by construction, and the mapping between CDFs is a simple lookup procedure. For future periods, systematic changes in both the mean and variance result in simulated data falling outside of the original training period CDF. In the case of temperature, the projected warming trend in the transient ECHAM5-WRF simulation causes a large displacement relative to the training CDF. Thus, for each 30-yr future period, the change in average temperature, relative to the training period, is removed before bias correction and added back to the bias-corrected value. Even after removing the mean temperature change, future simulated temperature or precipitation values frequently fall outside the range of the CDF derived from the training period, requiring extrapolation of the CDF. In general, this extrapolation could be accomplished by fitting a specific probability distribution to the observed and ECHAM5-WRF CDFs and then mapping between quantiles. Because of the relatively small samples used in this study, however, the tails of the distributions are not well resolved, and this approach frequently produced highly unrealistic results in practice. As a more robust alternative, values that fall outside the model CDF were mapped to an equivalent observed value via the multiplicative (P) or additive (T) anomaly from the mean. Note that these adjustments are equivalent to correcting the ECHAM5-WRF values for the multiplicative or additive bias in the mean.
The downscaled and bias-corrected WRF simulations are analogous to high-resolution “weather forecasts” for historical and projected future conditions. By construction, the downscaling process largely preserves the seasonal timing, spatial location, size, intensity, and sequence of events from the WRF simulations, while removing systematic biases in comparison with observations. Similarly, the simulated temporal correlation between temperature and precipitation are preserved by the downscaling procedure, an important consideration in accurately reproducing the timing and intensity of flooding events in mountain environments. For climate change scenarios, relative changes in key meteorological parameters are preserved by the bias correction.
4. Simulated flood statistics
Using bias-corrected data from the WRF simulations as input, the VIC hydrologic model simulates natural streamflows (i.e., streamflows that would occur without the influence of dams or diversions) at a daily time step, which will be referred to as ECHAM5-WRF-VIC simulations. Streamflows from ECHAM5-WRF-VIC simulations were generated for the three 30-yr periods: 1970–99, 2010–39, and 2040–69. These periods nominally represent the statistics for the decades of the 1980s, 2020s, and 2050s. Based on simulated daily flows, the annual daily peak flows and peak flow dates (the day when the peak flow occurred, starting from 1 October) were extracted for each water year for each of the 297 experimental watersheds.
Based on the annual peak daily streamflows (30 realizations in each sample period), we estimated the 100-yr flood using the generalized extreme value (GEV) distribution and L moments (Wang 1997; Hosking and Wallis 1997; Tohver et al. 2014). The 100-yr return interval flood was estimated from the fitted GEV distributions for each river basin and for the three 30-yr periods. The choice of the GEV distribution was based in part on studies showing that the performance of this distribution was superior when the true nature of the probability distribution is unknown, as here (Potter and Lettenmaier 1990; Hamlet and Lettenmaier 2007). The fit of the GEV distribution to the annual extremes was examined in spot checks of basins of different types and was found to be a close match with the simulated data in general, although issues related to the relatively small sample size in fitting to the extreme tails of the CDF were sometimes noted.
To evaluate the model performance, we compared streamflow from the OBS-VIC and bias-corrected ECHAM5-WRF-VIC simulations for 1970–99 to observed streamflow data from the Hydro-Climatic Data Network (HCDN; Vogel and Sankarasubramanian 2005). The performance of OBS-VIC results indicates how well the VIC model can simulate the observed hydrologic extremes when driven by gridded observations. Because the ECHAM5-WRF bias correction is based on the gridded observational temperature and precipitation dataset, by construction, the bias-corrected ECHAM5-WRF-VIC results should return similar statistics as the OBS-VIC simulation during the 1980s. Results are presented in the next section.
After validating the ECHAM5-WRF-VIC simulations, the impacts of climate change on flood risk were evaluated by comparing the 1970–99 ECHAM5-WRF-VIC simulations with future projections for the 2010–39 and 2040–69 periods. In contrast to statistically downscaled global climate models, the WRF Model is expected to simulate more realistic storms through an improved representation of terrain effects and mesoscale dynamics, which are important for projecting localized extremes in precipitation, temperature, and wind speed in a changing climate. To compare the new results with results from previous studies, we also compare the ECHAM5-WRF-VIC results with the same ECHAM5 simulation statistically downscaled using the hybrid delta (HD) method (ECHAM5-HD-VIC; Tohver et al. 2014; Hamlet et al. 2013). Additionally, streamflows using raw data from the ECHAM5-WRF simulation without bias correction are used to evaluate the bias correction process on simulated flood statistics. We expect relative changes in flood risks in the simulations without bias correction should broadly match those in the bias-corrected simulations.
5. Results and discussion
a. Evaluation of simulated peak streamflows
Results from OBS-VIC and 1970–99 ECHAM5-WRF-VIC simulations are compared with HCDN observed streamflow data for three representative watershed types: 1) a snowmelt-dominant watershed, the Similkameen River near Nighthawk (Fig. 5); 2) a transient, mixed rain–snow watershed, the White River below Tygh Valley (Fig. 6); and 3) a rain-dominant watershed, the Chehalis River near Grand Mound (Fig. 7). Note that the 1970–99 ECHAM5-WRF-VIC and HCDN data are only compared through CDF plots and date–magnitude peak flow scatterplots because the ECHAM5 time series is not expected to reproduce the historical timing of events, only the long-term statistics.
As expected, after bias correction, the CDF plots for peak streamflows from the ECHAM5-WRF-VIC simulations are very similar to CDFs from the OBS-VIC simulation (see lines marked by triangles in Figs. 5–7, left). This result indicates that the bias correction applied to simulated temperature and precipitation substantially improves bias in simulated peak streamflows. Flow magnitudes depend also on antecedent conditions resulting from the timing and sequence of weather events. Thus, simply reproducing the climatological CDFs of temperature and precipitation does not assure simulation of specific historic events and duplication of the CDFs of streamflow. Nevertheless, both ECHAM5-WRF-VIC and OBS-VIC simulations reproduce reasonably well the magnitudes and seasonal timing of the observed annual daily peak flows based on HCDN data. For the snow-dominant basin (Similkameen River near Nighthawk; Fig. 5), the flow magnitude with lowest exceedance probability (i.e., the highest flow event) from the OBS-VIC simulation is not represented in ECHAM5-WRF-VIC. Extreme flow events in this basin occur in the melt season and depend on concurrent antecedent snowpack, warm temperatures, and rainfall. This particular WRF simulation did not represent the sequence of events leading to the most extreme historical event, despite the bias correction. For the Chehalis River near Grand Mound, OBS-VIC and ECHAM5-WRF-VIC with bias correction simulations show a substantial low bias for flow magnitude (Fig. 7). Note that the bias occurs between OBS-VIC and the observed streamflow, while the ECHAM5-WRF-VIC simulation essentially matches the observation-driven VIC simulation. This bias is likely due to a deficiency in the gridded observed precipitation dataset that might be caused by an erroneous estimate of long-term mean precipitation at high elevation, the use of primarily low-elevation station data to provide the time series of precipitation, or spatial extrapolation across data gaps in complex terrain. Because of the deficiencies in the observed dataset, the ECHAM5-WRF-VIC without bias correction (Fig. 7, top) matches the HCDN data better than the bias-corrected simulation.
Figure 5–7 (bottom) show the HCDN observed and the OBS-VIC simulated annual peak flows for the water years 1971–88. In each case, the interannual variability is well represented in the peak streamflows derived from gridded temperature and precipitation observations. Correlations between the observed and simulated annual maximum flows are 0.91 in the snowmelt-dominant basin (Fig. 5), 0.72 in the transient rain–snow basin (Fig. 6), and 0.95 in the rain-dominant basin (Fig. 7). Note that, despite the deficiency in simulated peak flows for the rain-dominant basin, the temporal variability is well represented.
b. Projections of future flood risks
Figures 8 and 9 compare simulated current (1970–99) and future (2010–39 and 2040–69) streamflow statistics for six representative watersheds: three snowmelt-dominant watersheds (the Columbia River at Revelstoke Dam, the Kootenai River at Corra Linn Dam, and the Boise River at Boise; Fig. 8), a transient rain–snow watershed (the Skagit River near Mount Vernon; Fig. 9, top), and two rain-dominant watersheds (the Nisqually River at Alder Dam and the Chehalis River at Porter; Fig. 9, middle and bottom). To illustrate the coincident changes in the magnitude and timing of peak flows, the data are presented as scatterplots of annual peak flow date versus annual peak flow magnitude for each water year of the bias-corrected ECHAM5-WRF-VIC simulations, with results for the early twenty-first century (2010–39) on the left and mid-twenty-first century on the right (2040–69). The median date for peak flows is indicated by a vertical line. The direction and magnitude of changes in flood risk and peak flow date depend on a watershed’s location, elevation, and aspect. In very cold, strongly snowmelt-dominant watersheds, climate change is not expected to cause major timing shifts in annual peak flows, which will continue to occur because of spring snowmelt (Tohver et al. 2014). Changes in the peak flow magnitude may change, however, depending on projected changes in precipitation and temperature. The Kootenai River at Corra Linn Dam is a good example of a cold, strongly snowmelt-dominant watershed (Fig. 8, middle). For 2040–69, the simulation yields a modest shift toward earlier spring flooding and insignificant change in flood magnitude. By comparison, annual peak flows are projected to increase moderately for the Columbia River at Revelstoke Dam (Fig. 8, top) because of increasing winter and spring precipitation. For a warmer snowmelt-dominant watershed such as the Boise River at Boise, climate change is projected to cause a larger seasonal timing shift in annual peak flows and a greater increase in the magnitude of peak flows (Fig. 8, bottom) because of the greater sensitivity of the snowpack and freezing elevation to increased temperatures.
Transient rain–snow watersheds are more sensitive to warming because, as the freezing level rises with increased temperature, an increased basin area is subjected to rainfall and runoff production during storm events, especially in the early winter. Storms producing annual peak flows are projected to occur earlier in the water year in the ECHAM5-WRF-VIC simulations of future climate change, with more severe and earlier peak flows projected for transient rain–snow watersheds. The Skagit River near Mount Vernon is a good example of a transient rain–snow watershed (Fig. 9, top). For this basin, the seasonal timing of peak flows is shifted earlier in the water year, increasing October/November daily peak flows from about 10% of the total for the 1970–99 runs to about 30% for 2040–69. These findings are broadly consistent with previous studies evaluating observed changes in hydrologic conditions and hydrologic extremes in moderate elevation rain–snow basins in the PNW (Hamlet et al. 2005, 2007; Hamlet and Lettenmaier 2007; Elsner et al. 2010).
Results for the Nisqually River at Alder dam and the Chehalis River at Porter (Fig. 9) illustrate simulated changes in the annual peak flow date for relatively warm, rain-dominant watersheds. Because they lack substantial snowpack, rain-dominant watersheds are not very sensitive to warming (Elsner et al. 2010; Tohver et al. 2014), and the simulated changes in flood intensity are mostly due to projected changes in precipitation. When examining the projections of a single global model as we do here, the stochastic effects of decadal variability in precipitation become more apparent in the rain-dominant basins. Thus, while we might expect steadily intensifying flood risks projected from an ensemble of climate model projections, the ECHAM5-WRF-VIC simulation shows greater increases in annual peak flows for 2010–39 than for 2040–69. Although peak flows in both future time periods are higher than in the observed period (1970–99), decade-to-decade precipitation variability is strongly influencing the results as well. It is important to note that the sequence of events in any one model run do not match the time series observed in the past or that will actually occur in the future. Thus, the results should not be interpreted as a prediction that the flood risk for 2010–39 will be higher than for 2040–69, but rather that this is one possible sequencing of events that could occur within the combined influence of systematic climate change and decadal variability.
Following approaches developed by Hamlet and Lettenmaier (2007) and Tohver et al. (2014), the ratios of projected future 100-yr flood to the simulated contemporary (1970–99) 100-yr flood for 297 sites are plotted against average historical temperatures for December–February (DJF) for each basin in Fig. 10. The colors of dots for each station indicate the typical month of historical flood occurrence so that results may be differentiated according to basin characteristics. Note that values greater than one indicate an increase in flood risk and values less than one indicate a decrease. Overall, the bias-corrected ECHAM5-WRF-VIC simulations show increased flood risk for 2040–69 relative to 1970–99. Warm, rain-dominant basins with DJF temperatures larger than about 2°C and peak runoff early in the water year show exclusively increased flood risk while cold basins (DJF temperatures less than about −6°C, late peak runoff) show reduced flood risk at many sites (Fig. 10, top).
To illustrate the effects of different downscaling approaches, we compare 1) ECHAM5-WRF-VIC simulations with bias correction applied to WRF (Fig. 10, top); 2) ECHAM5-WRF-VIC simulations without bias correction (Fig. 10, middle); and 3) ECHAM5-HD-VIC, described in section 4 (Fig. 10, bottom). The results with and without bias correction are similar, but the bias-corrected results indicate more consistent increases in flood risk, especially for warm, rain-dominant basins (DJF temperatures >2°C). In contrast, the results from the hybrid delta statistical downscaling method (Fig. 10, bottom) are quite different from the WRF-based results. The statistical method is based on monthly time step data from the global model with daily temperature and precipitation constructed from historical time series (see Tohver et al. 2014). Hence, the statistical downscaling result has only a few degrees of freedom and does not represent the full range of changing storm and event-scale characteristics of the climate change projection in the detailed manner provided by the regional climate model. In particular, for the warmest, rain-dominant basins, the statistical downscaling results show modest changes in flood risk while the WRF results show much more substantial increases. More intense future storms simulated by WRF in relatively warm areas of the coastal lowlands and west slopes of the Cascades (see Fig. 2) contribute to the increased flood risk shown for these basins. The direct effects of temperature on snowpack, which are well represented by the statistical downscaling approach, do not have a substantial impact on flood risk in these warm basins. Likewise, many of the cold, high-elevation basins also do not show large changes in flood risk as they do not warm sufficiently to substantially alter snowpack dynamics. For the statistical method, it is only the transient rain–snow basins that show substantial changes in flood risk as warmer temperatures result in more rain-dominant events, increased contributing basin area, and less buffering of streamflow by deep snowpack.
The inclusion of realistic mesoscale precipitation processes in the regional climate model simulations results in a more complex projected response of flood risk to climate change across the region. The results from WRF (Fig. 10, top) show considerable variability from basin to basin, reflecting the sensitivity of flooding to highly variable mesoscale effects. Furthermore, differences in the circulation patterns associated with atmospheric river events can have a substantial influence on the watersheds experiencing high streamflows because of the complex orography of the region (Neiman et al. 2011). Since the diversity of potential changes in daily weather and associated orographic effects are more fully incorporated into the WRF simulations than the statistical downscaling, considerably more scatter is found in the climate change projections from WRF. This result, while providing a less clear-cut depiction of climate change effects, arguably better reflects the relative magnitudes of variability, uncertainty, and systematic trends (i.e., signal to noise ratio) in determining future changes in flood risk.
The sensitivity of these results to particular storm characteristics and decadal-scale precipitation variations points to the need to analyze a multimodel ensemble before drawing robust conclusions. Nevertheless, the results using a single RCM scenario support the hypothesis that changes in flood risks in individual basins depend on mesoscale processes and storm characteristics that are likely not well represented by statistical downscaling methods. Furthermore, future increases in heavy precipitation are consistently projected by many climate models, and the ECHAM5 model is not unique in this regard (Tebaldi et al. 2006). Our results suggest that in rain-dominant watersheds, flood risk is likely to increase much more than has been shown in earlier studies using GCM projections downscaled by the hybrid delta downscaling and similar techniques (i.e., the hybrid delta projections are very conservative for rain-dominant watersheds). The ECHAM5-WRF-VIC simulations are also much noisier than the hybrid delta runs for two reasons: 1) the hybrid delta results are based on 91 years of data in each period, which results in less noise in the GEV fits to the extreme data, and 2) by construction, the hybrid delta method replicates the same cycles of interannual variability for both the historic and future periods, which underrepresents the true number of degrees of freedom associated with potential changes in natural variability.
Initial projections of changing flood risk in the PNW from past studies using statistical downscaled GCM data suggest increasing flood risk in most areas of the region due to projected regional warming and increases in cool season precipitation (Tohver et al. 2014). Regional climate models based on WRF simulations offer a new research direction for improving our understanding of the spatial variations in changing hydrologic extremes across the region. These models can provide a more physically based representation of features such as atmospheric rivers and orographic precipitation that determine the local changes in heavy precipitation and the timing of flooding events. Our results using ECHAM5 A1B global climate change simulations as the large-scale forcing for WRF suggest more extreme storms in the early fall and general increases in flood intensity in the PNW by 2040–69. In particular, many sites show shifts toward flooding earlier in the water year and increased flood risk due to the combination of reduced snowpack and more intense precipitation events.
In contrast to changes in flooding estimated using statistically downscaled ECHAM5 GCM data, the WRF scenarios show more substantial increases and more variability in flood intensity, particularly in rain-dominant and mixed rain–snow basins, and an increased range of uncertainty in projected flood risk. These results follow from the more detailed representation of important storm characteristics such as position, size, and intensity of extreme precipitation that are related to flooding. Future changes in these storm characteristics, which are not well represented by statistically downscaling the global model, can have important implications for extreme streamflow in individual basins and can create additional pathways by which climate change can affect future flood risk at the local scale. Given that global climate model projections indicate a general widening of the tropical belt (Seidel et al. 2008) with a large-scale northerly shift in the midlatitude storm tracks (Ulbrich et al. 2009), it is plausible that the typical large-scale configuration of atmospheric rivers would be perturbed in the future climate.
The results of this study suggest that future changes in extreme weather systems may have a substantial effect on flood risk. These changes are likely a combination both of thermodynamic processes, which result in more moisture transport and heavier precipitation, and dynamic processes, which alter the terrain effects that control the local distribution of heavy precipitation. This added complexity due to mesoscale processes underscores the need for comprehensive multimodel ensemble analysis and is a topic of ongoing research using multiple regional climate realizations. In particular, a better understanding of the spatial and temporal variability of flood risk in both the current and future climate is essential to our ability to project the timing for the emergence of significant changes in the risk of flooding across the region.
This research was supported through funding from the U.S. Army Corps of Engineers. Dr. Hamlet also received partial support for this project from a grant from the National Science Foundation (Grant EAR-0838166). We thank Robert Norheim of the UW Climate Impacts Group for producing Figs. 1, 3, and 4.