Tracking an atmospheric river in a warmer climate: from water vapor to economic impacts

Atmospheric rivers (ARs) account for more than 75 % of heavy precipitation events and nearly all of the extreme flooding events along the Olympic Mountains and western Cascade Mountains of western Washington state. In a warmer climate, ARs in this region are projected to become more frequent and intense, primarily due to increases in atmospheric water vapor. However, it is unclear how the changes in water vapor transport will affect regional flooding and associated economic impacts. In this work we present an integrated modeling system to quantify the atmospheric–hydrologic–hydraulic and economic impacts of the December 2007 AR event that impacted the Chehalis River basin in western Washington. We use the modeling system to project impacts under a hypothetical scenario in which the same December 2007 event occurs in a warmer climate. This method allows us to incorporate different types of uncertainty, including (a) alternative future radiative forcings, (b) different responses of the climate system to future radiative forcings and (c) different responses of the surface hydrologic system. In the warming scenario, AR integrated vapor transport increases; however, these changes do not translate into generalized increases in precipitation throughout the basin. The changes in precipitation translate into spatially heterogeneous changes in sub-basin runoff and increased streamflow along the entire Chehalis main stem. Economic losses due to stock damages increase moderately, but losses in terms of business interruption are significant. Our integrated modeling tool provides communities in the Chehalis region with a range of possible future physical and economic impacts associated with AR flooding.

1 Introduction

On 3 December 2007, an atmospheric river (AR) event made landfall on the West Coast of the US. The resulting extreme precipitation event severely impacted the Chehalis River basin in western Washington and resulted in 6 h rainfall amounts close to the 100-year storm volume (NOAA, 2008). Two previous storms (on 1 and 2 December) brought heavy snow to the Oregon Coastal Range, the Olympic Mountains and the Cascades, while the third and strongest event brought mostly liquid precipitation. The hurricane force winds on 3 December produced wind damage with tree blowdowns, power outages, huge ocean swells and a record coastal storm surge. Eleven people lost their lives. Millions of people lost power throughout Washington and Oregon as a result of the storm. Portions of interstate 5, the major north–south freight corridor on the West Coast connecting the Puget Sound region of Washington with Oregon and California, were closed for 4 days, resulting in an estimated USD 47 million in economic losses (Ivanov et al., 2008). Major disaster declarations were issued in several counties in the states of Washington and Oregon, but most of the damages were concentrated in three counties in Washington: Grays Harbor, Lewis and Thurston. Lewis County, within which the most affected part of the Chehalis River basin lies, experienced the largest impact with USD 166 million in damages and 46 % of its agricultural land flooded (Lewis County, WA, 2009).

While this event was particularly extreme, more than 50 % of the total cool-season precipitation and more than 75 % of heavy precipitation (top decile) on the west coast of Oregon and Washington is related to AR events (Rutz et al., 2014). Water vapor transport during the winter season is often roughly orthogonal to the mountain ranges, which favors orographic precipitation enhancement (Hu et al., 2017; Neiman et al., 2008). Furthermore, ARs with anomalous warm, strong, low-level water vapor fluxes are responsible for nearly all of the extreme flooding along the Olympic Mountains and the western Cascade Mountains of Washington (Hu et al., 2017; Neiman et al., 2011; Warner et al., 2012).

Given the critical role of ARs for precipitation and flooding in the region, it is important to understand how these could change in a warmer climate. As tropospheric temperature increases, integrated water vapor transport (IVT) is projected to increase by 30–40 % by the end of the 21st century along the North Pacific storm tracks, including the West Coast of the US (Lavers et al., 2015; Salathé et al., 2015). In climate model projections, years with many AR storms are projected to become more frequent and water vapor content is projected to increase during intense AR events (Dettinger, 2011). The changes in IVT are driven mostly by thermodynamics through increased water vapor content of a warmer atmosphere, while changes in dynamics seem to have only a secondary effect along the northern West Coast of the US (Lavers et al., 2015; Payne and Magnusdottir, 2015; Salathé et al., 2015). Based on the analysis of IVT changes, it is tempting to conclude that the projected increase in intensity and frequency of AR events will lead to increased flooding in the region. However, to quantify the risk of inundation and its economic impact, it is important to understand the myriad of processes that happen between the impact of an AR in a watershed and the resulting flooding.

In this work, we present an integrated modeling system that quantifies the atmospheric–hydrologic–hydraulic and economic impacts of the December 2007 AR event. In addition, we use the modeling system to project physical and economic impacts under a scenario in which the same December 2007 event occurs in an atmosphere with increased greenhouse gas forcing. As opposed to a traditional approach that uses an ensemble of downscaled and bias-corrected climate model simulations, we use the regional model simulations of the December 2007 event in hypothetical future climate settings. We then use these high-resolution simulations in a warmer climate as forcing for the hydrology–hydraulic and economic loss models. Our work follows a similar procedure as the US Geological Survey (USGS) Multihazards Project, which used a synthetic but plausible California AR scenario to estimate the human, infrastructure, economic and environmental impacts for emergency preparedness and flood planning exercises (Porter et al., 2010). In our work, we focus on the Chehalis River basin in western Washington to provide an end-to-end model of severe weather, physical impacts and economic consequences of ARs in a warmer climate.

