Effect of Fluvial Discharges and Remote Non-Tidal Residuals on Compound Flood Forecasting in San Francisco Bay

: Accurate and timely ﬂood forecasts are critical for making emergency-response decisions regarding public safety, infrastructure operations, and resource allocation. One of the main challenges for coastal ﬂood forecasting systems is a lack of reliable forecast data of large-scale oceanic and watershed processes and the combined effects of multiple hazards, such as compound ﬂooding at river mouths. Offshore water level anomalies, known as remote Non-Tidal Residuals (NTRs), are caused by processes such as downwelling, offshore wind setup, and also driven by ocean-basin salinity and temperature changes, common along the west coast during El Niño events. Similarly, ﬂuvial discharges can contribute to extreme water levels in the coastal area, while they are dominated by large-scale watershed hydraulics. However, with the recent emergence of reliable large-scale forecast systems, coastal models now import the essential input data to forecast extreme water levels in the nearshore. Accordingly, we have developed Hydro-CoSMoS, a new coastal forecast model based on the USGS Coastal Storm Modeling System (CoSMoS) powered by the Delft3D San Francisco Bay and Delta community model. In this work, we studied the role of ﬂuvial discharges and remote NTRs on extreme water levels during a February 2019 storm by using Hydro-CoSMoS in hindcast mode. We simulated the storm with and without real-time ﬂuvial discharge data to study their effect on coastal water levels and ﬂooding extent, and highlight the importance of watershed forecast systems such as NOAA’s National Water Model (NWM). We also studied the effect of remote NTRs on coastal water levels in San Francisco Bay during the 2019 February storm by utilizing the data from a global ocean model (HYCOM). Our results showed that accurate forecasts of remote NTRs and ﬂuvial discharges can play a signiﬁcant role in predicting extreme water levels in San Francisco Bay. This pilot application in San Francisco Bay can serve as a basis for integrated coastal ﬂood modeling systems in complex coastal settings worldwide.


Introduction
Climate change is driving several processes that will lead to increased flood frequency in many parts of the world [1,2]. Global Sea-Level Rise (SLR) has been observed to be between 1.2 to 1.7 mm/yr Water Resources has funded the development of a state-of-the-art Advanced Quantitative Precipitation Information (AQPI) system to provide near-term precipitation and flooding forecasts using an integrated observation and modeling framework for the San Francisco Bay area [28,29]. The main goals of this collaborative project between NOAA, USGS, and several other stakeholders are to detect and provide advanced warning of storms, nowcast high-resolution precipitation with cutting-edge radar technology, and forecast watershed and coastal flooding up to 72 h in advance. As part of the AQPI framework, we developed an integrated modeling system, coupling HRRR and NWM to the hydrodynamic model of USGS Coastal Storm Modeling System (Hydro-CoSMoS, www.usgs.gov/cosmos) to achieve the ultimate goal of providing the public with an accurate coastal and inland flood forecast system.
Hydro-CoSMoS is unique in its inclusion of relatively small, steep, flashy watersheds. Much of the work on compound flooding has been focused on low-gradient coastal watersheds and regions susceptible to hurricanes [11]. The flooding in the San Francisco Bay Area on the coastal side has been more closely linked with El Niño, significant wave events, king tides, and coincident storm surge. Some of the larger oceanic storm events in the region have caused significant damage, such as the February 1998 El Niño-fueled storm [30,31]. Unfortunately, there is also a growing incidence of more common, nuisance-type flooding that causes minor damage, but substantial interruption to economic activity by impacting important transportation routes and other infrastructure. These nuisance flooding events are increasing in the region [7] and often require interventions to protect assets. An example of this is the flooding of roads and low-lying areas during the February 2019 storm in the San Francisco Bay area [32], caused by a combination of excessive rainfalls, low atmospheric pressure, and strong winds. Having both coastal and fluvial drivers of flooding made the February 2019 storm a perfect case study to validate Hydro-CoSMoS for extreme events and investigate its ability to capture compound flooding across the San Francisco Bay Area.
In this work, we explain how our coastal hydrodynamics model (Hydro-CoSMoS) imports meteorological, oceanic, and fluvial inputs from HRRR, HYCOM, and NWM in an operational setting. The focal point of this work is to show the effect of fluvial discharges and remote NTRs on compound flood forecasting by discussing our model's results from the February 2019 storm case study. Model performance, both with and without these larger scale inputs, is discussed to show and justify the need for incorporating these drivers in flood forecasting endeavors. The full details of the forecast system and hydrodynamic model setup will be discussed in parallel publications. Figure 1 shows the study area. Hydro-CoSMoS domain extends 20 km offshore of the Golden Gate, extending from Point Reyes in the north to Montara in the south. The model includes coastal zones of the nine counties within the San Francisco Bay area up to 10 m (NAVD88) elevation. Hydro-CoSMoS is powered by the Delft3D Flexible Mesh Suite (Delft3D Flexible Mesh (DFM); http://oss.deltares.nl/) platform and is an evolved version of the San Francisco Bay and Delta (SFBD) community model (http://www.d3d-baydelta.org/) which has been successfully used in several scientific and engineering projects [33][34][35][36]. The new version of the model is built together with Alameda County Flood Control and Water Conservation District for flood protection and forecasting purposes, with improved resolution using the latest topo-bathymetric sources [37]. We expanded upon the latest information on levees, berms, and other shoreline structures [38] with input from our local partners, especially on more upstream structures. Model resolution ranges from 2 km in the offshore to 200 m nearshore using an unstructured flexible mesh. It includes the four subembayments of Suisun Bay, San Pablo Bay, as well as South and Central San Francisco Bay. The new model applies the Delft3D-FM 1D/2D coupling, which includes the streams and rivers as 1D nodes. The San Francisco Bay Delta on the eastern extent of the model domain is primarily represented by 1D nodes. This 1D/2D approach shown in Figure 1 has resulted in more robust and efficient modeling than previous schematizations.

