{'fraction_non_alphanumeric': 0.0387186410464781, 'fraction_numerical': 0.024693947099378718, 'mean_word_length': 4.504761904761905, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 12, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 11, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The Study of the Hurricane-Induced Storm Surge and Bay-Ocean Exchange Using a Nesting Model

Lagoon systems are more heavily impacted by hurricanes, whereas the relevant storm surge modeling studies have been paid little attention to lagoon systems and the storm-induced exchange in lagoon systems is even less understood. To address this gap, a three-dimensional unstructured grid-based model was configured for the Maryland Coastal Bays, a typical lagoon system with two unique inlets (Ocean City Inlet (OCI) and Chincoteague Inlet (CI)), to investigate how Hurricane Sandy impacted inlet dynamics. A nesting model framework was applied to provide the necessary remote forcing from a large model domain and maintain the intricate shoreline and bathymetry of an inner model domain. Results indicated that the flux patterns varied in response to the change in wind direction and rising/falling high water levels from the coastal ocean, rather than a single flow pattern during the passage of Sandy. From October 29 05:00 to 17:00 UTC, mild (> 10 m/s) and strong (> 15 m/s) northerly winds accompanied by the rising high water level from the coastal ocean promoted a mean inflow pattern at the OCI and a mean outflow pattern at the CI. Strong southwesterly winds (> 15 m/s) dominated in the bays from October 30 03:00 to 15:00 UTC. Under strong southwesterly winds and falling high water levels from the coastal ocean, flux was transported landward at the CI and seaward at the OCI. Sensitivity experiments on various storm temporal scales showed that a net inflow pattern occurred in the bays, and the net exchange amounts became smaller in response to longer storm durations. Residual effect of relatively high river flow from Sandy could still influence the salinity at the OCI, whereas the CI salinity was not affected by river flow owing to a long distance between the CI and river locations.


Introduction
A storm surge is an abnormal water level change generated by a storm, which is closely associated with extreme atmospheric pressure and wind forcing. The surge can be potentially devastating along the coasts to properties and the safety of human life. Hence, a comprehensive understanding of storm surge dynamics is essential for the relevant policy makers to make better and more informed decisions, to conserve environmental economic and natural resource subject to storm surges, and to minimize coastal hazards as much as possible (Psuty and Ofiara 2002;Lewis and Babson 2018).
Hurricane-induced storm surge modeling has been a subject of continuous research interest, and enlarging the model domain is a critical issue for simulating storm surges (Morey et al. 2006;Xia et al. 2008;Orton et al. 2012). However, less effort on storm surge modeling has been paid to lagoon systems along the coasts (Beudin et al. 2017). A lack of relevant studies constrains our comprehensive understanding of the storm dynamics in lagoon systems. A grid refine technique, a model nesting method, can be applied to cover the necessary remote storm forcing from a large model domain (LMD) with a coarse spatial resolution. This model nesting method maintains the accurate shoreline and bathymetry characteristics of lagoon systems by a smaller inner model domain (IMD) for the study area with a fine spatial resolution. The open boundaries of the IMD are nested within the LMD, and the open boundary forcing of the IMD is extracted from the LMD. The model nesting approach could ensure the same hydrodynamic information at the nested boundary interface between the two model domains (Morey et al. 2006;Brown et al. 2010).
Maryland Coastal Bays (MCBs), a typical lagoon system (~453 km 2 ) with two unique inlets (Ocean City Inlet (OCI) in Maryland and Chincoteague Inlet (CI) in the south Chincoteague Bay in Virginia), is a restricted shallow coastal lagoon system (~1 m) located behind the barrier island of the Delmarva Peninsula (Figs. 1 and 2; Boynton et al. 1996;Wazniak et al. 2005). The bays include five sub-bays and one major river: Assawoman Bay, Isle of Wight Bay, Sinepuxent Bay, Newport Bay, Chincoteague Bay that occupy the largest portion of surface areas of the MCBs, and St. Martin River (Fig. 1b). Two inlets display different inlet orientations: the OCI follows a northwest to southeast orientation, whereas the CI represents a northeast to southwest orientation ( Fig. 2b and c). Compared to the CI with a width of 2000 m and various depths (maximum of 7 m), the OCI is narrower (155 m wide) and has a relatively flat bottom topography ( Fig. 2b and c). The bays can be divided into two regions (north and south MCBs: NMCBs and SMCBs) based on the location of the OCI (Duan et al. 2015). Tidal flushing is known to have a significant influence in the vicinity of the inlets (Boynton et al. 1996). Winds aligning with the shape of Chincoteague Bay have a great influence on the inlet exchange dynamics (Kang et al. 2017). The circulation in the interior of the bays is weak due to the limited freshwater discharge and constricted inlets (Lung 1994;Pritchard 1960;Kang et al. 2017). Relatively high salinity (> 25) can be observed in most parts of the bays except near creeks where freshwater flow dominates (Wang et al. 2013), and the bays are vulnerable to hurricanes because of their locations.
In general, the combined physically interactive forces, such as winds, tides, and river discharge, are the primary factors in governing the exchange dynamics of most estuarine systems (Xia et al. 2011;Umgiesser et al. 2014). Despite previous studies that illuminated the exchange dynamics in response to external forcing agents (i.e., winds, tides and river flow) on multiple temporal scales (Kang et al. 2017;Valle-Levinson et al. 2018;Liu et al. 2018), little is known about the storm-induced exchange dynamics in the MCBs as well as other lagoon systems. Although Li et al. (2009) documented the storm-induced saltwater flux through multiple inlets in Lake Pontchartrain estuary based on observations; the driving mechanisms of exchange dynamics still remain unclear. It is believed that storm durations impact the estuarine exchange dynamics, whereas the effect of storm duration on exchange dynamics in lagoon systems still remains unknown. Meanwhile, salinity dropped in response to the high river flow from hurricanes (Cho et al. 2012;Brown et al. 2014). Corresponding low salinity could alter the surrounding ecosystem, such as causing high mortalities for oysters (Munroe et al. 2013). Especially, lagoon systems are strongly influenced by relatively small, but intense river flow from hurricanes. Given that salinity is an important indicator for water quality and ecology, it is crucial to understand the salinity behaviors before and after a hurricane in lagoon systems. However, at present, the effects of the relatively high river flow during a hurricane on the salinity in the lagoon systems are not fully understood.
Hurricane Sandy was one of the most destructive hurricanes in US history. It developed in the Caribbean Sea on October 22, 2012, passed through the Mid-Atlantic Bight and made landfall as a Category 1 hurricane on the New Jersey coast on October 29 (Fig. 3a). Sandy brought huge impacts to the MCBs. According to the observations, the National Ocean Service gauge at the OCI measured a storm surge of 1.32 m with a maximum observed wind gust of  (Krantz et al. 2009). b Inlets, subbays, a river name and important creek/stream locations of the MCBs approximately 50 knots (Blake et al. 2013). However, the storm-induced flux and salinity behaviors in the MCBs during the passage of Sandy still remain to be elucidated.
To understand the exchange dynamics of a lagoon system (i.e., the MCBs) in response to a passing storm (Sandy), a state-of-the-art, finite-volume, and three-dimensional (3D) hydrodynamic model with the application of the model nesting method was applied. The study focuses on the following: (1) exploring the influence of the model domain size on storm surge simulations, (2) investigating the effect of strong hurricane winds and low atmospheric pressure on the hydrodynamics of the MCBs, (3) quantifying the flux in response to Sandy and various storm durations, (4) shedding light on the hurricane-induced salinity variations in the bays.

Model Description
The Finite Volume Community Ocean Model based on the MCBs (Kang et al. 2017) was used to simulate the hydrodynamics in the MCBs along with the adjacent coastal Atlantic Ocean for both LMD and IMD. With the applications of unstructured horizontal grids and vertical sigma coordinates, the model could accurately follow the irregular coastline and complex bathymetry of the MCBs. The Smagorinsky turbulent closure parameterization was used to solve the horizontal turbulence mixing, and the vertical turbulence mixing was determined by the Mellor-Yamada level 2.5 turbulent closure model.

Unstructured Mesh
The LMD with the extension of the entire North American East Coast was created by Surface-Water Modeling System in order to cover the remote forcing for storm surge simulations. The LMD encompasses the entire North American East Coast Shelf, stretching from Key West, Florida, to Nova Scotia in Canada, and extends out from the North American East Coast about 1600 km into the Atlantic Ocean with 30,247 nodes, 57,264 elements, and 113 open boundary nodes (Fig. 3c). The grid spatial resolutions range from about 40 km at the open boundary to about 20 m in the MCBs. The coastline data were derived from National Oceanic and Atmospheric Administration (NOAA) National Geophysical Data Center. The bathymetric data was derived from the 1 arc-minute global relief model ETOPO1 Global Relief Model from NOAA. Vertically, a total of five uniformly vertical layers were used for the LMD. The external mode step of the model integration was 1 s, and the internal to external mode ratio was 6.
In FVCOM, the bottom stress is calculated by a quadratic formula in the form of where u b and v b are the bottom current velocity along the axis of x and y, and C d is the bottom drag coefficient that is The open boundary nodes of the IMD for the MCBs were nested within the LMD ( Fig. 3d and e). The model nesting approach could allow the open boundary forcing including water elevations and currents in the IMD extracted from the LMD (Morey et al. 2006). The IMD with 7332 nodes and 12,428 elements was the same one used in Kang et al. (2017), and the grid resolutions are in the range from 1.60 km at the open boundary to less than 15 m in the bays (Fig. 3e). Five uniformly vertical sigma layers with a 4-s time step were applied for model simulations when applying the IMD.

External Forcing
The hydrodynamic model based on the LMD was forced by surface meteorological forcing (air pressures and wind fields) and equilibrium tides at the open boundary. Surface meteorological forcing from two sources was selected for model comparisons: NARR and the combined data from NARR and NAM (NARR+NAM; Warner et al. 2017). The NARR data obtained from North American Regional Reanalysis were at a three-hourly time interval in a 32-km spatial resolution. The NAM data from North American Mesoscale Forecast system were obtained from Rutgers University, and they were also at a three-hourly time interval but in a higher spatial resolution of 12 km. The NAM data could not be sufficiently applied for the LMD because the spatial coverage of the LMD extends beyond the spatial coverage of NAM data. Thus, the NARR data were used to provide the information for the regions that NAM data could not cover (Warner et al. 2017). The predicted tidal inputs were interpolated from the Oregon State University global inverse tidal model TPXO 7.2 (Egbert and Erofeeva 2002). Eight harmonic tidal components including four semidiurnal constituents (M 2 , N 2 , S 2 , and K 2 ) and four diurnal constituents (O 1 , P 1 , K 1 , and Q 1 ) were imposed at the open boundary of the LMD with the consideration of the inverted barometer effect.
Model applications with the IMD were driven by riverine inputs, surface meteorological forcing (air pressure, wind fields, long/short wave radiation, air temperature, and relative humidity), and the open boundary forcing. Only three observed riverine stations with daily freshwater discharge were available in the MCBs from the U.S. Geology Survey (USGS) (Fig. 2a). Due to numerous small creeks and streams in the bays, besides the three USGS stations, the related river flow was added through Grays Creek, Bishopville Prong, Shingle Landing Prong, St. Martin River, Manklin Creek, Turville Creek, Herring Creek (from the NMCBs), Ayer Creek, Newport Creek, and Marshall Creek (from the SMCBs), with the mean of the three observed river flows supporting the salinity simulation (Wang et al. 2013;Fig. 1b: purple diamonds from north to south). Surface meteorological forcing was obtained from two sources (NARR and NAM) over the entire IMD. Open boundary forcing in the IMD was extracted from the LMD based on the model nesting method.
Meanwhile, accurate wind fields are important to storm surge simulations. Previous studies applied the method from Holland (1980) to reconstruct wind and atmospheric pressure fields for storm surge simulations (Xia et al. 2008;Weisberg and Zheng 2006;Rego and Li, 2010;Ma et al. 2015). Here, we applied the original data from NARR and NAM coupled with the method from Holland (1980) to compare the model performance in simulating water levels during Sandy. Wind speed and atmospheric pressure relative to a storm center (Holland 1980) are addressed as below: where r is the radial distance from the hurricane center; V w and P are the wind speed and atmospheric pressure as a function of r; ρ a is the air density; P n and P c are the ambient and minimum central atmospheric pressures, respectively; V m is the maximum wind speed; e is the natural logarithm base; and A and B are the storm scale parameters related by A = (R max ) B , where R max is the radius of maximum winds. B determines the shape of a storm wind field. The working best track of Sandy was obtained from National Hurricane Center with six-hourly data subject to the modification during the life cycle of the cyclone. The data include storm center locations, atmospheric pressures, maximum 10 m wind speeds, and the 64 knots wind radius in four directions (NE, NW, SE, SW). The maximum wind radius (R max ) was the averaged maximum wind radial distance in four directions using the 64 knots wind radius. Further linear interpolations were done to provide the data every 3 h. However, the method from Holland (1980) could only construct an axisymmetric wind field, neglecting asymmetry in a real hurricane. This is one uncertainty in Holland (1980). From October 24 12:00:00 UTC to October 29 18:00:00 UTC 2012, the reconstructed wind and atmospheric pressure fields replaced the original NARR and the combined NAM and NARR data and forced the hydrodynamic model when applying the LMD (Table 1: RE-NARR-LMD and RE-NARR+NAM-LMD). Further experiments were conducted to compare the skills of two wind sources in simulating water levels along the US East Coast during Sandy. In a similar manner, the reconstructed wind and atmospheric pressure fields also replaced the wind and atmospheric pressure fields from NARR and NAM over the entire IMD during Sandy. Considering the MCBs is relatively small, to improve the accuracy of wind fields for the IMD simulations, the hourly observational data including wind speeds, wind directions, and barometric pressures at the OCI from NOAA were used to correct the wind and barometric fields that were computed by the method of Holland (1980) (Fig. 2b; RE-NARR-OB vs RE-NAM-OB; Shen et al. 2006a).
F(x, y, t) is a variable that was reconstructed based on the method from Holland (1980) at node (x, y).F x; y; t ð Þ is the new variable after corrections, expressed aŝ is the correction factor for observed variables at the ith observational station. F obs is the observed variable. W i (x, y) is the weighted function corresponding to the ith observational station.
Wind stress in our storm surge model is computed by the following: where ρ a is air density; V w ! is wind speed at 10 m height; C dw , a drag coefficient dependent on the wind speed, is given by Large and Pond (1981): Observational Data, Model Evaluations, and Experimental Designs The model using the LMD started running from September 21, 2012. Results from October 1 to November 10, 2012 were used for the Sandy study because the model reached an approximate dynamic equilibrium after running for 10 days. Fifteen observed water level stations along the US East Coast could be gained from NOAA Center for Operational Oceanographic Products and Services ( Fig. 3a and b). Model simulations by the IMD started running on October 1 and ended on November 10, 2012. In the MCBs, the observed station at the OCI, which is operated by NOAA, provided the observed water level and wind data (Fig. 2b). A USGS storm tide sensor at the CI also provided the water level information during Sandy from October 27 to November 3, 2012 (Fig. 2c). Data from four continuous observational stations located near the creeks were obtained from the Maryland Department of Natural Resources (MD-DNR) to exhibit the spatial salinity and water temperature variations in the MCBs during Sandy (Fig. 2a). Due to a lack of observed water level inside the bays, the data from the pressure sensor of YSI at the four continuous stations were used to convert to the estimated water level inside the bays (Fig. 2a). Correlation coefficients (Rs), mean absolute errors (MAEs), and root mean square deviations (RMSDs) were applied to appraise the model performance.
To investigate the influence of the model domain size on the storm surge modeling, model simulations based on a single usage of the IMD (Table 1: RE-NARR-OB-IMD and RE-NAM-OB-IMD) were conducted for comparisons of model simulations based on the usages of both the LMD and IMD with a model nesting method (Table 1: RE-NARR-OB-NEST and RE-NAM-OB-NEST). Considering the inverted barometer effect of a hurricane, the inverse pressure adjusting method were used at the open boundary of the IMD to account partially for the remote meteorological forcing combining with the tidal forcing obtained from TPXO 7.2 for storm surge simulations (Shen et al. 2006b;Shen and Gong 2009). One hectopascal drop in air pressure results in approximately 0.01 m rise in water level (Danard et al. 2003). The adjusted water level inputs for the IMD forced the model applications with a single usage of the IMD (Table 1: RE-NARR-OB-IMD and RE-NAM-OB-IMD). Model performances in simulating water levels based on different domain sizes and winds are described in "Results." Table 2 lists the experimental designs to identify the contributions of each external forcing factor (storm winds, high water levels from the coastal ocean) to exchange dynamics in the MCBs during Sandy. All the model settings in Table 2 are based on case RE-NAM-OB-NEST in Table 1.

Water/Salt Flux Calculation
A quantitative assessment linking water exchange and salt flux was applied to the inlet area of the MCBs to constitute a basis of understanding their behaviors in response to a passing storm ( Fig. 2b and c). Flux calculations were based on the experimental designs in Table 2.
Flux calculations were described in detail by Xia et al. (2011) and Kang et al. (2017).
where v is the current velocity that is perpendicular to the direction of the inlet; h is the water depth; η is the water elevation; S is the salinity; x 1 and x 2 are the left end and right end horizontal locations of the inlet; t 1 and t 2 are the starting and ending time for the integral period; dx is the grid width across the inlet; dh is the depth along the vertical direction and dt is the time step.

Results
The behaviors of NARR and NAM + NARR coupled with the method from Holland (1980) in simulating water levels along the US East Coast were compared during Sandy from October 27 to November 3, 2012 (Table 1: RE-NARR-LMD and RE-NAM + NARR-LMD). Through the application of the LMD, the model simulations captured the baseline of the observed water level at all 15 stations along the US East Coast (Fig. 4ao). Compared to RE-NARR, RE-NAM + NARR showed a better skill in water level simulations with a higher correlation coefficient and lower RMSD and MAE (Fig. 4p1). As Sandy made landfall on the New Jersey coast on October 29, the peak water level occurred in the Mid-Atlantic Bight including the Virginia, Delaware, New Jersey, New York, Connecticut, and Rhode Island coasts between October 29 and 30 ( Fig. 4fo). However, our model results underestimated the water level at Duck, NC and Swell Point, VA under RE-NARR ( Fig. 4e and f) and this discrepancy was partially attributed to the errors in the wind forcing and weaker winds of RE-NARR (Fig. 5). For the Virginia, Maryland, and Delaware coasts, our study areas of interest, RE-NAM + NARR produced a better performance in simulating water levels with higher correlation coefficients and lower RMSDs and MAEs than RE-NARR because wind speeds of RE-NAM + NARR were higher and closer to observed wind speeds than RE-NARR ( Figs. 4p2 and 5d).
Overall, RE-NAM + NARR showed a more satisfactory skill in water level and wind speed simulations. To fulfill the objectives of this study, RE-NAM + NARR-LMD were applied to provide the nested open boundary forcing for the IMD for Sandy simulations.
Through a single application of the IMD considering the inverted barometer method (Table 1: RE-NARR-OB-IMD and RE-NAM-OB-IMD), two wind fields (RE-NARR-OB and RE-NAM-OB) showed a similar performance in the water level simulation of the MCBs (Fig. 6a-f). It can be seen that the small IMD model considering the inverted barometer effect greatly underestimated the water level at both inlets and the other four stations in the bays (Fig. 6a-f; Table 3). This suggests the small IMD cannot account for the sufficient remote forcing and also highlights the importance of a LMD for storm surge simulations.
With the remote forcing from case RE-NAM + NARR-LMD, RE-NARR-OB and RE-NAM-OB still generated a similar performance in the water level simulation (Fig. 6a-f; Table 3). Although bias between observed and modeled water levels still existed, the model nesting method greatly improved the model results compared to the single application of the IMD without remote forcing (Fig. 6a-f) (Table 3). The underestimations may have resulted from the open boundary forcing and the errors in the wind fields (Fig. 6g, h1, and h2). This model nesting method for storm surge simulations covers the necessary remote forcing simulated by the LMD and keeps the complex shoreline and bathymetry of the IMD, and thus, it enhances the simulation accuracy. Given the better water level performance through the model nesting method, the relevant surge, salinity, and temperature discussions were based on cases RE-NARR-OB-NEST and RE-NAM-OB-NEST.
The maximum storm surge level and lasting periods of "highest surge" displayed distinct spatial and temporal variability in the bays (Fig. 7a1, a2, b1, and b2). The highest surge duration is defined as the durations that the surge exceeds 0.7 m (Beudin et al. 2017). The higher surge occurred around Sinepuxent Bay (~1.6 m) and the NMCBs (0.8-1.2 m), and the lower surge (0.6-1 m) mainly concentrated in the southern Chincoteague Bay (Fig. 7a1 and a2). Meanwhile, the highest surge duration lasted longer in Sinepuxent Bay (RE-NARR-OB-NEST: 27 h; RE-NAM-OB-NEST:~27 h) ( Fig. 7b1 and b2). Shorter highest surge duration concentrated in the southern Chincoteague Bay (RE-NARR-OB-NEST < 12 h; RE-NAM-OB-NEST < 12 h) ( Fig. 7b1 and b2). Our results were similar to the previous study in Chincoteague Bay under Sandy by Beudin et al. (2017). Due to the influence of different wind fields, the corresponding maximum surge level and highest surge durations exhibited significant differences (Fig. 7a1, a2, b1, and b2). Comparisons between RE-NARR-OB and RE-NAM-OB showed no significant variations in wind directions over the bays, whereas the wind speed statistics including the maximum, minimum and mean wind magnitudes of RE-NARR-OB were weaker than those of RE-NAM-OB in the bays (Fig.  6g, h1, and h2). Impacted by the stronger RE-NAM-OB, the corresponding maximum surge level was higher (0 .2 m) and highest surge durations were longer (5-15 h) than those influenced by RE-NARR-OB in the most parts o f C h i n c o t e a g u e B a y, e x c e p t i n t h e s o u t h e r n Chincoteague Bay (Fig. 7a3 and b3).
Most parts of the MCBs exhibit a fairly high-saline pattern (Kang et al. 2017); however, salinity could be low with freshwater present in the upper headwater regions, tributaries and streams near the creeks. The abilities on reproducing salinity forced by different wind fields (RE-NARR-OB and RE-NAM-OB) were alike (Fig. 8b1-e1; Table 3). Salinity at  Grays Creek, Bishopville Prong, and Newport Creek dropped to 0 in response to high river flow (Fig. 8a, b1-d1). However, Public Landing, located far away from river flow, had high and stable salinity (Fig. 8e1: > 25). Relatively high correlation coefficients with reasonable MAEs and RMSDs support that this model was reliable for the salinity simulation (Table 3). The model also captured the observed water temperature under RE-NARR-OB and RE-NAM-OB and water temperature slowly dropped as a result of Sandy (Fig. 8b2-e2). High correlation coefficients with low MAEs and RMSDs were found Fig. 6 Observed (red) and modeled data of cases RE-NARR-OB-IMD (solid cyan), RE-NAM-OB-IMD (solid blue), RE-NARR-OB-NEST (dashed black) and RE-NAM-OB-NEST (dashed gray) for water levels at a S1: Grays Creek; b S2: Bishopville Prong; c the OCI; d S3: Newport Creek; e S4: Public Landing; and f the CI. g Spatially average wind speed comparisons between RE-NARR-OB (dashed black) and RE-NAM-OB (dashed gray) with corresponding statistics analysis including the maximum, minimum and average as well as h1, h2 vector wind speeds time series from RE-NARR-OB and RE-NAM-OB, respectively between the observations and modeled water temperature (Table 3).
Overall, a reasonable agreement between observations and modeled water level, salinity, and temperature lent credence to this 3D hydrodynamic model, which ensured the accuracy of potential flux dynamics studies in the MCBs. Considering a small water level difference between cases RE-NARR-OB-NEST and RE-NAM-OB-NEST and a better agreement between the observed and modeled data in salinity and temperature impacted by RE- Fig. 7 Spatial distributions of a1, a2 the maximum surge and b1, b2 the durations of the "highest surge" of cases RE-NARR-OB-NEST and RE-NAM-OB-NEST, respectively. a3, b3 The surge and durations differences between case RE-NARR-OB-NEST and case RE-NAM-OB-NEST Table 3 Model evaluations with correlation coefficients (Rs), mean absolute errors (MAEs), and root mean square deviations (RMSDs) for water level, salinity, and temperature among cases RE-NARR-OB-NEST, RE-NAM-OB-NEST, RE-NARR-OB-IMD, and RE-NAM-OB-IMD during Sandy

Inlet Exchange Dynamics
Based on the flux at both inlets, the flux amounts at the OCI were smaller than those at the CI because of inlet width differences (Fig. 9). This is consistent with the findings from Mao and Xia (2018), who elucidated the inlet dynamics of the MCBs induced by Hurricane Irene (2011) and demonstrated that the inlet exchange amounts differences were due to the different inlet widths. Under a single factor of open boundary forcing or storm winds, the inlet flux also varied distinctly with each other, indicating that the inlet flux was sensitive to external factors, particularly from October 29 00:00 UTC to October 31 00:00 UTC (Fig. 9). During the passage of Sandy, the MCBs displayed explicit wind variations. Mild (> 10 m/s) and strong (> 15 m/s) northerly winds blowing from the north to south became prevailing in the bays from October 29 05:00 to 17:00 UTC and then wind directions switched ( Fig. 6g and h2). From October 30 03:00 to 15:00 UTC, the bays were dominated by strong southwesterly winds (> 15 m/s) ( Fig. 6g and h2). Thus, two phases (P1: October 29 05:00 to 17:00 UTC and P2: October 30 03:00 to 15:00 UTC) were defined based on the different wind fields to explain the exchange dynamics in response to storm winds and high water levels during Sandy. General flow patterns and spatial water level distributions varied distinctly during P1 and P2. The simulated salt flux shared the same flow pattern as the water flux for each inlet in the same phase (Fig. 10c1, c2, d1, and d2). During P1, mild and strong northerly winds accompanied by rising high water levels tended to pile water up in the southern bays ( Fig. 10a1  and a2), leading to an inflow pattern at the OCI and an outflow pattern at the CI (Fig. 10c1, d1, c2, and d2). Stronger wind effect (outflow) against the effect of rising high water levels from the coastal ocean (inflow) determined the outflow pattern at the CI, indicating the importance of wind fields and rising high water levels from the coastal ocean on the inlet flux ( Figs. 6f and 10a1, a2, c2, and d2). Conversely, higher water levels were present in the northern bays under the strong southwesterly winds during P2, resulting in an outflow/ inflow pattern at the OCI/CI (Fig. 10b1, b2, c1, c2, d1, and d2). In particular, stronger wind-induced flux (inflow) against the outflow by falling high water levels from the coastal ocean (outlfow) determined the inflow pattern at the CI during P2, indicating the importance of wind fields and fall high water levels from the coastal ocean on the inlet flux during Sandy (Figs. 6f and 10b1,b2,c2,and d2). Under the most intense wind effect and rising/falling high water levels from October 29 00:00 to October 31 00:00 UTC, flux was transported seaward at the OCI (water/salt flux − 256 m 3 /s and − 7614 m 3 /s) and landward at the CI (water/salt flux 542 m 3 /s and 17,596 m 3 /s). Overall, the bays gained flux from October 29 00:00 UTC to October 31 00:00 UTC (water/salt flux: 286 m 3 /s and 9982 m 3 /s; Fig. 11c and d).
The wind effect on the inlet flux patterns became more evident when the open boundary forcing of the IMD was turned off (Table 2: No_OBF). Under a single factor of winds, mild and strong northerly winds drew flux landward through the OCI (water/salt flux 60 m 3 /s and 1939 m 3 /s) and seaward through the CI (water/salt flux −2343 m 3 /s and −74,391 m 3 /s) during the early phase of Sandy (P1) (Fig. 10a1, c1, c2, d1, and d2). Owing to the change in wind direction blowing from the southwest during the later phase of Sandy (P2), the flow pattern at both inlets also became reversed with an outflow pattern at the OCI (water/salt flux − 506 m 3 /s and − 15,505 m 3 /s) and an inflow pattern at the CI (water/salt flux 4544 m 3 /s and 147,051 m 3 /s) (Fig. 10b1, c1, c2 d1, and d2). The more favorable wind directions paralleling the shape of Chincoteague Bay (P1: northerly vs P2: southwesterly) Fig. 10 a1, b1 Scattered wind rose plots and a2, b2 the corresponding modeled spatial mean water level distributions during P1 and P2, respectively. c1, d1 and c2, d2 are the mean water and salt flux transport at the OCI and CI of cases Sandy, No_OBF and No_W during P1 and P2. Positive/negative values for flux in c1, c2, d1, and d2 represent the flux onshore/offshore movement contributed to the larger amounts of transport at each inlet during P2 than those during P1 (Fig. 10c1, c2, d1, and d2). The findings of strong wind-induced inlet flux verified the flux dynamics based on sensitivity experiments with the selections of two typical wind directions (northwesterly and southwesterly) in the MCBs by Kang et al. (2017), indicating that the bay orientations aligning with the favorable wind directions play a vital role in determining the flux pattern at inlets.
Considering that Sandy only lasted 2 days in the MCBs, we artificially shortened the storm to 1 day (October 29 00:00 UTC to October 30 00:00 UTC) and prolonged the storm to 3 days (October 29 00:00 UTC to November 1 00:00 UTC) to investigate the exchange dynamics in response to various storm temporal scales (Fig. 11a1, a2, a3, b1, b2, and b3). It was found that the flow pattern changed from inflow to outflow at the OCI in response to longer storm durations due to the longer durations of southwesterly winds (Fig. 11a1, a2, a3, b1, b2, b3, c, and d). The amounts of mean flux at the CI became smaller in response to longer storm durations due to the longer durations of mild southwesterly winds and falling high water levels (< 15 m/s; Fig. 11a1, a2, a3, b1, b2, b3, c, and d). Overall, a net inflow pattern occurred in the bays and the net flux amounts became smaller under longer storm durations ( Fig. 11c and d).

Inlet Salinity
Here, another question was raised: can the relatively high river flow from the creeks affect the inlet salinity under the combined effect of strong winds and high water levels during a storm? To solve this question, we calculated the mean salinity at the inlet transects ( Fig. 2b and c) and found that the salinity at inlet transects displayed temporal variability ( Fig. 12b and c). The bays received the strong river flow on October 29 and 30, 2012 during Sandy and salinity across the OCI transect was still high and stable (> 30) ( Fig. 12a and b). Compared to the high and stable salinity at the OCI transect in case No_RF (Table 2), salinity at the OCI transect could drop to 28 after October 31 due to the residual effect of the relatively high river flow associated with Sandy ( Fig. 12a and b). However, high-saline water (> 32) occurred at the CI and there was no significant variation in the salinity at the CI during the passage of Sandy (Fig. 12c). This is because the CI is far away from river flow and the salinity at the CI transect was less influenced by river flow and mainly determined by salt water intrusions ( Fig. 12c; Fig. 1b).

Conclusions
This study investigated the inlet exchange dynamics of a typical lagoon system (MCBs) under Sandy using a 3D unstructured-grid hydrodynamic model (FVCOM). Through the application of a model nesting method, the simulated water level was greatly improved because the model nesting method provided the necessary remote forcing for storm surge simulations. However, uncertainties still exist in this study, such as the underestimation of the modeled water level at inlets. Nevertheless, the model was capable of reproducing Fig. 11 Time series of wind speeds for a1 1 day, a2 Sandy, and a3 3 days and vector wind speed for b1 1 day, b2 Sandy, and b3 3 days. c Water flux and d salt flux at each inlet and net flux for 1 day, Sandy and 3 days tests.
Positive/negative values for flux in c and d represent the flux onshore/ offshore movement the hydrodynamic characteristics reasonably within the study area and the major conclusions are as follows: (1) Flux between the bays and the adjacent coastal ocean was sensitive to major external forcing agents during a hurricane (storm winds and rising/falling high water levels from the coastal ocean). Storm wind fields and rising/falling water levels from the coastal ocean determined the flux patterns at inlets. In general, the bays gained flux from October 29 to 31, 2012 during Sandy. Sensitivity experiments on various storm temporal scales showed that a net inflow pattern occurred in the bays and the net exchange amounts became smaller under longer storm durations.
(2) Under the effect of high water levels, strong winds, and relatively high river flow, inlet salinities were high and stable (> 30) from October 28 to 31, 2012. Residual high river flow associated with Sandy still had a significant impact on the OCI salinity after October 31 and the OCI salinity even dropped to 28 subsequently. By contrast, the CI salinity was high and stable (> 32) before and after Sandy, indicating that the CI salinity was not affected by river flow because of the long distance between river locations and the CI.
Overall, this study demonstrated the inlet exchange dynamics of the MCBs (a typical lagoon system) under storm winds, rising/falling high water levels from the coastal ocean and river flow of Sandy. A full understanding of the effects of storm winds, rising/falling water levels, and strong river inflow in a lagoon system in response to a passing storm could help us to know the inlet flux pattern and even predict which area of the bay is more likely to suffer from the coastal flooding and inundation with the help of the potential storm surge inundation model. This potential prediction could help policy makers to identify the evacuation routines based on differentially affected areas more precisely and efficiently instead of evacuating all the residents nearby the coastal areas.
Besides, the health of an estuarine ecosystem depends on such estuary-ocean exchange processes, which are associated with water quality, pollutants, and primary production. This flux understanding of the MCBs could be beneficial to environmental managers to predict the movement of the polluted matter (e.g., oil spill, toxic substances), which would allow them to make decisions on water quality policy and other related policies.
More effort and work that look at the effects of hurricanes on turbulent mixing and coastal inundation in the bays are expected. Comparisons between the long term flux and short term flux need to be elucidated in the future. Similarly, this modeling method and work are useful for other shallow bays and lagoon systems, such as Great South Bay in New York, Barnegat Bay in New Jersey, Tampa Bay in Florida, and Perdido Bay in Alabama and Florida, which share similar lagoon characteristics with the MCBs for the study of the hydrodynamic response to hurricanes. Fig. 12 Time series of a river flow from the NMCBs and SMCBs; salinity at b the OCI transect and c the CI transect of cases Sandy (green) and No_RF (orange)