Simulating flood events at the Twin Ports of Duluth-Superior using a linked hydrologic-hydrodynamic framework

Generally, ports in the North American Great Lakes are not supported with navigational guidance (water level, water temperature, currents, ice) by NOAA’s Great Lakes Operational Forecast System (GLOFS). To examine the benefit of extending model coverage to this critical infrastructure, a linked hydrologic-hydrodynamic framework was developed for the Twin Ports of Duluth-Superior in western Lake Superior, and tested over three case studies of flooding due to storm surge and/or river flooding. Streamflow from 22 National Water Model (NWM) simulated and 3 gauged inflows within the domain was injected into a hydrodynamic model built on the Finite Volume Community Ocean Model (FVCOM), with a wet-dry grid that covered the harbor and surrounding floodplain. Model results from flood simulations compared well against time-series of water level and streamflow at gauges within the harbor, with root mean square error (RMSE) and bias values small relative to the typical fluctuations. Inclusion of NWM-simulated tributaries improved the accuracy of modeled water levels in the harbor, and increased simulated current speeds in navigational channels by up to 0.5 ms−1. Modeled inundation extent in the floodplain closely matched flood extent surveys conducted in response to a record flood event in 2012, demonstrating the capability of the modeling framework to provide flood guidance in the complex coastal setting.


Introduction
For thousands of years, the North American Great Lakes have been used as a commercial trade route. With over 3700 km of natural channels and man-made canals, the route connects the Atlantic Ocean through the St. Lawrence River to the largest freshwater system in the world. The route serves several regional and major international ports in the USA and Canada, including the high-traffic Twin Ports of Duluth-Superior. Located within the St. Louis River Estuary (SLRE) at the western end of Lake Superior, Duluth-Superior Harbor (DSH) is the farthest westward port in the Great Lakes for ocean bound vessels.
Having accurate, reliable, and timely short-term predictions of lake conditions (e.g., winds, waves, water levels, ice) are imperative for the safety and efficiency of vessels traveling the Great Lakes as well as emergency managers, coast guards, and port authorities. Since 2006, the National Ocean Service (NOS) has provided Great Lakes mariners with short-term forecast guidance on water levels, water temperatures, and currents using the National Oceanographic and Atmospheric Administration (NOAA) Great Lakes Operational Forecast System (GLOFS). The recently upgraded next-generation GLOFS  consists of four forecast modeling systems for Superior, Michigan-Huron, Erie, and Ontario, and is run four times per day to provide forecast guidance out to 120 h. GLOFS is forced by predictions from the National Weather Service (NWS) numerical weather prediction modeling systems and Responsible Editor :Joanna Staneva NOAA National Water Model CONUS Retrospective Dataset was accessed from https:// regis try. opend ata. aws/ nwm-archi ve.
1 3 NWS weather forecasts. In addition, GLOFS is forced by river discharge and water temperature observations from US and Canadian river gauges and discharge predictions from the NWS National Water Model at gauged locations.
GLOFS has two significant deficiencies. Firstly, the GLOFS model domain does not include ports and connecting channels and, therefore, does not provide forecast guidance for these small-scale features. For example, the Lake Superior Operational Forecast System (LSOFS) (Kelley et al. 2022) domain does not extend into the SLRE and does not provide guidance for DSH. Secondly, river inputs at ports are either not included or not resolved in GLOFS appropriately. This leaves a substantial gap in coverage between streamflow data and the lake model. These deficiencies can have serious impacts on the usefulness and accuracy of GLOFS's forecast guidance at significant ports, such as DSH, which are subject to flooding and fast-changing currents caused by heavy rain events, storm surge, or a combination of both. In recent years, a number of freight-carrying commercial shipping vessels have gone off-course while navigating DSH, leading to grounding or collisions with infrastructure (MPR News Staff 2014; Steil and Krueger 2018;Lovier 2020), demonstrating the need for navigational guidance in these areas.
Several modeling studies of the SLRE have previously been conducted. Smith (2020) developed a hydrodynamic model with the Environmental Fluid Dynamic Code (EFDC) that covered the SLRE and extended 7.5 m out into Lake Superior. The model was validated from April to November 2017 and found it was able to simulate water level, flow, and water temperatures on a daily time scale, but not at the hourly scale (Smith et al. 2020). Reschke and Hill (2020) applied the Finite Volume Community Ocean Model (FVCOM) for the SLRE and Lake Superior and focused on different scenarios of river discharge (low, moderate, high, and extreme flows) and Lake Superior water levels (low and high levels) for 2008 through 2018. The model's predictions of water temperatures and velocities at aquatic vegetation sampling sites for the different scenarios were used to refine aquatic habitat maps at sites in the SLRE. In addition, Anderson and Wu (2018) developed and implemented the Integrated Nowcast/Forecast Operation System for COastal, Riverine, and Estuarine environment (INFOS CORE) for Lake Superior and the SLRE. INFOS CORE is composed of SCHISM and WaveWatchIII models and provides experimental, real-time forecast guidance for Lake Superior and the SLRE out to 48 h.
This new study has three main goals: (1) Improve the hydrologic and hydrodynamic modeling framework to provide coverage for the SLRE and DSH; (2) include all streams and tributaries that flow through the SLRE to improve the water balance in this region; and (3) conduct a skill assessment during river and storm surge driven flood events. By addressing these goals, we seek to demonstrate the potential improvement to lake and flood forecasts for Great Lakes ports through enhancement to existing operational models.