Study Area and Model Description
Tides at the mouth of the bay (San Francisco NOAA Station in Figure 1) are mixed, semi-diurnal, with a tidal range near 2 m. Tidal fluctuations extend up the Sacramento over 100 km from Suisun Bay. On average the freshwater inflow from the Delta is less than one percent of the spring tidal prism [39]. The result of this is that during normal conditions tidal currents are much stronger than freshwater driven flows in the 2D domain ( Figure 1). The model has been calibrated and validated for the time period of 1950-2019. The formal introduction of the schematization, including a more detailed description of the setup, results and subsequent extreme value analysis of this 70-year long hindcast of water levels in the San Francisco Bay Area will be published in a separate publication [40]. This study's results showed that the model captured the water levels with less than 10 cm mean averaged error (MAE) at all of the tidal stations depicted in Figure 1. The model performed slightly better in Central and South San Francisco Bay than Suisun Bay, where the averaged MAE over 70 years of the simulation was less than 5 cm in San Francisco and was slightly less than 9 cm in Port Chicago. This slightly worse performance in the vicinity of Suisun Bay is likely due to the role of salinity gradients being slightly more important in this region of freshwater confluence. From San Pablo Bay eastward gravitational circulation and other stratification driven processes have been found to be dynamically important [41]. However, in order to keep this model computationally inexpensive we must simplify to two dimensions, although the errors are slightly larger, this simplification has a small impact on water levels.  Astronomic tides are forced along the offshore boundaries derived from the TPXO8.0 global tide data [42]. We used measured tidal constituents at San Francisco (NOAA station 9414290) and calibrated the amplitudes and phases in a couple of iterations.

Remote Non-Tidal Residuals
In addition to tidal forcing, remote NTRs are forced along offshore boundaries. The forecast system obtains remote NTRs from the sea surface anomaly data provided by HYbrid Coordinate Ocean Model (HYCOM) [12]. HYCOM uses the Global Ocean Forecasting System (GOFS3.1) outputs to force its global model with a resolution of 0.08 • providing hindcast and forecast data for sea surface anomalies, ocean currents, temperature, and salinity from 1993 to present. HYCOM ocean prediction system provides 180 h of forecast data with 3-h time steps each day. HYCOM sea surface height forecast data does not include the effects of atmospheric pressure. This allows us to use HYCOM directly as a boundary condition since Hydro-CoSMoS already captures the inverse barometric effect itself.

Discharge Boundaries
Fluvial discharges from 22 rivers and streams are included in the model setup. The model obtains forecast data from NWM [43,44]. The NWM operates in several configurations that provide analysis and assimilation of current conditions, a short range (18 h) forecast, a medium range (10 days) forecast and a long term forecast (30 days). Often investigations of fluvial flooding in the coastal zone use daily discharge values as that is what has been primarily available [45,46]. We are currently utilizing the 18-h forecast and when model latency requires it, the last hour is repeated to ensure a discharge time-series of 16 h with 1-h time steps. The Hydro-CoSMoS and NWM are 1-way coupled at upstream locations shown in Figure 1. These coupling points are mostly located upstream in areas with minimal tidal influence to minimize model instabilities at the boundaries and avoid the complications of 2-way coupling between two systems. For the minor tributaries the channel elevation was used to estimate the location of the head of tide, similar to the method used by San Francisco Estuary Institute [47].

