A cross-scale study for compound flooding processes during Hurricane Florence

We study the compound flooding processes as occurred in Hurricane Florence (2018) that was accompanied by heavy precipitation, using a 3D creek-to-ocean hydrodynamic model. We examine the important role played by barrier islands in the observed compound surges in the coastal watershed. Locally very high resolution is used in some watershed areas in 10 order to resolve small features that turn out to be critical for capturing the observed High Water Marks locally. The wave effects are found to be significant near barrier islands and have contributed to some observed over-toppings and breaches. Sensitivity results from applying each of the three major forcing factors (oceanic, fluvial and pluvial) separately are succinctly summarized in a “dominance map” that highlights significant compound effects in most of the affected coastal watersheds, estuaries and back bays behind the barrier islands. Operational forecasts based on the current model are being set 15 up at NOAA to help coastal resource and emergency managers with disaster planning and mitigation effort.


Introduction
Recently, more frequent occurrences of "wet" hurricanes (i.e. hurricanes accompanied by heavy precipitation) that stall near the coast (Pfahl et al., 2017;Hall and Kossin, 2019) have brought new challenges to coastal communities in the form of 20 compound flooding, which is defined as concurrence of flooding from same or different origins (river, storm surge and rainfall), especially in the coastal transitional zone that sits at the border between coastal, estuarine, and hydrologic regimes (Santiago-Collazo et al., 2019). Compound flooding highlights one of the major pitfalls of the current hurricane intensity scale, which is entirely based on wind speed, leaving the potential rainfall and flooding impacts to be glossed over in initial forecasts that emphasize hurricane category. The record-setting 2020 Atlantic hurricane season including a number of very 25 wet storms highlights the urgency, and exposes the current knowledge gap, for understanding compound flooding processes.
A recent example for compound flood events is Hurricane Florence that impacted a large area of North Carolina (NC) in September 2018. Hurricane Florence was the first major hurricane of the 2018 Atlantic hurricane season. Originated from a strong tropical wave near Cape Verde, west Africa, it acquired tropical storm strength on September 1, followed by a rapid intensification to a Category 4 status on September 4, with estimated maximum sustained winds of 130 mph, and eventually 30 reached its maximum strength on September 11. It made landfall south of Wrightsville Beach, near the border between NC and South Carolina (SC) as a Category 1 hurricane on September 14. The slow motion of the storm after the landfall brought heavy rainfall throughout the Carolinas for several days. Compounded by the storm surge, the rainfall caused widespread flooding along a large swath of the NC coast, and inland flooding in cities such as Fayetteville, Smithfield, Lumberton, Durham, and Chapel Hill. According to a USGS report (Fester et al., 2018), a new record rainfall total of 35.93 inches was 35 set during the hurricane in Elizabethtown, NC. Many other locations throughout NC and SC also set new records for rainfall ( Fig. 1). Florence is a quintessential example of major flooding that a slow moving, moisture-laden storm can bring, even if it does not have strong hurricane wind.
In this paper, we will study the response to the storm in the watershed rivers and estuaries and examine the processes and sources that lead to the compound flooding there. We will also examine the coastal responses to the event and the close 40 connection between watershed and coastal ocean. The existing modeling efforts on compound flooding (Chen et al., 2010;Cho et al., 2012;Dresback et al., 2013;Chen and Liu, 2014;Ikeuchi et al., 2017;Kumbier et al., 2018;Pasquier et al., 2018;Wing et al., 2019;Muñoz et al., 2020) often focus on a subset of the processes (storm surges, tides, waves, fluvial flooding, pluvial flooding, and potential baroclinic effects), leaving gaps in accurately representing the complex interactions among them (Santiago-Collazo et al., 2019). What distinguishes this study from traditional compound flooding simulations is a 45 holistic approach that solves interrelated processes in different regimes and on multiple temporal/spatial scales with the same hydrodynamic core (i.e., the same set of governing equations) in a single modeling framework. An overview of the processes studied in this paper is shown in Fig. 2. The primary tool used in this study is a proven cross-scale 3D baroclinic model designed for effective and holistic simulation of intertwined processes as found during this event. The trade-off between 2D and 3D setups was carefully weighed with the goal of operationalization in mind before the latter was chosen. In short, the 50 advantage of 2D is the speed (about 80 times faster than its 3D baroclinic counterpart) and the simplicity of the set up; the disadvantage is that it misses baroclinic effects and 3D processes that can vary in space and time. For example, the baroclinic effects during the adjustment phase after Hurricane Irene (2011) are discussed in detail in Ye et al. (2020), using a similar model setup as the one used here. Even though different setups (2D, 3D barotropic, and 3D baroclinic) were tuned to their best possible skills, the 3D baroclinic setup was shown to better capture the total elevation during the post-storm adjustment 55 phase. In addition, during the ongoing effort to operationalize the model, we found that including 3D processes greatly simplified the bottom friction parameterization at some coastal locations (e.g., NOAA Station 8447930 at Woods Hole, MA; Huang et al., submitted). A 3D model can also produce relevant 3D variables (e.g., 3D velocity and tracer concentration) that are important for safe navigation and ecosystem health. The 3D model presented in this paper is efficient enough for operational forecasts (see Section 3.1), which are being set up at NOAA (National Oceanic and Atmospheric 60 Administration).

3
The paper is structured as follows. Section 2 will review the study site and available observations collected during the event in the watershed, estuaries, and coastal ocean. Section 3 describes the numerical model used and its set-up. Section 4 presents model validation and important sensitivity test results; the validation is done in a cross-scale fashion from smallscale watershed areas to large-scale coastal ocean. Section 5 discusses the compound effects as revealed by the 3D model in 65 all regimes. Section 6 summarizes the major findings and planned follow-up work.

2 Study site and observation
The focus (high-resolution) area of this study is on the NC and SC coast and coastal watersheds that saw most of the impact from Florence ( Fig. 3b and Fig. 3f). Similar to what we did for other hurricane events, the spatial domain's landward boundary is set at 10 m above the NAVD88 datum, which is deemed sufficient to capture most backwater effects (Zhang et 80 al., 2020). A rich set of observations for physical and biological variables are available from satellites, autonomous instruments (e.g., gliders and Argo floats), in situ stations operated by NOAA and USGS's field estimates collected during after-event surveys (e.g., High Water Marks or "HWMs"; Fig. 2b and Fig. 2c). Analysis and quality control of these datasets have been done by the data distributors, together with uncertainty assessments. Some of the datasets will be presented in the context of model validation sections below to allow for a comprehensive and objective assessment of the model errors and 85 uncertainties. Assessment of compound flood models such as ours inevitably involves observation collected at disparate spatial and temporal scales of several orders of magnitudes contrasts, as illustrated in Fig. 2. To the best of our knowledge, this type of model assessment has rarely been attempted before even in 2D setting, due to the formidable challenges to numerical models (Santiago-Collazo et al., 2019), but is badly needed in order to gain a holistic understanding of the complex processes at play (Ye et al., 2020;Zhang et al., 2020;Huang et al., 2021). 90