The integrated modeling system allows us to incorporate different types of uncertainty, including (a) alternative future radiative forcings associated with different Coupled Model Intercomparison Project 5 (CMIP5) Representative Concentration Pathways (RCPs) – RCP4.5 and 8.5 scenarios, (b) different possible responses of the climate system to future radiative forcings as represented by 14 CMIP5 GCMs and (c) different possible responses of the hydrologic system as represented by two different hydrologic models. We do not account for possible changes or structural failures in the main channel hydraulics and we do not account for possible changes in private or public building infrastructure or trade flows. At each step in the modeling chain we provide an envelope of possible future responses of the system and present them as changes with respect to the historical control simulation. The modeling system is intended to provide decision makers with information about the range of physically plausible changes in flood-causing AR storms and floods, as well as a tool to quantify the related economic impacts.

2 Data and methods

The Chehalis River basin, with a drainage area of approximately 5400 square kilometers, is located in southwest Washington state (Fig. 1b). It heads in the Willapa Hills, flows east, then north and west into Grays Harbor. Most of the basin lies below 1000 m of elevation. Fall and winter precipitation mostly occurs as rain, with exceptions in small areas of the extreme northern and eastern portions of the basin. Floods in the basin generally occur in late fall and early winter and are associated with atmospheric rivers. The most significant floods in the observational period are January 1972, January 1990, November 1990, February 1996, December 2007 and January 2009 (USGS). We focus on the largest event recorded in the basin, the December 2007 event.

On 3 December 2007, an AR filamentary plume transporting more than 2000 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{s}}^{-\mathrm{1}}$ of water vapor at its core extended from the tropical Pacific, west of Hawaii, to the coast of Oregon and Washington (Fig. 1a). Selecting the cross section of the AR with the most intense transport and integrating IVT for all values exceeding 1500 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{s}}^{-\mathrm{1}}$, we can calculate the equivalent liquid water discharge. This AR carried approximately 847 000 m³ s⁻¹ of liquid water across its inner core, or the equivalent of about 50 times the average discharge at the mouth of the Mississippi River. Temperatures rose 17 ^∘C in less than 2 h ahead of the cold front (NOAA, 2008). Along this warm southwesterly tropical air mass, more than 70 % of the water vapor and precipitation that reached the coast was of direct tropical origin (Eiras-Barca et al., 2017). The catastrophic flooding along the Chehalis River basin was primarily due to unusually high and sustained hourly rainfall rates concentrated in less than 24 h, mainly on 3 December. The conditions were exacerbated by warm air advection into the region by the AR, which produced rain on snow conditions and partially melted the existing shallow, low-elevation snow. Ten USGS stream gauges experienced record flooding, including four on the Chehalis River or its tributaries (Grand Mound, Porter, Doty and South Fork Chehalis; see Fig. 1b for station locations). The peak discharge measured at Doty was a 500-year event – the only 500-year stream peak event ever recorded in western Washington.

Figure 1 (a) Integrated vapor transport (IVT) $\mathrm{kg}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{s}}^{-\mathrm{1}}$ on 3 December 2007 from ERA-Interim Reanalysis; dashed lines are the WRF outer and inner domain. (b) Chehalis River basin with topographical features and the largest urban areas (Centralia and Chehalis). The Chehalis main channel as represented in HEC-RAS is shown, along with the USGS gauging stations (red triangles) and precipitation stations (yellow circles) used in this study. Numbers correspond to the station information in Table 1.

Table 1 Streamflow and precipitation observations. Map ID corresponds to the locations on the map of Fig. 1b.

2.1 Data: observations

We used the 1∕16^∘ latitude–longitude daily gridded precipitation product derived from NOAA Cooperative Observer (COOP) stations by Livneh et al. (2013). In addition, we used hourly data from seven NOAA (four COOP and three HADS) stations in and around the Chehalis basin (Fig. 1b and Table 1). We used USGS streamflow observations from 15 gauges located throughout the basin (Fig. 1b and Table 1). During the flood event, the upstream-most gauge (Doty) measured streamflow up to approximately 60 000 cfs, but then malfunctioned during the time of peak flood (WSE, 2012); consequently, the peak discharge had to be estimated by the USGS. In addition, we used the European Centre for Medium-Range Weather Forecasts (ECMWF) Interim Reanalysis (ERA-Interim) (Dee et al., 2011) at 0.75^∘ resolution as lateral boundary conditions for Weather Research and Forecast (WRF) atmospheric model simulations. In terms of direct economic losses, we rely on infrastructure data and dasymetric dataset for buildings, which is embedded in the standard release of HAZUS-MH 3.0. To calculate their ripple effects throughout the local supply chain (also called indirect losses) we rely on the 2008 input–output tables from IMPLAN (2015). The sector-specific inoperability levels and sector-specific recovery rates are calculated using the inventories of finished goods. Input–output data contain information about trade flows across 16 different sectors that represent the economic structure of each of the counties within the state of Washington. They were obtained from IMPLAN (2015).

Figure 2 Diagram of the integrated modeling, including the models used and the input data for each model during the historical simulations (a) and the climate change simulations (b). Hydro-control represents both HEC-HMS and DHSVM-control simulations, while Hydro-PGW represents both HEC-HMS and DHSVM-PGW simulations.