Atmospheric Inputs
Wind and mean sea level pressure fields are obtained from NOAA's HRRR forecast system [21,22]. HRRR provides 18 h of forecast data with spatial and temporal resolutions of 3 km and 1 h. Figure 2 shows a flow chart of the Hydro-CoSMoS integrated structure. Hydro-CoSMoS is currently running 16-h forecast simulations every hour, with an average computation time of 24 min on a single processor. Table 1 shows a summary of the models used in this forecast system. Further information about the integrated system can be found at the AQPI website (https://psl.noaa.gov/aqpi/).

February 2019 Storm
In February 2019, California experienced a winter storm with recorded rainfall in excess of 200 mm in many parts of the State. The peak of the storm coincided with Valentines' Day on 14 February, which is why it has been referred to as the "Love Storm". A combination of flooding and landslides closed sections of several highways, such as Highway 1 on the Central Coast, Interstate 5 north of Sacramento, Highway 37, and California State Route 121 in Sonoma County [32,48]. Many rivers reached moderate to major flood stage, and flash flood watches were issued for most of the Bay Area. The storm was caused by a phenomenon known as "Pineapple Express", which is an atmospheric river system that draws warm and wet air from the tropics and brings heavy rainfall to the coast of North America [49,50]. A recent analysis of this atmospheric river event found that an event of this size has a return period of about four years [51,52]. The storms also brought a low atmospheric pressure system and strong winds to the San Francisco Bay, resulting in storm surge and, consequently, elevated coastal water levels. Considering the maximum recorded water level at the San Francisco tide station ( Figure 1) as an indicator of the coastal component of the "Love Storm", this event had a return period ranging from two to three years. The fact that this storm had both coastal and fluvial drivers of flooding made it a perfect test case to investigate Hydro-CoSMoS's ability to capture compound flooding within the San Francisco Bay Area.

Model Setup for Storm Experiments
We have modeled the 2019 February Storm using four different forcing combinations to study the effect of fluvial discharges and remote NTRs on flooding in San Francisco Bay, considering two different types of forcing for offshore and discharge boundaries. HYCOM sea surface height (Figure 3) forcing on the offshore boundaries was turned on and off to investigate the role of remote NTRs. Discharge boundaries were forced with measured real-time discharges to replicate circumstances with accurate high-temporal discharge forecasts available. In contrast, daily-averaged discharges were used to evaluate Hydro-CoSMoS performance under situations with no available discharge forecasts. It is common practice to use daily-averaged discharges when there is no available forecast data. Daily-averaged discharges are the mean of all measured values for each day of the year, creating a daily climatology. For example, the 16 February daily-averaged discharge for the Sacramento River is about 1034 m 3 /s (Figure 3), which is the averaged value over all years of measurement at that measurement site on 16 February [53]. The wind and atmospheric pressure for all four test cases were obtained from HRRR archived data. For this work, we used a two-day spin-up time to start the model; however, in the operational setting, the model starts using results from its previous forecast results (hotstart).  Figure 4 shows the comparison of the measured and calculated water levels at six NOAA tide gauges ( Figure 1) within SF Bay during the February 2019 storm. Without considering remote NTRs, Hydro-CoSMoS underestimated water levels by an average of 15 cm across the entire domain during the storm. It must be noted that water levels mentioned in this work are not what is known as "total water level" because they do not include the effects of wave runup. Figure 5 depicts the comparison between measured and calculated NTRs at six NOAA tide stations. The NTRs were calculated using a 33-h low-pass filter to remove tide from the water level signal. The measured Inverse Barometric Effect (IBE) is also shown in Figure 5. The IBE effect was calculated using measured atmospheric pressure at Oakland (NOAA station 9414763) and by considering a mean sea level pressure of 1020 mbar. The comparison between calculated NTRs and measured IBE shows that without forcing the remote NTRs on offshore boundaries, the coastal water levels are mainly influenced just by local IBE, and the model underestimates the peak of the storm by more than 30 cm. Moreover, Figure 5 shows that water levels are sensitive to the large fluvial discharges from the San Francisco Delta. At Martinez and Port Chicago, which are the closest measurement stations to the San Francisco Bay Delta, using real-time discharges improved the water level estimates by an average of 5 cm compared to the simulations forced by the daily-averaged discharges. Water levels in the rest of the San Francisco Bay are less sensitive to fluvial inflow.

