A New 1D/2D Coupled Modeling Approach for a Riverine ‐ Estuarine System Under Storm Events: Application to Delaware River Basin

Numerical simulations of three of the most severe historical tropical cyclones to affect the Delaware River Basin (DRB) are used to evaluate a new numerical approach that is a candidate model for the inland ‐ coastal compound ﬂ ood forecast. This study includes simulating interactions of tides/surges, freshwater stream ﬂ ows, winds, and atmospheric pressure for the DRB. One ‐ way coupling between the hydrologic (National Water Model [NWM]) and the ocean/wave (ADvanced CIRCulation model/WAVEWATCH III [ADCIRC/WW3]) models for the Delaware river ‐ estuarine system is developed. The links between the coastal processes and the NWM are provided by two different hydraulic and hydrodynamic models: (i) a well ‐ calibrated public ‐ domain 1D hydraulic solver model (Hydrologic Engineering Center's Results show Simulations show work


Plain Language Summary
The East Coast of the United States is prone to powerful winter nor'easters and tropical cyclones. Both weather patterns produce strong winds, large storm tides, enormous ocean waves, and potentially extreme precipitation. This makes the U.S. East Coast vulnerable to catastrophic damage and loss of life. The Delaware River Basin is part of the highly developed and densely populated megalopolis of the northeastern United States. The unique properties and movement of individual storms and the complexity of the coastal morphology require the implementation of the most skillful and comprehensive modeling approach for forecasting and planning purposes. In this paper, the appropriate modeling system is evaluated for the hurricanes that made landfall along the East Coast of the United States to investigate the compound inland-coastal flooding. Numerical results are compared with observational data from National Oceanographic and Atmospheric Administration (NOAA). In agreement with observations, simulations indicate that water level estimations depend on a precise representation of both river inflows and elevated sea levels. The findings of this work show that using the proposed modeling framework has advantages over existing river hydraulic models, particularly near a coastal plain estuary during storm events.

Introduction
The interactive riverine and coastal processes, especially during short-term changes like major storms or under long-term impacts such as climate change and sea level rise, are responsible for large cultural, economic, and biodiversity impacts in such highly populated estuarine systems (Moussa & Bocquillon, 2009;Neumann et al., 2015;Peduzzi et al., 2012;Zhang et al., 2020).
In river-estuary transition areas or low-gradient coastal watersheds, both inland hydrological and coastal ocean processes significantly interact. A river-estuary transition zone is subjected to tides, waves, saline water, freshwater inflow, rainfall, and sediment. River-bay-ocean systems are at risk to flooding exposures from combined influences of river inflows and storm surges produced by tropical cyclones (TCs; i.e., compound flooding) (Santiago-Collazoa et al., 2019;Zscheischler et al., 2018). Intense rainfall during extreme storms increases the riverine inflows and, as a result, produces overland flow and floodplain inundation. North Atlantic TCs are accompanied by strong winds and torrential precipitation capable of generating tremendous storm tide, enormous ocean waves, and large streamflow simultaneously. The culmination of complex interactions of forcing including surges, waves, strong turbulence, freshwater inflows, winds, and atmospheric pressure produces life-threatening flooding; catastrophic damage of infrastructure, homes, and commercial properties; and erosion of the coastlines (Bakhtyar, Orton, et al., 2018;Chang, 2001;Czajkowski et al., 2013;Di Liberto et al., 2011;Villarini et al., 2014).
The natural ecosystem, coastline's morphology, and water quality are often left in disrepair after the storm's impact. Besides, the cost of such events is estimated at about $28 billion per year and is expected to rise to $39 billion by 2075. Under extreme meteorological events, intense streamflow and storm surge occur at the same time or in temporal proximity. Sometimes, sustained precipitation and the combined effect of storm surge and flood tide obstruct the streamflow towards the bay/delta, which produce a disaster flood. Moreover, it is expected that the number of severe meteorological events will increase and their compounding influences will become stronger via sea level rise .
The forecast of flood inundation is vital for emergency action plans, risk assessment, flood insurance rates, and water resources management. Many of the existing flood inundation models disregard precipitation and river discharge (Santiago-Collazoa et al., 2019). As in these models, some important physical processes like momentum exchange are omitted , it is vital to simulate the combination of above-mentioned flooding mechanisms (i.e., interaction of streamflow and combined effect of storm surge and flood tide) and determine the effect of each one and their combined impacts on flood forecasts (Santiago-Collazoa et al., 2019).
Riverine-coastal hydrodynamics may be investigated either via field and laboratory experiments or numerical models. Attaining robust data through the laboratory and field measurements in extreme storm events is difficult and expensive. In addition, progress in numerical methods and computational power increase the attractiveness of numerical simulation of riverine/bay hydraulics and hydrodynamics (Chen et al., 2015;Neal et al., 2012). Furthermore, operational and probabilistic flood forecasting systems, which provide oceanographic conditions in nowcast and forecast time, need many simulations per hour (Bakhtyar, Orton, et al., 2018;Orton et al., 2016). Therefore, the advanced numerical models provide a tool for scientists to simulate the motions of air and water in the riverine and estuary systems.
The National Weather Service (NWS) is trying to efficiently minimize the flood risks and flood-induced losses in life and property. In order to achieve this goal, NWS needs a precise flood modeling framework. An accurate flood modeling system requires state-of-the-art numerical models, which appropriately consider the significant physical processes that affect the flooded areas (i.e., integrated storm surge/tide, water levels in rivers and tributaries). For many years, NWS has used 1D hydraulic models to forecast water levels on main rivers (http://www.nws.noaa.gov/ohd/hrl/hsmb/hydraulics). In coastal plain estuaries where rivers broaden and flow into bays, 1D modeling assumptions can be a limiting factor (Bakhtyar, Maitaria, et al., 2018Gong et al., 2007). The current kinematic wave hydraulic module in the National Water Model (NWM) is not able to resolve wind forcing, ocean tide, or baroclinic forcing, which can play dominant roles in changing total water level in coastal regions, especially during TCs. Additionally, the multifaceted interactions at tidal, subtidal, and shorter time scales create challenging assessment of riverine/bay models (Warner et al., 2005). Because of these multiple influences, a complex modeling scheme is necessary to fully understand and forecast the water behavior in the river-estuary transition zone. The model should provide considerable flexibility in resolving complex geometry/bathymetry/physics for flood forecasting and water management activities in the coastal-estuary zone.
Commonly, three main types of numerical models are used for hydrodynamic simulations of major storm events (Xie et al., 2016): (i) wave models (Bakhtyar, Orton, et al., 2018), (ii) ocean circulation models (van Heerden et al., 2007), and (iii) coupled wave and ocean models (Deb & Ferreira, 2018;Dietrich et al., 2013). The most accurate type of these models is the coupled ocean/wave model that includes the momentum exchange between currents and waves (Donelan et al., 2012). Although some of the open source ocean models have shown potential for inclusion within the NWM approach, the well-established hydrodynamic ADvanced CIRCulation (ADCIRC) model coupled with third-generation wave WAVEWATCH III (WW3) model has been broadly used in recent years to simulate storm surge, tides, and waves. ADCIRC (Luettich et al., 1992) is useful for studying storm impacts and flood forecasting in coastal zones (Kress et al., 2016;Orton et al., 2015). WW3 is also widely used in offshore and nearshore applications in high-resolution numerical domains. New developments in the WW3 model have improved the model accuracy, robustness, and scalability on various high-performance computing (HPC) environments compatible with community-based coupling infrastructures (NOAA Environmental Modeling System [NEMS]).
The main objective of this study is to investigate the combined effects of storm tide and riverine streamflow on flood/inundation modeling at the transition zone using an advanced ocean/estuarine/hydrologic modeling approach. There are some approaches to link estuary, ocean, and hydrologic models to create a river-estuary-ocean modeling system: (i) directly linking an ocean model with a hydrologic model; (ii) using a simplified 1D inland hydraulic model between a hydrologic and ocean model; and (iii) applying a more robust but more complex hydrodynamic model between a hydrologic and an ocean/wave model to simulate total water elevation caused by ocean storm tide, wave, freshwater flow, and wind forcing in both 1D and 2D domains. In the first approach, the inland hydrology model needs to be capable of receiving coastal-related signals (i.e., water surface elevation) and propagate it landwards such as backwater effect. This also requires that the coastal ocean model be able to accurately receive and resolve freshwater inflow and vertical exchanges in wet and dry parts of its computational domain (i.e., local generated runoff, evaporation, and infiltration). Most of ocean models cannot adequately simulate streamflows in the rivers and tributaries with the temporal and spatial resolution required for this analysis. Therefore, in this study, the first approach is not considered and the focus is placed on the development of the latter two approaches.
The one-way coupling between an ocean/wave model (ADCIRC/WW3) and a hydrologic model (NWM) under three major storms (Hurricanes Isabel, Irene, and Sandy) is validated by incorporating two well-established hydraulic (Hydrologic Engineering Center's River Analysis System [HEC-RAS]) and hydrodynamic (D-Flow Flexible Mesh [D-Flow FM]) models. HEC-RAS is a hydraulic model for water flow through natural rivers, and D-Flow FM is a modeling suite for hydrodynamic simulation on unstructured grids in combined 1D/2D. Either D-Flow FM or HEC-RAS can communicate with the NWM and ADCIRC by passing the output discharge time series from the NWM and water level/velocity time series from ADCIRC.
The study area here is Delaware Bay, which is one of the largest and most important bays along the U.S. East Coast. It is also a generalizable example of shallow coastal estuaries (Pareja-Roman et al., 2019). The hydrodynamic characteristics in the Delaware Bay are determined by the interaction of the tidal flow from the Atlantic Ocean and the freshwater flows from the mainstem of the Delaware River and its tributaries (Dupont, 2019). The results of the simulations conducted in the coupled modeling framework are compared with the available field observations of NOAA Center for Operational Oceanographic Products and Services (CO-OPS; https://tidesandcurrents.noaa.gov/). A sensitivity analysis is also carried out to assess the effects of streamflows and wind, which are the two variables that have a large impact on the model performance.
The description of the modeling approach is presented in section 2. Section 3 contains the model setup, model domains, and boundary conditions. Model validation for three major storms and different modeling approaches as well as sensitivity analysis are presented in section 4. Finally, the summary and conclusions are provided in section 5.

Description of the Model Couplings Tested
A new 1D/2D coupled modeling approach is introduced to predict total water levels and flood inundation in the coastal areas of the CONtinental United States (CONUS) aiming to improve the freshwater discharge forecasts of the NWM. The one-way coupling between the ADCIRC and the NWM at the transition zone between the inland and coastal zones is provided by riverine-estuarine models. In this study, two different riverine-estuary models are used: (i) the HEC-RAS hydraulic model and (ii) the D-Flow FM hydraulic/hydrodynamic model. The D-Flow FM is capable of seamlessly integrating 2D physics (e.g., estuaries, bays) with 1D hydraulics (e.g., rivers, tributaries). The effects of the wind-generated surface waves are modeled using the third-generation WW3 spectral wave model, which is two-way coupled with ADCIRC. Figure 1 shows a schematic representation of the proposed modeling approach. Summary details of all modeling components are provided in sections 2.1 and 2.2.

Riverine-Estuary Hydrodynamic/Hydraulic Models
The reliability of different 1D/2D hydrodynamic/hydraulic models in the Delaware River Basin (DRB) for simulating of inland-coastal compound flood propagation and inundation/water levels is examined. The performance of D-Flow FM and HEC-RAS models can be evaluated by comparing observed and corresponding simulated water levels, using different statistical measures.

Hydrologic Engineering Center's River Analysis System
HEC-RAS (https://www.hec.usace.army.mil/software/hec-ras) is a public-domain model that can be used to perform water flow simulations. This model was developed by the Hydrologic Engineering Center of the U.S. Army Corps of Engineers (USACE) and solves the dynamic Saint-Venant equations for unsteady flow (HEC-RAS, 2016). HEC-RAS models the hydraulics of water flow through natural rivers and other channels. HEC-RAS can be used to predict the extent of flood and create inundation maps. It may have some numerical instability problems in steep and/or highly dynamic rivers and streams.

D-Flow FM Model
D-Flow FM (http://oss.deltares.nl/web/delft3d/d-flow-flexible-mesh) is a new hydrodynamic model that solves the shallow-water equations on an unstructured mesh in a coupled 1D/2D fashion. Together with the curvilinear grids from Delft3D, the unstructured mesh can contain pentagons, triangles, and 1D river networks altogether in single mesh (D-Flow FM User manual, 2014). D-Flow FM is composed of Delft3D 4 and SOBEK 2 hydrodynamic models and contains unstructured grids. D-Flow FM is a very flexible mesh in the coastal transition zone and is suitable for tidal, estuarine, riverine, and inundation modeling by using the wetting-and-drying scheme (D-Flow FM User manual, 2014). The model achieves high performance by the use of multicore architectures, grid computing clusters, and subgrid methods. D-Flow FM is an open source code. A detailed conceptual description of the model and numerical approach was given in the D-Flow FM user manual.

Boundary Conditions and Atmospheric Forcing
In this study, an efficient couple system of a 2D ocean circulation model, ADCIRC (Luettich et al., 1992), and a third-generation spectral wave model, WAVEWATCH III WW3DG, 2019), is used to compute the ocean boundary condition. The radiation stresses from wave calculations are communicated to the ADCIRC ocean model and used in computing the water level. Flow field is then communicated back from ADCIRC to WW3. This exchange of data happens interactively . The simulations are conducted on an unstructured triangular mesh with~1.8 M node with lowest resolution of~200 m near the coast. The atmospheric forcings (wind and pressure) are provided to the system from the Hurricane Weather Research and Forecasting (HWRF) model. The water level and velocity time series at the offshore boundary are extracted and imposed to the hydraulic model as the ocean boundary condition (Figure 2). The inland boundary conditions at the river and tributaries are time-varying streamflow discharges that are computed by the NWM (http://water.noaa.gov/about/nwm). The NWM produces a variety of operational hydrologic analysis and prediction products for the CONUS (e.g., estimates of river water velocity and flow discharges for more than 2.7 million river reaches that include both data-rich and underserved locations). Figure 3b shows the study area and the model domain. The study area covers the Delaware River reach and its 10 primary tributaries from the upstream boundary (USB) to the Atlantic Coast ( Figure 3b). The USB of the river flow is downstream of Trenton, NJ, where the non-tidal flows meet tidal flows of the Delaware Bay. Each tributary is composed of a corridor under the influence of backwater effects from the Delaware River or Delaware Bay (see Table 1). Generally, the 11 tributaries contribute about 95% of the freshwater into the Delaware Bay. The Chesapeake and Delaware Canal is not considered in this study. The river discharge pattern in the Delaware River is seasonal. This is primarily from the snow melt in the Catskill and Pocono mountains. Based on a 27-year data set, the average discharge at Trenton, NJ, for June-October was 195 m 3 s −1 , for November-February was 334 m 3 s −1 , and for March-May was 510 m 3 s −1 (Sharp et al., 1986;Smullen et al., 1983). According to the U.S. Geological Survey (USGS) record at Trenton, NJ, station (https:// waterdata.usgs.gov/usa/nwis/uv?01463500), mean annual discharge range is between 140.93 to 677.91 m 3 s −1 from 1913 through 2020. The mean freshwater discharge varied from 330 m 3 s −1 at the upstream to 570 m 3 s −1 at the mouth of the bay. At the head of the tides, the freshwater discharge varies from 34 to 8,500 m 3 s −1 during the severe drought and wet periods, respectively (Partheniades, 2009). The question that arises is how to determine the minimum length of each tributary that is necessary to capture the essential physics in the model. This minimum corridor length (1D portion) must be longer than the maximum possible backwater extents in Table 1 ensuring that the actual spatial extent of flooding due to both tidal influences and freshwater fluxes is properly computed.

Study Area, Model Setup, and Model Domain
Using geospatial tools, a subset of the NWM river network coinciding with prominent river tributaries is selected. Next, the NWM hydrologic forecast locations that are collocated with the 11 USB limits of the tributaries are identified. A complete set of the retrospective streamflow data from the NWM for the respective named storm events is obtained. This data set forms the lateral boundary conditions for hydraulic/hydrodynamic models in DRB.

D-Flow FM Setup
The D-Flow FM model uses the unstructured mesh that covers the entire DRB. The mesh covers 250 km of the lower Delaware River (Figure 3a). It consists of a combination of 2D curvilinear/triangular mesh and a 1D river network including several cross sections. Two coupled 2D/1D setups are used: (i) 2D mesh for the bay and 1D network just for the main river (Delaware), and using discharge boundary conditions at other lateral rivers (Figure 3b), (ii) 2D grids for the bay and 1D networks for all tributaries (Figure 3c). The mesh includes sufficient resolution to represent realistic physics of rivers and the bay. The 2D mesh has 57,872 elements in total, with mesh element sizes varying from 50 m to 3 km. The mesh is refined in/near the river channel and gradually coarsens over the floodplains and towards the ocean. The upstream segment of the river, approximately 120 km, is represented with the 1D river network.
After the mesh is created, elevation data are assigned to each node ( Figure 3a). The primary source of the topo/bathymetric information is the USGS Coastal National Elevation Database (CoNED). The CoNED is a seamlessly integrated topo/bathymetric database from multiple topographic data sources with adjacent intertidal topo/bathymetric and offshore bathymetric sources (Danielson et al., 2016). Elevation data of the nodes located in the ocean were assigned using the ADCIRC bathymetric data set. The topo/bathymetry values have World Geodetic System (WGS84) as the horizontal datum and geodetic North American Vertical Datum 1988 (NAVD88) as vertical datum.
Spatially variable Manning's roughness coefficients are initially derived from the multiple studies that have been done by USACE, USGS, and Federal Emergency Management Agency (FEMA). D-Flow FM interpolates roughness values during the model initialization process, in order to achieve the smooth transition of roughness values over the model domain. Some sensitivity studies based on tidal analysis are done in order to achieve the best roughness configuration.

HEC-RAS Setup
The digital terrain model (DTM) is compiled in the form of a triangular irregular network (TIN) for use in HEC-RAS model development. The FEMA 500-year flood coverage plus a buffer zone around it is used to clip the TIN coverage. Stream cross-section information is generated using the HEC-GeoRAS/ArcInfo tools. Additional data not extracted by the HEC-GeoRAS/ArcInfo such as detailed channel roughness coefficient and ineffective areas are later defined in the HEC-RAS model. The stream network is composed of stream centerlines for the main stream Delaware River and all its major tributaries ( Figure 4). Stream centerlines identify the connectivity between different streams for flood routing within HEC-RAS.
The key inputs to HEC-RAS model are DRB geometric data, Delaware River floodplain data (length and elevation), the distance between consecutive river cross sections, Manning's roughness coefficient, and slope. The downstream boundary is located at a nearby ADCIRC point in the mouth of the bay. The USB of HEC-RAS is located at the USB of the river (i.e., downstream of Trenton, NJ). Ten internal boundary conditions (coincident with the NWM forecasts at respective tributaries) are also included.

Coupled Model Setup and Boundary Conditions
The ocean boundary condition for D-Flow FM is the time series of water elevations/velocities that are calculated by the two-way coupled ocean model ADCIRC with the wave model WW3. The inland boundary conditions at the river and tributaries are time-varying streamflow discharges that are computed by the NWM. Different wind and pressure reanalysis and hindcast data are applied as surface boundary forcing. The bottom Manning's roughness coefficient is spatially varying. D-Flow FM is run with hourly wind and pressure forcings on the unstructured mesh, while combined effect of storm surge and flood tide (tides/storms) and discharges are imposed as boundary conditions. As such, the simulations and model validations focus on estuarine and nearshore areas where tidal currents can significantly influence the water level, as well as upstream where the river inflow has an influence on water levels. The simulation period in all model runs is the same as ADCIRC simulation time frame (22, 23, and 14 days for Hurricanes Isabel, Sandy, and Irene, respectively). A spin-up period of 10-13 days prior to the storm landfall is used to remove initial condition effects. D-Flow FM implements a finite volume solver on a staggered grid. The simulations are run on the HPC cluster (a Linux-based computing system) at NOAA HPC Center (Theia) and personal desktop. In serial mode (desktop), the model run time per storm and for this domain using one CPU is about 1.5-2.5 h, while in an HPC parallel environment using eight cores, it takes about 15 min. We note here that D-Flow FM has been parallelized and it helps for much faster simulation times.

Atmospheric Forcing
The atmospheric forcings of the D-Flow FM simulations during the three named storms are supplied by two low-and one high-resolution wind/pressure products: (i) the Global Forecast System (GFS; https://www. ncdc.noaa.gov/data-access/model-data/model-datasets/global-forcast-system-gfs), (ii) the Climate Forecast 10.1029/2019JC015822

Journal of Geophysical Research: Oceans
System Reanalysis (CFSR; https://rda.ucar.edu/datasets/ds093.1), and (iii) the HWRF (https://www.emc. ncep.noaa.gov/gc_wmb/vxt/HWRF/index.php) (hourly with three inner nested moving grids with the finest resolution of 0.02°moving with the hurricane eye). All atmospheric products contain the data for wind speed components (ms −1 ) at 10 m above ground and atmospheric pressure (Pa) reduced to mean sea level (MSL). The domain of interest is confined between 39-41°N and 72-76°W. The three wind products are used to assess the sensitivity of the modeling system to different atmospheric forcing resolutions and physics and how these affect the simulations of the total water level. Details about the above wind products are presented in Table 2. For modeling the Hurricane Isabel and Super Storm Sandy, GFS, CFSR, and HWRF wind products were used, while for Hurricane Irene, we used CFSR and GFS wind products.

Observations
The calculated water levels are compared with the available field observations of NOAA CO-OPS at different sites from Atlantic City, NJ, to Newbold, PA. In order to select stations from areas with different dynamics (in the ocean, inside the bay, upstream the river) based on the availability of data sets, eight stations are selected for model validation (Figure 5a).

Tropical Cyclones
Simulations of three major and notable storms that have made landfall in/near the DRB (i.e., Hurricane Isabel, 2003;Hurricane Irene, 2011;and Super Storm Sandy, 2012) are used to evaluate the one-way coupled modeling approach. While physical processes in these areas are influenced by all three storms, each storm has unique and distinctive characteristics. Hurricane Isabel, a Category 5 major hurricane (

Results and Discussion
Different riverine-coastal modeling systems (NWM/D-Flow FM/ADCIRC/WW3 and NWM/HEC-RAS/ ADCIRC/WW3) are evaluated to confirm that models satisfactorily simulate DRB inland-coastal compound flooding during three major storms. Tidal analysis, time series comparisons of water levels at different stations for various storms, as well as sensitivity analysis of riverine discharges and atmospheric forcing for various locations and storms are presented.

Tidal Analysis
First, the Delaware River mainstem extent from upstream to the Atlantic Coast (approximately 480 km) was selected for the tidal analysis. Five USGS gauge stations (U1 to U5) and eight NOAA CO-OPS (C1 to C8) were considered in the study (see Figure 5a). Tidal analysis (T-TIDE, Pawlowicz et al., 2002) was performed on the water level data for the temporal window from 1 July 2016 to 30 November 2016. During this time frame, the river was subjected to neither extreme effects of wind nor large freshwater flows. Here, the extent of the tidal signature upstream was determined. The amplitudes of the significant tidal constituents were considered. The spatial variability of the dominant constituent, M 2 , is shown in Figure 5b. The magnitude of the M 2 component increases gradually from Delaware Bay (CO-OPS #8557380, C1) to a maximum value near Trenton, NJ (CO-OPS #8548989, C8). As the tide advances upstream from C8, it encounters the fall line near Trenton-Morrisvile Toll Bridge on Hwy 1, and at this point, the tidal propagation stops. Therefore, it is the proper location for transitioning from the 2D domain to the 1D network.
Subsequently, the D-Flow FM model was verified to test stability and general performance when it is forced with only tidal components. The D-Flow FM model was run for a 1-year period using known astronomical tidal components as the downstream boundary condition and minimal flow of 10 m 3 s −1 for the USB condition. For tidal calibration, the ADCIRC Western North Atlantic Tidal Database, the Eastcoast 2001 database, was used to extract tidal constituents. The Eastcoast 2001 database describes the calculated phase and amplitude for different tidal constituents (i.e., M 2 , S 2 , N 2 , K 2 , O 1 , K 1 , Q 1 , M 4 , and M 6 ) (Mukai et al., 2002).
The ability of D-Flow FM to reproduce tides at the NOAA CO-OPS stations was evaluated using the automated constituent selection algorithm, T-TIDE. Using T-TIDE, the harmonic constituents were added, giving the total tidal variation at a particular coastal location independent of the effects of wind, freshwater flows, and other deviations from the pure astronomical tide. In addition, a 95% confidence interval for the estimated tidal parameters was generated through non-linear mapping. The computed constituents were then compared with observed constituents at the collocated NOAA CO-OPS locations (see Figure 5a).
Then, some sensitivity studies based on tidal analysis were done in order to achieve the best roughness configuration. To obtain and fine-tune the acceptable roughness values, two factors were considered: (i) physical processes that cause energy losses such as turbulence and river meandering and (ii) deepening of the mainstem by USACE to ensure safe ship transit. Deeper channels have lower flow resistance/ roughness than the rest of the river. Spatially variable Manning's roughness coefficients were initially derived based on land use and land cover data (Figure 5c) from the multiple studies that have been done   (1) and (2). The fine-tuning started from the ocean boundary at the river mouth towards the transition zone (zone 1 of Figure 5d). According to tidal analysis, and in order to capture the most accurate tidal amplitudes and tidal phases for different stations in this zone, the roughness values were modified. Then, we started from very upstream (near Trenton) towards the bay and again modified the roughness in order to obtain the best tidal constituents. Finally, the spatial distributions of roughness in zone (3) were determined (Figure 5d). Therefore, multiple roughness values were specified along the river and the lateral cross sections. D-Flow FM interpolates roughness values during the model initialization process, in order to achieve the smooth transition of roughness values over the model domain. In subsequent modeling efforts, no attempt was made to recalibrate roughness values for changes in flow. Figure 6 shows the harmonic constituents of tidal amplitudes and tidal phases for different stations in the DRB using the best acceptable roughness values. Clearly, the D-Flow FM accurately captures the characteristics of tidal dynamics. Comparison of model outcomes and measurements show strong correlations at all locations. The small mismatch between computed and observed components is an artifact of non-tidal or residual "noise" in the data. The "true" values are within the confidence limits. Based on the verified model setup and fine-tuned roughness values, model validation is presented in the following section.

Model Validation for Three Extreme Storms and Different Riverine Models
In this section, the model results (i.e., detailed water level time series) are compared with field measurements of three major storms (i.e., Hurricanes Sandy, Isabel, and Irene). Table 3 shows the required input data for both models. HEC-RAS does not use flow velocity and atmospheric forcing.

Super Storm Sandy (2012)
Super Storm Sandy was an extraordinary TC in the New Jersey/New York area due to its large size, hurricane-strength winds, resulting large waves, and record-breaking flood height (Orton et al., 2016). The results of two different modeling approaches (NWM/D-Flow FM/ADCIRC/WW3 and NWM/HEC-RAS/ADCIRC/WW3) for simulating this hurricane are presented and compared.

D-Flow FM Model
The total simulation period was 25 days including the spin-up period of 13 days prior to the storm landfall that was used to remove initial condition effects beginning on 9 October 2012. Then, the storm simulation period was 12 days beginning on 22 October 2012. After developing D-Flow FM model with best available grid setup/size and wind product, roughness coefficient was the most important calibration parameter affecting the calculation of surface elevation and velocity distribution.   Figure 8 show comparisons of the NWM/D-Flow FM/ ADCIRC/WW3 modeling system's (see Figure 3b) results and field data for water level. In addition, Table 4 gives statistical measures (i.e., bias, R, root mean square error [RMSE], and skill) of observational and numerical results for Super Storm Sandy at four different stations. The equations for these measures are (Willmott, 1981):   In this study, the combined effects of river flow, tide, wind, and wave on total water levels were investigated. Furthermore, one of the main goals is to achieve accurate water level from the compound flooding via an accurate numerical framework under different weather conditions such as stormy and normal condition. In addition, lots of the peak water level events do not coincide with hurricane direct impact period. Therefore, simulation times before and after the direct impact of the hurricane were included in calculations of statistical measures, so that we can investigate the model's capability to produce total water levels based on all forcings available. As illustrated in Figure 8 and Table 4, the model results and field observations agree qualitatively. The temporal pattern of variation in water level is well predicted (with most R ranging from 0.92 to 0.97 and skill ranging from 0.82 to 0.85), which provides confidence for this coupled model setup. The measured phase and amplitude time series are well reproduced by the model. The model simulations during the peak of extreme storms are important for flood investigations. The present modeling approach is powerful for this aspect. The small discrepancies between modeled and observed values may be due to errors in the wind product, Manning's roughness coefficient, or topo/bathymetry. The middle panels of Figure 8, as well as Table 4, show the comparisons of the NWM/D-Flow FM/ADCIRC/ WW3 model results and field data for water level for 2D/1D/All-Rivers set up (according to Figure 3b). As stated in Table 4, both 2D and 1D setups have reasonable results with marginally better results for 2D/1D/All-Rivers. According to the 2D/1D/All-Rivers (Figure 3c) simulations, skill and R are larger and RMSE is lower than those of the 2D/1D/Delaware (Figure 3b), especially for upstream stations. This reflects a better ability of the 2D/1D/All-Rivers model than 2D/1D/Delaware to simulate the hydrodynamic characteristics.

Journal of Geophysical Research: Oceans
As the ADCIRC mesh enters the bay (Figure 2), the ADCIRC outcomes in the overlapped area were evaluated and compared with the D-Flow FM ones. Figure 9 depicts the water level prediction using ADCRIC model compared to the NOAA-observed data during Super Storm Sandy (2012). Table 4 and Figures 8 and  9 show that the accuracy of D-Flow FM models is greater than that of the ADCIRC especially in the overland area, with smaller values for the RMSE, while coefficients of determination and skill are closer to unity. The main D-Flow FM strength above ADCIRC and some of the well-known hydrodynamic/ocean circulation models (e.g., Regional Oceanic Modeling System [ROMS] and The Unstructured Grid Finite Volume Community Ocean Model [FVCOM]) is the ability to blend 1D computational domain with 2D and 3D within one modeling framework.

HEC-RAS Model
In order to compare the results of different riverine models, the results of NWM/HEC-RAS/ADCIRC/WW3 are presented in this section. The model is first configured and debugged using steady flow computations. The steady flow runs enable us to see potential structural problems in the model at particular flows. For steady flow modeling, USB conditions are input as discharges from the 10-, 25-, 100-, and 500-year recurrence interval storm events. The geometric information was adjusted in order to obtain a more robust hydraulic model.
The adjusted HEC-RAS model is then applied to the unsteady state simulation using the NWM forecast discharge as the upstream and lateral boundary conditions. Calibrated water elevation data are obtained from the NOAA CO-OPS. There are seven intermediate CO-OPS stations with stage information. If the initial runs do not closely agree with the observed data at the above-mentioned locations, one of the HEC-RAS parameters is modified within reasonable values of the parameter. A new run is conducted, and the results are again compared to the known data. This process is continued until the calibration is within a reasonable To evaluate the performance of the HEC-RAS model, the same criteria as for the D-Flow FM model are used. Comparison between simulated and observed water levels at various stations for Super Storm Sandy using HEC-RAS is presented in the lowest panels of Figure 8. Simulation results show that peak water level predictions are generally in good agreement; however, time of occurrence is not well predicted. Table 4 provides the statistics for model-observation comparisons for both the D-Flow FM and HEC-RAS models. Table 4 illustrates that the water levels can be calculated reasonably well using the D-Flow FM model as compared to HEC-RAS. In addition, the accuracy of D-Flow FM models is greater than that of the HEC-RAS, with smaller values for the RMSE, while coefficients of determination and skill are closer  Unbiased models are only superior to those with higher bias if the unbiased models also have superior precision and overall smaller variability. Generally, biased models have smaller total error than unbiased models. As bias is the average difference between the simulated and the observed values, clearly, it would not be the only measure for model effectiveness. Therefore, various criteria must be considered for a model to be designated superior to another. Furthermore, the HEC-RAS did not simulate tide well before and after hurricane made landfall (when water levels were mostly tidal), indicating that this model lacks the capability to be a good storm surge or tidal model. Based upon these results ( Figure 8 and Table 4), we may conclude that D-Flow FM is the better model. The improvement in the D-Flow FM over the HEC-RAS may be due to the wind-modeling capability used in the former, which might be crucial for the estuarine hydrodynamic simulation under extreme storms. Additionally, HEC-RAS did not perform well due to model's inability to simulate tide accurately. Mashriqui et al. (2014) reached a similar conclusion from their studies of estuarine river reach of the Potomac River.
In order to provide floodplain mapping, water surface profile data are exported from HEC-RAS simulations and processed by the HEC-GeoRAS. The HEC-GeoRAS is a set of procedures, tools, and utilities in the HEC-RAS for processing geospatial data. Next, the D-Flow outputs were used to calculate inundation extent over the terrain. D-Flow water surface elevation outputs were interpolated and extrapolated to generate a continuous water surface elevation profile. The water surface elevation profile was then subtracted by a 10-m digital elevation model with the positive values depicting inundation extent. Figure 10 depicts the inundation map (using both HEC-RAS and D-Flow FM models) that shows the areal extent of the flooding and the resulting water depth at each location. The inset in Figure 10 presents a snapshot of the simulated flooding along the Rancocas tributary, which shows the propagation of ocean forcing in the lateral rivers' streamflows. The Delaware Bay system is a riverine-estuarine transition zone where the water elevation is greatly influenced by the freshwater component, oceanic tides, and atmospheric forcings. During extreme events, these multiple influences are all important. Despite the advantages of the HEC-RAS implementation in normal events, it underpredicts flooding extents during extreme events, such as Super Storm Sandy (2012), because it lacks a wind forcing term in its momentum equation, which is evident in Figure 10a. The wind forcing term has been included in the D-Flow FM formulation (see Figure 10b) resulting in the larger flood inundation extent that matches high water marks. These simulation results suggest that implementing 1D is necessary because it provides vital water information in the dense urban corridors along the tributary during high-impact events.

Hurricane Isabel (2003)
Hurricane Isabel (2003) will be recalled because of its intensity, size, and impacted areas. Isabel was the most expensive storm to hit the U.S. coasts in 2003 (Preller et al., 2003). Extreme rainfall from the storm occurred in the East Coast. According to NWS reports, this rain caused flash floods over several tributaries and rivers. As in the previous section, two different modeling systems (i.e., NWM/D-Flow FM/ADCICR/WW3 and NWM/HEC-RAS/ADCIRC/WW3) are validated using the observed time series data at different stations in DRB.

D-Flow FM Model
The total simulation period was 25 days including the spin-up period of 12 days prior to the storm landfall beginning on 29 August 2003. Then, the storm simulation period was 13 days beginning on 10 September 2003. Figure 11 illustrates the water level prediction compared to the NOAA observed data during Hurricane Isabel (2003). Temporal variations of observational and predicted water levels at various stations in the DRB (2D/1D/Delaware) are shown in upper panels of Figure 11. Table 5 gives statistics for water levels. The model skill in all of the top panels is higher than 0.84, the overall magnitude and trend of the simulated water level time series match with observations, and inconsistencies are mostly small. It can be seen that the RMSE between the measured and predicted water levels are increased for stations farther upstream. In addition, the model overpredicted the water levels for Ship John Shoal (NJ), Burlington (NJ), and Newbold (NJ) (bias ¼ 0.02, 0.12, 0.10 m, respectively) and underpredicted the water level for Marcus Hook, PA (bias ¼ −0.015 m). The higher RMSE values for upstream stations may be due to the fact that in the narrower and shallower part of the river, different motions with various time scales and physics exist. The disagreement between modeled and observed water levels could also be related to the accuracy of wind and streamflow inputs, error in roughness estimation, or observations.
The middle panels in Figure 11 and Table 5 show the comparisons of the NWM/D-Flow FM/ADCIRC/WW3 model results and field data for water level for Hurricane Isabel (2003), for 2D/1D/All-Rivers set up (according to Figure 3b). Water level results have a typical bias around −0.02 to 0.12 m and skill ranging from 0.88 to

10.1029/2019JC015822
Journal of Geophysical Research: Oceans 0.93 inside the bay and in upstream stations. Albeit some differences can be seen between observations and numerical simulations, the model predictions are generally very good. Water level simulations were mostly accurate near the peak of each storm. Furthermore, for most of the stations, the model produces an overestimation of water level. According to Table 5, both 2D and 1D setups have reasonable results with marginally better results for 2D/1D/All-Rivers. For the 2D/1D/All-Rivers (Figure 3c) simulations, skill and R are higher, and bias and RMSE are lower than those of the 2D/1D/Delaware (Figure 3b), especially for upstream stations. This confirms a better ability of the 2D/1D/All-Rivers model to simulate the hydrodynamic characteristics than 2D/1D/Delaware. This modeling architecture (2D/1D/All-Rivers) promotes sufficient dynamic information exchange between the river and estuary and hence realistically addresses a wider range of physical processes.

HEC-RAS Model
Comparison between simulated and observed water levels at various stations for Hurricane Isabel using HEC-RAS is presented in the bottom panels of Figure 11. Simulation results show that peak water level predictions are generally in good agreement; however, the time of occurrence is not well predicted. Table 5 shows the statistics of observational data and numerical results for D-Flow FM and HEC-RAS models over the different regions of the domain (in the bay and upstream of the river). The results show good temporal correlation with the measured data for D-Flow FM and confirm that D-Flow FM is evidently the better model for water level estimation in the bay. This indicates that the model can be used as a better tool to predict future flooding during major storms. Comparison of the calculated statistical values shows that the accuracy of the D-Flow FM model is better than that of the HEC-RAS model; the RMSE of the water levels using D-Flow FM is less than 50% of that using HEC-RAS. In addition, by moving upstream, the differences between the two models are increased, which shows the ability of D-Flow FM to capture the physics of the rivers both near and far from the coasts.
The superiority of the D-Flow FM over the HEC-RAS model may be related to the exclusion of atmospheric components from the simulation, which is crucial for its performance in resolving the water movement during the storms near the coastal areas. Another possible reason may be other simplifications to the flow equation, which render the HEC-RAS less suitable for the estuarine-riverine system. Furthermore, HEC-RAS cannot accurately predict the storm surge and tide in areas like the bay and its tributaries. Figure 11 shows that when the water levels were mostly tidal (before and after the hurricane landfall), HEC-RAS-simulated tide was small in amplitude and out of phase.

Hurricane Irene (2011)
Although Hurricane Irene's maximum intensity did not go above Category 3, it was a very large storm. Due to the scale of its impact on people and property, Hurricane Irene (2011) has been the subject of earlier studies (Mooney et al., 2019;Yablonsky et al., 2015) for regional climate modeling framework. Hurricane Irene (2011) caused extensive and significant impacts such as severe flooding, widespread wind damage, and power outages along the U.S. East Coast (Klausmann, 2014). Hurricane Irene had a high moisture content that carried very heavy rainfall and streamflow rates to the U.S. East Coast.
The total simulation period was 15 days including the spin-up period of 7 days prior to the storm landfall beginning on 14 August 2011. Then, the storm simulation period was 8 days beginning on 21 August 2011. For Hurricane Irene, the high-resolution HWRF data for ADCIRC/WW3 simulations are not available. Therefore, wind and pressure fields for ADCIRC/WW3 model are produced based on a parametric wind model (Holland, 1980), which calculated barometric pressure and wind stress according to best-track estimation from National Hurricane Center (NHC).
Comparison between simulated and observed water levels at various stations for Hurricane Irene (2011) using D-Flow FM and HEC-RAS is presented in upper and lower panels of Figure 12, respectively. As can be seen from Figure 12, model results during the storm are generally reasonable, while after the storm peak, they are less accurate. Even though some differences can be seen between observations and numerical simulations, skills were usually higher than 0.80. As measured by skill, the water level results are less accurate than the results for other storms (Hurricanes Isabel and Sandy), likely due a more simplified wind product that does not capture the full physics. Statistical measures (Table 6) show that the model is moderately accurate in computing values of water levels but less accurate in capturing the variability (measured by skill). Figure 12 shows that by moving upstream, the capability of HEC-RAS decreases drastically. Furthermore,

Journal of Geophysical Research: Oceans
it does not show good temporal correlation with the measured data. Comparisons of the upper and lower panels in Figure 12, in addition to Table 6, show that the D-Flow FM is more appropriate for water level predictions in Delaware Bay than the HEC-RAS.
HEC-RAS is not a parallel model and is not able to resolve wind forcing and ocean tide, which can play a dominant role in changing total water level in coastal regions. Therefore, D-Flow FM would be a better candidate for riverine-coastal coupling. In addition, D-Flow FM has parallel computation capabilities and the potential for inclusion within the NWM approach. D-Flow FM can be implemented for a single river, estuary, or a coastal zone with ocean tide, freshwater, and wind forcing in both 1D and 2D domains.

Sensitivity Analysis of Atmospheric Products and Streamflows
Based upon the results for the three extreme storms in the previous sections, we concluded that D-Flow FM is superior to HEC-RAS. Therefore, in the rest of the paper, the former model is considered. In this section, the wind and streamflow impacts on water level are investigated.

Atmospheric Effects on Water Level
D-Flow FM is highly flexible and directly handles ASCII and NetCDF meteorological input formats allowing the use of different sources of meteorological data. Here, D-Flow FM is driven by three different wind products as discussed in section 3.4 and presented in Table 2 (i.e., GFS, CFSR, and HWRF). The quality of the predicted water levels, especially under the influence of significant wind variations, depends greatly on the quality and accuracy of the atmospheric forcing fields. The uncertainties introduced by atmospheric prediction models include errors related to the modeling of the physical processes, the spatial resolution, and the quality of the initial conditions (Krzysztofowicz, 2001).  In order to drive D-Flow FM with the best available wind forcing, sensitivity tests are conducted to determine the relative importance of the spatial and temporal resolution of the wind products on the model predictions. Figure 13 displays the comparison between the modeled wind speed and the observations for three different wind products for Super Storm Sandy (2012) and for the CO-OPS observation station Ship John Shoal, NJ, located at the mouth of Delaware River. The results show that among the three wind products, HWRF produces the best representation of the observations. From the three products, HWRF best coincides with the observations, especially during the passage of the storm. The other two low-resolution products capture the wind trend reasonably well, but do not capture the wind variations (peaks) during the passage of the storm as well as HWRF. Figure 14 depicts the comparisons between modeled water level predictions and NOAA observed data for different wind forcing products during Hurricane Sandy (2012), while Table 7 presents statistics for the predicted water levels during the Sandy simulations using different wind forcing products. Figure 14 and Table 7 show that HWRF is superior to GFS and CFSR for water level prediction during extreme storms.

Streamflow Effects on Water Level
Delaware estuary is a partially mixed estuary with salinity intrusion-controlled shoaling. The estuary length is 214 km (Partheniades, 2009). In order to examine the variations in simulated discharge along the river/bay, D-Flow FM is forced with inflow times series from the NWM retrospective model during the storm. The NWM is used to route the flood and determine expected flow hydrographs at the point of interest (11 available tributaries).
To see the influences of streamflow and oceanic fluctuations on the water levels, the cross-section discharge at different locations in both 1D and 2D sections during Hurricane Isabel (2003) was computed ( Figure 15). Importantly, the shape and timing of the cross-section discharge at different locations in 1D network (#2 to #4) are accurately simulated with nearly identical characteristics as the observed hydrographs, as illustrated in location (#1). From location (#3) towards the bay, the effects of ocean forcing on the streamflow can be seen. Specifically, from location (#4), some negative values that represent the effects of backwater in the river are visible. In location (#7), the water level and discharge are determined primarily by oceanic fluctuations.
These results suggest that total water depth from Trenton, NJ, to Marcus Hook, PA (#3 to #7), is determined by the interaction of the tidal flow from the Atlantic Ocean and the freshwater flows from the mainstem Delaware River and its tributaries. Downstream of Marcus Hook, PA (#7), the controlling factor for total water level is the oceanic forcing and associated meteorological conditions. Generally, the tidal flow is 300 times greater than the freshwater flow. Lastly, water depth upstream of Trenton, NJ (#1 to #3), is dependent on freshwater flow only.
In this study, salinity field was not computed by the model. A long-term analysis of mean axial salinity distribution and Delaware River discharge indicates that Delaware is a weakly stratified estuary with a typical vertical salinity variation of only 1 psu. In addition, the response of salinity to Delaware River discharge is weak, more enhanced during flood tide relative to ebb tide. The weak salinity response reduction with freshwater is due to the action of vertical shear flow dispersion in a tidally stirred regime and the action of lateral shear coupled to strong lateral salinity gradients which are driven by the lateral flows governed by the  (Figure 3b) 2D/1D/All tributaries (Figure 3c)
In order to see the effects of streamflow on the water levels, two different scenarios with two different streamflow time series are run: (i) Hurricane Isabel with actual streamflows and (ii) Hurricane Isabel with double streamflows from Delaware River. Then, the water level differences between these two scenarios are shown in Figure 16. At the upstream stations (Newbold, PA, and Burlington, NJ), the freshwater flux entering the tidal reach is significant relative to the contributions from tidal effects. At this river corridor, both the freshwater component and tidal forces influence the water level. Nevertheless, at Marcus Hook, PA (3), the maximum rate of tidal flow in the Delaware Bay proves to be very large, as compared to the freshwater The panels in Figure 16 show that among the stations, river streamflow has the most significant effect on water levels at upstream stations. These findings are consistent with the existing understanding of the relative importance of the river discharges on water level estimations. In addition, these results are consistent with the changes in the magnitude of discharges. In other words, significant changes in the discharges of rivers during/after storms show the significant effect on flooding for upstream areas and small changes for the areas near the coast. This test confirms that freshwater is a principal component of the total water depth at upstream river reaches. However, the water level in the bay is dependent on the tidal effects.
It is worth mentioning that in reality, lots of the peak water level events do not coincide with hurricane events. The discharges in Delaware River at Trenton are in the range of 113.3-1132.7, 84.6-1415.8, and 84.6-4247.5 m 3 s −1 for Hurricanes Sandy, Isabel, and Irene, respectively, while historical peak value from 1982 to 2020 is approximately 6,800 m 3 s −1 . However, since one of our main goals was to achieve accurate water level from the combined storm tide and freshwater inputs (compound flooding) via an accurate numerical framework, we compared our model numerical results with actual observations during extreme storms. These results indicate the strengths of the proposed modeling approach for simulations of flooding during the most extreme storms, including for use in future operational forecasting systems.

Summary and Conclusions
This paper demonstrated one-way coupling between NOAA's freshwater model (NWM) and a coupled wave/ocean model (WW3/ADCIRC) at the DRB using different inland hydraulic/hydrodynamic models (i.e., HEC-RAS and D-Flow FM) to determine the most appropriate modeling system (NWM/D-FLOW FM/ADCIRC/WW3 or NWM/HEC-RAS/ADCIRC/WW3). First, the modeling frameworks were verified by performing tidal analysis. The spatial variability of roughness was optimized, and the accuracy of     Journal of Geophysical Research: Oceans topo/bathymetry was tested. The characteristics of tidal dynamics were well captured by the D-Flow FM. Then, modeling system was validated for the cases of Hurricane Isabel (2003), Hurricane Irene (2011), and Super Storm Sandy (2012), and numerical results were compared to CO-OPS observational data.
In order to investigate flooding in tributaries, a 1D model was applied to the river/tributaries and linked to the 2D model for the estuary and ocean, with the system running as a single model (in 1D/2D coupling). The results showed that the 1D/2D hydrodynamic coupling is robust and the D-Flow FM model is accurate for water level and tidal predictions for all three storms, both inside the bay and tributaries. Water levels were generally accurate, including the phases and peak magnitude and time. Simulation results using HEC-RAS showed that peak water level predictions were generally in good agreement with observations; however, time of occurrence was not well predicted. Therefore, we may conclude that D-Flow FM is the better model.
A sensitivity analysis has been carried out to evaluate the effects of stream discharge and atmospheric forcing on the model performance. The analysis revealed that hydrodynamic predictions are dependent on the stream discharge, especially in the upstream regions. In addition, the numerical results confirmed that the wind velocities and air pressure are very important for predicting accurate water levels. Better quality wind products can further improve the numerical estimations of water level.
The findings of this work show that using NWM/D-Flow FM/ADCIRC/WW3 modeling system has advantages over 1D river hydraulic models, particularly during the storm events. Simulations show that water level predictions depend on a precise representation of both river inflows and elevated sea levels. These results indicate the strengths of the proposed modeling approach for simulations of flooding during the most extreme storms, including for use in future operational forecasting systems. These models can be implemented for a river, estuary, or a coastal zone to achieve accurate water level from the combined storm tide and

10.1029/2019JC015822
Journal of Geophysical Research: Oceans freshwater inputs in both 1D and 2D. These models have parallel capabilities that significantly reduce computational time. In the future, we will pursue the fully two-way coupled NWM/D-Flow FM/ADCIRC/WW3 model that we are currently developing.
Although a 3D modeling is not necessary for flood forecasts, literature highlights that the lateral circulation, solved only by a 3D model, affects the estuarine water exchange and the resulting estuarine outflow (e.g., Du et al., 2018;Xiao et al., 2019). The classification of the vertical stratification in the transitional zone using the Fischer's flow ratio or the gradient Richardson number working with 3D setting could be beneficial. As a part of the broader inland-coastal coupling and modeling capabilities, we are working with ROMS (Shchepetkin & McWilliams, 2005), FVCOM (Qi et al., 2009), and Semi-Implicit Cross-scale Hydroscience Integrated System Model (SCHISM) (Zhang et al., 2016) communities for extending this framework to couple NWM with 3D coastal ocean models. We are also working with ADCIRC developers to include baroclinicity effects in our inundation calculation by vertically integrating water density from global ocean models and include it as an extra forcing term.