What's streamflow got to do with it? A probabilistic simulation of the competing oceanographic and fluvial processes driving extreme along-river water levels

Extreme water levels generating flooding in estuarine and coastal environments are often driven by compound events, where many individual processes such as waves, storm surge, streamflow, and tides coincide. Despite this, extreme water levels are typically modeled in isolated open-coast or estuarine environments, potentially mischaracterizing the true risk of flooding facing coastal communities. This paper explores the variability of extreme water levels near the tribal community of La Push, within the Quileute Indian Reservation on the Washington state coast, where a river signal is apparent in tide gauge measurements during high-discharge events. To estimate the influence of multiple forcings on high water levels a hybrid modeling framework is developed, where probabilistic simulations of joint still water level and river discharge occurrences are merged with a hydraulic model that simulates along-river water levels. This methodology produces along-river water levels from thousands of combinations of events not necessarily captured in the observational records. We show that the 100-year still water level event and the 100-year discharge event do not always produce the 100-year along-river water level. Furthermore, along specific sections of river, both still water level and discharge are necessary for producing the 100-year along-river water level. Understanding the relative forcing driving extreme water levels along an ocean-to-river gradient will help communities within inlets better understand their risk to the compounding impacts of various environmental forcing, which is important for increasing their resilience to future flooding events.

1 Introduction

Coincident or compound events are a combination of physical processes in which the individual variables may or may not be extreme; however, the result is an extreme event with a significant impact (Zscheischler et al., 2018; Bevacqua et al., 2017; Wahl et al., 2015; Leonard et al., 2014). Flooding is often caused by compound events, where multiple factors impact both open-coast and estuarine environments. Storm events, for example, often generate concurrently large waves, heavy precipitation driving increased streamflow, and high storm surges, making the relative contribution of the actual drivers of extreme water levels difficult to interpret. Studies at the global (e.g., Ward et al., 2018), national (e.g., Wahl et al., 2015; Svensson and Jones, 2002; Zheng et al., 2013) and regional scale (e.g., Odigie and Warrick, 2017; Moftakhari et al., 2017) have evaluated the likelihood for variables such as high river flow and precipitation to occur during high coastal water levels, demonstrating that dependencies often exist between these individual processes.

Around river mouths, the elevation of the water level measured by tide gauges, or the still water level (SWL), varies depending on the mean sea level, tidal stage, and the nontidal residual contributors which may include the following forcings: storm surge, seasonally induced thermal expansion (Tsimplis and Woodworth, 1994), the geostrophic effects of currents (Chelton and Enfield, 1986), wave setup (Sweet et al., 2015; Vetter et al., 2010), and river discharge. Most commonly, estimates of nontidal residuals are determined by subtracting predicted tides from the measured water levels. However, residuals computed in this way often contain artifacts of the subtraction process from phase shifts in the tidal signal and/or timing errors (Horsburgh and Wilson, 2007). Another approach for extracting the nontidal residual is through the skew surge, which is the absolute difference between the maximum observed water level and the predicted tidal high water (de Vries et al., 1995; Williams et al., 2016; Mawdsley and Haigh, 2016). While this methodology removes the influence of tide–surge interaction from the nontidal residual magnitude, it does not differentiate between the many factors contributing to the water level, an important step for distinguishing when and why high water, and thus flooding, is likely to occur.

Hydrodynamic and hydraulic models have recently been used in attempts to quantify the relative importance of river- and ocean-forced water levels to flooding. The nonlinear coupling of wind- and pressure-driven storm surge, tides, wave-driven setup, and riverine flows has been found to be a vital contributor to overall water level elevation (Bunya et al., 2010). Furthermore, river discharge is often found to interact nonlinearly with storm surge (Bilskie and Hagen, 2018), exacerbating the impacts of coastal flooding (Olbert et al., 2017), which suggests that the extent or magnitude of flooding is often underpredicted when both river and oceanic processes are not modeled (Bilskie and Hagen, 2018; Kumbier et al., 2018; Chen and Liu, 2014). The computational demand of two- and three-dimensional hydrodynamic models, however, typically precludes a large amount of events to be examined. Therefore, while accurately modeling the physics of the combined forcings, researchers taking this approach are often limited to modeling only a select number of boundary conditions. On the other hand, statistical models allow for the investigation of compound water levels through the simulation of combinations of dependent events which may not have been physically realized in observational records (Bevacqua et al., 2017; van den Hurk et al., 2015). In addition, researchers have recently begun to generate hybrid models that link statistical and physical modeling approaches for understanding compound flood events (Moftakhari et al., 2019; Couasnon et al., 2018). Similar to the results solely from hydrodynamic and hydraulic models, statistical and hybrid modeling strategies show that simplifications of the dependence between multiple forcings may lead to an underestimation of flood risk.