2.2 Methods: models

Our atmospheric simulations of the December 2007 event used the Advanced Research version (ARW) of the WRF model (Skamarock et al., 2005), version 3.4.1, with two nested domains, one of 15 km and the inner domain of 3 km (Fig. 1a). The time period for our simulation is 30 November to 8 December 2007. The physics options used are the YSU planetary boundary layer scheme (Hong and Pan, 2009), subgrid-scale convection in the 15 km grid based on the Kain–Fritsch parameterization (Kain, 2010), WSM six-class microphysics (Hong and Lim, 2006) and the Noah-LSM V1.0 (Chen and Dudhia, 2001) land surface model. We tested other microphysics schemes, but we found that the WSM six-class yielded precipitation that was closest to observations.

Our hydrologic simulations used two different models to estimate the response of the Chehalis watershed to precipitation: the US Army Corps of Engineers (USACE) Hydrologic Engineering Center (HEC) Hydrologic Modeling System (HEC-HMS) and the University of Washington's Distributed Hydrology Soil Vegetation Model (DHSVM) hydrologic model (Wigmosta et al., 1994). Our goal in using the two models is to account for uncertainty in the physical representation of hydrologic processes. In HEC-HMS, we partitioned the watershed into 64 sub-basins with homogenous soil and land cover properties based on data from SSURGO (USDA-NRCS) and NLCD 2011 (Homer et al., 2015). HEC-HMS provides the streamflow response of each of the sub-basins that drain to the Chehalis main channel. We calculated base flow in three different ways: if there was a stream gauge, we used the USGS stream statistics; if the stream gauge was located downstream of a tributary, we calculated the initial base flow for the channel receiving from each sub-basin based on the fraction of the gauged area contributed by each sub-basin in the tributary; if there were no stream gauges available, we estimated the initial base flow through analogy with similar-sized sub-basins nearby. We used the Green and Ampt option in HEC-HMS to simulate infiltration in each sub-basin. Given the limited observations, we estimated the Green and Ampt parameters (saturated hydraulic conductivity, effective porosity and wetting front suction head) based on the values reported in the literature for each hydraulic soil group. For each sub-basin, we used the area-weighted properties. For the purposes of calculating soil infiltration rates, we estimated percent impervious area using the land use and land cover maps obtained from SSURGO. The runoff transform uses the Soil Conservation Service (SCS) lag time. The HEC-HMS simulated streamflow was compared to the observed streamflow at the USGS gauges listed in Table 1. The only parameter that was calibrated was the soil infiltration parameter, which was adjusted within the range of each soil type. In addition, the final model setup with 64 sub-basins of homogeneous soil and land cover types was found to be the optimum representation of the basin and it resulted in streamflow closest to observations. If the basin is represented with fewer sub-basins, the HEC-HMS simulated streamflow does not capture the timing or magnitude of the peak in the observed hydrographs.

DHSVM is an explicit, physically based, spatially distributed hydrological model developed primarily for use in regions with complex terrain. Unlike HEC-HMS, DHSVM uses a rectangular grid formulation, here with a spatial resolution of 150 m. DHSVM represents runoff primarily through the saturation excess mechanism using a representation of a shallow water table whose depth is modeled similarly to TopModel (Beven and Kirkby, 1979), with the exception that the spatial variation in depth to the water table is represented explicitly rather than statistically. At each grid cell, unsaturated moisture flow through the root zone is computed using a prescribed hydraulic conductivity that decays exponentially at the water table depth to the saturated hydraulic conductivity. Redistribution of moisture between pixels occurs (only) in the saturated zone where the hydraulic gradient is taken to be equal to the (computed) slope of the water table, following Wigmosta and Lettenmaier (1999). The model uses a linear storage scheme to route both overland and subsurface flow (which occurs at the intersection of the water table and the stream network) through a channel network identified using digital topographic data. We calibrated DHSVM using observed daily streamflow at the USGS stream gauges. To calibrate DHSVM for the 2007 storm we initially implemented a simple sensitivity analysis. DHSVM uses 18 different soil types, which the model links internally to soil hydraulic properties (e.g., saturated hydraulic conductivity, porosity, etc). We then determined sensitivity to the three dominant initial soil types (as suggested by Cuo et al., 2011) and other selected model parameters. We found that the soil maximum infiltration rate and Manning's roughness coefficient (for channel flow) were the most sensitive parameters. We then developed a Monte Carlo simulation approach that randomly picked these parameters (between the prescribed upper bounds and lower bounds defined by Cuo et al., 2011). We compared simulated flows with USGS gauge station observed streamflow (using RMSE) and identified the optimal parameter combinations within each sub-basin.