Discussion
We used four different forcing combinations to study the effect of fluvial discharges and remote NTRs on water levels in San Francisco Bay during the February 2019 storm. Figure 6 shows the Mean Averaged Error (MAE) for water levels, NTRs, and 95th percentile water levels. The model performs well at predicting measured water levels and NTRs throughout San Francisco Bay by using HYCOM-based remote NTRs and real-time discharges. Forcing remote NTRs along the offshore boundaries of Hydro-CoSMoS improved the model accuracy across the entire domain, correcting the water level bias by 15 cm (Figure 6) while reducing the error of highest NTRs by more than 30 cm ( Figure 5). Figure 6 shows that the water levels, NTR, and water level 95th percentile errors at every station were similar for each of the four forcing types used to simulate the February 2019 storm. In other words, the errors in estimating water levels were due to inaccurate NTR calculations. Therefore, it is safe to conclude that in order to forecast high water levels and their consequent flooding in San Francisco Bay, NTRs must be captured correctly, and for a storm similar to that of February 2019, that will not be the case if the remote NTRs are neglected. Although this is probably not the case for all of the storms that occur in San Francisco Bay, our results show that using the forecast of remote NTRs can be essential to predict water levels and flooding accurately, and that is why Hydro-CoSMoS imports remote NTR data from HYCOM in its operational setting. Moreover, the results shown in Figure 6 demonstrate that during the February 2019 storm, coastal water levels were more sensitive to fluvial discharges further away from the ocean inlet (Golden Gate) and closer to the delta. The water levels at Martinez and Port Chicago were underestimated by an average of 5 cm using the daily-averaged discharges. One main reason for this pattern is probably the fact that the inflow from the delta contributes to more than 95 percent of freshwater input into San Francisco Bay [27], which explains why the calculated water levels at Martinez and Port Chicago were sensitive to fluvial inputs, while the differences between the simulations forced by real-time and daily-averaged discharges were negligible in the rest of the SF Bay further from the delta. The impact of discharges should extend farther into the bay for a larger event, as it has in the past. In one of the worst storms recorded in the San Francisco bay area in January of 1862, water levels at the mouth of the bay were reported to be 17 cm higher than normal and tidal range was significantly reduced by large fluvial discharges [54]. A recent analysis of California's vulnerability to storms caused by atmospheric river incidents [55], similar to February 2019 storm, found that the 1862 storm was roughly a 500-year return period event that would greatly exceed the state's current flood control infrastructure. Although a rare occurrence, being able to provide some advance notice of water levels and flood extents throughout the bay for this type of storm will be useful for emergency response. Until now, no regional model of San Francisco Bay would be able to incorporate forecasts of the river runoff and its impacts on the entire bay system. As this is a level of storm that is likely to repeat itself [55] this new capability will be useful.
Although coastal water levels calculated by using real-time and daily-averaged discharges were similar for areas far from the San Francisco Bay Delta, their flooding extent was different in areas that are vulnerable to compound flooding, including estuarine zones. An example of this is the flooding extent on the west bank of the Petaluma River. By disregarding remote NTRs and using daily-averaged discharges, the flooded area on the west bank of Petaluma River (Figure 7) was underestimated by 15.5 km 2 compared to the simulation that used real-time discharges and remote NTRs from HYCOM. Compared to the simulation with the best water level results, forced by both real-time discharges and remote NTRs, the flooded area on the west bank of the Petaluma River was underestimated by 2.4 km 2 when using real-time discharges without remote NTRs, and it was underestimated by 4.8 km 2 when using daily-averaged discharges with remote NTRs. On the other hand, the differences on the east bank of the Petaluma River between different models are insignificant. This is because the west bank has a vast low-lying flat area, while on the east bank the area vulnerable to flooding is limited by steep slopes near the river. Although there is no available data to validate the flood extents shown in Figure 7, it is clear that the flooding extent in this estuarine zone is more sensitive to discharges than remote NTRs. Capturing these extents is vital for assessing hazards due to high water levels on roads, levees, and other nearshore infrastructure. Modeled with real-time discharges and remote NTR (blue), modeled with daily-averaged discharges and remote NTR (green), modeled with real-time discharges without remote NTR (red), modeled with daily-averaged discharges without remote NTR (orange).
The spatial pattern observed during the storm of Feb 2019, of fluvial discharge being important closer to the river and unimportant closer to the ocean, both in terms of water levels and inundation extent, is likely repeated for each of the hundreds of tributaries in the Bay, not just the Petaluma River example shown. Fluvial discharges are small relative to the tidal prism, where freshwater inputs contribute to less than 1% of spring tidal prism, reaching 19% during record high discharges [27]. This means the fluvial effects are not observed at the large open water tide gauges, and the small size of the majority of streams means that many are not monitored. However, water levels in the local tributaries are compounded by the confluence of upstream fluvial boundary conditions and downstream SF Bay water levels. Using accurate remote NTRs and fluvial discharge peaks and durations are essential in computing not only water levels but also the flooding extent. The fact that flood extent could be sensitive to discharge values during an extreme event, as well as better performance of the model near the delta in capturing water levels when using real-time discharges prove that having the Hydro-CoSMoS coupled to NWM is vital to forecast compound flood extent correctly in the San Francisco Bay Area.
Much of the work examining compound flooding events in the San Francisco Bay Area has been focused on the larger discharges like the Sacramento River [27,56] including our previous work that was focused on the Napa River estuary [25,28]. For these larger rivers, it has also been demonstrated that there is a correlation between high discharges and the occurrence of coastally driven high NTRs [56]. Unfortunately, the length of record and scarcity of data for all of the smaller streams (more than 450) makes it difficult to assess whether the same relationships hold up; however, it is clear from this modeling and the previous efforts [56][57][58] that the coastal water levels propagate throughout most of the bay, so inherently there will be interaction with the smaller streams. Many of the cities along the San Francisco Bay coast (e.g., San Jose, Hayward, San Rafael, Novato, Napa, and Petaluma) are centered on rivers near the head of the tide, which is indicative that water levels will be coastally influenced, especially in the future with higher sea levels [47]. In addition to the many towns, there is also a significant amount of infrastructure that crosses these small streams, such as fuel lines, sewer laterals, railroads, and highways. Anecdotally (via twitter), the storm we have modeled here produced flooding on several major roadways (Highway 121, Highway 37 and Highway 101). In order to forecast these road flooding events for future storms similar to February 2019 "Love Storm", it will be necessary to incorporate both river discharges and the coastal NTRs. A model, such as Hydro-CoSMoS that leverages the continental (NWM) and ocean basin scale (HYCOM) models to predict local compound flood hazards will be significantly useful. Providing early notice of flooding impacts before an event will allow resources to be allocated to protect the most at risk infrastructure and allow emergency-responders to divert people from hazardous areas. The developers of the forecast system are currently engaging in an iterative process with our stakeholders to establish appropriate alerting criteria that depends upon the role of the stakeholder. These notices will be based on forecasts of arrival time and the duration of the flood and inundation depth map, as well as other location-based hazard criteria such as levee over-topping potential. The settlement patterns of the San Francisco Bay region have created large population centers in places that are impacted both by coastal and fluvial water levels making this a location that can maximally benefit from a forecast model that accounts for both fluvial and oceanic drivers of flooding.