This study explores the influence of oceanographic and riverine processes on extreme water levels along a coastal river where there is a substantial river signal recorded in the tide gauge. In order to better understand the river- and ocean-forced water levels at this location, a hybrid methodology is developed for linking statistical simulations of ocean and river boundary conditions with a hydraulic model that simulates along-river water levels. First, river-influenced water levels are defined and removed from SWLs. Then, both river discharge and river-influenced water levels are incorporated into a nonstationary, probabilistic total water level model, which allows for multiple synthetic representations of joint ocean and riverine processes that may not have occurred in the relatively short observational records. Next, a one-dimensional hydraulic model is used to simulate water surface elevations along a 10 km stretch of river. Surrogate models are generated from the hydraulic model simulations and used to extract along-river water levels for each probabilistic joint occurrence of SWL and river discharge in a computationally efficient manner. Rather than determining the along-river return level from an equivalent return level forcing (e.g., the 100-year discharge event drives the 100-year water level), spatially varying along-river return levels are extracted and matched to the driving boundary conditions. This technique allows for a spatially explicit analysis of the ocean and river conditions generating extreme water levels.

The following sections describe the study area, present the hybrid modeling framework linking oceanographic and riverine systems, and evaluate the compounding drivers of along-river extreme water levels.

2 Study area

The Quillayute River is located in Washington state along the US west coast and drains approximately 1630 km² of the northwestern Olympic Peninsula into the Pacific Ocean (Czuba et al., 2010). The Quillayute River is approximately 10 km long; is formed by the confluence of the Bogachiel and Sol Duc rivers (Fig. 1); and enters the Pacific Ocean at La Push, Washington, home to the Quileute Tribe. The Quileute Indian Reservation is approximately 4 km² and the majority of community infrastructure sits at the river mouth, bordering the river and open coast. The Quileute Harbor Marina is also situated just inside the river mouth and is the only port between Neah Bay and Westport, Washington. Rialto Spit, which connects Rialto Beach to Little James Island, contains a rocky revetment which protects the marina and the community from ocean and storm wave impact.

Figure 1 Map of study area (a), which is denoted on the regional map (b) in the black box. The La Push tide gauge is represented as a red square while other regional tide gauges are represented as blue squares. The Calawah and Sol Duc River gauges are represented as black triangles, and USGS measurement sites from the May 2010 survey (see information in the Supplement) are depicted as yellow circles. Approximate river kilometers are denoted as black crosses on the study area map.

The Quillayute River is a natural, unstabilized river that is relatively straight at the confluence of the Bogachiel and Sol Duc rivers and increases in sinuosity moving towards the river mouth. Channel-bed materials are coarse (gravel and cobble) in the free-flowing channels and dominated by sand in the small estuary (Czuba et al., 2010). Upstream of river kilometer 3 there are numerous point bars and bends in the river. Between river kilometer 1.5 and 3, the Quillayute River is braided with several side channels, usually containing woody debris (Czuba et al., 2010). The channel is straight near the river mouth and is confined by the Rialto Spit revetment before draining into the Pacific Ocean.

The oceanic climate of the coastal Pacific Northwest (PNW) is cool and wet with a small range in temperature variation and the majority of rainfall between October and May. Four river basins, the Sol Duc, Bogachiel, Calawah, and Dickey rivers, feed into the Quillayute River and comprise the majority of the watershed. Streamflow in the region is primarily from storm-derived rainfall in the winter and snowmelt during the spring and summer (WRCC, 2017).

Oceanographically driven SWLs are generally comprised of mean sea level, tides, and nontidal residuals, which include storm surge. Regional variations in shelf bathymetry, shoreline orientation, storm tracks (Graham and Diaz, 2001), seasonality (Komar et al., 2011), and winds drive differences in storm surge along the US west coast. However, the US west coast's narrow continental shelf, in relation to broad-shelved systems, controls the magnitude of storm surge, which is rarely larger than 1 m (Bromirski et al., 2017; Allan et al., 2011). The PNW is also influenced by a unique interannual climate variability due to the El Niño–Southern Oscillation. During El Niño years, the PNW experiences increased water levels for months at a time, along with changes in the frequency and intensity of storm systems (Komar et al., 2011; Allan and Komar, 2002). In the PNW, tides are micro- and mesotidal, and at La Push the tidal range is mixed, predominantly semidiurnal, with a mean range of 1.95 m and a great diurnal range of 2.58 m (https://tidesandcurrents.noaa.gov/datums.html?id=9442396, last access: October 2017).

Global rise in sea level and local changes in vertical land motions result in significant longshore variations of relative sea level along the Washington coastline. The northern Washington coast is experiencing relative sea level rates of $-\mathrm{1.85}\pm \mathrm{0.42}$ mm year⁻¹ due to a rising coastline, while relative sea level in Willapa Bay in southern Washington is 0.94±2.14 mm year⁻¹ (Komar et al., 2011). Tide gauge records at La Push are too short to calculate robust trends in sea level; however, sea level is likely rising in this location rather than falling, partly due to local land subsidence (Miller et al., 2018).

Figure 2 (a) The joint relationship between storm surge and wave height for La Push, Washington (black); and Westport, Washington (pink). Example storm surge and discharge relationship at (b) La Push and (c) Westport, Washington.

3 Data