We used the output from the two hydrologic models as boundary conditions for the USACE River Analysis System (HEC-RAS) one-dimensional unsteady flow model to perform hydraulic simulations of water levels in the Chehalis River main stem and its largest tributaries. The calibrated HEC-RAS model was provided to our team by USACE. USACE and its contractor, Watershed Science and Engineering (WSE), updated previously existing hydraulic models of the Chehalis River based on data from a bathymetric survey performed by WSE and available lidar data. They then calibrated the updated model based on hydrologic observations in the watershed. The hydraulic model extends from the mouth of the Chehalis River to upstream of Pe Ell (173 km). The model includes portions of the following tributaries: Wynoochee River, Satsop River, Black River, Skookumchuck River, Newaukum River and South Fork Chehalis (Fig. 1b). HEC-RAS output includes river stage and streamflow calculations at each channel cross section, flood inundation extent and flood inundation depth. WSE calibrated the model to the February 1996 and January 2009 storm events and used the December 2007 storm event for validation. WSE adjusted channel and overbank values of Manning's n bottom roughness coefficient, flow roughness factors and the placement of ineffective flow areas in their calibration process. The HEC-RAS model provided by USACE used observed streamflow hydrographs as lateral boundary conditions; for this reason, we developed our own hydrologic models, as described above, to provide flexibility in our simulations of alternative storm scenarios.

We calculated the direct economic losses using HAZUS (HAZard USa), a software developed by the Federal Emergency Management Agency (FEMA, 2003), to calculate economic losses associated with different natural disasters, including floods (see, among others, Banks et al., 2014; Ding et al., 2008; Gutenson et al., 2015). We used HAZUS-MH to calculate how the HEC-RAS simulated flooding led to direct economic losses to agriculture (crops), buildings and public infrastructure such as telecommunication lines and roads. The dasymetric data embedded in HAZUS include information about the location and characteristics of the buildings and infrastructures (e.g., number of floors in a building, number of lanes in a road). These data allocate the use of land and buildings by economic sectors so that one can estimate how the direct economic losses result in direct production capacity constraints and losses by sector. Our HAZUS implementation contains several assumptions: as usual in the literature, production capacity constraints are based on the assumption of homogeneous productivity per square foot for each industry in a specific county and on the assumption that industries operated at full capacity before the disaster. As a result, we set the production capacity constraints based on the pre-disaster total output by industry. While HAZUS is able to calculate damages to crops and some crop areas were flooded during the event, crop losses are null because the event took place several months before the planting season. Buildings located on farmland were damaged, however, and their repair or reconstruction costs follow the same methodology as similar costs, as described further below.

Because each company or institution relies on a set of suppliers and purchasers to support its activities, they too will experience production losses as a result of the flood even though they have not been flooded themselves. These indirect economic losses are estimated from the 2008 input–output tables extracted from IMPLAN at a 16-sector aggregation level (Avelino and Dall'erba, 2016). In addition to production losses, the combination of HAZUS and input–output techniques allows us to quantify how local final demand decreases as a result of the employees suffering from labor income losses due to temporary closure of their workplace. We assume that the expenditure structure remains fixed in the post-disaster period and that demand decreases proportionally to the decrease in income. Reconstruction costs, on the other hand, correspond to a positive stimulus encompassing the total repair costs of buildings, infrastructure and vehicles that were destroyed or damaged during the flood. Since IO models are based on producer prices and HAZUS provides repair costs in purchase prices, we assume that manufacturing orders include margins split 20∕80 % between transportation and trade. Due to the small size of the economy of the affected counties, the model assumes that reconstruction efforts are supplied by companies located outside of the flooded area. The duration of the recovery phase is given by HAZUS (Tables 14.1, 14.5 and 14.12 of FEMA, 2015) and assumed to be linear in time. The total economic impact in the three affected counties and the rest of Washington is then estimated using the Inventory-Dynamic Inoperability Input–Output Model (Inv-DIIM) proposed by Barker and Santos (2010). In relation to other available input–output models, the Inv-DIIM offers a dynamic view of inoperability and recovery processes in addition to accounting for available inventories that can alleviate disruptions in the region (Avelino and Dall'erba, 2016). The inventory data for the DIIM are based on the December 2007 inventory-to-sales ratio for manufacturing reported by the Federal Reserve Bank of St. Louis in 2016. This ratio has been suggested by Barker and Santos (2010) and is equivalent to 1.23 for the period under study. We apply it homogeneously to all counties. Since the activities of wholesale and retail are recorded as margins, these sectors do not hold finished goods inventories.

Xiao-Ge et al. (2013)Arora et al. (2011)Gent et al. (2011)Voldoire et al. (2013)Rotstayn et al. (2010)Volodin et al. (2010)Dufresne et al. (2013)Watanabe et al. (2010)Watanabe et al. (2010)Zanchettin et al. (2013)Zhang and Yan (2012)Donner et al. (2011)Donner et al. (2011)Jones et al. (2011)

Table 2 CMIP5 GCM models used in this study, including the respective RCP scenario used.

2.3 Methods: climate change simulations