Summary and Future Developments
In this work, we introduced a flood forecast model (Hydro-CoSMoS) for the San Francisco Bay Area. We studied the 2019 February storm to show the effect of fluvial discharges and remote NTRs on compound flood forecasting in San Francisco Bay. Our results show that coastal water levels could be significantly altered by remote NTRs during a storm, highlighting the importance of accurate global and regional ocean forecast models to forecast flooding in the Bay Area. We also showed that having higher temporal resolution discharge data is essential to accurately compute high water levels near the San Francisco Bay Delta. Moreover, it was shown that fluvial discharges affect the flood extent during a compound flooding event in the estuarine zone of a smaller river.
In the future, we will add more fluvial tributaries to the model from the NWM. Given the large number of tributaries in the region, maintaining water level and discharge gauges at all of them is cost prohibitive. However, with the addition of all the ungauged tributaries to our model we can begin to explore the geographic variability in the compound flooding risk around the bay. This modeling framework, with 1D discharge channels into a 2-D coastal system allows for a lot of flexibility in investigating the risks of compound flooding, even in the smallest tributaries.
Moreover, we are in the process of studying the West Coast Operation Forecast System (WCOFS) [13] as an alternative to HYCOM to capture remote NTRs on offshore boundaries of Hydro-CoSMoS. WCOFS has a higher resolution than HYCOM, and our initial assessments show that it could perform better under our forecast system. We will also include waves in our forecast system to provide forecasts of wave runup and overtopping potential along the San Francisco Bay coastline. Moreover, we are conducting data assimilation efforts to use real-time measurements to improve our model results. We believe data assimilation will be a crucial step in generating a robust and reliable compound flood forecast system. Our ultimate goal is to create a seamless vulnerability forecast map based on the three processes that can drive compound flooding: pluvial, fluvial, and marine components. This operational system is planned to be online in early 2022 to provide emergency planners and decision-makers with critical forecast information, such as time of the first flood, water depth, flood duration, and many other products.

Acknowledgments:
We would like to acknowledge Rohin Saleh from Alameda County Flood Control and Water Conservation District for his critical eye while developing this new regional Delft3D model. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: