Current understanding of the combined effects of sea level rise (SLR), storm surge, and changes in river flooding on near-coastal environments is very limited. This project uses a suite of numerical models to examine the combined effects of projected future climate change on flooding in the Skagit floodplain and estuary. Statistically and dynamically downscaled global climate model scenarios from the ECHAM-5 GCM were used as the climate forcings. Unregulated daily river flows were simulated using the VIC hydrology model, and regulated river flows were simulated using the SkagitSim reservoir operations model. Daily tidal anomalies (TA) were calculated using a regression approach based on ENSO and atmospheric pressure forcing simulated by the WRF regional climate model. A 2-D hydrodynamic model was used to estimate water surface elevations in the Skagit floodplain using resampled hourly hydrographs keyed to regulated daily flood flows produced by the reservoir simulation model, and tide predictions adjusted for SLR and TA. Combining peak annual TA with projected sea level rise, the historical (1970–1999) 100-yr peak high water level is exceeded essentially every year by the 2050s. The combination of projected sea level rise and larger floods by the 2080s yields both increased flood inundation area ( 74%), and increased average water depth ( 25 cm) in the Skagit floodplain during a 100-year flood. Adding sea level rise to the historical FEMA 100-year flood resulted in a 35% increase in inundation area by the 2040's, compared to a 57% increase when both SLR and projected changes in river flow were combined.
There is a strong scientific consensus that human behavior is altering greenhouse gas concentrations and the global climate system, and projected changes in climate are expected to substantially impact human and natural systems in the twenty-first century (IPCC 2007). Projected changes in climate in the Pacific Northwest (PNW) and the Skagit River basin include losses of snowpack, changes in seasonal streamflow timing, increasingly intense hydrologic extremes (floods and low flows) and sea level rise (SLR) (Elsner et al. 2010, Hamlet et al. 2013, Lee et al. 2016, Tohver et al. 2014, Salathé et al. 2014, Mote et al. 2008, Nicholls and Cazenave 2010, National Research Council 2012). Natural climate variability also strongly influences PNW climate. The Pacific Decadal Oscillation (PDO) (Mantua et al., 1997) and El Niño Southern Oscillation (ENSO) are two climate phenomena that exert strong controls on seasonal air temperature and precipitation (Mote 2003, Mote et al. 2003), flood risk (Hamlet and Lettenmaier 2007), and sea level in the PNW.
Water levels in estuaries and floodplains can be particularly sensitive to changing climate. These areas are impacted from the marine side by tidal anomalies (TA) (storm surge) and SLR, and from the freshwater side by seasonal changes in river flow and hydrologic extremes. Understanding these changes is a crucial step in developing proactive strategies for the management of both the natural and human systems in the floodplain and near-coastal areas. This study quantifies future water levels in the Skagit River basin (Figure 1) caused by projected impacts on SLR, TA, riverine flooding, including the effects water resources infrastructure (dikes, levees, and storage dams). A principal goal of this research is to develop tools and approaches to make climate change projections that are useful for future planning in coastal and near-coastal areas.
Storms are important drivers of high estuarine water levels due to the combination of river flooding and coincident barometric pressure and wind effects on sea levels (Cayan et al. 2008). In this study the time series behavior of individual storms, simulated by a regional climate model (RCM), is kept intact within the TA and hydrology models. In this way, coincident TA and riverine flooding can be modeled as products of each individual storm, which facilitates a realistic pairing of SLR, tidal cycle, TA, and river flooding events in the projection of impacts to estuarine flooding.
The research described here has important practical implications for coastal and riverine floodplain management. Nearly all present-day floodplain management strategies and policies rely entirely on the historical records of TA and river flow, which is a questionable practice in a non-stationary climate (Milly et al. 2008). Using scenario-based approaches to manage floodplain environments, such as the one described in this paper, will likely result in more effective long-term planning and management strategies in response to a non-stationary and rapidly evolving environment.
Tidal Data—Observed hourly tidal data for Seattle (NOAA 9447130; 1899–2012), and Sneeoosh Point (NOAA 9448576; 05/2000-08/2000) were obtained from the National Oceanic and Atmospheric Administration's (NOAA) tide station network database (NOAA 2011).
Sneeoosh Point is the closest station to the Skagit estuary, but available data are extremely limited. Following the methodology used by the Federal Emergency Management Agency (FEMA) for the Skagit County Flood Insurance Study, the Seattle tidal time series was used as the historical observations. Although the observed time series at the two stations do not match perfectly in the overlapping period, the significantly longer period of record at the Seattle station and the geographic proximity of the study locations to Seattle justify the substitution.
Observed and Simulated River Flows—Mean daily flows for the Skagit River at Mount Vernon (USGS 12199000; 1908–2012) were obtained from the United States Geologic Survey (USGS) website (USGS 2011). Hourly flows for the Skagit River at Mount Vernon over a shorter period were used to disaggregate daily flood hydrographs for use in the hydrodynamic model.
Daily time step regulated flow scenarios based on statistically downscaled climate change scenarios (Hamlet et al. 2013) were simulated by Lee et al. (2016) for the Skagit River at Concrete and Mount Vernon using the SkagitSim reservoir operations model.
Sea Surface Temperatures and Climate Indices—Historical monthly values for the Niño 3.4 index were obtained from the NOAA's Climate Prediction Center (CPC 2011). The indices were calculated using the OISST.v3 gridded Sea Surface Temperature Data Set (Smith et al. 2008) between (5°N—5°S, 170°W—120°W) with a climatological base period of 1981–2010.
Climate Model Simulations—For this project the Coupled Model Intercomparison Project, Phase 3 (CMIP3) results from the ECHAM-5 global climate model (GCM), ensemble member one, forced by the A1B emissions scenario (Nakićenović et al. 2000), were used as the climate forcing for subsequent models and techniques. Model output was six hourly, which was also used as the forcing interval for the Weather Research and Forecasting (WRF) regional climate model (RCM) (See Salathé et al. 2014 for additional details on WRF implementation and model coupling).
In their evaluation of 10 climate models, Salathé et al. (2007) concluded that the ECHAM-5 model most accurately simulated the observed climate in the PNW relative to the other models. ECHAM-5 also ranked highly using updated criteria established by Mote and Salathé (2010). In particular, both studies found relatively small biases for ECHAM-5 relative to historical National Center for Atmospheric Research (NCAR) reanalysis data. The PNW temperature and precipitation changes predicted by the ECHAM-5 model also fall near the averages of the other global climate models (GCMs). Lastly, among the twenty-one coupled GCMs included in the CMIP3 and IPCC AR4 report, the ECHAM-5/MPI-OM was the only model to successfully replicate the observed variability of the ENSO pattern in the equatorial Pacific sea surface temperatures (Lin 2007), an important driver of interannual climate variability in the PNW. These qualities make the ECHAM-5 model a good selection for singlemodel consideration.
The A1B emissions scenario represents a “business as usual” progression of economic and policy responses, and technological development until the mid-twenty-first century, resulting in rapid increases in greenhouse gas concentrations and global temperatures in the first half of the 21st century, followed by more substantial greenhouse gas mitigation and reduced rates of warming in the second half of the century (IPCC 2007).
Overview of Numerical Models
Our study uses a suite of linked physically and empirically based numerical models to quantify changes in extreme high water levels expected in the Skagit River Basin under projected future conditions (Figure 2). Simulations from the ECHAM-5 A1B climate model scenario were used as the future climate forcing, and were dynamically downscaled using the WRF RCM (Salathé et al. 2014). Statistically downscaled results from previous studies using the Hybrid Delta approach (Hamlet et al. 2013, Tohver et al. 2014) were also used as a point of comparison and to perturb existing estimates of high-flow extremes used by FEMA in previous studies. Downscaled and bias-corrected outputs from the WRF simulations were used as the meteorological inputs for the Variable Infiltration Capacity (VIC) distributed hydrology model and a linear regression storm surge model based on a) ENSO forcing and b) pressure patterns simulated by WRF. Daily streamflows from the VIC simulations were used to force the SkagitSim reservoir operations model (Lee et al. 2016) to estimate regulated flows for both historical and future conditions. The final step was to create integrated scenarios by using disaggregated hourly flows and tides, plus daily TA, as inputs to a hydrodynamic model to determine the depth and spatial extent of inundation during flooding. The sum of hourly tide, daily tidal anomaly and base SLR were used as the lower river boundary condition for the hydrodynamic model in these simulations. The following sections outline the methods used in developing temporally consistent projections of streamflow, tides, and TA.
Climate Model Downscaling
Due to the low spatial and temporal resolution of GCMs, simulations of regional and local-scale climate have many shortcomings, especially in topographically complex regions such as the PNW. Mass et al. (2002, 2003), for example, show that a horizontal grid resolution of 15 km or less is necessary to successfully resolve orographic precipitation in the PNW. Thus in order to use contemporary GCM results for regional and local applications some form of downscaling and bias correction must be employed to recapture important local climate features in complex terrain.
Hybrid Delta Statistical Downscaling Approach—Hydrologic simulations based on the Hybrid Delta statistical downscaling method (Hamlet et al. 2013, Tohver et al. 2014) are used as one set of inputs in this study. This downscaling method produces realistic storms, and projections of future hydrologic extremes, based on observed storm patterns from the historical record. However, potential changes in the probability distributions of daily precipitation, seasonality, storm size, storm track position, and interarrival time of storms are not captured by this approach. We will refer to statistically downscaled climate scenarios for the 2040s and 2080s using this approach as ECHAM5-HD-2040S and ECHAM5-HD-2080S, respectively.
Dynamic Downscaling Using a Regional Climate Model—The dynamic downscaling technique used in this study is based on the WRF RCM implemented over the PNW at 12 km resolution (Salathé et al. 2010, 2014). For the historical run, this model was forced at the outer boundary by the NCAR Reanalysis (1970–1999), which is a largescale global simulation of the observed climate. ECHAM-5 GCM forced by the A1B emissions scenario (1970–1999, 2010–2039, and 2040-2069) provided a dynamically simulated historical baseline (“1970–1999”) and two future projections (2010–2039, and 2040–2069). We will refer to these climate simulations as ReAnal-WRF-1980s, ECHAM5-WRF-1980S, ECHAM5-WRF-2020s, ECHAM5-WRF-2050S, respectively.
In the context of understanding future hydrologic extremes, the use of a dynamic downscaling approach that better represents fine-scale topography (i.e. the Cascade Mountains) and atmospheric processes (i.e. orographic precipitation) provides a more realistic picture of climate change impacts in Western Washington. When forced by reanalysis data, and bias corrected, the model has been shown to successfully capture important storm characteristics related to flooding in the Skagit basin (Salathé et al. 2014). This approach maintains the explicit character of storms simulated by the 12-km WRF model (e.g. seasonality, location, size, intensity, interarrival time, etc.) by directly preserving the spatial and temporal patterns of simulated storms. Salathé et al. (2014) provide more details on the methods used in downscaling and bias correcting the WRF simulations.
The hydrologic model used to produce daily streamflows for this project was the Variable Infiltration Capacity (VIC) model (Liang et al. 1994) implemented over the PNW at 1/16th degree (Elsner et al. 2010, Hamlet et al. 2013). We will refer to historical VIC simulations based on gridded meteorological station observations as Hist-VIC simulations. Hybrid-Delta VIC simulations based on the ECHAM5 A1B scenario will be referred to as ECHAM5-HD-VIC simulations. These scenarios were extracted from simulations prepared by Lee et al. (2016) employing methods developed by Hamlet et al. (2013). VIC simulations based on downscaled and bias corrected WRF simulations were extracted from model runs carried out by Salathé et al. (2014). We will refer to these dynamically downscaled model runs as ECHAM5-WRF-VIC simulations.
Reservoir Operations Modeling
The SkagitSim reservoir operations model used for this study was developed by Lee et al. (2016) and adapted to accept the ECHAM5-WRF-VIC inflows. SkagitSim is implemented in the STELLA simulation modeling package (ISEE 2012) using reservoir simulation algorithms developed by Hamlet (1996) and further described by Hamlet and Lettenmaier (1999) and Lee et al. (2016).
Figure 3 shows the cumulative distribution functions (CDFs) of regulated peak annual flow at Mount Vernon from a) USGS daily streamflow observations, b) SkagitSim simulations forced by the Hist-VIC streamflow simulations, and c) SkagitSim simulations forced by ECHAM5-WRF-VIC-1980s streamflow simulations. Overall, the integrated model simulations do a good job of reproducing the observed probability distribution of regulated high flows at daily time step. Some of the bias in the CDFs in Figure 3 is due to bias in flow values inherited from the hydrology model, but also reflects some reservoir operations model bias in reproducing flood characteristics (See Lee et al. 2016 for additional model validation and discussion).
Sea Level Rise and Tidal Anomaly Modeling
Estuarine flooding can be caused by a combination of unusually high base sea levels (e.g. during warm ENSO years or as a result of climate change related SLR), tide levels, TA, or extreme river flow. A combined sea level and TA model (SLTAM) was developed using a least squares multiple linear regression approach to quantify the effects of sea level rise and TA on coastal flooding and backwater conditions that partly determine river stage in the tidally influenced zone. The following steps were taken to develop the model.
Deriving Tidal Anomalies—Using NOAA hourly water levels, a least squares approach was used to fit the principle harmonic tidal constants using the T Tide Matlab package (Pawlowicz et al. 2002). These constants were then used to predict hourly tides for the same time period (Figure 4). The daily mean tidal anomaly between the predicted and measured tides was calculated based on hourly data, and then binned by calendar month.
Developing Model Inputs—All inputs to the regression model were derived directly from the large-scale ECHAM5 GCM simulations or the ECHAM5-WRF output, with the exception of the historical Niño3.4 monthly time series, which was derived from the OISST.v3 gridded SST data (described above). Local surface pressure was extracted from ReAnal-WRF-1980 output at the grid cell centered over Seattle (48.3°N,-122.3°E). The simulated surface pressure time series was also aggregated from a 6-hour to 24-hour time step.
SVD Analysis—A singular value decomposition (SVD) analysis (also called Empirical Orthogonal Function (EOF) analysis) was performed on gridded standardized daily surface pressure anomalies from the ReAnal-WRF-1980 data (Figure 5). The principal components of the first and third EOFs, explaining 30.78% and 12.26% of the variance in the pressure field respectively, were used as explanatory variables in the TA regression model. For the climate change values, the gridded ECHAM5-WRF-2020 and ECHAM5-WRF-2050 data were projected onto the spatial EOFs and singular values of the reanalysis decomposition, yielding the principle component time series for future climate projections associated with EOF-1 and EOF-2.
ENSO Index Values—The historical monthly time series of the ENSO (Niño 3.4) index was also used to train the model (Figure 6B). For the future climate change scenarios for the 2020s and 2050s, the sea surface temperatures (SSTs) from the ECHAM5 GCM simulations were used to calculate the Niño 3.4 index value for each future month using methods described by Trenberth (1997). Figure 6A shows a scatter plot of the observed Nino 3.4 data and observed Seattle tidal anomalies. Figure 6B shows the observed time series of the Nino 3.4 Index (1950–2012) and the ECHAM5 simulated Nino 3.4 Index (1970–2069). Note that there is no expectation that the time series in the simulations should match observations, but the variability should be, and is, comparable.
Training the Tidal Anomaly Regression Model— An iterative approach was used to determine which variables best describe observed tidal anomalies. Local pressure (daily average and 3-day average), the 1st and 3rd principle component time series associated with the large-scale pressure SVD analysis, and the ENSO index were the only variables found to be statistically significant explanatory variables for the regression model. A separate regression model was constructed for each month to account for the seasonal nature of anomalies. Hence, the regression model was trained such that:
This approach yielded a good fit in winter (R2 > 0.75), which is the primary season of concern with regard to Skagit River flooding.
Predicting Future Tidal Anomalies—Using the regression parameters found for the observed training period and the input variables extracted from the ECHAM5 GCM and ECHAM5-WRF simulations, daily tidal anomalies were calculated for the ECHAM5 1970–2069 period.
Predicting Future Tides—An hourly time series of predicted future tides was created using the harmonic constants determined in step 1. The predicted hourly tides were then adjusted by adding the predicted daily tidal anomalies.
Sea Level Rise Scenarios—At the time of this writing, the science behind global and regional SLR projections is changing rapidly, and published projections of SLR have changed markedly even since the 2007 IPCC report. In particular, recent studies have pointed to higher estimates of global SLR and a greater range of uncertainty than those published in the 2007 IPCC report due to accelerated contributions from melting ice in Greenland and Antarctica and other factors (Cazenave et al. 2008). For example, Vermeer and Rahmstorf (2009) projected that global SLR could plausibly exceed 180 cm by 2100, far exceeding more conservative estimates associated with the IPCC AR4 (Cazenave et al. 2008). A recent summary report from the National Research Council (NRC 2012) projected mean global SLR for three time periods: 2030, 2050, and 2100 relative to observed global baseline sea levels for 2000. The NRC projections are approximately midway between the older AR4 projections and some of the more extreme values from the individual studies discussed above. Local estimates of SLR in Puget Sound show similar increases. For example, Mote et al. (2008) projected a median scenario of SLR for Puget Sound of 33 cm for the 2100 time frame based on the AR4 projections, whereas the NRC (2012) mean SRL projection for Seattle for this same time period is 62 cm (NRC 2012 Table 5.3). To construct quantitative scenarios we started with the NRC global projections (midpoints calculated from the range of values given in Table 5.2 from NRC 2012), fitted a second order polynomial to these primary data (shaded rows in Table 1), and then interpolated between the primary values to obtain values for 2025, 2035, 2045, 2055, 2085. Table 1 summarizes these calculations. Our intent was to produce SLR scenarios that would plot in the center of the range of the global values cited in the NRC report. The scenarios of sea level rise shown in Table 1 were then added to 2000 baseline sea levels, to which the hourly tidal cycle and simulated daily TA were added for each storm event.
Mean global SLR projections and interpolated values for future years based on a fitted second order polynomial. Shaded values are the central values from the ranges given in the NRC (2012) report (Table 5.2). Approach 1 value is from NRC (2012) (average of ranges, Table 5.2); approach 2 value is calculated based on best second order polynomial fit to NRC 2012 values.
It is important to note that these SLR scenarios are only intended to represent plausible values based on current projections of global SLR, and should not be confused with “predictions” of future conditions on specific dates, which for the reasons discussed above remain very uncertain. Based on the rapidly evolving scientific understanding of global ice dynamics, steadily improving measurements of local-scale vertical land motion, and the resulting relative SLR, future projections of SLR will undoubtedly continue to evolve, and could be substantially higher or lower than the values we have used here.
The NRC report (2012) also provided local estimates of SLR specifically for Seattle (in Table 5.3), however we were reluctant to use these numbers because they included positive vertical uplift estimates (i.e. reduced relative SLR numbers) for all areas north of Cape Mendocino, when in fact there are large differences between local sites in Puget Sound and on the Pacific Coast. GPS-based elevation records for the Skagit lowlands, for example, show no statistically significant uplift, whereas coastal GPS readings off the Olympic Peninsula show strong positive trends in land elevation (see Mote et al. 2008 for more discussion). In addition, given the very large uncertainties in global SLR estimates, any potential advantage associated with using a more detailed local SLR projection seemed minimal at this time.
Integrated Scenarios Pairing Tidal Anomalies with Flood Flows
For the statistically downscaled scenarios (based on ECHAM5-HD-2040S, ECHAM5-HD-2080s scenarios, which have essentially the same time series as observations), storm hydrographs are paired with SLTAM simulations forced by observed Nino 3.4 Index values and pressure patterns from the ReAnal-WRF-1980 simulations. This decision was based on the TA results (discussed later in the paper—See Figure 10) in which minimal changes in the future TA probability distributions were found. Using this approach, temporal consistency is maintained between the hydrologic model inputs and the TA model inputs.
It is worth noting that the SLTAM does not explicitly include the effects of wind-driven storm surge, changes in which could also affect future coastal impacts. Wind effects, however, do not appear to be the dominant factor controlling observed TA in winter since the SLTAM model explains about 75% of the variance without including this factor.
Hourly Disaggregation of Daily Flood Hydrographs
The hydrodynamic model used in the study (discussed below) required hourly flow inputs in order to achieve computational stability, which required the temporal disaggregation of simulated daily streamflows to hourly time step values. The Steepness Index Unit Volume Flood Hydrograph Approach for Sub-Daily Flow Disaggregation (Tan et al. 2007) was used to map daily flood flows to unit hourly hydrographs based on the steepness of the rising limb of the flood event. The assigned hydrograph is then scaled to match the four-day rising limb flood volume of the modeled daily flood. This method was chosen because it does not prescribe a certain hydrograph shape to all future flood events, rather it finds the closest approximation based on the characteristics of the simulated daily flood hydrograph, which could be substantially different in the future climate simulations. In testing the approach on observed Skagit River floods, a reasonably accurate fit between observed and disaggregated hydrographs was achieved (Figure 7A). Regulated daily flows at Mount Vernon simulated by the SkagitSim reservoir model were similarly disaggregated to hourly hydrographs for all peak flow events that exceeded a 1275 cms threshold (roughly equivalent to the mean annual flood).
To facilitate a more direct comparison between the floodplain impacts associated with the FEMA 100-yr flood (FEMA Q100) and future projections of Q100 using different downscaling methods, the FEMA Q100hourly hydrograph was simply scaled to match the fractional change in the Q100 values for each future time period (Figure 7B)(Tohver et al. 2014). For example:
Table 4 shows the adjustment factors for each future time period.
A 2-dimensional (2D) hydrodynamic model (Flo-2D) was implemented downstream of river mile 22.3 (below Sedro Wooley, WA). The Skagit model was developed as a combined effort between U.S. Army Corps of Engineers (USACE) and the Federal Emergency Management Agency (FEMA) (USACE 2007, 2008, 2009). The model was constructed with a grid cell size of 400 ft × 400 ft (121.9 m × 121.9 m). Input streamflow data were delivered at an hourly time step and the computation time step was allowed to vary between 0.1 and 10 seconds. Due to a wide range of levee conditions in the Skagit Basin, the floodplain model included seven levee failure scenarios used for composite flood risk mapping. Composite inundation maps were created by finding the maximum depth at each grid cell over the seven levee failure scenarios for specific zones as outlined in the USACE hydraulic reports (USACE 2007, 2008, 2009). The majority of the hydrodynamic model runs in this study were completed only for the “No Levee Failure” scenario, although the scaled USACE/FEMA 2040s and 2080s floods were also evaluated for all seven levee failure scenarios. The benefits of using the same model combination used by USACE/FEMA are a) that the results here can easily be compared to recent studies in the Skagit River and b) the methods used to develop model inputs will be compatible with future studies using similar tools.
It should be noted here that in the course of the USACE/FEMA modeling study the Q100 peak daily flood magnitude was adjusted based on a variety of historical factors. For this reason, the magnitude of the FEMA Q100 event (5772 m3/s) does not match the one shown in Table 3 based on a GEV distribution fit to the USGS historical record from 1970–1999 (6544 m3/s). For simplicity and in order to provide more direct comparisons with previous work, the scaled FEMA Q100 was not pre-adjusted to account for this difference. Instead, the scaling factor was directly applied to the FEMA Q100 value.
Integrated Hydrodynamic Scenarios
The key components of ten flooding scenarios are described in Table 2, and additional descriptions of each of the model runs are provided here.
The Hist- 1980s scenario establishes a simulated baseline from 1970–1999. Gridded meteorological observations were used to drive the hydrology model (Hamlet et al. 2013, Hamlet and Lettenmaier 2005), and ReAnal-WRF-1980 simulations drove the SLTAM. Sea level rise was assumed to be zero (i.e. the 2000 baseline condition) for this time period.
Two statistically downscaled scenarios (HD-2040s and HD-2080s) were used to perturb the FEMA flood. ReAnal-WRF-1980 simulations drove the TA portions of the SLTAM, which were combined with mean SLR projections for the 2040s and 2080s respectively.
The three dynamically downscaled scenarios were split into three 30-year time blocks (WRF-1980s, WRF-2020S, WRF-2050s). For each of these three time periods, all meteorological forcings came from the downscaled ECHAM5-WRF simulations that were used as inputs to the VIC model and SLTAM. VIC streamflow simulations were then used as inputs to the SkagitSim model to provide regulated daily flow, and these were disaggregated to hourly time step as discussed above. Sea level rise was uniformly added to the tidal anomaly simulations for each future time period according to the projected value at the center of each future period (i.e. 2035, 2045, 2055, 2085—see Table 1). The total TA was then added to the hourly tidal cycle to estimate the lower boundary conditions for the hydrodynamic model and its alignment with hourly flood hydrographs.
Description of integrated scenarios and associated model runs. A total of eight 30-year periods were analyzed using a combination of Hybrid Delta statistical downscaling (HD) and dynamic downscaling techniques using WRF. In addition, four scenarios based on scaled FEMA estimates of the historical 100-yr flood and SLR projections were carried out.
Changes in Regulated Peak Flows
A number of previous climate change impacts studies have shown that projected warming and precipitation change in the PNW will lead to larger flood events on the west slopes of the Cascades (Hamlet and Lettenmaier 2007, Hamlet et al. 2013; Lee et al. 2016, Tohver et al. 2014), and dynamic simulations using RCMs show larger increases in extreme precipitation events on the windward slopes of the Cascades than global scale projections (Salathé et al. 2010, Dulière et al. 2011, Salathé et al. 2014). These local-scale effects on precipitation extremes are included in the ECHAM-WRF-VIC streamflow simulations.
Lee et al. (2016) showed that regulation reduces the absolute value of peak flows in comparison with unregulated flows, and changes in regulated flows relative to historical baselines are similar to those for unregulated conditions, except projected percent changes in peak monthly flows are substantially larger under regulated flows than changes in unregulated conditions. Figure 8 shows historical and future cumulative distribution functions (CDFs) of peak daily regulated flow for the Skagit River at Mount Vernon, and Table 3 shows the peak flow values for three probability of exceedance values (0.2, 0.1, 0.01). The CDFs for the HD-2040S and WRF-2050s scenarios are similar, showing only small differences between the two downscaling approaches for the Skagit. Lee et al. (2016) also showed that because there are major unregulated tributaries which contribute to flooding in the lower basin, (particularly the unregulated Sauk River), modifying flood rule curves to create more flood storage provides relatively little buffering against predicted increases in flood risk.
5-, 10-, and 100-year regulated events for the Skagit River at Mount Vernon for each model run calculated using the generalized extreme value (GEV) distribution. The fractions in parentheses are the increases relative to historical period of that model (i.e. WRF-2050s is relative to the WRF-1980s). Grey highlighting indicates poor fit for the GEVD compared to Cunnane (1978) unbiased quantile estimates shown in Figure 8.
Sea Level and Storm Surge
In an earlier study, Stammer and Hüttemann (2008) argued that climate change will cause relatively little change, if any, in pressure-induced TA. Our results support this hypothesis. Figure 9A shows the maximum annual mean daily TA at Seattle for the four WRF model runs: ReAnal-WRF-1980s, ECHAM5-WRF-1980S, ECHAM5-WRF-2020s, ECHAM5-WRF-2050s. Three observations are immediately apparent from this figure, 1) there are only minor variations between the three ECHAM5-WRF model runs indicating no systemic shift in TA due to climate change, 2) there is a small positive bias in the ECHAM5-WRF-1980s runs for the largest events relative to the ReAnal-WRF- 1980s run, indicating that the nested ECHAM5-WRF model may produce somewhat more intense storms than have historically been seen, and 3) there is a relatively narrow range of tidal anomalies between the 1-yr and 100-yr events (about 45.7 cm). The PNW's lack of large tidal anomalies combined with little change in the TA probability distribution supports the argument that SLR, rather than storm surge, will more significantly alter the future probability distributions of water levels in Puget Sound (Tebaldi et al. 2012).
The relative importance of sea level rise impacts on tidal anomalies are clearly seen in Figure 9B, which adds the CDF of extreme TA to the SLR projection. In this figure, the ∼ 1-yr event from the ECHAM5-WRF-2050s simulation exceeds the ∼ 100-yr event of the ECHAM5-WRF-1980s simulation, meaning that the current-climate 100-year water level event is projected to be exceeded every year by the 2050s. Similar changes in the probability distributions are found when the hourly tidal signal is included in the analysis. In Puget Sound, the range of the tidal cycle is nearly an order of magnitude larger than the observed storm surge, and adding 60–90 cm of sea level rise to typical high tides (even without storm surge) is likely to dramatically increase extreme high water levels in the Puget Sound lowlands.
Coincident Occurrence of Flooding and Storm Surge
Combining temporally consistent riverine flooding and water levels in Puget Sound provides the ability to assess the likelihood of coincident events as they vary in time. This is especially important when considering the estuary and tidally influenced portions of the river channel. In some Puget Sound rivers, for example, backwater conditions created by TA or extreme tides can affect river flood stage as far as 20 km upstream (see e.g. Figure 11G, below).
Figure 10 shows the coincident relationship between regulated peak annual streamflow in each year and the paired tidal anomaly. Although there is a wide range of response shown in the scatter plot, numerous large flood events in the future model runs occurred during periods of depressed sea levels, an association not apparent in the historical run. Attributing a cause to this observation is difficult but may be the result of increased influence of the ENSO phenomenon on flooding in the future simulations (e.g. more large floods in La Nina years which are associated with depressed sea levels) or altered storm dynamics (e.g., more frequent high pressure systems following large low pressure systems).
A positive trend relating these two variables is found in each model run (i.e. high TA tends to be paired with high flow), but the trend is not statistically significant at the 95% confidence level. This suggests that the TA during peak flow events is essentially decoupled and primarily reflects a random pairing of events. Although streamflow is systematically higher in the future scenarios, any relationship between the two variables remains weak in all cases. This is an important result, because it demonstrates that base SLR, which will dramatically increase water levels in most if not all events in the future, is likely to be the dominant factor in comparison with tidal anomalies that vary widely from event to event. These results also support the argument that a random pairing of flooding and observed storm surge events in Monte Carlo simulations will successfully capture the joint probability distribution of these events and their uncertainties due to sample size and other factors.
Hydrodynamic Modeling Results
Two approaches were used to model the physical dynamics of future flooding in order to consider coincident changes in flood volume, tidal cycle, tidal anomaly, and SLR described above. First, in order to assess some of the dynamics of future flooding given paired riverine flooding and sea level/storm surge, an ensemble of the three largest storms experienced in the Skagit River basin during the ECHAM5-WRF-2050s run were evaluated using the hydrodynamic model. These storms, paired with their coupled hourly tidal and daily TA signal (discussed above), were used as inputs to the hydrodynamic model for the “No Levee Failure” scenario (Figure 11D–F). The largest flood of the 2050s run occurred on 1/30/2069 and had a peak daily flow of 5441 m3/s. Although the regulated daily flow is comparable to the FEMA Q100 event, the inundated area created by the flood when combined with SLR is approximately 22% greater than the FEMA base case. The other two storms produced lower simulated inundation areas.
To eliminate any potential bias in regulated flows from the WRF scenarios and to facilitate a more direct comparison with the FEMA analysis, the second approach to modeling future floods scaled the FEMA Q100 storm to reflect projected increases in the Q100 flood for the HD-2040s and HD-2080s runs. These future runs used the idealized tidal cycle from the historical FEMA study, but with SLR projections for each time period added to the time series. Table 4 shows the amount of inundation for each of the scaled FEMA storms coupled with projected SLR estimates. Figure 11A–C shows the extent of inundation from the FEMA Q100 “No Levee Failure” scenario. The HD-2040 scenario resulted in a 57% increase in inundated area relative to the historical run (Figure 11B), and the HD-2080 scenario resulted in a 74% increase (Figure 11C).
A sensitivity analysis was performed to isolate the effects of SLR alone on the results for the 2040s. Combining the historical Q100 flood FEMA hydrograph with projected 2040s SLR (28.1 cm), produced increased inundated area of 35% (from 171 km2 to 230 km2) (Figure 11 G, Table 4). Combined with increases in flooding for the HD-2040s scenarios, inundation increases by 57% to 268 km2 (Figure 11 B, Table 4). Thus inundated area increases by an additional 22% due to an increase in peak flow during the event.
Summary of hydrodynamic modeling results for the “No Levee Failure” Scenario.
Figure 12 compares the composite flood maps for the historical, HD-2040s, and HD-2080s Q100 flood events involving FEMA's 7-levee failure scenario. The model simulations project average Q100 flood depths in the lower Skagit River 12.7 cm (5 inches) higher in the 2040s and 25.4 cm (10 inches) higher in 2080s compared to the historical base case.
As expected, larger peak flows and SLR produce more extensive inundation in the simulations. A summary of peak flow and inundated area is shown in Table 4. Increases in Q100 increase inundated area, but SLR alone is also shown to be an important factor. For example, peak flows comparable to the FEMA Q100 value, such as the largest WRF-2050s event (Figure 11 D) and the SLR sensitivity run (Figure 11 G), can create substantial increases in inundation area when combined with mid-21st-century SLR projections.
This study uses numerical modeling approaches to establish future flood risk in response to projected coincident changes in sea level, tidal anomaly, and regulated peak flows in the Skagit River projected for the 21st century. These changes are shown to dramatically increase extreme water levels and inundated area in the Skagit floodplain during high flow events (Table 4, Figures 11–12). It is important to note that changes happen not only with the extreme 100-year flood events, but also with smaller, more frequent events that regularly impact communities. For example, the HD-2040s 5-year flood is projected to carry 30% more volume than historically (Table 3), which could substantially increase impacts to communities in the floodplain.
A dynamic storm surge model based on season, ENSO, and simulated barometric pressure effects shows little systematic effect of climate change on storm surge, and instead the impacts of SLR dominate increases in near-coastal flooding forced from the marine side. Sensitivity analysis showed that changes in flooding due to SLR alone will be focused primarily in near-coastal areas, whereas increases in river flooding combined with SLR will affect the entire floodplain in the Skagit lowlands.
Analysis using dynamic downscaling demonstrated that peak flow events and tidal anomalies were essentially decoupled in the historical period, and remained so in the future scenarios. Furthermore, the probability distributions of storm surge were not meaningfully altered in comparison with historical observations. These findings supports the conclusion that base SLR will be the dominant effect on increased sea levels during flooding events, rather than storm surge.
As we move further into the twenty-first century, global climate change is expected to continue to alter the climate of the PNW. Assumptions of stationarity in flood plain management applications will become increasingly problematic under these conditions, and the results of this study highlight the need to develop forward-looking floodplain management strategies based on modeling rather than historical observations. The new approaches developed in this study are intended to help provide a detailed “road map” for future studies addressing the combined effects of future sea level rise, storm surge, and river flooding on extreme high water levels in coastal areas.
A single large-scale forcing from the ECHAM5 GCM was used in this study. More realizations of PNW climate using regional scale climate models will be needed to better characterize the changing risks in the Puget Sound lowlands and to better understand the effectiveness of RCMs in simulating both local-scale climate and global teleconnections (e.g. ENSO) affecting extremes. Increased sample size is also needed to improve GEV fits and better characterize high-flow extremes.
This study has been funded in part by the United States Environmental Protection Agency under assistance agreement 00J30401-0 to The Nature Conservancy. Thanks to Ted Perkins at FEMA for assistance obtaining and running the hydrodynamic model for the Skagit. Thanks also to FLO-2D Software, Inc. who provided free access to the FLO-2D hydrodynamic modeling software.