To understand how the December 2007 event would change if it occurred in a warmer climate, we used a “pseudo-global warming” (PGW) approach (Kawase et al., 2009; Lackmann, 2013, 2015; Lynn and Druyan, 2009; Rasmussen et al., 2011; Sato et al., 2007; Schär et al., 1996). The PGW can provide complementary information to the traditional downscaling approach as it gives more physical insight into detailed spatial processes and potentially a better way of communicating with regional stakeholders, as argued by Hazeleger et al. (2015). In this approach, the lateral and initial boundary conditions used in the WRF-control simulation are modified by adding a perturbation “delta” to reflect future changes in temperature as simulated by global climate model (GCM) projections for the future. We only modified vertical and surface temperature and SSTs, while increasing the specific humidity to maintain constant relative humidity. In this way, we ensured that the storm dynamics remain unchanged (Schär et al., 1996). It is important to emphasize that this method does not account for possible changes in large-scale dynamics, such as changes in the storm track. However, it has been shown that the changes in future AR events in this region are dominated by thermodynamic (changes in humidity) as opposed to dynamic processes (changes in wind) (Lavers et al., 2015; Payne and Magnusdottir, 2015; Salathé et al., 2015). For this reason, the PGW method provides useful information about possible future AR changes in the Chehalis basin.

The 14 different CMIP5 global climate models used to calculate the changes in temperature over the region (WRF model outer domain) are listed in Table 2. Based on one simulation from each model for two different Representative Concentration Pathway scenarios (RCP4.5 and RCP8.5), we obtained an envelope of possible changes in temperature between the future (2071–2098) and the historical (1980–2004) mean December–January–February (DJF) temperatures (Fig. 3). We denote “lower” as the smallest change in temperature and “upper” as the largest. Surface temperature changes range between approximately 1 and 4 K, increase to between 2 and 6 K around 350 mb and then decrease sharply to approximately −1 to 2 K at 50 mb. These patterns are similar to the global-averaged changes in temperature, which have maximum warming in the upper troposphere and cooling in the stratosphere (IPCC, 2013).

Figure 3 Upper (red) and lower (green) bounds of the area-averaged temperature changes, as represented by the 14 CMIP5 models listed in Table 2, using the RCP4.5 and RCP8.5 simulations.

We interpolated the domain-averaged changes in temperature from the upper and lower scenarios to the same 26 vertical levels of ERA-Interim. Then, we added these deltas to the ERA-Interim forcing to perform two simulations, one for the upper scenario and one with the lower scenario. In this way, we are only evaluating the change in precipitation due to horizontally homogeneous changes in temperature – all other variables remain exactly the same as in the control simulation. This ensures that the AR's path and orientation do not change due to changes in atmospheric dynamics (see mathematical derivation in Schär et al., 1996). This is important because AR precipitation is strongly influenced by the angle of impingement on regional topography (Hu et al., 2017).

“Delta method” for model simulations

Each model is sensitive to its input data. In particular, the socioeconomic evaluation requires precise information about the spatial location and depth of inundation. For this reason, in each part of the model chain, we decided not to use the raw model data but rather the changes in total water flux as simulated by the different models (see Fig. 2). Our strategy for each model simulation was as follows.

1. We performed control simulations of each model forced with observed or reanalysis data (WRF is forced by ERA-Interim, HEC-HMS and DHSVM are driven by observed precipitation, HEC-RAS is forced by observed streamflow). Due to a lack of observed maximum flood extent, we forced HAZUS with the inundation depth and extent as modeled by the control HEC-RAS simulation.

2. We calibrated each model so as to best simulate the relevant observations.

3. We ran WRF with the PGW conditions, both the upper and lower scenarios, and obtained changes in precipitation (WRF-PGW).

4. Based on the ratio of WRF-PGW and WRF-control precipitation, we obtain a percent change in precipitation over the entire 1–4 December period. We modified the observed precipitation by this percent change and then ran the hydrologic models with modified precipitation (HEC-HMS-PGW and DHSVM-PGW).

5. Based on the ratio of HEC-HMS-PGW and HEC-HMS-control (and DHSVM-PGW to DHSVM-control) streamflow, for each type of inflow into the main Chehalis channel we obtain a percent change in total streamflow volume for the 1–7 December period. We modified the observed streamflow by this percentage change and then ran HEC-RAS with modified streamflow (HEC-RAS-PGW).

6. Based on the new HEC-RAS-PGW inundation extent and depth, we ran HAZUS and our input–output model to obtain new economic loss estimates.

Figure 4 (a) Observed daily precipitation (mm day⁻¹) averaged for 1–4 December 2007 from Livneh et al. (2013), (b) WRF-control simulated precipitation for the same period and (c) bias in simulated precipitation for each of the HEC-HMS sub-basins within the Chehalis basin.

Figure 5 USGS observed (solid black), HEC-HMS simulated control (dashed blue) and DHSVM simulated (dotted blue) discharge for four representative sub-basins within the Chehalis.

3 Results: historical simulations

The WRF-control simulation captures the observed extreme precipitation over the Oregon Coastal Range and Olympic Mountains with precipitation on the order of 80 mm day⁻¹ over some areas (Fig. 4a and b). However, the simulation overestimates precipitation over the Cascades and underestimates precipitation over most of the Chehalis basin by about 30–40 % (Fig. 4c). The simulation overestimates precipitation over the Willapa Hills in the southern part of the basin.