Study area
The study area is the SLRE and DSH. The upper part of the SLRE has numerous shallow bays (e.g., Saint Louis Bay, Allouez Bay, Spirit Lake), wetlands, and forested areas along its shore. The lower part of the SLRE is highly industrialized with commercial and residential areas (cities of Duluth and Superior), culminating at DSH and discharging into Lake Superior. The major rivers flowing into the SLRE include the St. Louis River, the largest US tributary to Lake Superior, the Nemadji River, and the Pokegama River. The average yearly flow of the St. Louis River at the Scanlon (MN) Dam is 64.7 cms (2,284 cfs) with average yearly high and low flows of 414 cms (14,617 cfs) and 13.2 cms (465 cfs), respectively (Minnesota Pollution Control Agency 2022). The SLRE is affected by river flooding (e.g., Aug. 1972, and also back flows from Lake Superior during storm surges (e.g., Oct. 2018, Oct. 2019), and impactful seiche events (e.g., Nov. 27, 2001).
DSH is North America's farthest-inland freshwater seaport and consists of 79 km of harbor frontage. DSH is sheltered by a 14.5-km natural breakwater/barrier beach from the Duluth entrance channel to the Superior entrance channel (Fig. 1), the longest fresh-water barrier beach in the world (Peek 1914). There are several shipping channels in DSH (North, South, Upper, Minnesota, Superior Front, Allouez Bay) and harbor basins (Duluth, East Gate, and Superior). An estimated 35 million tons of cargo (iron ore, coal, grain, cement) move through DSH each year (Duluth Seaway Port Authority 2022). Shipping channels within DSH are maintained through dredging, with the entrance channels between the Lake Superior and DSH approximately 9-m deep and channels in the upstream portions of DSH approximately 7 m deep. The maximum depth within DSH is approximately 12 m.

Model setup
The numerical lake model used by GLOFS is the FVCOM. FVCOM is a prognostic, unstructured-grid, finite-volume, free-surface, three-dimensional (3D) primitive equation coastal ocean circulation prediction model developed by the researchers at the UMASS-Dartmouth and Woods Hole Oceanographic Institution (Chen et al. 2003(Chen et al. , 2013 and adapted for freshwater application to the Great Lakes (Anderson et al. 2010). The model consists of momentum, continuity, water temperature, salinity, and density equations and is closed physically and mathematically using turbulence closure sub-models. The horizontal grid is composed of unstructured triangular cells with a generalized terrain-following vertical coordinate system. The Mellor Yamada 2.5 turbulence closure scheme (TCS) was used for the vertical and the Smagorinsky TCS was utilized for the horizontal turbulence parameterization. FVCOM is solved numerically by a second-order-accurate discrete flux calculation in the integral form of the governing equations over an unstructured triangular grid. The 3D model solution is determined using a mode-splitting technique by which a two-dimensional external mode is updated at frequent intervals, while the more slowly evolving internal mode is obtained less frequently. The free surface, defined as the external mode, is integrated by solving vertically averaged equations with a smaller time step, while the 3D momentum and tracer equations, defined as the internal mode, are integrated with a larger time step. Following every internal time step, an adjustment is made to maintain numerical consistency between the modes (Chen et al. 2013). Simulations for this study were carried out using FVCOM version 4.3.1 in a 3D baroclinic mode, with the Coupled Ocean-Atmosphere Response Experiment (COARE) heat flux algorithm Version 2.6 (Fairall et al. 1996(Fairall et al. , 2003, and with dynamic wetting and drying for floodplain areas. The NOAA National Water Model (NWM) is a hydrologic modeling framework (Gochis et al. 2018) that simulates different physical processes of the water cycle in order to provide high-resolution forecasted streamflow across the continental United States (NWS 2022a; NWS 2022b). Simulated streamflow data was obtained from a 42-year retrospective run using the latest 2.1 version (NOAA National Water Model CONUS Retrospective Dataset n.d.). The retrospective run used meteorological forcing data from the National Weather Service (NWS) Analysis of Record for Calibration (AORC), which incorporates several near-surface meteorological products such as the North American Regional Reanalysis (NARR) data, as well as radar data and precipitation observations (NWS 2021). It is important to note that while the real-time NWM products assimilates streamflow gauge observations, the retrospective model run does not.
In order to link the hydrodynamic and hydrologic models, a high-resolution grid was created to cover the SLRE and floodplain from Lake Superior to the Fond du Lac Dam (Fig. 1). The LSOFS operational grid (purple) only extends to the barrier beach at the western end of Lake Superior and ranges from 200 m near the coastlines to 2.5 km in the open water. The extended grid (red) has a much higher resolution Orange dots represent the three locations of the USGS sites with observational streamflow and the green dot represents the NOAA CO-OP site with observational water levels. The light blue lines represent the streams and tributaries from the operational version of the NWM that feed into the estuary, and the bold blue lines represent the 24 rivers whose streamflow was used in the linked model framework. Red outlines the FVCOM extended grid developed and used for this project. The purple outline is the extent of the operational LSOFS grid. The black rectangles on the map capture the locations of the two inserts showing the high-resolution extended grid across the Duluth entrance channel (top) and the Superior entrance channel (bottom) 1 3 of 5-10 m inside the SLRE and DSH and a lower resolution of 10 km in the open water of Lake Superior. This configuration directly simulates the influence of Lake Superior on the SLRE by including a coarse representation of the lake in the domain, while focusing computational resources on the high-resolution harbor-estuary portion of the domain that is the focus of this study.
There are 24 streams that flow in and through the SLRE and one outflow stream at the far eastern point of Lake Superior, the St. Mary's River. Simulated streamflow data for all 25 streams was extracted from the NWM 2.1 Retrospective Dataset. The NWM simulated streamflow for the St. Louis River, Nemadji River, and St. Mary's River were compared to observational data from USGS gauges (USGS Water Data for the Nation n.d.). The simulated flows were grossly underestimated which is why the operational NWM assimilates observations into their streamflow forecasts. To better represent the streamflows in this experiment, we replaced the simulated flows at those three rivers with the observations. For the remaining 22 tributaries, simulated streamflow from the NWM retrospective run was used. All streams were applied in the lake model at the nearest FVCOM node to their endpoint (see Table S1). For larger rivers with channels that span multiple nodes in the model grid, flows were evenly distributed across all nodes comprising the river.
Topobathymetric data used to develop the model grid is based on the National Centers for Environmental Information's Lake Superior bathymetry dataset for the Lake Superior, nautical charts from the NOS Office of Coast Survey for the SLRE, and was obtained from the NOS Office for Coastal Management for floodplain regions. Figure 2 shows the bathymetry used in the model simulations. The bathymetry was not smoothed when developing the model grid. While the nautical charts are believed to provide the bestavailable representation of the SLRE bathymetry, there is some uncertainty in the depth information due to the age of the data and the dynamic nature of the domain. Processes such as scouring, sedimentation, and dredging are likely to alter the bathymetry over time, particularly at critical locations such as the harbor entrances, and such changes would not be accounted for in the bathymetry.
A 1-month spin-up period was used for each of the case studies using all 24 inflows and 1 outflow. Simulated data was output hourly, with the exception of water levels which output 1-min probe data for direct comparisons to 6-min observational data during the skill assessment.
Water level tracking was used to prevent drift based on a 3-day moving average comparison between observed and modeled water levels at 4 water level stations spaced roughly equidistant across the south shore of the lake, which was calculated on a daily basis. Extra precipitation or evaporation was evenly applied over the next 24-h run to counter the difference between observed and modeled water level.
One-second internal and external time steps were used for all of the model runs except for calendar days where the St. Louis River exceeded 800 cms. During those days, 0.2-s internal and external time steps were used. It was found that the shorter 0.2-s time step was necessary to satisfy the Courant-Friedrichs-Lewy (CFL) condition at the inflow location of the St. Louis River during high flow events, but that the Fig. 2 The bathymetry [m] used in the FVCOM case study experiments. Negative values represent land. Color scale is saturated at 12 m in Lake Superior for better representation of smaller depths in DSH 1-s time step could otherwise be used to improve computational efficiency. Running on 224 processors, computational efficiency was 3-h wall time per model day when using the 1-s time step and 12-h wall time per model day when using the 0.2-s time step.

Case studies
This study focuses on two record breaking flood events in DSH: the June 2012 heavy rain event, the October 2019 storm surge event, as well as a third case combining the first two to create a combined hydrologic-storm surge, or "worst case scenario" (WCS) flood event. These flooding events were chosen to test the model framework because of their different meteorological causes (rainfall versus wind and pressure).
The first case study was based on a flood event caused by heavy rainfall in June 2012. The week prior to the record event, several weather systems moved across Minnesota, priming the region with several inches of rain and saturating the soil. From June 19 to 20, the official rainfall total at the Duluth International Airport weather station was 18.4 cm (7.25 in) with locally higher amounts along the northern Lake Superior coast (NWS 2022c). During those 48 h of heavy rainfall, the water level rose by 3.35 m (11 ft) at the Scanlon Dam on the St. Louis River, breaking the previous record from 1950. The heavy rains along with steep terrain left the region with significant damage, which included roads being washed away and the death of several zoo animals at the Lake Superior Zoo.
The second case study focuses on the flooding event in October 2019. This event was different from the 2012 event in that the flooding was mainly driven by a rapid drop in air pressure, strong winds, and storm surge. On October 21, the meteorological station at the Blatnik Bridge recorded wind gusts up to 61.7 knots (31.7 ms −1 or 71 mph). The water level at the NOS DSH gauge reached 184.3 m (604.75 ft), breaking the previous 184.2 m (604.42 ft) record set in 1985 (NWS 2022d).
Lastly, a third case study was developed to look at a hypothetical combined hydrologic-storm surge WCS that combines the strong streamflow caused by heavy precipitation from the 2012 event with the storm surge and high winds and air pressure gradient of the 2019 event. This scenario was designed to look at the future flood possibilities, surface water currents, and how useful accurate flood guidance can be in this region.

Model forcing
For both the 2012 flood event and the 2019 flood event, two different simulations were run with different configurations but the same linked framework ( Table 1). The "Base" model run is considered the baseline simulation and includes the hybrid streamflow dataset consisting of NWM simulated streamflow from the 22 tributaries as well as observational streamflow data from the three gauged streams as described in the previous section. The other model run, the "USGS Rivers Only," removed the NWM simulated streamflow and only included observations from the three USGS river gauges.
The 2012 simulations used meteorological forcing from interpolated observations from the NOAA Great Lakes Environmental Research Laboratory (GLERL) Marine Observations Network (MAROBS) archive. Precipitation data was from the NWS/National Centers for Environmental Prediction's Environmental Modeling Center (NCEP/EMC) 4 km gridded product known as Stage IV data, which is based on multi-sensors and manual quality control (Du 2022).

3
For the 2019 event, NOAA High-Resolution Rapid Refresh (HRRR) (Dowell et al. 2022;James et al. 2022;NWS 2020) data was used for meteorological forcing including precipitation. The HRRR output was available hourly at 3-km spatial resolution. Figure 4 shows HRRR wind speeds and air pressure over the region during the time of peak water level at DSH (October 21, 2019, 2100 UTC). A spinup period from September 1, 2019, to October 18, 2019, was completed using the base simulation configuration. Both the Base and USGS Rivers Only simulations were initiated from the same spin up and were run from October 18 to 25, 2019.
The WCS simulation was configured by shifting the timing of forcing inputs such that the hydrologically driven water level peak from the June 2012 event coincides with the wind-and pressure-driven water level peak from the October 2019 event. Because the initial peak in 2012 water level was influenced by the movement of the low-pressure system across the domain, an interim simulation was conducted for the June 2012 event that was forced with a constant air pressure of 1000 hPa. This run removes the influence of pressure on water levels and isolates the hydrologically driven component of the water level rise. After reviewing results from these scenarios, river and precipitation forcing from the 2012 event were shifted forward in time by 2678 days and 12 h, calculated as the time difference between the peak water level from the 2012 constant air pressure simulation (June 21, 2012 0900 UTC) and the peak water level from the 2019 base simulation (October 21, 2019 2100 UTC). This results in a WCS configuration that is forced with the wind and pressure data from October 2019 coincident with the inflow and precipitation data of June 2012. Additionally, in shifting the 2012 forcing forward to align with the 2019 storm surge, the scenario is run under the high lake level conditions of 2019, thus exacerbating the flood severity. Figure 5 shows a time-series of air pressure, accumulative precipitation, surface wind speed and direction, and water level at DSH during the June 2012 simulation period. The weather system begins to move through the region on June 18 as seen by the air pressure and wind shift (Fig. 5a, c). The heaviest rainfall occurs on June 20 (Fig. 5b) which is when the surge in water level is observed at the harbor. Both of the model simulations do well at capturing the short-period oscillations in water levels prior to the event as well as the increase on June 20 (Fig. 5d). The model has an overall root mean square error (RMSE) and bias of 6.4 cm and − 0.4 cm, respectively, at the NOS DSH gauge 9,099,064 (Table 2), which is small compared to sub-daily water level fluctuations of approximately 20 cm observed at the location (Fig. 5d). The low bias in water level increases in magnitude from − 0.4 to − 1.7 cm when only observations from USGS where the black line is the 6-min observations, blue line is 1-min output subset to 6-min data from the Base simulation, and red line is 1-min output subset to 6-min data from the simulation with only USGS gauged river inflows A USGS-published scientific report (Czuba et al. 2012) used surveys of high-water marks combined with high-resolution elevation maps to create flood extent and inundation maps for six regions that were highly affected by the June 2012 flood. One of these areas, the Fond du Lac neighborhood of Duluth, is inside the study area, providing information on flood extent that can be used to validate modeled flooding in the wet-dry floodplain. The river morphology within the Fond du Lac region is characterized by several small islands and a narrow peninsula, with houses, parks, a highway, and a boat launch along the northern river bank, and with the residential neighborhood extending out onto the peninsula in the SLRE.

Hydrologic event
Modeled inundation generated through dynamic wetting and drying in FVCOM from the two simulations were compared to the inundation survey published in the USGS report, with regard to both maximum flood extent (Fig. 6) and maximum inundation depth (Fig. 7). Impacts from flood water encompassed all islands and the peninsula within the SLRE, and extended to portions of the northern river bank. Both model simulations agree closely with the survey data, demonstrating the capability of the FVCOM model to provide accurate guidance of flood extent in floodplain areas. The largest deviations between modeled and observed inundation occur near the coast in the Base run at locations where the NWM streams enter the hydrodynamic model. The model can produce localized unrealistic flooding in areas where smaller NWM inflows have poorly-resolved thalwegs (deepest channel in a river) in the model grid. This can occur where the width of the inflow channel is narrower than the resolution of the model grid, or where the channel passes underneath roads or bridges, and hydrologic connections (e.g., culverts) are obscured in the floodplain topography.
Spatial analysis of the SLRE and DSH showed differences in water levels and surface currents between the simulations (Fig. 8). Maximum modeled water levels were highest in the Fond du Lac region, which saw significant flooding during the event. The additional inflows in the base simulation increased the water volume during the 2012 flood event causing a 10-20-cm increase in water level in the SLRE and 5 cm in the DSH (Fig. 8c). Circulation in the SLRE and DSH is strongly impacted by Lake Superior's 8-h period seiches leading to periodic flow reversal in both entrance channels between the lake and DSH. Maximum currents occur in channels and also along the thalweg up to the most upstream location. Before the June 20 flood event, maximum surface current speed reached ~ 1 ms −1 in both entrance  Czuba et al. (2012). The red line represents the flood extent from the base simula-tion, purple represents the USGS Rivers Only simulation, dark blue lines represent the rivers as used from the NWM, and the gray lines show the natural shoreline as observed when the region is not flooded channels and ~ 0.5 ms −1 elsewhere within DSH. During the flood event, maximum surface current speeds occurred far upstream, reaching ~ 3 ms −1 just below the most upstream location, and increased up to ~ 1.5 ms −1 in the DSH entrance channels. Most of the SLRE flow exits DSH through the Duluth entrance channel, while Nemadji River flow exits primarily through the Superior entrance channel. Additional flow from smaller tributaries in the NWM led to an increase in maximum surface current speed up to 0.5 ms −1 in entrance channels and up to 0.25 ms −1 elsewhere in DSH (Fig. 8f).
Spectral analysis was conducted on water level timeseries to quantify the model's representation of the Lake Superior seiche signal (Fig. 9). Spectral analysis was performed using the Welch (1967) time-averaging method on 1-min water level output from the model at the NOS DSH water level gauge 9,099,064 (Fig. 1), with 72-h segments and a 24-h overlap. A dominant 8-h period is shown in both the gauge observations and the model output, which is consistent with the known 7.9-h period of the Lake Superior seiche (Jordan et al. 1981;Sorensen et al. 2004). The ability of the linked model to reproduce this key component of the lake-estuary circulation illustrates the need to incorporate the SLRE, including streamflow, into a linked hydrologic-hydrodynamic domain when simulating water levels within DSH. Although Lake Superior is downstream of the SLRE, it has a critical influence on estuary dynamics that is not represented in a model that is based on streamflow alone. The modeled energy at the 8-h period is higher in the model than the observations, indicating that the model is overestimating the seichedriven water level fluctuations in the DSH. In preliminary analyses, it was found that exchange between the lake and the harbor is highly sensitive to model bathymetry within the entrance channels. Given that the influence of the Lake Superior seiche signal at the water level gauge is controlled by exchange through the channels, any inaccuracies in the modeled channel depths would negatively impact the accuracy of modeled seiche-driven water level fluctuations in the harbor.

Storm surge event
The October 2019 event was different from the June 2012 event in that it was triggered by a strong low-pressure system that brought high winds to the region. DSH saw a rapid decrease in air pressure of 26 hPa from October 21 0300 UTC to 2300 UTC (Fig. 10a). The Duluth International Airport recorded several centimeters of precipitation on October 21, 2019; however, when total precipitation from HRRR forcings was analyzed, little to no rainfall was simulated at DSH during this period. Wind speeds reached 20 ms −1 (Fig. 10b), with even higher wind gusts, over DSH on October 21, 2019. These high winds combined with the strong air pressure gradient across Lake Superior induced a storm surge along coastal regions of western Lake Superior that impacted DSH. Fluctuations of streamflow in and out of DSH (Fig. 10c) and water levels within DSH (Fig. 10d) are dominated by the Lake Superior seiche for most of the simulation period. However, the storm surge produced an acute  influx of water from Lake Superior to DSH during the latter half of October 21 and caused water levels to rise by 30 cm in less than 24 h. This was followed by a period of high flow from DSH back into Lake Superior, which continued into October 22, as water levels in DSH receded to pre-surge levels.
While both the Base and USGS Rivers Only model simulations have a slight delay in the timing of the water level increase, they capture the peak water level in DSH on October 21 with a slightly higher peak water level in the Base run (blue) than the run with only USGS rivers (red). The 8-h seiche fluctuations in water level are also well-captured in the model simulation, both in the model time-series (Fig. 10d) and spectral analysis (Fig. 9b). The modeled water level in the Base run has an overall RMSE and bias of 4.6 cm and − 0.2 cm, respectively, relative to observations ( Table 2). The low bias increases in magnitude from − 0.2 to − 0.9 cm when only USGS-gauged rivers are included in the model, which demonstrates that incorporating smaller tributaries in the domain benefits water level simulations even under events for which river inflows are not the primary driver of the water level surge. The higher skill in water level results from the 2019 case, in comparison to the 2012 case, may be due to the use of more sophisticated HRRR forcing, which was not available for the time period used in the 2012 case, or it may be due to the surge-driven nature of the 2019 water level rise, which is a result of fundamentally different processes than the hydrologically driven 2012 flood.
Because the 2019 flood event was dominated by wind, not precipitation, the differences in maximum water level and maximum surface water currents between the two simulations are smaller than during the 2012 flood (Fig. 11). The additional inflows in the Base simulation raised the overall water level in the SLRE and DSH by only 5 cm (Fig. 11c). The strongest surface currents (up to ~ 0.6 ms −1 ) were again along the thalweg up to the most upstream location where they reached ~ 1 ms −1 (Fig. 11d, e). Very strong coastal flow is also seen just outside of DSH, east of the Superior entrance channel, with surface current speeds reaching 1 ms −1 . The model captures the full cycle of the 8-h seiche but mute the high-frequency oscillations, most likely due to the low-frequency of the simulated streamflow output (Fig. 10c). Seiche-related oscillations are still present on October 21, but are masked at the peak of the event by strong outflow which lasts in both observations and model runs for about 6 h, i.e., close to the 8-h seiche (Fig. 10c). Currents  [cms, lakeward] where the black line represents the 6-min observations (464646092052900), blue line is hourly output from the Base simulation, and red line is hourly output from the simulation with only USGS-gauged river inflows, d water level [m] where the black line is the 6-min observations, blue line is 1-min output subset to 6-min data from the Base simulation, and red line is 1-min output subset to 6-min data from the simulation with only USGS-gauged river inflows 1 3 were somewhat stronger in the Base simulation but the overall pattern in the difference plot is patchy; in some areas the Base simulation's currents are stronger, in others the USGS Rivers Only simulation's currents are stronger (Fig. 11f). At most locations, the difference in maximum surface current speed is less than 0.25 ms −1 .

Combined storm surge-hydrologic event
Time-series of the combined storm surge-hydrologic WCS event (Fig. 12) show simulated WCS water level (green) at DSH during the 2019 flood event compared to the water level from the 2019 Base simulation (blue). By combining the strong inflows from the 2012 event with the low-pressure system of the 2019, DSH sees a significant increase in the water level. The peak of the high-water level was seen at 2200 UTC on October 21 and the WCS saw a 0.11-m-higher simulated water level than the Base simulation of the 2019 Storm Surge Case.
Water levels and surface current patterns in the WCS (Fig. 13) generally resemble that of 2012 (e.g., strong currents along the thalweg and in both Duluth and Superior entrance channels) with additional features originating from 2019 event (e.g., strong coastal currents outside of the Superior entrance channel). Additional streamf low provided by the 2012 event produced higher water levels (Fig. 13a) and stronger current velocities (especially along the thalweg and in entrance channels) with surface currents increasing up to 1 ms −1 (Fig. 13f). At the same time, there are areas in DSH where the 2019 Base simulated currents that are stronger, although the difference is not as dramatic (only up to 0.25 ms −1 ).

Summary and conclusions
The FVCOM grid of the LSOFS was extended to cover the Twin Ports of Duluth and Superior including its harbor entrances and shipping channel, the SLRE and adjacent floodplain. The spatial resolution of the grid inside the SLRE and DSH ranged from 5 to 10 m. The NWM endpoints were associated with the nearest FVCOM grid node or multiple nodes for large rivers. The linked hydrodynamic and hydrologic modeling system was run for three scenarios: hydrologic/river flooding (June 2012), storm surge (October 2019), and combined hydrologic-storm surge event. For each scenario, two simulations were run: the Base run which included NWM simulated streamflow to 22 tributaries and observed streamflow at three USGS gauges and USGS Rivers Only which only used observed streamflow at the three gauges.
The model simulations during the hydrologic event demonstrated that the linked hydrodynamic-hydrologic modeling systems captured the short-period oscillations in the water levels before and during the event. They also indicated that incorporating streamflow from smaller tributaries improved the water budget in DSH. The simulations depicted the inundation well, both maximum flood extent and depth. However, the model produced localized unrealistic flooding in areas where smaller NWM inflows have poorly defined thalwegs in the model grid. This indicated that while the inclusion of NWM tributaries improved estimates of water level within the harbor, simulations of flood extent would further benefit from advanced spatial treatment of rivers in the model grid.
The simulations for the storm surge event accurately predicted the influx of water from the lake into DSH, which caused water levels to rise rapidly in the harbor and then the high outflow from DSH to Lake Superior. The simulations demonstrated that incorporating streamflow from smaller tributaries in the model domain during a non-river flooding event improved the estimation of water levels in DSH. While the model provides skillful simulations of water level and flow fluctuations between the Lake Superior and DSH, results may further be improved through improved representation of the bathymetry, particularly around the harbor entrance channels. Hydrodynamic modeling at the entrance channels is particularly challenging, due to the complex and dynamic nature of the lake-estuary system, and a bathymetry that is constantly evolving as a result of scouring, sedimentation, and active dredging.
The simulations for a hypothetical combined storm surgehydrologic event which combined the high river inflows from the June 2012 flooding event with the low atmospheric pressure and high winds of the October 2019 event resulted in water levels being significantly higher than water levels caused by high winds alone. This implies that higher water levels could occur in DSH by the simultaneous occurrence of both river flooding and storm surge when an extra-tropical cyclone passes through the western Great Lakes region, especially when overall lake-wide water levels are higher than normal.
Simulations showed high surface currents within the shipping channels linking DSH to Lake Superior. In particular, surface currents reached 1.5 ms −1 during the hydrologically driven 2012 flood event, due in part to additional flow from the NWM tributaries that were incorporated into the model configuration as a part of this study. The additional flow increased surface currents by approximately 0.5 ms −1 during the flood event, compared to the run with only USGS-gauged rivers. This strong influence of tributary flows on currents within the DSH entrance channels is consistent with past assessments of DSH hydrodynamics. The ability to accurately simulate and predict such changes in currents within DSH is important toward improving the safety of commercial shipping vessels and other watercraft. Grounding of commercial shipping vessels and collisions with the infrastructure illustrate the need for improved navigational guidance in Great Lakes ports and harbors, which are often characterized by a high volume of shipping traffic and fast-changing hydrodynamics, and which are currently a blind spot in GLOFS forecast guidance. This study demonstrates that hydrodynamic models in the Great Lakes can be extended to produce accurate simulations of dynamics within these coastal features, and is a step toward filling this critical gap.
Funding This project was funded by the NOS Water Initiative. Funding was awarded to the Cooperative Institute for Great Lakes Research (CIGLR) through the NOAA Cooperative Agreement with the University of Michigan (NA17OAR4320152). This is CIGLR contribution number 1216 and NOAA GLERL contribution no. 2021.

Data availability
The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.

Competing interests The authors declare no competing interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.