Model setups
To capture the storm surge effects we use a large study domain that encompasses the north Atlantic west of 60 o W (Fig. 3a).
An added benefit of using such a large domain in conjunction of a 3D baroclinic model is that the interaction between largeand small-scale processes can be organically examined in a single model. The disruption and oscillation of Gulf Stream by 110 storms can directly affect the coastal inundation (i.e., the fair-weather flooding reported by Ezer (2018)); our results suggest that the converse is also true, as watershed processes can also affect Gulf Stream and other coastal processes (Ye et al., 2020). Therefore, a seamless creek-to-ocean model is advantageous for compound flood studies. SCHISM (schism.wiki) uses efficient semi-implicit solvers to solve the hydrostatic form of the Reynolds-averaged Navier-Stokes equations and transport equation (Zhang et al., 2016), which govern all flow movements (including overland flow in the watersheds as well as estuarine and ocean circulations) inside the 3D model domain. Major characteristics of the model that ensure a balance of accuracy, efficiency, robustness and flexibility include: hybrid finite-element/finite-volume methods and a highly flexible 3D gridding system ("polymorphism") that combines a hybrid triangular-quadrangular unstructured 120 grid in the horizontal dimension; and localized sigma coordinates with shaved cells (dubbed as LSC 2 ; Zhang et al., 2015) in the vertical dimension. The polymorphism allows a single SCHISM grid to seamlessly morph between full 3D, 2DH (2D depth-averaged), 2DV (2D laterally averaged), and quasi-1D configurations. The employment of shaved cells near the bottom in particular faithfully preserves the original bathymetry without any smoothing required. The detrimental effects of bathymetry smoothing on important physical and biological processes (e.g., residual transport, lateral circulation, nutrient 125 cycling etc) have been documented in Ye et al. (2018) and Cai et al. (2020). Similar to a recent compound flooding study using SCHISM for Hurricane Harvey (Huang et al., 2021), the current model domain covers the entire US East Coast and the entire Gulf of Mexico, with all major bays/estuaries and coastal watersheds resolved (Fig. 3). The horizontal grid, generated using the software SMS (aquaveo.com), has 2.2 million nodes and 4.4 million elements (Fig. 3). About 50% and 40% of the elements have resolution finer than 300 m and 220 m respectively (Fig.  130 3e). The grid bathymetry is interpolated from a combination of DEM sources from coarse (ETOPO1, 90-m Coastal Relief

Model
(https://ngdc.noaa.gov/mgg/coastal/crm.html)) to fine resolution (1/9 arc-second CUDEM (https://www.ncei.noaa.gov/metadata/geoportal/rest/metadata/item/gov.noaa.ngdc.mgg.dem:999919/html), 1-3 m CoNED (https://www.usgs.gov/core-science-systems/eros/coned)), and the vertical datum used is NAVD88, with appropriate conversion between datums done by VDatum tool (vdatum.noaa.gov). Note that NAVD88 is a more convenient datum to use 135 in the model as most of recent observation data refer to this datum, and therefore we use this datum in the model setup and allow the model to automatically set up the sub-tidal surface slope from coastal ocean into watershed (due to the friction effects). Altogether close to 400 DEM tiles are used to cover such a large region. The horizontal grid resolution ranges from depth of coastal streams (e.g., in the South Carolina watersheds), which is a potential error source of our model. In the vertical dimension we use 1-43 layers, with 43 layers being applied in the deep ocean and 1 layer in most of the watershed ( Fig. 4), thus effectively rendering the model 2DH there, which is sufficient for processes like overland flow and inundation. Table 1 shows the setups for "baseline" and important sensitivity simulations used in this paper. For the baseline, imposed at the ocean boundary (60 o W) are tidal elevation and barotropic velocity of 8 tidal constituents (S2, M2, N2, K2, K1, P1, O1, 155 and Q1), extracted from the FES2014 database (https://datastore.cls.fr/catalogues/fes2014-tide-model/). The baroclinic components for the elevation and velocity are derived from the daily outputs of the "HYbrid Coordinate Ocean Model" (HYCOM; https://www.hycom.org/). The initial condition for the water elevation is set to be 0 in all areas with positive grid depths ("wet") and to be equal to the bottom elevation at "dry" watershed locations (0 initial water depth). The initial condition for the horizonal velocity is zero in the watershed and is interpolated from HYCOM elsewhere. Commensurate with the non-zero velocity (fully dynamic state) are salinity and temperature values interpolated from HYCOM; however, in the nearshore areas where HYCOM results are less accurate, the initial conditions for salinity and temperature are interpolated from the sparse observation at several USGS gauges in order to speed up the dynamic adjustment process. The surface meteorological forcing applied in the model is a combination of two products: (1) a high-resolution ERA re-analysis product from the European Centre for Medium-range Weather Forecast (ECMWF) with a ~9 km horizontal resolution (see 165 Acknowledgements), and (2) a 3-h time interval with NOAA's High-Resolution Rapid Refresh (HRRR; https://rapidrefresh.noaa.gov/hrrr/), which is a cloud-resolving and convection-allowing atmospheric model with a 3-km horizontal resolution and a 1-h time interval. The friction of the baseline model was tuned in the wet area (river, estuary, ocean; lower than 1 m, NAVD88) and on higher grounds (higher than 3 m, NAVD88) separately. In the wet area, drag coefficients within a range of 0.001-0.01 were tested. Commonly accepted default value of 0.0025 gave good error statistics 170 near the landfall site; but the sensitivity is small in the range of 0.001-0.005. In the watershed, drag coefficients within a range of 0.01-0.5 were tested. The optimal value was chosen based on the High Water Marks (HWMs) comparisons (Section 4.3) at 276 locations collected by USGS. A small friction value within this range tended to under-predict the elevation at HWMs, and a large value led to over-prediction. Values within the range of 0.02-0.05 gave good error statistics. We chose 0.025 because it gave slightly better results in the Cape Fear River watershed near the landfall site. To sum up, bottom 175 friction coefficient is set at a constant value of 0.0025 at all "wet" locations, then linearly increases to 0.025 as the ground elevation increases from 1 m to 3 m (NAVD88), then a constant value of 0.025 is used for higher grounds where the bed texture is generally rougher than the riverbed. Note that this is the parameterization based on the region influenced by Hurricane Florence. Spatially varying parameterization of bottom friction for different systems is an on-going effort as we study more recent hurricanes and operationalize the model along the East Coast and Gulf Coast. However, as we presented 180 in Ye et al. (2020), Zhang et al. (2020) and Huang et al. (2021), the choices described above seem to work fine in general for other systems as well. The method used to impose the river flow in the model is described in the next subsection.
Choices of the "baseline" model parameters are similar to those used for Irene (Ye et al., 2020). The time step is 150 seconds (sensitivity tests using 100-150 s gave very similar results). The level-2.5 equation turbulence closure scheme chosen is from the generic length scale model of k-kl (Umlauf and Burchard, 2003). The simulation starts from 2018-08-24 185 00:00 (UTC) and lasts for 36 days to cover the hurricane and ensuing restoration period. Although the model covers a large domain, most of the elements (those in the watersheds) are quasi-2D, making it efficient enough for operational forecast. For the baseline run, the real time to simulation time ratio is 80 with 1440 cores on TACC's Stampede2 cluster and 30 with 480 cores on W&M's SciClone cluster. Intel Skylake cores with a nominal clock speed of 2.1 GHz were used on both clusters. This means a 3-day (typical forecast duration) simulation will take about 0.9 hours using 1440 cores or 2.4 hours using 480 190 cores.

Scenario Description Baseline
With forcing from tides, atmosphere, rivers (from NWM) and precipitation, initialized with HYCOM GOFS 3.1 Baseline_Wave Baseline with added wave effects Ocean Baseline forced by ocean and atmosphere only (i.e., tides and storm surge) River Baseline forced by rivers (freshwater inputs from NWM) only Rain Baseline forced by precipitation (directly on top of the domain) only

Coupling with NWM 195
River discharges are introduced into our model at its land boundary. About 6752 intersection points are identified between the NWM river segments and SCHISM's land boundary, where the freshwater is injected as volume sources. Inside the model domain, streamflow, overland flow and precipitation are directly handled by the hydrodynamic core of SCHISM. This fully coupled configuration is rare in the existing compound flooding simulations (Santiago-Collazo et al., 2019). To ensure the accuracy and robustness of SCHISM in simulating hydrological and hydraulic processes including the overland flow, we 200 already examined the model's performance in both lab-scale and field-scale tests in a previous study (Section 2.2 and 2.3 of Zhang et al. (2020)) and applied the model in the Delaware Bay watershed including the Delaware River with a hydraulic jump up to 40 meters above the NAVD88 datum ( Fig. 14 in Zhang et al. (2020)). The NWM segments explicitly reproduced in our grid (Fig. 3c) during the mesh generation stage help capture the bathymetry of main flood pathways (thalwegs).
However, flow is not restricted to these 1D segments; in fact, precipitation may generate overland flow on any 2D horizontal 205 grid elements in the watershed. Note that the river flows injected the land boundary have indirectly incorporated the precipitation that occurred outside (but not inside) the model domain, and therefore, the addition of direct precipitation onto the model domain is appropriate and is an integral component of the compound flood processes. To accurately simulate the initial movement of the very thin layer of rain water on the dry land, which is dominated by friction, a very small threshold of 10 -6 m is used to differentiate between wet and dry states (Zhang et al., 2020). Since we have no information on the scalar 210 concentrations for river inflows and rainfall, we applied 0 PSU for salinity and ambient water temperature (i.e., the temperature at the local receiving cell calculated without accounting for the rivers or raindrops) for the injected river water and also for the rainfall. Obviously, the latter represents a source of uncertainty for the modelled temperature results. As explained in Huang et al. (2021), heat exchange between air and water would misbehave on such a thin layer of water on watershed, so a threshold of 0.1 mm is set for local water depth, below which the heat exchange is turned off. As a model 215 limitation, infiltration is neglected in this work. In the case of Hurricane Florence induced flooding, we expect the effect of infiltration to be minor. According to NOAA's weather map (https://www.wpc.ncep.noaa.gov/dailywxmap/), there was continuous rainfall along the US east coast from Sep 11, 2018 to the date of Florence's landfall (Sep 14, 2018), so the infiltration capacity of the soil was already reduced. Moreover, "wet" storms like Hurricane Florence (2018) and Hurricane Harvey (2017) tend to dump large amount of rain fall at a location for days because of the slow movement of the storm, so 220 most of the rainfall is on saturated soil. As another model limitation, the drainage in urban settings is not included in our model. This may have led to some occasional big errors in the predicted elevation on high water marks, for example the one large error in the urban area in Fig. 13d. We do have a plan of explicitly accounting for infiltration and drainage as volume sinks based on NWM (or other hydrologic models). However, considering the additional uncertainty this would bring, for now we choose to continue improving more important aspects of the model for operational use; the focus is on the quality of 225 model grid, which is very likely responsible for most of the existing large errors.
The assessment of NWM calculated flow against observed flow at the two largest rivers in the region is shown in Fig. 6.
Similar to our findings for other storm events, the flow produced by this particular version of NWM (v2.0) is generally consistent with USGS observation but tends to show narrower and higher peaks, with roughly the same total volume of water during the event (Fig. 6). The observation indicates the peak streamflow occurs about 7 days after the landfall, which is the 230 time it takes for the rainfall induced flood to reach the coastal rivers. Note that there is typically a time lag of 1-2 days between the peak flow in NWM and the gauged flow (Fig. 6). The forcing errors in the magnitude and timing of NWM's peak flow should explain part of the model errors especially in the watershed. For example, we found that replacing the NWM streamflow with the gauged flow at USGS Station 02109500 (Waccamaw River at Freeland, NC) improves the model skill locally. However, this is not cost-effective for our goal of operationalizing this compound flood model along the US 235 East Coast and Gulf Coast. The developers of NWM (Gochis et al., 2018) showed that NWM's model skill was improved by each version update, with 44% of the gauges having bias < 20% in the latest version (NWM v2.0). We will adopt the newest and best NWM version in our ongoing study and operational forecast as soon as it is available. And we are open to using any other hydrologic sources to drive our model.

Model validation and sensitivity
In this section we will assess the model results for elevation, inundation and flow in the watershed and estuary. The spatial scales covered by the validation vary from O(10 km) to O(1 m). The model validation for non-storm period will not be 250 discussed here; in short, the averaged amplitude error for the major constituent (M2 in east coast and K1 in northern Gulf of Mexico) for non-storm period is 3-4 cm (see Huang et al., 2021). We will start by looking at the wave effects nearshore.

Wave effects
Multiple breaches and over-toppings were reported across several NC barrier islands during the event, including Surf City, North Top Sail Beach, and New River inlet (https://www.wusa9.com/article/weather/before-and-after-hurricane-florence-255 changes-north-carolina-coastline/65-595918389). Some of these breaches may be related to significant wave activities, for example the maximum wave height at buoy 41025, ~30 km offshore from Cape Hatteras, reached ~10 m. Interestingly, the maximum wave heights are more modest nearer to the landfall: 4-6 m at Buoys 41159 and 41013 (see Fig. 3b for their locations). Therefore, to investigate this possibility, we have also conducted a simulation with the wave model inside SCHISM activated ("Baseline_Wave" in Table 1). The details of the wave model (Wind Wave Model) have been described 260 in Roland et al. (2012) and the coupled model has been applied to other systems (Guerin et al., 2018;Khan et al., 2020). The wave model was initialized and forced at the ocean boundary by a global Wave Watch III simulation (ftp://ftp.ifremer.fr/ifremer/ww3/HINDCAST), and used a spectral resolution of 36 directional bins and 24 frequency bins to cover a frequency range of 0.04 to 1 Hz. The coupling time step between the two models (at which information was exchange including surface elevation, velocity and wave radiation stress) was 600 seconds. 265 Our results indicate that the barrier islands near Surf City, North Top Sail Beach and New River inlet were indeed overtopped with 1-2 m of water ( Fig. 7b; the locations of the islands can be seen in Fig. 8b). On the other hand, large portion of the island north of Cape Hatteras (Outer Banks) was spared (Fig. 7c), likely due to its N-S shoreline orientation, even with the 10 m wave approximately 30 km offshore from there. More quantitative validation for the breaching processes is out of scope here because (1) we do not have accurate and update-to-date bathymetry just before the event, and (2) more 270 importantly, a sediment transport study is required to simulate the bathymetric changes.
The wave effects on the surface elevation are further quantified in Fig. 8, which suggests that the effects are most pronounced (with 30 cm or larger differences) inside the estuaries and Albemarle-Pamlico Sound (APS) due to large wave breaking nearby. In the intermediate and deep water, however, the wave effects are on the order of a few centimetres and negligible (Fig. 8a). The Baseline results (without wave effects) also showed over-topping of the barrier islands similar to the 275 Baseline_Wave (not shown).  In summary, the wave effects are significant nearshore and have contributed to the observed breaching and over-toppings of barrier islands. However, because the current study does not focus on the breaching processes and because of the 285 significant computational overhead introduced by adding the wave model (>50%), we will proceed in the following by using the run without waves as the "baseline".

Bays and estuaries
We assess the calculated water levels at six NOAA tide gauges near the impact area (see locations in Fig. 3b). Three gauges are facing the open ocean (Charleston, Springmaid, and Wrightsville), and the other three gauges are either sheltered 290 inside a bay (Beaufort) or behind barrier islands (Hatteras and Oregon Inlet). The responses to the hurricane are different at those gauges, as seen from the total elevation (Fig. 9a) and the sub-tidal signals (Fig. 9b). The latter applies a low-pass Butterworth filter (Butterworth, 1930) only preserving longer-period (longer than 2 days) components. The observation shows see set-downs of ~0.3-0.5 m at Charleston, Springmaid, Hatteras and Oregon Inlet followed by surges of ~0.2-0.5 m; and surges of ~0.5 to 1 m at Wrightsville and Beaufort (Fig. 9). The different responses at these gauges are due to the wind 295 curl of Florence that led to different dominant wind directions between southern and northern stations, and are also due to specific geographic settings of each gauge. Most intriguing are the prominent set-downs observed at Hatteras and Oregon Inlet, which are explained by a combination of wind direction and blocking effects of barrier islands. Fig. 10 demonstrates that around the time of landfall of the hurricane, the predominantly westward wind felt in the Pamlico Sound has pushed water away from the barrier island. Meanwhile, the surge that propagated from the ocean side is effectively blocked by the 300 barrier island chain, thus creating a ~70 cm elevation difference between the waters immediately outside and inside the islands (Fig. 10), and the difference lasted for ~2 days (Fig. 9). The mechanism causing the water level set-downs at the two South Carolina stations is similar to that causing the set-downs behind the barrier islands in North Carolina. The two South Carolina stations (Charleston and Springmaid) are located to the south of the landfall site and the wind direction is from the land to the ocean, pushing water away from shore. The model captured the different regional responses; overall, the averaged 305 MAE (mean absolute error) for elevation is 11 cm. The averaged MAE for the subtidal comparison is 8.6 cm and averaged correlation coefficient is 0.92. The peak errors at different stations occur around the storm surge, with a maximum overprediction of 0.64 m for the peak surge at Spring maid Pier, SC. The overpredicted peak surge can lead to overpredictions in elevation on coastal high water marks (HWMs). In addition, there is a maximum under-prediction of 0.66 m for the set down at Hatteras, NC, mainly due to the mismatch in the set down timing. The uncertainties in wind forcing may be the main 310 cause of the error, which is predominantly from subtidal signals. The grid quality near the barrier islands may also contribute to the error. Adding wave effects slightly increases the surge and rebounding waves at the last three gauges, resulting in slightly better model skills there.  Table 1). The stations in (a)-(d) are exposed to the ocean, showing small differences among scenarios; the two stations in (e) and (f) are on the land side of the barrier islands, showing larger differences between the "Baseline" and the "Ocean" scenario in the post-storm period. Note the plots have different y-axis ranges.

Obs
Ocean Baseline Baseline_Wave

Watershed
High Water Marks (HWMs) are collected by USGS experts more than two weeks after Hurricane Florence's landfall on Sep. 14. They are derived from small seeds or floating debris carried by floodwaters that adhere to smooth surfaces or lodge in tree bark to form a distinct line, and also by stain lines on buildings, fences, and other structures. Therefore, HWMs are 330 time sensitive and usually have vertical uncertainties of ± 0.3 feet, or equivalently ± 0.09 m (Koenig, et al., 2016;Austin, et al., 2018).
The simulated elevation on 276 HWMs in the NC and SC watersheds are compared with field estimates (Fig. 11). The model is able to capture the transition from estuarine to riverine regimes; note that the averaged bottom elevation for all observation points is 3.8 m based on NAVD88, and about 70% of the points are located above 2 m (NAVD88), beyond the 335 reach of storm surges. Overall, the averaged MAE for all HWMs is 0.73 m, with a correlation coefficient of 0.92 and a positive mean bias of 0.09 m. There is a slight positive bias on near-shore HWMs (Fig. 11), corresponding to the overprediction of peak elevation at coastal stations (Fig. 9a-c). These skill scores are similar to what we obtained for Hurricane Harvey (Huang et al., 2021).

Fig. 11: Comparison of HWMs. (a) Field estimates from USGS; (b) model prediction; (c) model minus observation; (d) regression between model prediction and observation; (e) histogram showing the error distribution.
Some of the HWM estimates are derived from the measured time series at USGS in situ gauges. Fig. 12 shows the comparison for both elevation (gauge height) and discharge in a large river in the study region. The gauge is in the interior of our grid, near "Freshwater Source 2" in Fig. 6. Because the observation's vertical datum is NAVD29 and the model's datum 345 is NAVD88, we have adjusted the mean model elevations to match the observed mean in the elevation comparison (Fig.   12b); in other words, only the elevation variability is compared. Our model over-predicted the flow and under-predicted the flood-induced surges. Using more accurate fresh water source at the land boundary, improving the channel representation in the model grid, and locally adjusting the bottom friction should help improve the skill.

355
Our tests show that the simulated elevation on the High Water Marks (HWMs) in the watershed is sensitive to: grid resolution, precipitation, river inputs through the land boundary, and bottom friction. Grid resolution and quality is the most important factor. Misrepresenting flood pathways can easily lead to errors of a few meters near some very localized features like ditches, highways etc. Fig. 13 illustrates such an example around the Burnt Mill Creek in the city of Wilmington, NC. 360 Large HWM errors were found in the preliminary setup, because the computational grid did not resolve the small creeks that served as the main conduit in draining out the storm water after the flood. This resulted in stacking of water locally and thus large over-prediction of HWMs there. The channel of the creek is about 6-10 m wide, and once resolved using 2 rows of quadrangles as done in the baseline (Fig. 13c), the model skill was greatly improved (Fig. 13d). The only remaining large error point in the "baseline" occurs in an urban area away from the river (Fig. 13d), likely due to the building or drainage 365 effects that have not been incorporated in the model. The defects in grid quality can lead to large errors that are not likely to be rectified by tuning other parameters. To fix the remaining few large errors away from the landfall site, grid quality should be examined first. The continuous improvement on this model grid is part of an ongoing effort of operationalizing the model along the US East Coast and Gulf Coast, and we will report this in future studies. Resolving small-scale flow routing features in a national scale requires automated tools such as Pysheds (https://github.com/mdbartos/pysheds) that can detect and 370 delineate the channels automatically. Initial tests showed very promising results from this package. We remark that it is feasible to resolve these features efficiently without significantly increasing the grid size due to SCHISM's flexibility and robustness in handling poor-quality meshes. Afterwards, the inclusion of urban drainage should reduce the occasional large errors there. Other factors such as uncertainties in DEM, precipitation and the river flow through land boundary also play minor roles. 375

Compound effects
A carefully validated 3D model such as the one presented here can effectively separate out compounding factors from different sources: coastal surge, river flooding and precipitation. In this section we apply this approach to examine the 385 contributing factors to the observed total flooding during Florence. The design of the numerical experiments is such that we selectively turn on/off forcing from ocean, river, and precipitation to examine their individual effects in isolation (Table 1).
As an overview, the conditions of maximum inundation extent from all scenarios are listed in Table 2. To facilitate the comparison of inundated area, a practical value (1 foot) on the same order of the mean inundation depth is used as a threshold. For the two indices shown in Table 2, the baseline values are always larger than those from a single sensitivity 390 test. This confirms the existence of compound regions in the two states (North Carolina and South Carolina) during the hurricane. More details of each forcing's effect and the compound effects are discussed below. Turning off both rivers and precipitation (i.e., storm surge only) is expected to have major impact on flooding in the watershed. This is confirmed in Fig. 12 in the previous section. Not surprisingly, without rivers and precipitation, watershed is mostly dry as the storm surge cannot propagate over steep terrains. As a result, the predicted HWMs are biased too low ( Fig. 14) as the steep topography quickly damped out any surges brought in by the ocean. This leads to systematically underpredictions in the watershed and a 64% increase in MAE compared to the baseline (Fig. 11). On the other hand, Zhang 400 et al. (2020) demonstrated that the storm surges can propagate much further into watershed if watershed rivers are included.
There is no apparent deterioration of model skill on the near-shore HWMs, because those locations are predominately affected by oceanic processes.

405
Fig. 14: Simulated HWMs from the sensitivity run "Ocean", i.e., without the freshwater inputs from NWM or precipitation. Note the underpredictions in the watershed and worse model skill compared with "baseline" (Fig. 11).
Less obvious are the effects of rivers and precipitation on the observed surges in the coastal bays. Fig. 9 indicates that the effects are negligible at the 3 coastal stations away from barrier islands (Fig. 9a-c), as the large amount of freshwater from 410 the watershed directly drains into the coastal ocean (which has much larger volume of water). On the other hand, the impounding effects are clearly seen at the 2 stations behind the barrier islands (Fig. 9ef), where the discharged water is trapped for almost a week; the difference reaches ~30 cm at times there (see Fig. 8). As we showed in Section 4.2, the barrier islands are effective in creating separation and thus large elevation gradient between water immediately outside and inside (see Fig. 10). 415 To assess the contributions from each of the 3 forcing factors to the total sum, we follow Huang et al. (2021) and use the concept of "disturbance". We recognize that for compound flooding processes that involve both ocean and watershed, neither the water surface elevation nor the water depth is a satisfactory metric, because the nominally large water elevations on the high ground of watershed are dominated by the high bottom elevation there, and the large water depths in the bays and ocean are dominated by the local bathymetry. Therefore, we adopt the concept of "disturbance", defined as: 420 where is the water surface elevation and ℎ is the bathymetry (positive downward; ℎ > 0 for ocean and ℎ < 0 for high grounds in watershed), so ( + ℎ) is water depth. Basically, represents the departure from "initial condition" (either the initial water surface or bottom, whichever is higher). Note that is continuous across ℎ = 0. On the initially "dry" ground in watershed, represents the local water depth; whereas at initially "wet" locations, is simply the surface elevation. is 425 also a smoother metric to measure the compound effects as one transition from oceanic into watershed regimes.
The comparison between the maximum disturbances from the three experiments using only one of the three forcing factors and the total sum helps shed lights on the contributions from each. The results are presented in Fig. 15, in terms of proportion of total maximum disturbance as explained by each forcing factor at a given location (therefore, the sum of all proportions equals unity). Our results clearly indicate that: (1) the ocean (and atmospheric) forcing dominates in the open 430 and part of southern Pamlico Sound in the APS system; (2) river forcing is most dominant in the river network of watershed; (3) precipitation effects are dominant in other parts of the watershed away from the river network. However, the presence of the barrier islands significantly complicates the interaction among different forcings (e.g., the significant contribution from precipitation in Pamlico Sound as shown in Fig. 15c). On the other hand, the inclusion of the wave effects is not expected to alter the findings here because their contribution to the total elevation is relatively minor as compared to the atmospheric 435 effects (see Fig. 9).
The competition among different forcing factors in different regions can be succinctly summarized as a "dominance" map as shown in Fig. 16: a factor is deemed dominant if it explains at least 80% of the total disturbance; if, however, none of the 3 factors contribute to 80% or more at a particular place, the nonlinear compound effects are expected to be significant there.
The ocean response is overwhelmingly dominated by the oceanic and atmospheric forcing, the response in the watershed 440 rivers by the river flow, and the response in large portion of highly-elevated watersheds by the precipitation, as seen in Fig.   16. Most of the response in the southern Pamlico Sound is of oceanic origin, because of the wider openings to the south (e.g., Ocracoke inlet). On the other hand, there are only a few narrow inlets to the east (e.g., Oregon Inlet), thus effectively blocking off the oceanic influence there ( Fig. 16b; also see Fig. 10). It is the "grey areas" (Fig. 16) of compound flooding zones that are most intriguing. These include most of the APS and coastal rivers (Fig. 16b), where the weakened oceanic 445 influence competes with river flow and large rainfall there (Fig. 15). In the estuaries with large river discharges with limited openings to the coastal ocean ( Fig. 16c and Fig. 16e), the compound flooding zone is a result of competition among all 3 factors. On the other hand, for the estuaries to the south of the landfall in SC that are less influenced by the heavy precipitation, and not protected by barrier islands and have smaller river discharges, oceanic effects can penetrate very deep into the coastal watershed (Fig.16d). Note that the ocean dominance near Wrightsville Beach, NC (NOAA Station 8658163) 450 and Springmaid Pier, SC (NOAA Station 8661070) may be exaggerated, considering the overestimated peak elevation there (Fig. 9). The compound map like the one shown in Fig. 16 clearly demonstrates the urgent need for a holistic management approach in the planning of mitigation effort for the flood hazard during and after hurricane events.  Table 1); (b) 455 river force ("River" in Table 1); and (c) precipitation ("Rain" in Table 1). The value is the proportion of a factor's individual "disturbance" (see definition in Section 5) to the sum of the disturbance from all factors. The colors from blue to red represent increasing importance of a factor at a specific location.

Conclusions
We have successfully applied a 3D cross-scale model to examine the compound flooding processes that occurred during Hurricane Florence. The model is fully coupled in the sense that the hydrologic and hydrodynamic processes are solved by 465 the same set of governing equations. Limitations of the model include the neglection of infiltration and urban drainage, which will be implemented soon. The model was validated with observation data collected in the watershed and coastal ocean. The mean absolute errors for major variables are: 11 cm for coastal elevation, 72 cm for High Water Marks (HWMs).
Locally very high resolution was used in some watershed areas to resolve small features that were critical for a good model skill for the HWMs. The wave effects were found to be significant near barrier islands and had contributed to some over-470 toppings and breaches there. The validated model was then used to reveal significant nonlinear compound effects in most parts of coastal watersheds and behind the barrier islands. The barrier islands were shown to be particularly effective in separating the processes in the water bodies inside and outside.
The results of the current study, especially the regional compound zone map, filled in a critical knowledge gap in our understanding of compound flooding events. In fact, operational forecasts based on the current model are being set up at 475 NOAA to help coastal resource and emergency managers with disaster planning and mitigation effort. The model can also be used to facilitate new scientific discovery of novel coastal processes; for example, preliminary results for the fate of pollutants discharged from the watershed suggest that the large watershed outflow resulted from heavy precipitation played an essential role in exporting pollutants far into the ocean through the large and long-lasting freshwater plumes that occurred after the event. 480

Data availability
The model source code is freely available at https://github.com/schism-dev.