HEC-HMS captures the timing of peak stage and flow; however, it has problems with underestimation of peak flow and more generally underestimates discharge throughout most of the basin (Fig. 5, dashed lines). DHSVM, on the other hand, adequately captures peak flow in the upper basin (Doty and Newaukum), but overestimates peak discharge in the lower basin and underestimates recession flows (Fig. 5, dotted lines). As explained above, the hydrologic models are different and they both have strengths and weaknesses in simulating different parts of the hydrograph at different locations. It is important to note that we used a combination of Livneh precipitation data (daily timescale) with hourly data from five NOAA stations (shown in Fig. 1) to partition the Livneh daily totals. Hence, while the total daily volumes match the Livneh product, the hourly variability comes from the station data. There is considerable uncertainty in the Livneh precipitation product daily totals for this storm and even more uncertainty as to the hourly precipitation throughout the basin. Errors in the hydrologic response are largely due to error in the precipitation estimates. Since the 2007 flood, an NWS precipitation radar has been installed (at Langley Hill) and the number of HADS stations has increased, helping to better resolve the space–time distribution of precipitation over the basin. These assets were not, however, available during the 2007 storm.

The calibrated HEC-RAS hydraulic model, driven by observed streamflow from the USGS stations (see station locations in Fig. 1), performs very well (Fig. 6). The differences between the observed and simulated stage along the Chehalis main stem range from −0.54 to 0.65 m, while the difference in peak flow magnitude ranges from about −1.4 % at Doty (upstream) to −16.9 % at Porter (downstream). The resulting inundation depth and extent are shown in Fig. 6. Large areas around the cities of Chehalis and Centralia (see Fig. 1b for location) were inundated.

We used the inundated areas and depths from HEC-RAS to calculate the local damages to arable land, buildings and content, infrastructure and vehicles using HAZUS. Then the net loss in local production is calculated using the Inv-DIIM. The total physical damages for Lewis, Thurston and Grays Harbor combined were estimated at USD 678 million with business disruption losses of USD 51 million (Table 3, “Base” rows), most of which was in Lewis County (Avelino and Dall'erba, 2016). While reported loss estimates are difficult to obtain, the Department of Commerce estimated that losses for the states of Washington and Oregon combined for this flooding event were approximately USD 1 billion, so our estimates for the three counties seem reasonable. In addition, the official building and inventory damages in Lewis county were estimated at USD 166 million (Lewis County, WA, 2009), which is close to our estimate of USD 151 million for the same categories. It is important to clarify that we do not have a counterfactual that can be used to calibrate the economic model in the same way that we calibrate the physical models.

Overall, we find that individually the models of the integrated system realistically capture the dominant physical and economic processes. However, it is clear that there are problems with some variables, particularly precipitation and the associated hydrologic response. For this reason, we decided not to use the raw model output (from WRF, HEC-HMS, DHSVM or HEC-RAS) to drive the subsequent model in the historical simulations (Fig. 2). Instead, we use the individual historical model simulations forced with observations. To simulate the climate change response, the observations are then multiplied by a factor that accounts for the changes projected by the models in a warmer climate (as described in Sect. 2.3). The underlying idea is that the models cannot provide precise spatiotemporal values of the different variables; however, because their representation of the dominant processes is realistic, we trust they are able to capture the changes between the past and the future. This is the rationale behind the “delta method”.

Figure 6 (a, c, d) USGS observed (solid) and simulated (dashed) stage for three cross sections of the Chehalis River main stem as represented by HEC-RAS. (b). Flood extent and depth map as simulated by HEC-RAS.

Figure 7 (a) WRF-control simulated IVT ($\mathrm{kg}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{s}}^{-\mathrm{1}}$) for 3 December 2007, (b) WRF-PGW simulated IVT ($\mathrm{kg}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{s}}^{-\mathrm{1}}$) for 3 December 2007 for the upper scenario and (c) absolute change in IVT between the WRF-PGW upper scenario and WRF-control.

Figure 8 Changes between WRF-PGW for the upper scenario and WRF-control (inner WRF domain) averaged for the 1–4 December period for (a) water vapor mixing ratio percent change at 800 mb, (b) relative humidity percent change at 800 mb, (c) absolute change in cloud water mixing ratio at 800 mb and (d) percent change in precipitation; (e) percent precipitation change area averaged over all Chehalis sub-basins of the HEC-HMS model.

Figure 9 (a) Percent change in precipitation for the lower (green) and upper (red) scenarios, as simulated by WRF for all Chehalis sub-basins used in the HEC-HMS simulations. (b) Percent change in streamflow for the lower (green) and upper (red) scenarios, as simulated by HEC-HMS for all Chehalis sub-basins

Figure 10 Streamflow hydrographs for the HEC-HMS-PGW upper scenario (dashed) and HEC-HMS-control (solid) for select sub-basins in the Chehalis.

4 Results: climate change simulations