Observational records in the region are generally sparse; one tide gauge exists in the marina near the river mouth and only two of the four rivers which feed into the Quillayute watershed are gauged (Fig. 1). The Sol Duc River gauge (WA Dept of Ecology 12A070) is located 12 km upriver from the Quillayute River and measures hourly discharge and stage observations from 2005 to 2014. The second river gauge is located on the Calawah River (USGS 12043000), approximately 25 km upriver from the Quillayute River. The Calawah River flows into the Bogachiel River and has hourly discharge and stage measurements from 1989 to 2016. The hourly record of discharge measurements from the Sol Duc River is 100 % complete, while the Calawah River is 99 % complete. An area scaling watershed analysis (Gianfagna et al., 2015) is undertaken to rectify the discharge by the amount of ungauged watershed. The watershed delineation shows that the Bogachiel, Calawah, Sol Duc, and Dickey rivers account for 24 %, 22 %, 37 %, and 17 % of the total Quillayute River watershed area, respectively. Noting the similar watershed characteristics and proportional watershed areas, the contribution of the Bogachiel River is estimated by scaling the Calawah River discharge measurements by a factor of 2.09. This scaling factor for estimating Bogachiel River discharge is validated by comparing to eight discharge point measurements taken during a US Geological Survey (USGS) survey in 2010 (see Supplement). Discharge for the Quillayute River is estimated by adding together discharge from the Sol Duc and Bogachiel rivers.

Hourly measured SWLs at the La Push tide gauge (NOAA station 9442396, 2004–2016) relative to Mean Lower Low Water (MLLW) are downloaded, transformed into NAVD88, and decomposed into mean sea level (η_MSL), tide (η_A), and nontidal residual (η_NTR). The η_NTR is further decomposed into monthly mean sea level anomalies (η_MMSLA), seasonality (η_SE), and storm surge (η_SS), using methods described in Serafin et al. (2017). Peak η_SS events at La Push are found to be the highest on record compared to all US west coast tide gauge stations (Serafin et al., 2017). Upon further investigation of the η_SS record, a large portion of extreme η_SS events occur during low-wave events (Fig. 2a) and high-river-discharge events (Fig. 2b). This is inconsistent with η_SS in Westport, Washington (Fig. 2a and c), just south of La Push, and with other tide gauges along the US west coast (not shown). It is therefore hypothesized that the anomalously large signal seen in the η_SS is river-induced.

To further investigate the anomalously large η_SS at the La Push tide gauge, the hydrodynamic model ADvanced CIRCculation (ADCIRC, Luettich Jr. et al., 1992) and Simulating Waves Nearshore (SWAN, Zijlema, 2010) model (ADCSWAN; Dietrich et al., 2011) is used to simulate water levels at the tide gauge during a storm event corresponding with the peak river discharge on record occurring on 8 January 2009. ADCIRC is run in 2-D depth-integrated barotropic mode which performs well for calculating water surface elevations during storm events (Weaver and Luettich Jr., 2010). SWAN is run in nonstationary mode on an unstructured grid, allowing for tight coupling to ADCIRC. The model is run with two forcing implementations: one including a full forcing of waves, wind, pressure, streamflow, sea level anomalies, seasonality, and tides and one including only streamflow and tides. Once the river-influenced water level is validated, it is removed from the η_SS signal and saved as a sixth geophysical variable (η_Ri; see Supplement for removal technique).

Because of the short length of the La Push tide gauge record, decomposed water levels from the La Push tide gauge are merged with decomposed water levels from the Toke Point tide gauge (NOAA station 9440910) to create a combined water level record with a length of 36 years. Details of this methodology are explained in the corresponding Supplement, as well as in Serafin et al. (2019). Once the two tide gauges are merged, the combined hourly tide gauge record extends from 1980 to 2016 and is 97 % complete.

4 Methods

Return level flood magnitudes, such as the 100-year event, are typically assumed to be driven by a specific forcing event, such as the 100-year rainfall or storm surge. However, for processes driven by multiple dimensions, different sizes and combinations of forcing conditions could potentially generate extreme flood magnitudes. To explore the role of compounding forcings in generating extreme water levels, a hybrid modeling framework is developed by merging a hydraulic model simulating river flow with probabilistic simulations of jointly occurring boundary conditions, in this case SWL and river discharge (Fig. 3). Statistical simulations allow for long, synthetic records of joint forcings that may not have occurred in the short observational records but are physically capable of co-occurring. Modeling all of the statistically simulated boundary conditions in a hydraulic model to output along-river water levels would be prohibitively expensive. As an alternative to time-consuming simulations, surrogate models (Razavi et al., 2012) are developed to approximate the response of a hydraulic model simulation at each along-river location. This technique allows for the analysis of along-river water levels driven by a variety of boundary conditions. Long synthetic records on the order of 500 years allow for the direct empirical extraction of water level return levels rather than an extrapolation from historic observational forcing conditions. In addition, the large sample space of simulated variables permits a comparison of event-based return levels, where the 100-year water level is determined by the 100-year forcing, to response-based return levels, where the 100-year water level is derived and then mapped to its respective forcing conditions. This novel framework is flexible for input of any statistical or hydraulic model. In this application, we use the Serafin and Ruggiero (2014) full simulation total water level model and the US Army Corps of Engineers' (USACE) Hydrologic Engineering Center's River Analysis System (HEC-RAS; Brunner, 2016), which are described in more detail below.

Figure 3 Schematic of hybrid statistical–physical modeling technique. Models are portrayed as squares, while circles portray model outputs.

4.1 Probabilistic simulations of boundary conditions

