A numerical investigation of the mechanisms controlling salt intrusion in the Delaware Bay estuary

.


Introduction
Estuaries are important coastal systems whereby fresh river water mixes with saline ocean water to create critical habitats for marine animals while also providing natural resources and safe harbors that support human civilization.The distribution of organisms and habitat is largely determined by the temporal and spatial distribution of salinity (Gunter, 1961).Industrial and municipal infrastructure, including drinking water intakes, rely on freshwater from these systems, and are damaged when salinity intrudes into these fresher reaches.Salt intrusion occurs when the normally fresh parts of an estuary are inundated with denser, saltier water from the coastal ocean.Salt intrusion is influenced primarily by the balance of river discharge and tides (Hansen and Rattray, 1965;Fischer, 1972;Lerczak et al., 2006) and secondarily by interactions with bathymetric and topographic features (Ralston et al., 2012) and meteorological conditions (Chen and Sanford, 2009;Aristizabal and Chant, 2015;Gong et al., 2018).In general, the distribution of salinity in an estuary is primarily controlled by mixing between freshwater flow from rivers and saltwater from the coastal ocean.However, during periods of low discharge, low frequency subtidal motions and wind driven transport can push salt further upstream, threatening ecological habitats and industrial and civil infrastructure (Zhu et al., 2020).Bathymetric and topographic features also affect salinity intrusion by changing the momentum balance in the along and across channel directions (Fischer, 1976).The interaction of these mechanisms ultimately determines the distribution of salt in an estuary and has important implications for economic and ecological health.
The Delaware River Basin (DRB) is located on the eastern seaboard of the United States, with over 13 million people, including New York City, New York, relying on this body of water for freshwater for drinking, recreational, agricultural, and industrial use.Severe drought in the 1960s led to dangerously low drinking water reserves coinciding with severe salt intrusion that led to heavy corrosion of industrial and civil infrastructure (DRBC, 2019).Response to the drought included managing the Delaware River to meet a minimum discharge target, hoping to prevent the salt front from intruding upstream into ecologically and economically sensitive areas.In this system, the salt front is the approximate location of the 250 mg Cl − L − 1 or approximately 0.52 psμ averaged over a 7-day period based on drinking water quality standards (New York State, 1991).Minimum discharge flow objectives are imposed to limit the movement of the salt front upriver through managed reservoir releases upriver.These flow objectives are modified during drought conditions to manage both the salt front and low reservoir resources (Flexible Flow Management Program, 2017).Managing reservoir releases to meet a minimum riverine flow target without an understanding of other salt flux mechanisms and their interactions could lead to wasted water resources.
Salinity distribution in estuaries is a balance between advection downstream from river discharge and upstream transport of dispersive processes due to tidal and baroclinic flows (Hansen and Rattray, 1965;Fischer et al., 1979).Generally, increased river discharge drives salt seaward and increases stratification in parts of the estuary.Empirical correlations have shown that salinity intrusion length (L x ) generally scales with river discharge as Q n r and when baroclinic flows dominate the salt transport, n is approximately − 1/3, or L x ~ Q r − 1/3 (Hansen and Rattray, 1965;Monismith et al., 2002).In the Hudson River, L x ~ Q r − 0.38 for low river discharge and L x ~ Q r − 1.19 for higher discharge (Abood, 1974;Wells and Young, 1992).This relationship is weaker in San Francisco Bay, where L x ~ Q r − 1/7 due to basin geometry and the stratification-induced reductions in vertical mixing (Monismith et al., 2002).They concluded that during high river discharge the stratification increases, vertical eddy viscosity decreases, allowing the salt field to move upstream.They found that the along-channel dispersion coefficient was proportional to the tidal velocity and is enhanced during periods of high stratification (i.e., neap tide), and that the landward salt fluxes were confined to the channel.In Delaware Bay, salt intrusion also responds weakly to river discharge (Garvine et al., 1992), likely due to the widening of the bay in the down-estuary direction (Whitney and Garvine, 2006;Aristizabal and Chant, 2014;Geyer et al., 2020).Using a three-dimensional numerical model of Delaware Bay, Aristizabal and Chant (2013) found similar results to Monismith et al. (2002), namely that increased discharge suppresses vertical mixing and promotes landward salt flux that is primarily confined to the main channel.Bathymetry and topography can plan a large role in salinity intrusion, and although landward salt flux is mainly confined to the main channel, Delaware Bay consists of large shoals in the lower bay and bends and islands in the upper reaches that have an important effect on landward salt flux (Fischer, 1972).Salt can be stored in wide embayments during different phases of the tidal cycle, and then reintroduced back into the channel during the next phase of the tide (Okubo, 1973;Fischer et al., 1979).Observations from Aristizabal and Chant (2015) show enhanced salt flux during neap tide in regions where there is an along-channel phase shift between velocity and salinity less than 3.1 h due to the lateral flow structure.Geyer et al. (2020) examined how the channel-shoal geometry in the lower Delaware Bay affects the strength and mechanisms of exchange flow by scaling the estuarine circulation with an "active" cross-section that includes the dynamically relevant geometry.They found net exchange flow during neap tide was 5-6 times that of the spring tide, and that the wide shallow shoals essential for maintaining stratification only make minor contributions to the overall exchange.However, these studies were focused on the lower portion of the Delaware, and not the upper reaches.We know channel-shoal geometry plays a role in landward salt flux, but what about large bathymetric or topographical changes in the upper reaches?
Varying bathymetry and topographical features can contribute to circulation patterns that drive localized maxima in density gradients that can create localized fronts (Simpson and Linden, 1989;Ralston et al., 2008).Geyer and Ralston (2015) demonstrated how fronts can also form in the presence of uniform density, and all that is required is a velocity convergence.In the Hudson River, they found that frontal propagation is greatest during neap tides, propagating landward more than 60 km over several days.This front formation was found downstream of lateral constrictions, formed during late ebb and advected landward during the flood, and primarily during neap tides.Warner et al. (2020) found distinct zones in the Hudson River where mixing occurred, with more mixing occurring predominately during ebb near channel expansions, with mixing during flood at the top of the bottom mixed layer.Several channel constrictions, expansions, and islands exist in the upper reaches of the Delaware Bay but have not been studied in the context of salt intrusion.In this study, we will investigate the role of these features on controlling salinity intrusion up-estuary, particularly in reaches close to important infrastructure.
Meteorological effects, including local wind-driven circulation and barotropic changes in water levels due to offshore storm systems can also have a dramatic effect on salinity intrusion.Longer period fluctuations in sea level (>33 h), hereafter referred to as subtidal motions, due to meteorological systems can produce a net set up/set down of sea level, either locally or from offshore and impact the net salt flux into and out of the estuary.Wang and Elliott (1978) showed that dominant subtidal sea level fluctuations in Chesapeake Bay are the result of the propagation of sea level fluctuations generated by alongshore winds offshore.In San Francisco Bay, Ryan and Noble (2007) showed that approximately 40% of the subtidal fluctuations at the mouth of the estuary are due to large-scale regional wind fields affecting the sea level of the adjacent coastal ocean.In Delaware Bay, Wong and Garvine (1984) found that large subtidal fluctuations at the mouth were forced by winds parallel to shore, indicating Ekman transport into the bay.Wong and Moses-Hall (1998) used month-long sea level and current measurements to further assess the remote and local wind effect on the observed subtidal variability in Delaware bay.They concluded that the remote wind effect was more important than local wind on producing sea level fluctuations within the estuary.However, local wind effects were more important in producing subtidal currents at any given point within the estuary.Both local and remote wind effects can be important drivers of salinity intrusion.The Wedderburn number (Monismith, 1986;Chen and Sanford, 2009) is a ratio of the wind stress to the baroclinic pressure gradient and is a useful parameter for understanding the role of wind on landward salt flux.
In this paper, we use a coupled ocean-wave model to quantify four major mechanisms responsible for salinity intrusion in the Delaware Bay estuary: river discharge, tides, interactions with bathymetric and topographic features, and meteorological events.Previous studies have not looked at the interaction of these dynamics altogether on daily, seasonal, and yearly time scales.This work focuses on the combination of these mechanisms driving the movement of salinity landward and seaward, with particular attention to the salt front metric during periods of low discharge.In Section 2, we describe the model setup and calibration.In Section 3, we describe the results of the model runs.In Section 4, we provide a detailed description of the mechanisms controlling the spatial and temporal patterns of salinity and salt flux in Delaware Bay using a numerical model.Final conclusions and summarizing remarks are made in Section 5.

