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

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 2Diagram 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.

Download

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.

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

Download Print Version | Download XLSX

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 3Upper (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.

Download

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.

Download

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

Download

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.