The nonstationary, probabilistic simulation model of Serafin and Ruggiero (2014) (hereinafter SR14) was developed to produce synthetic time series of daily maximum total water levels (TWLs), the combination of waves, tides, and nontidal residuals, on open-coast sandy beaches. SR14 simulates the individual components of the TWL in a Monte Carlo sense, while appropriately accounting for any dependencies existing between the variables. This modeling technique is able to include nonstationary processes influencing extreme and nonextreme events, such as seasonality, climate variability, and trends in wave heights and water levels. SR14 outputs a number of synthetic records of all variables driving TWLs that produce alternate, but physically plausible, combinations of waves and water levels along an identified stretch of coastline. This technique is flexible to allow for both (i) the simulation of the present-day climate for computing robust statistics on extreme TWL events and (ii) the simulation of future climates and their impact on extreme TWLs. Because SR14 was developed for use in open-coast environments, it does not include a procedure for simulating estimates of river discharge, which is present in the local tide gauge at the La Push study site. SR14 is therefore modified to produce synthetic time series of river discharge as well as a river-induced water level.

Figure 4 (a) Joint relationship between wave height (H_s) and discharge (Q) for the observational record (black) and one example 500-year simulation (red). (b) Seasonal model fit for the probabilistic simulation of the Sol Duc River Q in relation to the Bogachiel River Q. The inset displays the model fits for discharge less than 100 m³ s⁻¹.

High-discharge events on the two gauged rivers in the watershed, the Sol Duc and Calawah rivers, tend to occur within hours of peak wave events recorded in offshore wave buoy records and water level events recorded in the tide gauge data. Due to the interrelated nature of these forcings, daily maximum estimates of Calawah River discharge (Q_C) are compared to all variables simulated in the SR14 model (e.g., wave height, η_SS, η_NTR, η_MMSLA) to capture any dependencies between these processes. The variable with the highest monthly correlation to Q_C is wave height (H_s). Extreme Q_C events are simulated using a bivariate logistic model, which is the same technique used to simulate η_SS. The bivariate logistic model preserves the dependency and frequency of occurrence of joint H_s−Q events in extreme and nonextreme space. This technique generates a synthetic record of Q_C that is seasonally varying, related to larger-scale climate variability through wave height (essentially as a proxy for storms), and carries the same dependency between variables as the observational record (Fig. 4a). Q_C is then multiplied by 2.09 to represent inflow from both the Bogachiel and Calawah rivers.

Discharge measurements at the Sol Duc River are highly correlated with the discharge measurements at the Calawah River (ρ=0.9, τ=0.83); thus Sol Duc River discharge (Q_SD) is modeled based on a relationship with the scaled Q_C, representing the Bogachiel River (Q_B). Estimates of Q_SD are related to Q_B during the summer and winter seasons. First, daily maximum Q is split into summer (May, June, July, August, September, and October) and winter (January, February, March, April, November, and December) seasons. Next, models are fit to the joint relationship between the Q_SD and Q_B each season, such that for the summer season,

$$\begin{array}{}\text{(1)}& Q_{\mathrm{SD}}=\mathrm{1.186}{Q}_{\mathrm{B}}+\mathrm{0.226}\end{array}$$

is used when Q_B falls between 0 and 10 m³ s⁻¹, and

$$\begin{array}{}\text{(2)}& Q_{\mathrm{SD}}=-\mathrm{1.0}\times {\mathrm{10}}^{-\mathrm{4}}{Q}_{\mathrm{B}}^{\mathrm{2}}+\mathrm{0.38}{Q}_{\mathrm{B}}+\mathrm{14.07}\end{array}$$

is used when Q_B falls between 10 and 700 m³ s⁻¹ (Fig. 4b). When Q_B is greater than 700 m³ s⁻¹, Q_SD is determined using

$$\begin{array}{}\text{(3)}& Q_{\mathrm{SD}}=\mathrm{0.216}{Q}_{\mathrm{B}}+\mathrm{61.25}.\end{array}$$

For the winter season,

$$\begin{array}{}\text{(4)}& Q_{\mathrm{SD}}=\mathrm{0.816}{Q}_{\mathrm{B}}+\mathrm{1.168}\end{array}$$

is used when Q_B falls between 0 and 25 m³ s⁻¹, and

$$\begin{array}{}\text{(5)}& Q_{\mathrm{SD}}=-\mathrm{1.0}\times {\mathrm{10}}^{-\mathrm{4}}{Q}_{\mathrm{B}}^{\mathrm{2}}+\mathrm{0.46}{Q}_{\mathrm{B}}+\mathrm{16.11}\end{array}$$

is used when Q_B falls between 25 and 2300 m³ s⁻¹ (Fig. 4b). When Q_B is greater than 2300 m³ s⁻¹, Q_SD is determined using

$$\begin{array}{}\text{(6)}& Q_{\mathrm{SD}}=\mathrm{0.075}{Q}_{\mathrm{B}}+\mathrm{500.42}.\end{array}$$

Summer and winter Q_B is binned and residuals of Q_SD from the above model fits are generated. Normal distributions are fit to Q_SD residuals in each bin, except for low bins (less than 25 m³ s⁻¹) where residuals are fit to exponential distributions. Q_SD is then directly related to simulated estimates of Q_B; Q_SD is first determined by fitting the prescribed model to each estimate of Q_B, and then a random sample is taken from the residuals per binned Q_B and added to the model. This technique captures the joint peaks of the river systems visible in the observed dataset, while allowing for variability between the simulated estimates (Fig. 4b).

Figure 5 (a) The relationship between the river-influenced water level (η_Ri) and Bogachiel River discharge on a log-linear scale. The solid black line represents the model fit to the observational records (black dots). (b) The percentage of time η_Ri occurs in the record during a specific Q_B. In both panels, black represents the observational record and red represents one example 500-year simulation.

4.1.1 Modeling the river-induced water level

At tide gauges along the US west coast, the maximum daily SWL generally occurs during, or close to, the daily high tide (Serafin and Ruggiero, 2014; Serafin et al., 2017). Modeling peaks in η_Ri that occur during low tide would therefore erroneously increase simulated estimates of the SWL occurring during high tide. Thus, instances of η_Ri occurring approximately during high tide are retained and all other η_Ri peaks are discarded, resulting in 155 η_Ri events.

Synthetic estimates of η_Ri are developed by relating Q_B and η_Ri. This relationship is modeled using

$$\begin{array}{}\text{(7)}& {\mathit{\eta }}_{\mathrm{Ri}}=\mathrm{0.039}{Q}_{\mathrm{B}}+\mathrm{0.854}\times {\mathrm{10}}^{-\mathrm{3}}\end{array}$$

when Q_B is below 190 m³ s⁻¹ and

$$\begin{array}{}\text{(8)}& {\mathit{\eta }}_{\mathrm{Ri}}=\mathrm{0.093}{Q}_{\mathrm{B}}+\mathrm{0.284}\times {\mathrm{10}}^{-\mathrm{3}}\end{array}$$

when Q_B is above 190 m³ s⁻¹ (Fig. 5a). Next, coarse bins ranging from 100 to 4000 m³ s⁻¹ are created and the standard deviation (σ) of η_Ri within each bin is saved. For bins that contain less than 10 observations, observations from the previous bins are included until there are more than 10 observations per bin for σ calculations. Finally, a two-point running average is used to smooth σ from each bin to ensure continuous transitions and to avoid the edge effects from binning a sparse dataset.

There are times of high Q_B without a distinguishable η_Ri in the tide gauge record; thus a model is also developed to simulate the frequency of occurrence of η_Ri during daily maximum SWLs. The frequency of occurrence of η_Ri is defined as the percentage of time η_Ri occurs in the observational record, which is less than 10 % of the time when Q_B is less than 840 m³ s⁻¹ and 10 %–25 % of the time when Q_B is between 840 and 2090 m³ s⁻¹ (Fig. 5b). For Q_B greater than 2090 m³ s⁻¹, η_Ri occurs approximately 50 % of the time. The frequency of occurrence of η_Ri is modeled using a best-fit cubic function, where the frequency of occurrence is a function of Q_B based on the percentage of time the values have occurred in the record. Because there are no events greater than 2500 m³ s⁻¹ on record, we represent the percentage of occurrence over this value as 100 % (Fig. 5b).

Once Q_B is simulated using SR14, η_Ri is simulated for every day in time by selecting the synthetic daily estimate of Q_B and randomly sampling from a normal distribution for each Q_B bin, where μ is the regression model and σ is the standard deviation from each bin (Fig. 5a). The frequency-of-occurrence model is then used to select the correct proportion of η_Ri events to retain for each synthetic simulation. These techniques capture both the spread of η_Ri related to Q_B and the percentage of time of occurrence (Fig. 5).

Figure 6 Digital elevation model (DEM) used for the HEC-RAS simulations of the Quillayute River. HEC-RAS cross sections are depicted as grey lines. Approximate river kilometer and the location of the tide gauge are depicted as diamonds and a square, respectively.

4.2 Hydraulic model for along-river water levels

While a variety of hydraulic models can be used for determining the elevation of along-river water levels, we employ the Hydraulic Engineering Center's River Analysis System (Brunner, 2016). HEC-RAS is used to estimate water surface elevations in rivers and streams in both steady and unsteady flow and under subcritical, supercritical, and mixed flow regimes (Goodell, 2014). HEC-RAS has been previously used to model water surfaces for a range of applications including, but not limited to, floodplain mapping (Yang et al., 2006), flood forecasting (Saleh et al., 2017), dam breaching (Butt et al., 2013), and flood inundation (Horritt and Bates, 2002). HEC-RAS computes water levels by solving the 1-D energy equation with an iterative procedure, termed the step method, from one cross section to the next (Brunner, 2016). For subcritical flows, the step procedure is carried out moving upstream; computations begin at the downstream boundary of the river and the water surface elevation at an upstream cross section is iteratively estimated until a balanced water surface is obtained. Energy losses between cross sections are comprised of a frictional loss via the Manning equation and a contraction/expansion loss via a coefficient multiplied by the change in velocity head (see Brunner, 2016, for more details).