In the WRF-PWG simulation, we added the changes in temperature shown in Fig. 3 (both upper and lower scenarios) to each level of the ERA-Interim boundary conditions used in the control simulation, while maintaining constant relative humidity. This necessarily implies an increase in the specific humidity, as higher temperatures increase the saturation specific humidity. These changes induce variations in the IVT of the projected AR event, which increases by 12.6 % in the lower scenario to 38.5 % in the upper scenario for the WRF outer domain (Fig. 7 shows the spatial changes for the PGW-upper scenario). The increase approximately follows the Clausius–Clapeyron scaling of about 7 % per degree of warming. The increase in IVT can be as large as 500 $\mathrm{kg}{\mathrm{m}}^{-\mathrm{1}}{\mathrm{s}}^{-\mathrm{1}}$ throughout the AR corridor. IVT also increases within the inner WRF domain by 12.4 to 42.3 % for the two scenarios. The water vapor mixing ratio increases everywhere, but not homogeneously in space (Fig. 8a), with a clear structure of changes above 40 % at the 800 mb level. However, due to the differences in temperature, the relative humidity can increase or decrease in the PGW-upper simulation, and this leads to both positive and negative changes in the cloud water mixing ratio (Fig. 8b and c). In Fig. 8 we show these results at the 800 mb level, but these heterogenous changes in relative humidity and cloud water can be seen throughout the lower troposphere. As a consequence, precipitation shows both areas of significant increase and decrease throughout the WRF inner domain. The inner domain area-averaged precipitation change is 8.2 % for the lower scenario and 17.8 % for the upper scenario – significantly below the Clausius–Clapeyron scaling. On the basin scale, precipitation increases significantly (exceeding 30 %) in the northern part of the watershed and deceases significantly (below 30 %) in the southeastern Chehalis basin (Fig. 8e). We calculated the fractional changes in precipitation for each sub-watershed as the total precipitation that accumulated between 1 and 4 December of the WRF-PGW simulation divided by the WRF-control accumulated precipitation for the same period (Fig. 9a). The upper basin (lowest sub-basin numbers) clearly shows precipitation increases, the eastern part of the basin shows decreased precipitation and the lower basin shows increased precipitation.

Figure 11 Streamflow hydrographs for observed (black) and simulated using HEC-HMS (dashed) and DHSVM (dotted) for the lower scenario (green) and higher scenario (red).

Figure 12 Change in streamflow and flood depth along main channel.

Table 3 Projected economic losses for the historical simulations (Base) and the upper and lower scenarios for the two hydrologic models. Values are in millions of US dollars (at 2008 rates). Values in parentheses represent losses.

We multiplied the observed precipitation by the fractional change in precipitation (shown in Fig. 9a for each HEC-HMS sub-basin) and used the result to force the HEC-HMS and DHSVM PGW simulations. There are two different scenarios that result in four different hydrologic simulations (HEC-HMS-lower, HEC-HMS-upper, DHSVM-lower and DHSVM-upper). The results show that some regions generate significantly more runoff due to increased precipitation, while the southeastern part of the basin generates less runoff (Figs. 9b and 10). Notably, the Doty station in the headwaters of the basin shows an increase in peak runoff that ranges from 13 % in DHSVM-lower to 44 % in HEC-HMS-upper. The use of the two hydrologic models provides an envelope of uncertainty in the numerical representation of the hydrologic response (Fig. 10). We find that the sharp increase in streamflow in the headwaters dominates the response in the main channel, as simulated by HEC-RAS (Fig. 11). There is an increase in both stage and flow throughout most of the channel, with increases that range from about 12–42 % in the headwaters (depending on the scenario), to −6 to 5 % in the eastern part of the basin and then about 10–30 % at the outlet into Grays Harbor (Figs. 11 and 12a). Only the DHSVM-lower scenario shows small decreases in the eastern part of the basin.

The associated socioeconomic losses, as simulated by HAZUS and Inv-DIIM, show an increase in physical damages of 2–33 % in Grays Harbor County, 9–171 % in Lewis County and −1–10 % in Thurston County. The results are sensitive to the scenario and the hydrologic model used (Table 3). Our results indicate a larger loss in Lewis County because it is where Centralia and Chehalis, two of the most populated cities of our watershed, are located and they hold the largest stock of private and public buildings and infrastructures. Interestingly, in terms of business interruption losses (Table 3 lower), the increases are substantially higher and can be very different from the changes in physical damages (27–250, 14–314 and 46–619 %, respectively). The economy outside of these three counties (“Rest of WA”, Table 3) is positively impacted as reconstruction and recovery efforts stimulate production in the rest of Washington. As a result, the net impact on local statewide production and internal trade is positive.

Figure 13 Methodology to calculate the historical and future expected annual losses using only HAZUS and streamflow observations. (a) Flow duration curve for the Porter gauge and the fitted lognormal distribution. (b) Fitted streamflow for different return periods for the historical period (blue) and the future (red). The changes in streamflow in the future are calculated by assuming a 15 % increase in streamflow in the future. We then calculate the changes in return period. (c) Economic loss probability curve for the current and future period.

4.1 Interpretation

Despite the fact that some sub-basins experience lower streamflow in the climate change simulation (see Skookumchuck and Newaukum in Fig. 10), streamflow throughout the main stem of the Chehalis increased. This implies that the dramatic increases in flooding of the headwaters (see Doty in Fig. 10) dominated the system response and caused flooding in populated downstream areas along the main stem of the river, including Centralia and Chehalis (the largest population centers in the basin). Our results highlight the fact that the economic impacts are very sensitive to the geographical location of inundated area and depth. The parts of the basin with large population centers are most vulnerable to direct economic losses and account for most of the stock damages (Table 3). But this is not the only factor. Indeed, Thurston County has strong trade linkages to other regions (such as the Seattle metropolitan area) and for this reason, despite modest changes in direct impacts, the net impact on trade increased significantly in the climate change simulation (480–619 %; Table 3). This indicates that, depending on the hydrological impacts, the simulated economic scenarios can lead to flooding patterns that impact key interconnected sectors of the economy, significantly increasing negative spillover effects.

