Data Resources

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).

Figure 2.

Flow chart for the chain of numerical models used in this project.

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.

Hydrologic Modeling

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).

Figure 3.

Probability of exceedance of regulated peak annual flow for the Skagit River at Mount Vernon. The figure compares three data sets: 1) USGS observations (1970–1999), 2) Simulated regulated flow based on HIST-VIC-SkagitSim simulations (1970–1999), and 3) Simulated regulated flow based on ECHAM5-WRF-VIC-SkagitSim simulations (1970–1999).

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.

Figure 4.

Observed vs. predicted (harmonic constants only) water levels at the Seattle Tidal Station during storm event from December 31, 1996 to January 3, 1997.

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:

Figure 5.

First and third EOFs derived from gridded ReAnal-WRF-1980s output of daily surface pressure anomalies.

This approach yielded a good fit in winter (R² > 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.

Figure 6.

A) Relationship between the Nino 3.4 Index and mean monthly tidal anomaly for Seattle, WA (1970–1999). B) Observed Niño3.4 index (black line), 1950–2012, and ECHAM5 simulated transient Niño3.4 index values (gray line), 1970–2069.

TABLE 1.

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).

Figure 7.

A) Example of disaggregated hourly flood event from December 5, 1989. Solid lines represent hourly hydrographs, dots represent mean daily flow. B) Example of the Q₁₀₀ FEMA hydrograph (gray line) scaled by a factor of 1.32 for the 2080s (black line).

To facilitate a more direct comparison between the floodplain impacts associated with the FEMA 100-yr flood (FEMA Q₁₀₀) and future projections of Q₁₀₀ using different downscaling methods, the FEMA Q₁₀₀hourly hydrograph was simply scaled to match the fractional change in the Q₁₀₀ 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.

Hydrodynamic Modeling

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 Q₁₀₀ peak daily flood magnitude was adjusted based on a variety of historical factors. For this reason, the magnitude of the FEMA Q₁₀₀ event (5772 m³/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 m³/s). For simplicity and in order to provide more direct comparisons with previous work, the scaled FEMA Q₁₀₀ was not pre-adjusted to account for this difference. Instead, the scaling factor was directly applied to the FEMA Q₁₀₀ 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.

TABLE 2.

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.

Figure 8.

Probability of exceedance plots for simulations of regulated peak annual streamflow for the Skagit River at Mount Vernon. A) shows statistically downscaled scenarios and B) shows dynamically downscaled scenarios.

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.

TABLE 3.

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 9.

A) Probability of exceedance of peak annual storm surge (i.e. the largest annual value predicted by TA regression model). B) Probability of exceedance of peak annual storm surge including sea level rise (i.e. SLR plus SS predicted by regression model).

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 m³/s. Although the regulated daily flow is comparable to the FEMA Q₁₀₀ 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.

Figure 10.

Relationship between regulated peak annual streamflow (y-axis) and coincident tidal anomaly (x-axis) for the Skagit River.

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 Q₁₀₀ storm to reflect projected increases in the Q₁₀₀ 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 Q₁₀₀ “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 Q₁₀₀ flood FEMA hydrograph with projected 2040s SLR (28.1 cm), produced increased inundated area of 35% (from 171 km² to 230 km²) (Figure 11 G, Table 4). Combined with increases in flooding for the HD-2040s scenarios, inundation increases by 57% to 268 km² (Figure 11 B, Table 4). Thus inundated area increases by an additional 22% due to an increase in peak flow during the event.

Figure 11.

Results from the hydrodynamic model for the no levee failure scenarios: (A) is the historical Q₁₀₀ event, (B–C) are the scaled FEMA Q₁₀₀ events for the HD-2040s and HD-2080s scenarios respectively, (D–F) are the three results for the 3 largest floods generated from the WRF-2050s scenario, and (G) shows the results from the SLR-Sensitivity run. The red outlines represent the bounds of the historical (FEMA) Q₁₀₀ event and are provided for spatial reference.

TABLE 4.

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 Q₁₀₀ flood events involving FEMA's 7-levee failure scenario. The model simulations project average Q₁₀₀ 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 Q₁₀₀ increase inundated area, but SLR alone is also shown to be an important factor. For example, peak flows comparable to the FEMA Q₁₀₀ 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.