In this application, HEC-RAS is used to model 1-D water levels under gradually varied, steady-flow conditions at transects along the Quillayute River. While a simplification of flood processes, the 1-D application is commonly used to create flood hazard maps. A detailed digital elevation model (DEM) is developed for the river network, including bathymetry and topography for the floodplains of interest (Fig. 6). Model domain boundary conditions are chosen as the SWL at the tide gauge (m; downstream boundary) and river discharge from the Sol Duc and Bogachiel rivers (m³ s⁻¹; upstream boundary). The HEC-RAS model is validated using water surface measurements from a 2010 survey. Details of the HEC-RAS model validation and calibration procedures are documented in the Supplement.

4.3 Hybrid statistical–physical modeling

The modified simulation technique of SR14 is used to produce 70 500-year-long synthetic records representing present-day climate for the time period of 1980–2016 of daily maximum SWL and Q for both the Sol Duc and Bogachiel rivers. Rather than run the ∼13 million simulated conditions through a numerical model, a limited set of joint boundary conditions of SWL and Q (at the Bogachiel and Sol Duc rivers) are run through HEC-RAS, outputting the elevation of the along-river water level at each HEC-RAS transect. Surrogate models are generated from the HEC-RAS runs for each transect using a scattered linear interpolation of the 3-D surface of boundary conditions. The number of combinations of SWL and Q used to develop the surrogate models are chosen to minimize interpolation errors during validation runs. A daily estimate of water level elevation at each transect is produced by inputting all daily maximum SWL and Q conditions into the surrogate models, which efficiently extract along-river water levels for any set of SWL and Q inputs. Using the count-back method, where for example the fifth largest event for each synthetic record would be the 100-year event, water level return levels are extracted for all 70 500-year synthetic records for the (1) along-river water levels at each transect, (2) SWLs, and (3) Q. This methodology provides both an estimate of the return level magnitude (e.g., the average of the 70 100-year events) and the uncertainty around that magnitude (e.g., the distribution of the 70 100-year events). It also provides a technique to compare the response-based return level (e.g., the 100-year water level) to the event-based return level (e.g., the water level driven by the 100-year SWL or 100-year Q event).

5 Results

The following section first validates the presence of a river-induced water level within the tide gauge signal and then demonstrates the effectiveness of the surrogate models in representing along-river water levels for unmodeled HEC-RAS boundary conditions. Next, the spatial and temporal variability of the magnitude of along-river water levels and their driving conditions are examined. Finally, low-probability water levels, like the 100-year event, are extracted from the simulated records of along-river water levels, and the dominant drivers are evaluated.

5.1 River-induced water level validation

Results from ADCSWAN modeling of the 8 January 2009 storm event show that the simulation including only river discharge and tides is nearly able to recreate the measured peak η_SS signal at the tide gauge (Fig. 7a). The addition of wind, pressure, waves, sea level anomalies, and seasonality is found to have minimal impact on the peak observed η_SS. Furthermore, the maximum η_SS occurs during low tide (Fig. 7b), which indicates a potential relationship between water surface elevation, tidal level, and river discharge. While the ADCSWAN runs only explore one instance of this phenomenon, it provides physics-based evidence that anomalously high η_SS at the La Push tide gauge is likely being driven by large discharge events.

Figure 7 Resulting storm surge (a) and still water level (b) at the La Push tide gauge modeled using ADCIRC for a simulation including full forcing (red) and a simulation including only discharge and tides (blue) compared to the observational record (black). The ADCIRC simulation was run for the maximum discharge event on record occurring on 8 January 2009.

5.2 Surrogate model validation

A number of validation scenarios are modeled in HEC-RAS to determine whether the combinations of Q and SWL boundary conditions used to develop the surrogate models represent a large enough sample space of forcing conditions for the interpolation of along-river water levels. The validation scenarios are chosen to cross through both HEC-RAS modeled and unmodeled conditions (Fig. 8a). Across all validation scenarios, the average root mean square error (RMSE) between the HEC-RAS directly modeled and the surrogate model-interpolated water levels is 1 cm. Only 1.5 % of the validation scenarios have a bias greater than 10 cm, and the largest RMSE at any transect is 20 cm across all scenarios (Fig. 9). The validation scenario with the worst performance occurs during high Q_B and low Q_SD paired with low-SWL events. However, even during the worst-performing case, the difference between the HEC-RAS directly modeled water level and the surrogate model-interpolated water level is small (Fig. 8b). The main research focus here is extreme water levels, and the conditions driving low-probability return level events rarely fall around the scenarios with the highest bias.

Figure 8 (a) Modeled HEC-RAS Q boundary conditions used to generate the surrogate models (red-dotted lines) compared to the simulated conditions used for surrogate model validation (green dots). The black dots represent the observational daily max conditions, while the colored circles represent the worst performing of the validation tests. The red and blue colored circles represent the scenarios where the interpolated water surface had a bias of over 10 cm lower than the model. (b) Example along-river water level for the worst-performing condition in the validation tests.

Figure 9 (a) Average root mean square error (RMSE) and (b) bias for all 197 discharge validation scenarios across 4 out of the 15 SWL scenarios. The worst-performing model is discharge scenario 153.

5.3 Hybrid modeling of along-river water levels

5.3.1 Temporal variability