Interestingly, despite general increases in streamflow in the climate change simulation, the changes in inundation extent are minimal (Fig. 12b). The reason for this is that the December 2007 event was so large that the flooding extended throughout much of the flood plain to the bounding and steeper hills. As a result, the changes in economic impacts might not be very large for an event of such low probability of exceedance. Smaller events may well be (proportionately) more affected under climate change in this river basin (clearly, the extent of the flood plain and the characteristics of the bounding topography are basin specific). We were able to get some insight into the nature of the basin's response to changes in more modest floods using a simplified method (in contrast to the full chain of model calculations that underlie our estimates for the 2007 flood) by using the default data for flood extent and depth for different return periods from HAZUS (without performing the atmospheric, hydrologic or hydraulic analysis) and applying the changes to gauge observations.

The Porter stream gauge (gauge 10 in Fig. 1) provides representative data for the entire watershed and allows us to identify the streamflow for different return periods (Fig. 13). Assuming that climate change will result in 15 % more streamflow for all return periods (an assumption based on our PGW results and results from Hamlet and Lettenmaier, 2007; see their Fig. 10), we used HAZUS and a method similar to Velasco et al. (2015) to evaluate the losses for historical and future events. We then calculated expected total losses for the historical period as the integral under the blue curve in Fig. 13c (USD 6.2 million) and expected total losses for the changed climate condition as the integral under the red curve (USD 8.6 million) for a total increase in expected losses of 39 %. In the future, we plan to repeat this analysis using the full integrated model chain to obtain more realistic values for the changes in streamflow, which would replace the assumed 15 % increase in streamflow independent of return period.

5 Conclusions

ARs are responsible for most of the extreme winter flooding events in the western US. As the climate warms, the thermodynamic response of these atmospheric structures will likely lead to significantly more water vapor content and fluxes. Others have hypothesized that a warmer climate will lead to more intense AR-related flooding events and societal impacts. However, the way that the water vapor carried by an AR is transformed into precipitation, runoff and streamflow along a channel is highly nonlinear and depends on a myriad of fine-scale processes both in the atmosphere and on the land surface. Furthermore, the economic impacts depend on the human footprint, the economic structures in the affected areas and their trade linkages with other regions. Because of the risk associated with these events, we need appropriate tools to assess the physical and economic impacts of ARs in a warmer climate.

We have presented an integrated modeling tool that tracks an AR from its atmospheric development to the economic impacts related to inundation and flooding. We have used this tool to understand how the ARs and their impact could change in a warmer climate using a PGW approach. As argued by Hazeleger et al. (2015), this type of approach is particularly useful for the affected communities because it uses high-resolution models to simulate an extreme hydrologic event that occurred in the past and that the community can remember. The method is flexible enough to tailor the projections to a narrative; in this case, how would this extreme event change in a warmer climate? Furthermore, the method takes into account three types of uncertainty: (a) uncertainty in future radiative forcing, (b) uncertainty in the climate system response to this radiative forcing and (c) uncertainty in the hydrologic response of the system. In this way, we provide the community with a range of uncertainty of possible future conditions.

In our application to the December 2007 AR flooding event over the Chehalis River basin, we found that while there is a clear intensification of AR specific humidity and integrated vapor transport for both the lower and upper PGW scenarios, these changes do not translate into generalized increases in precipitation throughout the basin due to spatially heterogeneous changes in relative humidity and water vapor mixing ratio. For this reason, some parts of the basin receive more precipitation, while others receive less. These changes in precipitation translate into amplified changes in sub-basin runoff (in terms of percent change in water mass). But, because the upper basin runoff increases substantially, the streamflow along most of the Chehalis main stem increases in the warming scenarios. Interestingly, this event was so large that even in the control simulation most of the inundated area was occupied. As a consequence, while the PGW simulation resulted in significant changes in inundation depth, changes in the inundated area were minor. However, these changes in flood depth resulted in economic losses due to stock damages that ranged between −1 and 171 %, while losses in local production and trade within the three impacted counties were between 14 and 619 % (depending on the affected county, PGW scenario and hydrologic model). The economy outside of these three counties actually benefited from the event as it provided the entirety of the reconstruction efforts after the flood. Because the 2007 event was so rare, we also offer a simplified way to estimate the economic losses associated with floods with a shorter return period and calculate changes in expected annual losses.

The meteorology and hydrology combined with public policy and mitigation cost–benefit considerations will remain a difficult challenge in the future for the Chehalis basin. Flooding potential may need to be reconsidered in light of possible changes in atmospheric rivers in a warmer climate. Our integrated modeling tool provides communities in the Chehalis region with a range of possible future physical and socioeconomic impacts associated with AR flooding. The framework takes into consideration several important sources of uncertainty. It can be applied to other intense flooding events that perhaps affected other parts of the basin. Furthermore, the tool can be modified to understand different future scenarios, including the failure of hydraulic structures and changes in land use and land cover. In this way, communities in the region will be better prepared to mitigate the losses and improve disaster relief efforts associated with likely changes in precipitation and flooding that a warmer climate will bring.