Study region
Delaware Bay is a large coastal plain estuary that begins at the Atlantic Ocean, expands to a large bay, and then converges into a narrow navigational channel deepened by maintenance dredging to accommodate shipping vessels since the late 19th century (Fig. 1; Moore, 1908).The navigational channel supports the Delaware River Port Complex, one of the busiest ports on the eastern seaboard of the United States and the largest freshwater port in the world (DRBC, 2018).The mean depth of the estuary is 8 m with the channel maximum depth reaching 45 m.The main tributaries are the Delaware River and the Schuylkill River, which account for 76% and 24% of the river discharge above Philadelphia, Pennsylvania, respectively, and ~60-70% of the total freshwater input into the system.The Delaware Bay estuary is connected tidally to the Chesapeake Bay by the Chesapeake and Delaware Canal (C&D Canal, hereafter), where the general inter-basin transfer results in a net transport of fresher upper Chesapeake Bay water eastward toward Delaware Bay (Hsieh et al., 1993).Based on a 41-year discharge record at Trenton, New Jersey [U.S. Geological Survey (USGS) site no.01463500] (location Fig. 1, observed discharge Fig. 2a) the mean annual discharge is 365 m 3 s − 1 , the low monthly average in August is 200 m 3 s − 1 and event scale flows greater than 6000 m 3 s − 1 (USGS, 2022).When discharge drops below approximately 85 m 3 s − 1 at Trenton, NJ, water resource managers release freshwater from upstream reservoirs to manage the salt front location (FFMP, 2017; see Section 2.4).
Tides in the Delaware Bay are mainly semi-diurnal with a tidal range of 2 m and a tidal excursion of approximately 16 km in the upper reach of the bay.The tides are predominantly semidiurnal, dominated by the M2 tidal frequency, with other contributions from the S2, N2, O1 and K1 constituents.The head of tide is located at Trenton New Jersey, approximately 214 km (upstream of the entrance of Delaware Bay, located at Lewes, Delaware (Fig. 1).The Delaware Bay estuary is generally weakly stratified to well mixed (Garvine et al., 1992) with periods of strong stratification in the main channel during neap tides (Geyer et al., 2020).The median location of the salt front (~0.5 psμ; see Section 2.4) is between 107 and 122 km from the entrance to the bay (Fig. 1; DRBC: https://www.nj.gov/drbc/programs/flow/salt-front.html).Average wind speed during 2019 is approximately 6 m s − 1 [NOAA Station ID: Ship John Shoal, NJ 8537121; https://tidesandcurr ents.noaa.gov/stationhome.html?id=8537121], with the strongest winds being 20 m s − 1 coming from the north-northwest, aligned with the channel axis of the bay (Fig. 3).

Modeling framework and model setup
The Coupled Ocean-Atmosphere-Wave Sediment Transport modeling system (COAWST; Warner et al., 2008Warner et al., , 2010) ) was used to simulate the water levels, currents, salinity, and temperature transport in the Delaware Bay.COAWST includes the Regional Ocean Modeling System (ROMS; Haidvogel et al., 2008;Shchepetkin and McWilliams, 2005).ROMS is a three-dimensional, free surface, terrain following hydrodynamic model that solves the Reynolds-averaged Navier Stokes equations assuming hydrostatic equilibrium and Boussinesq approximations, using a finite differences method on an Arakawa C grid.The coupled system also includes the Simulating WAves Nearshore (SWAN).SWAN is a third-generation spectral wave model that simulates wave generation, dissipation, wave-wave interactions, shoaling, refraction, and depth-limited breaking processes (Booij et al., 1999).When run together, data fields are exchanged between ROMS and SWAN using the Model Coupling Toolkit (Warner et al., 2008).Depending on the model case (Table 1) either the ocean circulation model was used alone, or two-way coupled to the wave model.ROMS and SWAN simulations were performed over the same grid.Results are analyzed to determine the dominant physics controlling salinity intrusion.The model simulations that support the findings in this study are available from Cook and Warner (2023).

Model grid
The model grid was based on a version used in previous studies (Chen et al., 2018;Geyer et al., 2020) and was modified to include the C&D canal and extended northward to the head of tide at Trenton, New Jersey (Fig. 1).Bathymetry was based on the CONED dataset (Coastal National Elevation Database (CoNED) Project -Topobathymetric Digital Elevation Model (Thatcher et al., 2016)).The grid extends approximately 100 km in the along channel direction and typical cell spacing is 450 m in the main bay region down to 5 m in the upper reaches.The model was simulated with 16 vertical levels and run with a 2 s-10 s time step depending on the run.SWAN was simulated on the same grid with time step of 5 min.The model runs were simulated from Jan 1, 2019, to Dec 31, 2019.A uniform bed roughness was set to z 0 = 1.5 mm and boundary layer was logarithmic.The k− ε turbulence model is implemented following the GLS method (Warner et al., 2005).

Model forcing
This study improves upon past modeling work (Aristizabal and Chant, 2013;Pareja-Roman et al., 2019;Geyer et al., 2020) by forcing the model for longer duration (1 year) and variable forcing based on observational time series.The coupled modeling system was forced from observations and other model data.At the northern boundary, river discharge and salinity were prescribed from the USGS tidal gage data at Trenton, New Jersey [Fig.2a; USGS site no.01463500] (USGS, 2022).Flows at the C&D canal were estimated from water level and velocity measurements, and salinity was prescribed from the NOAA station at Chesapeake City, Maryland [Fig.2b; NOAA Station ID: 8573927; htt ps://tidesandcurrents.noaa.gov/stationhome.html?id=8573927].
At the southern boundary, tidal forcing was generated from the Advanced CIRCulation model for oceanic, coastal and estuarine waters (ADCIRC) tidal data base (Luettich and Westerink, 2004) and water level as well as depth-integrated transport was imposed across the open boundary.Additional subtidal water level forcing was created from a 33-h low-pass filter of Lewes, Delaware, tidal data [Fig.2c; NOAA Station ID: 8557380; https://tidesandcurrents.noaa.gov/stationhome.html?id=8557380] and Cape May, New Jersey, gages [NOAA Station ID: 8536110; htt ps://tidesandcurrents.noaa.gov/stationhome.html?id=8536110].Data from both sites were interpolated to create a low-pass storm frequency component of the boundary forcing.COAWST forecast data were used for salinity and temperature at the ocean boundary (Warner and Kalra, 2022).
Surface forcings for model runs that included atmospheric dynamics .Shorelines from GSHHG 2.3.7 (Wessel and Smith, 1996).
(Table 1) were driven by the European Center for Medium-Range Weather Forecasts (ECMWF) Reanalysis v5 (ERA5; H Hersbach et al., 2020) with a spatial resolution of 30 km and temporal resolution of 1 h.Data included wind speed and direction, long-and short-wave heat fluxes, relative humidity, and air temperature.Air-sea boundary conditions were imposed using a bulk parameterization approach.Different model configurations generated an increasing level of complexity (Table 1).Simulation A was for ROMS with river and tides, B added the subtidal component, C added surface atmospheric forcings, and D was a coupled ROMS and SWAN simulation.The Model D ROMS-SWAN setup used vortex force formulation (Kumar et al., 2012) with wave dissipation provided by SWAN, and surface roughness was estimated from the wind stress (Charnok, 1955;Craig and Banner, 1994).For more details on the wave model setup and parameters see Chen et al. (2018).No offshore waves were imposed at the boundary.

Model skill
The model skill is based on the method presented by Wilmott (1981), where X is the variable of interest, subscript denoting model or observed quantity, and the overbar represents the mean.A skill of 1.0 represents exact agreement between the model and observations, where a skill of 0 represents complete disagreement.Model skill is evaluated for water levels, currents, and salinity.We also include the Brier Skill Score,  variance in the observations (denominator) (Brier, 1950;Murphy and Epstein, 1989).The maximum Brier score is 1, and a Brier score of zero represents a mean squared model error equal to the variance in the observations.

Salt front 2.4.1. Observation-derived salt front
The salt front location is a derived position determined by the Delaware River Basin Commission (DRBC: https://www.nj.gov/drbc/p rograms/flow/salt-front.html).The front location is defined as the location of the 7-day average 250 mg L − 1 chloride concentration based on a water quality standard for public water supplies (New York State, 1991).The salt front is calculated using specific conductance data measured from four stations: Reedy Island Jetty, DE [USGS site no.01482800], Chester, PA [USGS site no.01477050], Fort Mifflin, PA [USGS site no.01474703], and Ben Franklin Bridge, PA [USGS site no.01467200] (USGS, 2022).The specific conductance data are then converted to a chloride concentration using an established relationship specific to the Delaware River (Preucil and Reavy, 2020).At low salinities, the relative concentrations of ions in brackish water are more variable than in the open ocean.The salt front location is then determined using a log-linear interpolation between the two gages greater than and less than the 250 mg L − 1 chloride isochlor line (Preucil and Reavy, 2020).Although the 7-day is used for water resource management, we use the 1-day average to include the influence of higher frequency events (i.e., storm surge).The salt front is not tracked below river-kilometer (rkm) 86.4 at Reedy Island Jetty, Delaware (Fig. 1).Fig. 2 shows the location of the published salt front (Fig. 2d) compared with discharge at Trenton, NJ [USGS site no.01463500; Fig. 2a].

Modeled salt front
The modeled salt front is determined using the daily average salinity from the bottom vertical cell equal to approximately 0.52 psμ similar to other studies (Philadelphia Water District; April 9, 2019; https://water.phila.gov/pool/files/presentation-to-rfac-2019-04.pdf).

Table 1
Forcing conditions for four model runs using COAWST.

Model
River Discharge Tides Subtidal water levels Winds Waves X, included in model; -, not included in model.

Model-observation comparison
Hourly Water levels vary on tidal and subtidal time scales (Fig. 4).The neap and spring tidal cycles are more apparent in stations closer to the mouth, near Lewes, DE, and Ship John Shoal, NJ, and less obvious in upstream stations near Philadelphia PA.Several wind events and storm systems are apparent in the water level record, specifically in early October, and can be seen propagating throughout the system (Fig. 2b).The model captures these variations in water level well throughout the system.Observed current magnitudes are between 0.5 and 1 m s − 1 and vary slightly between the Delaware River near Philadelphia, PA (upper panel, Fig. 5) and the navigational channel in Delaware Bay (bottom panel, Fig. 5).The model (Model C) represents the variations in current magnitude well, with an overall slight underestimate at both locations.
Time series of salinity vary from 0 to 15 near Reedy Island Jetty, NJ (bottom panel, Fig. 6a) and decrease landward towards Trenton, NJ (top panel, Fig. 6b) where salinity is on average 0.15 psμ.Variability in salinity decreases from Reedy Island Jetty, DE to Trenton, NJ as tidal processes become less important and river discharge dominates.Salinity intrusion is apparent in the fall months as the salinity closer to Trenton, NJ increases above 1 psμ at Chester, PA and around 0.5 psμ at Fort Mifflin, PA.Fig. 6b shows a zoomed in time series of salinity intrusion from October 2019.All the model solutions underestimate salinity slightly (<0.05 psμ) for Fort Mifflin, PA and Ben Franklin Bridge, PA (Fig. 6b, top two panels).
Both the Wilmott and Brier skill scores indicate a general increase in model skill with added forcings for both water level and salinity (Table 2).Water level Wilmott skill scores are all above ~0.95 and Brier scores above ~0.8, with the latter demonstrating that variability in the differences between the model and observations is similar to the variability of the observations.The skill scores for salinity also generally improve with added forcings, however the Brier skill score is much lower in the upper reaches of the Delaware.The observed salinity at Fort Mifflin, PA and Ben Franklin Bridge, PA are relatively small (<0.2 psμ) and could be affected by processes not included in the model (i.e., runoff from winter road salt (Philadelphia Water District (2020)).Therefore, any relatively small mismatch in the model-observation comparison and any increased variability in the observations in this portion of the estuary will lead to smaller Brier skill scores.The Wilmott and Brier skill scores don't show large differences between the runs, and the root-mean square-error decreases from Model A to Model D.

Model salt front comparison
Fig. 7 compares the daily average modeled salt front line (solid lines) with observationally derived salt front line (black dots) for the full year of 2019.Predicted water levels and 33-hr low-pass filtered water levels at Lewes, DE are included at the top of Fig. 7 to demonstrate the springneap tidal cycles and subtidal motions, respectively (also in Fig. 2b).River discharge from Trenton, NJ (also in Fig. 2a) is also shown and included in all model simulations.Note that caution should be taken in comparing the model to salt front locations <rkm 86.4, as the observation derived estimate of the salt front in this region is based on one observational station at Reedy Island Jetty, Delaware (see Section 2.4).Model A (Tides+Discharge) underestimates the magnitude of the movement of the salt front during the winter months Jan-April, shows a similar estimate with Models B, C, and D in May-Sept, and deviates during two storms in early October 2019 (Fig. 7 almost twice as much computing power as Model C but does not increase model performance in predicting salinity upstream (Table 2).

Tides and river discharge
River discharge exerts a leading control on the location of the salt front in the estuary.The spring of 2019 (March-May; Fig. 2) was particularly rainy and accompanied by large discharge events.These fluctuations of river flow drove corresponding large variations in the salt front location as demonstrated in all the model runs (Fig. 7) and obscure the effects of tides and meteorological events.Conversely, during the summer (June-August) and fall (September-November), river discharge decreases to approximately 100 m 3 s -1 and the salt front location moves upriver at a constant rate (Fig. 7; black dots), maxing out at rkm 137.7.Model A (Fig. 7; blue line) diverges from other runs slightly in January-April, and more severely in October, demonstrating the effects due to winds and waves, particularly low frequency water level changes from off-shore events.
Tidal oscillations generate currents that can enhance mixing and transport of salt in the estuary.Fig. 8 shows the hourly location of the modeled salt front (thin blue line) for Model-A (Tides+Rivers), for May through October 2019.Spring-neap tidal modulation is seen in the Lewes, DE predicted water level signal (Fig. 8; top), with vertical gray lines representing the transition from neap to spring tides.The salt front moves up-estuary during the transition from neap to spring tides and tends to stall during subsequent spring tides.An example of this can be seen during end of June to early July (Figs.7 and 8).The greatest landward movements occur during transitions from weak neap tides and are 15, 12 and 5 km in July, August and September, similar to the Hudson River (Geyer and Ralston, 2015).Modeled stratification, the bottom salinity minus the surface salinity [Fig.8, top panel], at the Delaware Memorial Bridge (rkm 112) grows to 0.5 psμ following the apogean neap tide in late June 2019 and stays stratified through September.Throughout the month of July, the salt front "stalls out" around rkm 112 and does not progress upstream to the next bathymetric control point at rkm 116 until the next neap-spring transition that the salt front is able to progress upstream to the next bathymetric control point at rkm 116.Once past rkm 116, the salt front propagates at a constant rate of 0.41 km-day − 1 .

Bathymetric control points
Bathymetric features can affect salinity intrusion by modifying the along-estuary salinity gradient, ds/dx (where ds=(s bot -s top )), through localized convergences in bottom velocities driving frontal zones (Geyer and Ralston, 2015).Model results indicate that there are several frontal zones in the Delaware Bay (Fig. 9).The large salinity gradients between rkm 90-95 are due to the interactions of the main stem of the Delaware Rriver with the C&D Canal (Fig. 1; Fig. 9).This study however is focused on the region north of rkm 100 (Fig. 10a) where several bathymetric features stand out.The first is located between rkm 100-102 where the Delaware River takes a 90 • turn around a headland-like feature to the northeast and includes a deepening of the main channel (Fig. 10b).The salinity gradient here is generally >0.25 psμ-km − 1 with a local minimum just north of rkm 102 where there exists a shallowing in channel depth from 18 m to 14.5 m over 1.5 km (Fig. 10b; top panel; red line).The second feature is a mid-depth lateral expansion of the channel at 7.5 m depth around rkm 106 (Fig. 10b).The third is a constriction in the river at rkm 111, where just upriver of this constriction exists a second deepening of the channel around rkm 112 (Fig. 10c).This is where the Cristina River joins the main channel (not considered in this study).Upriver of this deepening there is a shallowing in the channel depth from 19 m to 14 m over 1 km, (Fig. 10b; top panel; red line).The last notable feature considered in this study is between rkm 113 and 115, where another mid-depth expansion occurs at the 14 m isobath and an expansion of the river which now includes a series of flats to the east of the main navigational channel (Fig. 10c; top panel).
Model results indicate that there are several frontal zones of salinity in the Delaware Bay (Fig. 9).The large salinity gradients between rkm 90-95 are due to the interactions of the main stem of the Delaware River with the C&D Canal (Fig. 1; Fig. 9).The frontal zones in March through May correspond to large discharge events that push salt oceanward and remove the stratification in this part of the estuary.In late May, the salinity gradient grows due to the decrease in river discharge coinciding with a neap-to-spring tidal transition.The following discharge event in early June creates an exchange flow that drives a dramatic increase in ds/dx in this part of the estuary.In late June, the salt front moves on average 1.5 km-day − 1 until it reaches rkm 106 where it stays until late July.The stall of the front at this location is due to the lateral expansion at mid-depth.
In late July-early August the salinity gradient propagates from rkm 104 to rkm 111 at a rate of 0.7 km day − 1 (Fig. 11c).The modeled salt front shows similar behavior between Model A (Tides+Discharge) and Model C (Tides+Discharge+Subtidal Water Levels+Winds) suggesting this is the result of tidal action, since river discharge is relatively low during this period.The second propagation period from rkm 111 to rkm 113.5 at a rate of 0.6 km day − 1 , has lower tidal currents, low discharge, and low wind activity (Fig. 11), with similar behavior between the model runs, suggesting again this is due to tidal action.The third propagation period from rkm 114 to rkm 118 at a rate of 1 km day − 1 occurs during similar tidal and river conditions as the second propagation period, however, there are moderate winds (3 m s − 1 ) out of the north (Fig. 11a).There is also a deviation between Model A and Model C suggesting this propagation is due to wind action.

Meteorological effects
Meteorological storm systems affect the movement of the salt front in Delaware Bay in two ways.We define a "remote" meteorological intrusion event as those from low-pressure systems that are generally offshore and do not occur over the estuary, and a "local" meteorological intrusion event as those from low-pressure systems that occur within the estuary.There were two major storm events in the fall of 2019, the "remote" subtropical storm Melissa between October 8-12, and a "local" S.E.Cook et al. bomb cyclone between October 16-17.
Subtropical storm Melissa developed from a cold front in the southwestern Atlantic on October 6, then evolved to a nor'easter on October 9 and upgraded to a tropical storm on October 11 and finally dissipated on October 14 (Berg, 2019).Because the storm track stalled off the mid-Atlantic coast, water levels increased along the coast (Fig. 2b) due to winds and barometric effects, leading widespread coastal flooding and a salinity intrusion event in Delaware Bay (Fig. 2a).During this remote event, the propagation of the salt front landward averaged 2.5 km day − 1 (Fig. 7).Model results for Model C and Model D match observed salt front values, are overestimated by Model B and completely missed by Model A. A few days later a powerful nor'easter occurred with wind gusts up to 40 m s − 1 (Fig. 2b) that set new low-pressure records in the region with precipitation accumulation between 2 and 4 inches.This storm system was designated a bomb cyclone because the central pressure in this system dropped by at least 24 mbar in 24 h or less.This local event caused a large drop in water level and corresponding migration downstream of the salt front, which receded on average 4 km day − 1 .Model results for Model C and Model D match observed salt front values, are again overestimated by Model B and completely missed by Model A. For both events described above, results indicate that to capture meteorological based intrusion events it is important to include winds, and subtidal water levels alone are not enough.

Salt intrusion length and river discharge
Theoretical predictions for salinity intrusion length, L x , as a function of river discharge, Q r , follow a power-law relationship, typically L x ≈ Q n r , with n ~ − 1/3 (Hansen and Rattray, 1965).However, this scaling applies to gravitation circulation with a constant cross-sectional area, where the "diffusive fraction" of salt transport is assumed to be zero.Aristizabal and Chant (2013) used a three-dimensional model in the Delaware Bay to determine the salt intrusion length for different constant river inputs: 350, 650, 1,000, 1,300, 1700 and 3000 m 3 s − 1 , and found that for the isohaline s = 0.1, L x (t, s = 0.1) = 206.1 Q − 0.113 r .In 2019, we find a constant discharge of 100 m 3 s − 1 existed for approximately 40 days from August 30th to October 8th (Fig. 12a).During this time the salt front propagated upstream at a constant rate of 0.4 km day − 1 .If we apply a linear regression to this observational subset, we find that L x (t, s = 0.5) = 421 Q − 0.27 r (Fig. 12b and c; r 2 = 0.55), close to the theoretical value, n ~ − 1/3 and greater than model estimates from Aristizabal and Chant (2013).Looking at discharge data for 2019 it is rare to find steady-state river discharge values greater than 100 m 3 s -1 therefore this scaling most likely only applies to seasonal low flow conditions in the Delaware, when the salt front is north of rkm 100.Applying least squares linear regression to the entire derived salt front dataset , we predict the salinity intrusion length, L x (t, s = 0.5) = 260 Q − 0.16 r (r 2 = 0.55), which deviates from n ~ − 1/3, which is not unexpected as the period of record data include a substantial amount of non-steady state conditions.If we apply a filter to the data for time periods where discharge is less than 100 m 3 -s − 1 and the salt front is greater than rkm 112, we find L x (t, s = 0.5) = 333 Q − 0.21 r (r 2 = 0.27).The low r 2 value is likely due to the large spread in the data over the almost 60-year period.In this case the n-value is closer to − 1/5, which aligns with studies in other estuaries for lower flow periods.Studies using long term datasets in the Hudson River have found n to be between − 1/5 (Oey, 1984) and − 1/3 (Abood, 1974) for low flows, and ñ -1 for high flows (Abood, 1974).Using almost 20 years of data in San Francisco Bay, Monismith et al. (2002) found that n is closer to − 1/7 than − 1/3.

Bathymetry
Several bathymetric control points exist in the upper reaches of the Delaware Bay estuary (see Section 3.2.2) and are apparent in the time series location of the salt front (Fig. 7).The salt front spends most of the year south of the first feature located at rkm 100-102 (Figs. 7 and 10b).The second is a lateral expansion at 7.5 m depth around rkm 106.In mid-October at the beginning of ebb tide (Fig. 13c1), there is a pocket of low velocity in this region along with a sharp along-channel gradient in salinity (Fig. 13c2).This lateral expansion is the cause of the increased ds/dx in this region, also seen in Fig. 9.The next control point at rkm 112 (Fig. 10c) is of particular interest as after crossing this point the salt front continues to advance till a large discharge event pushes salt south (Fig. 7).Between rkm 111 and rkm 114 there is a shallowing of the main channel and a lateral expansion which drives a large salinity gradient and a front form (Fig. 13b2-g2).This frontogenesis (Simpson and Linden, 1989;Geyer and Ralston, 2015) is common in other estuaries and is an important process responsible for salinity intrusion.This behavior exists for all model runs (Fig. 7) and therefore regardless of forcing is an important control on salt intrusion for Delaware Bay estuary.

Meteorological
The modeled estimate of the location of the salt front improves when including the influence of winds (Model C; red line Fig. 7).These atmospheric effects are particularly strong in the winter-spring (Jan-April) and during the intrusion event in October.Fig. 14 focuses on three time periods in January, March, and late September-October to illustrate the effects of local vs.remote meteorological systems on salinity intrusion events.

Local effects
When the winds are >10 m s − 1 and direction from the northnorthwest (~270-300), there is a reduction in the subtidal water level at Lewes, DE, and a corresponding salt front seaward migration (for example: January 10, January 21, March 23, October 16-17; Fig. 14 gray boxes).The drop in water level, and salt front location is deemed a "local effect", where local winds and a low atmospheric pressure system over the estuary drive water and salt out of the estuary.In both January (Fig. 14a) and March (Fig. 14b), including meteorological effects -Model C (red line) improves the estimated salt front location, and salt intrudes further into the estuary.Whereas in October including winds suppresses the movement of salt up estuary and still improves the model prediction of salt front location (Fig. 14c -red line).
Between October 16-17, a mid-latitude cyclone moved across the United States and developed into a bomb cyclone, with an atmospheric pressure drop and winds equivalent to a category 1 hurricane.Corresponding strong local winds from the north-northwest, and severe drop in water levels resulted in a large flux of water and salt leaving the Delaware Bay (Fig. 14c).The salt front moved seaward on average 4 km d − 1 over the following week.Precipitation and increased discharge at Trenton, NJ (Fig. 2a) contributed to this seaward migration of the salt front.In order to tease out the role of wind stress versus scaled baroclinic pressure gradient force we can use the Wedderburn number (Monismith, 1986;Chen and Sanford, 2009), where t w is the along-channel wind stress (positive is up-estuary; nega-tive is down-estuary winds), L is the length of the estuary (215 km), Δρ is the density change over L, g is the gravitational acceleration, and H is the averaged depth.The Wedderburn number over the entire estuary for Model C for the month of October (Fig. 15c) indicates that for the bomb cyclone event, the wind stress is dominant over gravitational circulation (W > 0).This is also true for the other local wind events in January (Fig. 15a) and March (Fig. 15b).

Remote effects
When subtidal motions at Lewes, DE show a large rise in water levels, that generally corresponds to an offshore low-pressure system, or a "remote" effect (Fig. 14a-c).Two examples of this effect are the event on January 14 and subtropical storm Melissa October 8-October 12. On January 14, there was an offshore low-pressure system that setup a pressure gradient leading to an increase in water level and corresponding flux of water and salt at Lewes, DE (Fig. 15a; bottom panel).If the increase in water level due to subtidal fluctuations was the only mechanism driving salinity intrusion, then the model run with subtidal motions Model B (yellow line, bottom panel, Fig. 14a) would be substantially different from the tidal only model, Model A (blue line, bottom panel, Fig. 14a).However, including both winds and subtidal motions, Model C (red line, bottom panel, Fig. 14a) matches closer to the derived salt front location (black dots, bottom panel, Fig. 14a).This could be due to the moderate winds (~5 m s − 1 ) from the south-southeast (~100-180 • ), driving more water and salt into the estuary.Overall, the January 14 event did not lead to a large intrusion event because of the location of the salt front (~rkm 90) was further upstream than any water movement from the lower estuary.Similar dynamics are seen in March with a more consistent match of Model B (red line, bottom panel, Fig. 14b) to the derived salt front.In this case, there are consistent winds between 5 and 10 m s − 1 from the south-southeast (~100-180 • ) over the estuary, driving water and salt into the upper estuary.
Between October 8-12, subtropical storm Melissa stalled offshore of the east coast of the United States for several days leading to a large pressure gradient offshore, higher subtidal water levels at Lewes, DE, and a salinity intrusion event in the upper Delaware Bay (Fig. 14c).The daily average salt front moved on average 2 km d − 1 , with the highest location at rkm 137.This was the only significant remote effect in 2019.Aristizabal and Chant (2015) found similar results in their observational study from 2011.During this period, the advective salt flux, mainly driven by changes in sea surface height, dominated over the steady shear dispersion and tidal oscillatory salt flux.They found that changes in sea surface height coincided with named tropical storms in the North Atlantic, driving low-frequency flux variability that overpowered the seaward flux associated with river discharge.The present study confirms their findings and therefore, to determine transport and distribution of tracers such as salinity, in the estuary, both local and remote effects should be considered.Overall, including the local meteorological winds improves the model ability to predict the salt front location, and the wind either enhances or suppresses the movement of the salt front.

The impact of waves on salinity intrusion
Observational studies of waves in the Delaware Bay estuary are difficult, and therefore we employed a coupled ocean-wave model to study the effects of wave-driven transport on salinity intrusion.Past studies show that offshore swell does not propagate into the estuary due to bathymetric refraction (Kukilka et al., 2017), and most waves are locally generated and can enhance transport depending on wind direction (Pareja et al., 2019) and bathymetric features (Chen et al., 2018).We follow Pareja et al. (2019) and use ROMS+SWAN to determine the effect of wave-driven transport on salinity intrusion.Model results for wave height and direction for 2019 indicate that larger waves >1 m are primarily from the south and north-northwest, which correspond to local wind events (Fig. 16).
Our results (Model D) indicate that including wave-driven processes does not improve the model ability to predict salinity intrusion in the upper reaches of the Delaware Bay estuary (Figs. 7 and 14), particularly upriver of Philadelphia, PA (Table 2).Waves and wave-driven transport occur during large wind events, primarily during storms.Although storm events can drive salinity intrusion events (see Section 4.3), they do not always include large waves, like during remote event subtropical storm Melissa, where winds were <10 m s − 1 in Delaware Bay (Fig. 14c) and waves were <1 m (Fig. 16).

Summary and conclusions
A coupled wave and circulation model (COAWST) was utilized to evaluate the role of river discharge, tides, subtidal water levels, winds, and waves on salinity intrusion in the Delaware Bay in 2019.Modeled water surface elevation and salinity compared well with observations from several NOAA and USGS gauges throughout the domain, with model skill for salinity showing improvement by adding atmospheric forcing, specifically the local wind fields.This is particularly true during the winter months (nor'easters) and fall months (tropical and subtropical storm systems) when salt fluxes are affected by local wind field.Local storm systems (e.g., bomb cyclone) generally led to salt flux out of the system and a seaward migration of the location of the salt front.Remote storm systems (e.g., subtropical storm Melissa) tended to create a pressure gradient offshore that increased salt flux into the Delaware Bay estuary and led to salinity intrusion upstream.Although including wave-based physics (Model D) affected the distribution of salinity in the lower Delaware Bay estuary, it did not improve the performance of the model in predicting the location of the salt front.
The role of discharge, tides, subtidal water levels, winds, and waves during dry conditions is not well understood.Our study year (2019) was a relatively high discharge year, with large variability in discharge during the winter months (Jan-April).Future work could focus on how these dynamics change during a drought year like 2016-2017.Future models could also include the effects of the Christina River, the third largest input of freshwater to the system near rkm 112, which is also the location of one of the bathymetric control points.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 1 .
Fig. 1.Model domain with topography and bathymetry (meters) and alongestuary distance from the mouth (Lewes, DE) of the bay to the head of tides (Trenton, NJ).Shorelines from GSHHG 2.3.7(Wessel and Smith, 1996).

Fig. 3 .
Fig. 3. Wind statistics in Delaware Bay based on observations at Ship John Shoal, NJ (NOAA Station ID: 853712) from January 2019 to December 2019.Wind rose with wind magnitude at 10 m above the surface (U10) and direction distribution.The wind rose highlights the predominance of down-estuary and up-estuary wind direction for 2019.

Fig. 4 .Fig. 5 .
Fig. 4. Comparison between time series of Model C modeled water level (m) for NOAA gages for sample time period August 15, 2019 to September 15, 2019.

Fig. 6 .Fig. 8 .
Fig.7compares the daily average modeled salt front line (solid lines) with observationally derived salt front line (black dots) for the full year of 2019.Predicted water levels and 33-hr low-pass filtered water levels at Lewes, DE are included at the top of Fig.7to demonstrate the springneap tidal cycles and subtidal motions, respectively (also in Fig.2b).River discharge from Trenton, NJ (also in Fig.2a) is also shown and included in all model simulations.Note that caution should be taken in comparing the model to salt front locations <rkm 86.4, as the observation derived estimate of the salt front in this region is based on one observational station at Reedy Island Jetty, Delaware (see Section 2.4).Model A (Tides+Discharge) underestimates the magnitude of the movement of the salt front during the winter months Jan-April, shows a similar estimate with Models B, C, and D in May-Sept, and deviates during two storms in early October 2019 (Fig.7).Model B (Tide-s+Discharge+Subtidal Water Levels) performs better than Model A in Jan-April but overestimates the salt front migration upstream in October 2019.Model C (Tides+Discharge+Subtidal WaterLevels+- Winds)   and Model D (Tides+Discharge+Subtidal Water

Fig. 9 .
Fig. 9. [top] modeled depth averaged velocity m s − 1 at mouth of estuary (gray) low pass filter depth averaged velocity m s − 1 (black) and Trenton discharge m 3 s − 1 (blue).Yellow box indicates subtropical storm Melissa, and the orange box indicates the bomb cyclone.[bottom] spatially averaged salinity gradient, ds/dx psu km − 1 for Model C (Tide-s+Discharge+Subtidal Water Levels+Winds) with modeled salt front.Horizontal dashed lines are at rkm 100 and rkm 116.(For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Fig. 10 .
Fig. 10.(a) Zoomed in grid from rkm 80-120.(b) Feature #1 (rkm 100-102), red arrow and line indicate a shallowing in the channel depth, and Feature #2 (rkm 105).Solid white and dotted white are the 7.5-m and 16-m contours, respectively.(c) Feature #3 (rkm 110-113), red arrow and line indicate a shallowing in the channel depth, and Feature #4 (rkm 113-115).Solid white and dotted white are the 14-m and 18-m contours, respectively.Colorbars indicate topography and bathymetry in the model domain.(For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Fig. 12 .
Fig. 12.(a) Observed salt front location (black circles) and discharge at Trenton, NJ (blue line) and subset data for constant discharge from August 30th -October 8th (yellow).(b) discharge vs. salt front location for 2019, subset data in yellow, and nonlinear regression in red (r 2 = 0.55).(c) zoom of subplot (b).(For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Fig. 13 .Fig. 14 .
Fig. 13.Model C results for the second ebbing tide during October 11, 2019.(a1, a2) modeled water levels at rkm 107 with observed discharge (blue line) at Trenton, NJ [USGS site no.01463500; USGS, 2022] with circles corresponding to the panels below (b1-g1) along-channel velocities (b2-g2) salinity, black lines are isohalines.Vertical gray line indicates rkm 113.(For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

Fig. 15 .
Fig. 15.Wind speed and direction from Ship John Shoal station (a,b,ctop panel) and along-channel stress (black line) and Wedderburn number (a,b,cbottom panel) for (a) January, (b) March, and (c) October 2019.

Fig. 16 .
Fig. 16.Wave rose for Model D from the middle of Delaware Bay for 2019.

Table 2
Model skill and root mean square error (RMSE) for water level, salinity, and currents.