Seasonal variability exists in the elevation of along-river water levels. The highest elevation water level occurs during the winter (here defined as December, January, February), while the lowest elevation water level occurs during the spring (March, April, May) (Fig. 10a). The spring along-river water level is on average 0.50 m lower than the winter along-river water level, 0.33 m lower than the fall (September, October, November) along-river water level, and 0.03 m lower than the summer (June, July, August) along-river water level (Fig. 10b). The difference between seasonal along-river water levels is nonlinear upstream, and certain sections of the river have larger changes in elevation between months (Fig. 10b). However, this variation becomes relatively linear downstream of river kilometer 3.

Figure 10 (a) Variability of along-river water levels averaged over spring (MAM), summer (JJA), fall (SON), and winter (DJF). (b) The difference between the spring and summer, fall, and winter along-river water levels.

The seasonal variability of the along-river water level is driven by the seasonality of the forcings, which are well represented in the simulations compared to the observations (Fig. 11). The monthly median SWLs and η_NTR's are higher in the winter than in the summer (Fig. 11a and b). This cyclical variability is also depicted in the monthly median river discharge from the Quillayute River (combined Sol Duc and Bogachiel Q) and is approximately 200 m³ s⁻¹ higher in winter months than summer months (Fig. 11c). The 98th percentile values of SWL, η_NTR, and Q have a similar seasonal variability to the median conditions (Fig. 11d–f).

Figure 11 (a–c) Observational (black) and simulated (red) monthly median still water level (SWL), nontidal residual (eta_NTR), and discharge (Q). (d–f) Observational (black) and simulated (red) monthly 98th percentile of the SWL, η_NTR, and Q. Red shading indicates the bound from each simulation.

5.3.2 Spatial variability

The large number of joint SWL and Q conditions allows for the direct extraction of water level return levels and the corresponding univariate or multivariate drivers along each HEC-RAS transect. The magnitude of the 100, 25, 10, and annual return level water levels is between 3 and 17 m (NAVD88, Fig. 12a). While the peaks in return level events occur at similar locations, the difference between return level events varies spatially moving upriver. For example, at river kilometer 1, the difference between the average (of all simulations) annual and 100-year event is approximately 0.9 m, whereas at river kilometer 8 and upstream, the difference between these two events is closer to 2 m (Fig. 12b).

Figure 12 (a) The water level return level at each transect for all 70 probabilistic simulations. Each return level event displays the average of the simulations (solid line) as well as the range around the average (shaded). (b) The along-river difference between the annual and 100-year event, averaged over 70 simulations.

The dominant forcing conditions driving water level return levels vary along-river. At the river mouth, the annual water level event (e.g., the event that is expected every year) in each simulation occurs during Q ranging from 40 to 2600 m³ s⁻¹ and SWLs around 3.3 m, which corresponds with the annual SWL event (Fig. 13a). Moving upstream to river kilometer 1 and 2, the annual water level event is driven by both high SWL occurring during low Q and low SWL occurring during high Q. At river kilometer 4, the annual water level event occurs during the annual Q event coincident with SWLs that range from 1.8 to 3.9 m (Fig. 13a). These results are similar, albeit events are larger magnitude, for the 100-year water level event. Downstream 100-year water levels are driven by SWLs, upstream 100-year water levels are driven by Q, and the 100-year water level between kilometer 1 and 2 is driven by different combinations of high- and low-SWL and Q events (Fig. 13b).

Figure 13 The individual Q or SWL condition driving the (a) annual and (b) 100-year water level event at specific along-river locations for each 70 500-year simulation. In both figures, the black lines represent the annual and 100-year return level magnitude for Q and SWL.

The relative importance of both oceanic and riverine forcing to extreme water levels emerges when averaging the magnitude of the drivers of the water level return levels at each transect from all 70 500-year-long simulations (Fig. 14). The magnitude of the average Q driving water level return levels gradually increases by approximately 1000 m³ s⁻¹ over river kilometer 0–2 and then is consistent from river kilometer 2 to 10 (Fig. 14a). Downstream, between river kilometer 0 and 0.25, the magnitude of the average SWL driving water level return levels is consistent and then gradually decreases over a 1 km zone (Fig. 14b).

Figure 14 The average forcing condition driving along-river return levels at each transect where (a) displays the Quillayute Q conditions and (b) displays the SWL conditions. The dashed lines depict the univariate forcing conditions, where the along-river return level is assumed to be driven by either Q or SWL. Red, orange, blue, and black lines represent the 100, 25, 10, and annual return level event. The grey shaded area represents a transition zone, where the water level is driven by a combination of SWL and Q events.

When comparing to water level return levels driven by a univariate forcing or event return level (e.g., along-river water levels modeled from the 100-year Q or SWL event), we find that the stretches of river driven by a consistent SWL or Q forcing approximate the univariate return level event. Therefore, the 100-year SWL does indeed cause the 100-year water level downstream, between river kilometer 0 and 0.25, while the 100-year Q event drives the 100-year water level upstream, between river kilometer 2 and 10 (dashed lines, Fig. 14). However, between river kilometer 0.25 and 1.75 a flood transition zone is present, where neither the SWL return level or the Q return level events drive the water level return level. This is consistent across all return level events.

6 Discussion

The hybrid model developed in this study, which combines statistical simulations with a physics-based model, provides an approach for probabilistically evaluating the conditions that drive extreme water levels not only in an open-coast setting, but also kilometers upriver. The ability to simulate millions of combinations of Q and SWL events allows for a robust estimate of resulting along-river water levels, which numerical models alone are unable to consider due to large computational expenses. While some of our modeling techniques are specific to this location, the overall framework for combining statistical and physics-based models is general enough for use in coastal locations throughout the globe where flooding arises from compounding processes.

The decomposition of the SWL into low- and high-frequency signals, including a river-influenced component, helps identify the importance of physical processes for generating high water levels across various regional settings. This is especially important in locations like the US west coast, where the steep, narrow continental shelf prevents wind- and pressure-driven storm surge from being overwhelmingly large (Allan et al., 2011). The influence of the river signal in the tide gauge is directly related to the setting of our study site. The estuary is relatively small and narrow with the river discharging directly into the ocean. This is dissimilar to other tide gauges in the region which are located in larger estuaries, situated away from river input. Estuaries typically exhibit wave, tide, or river-dominant morphology, based on the relative energy of each process (Dalrymple et al., 1992). The Quillayute River outlets directly to a high-wave-energy environment and has a small estuary volume compared to its river input volume. The steep catchment of the mountainous environment means a short response time for rainfall, therefore producing peak discharges temporally similar to peak storm-induced still water levels, allowing for interaction between the two. In contrast, water level elevations with large estuary volume compared to river discharge are less influenced by fluvial processes. Furthermore, a larger estuary may experience variability in the water surface elevation due to wave-induced setup and/or other local storm-induced processes (Cheng et al., 2014; Olabarrieta et al., 2011), which may further dampen the influence of a river signal.

This research confirms the presence of an oceanographic–fluvial transition zone, where traditional, univariate methodologies for defining return level events are insufficient for defining water level return levels. Between river kilometer 1 and 2, we find that a range of SWL and Q conditions drive all return level events, and water levels are driven by neither the univariate SWL or Q return level event. A similar flood zone transition was recently modeled numerically and, albeit for a single event, physically demonstrates the importance of including multiple variables to reproduce accurate flooding (Bilskie and Hagen, 2018). Thus, flood hazard assessments on systems with multivariate forcings may misrepresent water level elevations for low-probability events if only univariate variables are modeled. This has large implications for characterizing the risk to flooding, especially in the context of mapping flooding hazards. Furthermore, we show that return level water levels can occur over a range of combined extreme and nonextreme forcing in the flood transition zone. This illustrates that, in order to properly understand the impacts of compounding flooding, more than just design scenarios need to be considered for the proper assessment of risk.

Many of our results can be explained by dynamics that occur during interacting ocean and river flows. For example, a coincidence of high SWL and peak river discharge may induce blocking, where river-induced water levels are trapped upstream and either flood overbank or outlet to the ocean as the tide recedes (Kumbier et al., 2018; Chen and Liu, 2014). The outletting to the ocean as the tide recedes artificially inflates SWLs at the tide gauge, increasing water levels for days at a time and prolonging exposure to flooding. When subtracting a tide time series from this signal, storm surge would appear to be elevated at low tide. While the ADCSWAN simulation confirms the presence of this effect by matching the peak storm surge at low tide, our hybrid methodology only models steady-flow scenarios. Thus, with co-occurring daily maximum SWL and discharge, our model may miss certain dynamics important for flooding over unsteady conditions. Furthermore, interactions between storm surge and river discharge may increase the overall elevation of the residual (Maskell et al., 2013). While beyond the scope of our present study, these unsteady characteristics are important to consider in future research.

Because sea level rise, along with other changes to the climate, will exacerbate the compounding effects of flood drivers (Moftakhari et al., 2017; Wahl et al., 2015), it is also important to consider the impact of changes to processes driving flooding events in the future (Zscheischler et al., 2018). By 2100, the likely range of relative sea level rise in the La Push area is projected to be between 18 and 80 cm, considering vertical land motion and various emissions scenarios (Miller et al., 2018). The western Olympic Peninsula is projected to experience increased winter precipitation (Mote et al., 2013; Halofsky et al., 2011), which could subsequently increase either the frequency or intensity of high-Q events along the Quillayute River. While we have characterized the spatial variability in extreme water levels in the present day, there is a high likelihood changes in the future climate will shift the importance of these interacting processes.

7 Conclusions

This research illustrates the importance of considering a large number of forcing conditions to model compounding processes when evaluating extreme water levels. Here we find that, in coastal settings, river discharge can be an important driver of high water levels measured in a tide gauge. We also find that the univariate, event-based return level event, like the 100-year discharge, does not always match the response-based return level, like the 100-year water level. Furthermore, when processes compound, along-river return levels may be driven by events that are not considered extreme themselves. Probabilistic techniques allowing for the analysis of thousands to millions of combinations of events not captured in the observational record provide a characterization of where river, ocean, or the combination of the two may be important for generating extreme events.

Overall, the hybrid merging of a statistical and numerical model provides a methodology for better understanding the drivers of flooding along the length of a river. While our model does not actively resolve the physical interaction of river and oceanographic flow, it develops an approach for characterizing and extracting river-influenced water levels measured at tide gauges while robustly modeling the drivers of extreme along-river water levels. Understanding the dominant, spatially variable drivers of flooding events will help coastal communities better understand their risks, which is important for increasing resilience to future events.