Circulation & Exchange Within Shelf & Estuarine Waters Driven by the Atmosphere, Tides and Buoyancy

The proximity of human populations to the coast renders the studies of the coastal ocean important due to impacts by natural hazards, environmental management and material transport. This dissertation will quantify the impacts of the atmosphere, tides and buoyancy forces in Rhode Island Sound (RIS) and Narragansett Bay (NB) at a variety of scales. Observational analysis and numerical modeling provide comprehensive tools to study resonance between the atmosphere and ocean beneath synoptic storm systems; to describe circulation in RIS; and to study the effect of tides on stratification in NB. The characterization of the shallow shelf circulation under various forcing mechanisms provides information useful for hazard monitoring. The study of momentum propagation over the East Coast continental shelf is discussed in Chapter 1. We determine the role of atmospheric pressure in creating high frequency surface gravity waves extending over a large geographic range. We find that shallow water waves are generated near the continental shelf break under long-lived squall lines. Buoyancy, tides and wind driven circulation influencing RIS circulation and hydrography are examined using observations in Chapter 2. We find that the largest contributor to circulation are tides and the largest controller of stratification is solar insolation. At monthly time-scales we find a cyclonic coastal current present along the periphery of RIS, enhanced during the summer as a result of baroclinic and barotropic pressure gradients. During the winter the coastal current is reduced and, in some areas, not detectible. The reduction is thought to be caused by smaller baroclinic forces and barotropic forces that oppose the cyclonic circulation. Chapter 3 expands beyond the spatial-temporal limitations of data moorings within NB by using numerical hydrodynamic models to reveal 4-D aspects of estuary dynamics. We find that much of NB has a maximum stratification during either low or high tide resulting from a combination of straining, advection and mixing. The predicted tidal change in stratification is confirmed by observations from buoys located in NB. The coastal area of New England, an area of great economic, recreational, and environmental importance, is an ideal area to study contemporary hydrographic and dynamic processes. This dissertation applies innovative analytical and numerical techniques to investigate the relative roles of atmospheric, buoyancy, and tidal forcing in determining circulation in the shallow shelf sea of RIS and in the NB estuary.


Introduction
Tsunamis are most often caused by sudden movement of the seafloor due to submarine earthquakes, landslides, or volcanic activities. Tsunami-like events created by disturbances on the ocean surface are less well known, but their existence is welldocumented in areas such as the Balearic Islands, the Adriatic Sea, South Japan, New Zealand, northeastern North America, the United Kingdom, and northwestern North America 1,2,3,4,5 . Common to all of these ocean surface-generated tsunami-like events is forcing by an atmospheric surface pressure anomaly moving over a relatively shallow body of water 3 . Hence, this type of tsunami-like event is typically classified as a "meteotsunami." A large sea level anomaly associated with a meteotsunami can occur through a variety of atmosphere-ocean resonance mechanisms, including Greenspan, Shelf, and/or Proudman resonance 3 ; for the meteotsunamis examined in the current study, Proudman resonance has the largest influence on the sea level anomaly. Proudman resonance exists when the ground speed of the atmospheric pressure anomaly (U) matches the phase speed of meteotsunami waves (C), which travel as shallow water waves so that C = (gh) 1/2 , where g is the gravitational acceleration (9.8 m s -2 ) and h is the water depth 6 . For an atmospheric pressure anomaly propagating at a common ground speed of 20-40 m s -1 , the corresponding resonant water depth is 40-160 m, within the depth range of most continental shelves.
Although, not widely known, meteotsunamis along the east coast of North America are not rare events. There have been at least two documented meteotsunami events each year with sea level oscillations of 0.1-1 m along the U.S. East Coast from 2006 to early 2012, although the atmospheric forcing for these meteotsunamis has not been investigated in detail 7,8 . Mercer et al. 2 describes two tropical cyclones, Helene (2000) and Jose (1999), which generated meteotsunamis (with wave heights up to 3m) as the low-pressure anomalies at the center of the tropical cyclones rapidly propagated across the Grand Banks. These meteotsunamis then reflected back toward the coast once the storms crossed the shelf break into deeper water. Since limited observations of the pressure distributions in the tropical cyclones were available, Mercer et al. 2 estimated the pressure distribution using a simple analytical model.
Tropical cyclones do not appear to be the most common cause of meteotsunamis along the U.S. East Coast. Churchill et al. 9 describes a mesoscale convective system (MCS), which propagated southward along the east coast of Florida, generating a 3-m meteotsunami (recorded at Daytona Beach) that was forced by a high surface pressure anomaly under the squall line portion of the MCS. They used hourly radar reflectivity and barometric pressure readings from sparse stations to illustrate the relationship between radar reflectivity and the high atmospheric pressure anomaly. Churchill et al. 9 recorded a positive pressure anomaly of ~2 hPa, and through their estimations of atmospheric forcing, they concluded that, in addition to Proudman resonance, bathymetric effects (including wave refraction off an underwater ridge) may have played a large role in generating such a high amplitude meteotsunami.
Some MCSs include a fast-moving, long-lived, quasi-linear squall line (i.e., derecho), which produces strong winds and has a well-defined surface pressure anomaly signature (Fig. 1a). Derecho-producing MCSs are not uncommon in the interior of the U.S., numbering on average ~20 per year with a possible upward trend in frequency 10 . Most common in the months of May, June, and July, derechos frequently occur in groups along the eastern half of the U.S. 10,11 . On land, the most destructive impact of a derecho-producing MCS is typically straight-line wind damage, but as the MCS passes over the ocean, the potential initiation of a meteotsunami creates a different hazard.
On June 13, 2013, two high-frequency sea level oscillation (i.e., meteotsunami) events hit the U.S. East Coast. The maximum sea level oscillations recorded at National Oceanic and Atmospheric Administration (NOAA) tide gauges for the first and second meteotsunamis were 0.59 m at Providence, RI and 0.44 m at Atlantic City, NJ, respectively (i.e., stations 27 and 19, respectively, in Table 1 and Supplementary   Fig. 1). Earthquakes and landslides can be dismissed as the likely causes of these events because no earthquake was detected near the coast around the times of the events, and a subsequent offshore survey near Hudson Canyon by NOAA found no evidence of a significant submarine landslide (Jason Chaytor, personal communication  13 . These measurements provided the first detailed account of the magnitude, dimension, and duration of the atmospheric pressure anomalies that were ultimately responsible for the meteotsunamis along the U.S. eastern seaboard, once the pressure anomalies moved offshore. Along the coast, NOAA tide gauges are used to monitor sea level, but many also measure coincident atmospheric sea level pressure, with 0.1 hPa precision and a six-minute sampling frequency 14 (Table 1 and Supplementary Fig. 1). All barometric pressure anomalies are calculated by demeaning the recorded pressure values and then high-pass filtering above five hours to highlight the high-frequency pressure anomalies (Fig. 1).
A derecho-producing MCS typically includes a quasi-linear group of thunderstorms, consisting of well-developed convective and stratiform cloud regions.
Since radar measures reflectivity off water and ice particles, areas with the heaviest precipitation, such as beneath the convective towers (i.e., squall line) in a MCS, have the highest reflectivity, as shown in Fig. 1. In a well-developed MCS, the downdraft is often directly below this area of high reflectivity, creating a mesohigh 15 . In addition to high pressure under the convective towers, warm air preceding the convective towers is forced upwards in the gust front, creating a precursory mesolow ahead of the mesohigh 15,16 . Also, another mesolow often follows the convective line in the stratiform region of the MCS (i.e., wake low), resulting from adiabatic warming of the descending air mass in the rear of the MCS 17 . This mesolow/mesohigh/mesolow pressure combination is well illustrated by TA station P61A as the first MCS passed over the coastline (Fig. 1). Here, the first mesolow pressure anomaly of -0.5 hPa corresponds to a narrow area of stratiform precipitation, indicated by moderate radar reflectivity. Immediately behind this first mesolow is an area of strong reflectivity, up to 45 dBZ, with a corresponding pressure anomaly of 2.8 hPa. Finally, as the reflectivity decreases to 20 dBZ, the pressure anomaly associated with the wake low reaches a minimum value of -1.8 hPa. The MCS creates a peak-to-trough anomaly of 4.6 hPa. Radar images from station KDOX show that the cross-sectional reflectivity of the first MCS extends about ~90 km front-to-back, roughly east-to-west, as the MCS passes over TA station P61A (Fig. 2e). While the highest reflectivity occurred concurrently with the highest positive pressure anomalies ( Fig. 2a-d), the mesolow/high/low sequence may further enhance meteotsunami generation relative to a scenario whereby the precursor mesolow and wake low were not present. By examining many TA stations and NOAA tide gauges simultaneously, the highest pressure occurs beneath the MCS's convective squall line (Fig. 3), though interpolation of pressures at the stations broadens the apparent pressure anomaly along the convective line. The second MCS passed further south ~6 hours after the first MCS, arriving at Bishops Head, MD with an associated peak-to-trough pressure anomaly of ~5 hPa (Fig. 2d). The TA station O61A (Fig. 2a), which is located on the extreme northern fringe of the second MCS near the synoptic scale low pressure center, did not show significant atmospheric pressure anomalies and did not have reflectivity measurements above 40 dBZ.
The maximum peak-to-trough atmospheric pressure anomalies during the two MCSs at all of the TA stations and tide gauges are shown in Supplementary Fig. 4 Fig. 7 and described in the Method section. As the first MCS on June 13 th ,2013 moves across the interior of the United States, the position of the center of the atmospheric pressure anomaly is never more than 150 km away from the center of the radar reflectivity values. Second, to confirm that the storm moves at the same velocity as the atmospheric pressure anomaly, the propagation velocity derived from both radar reflectivity and atmospheric pressure anomalies is shown in Supplementary   Fig. 8 and described in the Method section. These velocity estimations are within 5 m/s of one another. The atmospheric pressure anomaly velocity was estimated by temporal averaging the atmospheric pressure anomalies observed on the 90 km spaced TA array. We believe that the discrepancy in the velocity estimations is due to the station spacing of the TA array, which spatially under samples the atmospheric pressure anomaly with an average wavelength of 150 km at any given time. The estimation of the atmospheric pressure anomaly velocity would be improved with the increase of station spacing or modeling of the phenomena. Finally, over a given continental location (in this case at TA stations), it is shown that the arrival of the maximum observed atmospheric pressure anomaly and the arrival of the maximum observed radar reflectivity occur at the same time ( Supplementary Fig. 9), verifying that the average velocities of the MCS and pressure anomalies must be similar. These observations of location, speed, and timing provide at least to first order the basis to infer MCS pressure anomalies offshore using radar reflectivity data, where direct atmospheric pressure measurements at sea level are currently lacking.
Once offshore, the propagation speeds and back-azimuths of the MCSs are estimated from radar reflectivity (Method). The first MCS maintained an average speed of 22 m/s as it propagated over the Atlantic Ocean at ~1500 UTC and dissipated by ~2000 UTC (Fig. 4) With an average wavelength on the order of 150 km, the pressure anomaly may be too wide to exert its force fully on Delaware Bay, which has a maximum width of 60 km.
The maximum water depth in the bay is 30 m 20 , corresponding to a maximum shallow water wave speed of ~17.1 m s -1 . This speed is lower than the forward speed of the atmospheric pressure anomaly during the first MCS (~22 m s -1 ) over the bay, implying a supercritical flow regime, in which a positive sea level anomaly corresponds to a positive atmospheric pressure anomaly (Method), as observed at the tide gauges in the bay (Fig. 5, Supplementary Fig. 10, Supplementary Fig. 12). Further investigation, perhaps using a numerical model, is required to determine the relative importance of atmospheric pressure and wind forcing within Delaware Bay.  (Fig. 4) as it propagated towards the shelf break, which limited the region of Proudman resonance along the shelf and reduced the amplitude of the second meteotsunami in areas decoupled from the atmospheric pressure forcing.

Discussion
The observed correlation between radar reflectivity and atmospheric pressure anomalies recorded by TA and tide gauge stations (Figs. 2 and 3) indicates that the largest positive pressure anomalies occur where reflectivity is >40 dBZ. These areas also experience mesolows before and after passage of the convective line of the MCS.  1,2,3,4,5 .
For near-term meteotsunami prediction along the U.S. East Coast, three stages in the forecast process are suggested: (1) using land-based pressure measurements to monitor the magnitude of the atmospheric surface pressure anomalies as a radarindicated MCS propagates towards the coast from the interior of the U.S., (2) monitoring radar reflectivity as the MCS propagates off the coast to the potential geographical area where Proudman resonance may occur, and (3)  Developing and then operationally implementing accurate numerical models of meteotsunami waves under atmospheric pressure (and perhaps wind) forcing is the next step towards understanding and predicting when and where MCS-generated meteotsunamis may occur. Ports, harbors, and bays may all have varying risk, and characterizing areas that may have extreme resonant mechanisms is also required.
Places like Atlantic City, NJ, documented both in the present study and in others 7 , appear to have an especially large ocean response to MCS forcing. Along with numerical modeling, monitoring high-frequency atmospheric pressure anomalies in the interior of the continental U.S. and along the coastline is essential for early warning of potentially hazardous meteotsunamis.

Method
Pressure and Radar Anomaly Location. Estimates of the center of the atmospheric pressure anomalies and radar reflectivity of the first MCS front were found using pressure data from the TA stations and NOAA/WSR-88D radar station KPBZ. Station KPBZ was used because of the dense coverage of TA stations surrounding it.
Atmospheric pressure time series were demeaned, high pass filtered over 5 hours, and squared so positive and negative anomalies did not cancel. Pressure data were averaged a half an hour before and after a given time step and then linearly interpolated to a 10 km grid. Radar data were down-sampled to a 10 km by 10 km grid. Radar reflectivity >40 dBZ was used to isolate the MCS front. The centers of the pressure anomalies and radar reflectivity measurements were calculated at fiveminute intervals using a gray-level-weighted average of values. Results are displayed in Supplementary Fig. 7.

MCS Propagation
where Froude number Fr = U/C, g is gravitational acceleration, and ρ is the water density. From equation (1), a positive (negative) atmospheric pressure anomaly yields a positive (negative) sea level anomaly under the supercritical condition (Fr > 1), where the pressure anomaly is moving faster than the tsunami waves. In contrast, a positive (negative) atmospheric pressure anomaly yields a negative (positive) sea level anomaly under the subcritical condition (Fr < 1), where the pressure anomaly is moving slower than the tsunami waves. As Fr approaches 1 (the critical condition), Proudman resonance occurs, resulting in a large sea level anomaly. The steady-state sea level anomaly becomes infinite when Fr =1 in the idealized situation described by equation (1).

Reflected Travel Times.
To estimate arrival times at tide gauges [19][20][21][22][23][24][25][26][27][28][29][30] ( Supplementary Fig. 14 where TR is the travel time from the continental shelf break to the shore and back to the buoy, and TT is the travel time from the continental shelf break to the buoy. Here, the continental shelf break and the shore are defined as 140 m depth and 0 m depth, respectively.

Figure 4| Estimated (a) propagation speed and (b) back-azimuth for both MCSs.
The speeds of the first MCS (black) and second MCS (gray) are calculated starting when the MCS moves over Chesapeake Bay and ending when the MCSs dissipate or are out of the range of land based radar. Time is in hours from 0000 UTC on June 13. The meteorological convention is used here, whereby the back-azimuth is the direction from which the disturbances originate.  Table 1 and displayed here. Maps were created with software Generic Mapping Tools (GMT v4.5.12; http://gmt.soest.hawaii.edu). Maps were created with software Generic Mapping Tools (GMT v4.5.12; http://gmt.soest.hawaii.edu).          Table 1 and displayed in Supplementary Fig. 1.      gradients. In addition, we find that there is a seasonal difference in the driving mechanisms. During the summer, baroclinic gradients enhance the cyclonic flow that is already established by tidal rectification and possibly wind-driven setup. During the winter, baroclinic and barotropic gradients are both small.

Introduction
Coastal currents are responsible for the transportation of water along the inner shelf. In Rhode Island Sound (RIS), the coastal current transports water in a cyclonic motion primarily parallel to local isobaths. This connects Vineyard Sound, Buzzards Bay, Narragansett Bay and Block Island Sound with the continental shelf. The summer intensification and winter dissipation of this current is thought to be a result of three general mechanisms: tidal rectification, density gradients, and sea-surface setup. Luo et al. (2013) and Liu (2015) numerically estimate that cyclonic circulation is doubled due to the development of a bottom cold pool during summer months. This 52 drives a baroclinic flow that adds to already present tidally rectified flow. Studies done on either side of RIS in Block Island Sound (Ullman and Codiga, 2004) and south of Martha's Vineyard (Fewings and Lentz, 2010), find that a wind driven sea-surface setup plays a dominant role at various times of the year in driving barotropic flow.
The goal of this study is to determine what the relative contributions of the tidally rectified to the depth-average coastal current in RIS are. To address this question, we will quantify coastal circulation in RIS, document temporal and spatial distributions and estimate the mechanisms that contribute to the sustained cyclonic circulation.
Since the 1970s, many studies have been conducted along the North American continental shelf, which are adjacent to our study area (e.g. Bumpus, 1973;Brown et al., 1985;Beardsley et al., 1985;Chapman et al., 1986;Beardsley and Winant, 1979).
The water along the middle and outer shelf, southeast of RIS, originates from the Scotian shelf (Beardsley and Winant, 1979). Large scale wind stress and heat fluxes over the continental shelf generate an alongshore pressure gradient (Beardsley et al., 1985). The alongshore pressure gradient drives persistent flow to the southwest. In contrast in the inner shelf we expect the circulation to vary significantly over seasonal times-scales due to stratification changes, sea-surface setup and wind-forcing.
Previous observational studies of RIS have focused on the description of seasonal variations of hydrography and to a limited extent of coastal jets. We define the coastal jets in this area as specific sub-regions of a coastal current, often 10's of km long. These studies have a scattered temporal and spatial distribution of observations. Ship-based measurements made over several days throughout the year indicate seasonal stratification develops in the summer and dissipates in the fall and 53 winter months (Shonting and Cook, 1970;Armstrong, 1998;Hicks and Campbell, 1953). The relatively strong current flows to the west at the mouth of Narragansett Bay. The current was observed with acoustic Doppler current profiler (ADCP) surveys .  attribute this development to lateral density differences that extend from central RIS to the shoreline, resulting from differential heating and mixing. Ullman and Codiga (2004), Edwards (2004) and Ullman and Cornillon (2001)  General studies of circulation along inner shelf regions suggest generation by tides, buoyant plumes, surface gravity waves, across-shelf wind stress, along-shelf wind stress and pressure gradients (Lentz and Fewings, 2012). Two recent numerical studies of RIS, have illustrated the importance of both tidal rectification and buoyancy driven flow in the RIS basin wide cyclonic circulation. Liu (2015) and Luo et al. (2013) find that in addition to tidal rectification, circulation is enhanced by density gradients developing primarily in the summer. They also find that surface currents are affected by seasonal wind forcing in the winter.
Despite the early hydrographic studies and recent modeling, the response of the whole water column to wind forcing, variability of sub-tidal currents across RIS and observational confirmation of tidal rectification, are not well measured and 54 lacking in this inner shelf region. From late 2009 to late 2011, multiple moored ADCP and hydrographic instruments were deployed in RIS over multiple months, providing an unprecedented view of the sub-mesoscale circulation and fluctuations. In addition, spatially dense hydrographic surveys were performed during multiple seasons to further illuminate the regional variability of the hydrographic environment.
The resulting time series and hydrographic distribution allows us to explore the variability and dynamics of the coastal current.
We intend to improve on the knowledge in RIS by analyzing hydrographic and current measurements over seasonal time periods. We estimate the depth-averaged

Regional Characteristics
RIS is an inner shelf environment, connecting surf zones with the deeper middle and outer shelf. It is the region where the shelf circulation adjusts to the presence of the coastal boundary conditions, in both the horizontal and vertical (Lentz, 1995). The inner shelf extends to a depth of 10's of meters depending on wind strength, waves and vertical stratification. Often the inner shelf is defined as the region where surface and bottom boundaries interact (Lentz and Fewings, 2012).
Assuming a typical eddy viscosity (A) of 0.01 m 2 s -1 , a characteristic Ekman boundary layer ( 2 = 32 / ) at the latitude of 41 o would extend 15 m. Much of RIS is less than double this depth, indicating that the surface and bottom boundary layers should overlap. Therefore, RIS is classified as an inner shelf region.
Bathymetry and the coastal configuration plays a part in governing circulation in RIS. RIS is located outside the mouth of Narragansett Bay and bounded to the east and west by Buzzards Bay and Block Island Sound, respectively. RIS makes a bowl shape with deeper bathymetry in the center and southern edges, with depths averaging between 30-40 m. In the center, there is a long trough running northeast-southwest reaching depths of up to 55 m. RIS is open to the continental shelf to the south. We next discuss the major known mechanisms driving circulation and hydrography.

Tidal Flows
The tidal velocities in RIS range from 0.1 to 1 m s -1 (Luo et al., 2013), with the largest velocities observed around Block Island (Edwards, 2004). Tides are primarily semidiurnal and corresponding tidal ellipses are oriented along the northwestsoutheast direction (Shonting and Cook, 1970). The M2 is the largest constituent , accounting for over 85% of the tidal constituent energy (Liu, 2015).
Modeling done by He and Wilkin (2006) found tides in RIS are co-oscillatory with the open ocean. In this area, ocean tides in deep water force tidal amplitudes onto the shelf. The resulting tidal amplitudes increase with distance away from the shelf break due to the onshore progression of the tidal wave (He and Wilkin, 2006;Moody et al., 1984). The resultant orientations of the major tidal axes are primarily perpendicular to the shelf break and do not vary substantially across RIS. For any constituent, the ensuing high and low tides occur approximately uniformly across the shallow sound.

Residual Flows and Forcing Conditions
Circulation longer than tidal periods (e.g. > 33 hours) occurs primarily parallel to the shore. Large vertical changes in bathymetry over a relatively short distance cause cross-shelf circulation to decrease toward the shore. Any induced cross-shelf exchange creates vertical velocities, i.e. upwelling and downwelling. Vertical velocities are key to transporting larvae, nutrients, sediment, pollutants, and oxygen in the inner shelf (Lentz and Fewings, 2012). As such, we focus on the local forcing conditions in RIS in an aim to understand the drivers that influence the residual circulation.

Winds
Winds influence circulation through direct shear stress, induced sea-surface height gradients, and through mixing. During the winter season, winds blow to the east-southeast averaging less than 12 m s -1 . However, strong variation in both magnitude and direction are observed as storms pass through the area reaching up to 25 m s -1 (e.g. . Summer winds are less variable, with average wind directions to the northeast and magnitudes averaging less than 7.5 m s -1 . At shorter time-scales, the New England inner shelf is susceptible to sea breezes, on the order of 2-5 m s -1 , created by differential heating of the land and sea (Fisher, 1960).
Winds in the area have an average eastward component during both the winter and summertime, commonly generating upwelling conditions in RIS . Although, wind magnitude and direction have been documented for the southern coast of New England, the effects on circulation have not been well reported in RIS.

Tidal Rectification
Tidally rectified flow is thought to be more important than wind influenced circulation for periods longer than synoptic scales (3-5 days) in RIS (Liu, 2015). The resulting residual flow can be produced by tidal rectification created by exchanging momentum from tidal frequency motions to lower frequency motions in areas with sloped bathymetry . Tidally rectified flow is usually along isobaths with shallow bathymetry to the right of flow. When looking at the shore from RIS, circulation is counterclockwise or down-shelf. The resulting tidal rectification is on the order of 3 cm s -1 and predominantly along the periphery of RIS (Luo et al., 2013). The rectified flow is cyclonic with negligible flow in the center of RIS (Luo et al., 2013).

Solar Insolation
In the summer, numerical models suggest RIS residual flow is intensified and is almost doubled along the periphery (Liu, 2015). The increase is thought to be created by the presence of a bottom thermal front. The front starts when solar insolation overcomes tidal mixing in warmer months, leading to a well-developed pycnocline (Rosenberger, 2001;Shonting and Cook, 1970). The horizontal variations in vertical mixing creates a tidal mixing front and induces strong lateral density gradients (e.g. Holt and Umlauf, 2008). A numerical model supports stronger mixing along the shallow edges of RIS, and a resulting summertime bottom thermal front 58 (Liu, 2015). In addition to numerical models, measurements near the mouth of Narragansett Bay  have confirmed the intensification of residual flow associated with increased stratification during summer months.

Non-Local Forces
Non-local forces are thought to be small but could have limited effect on residual flow in RIS. Large-scale pressure gradients along the continental shelf remotely forces flow from the northeast to the southwest (Luo et al., 2013). However, idealized studies such as Chapman et al. (1986), indicate that along-shore flow remains mostly near the shelf break, limiting the effect on our region of interest. It is unclear the magnitude of the non-local forcing in this area.

Data
In order to further assess the forcings and local circulation properties in RIS, we focused on five key measurements: temperature, salinity, water velocity, tidal height, and wind velocity. We used moored acoustic Doppler current profilers (ADCPs), thermistors, conductivity, temperature and depth (CTD) instruments, underway hydrographic surveys, tide gauges, and weather buoys to best characterize these properties in RIS.
Moored instruments, measuring water velocity and hydrography were placed in and around RIS. Stations of moored instruments are described with a three-letter abbreviation given in Table 1. All stations had a moored ADCP (Fig. 2). Most stations had a chain of thermistors or CTD sensors located within 100-200 meters of the ADCP location. Therefore, the location of the hydrographic measurements is also referred to as the three letter station abbreviations of the ADCPs. Four deployments 59 were undertaken, each lasting at least a month (Table 1). Finally, four dense underway hydrographic surveys were collected in RIS. Location of the hydrographic measurements are illustrated in Fig. 2.

Velocities
Observations of velocity were made using RDI broadband 300 kHz and 600 kHz ADCPs moored on the seafloor. The ADCP data consisted of four deployments

Temperature and Salinity
To measure hydrography, we used HOBO Water Temperature Pro instruments and HOBO Pendant Temperature thermistor chains as well as Seabird MicroCat CTD chains. Exact locations are shown in Fig.1 and the deployment times are listed in Table 1. Seven to nine thermistor instruments were evenly spaced along a chain with one pressure sensor below the surface float. These pressure sensors were used to estimate the depth of the thermistors during the deployments. The known distance between pressure sensors and thermistors was used to estimate the depths of each 60 thermistors throughout the deployment. Similarly, five or six CTDs were attached to chains near moored ADCPs. The duration of the thermistor and CTD deployments are listed in Table 1. In addition to measurements in the water column, bottom temperatures were recorded during the entire experiment, with thermistors attached to the bottom mounted ADCPs.
Higher resolution spatial measurements of hydrography were obtained with CTD surveys. Multiple vertical measurements are averaged to 1m depth intervals on September 22-24, 2009;December 7-8, 2009;March 9-11, 2010;and June 16-18, 2010. The locations are displayed in Fig. 2, and are the same for the four surveys.

Tidal Height
We use tidal range to compare variations of observations with the spring-neap tidal cycle. Our reference tidal range is from tidal height at Newport, RI from NOAA tide gauge 8452660 ( Fig. 1). This is the closest tide gauge to our study area. In addition, we use tide gauges located at Newport, RI; Montauk, NY; New London, CT and Woods Hole, MA to estimate sea-surface height gradients over RIS.
Sea level heights were recorded every 6 minutes during this experiment and measured from mean lower low water. Data is available from the start to the end of  our experiment. The tidal range, the difference between high and low tide, varied from 0.7 m to 1.5 m at the tide gauge in Newport, RI.

Wind Velocity
Wind data was collected at station BUZM3, a NOAA weather station, located in the northeast corner of RIS (

Results
Analysis of the moored ADCPs, moored hydrographic instruments and underway hydrographic surveys provides information at a variety of temporal and spatial scales in RIS. Due to the sampling rate of our instruments and the multiple year duration of the experiment, we analyze RIS data at time-scales from hours to seasons. The spatial resolution of velocities is limited to the seven locations of the moored instruments. These instruments are 5 to 10 km apart and therefore sample a wide range of environments in and around RIS. The hydrographic surveys provide more detailed observations of hydrography between the moored stations but limited temporal evolution. Below we discuss in detail our observations.

Tidal Circulation
Temporal and spatial variations within RIS, are first analyzed with a spectral analysis of the velocity time series data to determine the most prominent frequencies of interest. Data was linearly detrended and we use the Welch method, utilizing the 63 Fast Fourier transforms of the auto-correlated time series, to create a power spectrum . Overlap of the time series is set at 50% to reduce variance. Gaps in the data are filled with station averages . It is important to note that padding the data may bias some of the high frequency signal close to the sampling frequency. We produce power spectra of our depth-averaged velocities and depth-averaged shear (Fig. 4 & 5).
Strong power spectral density peaks are observed for depth-averaged flow in the diurnal (1 cpd In addition to depth-averaged velocities, we also explore depth-averaged shear. Depth-averaged shear is calculated by differencing the near-surface and near-bottom velocities and dividing by the separation distance of the two measurements. Power spectral density analysis of depth-averaged shear flow indicates strong peaks at the M2 and diurnal frequency (Fig. 5). This is true for both the along and across-shore directions. The M2 has the strongest peak at all stations. Some stations such as BIN, WBF, and MV7 have peaks at the M4 tidal frequency, but this is not consistent across all stations.
In addition to spectral analysis described above, we utilize the MATLAB toolbox t_tide to quantify the tidal components in our velocity measurements (Pawlowicz et al., 2002). T_tide provides estimates of the significant tidal constituents, the amplitude of the major and minor axis of the tidal ellipse, inclination and phase. The tidal ellipse is a way to describe the current vector for tidal velocities, with the current vector tip tracing out the path of the ellipse. Inclination and phase describe how and where the vector rotates around the ellipse. Specifically, the inclination is the direction of the semi-major axis and is measured counterclockwise from due east. Phase is the orientation of the vector relative to a specific time (Pawlowicz et al., 2002).
Depth-averaged analysis of tidal constituents shows the M2, the principal lunar semidiurnal constituent, is the largest constituent by a factor of 4 for all stations.
Results are listed in Table 2. M2 tidal ellipses are primarily oriented perpendicular to the shelf break, except for station BIN (Fig. 6). Magnitude of the major axis primarily increases as the depth of the stations decreases. Consistent with the theory that momentum must be conserved as the tidal wave approaches shore. The second two largest semidiurnal frequencies are the N2 and S2, the larger lunar elliptic and principal solar semidiurnal constituents, respectively. These constituents have comparable magnitudes to one another at each station ( Table 2). The two largest diurnal frequencies are the K1 and O1, both lunar diurnal constituents. The major axis of these diurnal frequencies is primarily oriented parallel to isobaths and the shoreline ( Fig. 7). Our observations support the analysis of He and Wilkin (2006) who modeled the barotropic amplitude of tidal constituents and found that the semi-diurnal constituents dominate this area over diurnal frequencies. Table 2 . Bottom ellipses tend to be about 10 % smaller in magnitude and more circular than the overlying layers. The reduction in tidal amplitudes is indicative of bottom friction (e.g. Edwards and Seim, 2008).
The bottom layer is also dissimilar in that the phase marker for the largest tidal component, the M2, (Fig. 7) is rotated relative to the overlying layers. All stations except BIN and BIE have bottom tidal ellipse phases that are rotated clockwise from the surface phase. This indicates the bottom layer velocities slightly precede the overlying layer velocities, as the rotation at these stations is clockwise.

Residual Circulation
Depth-and deployment-averaged horizontal velocities are largest closest to shore. Deployment-averaged moored ADCPs provide a mean current field for RIS ( Fig. 8). Central stations of CLC, CRS, SAK and WBF have residual velocities on the order of 1 cm s -1 , that are within the corresponding variance ellipses. Stronger flow was measured closer to shore on the order of 5 cm s -1 at stations BIN, BIE and MV7, primarily parallel to the local shoreline (Fig. 8).
Inflow into RIS, comes in around the east, near Martha's Vineyard as indicated by deployment-averaged velocities at station MV7 (Fig. 8). We plot depth-averaged velocities filtered with a 33-hour low-pass filter for all four deployments to explore shorter time variations ( Fig. 9-12). The strongest depthaverage velocities are observed around Block Island at stations BIE and BIN. Two major time-scale variations are apparent from the residual velocities ( Fig. 9-12). One is that there is a lot of variation across all stations at synoptic time-scales. However, velocities at stations are not necessarily coherent or in the same direction as the wind measured at BUZM3 during our experiment ( Fig. 9-12). The second time-scale is on the order of months, indicated by an intensification of velocity magnitudes at several stations during different seasons. For example, during deployment 2 ( Fig. 10), throughout summer months, the velocities at stations BIN and BIE are qualitatively larger than velocities measured at the same stations during the 1 st deployment (Fig. 9), i.e. during winter months.

68
Subtidal velocities are variable across stations. We quantify of how coherent the subtidal circulation is across RIS with correlation coefficients of the complex depth-averaged velocities. The complex depth-averaged velocity vectors lead to a complex correlation coefficient. The magnitude (Fig. 13a) is used to quantify the linear similarity of the two-time series. An angle (Fig. 13b) is calculated from the resulting complex value, representing the rotation of one station relative to another.
For all correlations we consider a method that utilizes effective degrees of freedom to determine statistical significance ( Fig. 13c) described by MacKenzie and Schiedek (2007). Using the Student's t-distribution we determine if our correlation is significant. If the computed t-value based on the correlation coefficient (C) is greater than a critical value of a Student's t-distribution for a probability of 95% with a specific effective degrees of freedom ( ), than the correlation is considered significant (Hald, 1976): Effective degrees of freedom are based on the autocorrelation of the time series of interest and accounts for the fact that filtering data reduces independence of samples.
Although, many of the correlations are significant between stations, (Fig. 13c) we find stations MV7 and SAK have limited significant correlations. MV7 is not significantly correlated with any station and SAK is not significantly correlated with stations BIE and BIN, around Block Island. In addition, correlation coefficients are relatively low. Only correlations between stations WBF and CLC are greater than 0.6 ( Fig. 13a). Low and insignificant correlations between stations, suggests smaller scale forcing is important to circulation over our stations at synoptic time-scales. We explore longer time periods, starting with seasonal averages, in hopes of finding larger scale patterns consistent over all stations.
Seasonal averages indicate stronger flow along the periphery of RIS during summer. We divide data into two major time intervals (Fig. 14). Periods during well- Mean monthly depth-averaged flow indicates a persistent negative circulation in the along-shore direction at all stations ( Fig. 15b &e). Along-shore flow is maximized during June and July of 2010 and minimized during November and December for most stations (Fig. 15b & e). We discuss the possible mechanisms for the along-shore flow in section 5.

Regional Hydrography
The four CTD surveys made during the experiment illustrate the seasonal changes in density and stratification in RIS. Depth-averaged density increases from north to south during all four deployments, mostly a function of the increasing depth Temperature and salinity measurements from CTD deployments indicate that density distribution is dependent on both properties. During the September, December and March surveys, temperature is fairly homogenous (Fig. 18a, b, c), indicating that 71 salinity is controlling the density differences observed. In December 2009, interestingly the survey measures a warm salty bottom (Fig. 19b), with a temperature inversion (Fig. 18b). Ullman et al. (2014) have documented the temperature inversion in this area as a warm salty intrusion from the shelf break. The survey in June 2010 indicates a spatial salinity gradient (Fig. 19d) that coincides with spatial temperature gradients ( Fig. 18 d), creating the largest spatial density gradients. Strongest density differences both vertically and horizontally are measured during this survey.
We quantify the vertical distribution of density with stratification measurements. Stratification can be measured as a potential energy anomaly ( ) of a given water column.  describe potential energy anomalies as the amount of work required to completely mix the water column: where H is the water depth, ρ is the density, g is 9.8 m s -2 and z is depth below sea level. When CTD data are available we can directly estimate the density at various depths and calculate ϕ.

Wind and Freshwater Influences
Before the discussion of the depth-averaged momentum balance in section 4.5, we explore two regional environmental parameters that impact the momentum balance. Winds and freshwater can directly influence the inner shelf dynamics through wind stress, sea-surface setup and baroclinic effects. We describe the first order trends of monthly averages of wind and river runoff during our experiment.
The seasonal cycle of monthly mean wind has peak magnitudes during November through January ( The seasonal cycle, in monthly mean freshwater river inputs into RIS, can be qualitatively viewed from the transport of two major rivers in the area. Transport of water masses comes into RIS from Long Island Sound, Narragansett Bay, Buzzards Bay and Vineyard Sound. The largest transport, by a factor of 10, is the transport from Long Island Sound followed by Narragansett Bay . We use the Connecticut river, the major source of freshwater into Long Island Sound and the Blackstone river, a major source for Narragansett Bay to illustrate the seasonal variation in freshwater input into RIS (Fig. 22). River transport into Buzzards Bay and Vineyard Sound is tiny and therefore not highlighted in this study.
Both the Blackstone and the Connecticut rivers have a peak in river discharge during February and March (Fig. 22). Although, the Blackstone is the largest of several tributaries into Narragansett Bay, the magnitude is much smaller than the Connecticut river. The size discrepancy illustrates that freshwater influence of Narragansett Bay is of secondary importance compared to freshwater coming into RIS from Long Island Sound. Similar to Long Island Sound and Narragansett Bay, the freshwater input into RIS likely increases during spring. Studies of Long Island Sound suggest that much of the freshwater exits Long Island Sound in-between Block Island and Montauk Point, not entering RIS (Edwards, 2004). Therefore, we do not expect freshwater plumes to play the largest role in governing circulation through the entire RIS.

Depth-averaged Momentum Balance
The monthly depth-averaged flow provides insight into what is driving circulation in RIS. The relationship between flow and forcing conditions is observed through analysis of the depth-averaged momentum balance. The resulting momentum balance in each direction is: where (Y, ̅ ) are the residual depth-averaged velocities, ( , ) are the tidal velocities, f is the Coriolis parameter, _ is a reference density, is pressure, H is the depth, ( ] , ] l ) are the surface wind stresses, ( â , a n ) are the bottom friction stresses and 〈 〉 denotes averaging over a tidal cycle. The o (eq. 3) direction is across-shore and the o (eq. 4) direction is along-shore. The first two terms in the above equations are the local acceleration and the Coriolis term. These terms can be estimated directly for the moored ADCP velocity measurements. The next four terms of the eq. 3 & 4 describe the wind stress, bottom stress terms, advective tidal stresses and the pressure gradient respectively (Ullman and Codiga, 2004;Visser et al., 1990). The advective tidal stress describes the momentum exchange from tidal frequencies to subtidal frequencies in areas where tidal ellipses change shape.
Monthly averages of the momentum terms calculated from deployment velocities and wind speed are displayed in Fig. 23. These include the acceleration term, the Coriolis term, wind stress and bottom stress. The magnitude illustrates the terms that are of primary importance in the momentum equation.

Local Acceleration and Coriolis terms
We calculated the local acceleration from a forward difference method of 33hour low pass filtered velocities. In both the across and along-shore directions, the acceleration term has small magnitudes measured over all stations, relative to other momentum terms.
The Coriolis term, found by multiplying the observed depth-average velocities by the Coriolis parameter, has magnitudes in both eq. 3 and 4 on the order of 5 x 10 -6 m s -2 . Monthly averages of the Coriolis term indicate there is always a negative acrossshore ( ) term measured at all stations. This reveals a cyclonic flow around RIS. The monthly time series shows a maximum magnitude near June and July and a minimum across stations during December and January (Fig. 23). Station BIN does not fit this pattern as well as the other stations. The along-shore monthly-averaged time-series of eq. 3 shows the Coriolis terms also has a large magnitude but smaller than the acrossshore direction. Stations in RIS (BIE, CLC, CRS, MV7, SAK and WBF) have values close to zero while the station north of Block Island (BIN) is positive.

Wind Stress
The wind stress can be estimated using the depth of each station and drag coefficient estimated from the wind speed. The wind stress is calculated using where is the wind vector at 10 m above sea level and is the density of air. We use a drag coefficient as defined by Garratt (1977) Station BUZM3 measures wind at 24.8 m above sea level. We calculate the U10 by using a logarithmic relationship: where ( ) is the wind velocity vector measured at height z, is the Von Karman constant and U is the surface roughness length scale. To calculate the wind stress, we start with an initial guess in eq. 6 to calculate Cd, found from approximating U10 with U25. We calculate z0 at 10 m height, where U = (10 )/ ƒ/3A ‚ , found from rearranging eq. 8. Next, we solve for a new U10 by dividing eq. 8 solved at z = 10 m, by eq. 8 solved at z = 25 m, providing the relationship between U(10 m) and U(25 m): 77 Finally, we update Cd with the new U10 velocity. We iterate this process until a solution for U10 and Cd converge, and changes are less than 0.01 %, usually 5 iterations. The wind stress is calculated instantaneously and averaged over a month.
The wind stress terms have small magnitudes, less than 1 x 10 -6 m s -2 (Fig. 23), when averaged over a month. This is true for both the along-and across-shore directions. Magnitudes do increase during winter months when wind speed is largest.
However, the wind stress is not large enough to balance the observed Coriolis term at any station over monthly time-scales.

Bottom Stress
The bottom stress is highly uncertain due to poorly estimated drag coefficients.
We, therefore, use two estimates to calculate bottom friction and provide a range of possible estimates. First, we use a quadratic bottom stress, defined as: where is a near bottom velocity vector. Bottom velocities are approximated with the deepest measurements recorded at our ADCPs. r is set to 2.5 x 10 -3 , a commonly used value along the continental shelf (e.g., Hofmeister et al., 2009;Lund-Hansen et al., 1996;Edwards, 2004;Pu et al., 2015;Sylaios et al., 2013).
The second method we use is the linear drag bottom stress, defined as: where " is the depth-averaged velocity and we use a value of 4 x 10 -4 m s -1 for , similar to Liu and Weisberg (2005). Both methods calculate bottom stress instantaneously and values are averaged monthly.

78
Monthly-averaged bottom stress term in eq. 3 &4 is small on the order of 1 x 10 -6 m s -2 for most stations, around the same order of magnitude as the wind stress.
The one exception to this is station BIN. This station has positive bottom stress on the order of 2 x 10 -6 m s -2 throughout the year in both the across-and along-shore directions. In Fig. 23, we only display the quadratic estimate of bottom stress as there was little difference between the quadratic and linear bottom stress estimates.

Advective Tidal Stress
We use the Regional Ocean Modeling System (ROMS) to estimate advective The numerical experiment is verified against the measured tidal components at each moored station. Fig. 25 displays the depth-average tidal ellipse at each of the seven moored ADCPs. In general, the ROMS numerical model predicted tidal ellipses 80 that have a slightly larger eccentricity but the same general magnitude, rotation direction and phase as the ADCP measurements. The one exception is station BIN. At this station, the amplitude predicted by the ROMS model is comparable to observations, but the orientation and phase are off by about 20 degrees (Fig. 25). The rotation at this station is also not well modeled, as ROMS predicts a counterclockwise rotation instead of the observed clockwise rotation. Despite the discrepancy at station BIN, the rest of the stations agree well with observations and we precede to use the model output to estimate the advective tidal stress.
Studies such as Ullman and Codiga (2004) and Visser et al. (1990) have calculated the advective tidal stress using depth-average tidal velocities. We improve upon this estimate by depth-averaging the depth-dependent advective tidal stress:

Advective Tidal Stress Sensitivity
We test the sensitivity of the tidal advective stress calculation by changing zo The change in magnitude of the tidal stress is found to be -19%, over the RIS domain.
The negative value indicates the magnitude of the tidal advective stress is smaller for the experiment with _ = 0.01 .
Reducing _ 0.0001 , we find ∆ and ∆| | to be 18 % and 13%, respectively. With decreased bottom roughness the tidal advective stress magnitude increases slightly. Therefore, we find that changing the bottom roughness by orders of magnitude from 0.01 to 0.0001 m, changes the tidal advective stress by + 20 %.

Pressure Gradient
The last term in the depth-averaged momentum equation changes the momentum balance through the addition or reduction of mass and changes in density distributions. Particularly, pressure gradients can create depth-average momentum force, through changes in volume (steric) or through changes in mass (non-steric). We use the method described by Ullman and Codiga (2004), who separate the pressure gradient into a steric and non-steric contributions, which we will call baroclinic and barotropic respectively:

83
where is defined as = _ [1 + ( , , , )], _ is the maximum density, and ¸ is a reference depth. The steric ( ] ) and non-steric ( >] ) sea-surface heights are defined as: ¸ is defined as 52 m, the deepest survey position in our study. The last term in eq.
14 and 15 are found by integrating along the bottom from the deepest point (−¸) to the bottom of the survey (− ) (Csanady, 1979). We use the four hydrographic CTD surveys to estimate the steric contributions of the pressure gradients for the month the survey was taken.
If we assume only a geostrophic balance, we can calculate the depth-mean velocity induced by baroclinic pressure gradients. The results of these calculations are displayed in Fig. 28 for all four CTD surveys assuming no flow normal to the bottom.

84
The baroclinic pressure gradient is interpolated to each moored station to estimate the contribution to the depth-averaged momentum balance. We apply the baroclinic pressure gradients to the months the survey deployments was taken in. In the across-shore direction, we find that most of the estimates of the baroclinic pressure Sound from the Connecticut river. Negative gradients indicate there is dense water down-shelf. This is the case for station BIE, which indicates denser water is located closer to the middle shelf and less dense water located closer to the mouth of Narragansett Bay.
The last component of the momentum balance is the barotropic pressure gradient. Estimations were made with tidal height measurements from NOAA tide gauges surrounding our study area: Newport, RI; Montauk, NY; New London, CT and Woods Hole, MA (Fig. 29a). Sea surface height (SSH) relative to NAVD88 was used without correcting for the inverse barometer effect. Collocated with the tide gauges 85 were atmospheric pressure measurements available from NOAA. There are many gaps longer than several days in the pressure measurements, especially during the summer of 2010. The advantage to not using the inverse barometer correction is that the SSH can be used to find the total sea surface set up and does not limit the amount of SSH data available.
A 4 th order Butterworth 33-hour low pass filter was used to remove the tidal signal of the SSH and calculate the difference between stations. The sea surface slope was calculated using a linear least-squares fit to the four tide gauge stations SSH and results are shown in Fig. 29b.
Comparing the monthly averaged SSH gradient across the four tide gauges with wind indicates a similarity between the two data sets (Fig. 30). As wind changes direction from southeastward blowing in winter to northeastward blowing in summer ( Fig. 30 b), the SSH gradient changes from increasing heights to the south (positive slope) during winter periods to increasing heights to the north (negative slope) during the summer (Fig. 29b). The change in SSH gradients shows a seasonal variation contemporaneous with wind direction changes across RIS.
The SSH gradient across tide gauges was used to calculate the barotropic gradients at each station. Specifically, the SSH gradient (∇ M ) is the combination of baroclinic and barotropic gradients: where ∇ M is the horizontal gradient. ∇ M ] is calculated using spatial distribution of ] found from the CTD surveys. Due to the fact that we only have a relative ] we subtract the averages from both and ] to find the relative >] shown in Fig. 31. magnitudes as large as the Coriolis, baroclinic and barotropic terms (Fig. 27). It is likely that the barotropic term is over estimated. We have approximated a linear seasurface setup created by wind. Most likely the SSH gradient has steeper slopes closest to land. It is likely we have overestimated the barotropic contribution at each station, as our stations are in deeper water far from the shoreline.
In the along-shore directions residual estimations have the same order of magnitude as other terms in eq. 4. However, there are both positive and negative residuals, unlike the across-shore residuals estimates (Fig. 27 h). It should be noted that residuals in the along-shore direction all have the same sign as the barotropic estimates, except BIN. This again is consistent with the barotropic term being over estimated.

Discussion
The coastal boundary causes a restriction in cross-shelf flow and forces subtidal, depth-averaged flows to be stronger in the along-shore direction (Lentz and Fewings, 2012). The sloping bathymetry as well as the geographic shape governs how and where circulation can occur. We discuss the impacts the three-dimensional shape of RIS has on circulation, in the following sections.

Tidal Rectification
Bathymetry directly influences the depth-averaged dynamics through tidal rectification. In order for tidal rectification to occur, the tidal excursions must be on the same order of magnitude as the bathymetric features. Tidal excursion is the total distance a parcel of water travels from flood to ebb tide. The distance captures the horizontal Lagrangian motion (Parsa and Shahidi, 2010). Tidal rectification studies characterize the topographic length scale as ( º ‚» ‚¢ ), where h is the water depth and rº r^ is the slope of the bathymetry Loder, 1980).
We calculate the tidal excursion with the M2, the largest constituent across RIS. Defined here as l=2V/ω, where V is the amplitude of the M2 major axis of the depth-averaged tidal velocity and ω is the M2 angular frequency. This results in excursions ranging from 1.1 km at SAK up to 5.5 km at BIN ( Table 3). The average excursion over the 7 stations is 2 km, slightly smaller than the tidal excursions of Narragansett Bay . We calculate the topographic length scale near our moored stations by using the depth and local slope. The local slope is found by finding the change in height over 1 km in the across isobath direction, centered at the station of interest (  Table 4: Topographic slope and length scale calculated from the depth of stations and bathymetric slope. Liu (2015) modeled the effects of tidal rectification in RIS also using ROMS and only the M2 tides. They found a persistent tidally rectified flow over most of RIS with stronger velocities towards the shore. Compared with our ROMS experiment, both Liu (2015) and our model predict tidal rectification generates a cyclonic flow on the order of 1 cm s -1 , with higher velocities measured in shallower waters, as well as a clockwise flow around Block Island (Fig. 32). Predictions of tidally induced flow from our numerical model should be taken with caution as the boundaries of our model are close to the domain of interest. However, the fact that our ROMS results are similar to Liu (2015), provides confidence that our estimates are reasonable.
Our model deviates from Liu (2015) in that our coastal current is narrower on the order of 5-10 km wide instead of around 20 km. This agrees with our analysis of tidal excurison and topographic length scale, which suggests only BIN and MV7 should be influenced by tidal rectification. We predict stronger residual flow on the order of 1-2 cm s -1 at the mouth of Narragansett Bay, not present in the Liu (2015) model (Fig. 32). In addition, our model predicts two headland eddies north of Block Island and a horizontal detachment of the coastal current from the eastern side of Block Island also found by Sun et al., (2016). The differences between our and Liu (2015) ROMS models are thought to be a result of resolution, since tidal rectification is strongly dependent on the bathymetric gradients. Liu (2015)

Pressure Gradients
Bathymetry also affects pressure gradients through mixing. With increased velocities and decreased water depths, more mixing can occur. Mixing specifically effects baroclinic gradients by changing the horizontal density gradients. Therefore, mixing influences the pressure gradients in the momentum term. Often the boundary of the tidal mixing front is defined with the Simpson-Hunter criterion (h/U 3 ). The water depth (h) and the tidal velocities (U) are used define the position of tidal mixing fronts and the degree of mixing in shallow coastal waters (Simpson and Hunter, 1974). Liu (2015) estimated this parameter and did not expect a tidal mixing front to be extended to the surface in RIS. Their numerical analysis however, did reveal a bottom thermal front.

Presence and Modification of the Bottom Thermal Front
Bottom thermal fronts in RIS have been modeled to estimate the influence of buoyancy driven flow in RIS driven by mixing due to the M2 tidal constituent. Liu RIS is subject to not only the M2 but other tidal constituents that create variable tidal mixing over fortnightly periods. The fortnightly intensification of tidal mixing fronts has been theorized (Dong et al., 2015;Sharples, 2007) but few studies have documented this phenomenon. We look for direct evidence for the enhancement of the coastal current due to either tidal rectification or tidal mixing within our dataset.
Over many tidal periods, we expect tidal rectified flow to have a spring-neap cycle as the M2 is not the only tidal constituent present in RIS. To test this variation, we quantify the correlation of the along and across-shore velocities with tidal range measured at Newport, RI. Particularly, we expect the along-isobath at BIN and MV7 flow to be correlated with the spring-neap cycle as tidal.
We use a 7 day to 3-month bandpass to filter the tidal range and depthaveraged velocities to identify the spring-neap signal. Filtering aids in the removal of the semidiurnal and diurnal tidal frequencies as well as synoptic ( Correlations are performed by using the MATLAB function xcorr. Significance is determined with a Student's t-distribution, comparing a calculated tvalue (eq. 1) with a critical value found by using the effective degrees of freedom and a 95% confidence level.
For this correlation analysis, we use all depth-averaged velocity data available from November 2009 to December 2011. All but the two deepest stations (WBF and CLC) have significant correlations with tidal range in the along-shore direction ( Table   5). Stations SAK, BIE, MV7 and CRS all have negative correlations, indicating an increase in cyclonic flow in the negative o direction, during spring tides. For these stations, we found that our numerical estimates of the advective tidal stress are small.
Therefore, the correlation is hypothesized to be a result of baroclinic pressure gradients changing as tidal mixing increases during spring tides. Table 5: Correlation coefficients and lag between tidal range measured at Newport, RI and velocity measurements. Velocity measurements are broken up into along and across-shore depth-average velocities and velocity shear. Data is used from all available periods and bandpass filtered from 7 days to 3 months. Significant correlations marked with * and lag (days) of velocity behind tidal range is indicated in parenthesis. Station BIN has depth-averaged along-and across-shore velocities, positively correlated with tidal range, resulting from onshore flow being minimized during spring tide and maximized during neap tide. As BIN is likely located in an area influenced by a headland eddy, we did not necessarily predict a negative correlation with spring neap tidal cycle. It is unclear if this change in velocity is consistent with a headland eddy as we would expect velocity to increase during spring tide. However, if migration of the eddy occurs as tidal velocities increase, there could be a reduction of velocity over the station as the eddy moves further offshore. This station has large advective tidal stress (|1|x10 -5 ms -2 ) and the change in velocities is likely a combination of tidal rectification and tidal mixing.

Depth-averaged Velocity Depth-averaged Shear
The increase in tidal velocities creates more mixing in shallow areas and we look for evidence of this in our hydrographic data. The horizontal resolution of hydrographic measurements is limited in our study and we therefore, look for temporal changes at moored stations. We use the moored hydrographic measurements collocated at the ADCP to examine fortnightly trends and potential mixing. Three days after spring tide, local minimums in potential energy anomalies (eq. 2) are observed at BIN, WBF, SAK and MV7 (Fig. 33d). This is a reduction of stratification after the spring tide.

Wind Driven Sea-Surface Setup
Nearby studies of depth-averaged momentum have documented along-shore flow created by wind generated sea-surface setup, as well as baroclinic gradients. To the west of our study area, at the outflow of Long Island Sound, a coastal jet was found to be primarily in geostrophic balance with an across-shore pressure gradient (Ullman and Codiga, 2004). They define the coastal jet as a strong persistent circulation located close to shore. This pressure gradient, which was intensified during summer months, resulted from buoyancy driven flow. During winter months, the coastal jet was reduced as a result of sea-surface setup by upwelling-favorable winds, i.e. the surface tilt created a pressure gradient that opposed the jet flow.
East of our study area, off the southern coast of Martha's Vineyard, a similar depth-averaged momentum balance was performed by Fewings and Lentz (2010).
These authors found local wind generated sea-surface tilting was an important mechanism along the inner shelf. Upwelling winds in this area, hindered the alongshore westward flow and downing-welling winds enhanced the along-shore westward flow. It is likely, given the influence of pressure gradients in nearby areas, that both buoyancy and sea-surface tilt play an important role in governing circulation in RIS.
We explore the effect of wind generated sea-surface setup by quantifying the correlation between wind and velocity. All available data, both depth-averaged velocities and depth-averaged shear measurements, is bandpass filtered at 33-hours to 3 months, exploring the synoptic timescale to compare to wind. At each station, for 96 both the across-and along-shore directions, we calculate the correlation coefficient and lag time for 5 o increments in wind directions. In Table 6 Although, wind stress and sea-surface tilt effect all stations, each station may be more sensitive to one of these forcings. We suspect that wind blowing to the northeast forcing water directly at stations BIN, BIE and SAK through shear stress, causes flow in the same direction as the wind (Table 6), has a larger impact on depthaveraged momentum at this synoptic timescale. This hypothesis is supported by high correlations for winds blowing in the southwest direction (Table 6). Along the eastern side of RIS, northwestward blowing winds create a pressure gradient through tilting sea-surface that induces flow in the negative o direction, that is dominant over wind shear at these timescales. This hypothesis is supported by high correlations for winds blowing in the northeast direction for stations CLC, CRS, MV7 and WBF (Table 6). Wind directions found for maximum correlations with the across-shore direction have a greater variation than the along-shore estimates. Strong variations in wind direction were also found for both the along-and across-shore depth-averaged shear comparison with wind.

97
We further postulate on the relationship between wind and sea level during the summer season. Longer than synoptic timescales, we have found wind stress to be small in the depth-averaged momentum equation (Fig. 23). Fig. 31 displays the barotropic estimation of sea-surface setup as well as the monthly average wind magnitude and direction. Throughout all seasons the wind is parallel to the coastline in this area and should create upwelling conditions. For example, the SSH are low near shore for all estimations except March 2010 (Fig. 31). At our moored stations this creates a negative across-shore barotropic gradient that indicates sea level increases offshore (Fig. 27e). However, RIS is more geographically complicated than a northeast trending coastline and is bounded to the east with a bend at the northeast corner. We hypothesize the bend on the eastern side of RIS retains water, creating a sea-surface slope that enhances the coastal current and reduces the barotropic across-shore term during the summer. The June 2010 barotropic estimate (Fig. 31d) illustrates the northeastward blowing average summertime wind. The June estimate shows a SSH high along the eastern side of RIS and a reduced southward gradient. This is illustrated through more north-south trending contours of constant SSHs. Further observations are required to verify this hypothesis.

Seasonal Setup
To review all the contributions to the cyclonic coastal current, we summarize the relative size of the momentum terms averaged over stratified (summer) and wellmixed (winter) periods in various regions of RIS. We calculate momentum term averages during stratified periods, by taking the mean magnitude of the terms calculated by eq. 3 and 4. Well-mixed estimates were averaged over November through March. We characterize the magnitude of the momentum terms into three categories: strong (>|0.5 x10 -5 m s -2 |), moderate (|0.1-0.5 x10 -5 m s -2 |) and weak (<|0.1 x10 -5 m s -2 |). The pressure term is divided into the baroclinic contribution, estimated from the CTD deployment in summer and winter months, and barotropic estimates.
Results for each term in various regions of RIS are listed in Table 7 and shown in Fig. 34. In regions 1, 2 and 3, seasonal changes in the barotropic component of pressure gradient was observed through the use of tide gauges. When the monthly average wind was blowing to the southeast (Fig. 31 b), the wind-driven sea-surface gradients impose a force that inhibits the cyclonic flow. Similarly,  found that the coastal current was minimal or non-existent during winter seasons, also supporting our theory that the wind set-up of the sea-surface inhibits the coastal current. During stratified months (Fig. 31 d), the sea-surface slope, which normally increases to the south, is reduced. This reduces the effect of the barotropic force on the cyclonic coastal current.
The last area considered is region 4, which is characterized by weak momentum terms. The three moored instruments in this area (Fig. 34) showed the least 101 variation throughout the year. The flow in region 4 is sluggish and weak at seasonal timescales (Table 7).

Conclusion
Circulation in RIS is dominated by tides. Tidal velocities range from 0.       Frequencies that are not significant are not plotted. The S2 frequency was not significant for WBF surface velocites. The O1 frequency was not significant at stations BIN and SAK. The O1 frequncy is also not significant for BIN, MV7 and WBF bottom velocities and the K1 frequency is not significant for BIE bottom velocities.                           To be submitted to Journal of Marine Systems

Abstract
Since 1999, seasonal time series of temperature and salinity measurements have been collected from multiple buoys in Narragansett Bay. These in situ data have provided hydrographic information about the local water column and stratification. In this study, we utilize these data along with numerical modeling to assess the influence of particular environmental factors on stratification. Specifically, we find that tidal height directly correlates with stratification changes at timescales on the order of hours. We set up an idealized simulation with constant freshwater input and no surface fluxes to determine the relative contributions of advection, straining, and differential diffusion from tidal velocities on stratification. We find that, due to the local geometry and bathymetry of Narragansett Bay, the stratification varies both temporally and spatially. Temporal changes are controlled by advection and straining. These common stratification variations can further be characterized into two distinct advection-driven regimes. The location of these regimes depends on the local stratification gradient.
The first regime is characterized by advection of more highly stratified water from the north, which produces a maximum in stratification occurring during slack low tide.
The second regime results from advection bringing more stratified water from the south into relatively shallow, well-mixed areas. This regime is characterized by the occurrence of a stratification maximum preceding slack high tide.

Introduction
Water column stratification plays a major role in biological, chemical and physical processes within estuaries. The degree of stratification regulates the vertical exchange within the water column. Weaker stratification allows for increased vertical 132 mixing, which leads to more uniform distribution of nutrients, dissolved oxygen and plankton. The vertical exchange between oxygen-depleted deep water and the oxygenrich surface is reduced as a result of strong stratification, potentially leading to deleterious conditions such as hypoxia and anoxia (Lin et al., 2006). For estuaries, such as Narragansett Bay (NB), stratification has been noted to play an essential role in the regulation of hypoxic events created by algal blooms (Deacutis et al., 2006).
Stratification patterns in estuaries vary over a range of temporal scales.
Previous studies indicate that changes in stratification, defined as the vertical difference in density, can be affected by surface heat fluxes (Lund-Hansen et al., 1996), precipitation and evaporation (Nahas et al., 2005), mixing (Burchard and Hofmeister, 2008), straining (Rippeth et al., 2001), and advection . Timescales for changes in stratification range from less than an hour, due to mixing and turbulence, to over a year, due to seasonal changes in freshwater input and solar heating.
NB is a unique estuary to study changes in stratification because of its low freshwater input and complex geometry. We use observations and numerical models to answer three questions. First, when and how much does stratification change in NB? The available time series, over 10 years of hydrographic data, allows for the direct characterization of spatial and temporal stratification patterns. The second is why does stratification change? We employ the statistical techniques of spectral and coherence analysis on in situ stratification measurements to determine which environmental parameters contribute to the patterns observed in NB. This analysis focuses on the effects that environmental factors such as solar radiation, river flow, 133 tides, and wind have on stratification. In addition, we use numerical model simulations to investigate a subset of these external forcing factors in greater spatial and temporal detail than can be accomplished using direct observations. Our last question is how does the stratification change compared to that in other estuaries? We relate our analysis of stratification to larger scale variations in estuarine processes.

Narragansett Bay Background
NB is a partially mixed estuary classified geologically as a drowned river estuary (McMaster, 1960). The bay is 40 km long and 15 km wide with an overall north-south orientation. The southern limit of the bay is delimited by Rhode Island Sound and the northern reaches are characterized by the inflow from the Taunton and Blackstone rivers. The bay is littered with embayments and islands resulting in an average channel width of around 5 km. Despite the relatively small channel widths, NB has a Kelvin number of order 1, meaning rotational effects are often important when describing fluid flow (Codiga, 2012;Pfeiffer-Herbert et al., 2015). The average depth is relatively shallow at approximately 10 m relative to mean sea level (e.g. . The deepest part reaches up to 50 m in the lower East Passage. Studies of topographic effects in estuaries are wanting, as many numerical models simplify geometry and bottom bathymetry. We add to the understudied area of bathymetrically complex estuaries by examining stratification in NB.  2016). Precipitation is a source of freshwater applied to the surface of NB that increases stratification. In the summer, evaporation exceeds precipitation due to the increase in solar insolation reducing stratification (Fan and Brown, 2003;Pilson, 2008). Surface heat flux also changes the thermal structure of the water column, tending to increase stratification when heat is applied to surface waters.

Sources of Stratification
Another component of the climate in Rhode Island is the seasonal change in wind patterns. During the summer, winds are predominately to the northeast (e.g., . During the winter, wind direction is more variable, but primarily blows to the east or southeast. Synoptic storms are more common in the winter and blow from the northeast (Kincaid et al., 2008). In addition to direction, winds vary in magnitude. Stronger winds occur during the winter, averaging 10 m s -1 .
Synoptic weather patterns in the winter tend to be even larger with wind magnitudes up to 25 m s -1 (e.g.  and often lasting 2 to 3 days (e.g. Weisberg and Sturges, 1976). Conversely, the summer season is often characterized by weaker winds, averaging 5 m s -1 (e.g. Pilson, 2008). In addition to the seasonal variation in wind magnitude, the summertime winds also exhibit a diurnal periodic variation, known as a sea breeze. The sea breeze effect, producing landward (northward) winds during afternoon hours, is created by diurnal fluctuations in the temperature difference between the land and water (Spaulding and White, 1990).
Not only do winds throughout the year play a role in vertical mixing, but they also induce changes in flow within the bay. Pressure gradients caused by along-estuary winds (axial winds) are an important driver in estuarine circulation, as the winds force flow into or out of an estuary (Geyer, 1997;Li, 2011, 2012;Scully et al., 2005).
The axial winds create a buildup of water at the head or mouth of estuaries, creating a pressure gradient resulting from the tilting of the sea-surface. The pressure gradient drives a return flow in the opposite direction of the wind, along the bottom. The resulting circulation induces straining, which changes stratification through the interaction of velocity shear and horizontal density gradients. Often, down-estuary winds create an increase in stratification. The resulting pressure gradient increases estuarine two-layer flow and drives less saline near-surface water over more saline water at depth. Up-estuary winds tend to reduce stratification as pressure gradients oppose estuarine two-layer flow.
Another driver of stratification is the rivers that bring freshwater into NB, primarily from the north. The average freshwater input is estimated to be 93 m 3 s -1 (Pilson, 1985). and Delaware Bay (570 m 3 s -1 ) (Gay and O'Donnell, 2009;Whitney and Garvine, 2006).
River discharge at the head of NB creates a spatial density gradient, with density increasing towards the south (Pilson, 1985). Vertical salinity differences account for approximately 80% of the stratification in the bay and are directly dependent on river discharge (Codiga, 2012). NB usually has some amount of stratification throughout the year, classifying it as a partially mixed estuary.
The highest river discharge is in the spring after seasonal snow melts (Spaulding and Swanson, 2008). Synoptic storms also contribute to the periodic nature of the river flow on shorter time-scales. Storm runoff generates a large peak in river discharge followed by a slower decline in river flux for several days.

Mechanisms of Changing Stratification
Tides are perpetually a driving force acting on NB. The characteristic tidal range at Newport, RI near the mouth of the bay is 1 meter and increases to 1.3 meters at the head of the bay in Providence, RI (Tidal Current Tables 2016-Atlantic Coast of North America, 2016). Hicks (1959) found that the largest tidal component in NB is the M2, the principal lunar tidal constituent, followed by the M4 and M6. The M2 tidal flow interacting with bottom friction and momentum advection produces the overtides at higher frequencies, i.e. M4 and M6.
The observed tides and tidal velocities are representative of a standing wave, with high and low water approximately corresponding to slack velocities (Hicks, 1959). The tidal ellipses, a way of representing the magnitude and direction of tidal flow, are elongated and oriented in a north-south direction. This indicates the influence of forced flow along the axis of NB. Tidal ellipses become more circular with depth, with reduced velocity magnitude due to bottom friction. Changes in tidal 137 current speed with depth creates a velocity shear during the tidal cycle. This is important as it results in tidal straining.
On time-scales longer than the tides, the general circulation, often called the residual flow, affects stratification through fluid motion. Residual flow is driven by the buoyancy and Coriolis forces in NB (Rogers, 2008). The horizontal density difference, created by the influx of fresh water to the north and salt water to the south, drives a gravitational circulation common in many partially mixed estuaries (Weisberg and Sturges, 1976). Often described as a two-layer system, the fresher upper layer moves seaward to the south, while the saltier lower layer moves landward to the north.
Tidal velocities are not strong enough to completely mix away vertical density differences and the gravitational flow acts to stabilize the pycnocline. Average residual velocities in NB are on the order of 10 cm s -1 .
The gravitational flow observed in NB is altered by the rotation of the Earth, the Coriolis force. Rotation plays a role in fluid motion in NB as the widths of the channels are on the same order of magnitude as the internal radius of deformation, the internal Rossby radius (2-4 km) . The result of the Coriolis force is a deflection of water to the right in the northern hemisphere. This causes a general counter-clockwise flow around NB and the large islands. A result of this counterclockwise flow is the deflection of isopyncals horizontally. Spatially, the isopyncals trend northeast to southwest, causing the East Passage to be denser than the West Passage (Hicks, 1959).
The combination of tidal and residual flow results in circulation that changes stratification. Flow changes stratification through advection, straining and mixing.
Two end-members of tidally driven stratification change are pure advection and pure straining. In areas where stratification gradients are large over the length of water parcel tidal excursions, advection primarily modifies stratification (Whitney et al., 2012). This scenario occurs when there are strong tidal velocities and horizontal stratification gradients. Often freshwater input and mixing creates a stratification gradient whereby stratification increases toward the head of the estuary. Advection of the stratified water column, during ebb tides, increases stratification. During flood tides, advection decreases stratification, bringing less stratified water into the estuary.
We illustrate this schematically in Fig. 1a. An example of this would be Conway estuary in north Wales (Turrell et al., 1996).
Straining becomes important in estuaries with strong velocity shear and horizontal density gradients. For an estuary with increasing density toward the mouth, differential advection, or straining, created by the sheared ebb current brings freshwater over saltier water, increasing stratification. During flood tides, the differential advection works in the opposite sense, bringing saltier water over fresher water, reducing stratification (Whitney et al., 2012). We illustrate this schematically in Fig. 1b. Along with turbulent mixing, tidal straining is the dominant mechanism in estuaries such as the Chesapeake Bay, the York River, and Liverpool Bay (Li and Li, 2011;Scully and Friedrichs, 2007;Simpson et al., 1990). As most estuaries have spatial variations in stratification, density, depth-averaged velocity and velocity shear both advection and straining participate in changing stratification. We find that in NB straining and advection both play a major role in controlling stratification over a tidal cycle, which we explore in section 6.3.
In areas with strong velocity shear, mixing becomes important. Mixing and diffusion tend to reduce gradients. Vertical mixing, often driven by the bottom friction interacting with a flow field, will homogenize the water column vertically, reducing stratification. Vertical mixing is thought to often dominate mixing in the open ocean and estuaries, especially on shorter time-scales (e.g. Rueda and Schladow, 2009).
Typical vertical diffusivity values in the literature range from 10 -6 -10 -1 m 2 s -1 (e.g. Bowden, 1967;Geyer and Signell, 1992;Li et al., 2005). In our numerical model the ). In our model, with a time step of 20 s, vertical diffusion is more important at these shorter time-scales.

Data
In this study, we set out to compare stratification with environmental factors in order to assess the major driving mechanisms. High frequency data from moored stations in NB allows for comparison with environmental variables at various timescales, and is ideal for our analysis. In particular, we focus on how tides are responsible for the periodic changes in stratification. We do this using our idealized numerical experiments and find that our models are consistent with observations.
Below are descriptions of the data used, including buoys, environmental observations and our numerical model results.

Buoy Data
The Narragansett Bay Fixed Site Monitoring Network (NBFSMN) is a collection of sensors that measure temperature, salinity, dissolved oxygen and 140 chlorophyll in NB Stoffel and Kiernan, 2009). These buoys have a sensor located 1 m below the surface and another located at 0.5 m above the sediment-water interface (Stoffel and Kiernan, 2009). The buoys are primarily positioned in the northern half of the bay, with most of the stations located in the center or western side of channels and embayments (Fig. 2a).
Data is primarily collected in the warmer months from May through October, with sampling every 15 minutes. The exceptions are stations GB and TW, which have predominantly continuous records during all times of the year (Fig. 3). We choose to perform our numerical experiments during summer months, described below, when we can compare to the majority of the NBFSMN station observations. The NBFSMN data began in 1999, however, the period from 1999 to 2003 is not particularly useful, as only two stations were recording data. Half of the stations were online by 2003, and by 2008 all 10 stations were recording observations. Timelines of data availability, from 1999 to 2014, are shown in Fig. 3.

Environmental Data
We investigate four environmental factors that contribute to stratification: surface heat flux, river runoff, tides and wind. Surface heat flux is a combination of radiative fluxes, latent heat and sensible heat fluxes. Estimates, for the region around Narragansett Bay, are obtained from NOAA National Center for Environmental Prediction North American Regional Reanalysis (NARR). The third environmental factor, tides, are measured at tide gauge station 8452660 at Newport, RI (Fig. 2a). We use data from the time period 1999 to 2014.
Tidal heights, referenced to Mean Lower Low Water, are sampled at 6-minute intervals.
Finally, wind velocities, measured 10.6 m above sea level at hourly intervals, are obtained from the National Oceanographic and Atmospheric Administration (NOAA) tide gauge station 8452660 located in Newport, RI. Data is available beginning in October of 1999, providing over 14 years of comparison with stratification data.

ROMS Model
In addition to water column observations, we use the Regional Ocean Modeling System (ROMS) to characterize spatial variations and quantify changes in stratification to compare with our observations. ROMS is widely used in the coastal modeling community for simulating estuarine flows and is particularly well suited for modeling estuarine exchange, circulation and dynamics (e.g. Kremer et al., 2010;Lerczak and Geyer, 2004;Li and Li, 2012). We simulate estuarine processes in NB with ROMS due to the spatial limitations of the NBFSN buoys. We use a numerical model to constrain the horizontal variability within NB.
ROMS is a hydrostatic, primitive equation model that uses a structured curvilinear grid. The model domain includes NB and extends into Rhode Island Sound, capturing exchange between the continental shelf and the estuary (Fig. 2b). We perform two numerical experiments for this study. The first is a spin-up experiment to obtain reasonable hydrographic properties, and the second is to test how tides change stratification. In the first experiment, for initiation and validation, the model is forced by tidal constituents, horizontal boundary conditions, river inflow and meteorological surface conditions. Nine major tidal constituents are obtained from the ADvanced CIRCulation model (ADCIRC) (Mukai et al., 2002). Boundary conditions are forced with the results of regional model runs of the northeast U.S. shelf which were made available by the University of Massachusetts, Dartmouth and Woods Hole Oceanographic Institution (Finite Volume Coastal Ocean Model [Internet]., 2013).
The time series of boundary values are filtered with a 33 hour 4 th order Butterworth low-pass filter to remove tidal signals. The model is run from January 1, 2010 until July 31, 2010. The temporal evolution of the depth-averaged density, at NBFSMN locations, is displayed in Fig. 4. The depth-averaged density, in NB, decreases in the spring and summer.
Model surface forcing, for the spin-up experiment, includes the net surface radiation, surface air temperature, relative humidity, air pressure and wind speed.
Estimates of meteorological forcing are obtained from the Weather Research and We compare output from the spin-up numerical experiment with realistic forcing from January 1, 2010 to July 31, 2010 with NBFMSN buoy data. The purpose of this model run was to initialize summer stratification in the bay, not to represent the buoy stations exactly. We estimate how well our model predicts 2010 buoy data using the Willmott Skill (Willmott, 1982): where Tmod is model prediction of a variable, such as temperature or salinity, and Tobs are the NBFSMN observations. We calculate the Wilmott Skill for surface and bottom measurements of temperature and salinity as well as the vertical density difference measured at the buoy locations.
Over the time period modeled, we find temperature agrees the best with observations, having the highest Willmott Skill. Both surface and bottom temperature have a station average skill over 0.90. Skills for salinity and vertical density differences are on average greater than 0.55. Codiga (2012) found that the vertical density difference in NB was driven by salinity, therefore, the salinity and the vertical density difference skill should be similar.
One reason the salinity agreements might be low could be our choice of river forcing. We only used 4 rivers, and although, they are the largest, the rivers we chose do not account for the total freshwater input into the bay. Another potential discrepancy might come from the complex bathymetry of NB. Our model only resolves bathymetry variations down to 40 m and small-scale structures (~1 m) could have an influence on fluid flow especially at stations BR and CP, which are located adjacent to channels.  We are most interested in the effects that tides have on stratification, as rivers, wind and solar forcing have been well studied in NB and other estuaries (e.g. Codiga, 2012;Scully et al., 2005;Rippeth et al., 2001;. Therefore, the second experiment has no sources and sinks of stratification, except for riverine input.
Specifically, the ROMS model is run with no boundary forcing, no surface fluxes, no wind stress and only the M2 tidal constituents. We call this our idealized experiment.
The idealized experiment is started during a period when known stratification is present in the bay; i.e. during the early summer, at the end of July. The two major tributaries include the Blackstone and Taunton rivers. The constant freshwater flux is specified as 22 m 3 s -1 for the Blackstone and 14 m 3 s -1 for the Taunton. Our model provides a valid first order estimation of stratification as these two rivers account for about 55% of the freshwater flux into NB (Ries III, 1990). In addition, all but the M2 tidal constituents are removed to limit the change in stratification due to spring-neap variations in tidal velocities. We let the model come into periodic steady-state. This is determined by running the simulation until velocity, density and stratification are not changing at periods longer than a tidal cycle. In this way, we can isolate the contributions from the tides.
We confirm that our idealized experiment has a reasonable temperature and salinity distribution by comparing buoy locations in the model with real observations. After our idealized experiment comes into steady state, we plot near surface and near bottom estimates of temperature and salinity, averaged over a tidal cycle (Fig. 5 & 6).
We find that our model conditions are closest to June and July averages of NBFSMN temperature and salinity measurements (Fig. 5 & 6). The temperature of our steadystate model is cooler than the July averages as a result of no heat exchange with the surface as well as the interaction with a colder Rhode Island Sound.
Observational data was used from all available years at the buoy stations (Fig. 5 & 6). Observational data tends to be fresher than our model, but our numerical results fall within the standard deviation of the buoy observations for June and July. We expect the model to be saltier than the observations, as our idealized experiment only introduces freshwater at the Blackstone and Taunton, a subset of the total freshwater influx coming into the Bay. At stations in Greenwich Bay, GB and SR, vertical density differences do not agree well with observations. This is likely because our model does not include a freshwater source in or near Greenwich Bay.

Analysis
A goal of this study is to better quantify the variability of stratification both temporally and spatially. We do this through a combination of data exploration and comparison with known forcing variables. Our analysis of observations helps identify what scales and forcing conditions are important to study in our numerical model.

Analysis of Observations
We analyze the time series of both stratification observations and numerical modeling results. Stratification in the observational data specifically refers to the density (sigma-t) difference between the top and bottom NBFSMN sensors. The NBFSMN sensors record temperature and salinity, which is used to calculate sigma-t, with the MATLAB SEAWATER Library. Similar to Codiga (2012), we find that using the Brunt-V̈̈̈ frequency, or buoyancy frequency, would be misleading with this data set. This is because the stations are primarily fixed above and below the pycnocline. The Brunt-V̈̈̈ frequency could change solely based on the depth at the various stations, making inter-station comparison challenging.
Through temporal analysis of the observed stratification, we characterize the frequency content of the dataset. Since the power spectral densities are calculated with continuous signals, we pad times with missing data with station averages. This may bias some of the high frequency signals but does not affect our analysis. We use the Welch method to estimate the power of a signal at different frequencies. Specifically, the power spectral density is found from using Fast Fourier transforms of the autocorrelated stratification time series . We use a Hamming window of 2048 points and each time step is 15 minutes. The longer length of the time series window provides us with higher frequency resolution but increases the variance in the spectral density estimate. Overlap of the time series is set at 50%, to reduce the variance .
Calculating the magnitude squared coherence of stratification with environmental observations helps establish which environmental factors are important.
Magnitude squared coherence determines the frequencies at which two signals are most related. A value of 0 indicates no relation of the two signals, while a value of 1 indicates that the two signals are perfectly and linearly correlated. We use the mscohere function in MATLAB to calculate the coherence, which utilizes Welch's averaged modified periodogram method. The method utilizes the cross-spectral density of two-time series as well as the auto-spectral density of each time series to calculate coherence at each frequency of interest. Welch's average method uses multiple segments of equal length to improve results and reduce variance. The crossspectral density is a complex number and therefore, can be used to calculate the phase lag between two signals. For each station, we compare and calculate the coherence of the vertical density difference with four data sets; surface temperature, river discharge, tides and northward wind.
We determine the significance of each magnitude squared coherence with a method developed by Bendat and Piersol (2010): where Cxy is the coherence at a given frequency and df is the degrees of freedom. We calculate the degrees of freedom as = .Ç >ÈÈG , where L is the total length of the time series and nfft is the length of the Hamming window. If the magnitude squared coherence is greater than , then the coherence is considered significant at that frequency.

Analysis of ROMS
We use the numerical output of our two ROMS experiments to develop a stratification scenario and to explore possible causes of stratification change at tidal time-scales. Our spin-up experiment is used for direct comparison with observations for the year 2010. We calculate the stratification as the density difference between near-surface and near-bottom levels at NBFSMN station locations in the model output.
Additionally, we explore frequency analysis on the output from our spin-up numerical model to assess the similarity between the observations and model. It is important to note that the numerical spin-up experiment does not produce the exact results observed during the modeled year 2010. What is more important is that the model reproduces the same scale of variation that we observe at NBFSMN stations, which is why spectral analysis is useful.
We use our idealized experiment to predict the change in stratification, to within the capability of our experiment, at a finer spatial scale than our observations. We choose to use a definition related to the vertical gradient of density to define stratification. Spatial variability of stratification is explored with the numerical model, which carries out precise calculation of density and velocity. The study of the stratification balance allows for the evaluation of mechanical changes by analyzing advection, straining and differential diffusion (Chen and Sanford, 2009;Li and Li, 2011 (1) where is density, is the velocity vector and M and 9 are the eddy diffusivity for the horizontal and vertical, respectively. The rate of change in density is a function of advection, the first three terms on the right-hand side, and diffusion, the second three terms on the right-hand side of eq. (1). We assume that there are no additional sources from heat or freshwater. This assumption is valid for our numerical idealized experiment in which surface fluxes are turned off and rivers are held constant. The exception is at the river nodes, where freshwater comes into the model. We exclude these nodes from our analysis. To find the stratification balance we multiply eq. (1) by w−« where . is defined as −« • -« - † , _ is a reference density and g is the gravitational acceleration. In addition to an advection and differential diffusion term, the stratification balance now includes a straining term. Stratification will change due to straining when a vertical shear is applied to a horizontal density gradient.
ROMS uses the structured Arakawa C-grid and therefore, careful attention must be used when calculating the N 2 balance with a finite difference scheme.
Calculations need to consider the volume of the water parcel and the volume of material being advected, strained and diffused into or out of the parcel of interest. Our analysis is based on methods by Li and Li (2011), who describe in detail how to The next two terms in eq. 4 are the depth-averaged advection terms calculated in both horizontal directions. The following two terms are the depth-averaged straining terms. The last three terms on the right-hand side of eq. 4 are the differential diffusion terms. Note that these terms are the vertical integration of a vertical gradient.
Therefore, the last three terms can be evaluated at the surface and bottom sigma layers. An error estimate is found by subtracting all the calculated terms on the righthand side of eq. 4 from the estimate of -Ã B YYYY -G , the left-hand side of eq. 4. Instantaneous values are used to calculate each term on the right-hand side of eq. 4. However, calculated . using output every 6 minutes. Therefore, there is a slight error associated with our calculations as -Ã B YYYY -G is calculated between output intervals.

Temporal and Spatial Patterns
Periodic temporal variations in water temperature, salinity and stratification are observed at each NBFSMN station. For discussion, we look at one representative station, CP, and one representative year, 2010 (Fig. 7). The most obvious signal is the seasonal cycle observed in the temperature time series, with the warmest temperatures occurring in late summer and early fall (Fig. 7a). Peak temperatures occur at different times, with the surface maximum in July (26 o C at CP) preceding the bottom temperature maximum in mid-August (24 o C at CP). Vertical temperature differences that help drive stratification are present in the spring and are negligible by early fall (Fig. 7b). The maximum vertical water temperature difference was recorded at 10˚ C in July. The water column at CP is nearly uniform in temperature during the fall and winter months. This is true for other stations in NB as well. BR, TW and GB recorded data during winter months and confirm that vertical temperature differences were negligible during the winter.
Salinity is not uniform throughout the year. At station CP, salinity has an annual periodicity with the highest salinity values recorded in September (Fig. 7a).
The peak is most likely due to reduced river flow and increased evaporation rates.
Because of the importance of salinity on stratification, we look at the vertical salinity differences (Fig 7b). Vertical salinity differences, at station CP in 2010, show a seasonal periodicity with a minimum occurring in September. The largest salinity differences are in the spring, with magnitudes of approximately 5 PSU.
Density exhibits less variability than temperature and salinity, but still shows an annual periodicity. Bottom density at CP is much less variable when compared with the surface density (Fig. 7c). The station CP data indicate a seasonal trend in stratification that is maximized in late spring and mid fall and minimized during summer months (Fig. 7d & 8). When averaged over individual months, stratification at other stations shows similar patterns, with maximized vertical density differences during the spring months of May and June and minimum values in late summer (Fig.   8).
We calculate the percent contribution of salinity to vertical density differences by dividing the vertical density differences computed with the average temperature at each station with the actual vertical density differences as described by Codiga (2012): where ∆ is the vertical density difference, S is the in situ salinity, T is the in situ temperature and 〈 〉 is the temporally averaged temperature measured at a given station. Results are listed in Table 2. Contributions of salinity dominate the vertical density differences for warmer months, with percentages ranging from 70 to 88%.
Stations in the north such as BR and GB have the highest contribution to stratification from salinity and stations in the south like QP and TW exhibit the lowest. This illustrates the influence of freshwater input from the northern tributaries. General spatial trends are observed when plotting the average vertical density difference against latitude (Fig. 9). Temporally-averaged vertical density differences at each station are calculated from data collected May through October from all available years. Average stratification increases from under 1 kg m -3 at the most southern station, TW, to almost 4 kg m -3 at the northernmost station, BR. Stratification increases with latitude. We compared sensor depth difference with average stratification and latitude to confirm that there is not a bias due to individual station configuration. Station sensor depth differences range from 1 m, at GB and SR, up to 9 m, at NP. No trend was found when we compared sensor depth differences with either stratification or latitude. This indicates that the correlation of vertical density differences with latitude is not a result of instrument setup. We conclude that the latitudinal increase in stratification is likely due to the freshwater influx from northern river input.  (Fig. 9). These stations could be less stratified compared to PP due to greater differential diffusion and/or the reduced transport of fresh or salty water.
Station PP is likely more stratified due to the connection with the mouth of NB and the influx of freshwater from the Providence River Estuary. Comparison of stations at similar latitudes indicates, that of the three passages connected to the mouth of the bay, the East Passage is the most stratified. This channel is flanked on either side by less stratified areas of the West Passage and Mt. Hope Bay.

Power Spectral Analysis
To look at higher frequency temporal trends, we perform power spectral analysis on vertical density differences from NBFSMN stations. We plot results in a variance preserving power spectral density to highlight the higher frequencies ( Fig.   10). Power spectral density of the stratification reveals several frequencies that have strong influence over the observed signal. Strongest peaks present, from lowest frequency to highest, include a diurnal (1 cpd Similar to the observations, the four most prominent frequencies in the spin-up ROMS experiment are diurnal, the M2, the M4 and M6. Spectra computed from the ROMS output sampled at the NBFSMN buoy locations are shown in Fig. 11. The model is resolving mechanisms that change stratification at tidal frequencies. The spin-up experiment has less variance compared to observations (Fig. 12). This spectral analysis does not conserve variance. Amplitudes of at the M2 frequency are on the O (10 -2 -10 -1 ) kg 2 m -6 cpd -1 , averaging 0.05 kg 2 m -6 cpd -1 for observations. The average of the M2 peak amplitude is also 0.05 kg 2 m -6 cpd -1 for the spin-up experiment output. The amplitude of the higher peak frequencies (M4 & M6) are lower than observations.

Coherence
To explain the variations observed in the vertical density differences, we compare observed stratification measurements to four key environmental parameters.
Two parameters are sources of stratification (i.e. net heat flux and river discharge) and two are forcing mechanisms (i.e. tides and wind). We assess the likelihood that these mechanisms affect stratification by evaluating their coherence with stratification.
Coherence between vertical density differences and net heat flux, Blackstone River input and northward wind measured at Newport, RI is not statistically significant for most stations. Coherence is determined over periods of 30 minutes to 11 days to encompass the frequency range of the observed spectral peaks. The significance level is calculated by taking into account the degrees of freedom in our analysis (Bendat and Piersol, 2010). The river discharge data is collected daily and linearly subsampled to 15 minutes. While these parameters do not influence vertical density difference over the period band between 30 minutes and 11 days, they likely influence NB stratification on longer time scales. For example, Codiga (2012) has shown the importance of for riverine input at monthly timescales.
The strongest coherence, between vertical density differences and any variable, is obtained when comparing stratification with tidal heights measured at Newport, RI.
We highlight the frequencies that are considered significant for each station (Fig. 13).
The strongest coherences are measured at the diurnal, M2, M4 and M6; i.e. tidal frequencies. Seven of the ten stations all have their strongest coherence at the M2 frequency and all, but station MH, had significant coherence at this tidal frequency ( Fig. 13).
A combination of varied magnitudes across straining, advection and differential diffusion mechanisms can create a stratification maximum at any time during the tidal cycle (Whitney et al., 2012). We explore the timing of stratification by calculating the phase lag between observed tidal heights and vertical density differences at various frequencies (Fig. 14).  (Fig. 15a). The phase lag increases at the northern stations. At these stations, stratification lags high tide by a half a period (6 hours), around low tide (Fig. 15a).
The tidal wave that proceeds up NB could hypothetically explain the phase lag between southern station and northern stations. For this to be true the tidal height phase lag must be on the order of hours. However, there is only a 15-minute lag in tidal height from the tide gauge stations of Providence and Newport. Therefore, the tidal wave travels too fast to account for the lag of several hours in stratification observed across stations.
A possible explanation for the variations in stratification timing involves straining. Straining-induced stratification suggests that with an increasing density towards the mouth of an estuary, differential advection should maximize stratification during ebb. In this scenario, stratification should lag by half a period behind tidal heights, as stratification should be out of phase with tides (Scully and Friedrichs, 2007). In addition to straining, pure advection could produce the same result, if stratification increases from the mouth of an estuary to the head and mean velocities are strong (Whitney et al., 2012). Therefore, we postulate that either straining or pure advection could be responsible for the stratification maximum near low tide found at the northern stations of BR, CP, GB, PP and NP. We further explore the combination of mechanisms that causes the phase lags for all the NBFSMN stations in following sections.
Power spectral analysis shows the dominant frequencies present in our idealized experiment. We identify a peak at the M2 frequency. After calculating the coherence of the vertical density differences with tidal height at Newport, RI, we plot the phase lag for the M2 frequencies in Fig. 15b. This characterizes the lag of stratification behind the tidal height. We plot the phase lag of the same stations as the observational analysis (Fig. 15a) at the M2 frequencies.
Station phase lag between vertical density differences and tidal height at Newport, RI are similar to the observations in Fig. 15a Since the coherence of stratification with tidal heights is present at many stations and at tidal frequencies, we conclude that the tidal flow is likely causing the change in stratification observed at the M2 frequency. The variability in phase lag exhibited across the stations suggests that a combination of tidally-driven advection, straining, and differential diffusion play a role in changing stratification. We explore these mechanisms further in following sections.

Tides
Dynamic changes in stratification are created by water movement and turbulence. Observations of NBFSMN stratification indicate changes on the order of 0.20 kg m -3 per tidal cycle at all stations (Fig. 12). Coherence between stratification and tidal heights at Newport is the strongest amongst the four environmental variables (tides, heat flux, river input and northward wind) compared with stratification. The coherence, along with previous studies on tidal straining (e.g. Scully and Friedrichs, 2007), suggest tidal currents have a large impact on stratification changes in NB.
To further explore the mechanisms that control stratification at tidal frequencies, we use our idealized numerical model to investigate the stratification balance. We analyze the terms in eq. 4 computed for all nodes located in NB.
The largest 〈N . YYYY 〉, on the order 10 -2 s -2 , is measured near the Blackstone and Taunton Rivers. This is where the freshwater is introduced into the model. Stratification decreases by an order of magnitude towards the mouth of the bay. Shallow embayments, like Greenwich Bay and the northern part of Mt. Hope Bay, also display low stratification.
We characterize the change in depth-averaged stratification (∆N . YYYY ) over a tidal cycle by taking the difference between the maximum modeled N . YYYY and 〈N . YYYY 〉, similar to Whitney et al. (2012). Results are displayed in Fig. 16b. ∆N . YYYY changes over a tidal cycle from 10 -4 to 10 -2 s -2 , with the largest changes occurring coincident with elevated stratification. We normalize the ∆N . YYYY by the 〈N . YYYY 〉 (Fig. 16c). The average percent change is between 10 and 30 %. Stronger variations, up to 90%, are measured in shallow areas such as the periphery of Ohio Ledge or the northern area of Mt. Hope Bay. The strong changes (red in Fig. 16c) around Greenwich Bay are likely not modeled properly due to the lack of local riverine sources in this area.
An N . YYYY maximum during the tidal cycle is measured at every position in NB, illustrating the temporal variability throughout the bay (Fig. 16d). To quantify when the maximum N . YYYY occurs, we take the cross correlation between N . YYYY and the depthaveraged up-bay velocity. The M2 is the largest contribution as seen with the variance preserving power spectral density for most stations (Fig. 11). Therefore, we evaluate the phase lag in the time domain. The phase lag at which the correlation is maximized provides the lag between the two-time series. In our experiment, we divide the phase lag into four categories each 3 hours long, i.e. a fourth of the tidal period. Most of the bay has peak stratification during slack low tide but there is a wide range of phases (Fig. 16d).
The large-scale evolution over a tidal cycle can be described through the distribution of salinity, temperature and stratification throughout the four tidal cycle epochs (Fig. 17, 18, and 19): a) From high tide to ebb, velocities are primarily down bay, bringing freshwater towards the mouth of the bay. b) During slack low tide, the velocities weaken, and for many areas, this is coincident with a salinity minimum (Fig. 17c). For example, just south of Greenwich Bay along the western side of the west passage the fresh-water tongue reaches it southernmost extent during slack low tide (Fig. 17c).
c) During maximum flood, the tide enters the bay, up-bay velocities are maximized, and the salt front starts to move up-bay (Fig. 17d). Throughout this phase we see stratification start to retreat or dissipate, especially within channels (Fig. 18d).
Particularly looking at Ohio Ledge, we see the river plume that was apparent during maximum ebb and low tide start to dissipate. d) Once NB reaches high tide, most surface salinities are maximized as a result of the intrusion of Rhode Island Sound water and stratification is generally at a minimum ( Fig. 17). Only slight temperature differences were observed throughout the tidal cycle ( Fig. 19). Note the high temperatures towards the north are a result of warmer input in rivers, while the high temperatures in Greenwich Bay are a result of the model configuration. Greenwich Bay is an area that has a longer residence time and is shallow, allowing heat to be absorbed in this area during the spring and early summer modeled spin-up period.

Stratification Driving Mechanisms
We next explore what causes the temporal variations in N . YYYY by estimating contributions from advection, straining and differential diffusion. We start with an illustrative location, station Quonset Point (QP), and evaluate the terms from eq. 4.
We plot the station's tidal height relative to mean sea level and the up-bay ( ) depthaveraged velocity in Fig. 20a. The up-bay velocity precedes high tide by a quarter of a period (3 hours). This is characteristic of a standing wave. Our model suggests that most of NB behaves like a standing wave. Velocity and tidal height measurements from Quonset Point and Fall River tidal and current gauges also suggest a standing wave in Narragansett Bay from observations (Tidal Current Tables 2016-Atlantic Coast of North America, 2016). Furthermore, our numerical model indicates that stratification at QP is in phase with the tidal velocity, maximized during flood and minimized during ebb (Fig. 20b).
The mechanisms responsible for -Ã B YYYY -G at QP are plotted in Fig. 20c. The major terms, advection, straining and differential diffusion, govern -Ã B YYYY -G and vary over a tidal cycle. At station QP, advection has the largest amplitude and dictates when . YYYY is 164 maximized. When advection is maximized around 0.38 and 0.9 days, . YYYY increases (Fig. 20). When advection is negative around 0.1 and 0.5 days, . YYYY decreases (Fig.   20). It should be noted that straining and differential diffusion play a role in the stratification balance. Straining both increases and decreases . YYYY . It also tends to act against advection. Differential diffusion has a smaller variation but always reduces . YYYY .
We explore advection and straining terms further at station QP, by plotting the horizontal components of the two terms in Fig This implies advection determines when the change in stratification occurs. Also similar to QP, most of the bay has a larger variance in the direction for both advection and straining. We explore the size of the three mechanisms in the rest of the bay by determining the variance of terms in the stratification balance.

Variance
Unlike Simpson et al. (1990), who found that tidal straining in many estuaries was the primary cause of increased (decreased) stratification during ebbing (flooding) tides, we find that a combination of mechanisms is responsible for observed temporal trends in the bay. This is suggested from the modeled . YYYY time series described in the previous section, and further supported here by our variance analysis.
The major contribution to -Ã B YYYY -G is from advection, straining and diffusion.
Convergence/divergence plays a minor role in changing stratification and our error 165 estimates are also small. We proceed with our variance analysis by normalizing the major three terms. Normalized variance is calculated by: where ? is the variance of the advection term, . is the variance of the straining term and x is the variance of the differential diffusion for any location in NB. By comparing the relative strength of the variance of one component to the total variance, we can determine the relative contribution of each component to the stratification balance over a tidal cycle.
We use ternary diagrams to plot the normalized variance of the three variables.
The total variance sums to a constant, in this case S1+S2+S3=1. Each component, advection, straining and differential diffusion, is plotted along one side of the triangle.
This allows one to represent three variables in two dimensions. Each corner represents 100% of the variance contributed to -Ã B YYYY -G . by that component. For example, the top corner represents samples that have only variations due to the differential diffusion component; where S3 = 1 (100%), S1=0 is 0% variation in advection and S2=0 is 0% variation in straining. The percentage of any component decreases linearly away from that corner. Therefore, the bottom horizontal edge of the triangle represents 0% variance contributed by differential diffusion.
After calculating the normalized variance for all positions in NB, we plot them on a ternary diagram in Fig. 21a. Warmer colors indicate a higher density of samples.
Each sample is a location in NB ROMS model. The ternary diagram illustrates that the largest contribution of the variance of the stratification change is from advection and 166 straining. Although, on average advection contributes a higher percentage to the variance, the straining contribution is comparable in magnitude.
Variance analysis shows that advection, straining and differential diffusion contribute differently to -Ã B YYYY -G at a given location. The ternary diagram (Fig. 21a) illustrates that advection and straining both contribute more than differential diffusion.
Our analysis shows that in only about 5% of NB, differential diffusion dominates the balance. The other 95% has either advection (65%) or straining (30%) as the largest contributor to -Ã B YYYY -G variance.

Phase
The To determine what causes the . YYYY maximum we need to determine the phase of advection and straining. We take the depth-averaged advection term and crosscorrelate it with the depth-average velocity in the direction. Correlations indicate the 167 time during the tidal cycle when advection is maximized. There are two peaks in the 2-D histogram of these results (Fig. 22e). One peak occurs at maximum flood, in phase with the tidal velocities and the other peak is right after max ebb. This indicates advection is out of phase with the tidal cycle (Fig. 22e).
We take the depth-averaged straining term and cross correlate it with the depth-averaged velocity in the direction. Correlations indicate most of the straining is either maximized at maximum flood or just after maximum ebb (Fig. 22c). The maximum flood peak is larger than the max ebb peak indicating more of the bay has a straining term that is directly in phase with the tidal velocity cycle in the direction.
The two peaks in advection and straining provide valuable information about when the term separately is maximized but do not explicitly explain the distribution of stratification phase lags observed over the bay (Fig. 16d).
Each grid cell in our model represents a sample where we have estimated the timing of maximum . YYYY , advection and straining. The comparison of when the . YYYY is maximized against when advection is maximized provides more detail about the dynamics occurring in NB (Fig. 22b). A 2-D histogram of NB timing of maximum . YYYY versus the timing of maximum advection, colored by the density of the number of samples, reveals there are primarily two regimes. One regime occurs when advection is maximized at maximum ebb, resulting in . YYYY being maximized during slack low tide (Fig. 22b). The other regime occurs when advection is in phase with the tide and maximized during maximum flood. . YYYY is maximized just after high tide for this regime (Fig. 22b). Potentially a third regime is noticed where advection is still in phase with the tidal cycle, but . YYYY is maximized during slack low tide.

Discussion
Estuaries, such as Liverpool Bay, are dominated by straining (Simpson et al., 1990), while other estuaries have found to have a balance between straining and advection (Whitney et al., 2012). Our study illustrates that there can be multiple regimes within one estuary. In the following sections, we focus on two major regimes within NB.

Two Advection Dominated Regimes
The two dominant peaks in Fig. 22b suggest there are two coherent regimes controlling -Ã B YYYY -G in NB related to advection. The high density of samples indicates that many parts of the bay behave similarly. We explore how advection and straining changes stratification over a tidal cycle by defining two regimes. The first regime is defined as any position in NB with a maximum advection around max ebb and maximum . YYYY occurring around low tide. We call this regime 1 (Fig. 23a). Regime 2 is defined as any position that experiences a maximum advection around max flood and maximum . YYYY around high tide (Fig. 23a). Both regimes use a cut-off sample density of 0.25 counts/deg 2 to define the region. These contours are outlined in the blue and red for regime 1 and 2, respectively (Fig. 23a). It appears that much of the bay's channels are included in regime 1, while embayments like Greenwich Bay or inter-channel areas tend to fall into regime 2 (Fig. 23b).
In NB straining and advection tend to be anti-correlated. This is illustrated in Fig. 22d, where straining and advection are largely out of phase. Therefore, straining works against advection in both regimes. As a result of straining being consistently anti-correlated with advection, we can reclassify major regimes based on one term.

Slack Low Tide Maximum Stratification, Regime 1
Ohio Ledge is an example of regime 1 dynamics. Tidal evolution at Ohio Ledge is shown in Fig. 24. Density profiles are shown at various locations and indicate the steepest gradients during maximum ebb, when the river plume emerges into Ohio Ledge (Fig. 24 c & d). Spatially averaged terms (eq. 4) over Ohio Ledge, indicate that . YYYY is maximized in-between maximum ebb and maximum low tide ( Fig.   25 a & b). The contribution from diffusion is always negative at around -5x10 -8 s -3 and convergence and divergence has almost no net effect in eq. 4 (Fig. 25c). Straining balances diffusion in this area and has a positive value with average 5x10 -8 s -3 .
Straining is positive in this area because of the induced estuarine two-layer flow. The vertical gradient of velocity and the along-estuary horizontal gradient in density are negative, resulting in a net positive straining term.
-Ã B YYYY -G is thus a result of the advection variation, as straining and diffusion balance one-another. This can be seen in the similarity in phase, magnitude and shape of the advection term (Fig. 25 c) and Fig. 25b).
A cross-section through Ohio Ledge (Fig. 26) illustrates that during maximum flood stratification increases up-bay (right of figure) and density increases down-bay ( Fig. 26). Maximum flood results in minimum stratification as up-bay velocities move unstratified waters toward the head of the bay (Fig. 26 a & c). During maximum ebb, stratification increases due to the advection of the freshwater plume into Ohio Ledge ( Fig. 26 b & d). Although, straining is positive throughout the tidal cycle, it tends to be balanced in this area by differential diffusion, not influencing -Ã B YYYY -G greatly.

Slack High Tide Stratification Maximum, Regime 2
The region between Jamestown and Prudence Island, referred to here as interchannel, represents an area that is characterized as regime 2. In general, . YYYY is maximized at slack high tide as more stratified water advects from the East Passage to the West Passage (Fig. 27). Vertical density differences are larger during flood and slack high tide than during parts of the tidal cycle (Fig. 27 b Fig. 28 b) is smaller than the advection term ( Fig. 28 c), the phases are similar.
Straining, in this location, reduces the effect of advection as it is out of phase with advection ( Fig. 28c). Diffusion always reduces . YYYY (Fig. 27c).
A schematic representation of regime 2 illustrates a bottom intrusion in Fig. 29. An imaginary cross-section is drawn northwest (right of figure) to southeast (left of figure), in between the two islands. In this area, the pycnocline intersects the bottom topography. Regime 2 is controlled by the pycnocline moving up-bay during flood which increases . YYYY and similarly, as the pycnocline moves down-bay during ebb . YYYY decreases.

Advection & Straining
Together, regime 1 and 2 account for 35 % of the bay. In the other areas, the sign and phase of advection is controlled by the direction of the horizontal stratification gradient because tidal velocities are directionally invariant within the model. However, the variance amplitude of the straining term becomes more important in these other areas. Peak . YYYY can occur at any time during the tidal cycle 171 depending on the magnitude of both straining and advection. We identified the most coherent regimes present in NB, through the phase of . YYYY , advection and straining, and note the rest of NB depends on the balance of all terms in eq. 4.

Comparison to Other Estuaries
Estuaries can be driven by pure advection or pure tidal straining. For example, the Chesapeake Bay, is thought to be dominated by straining (Li and Li, 2011).
Advection likely plays more of a role in NB because the stratification gradients change faster over a shorter distance. NB is different from larger estuaries like the Chesapeake, due to its length and shape.
The length of NB, is only a few 10's km, and it has multiple embayments and relatively shallow areas. These factors contribute to juxtaposing regions of high stratification next to areas of low stratification, making advection as important as straining. This is a result of strong horizontal stratification gradients.
Numerical modeling of estuaries has started to explore spatial variations in mechanisms, however there is a lack in documentation through observations (Burchard and Hofmeister, 2008;Giddings et al., 2011). As more estuaries are studied with numerical techniques, variations in mechanisms have been found within one estuary (Whitney et al., 2012). One such case is the Rhine river outflow region (Rijnsburger et al., 2016). Rijnsburger et al. (2016) found that advection and straining are opposing and lead to various timings of the maximum stratification, similar to our NB modeling results.
The analysis of the NBFSMN network observations and NB numerical model has found NB to experience multiple regimes. More detailed studies of other estuaries may reveal more spatially variant regimes as well. We suggest that it is important, even within the same estuary, to consider different stratification regimes as it effects turbulent mixing, physical, biological, geological and chemical processes.

Impact on Mixing
Mixing is dependent on stratification and shear (e.g. Whitney et al., 2012;Peltier and Caulfield, 2003). Shear increases turbulence and is commonly maximized during flood or ebb tides, when velocities are often greatest. Stratification can limit turbulence, as it affects the gradient Richardson number Ri (N 2 /S 2 ) and eddy viscosity (K) (Whitney et al., 2012). S 2 is defined as (|du/dz|) 2 . In addition to tides, wind can also create shear and induce additional mixing. Our numerical model predicts stratification varies over a tidal cycle throughout NB. We postulate that mixing efficiency from wind events lasting a couple hours would result in different amounts of mixing depending on the phase of the tidal cycle.
Ohio Ledge is a good example of an area that experiences large changes in stratification. This area would be susceptible to variable mixing that would depend on the time a wind event occurs in relation to the tidal cycle. During ebb tide we predict most of Ohio Ledge to have a maximum in stratification (Fig. 24b). However, during maximum flood, a small part of the periphery reaches its maximum stratification ( Fig.   24d & Fig. 16 d). A strong wind event applied to our model, lasting on the order of hours would be more effective at mixing the periphery during ebb tides than during flood tides. During such a wind event, we would expect stronger residual flows around the periphery due to baroclinic gradients, as the boundary between well mixed and stratification is enhanced.

Conclusion
Spatial and temporal stratification patterns in NB are primarily driven by seasonal, synoptic and semi-diurnal variations of environmental parameters including winds, run-off, solar heating and tides. The analysis of the multiple year buoy network in NB and hydrographic numerical modeling of NB have facilitated the characterization of stratification changes. Spatially, there is a latitudinal increase in stratification resulting from the freshwater input from the northern tributaries. This spatial trend is temporally changed due to several environmental factors. Over yearly time-scales, freshwater fluxes control the general stratification with more freshwater increasing stratification in the spring. These observations are in agreement with Codiga (2012), who found that riverine input was the most important source of stratification on monthly time-scales in NB.
On shorter time-scales, we find stratification varies at tidal frequencies.
Periodic stratification changes are observed at NBFSMN buoys at diurnal, M2, M4 and M6 frequencies. Power spectral densities of observations reveal stratification changes on the order of 0.5-1 kg m -3 at the M2 frequency. Through numerical modeling we found that advection and straining play a key role in controlling stratification changes.
Our model results indicate that in NB stratification is not spatially uniform. NB is a relatively shallow estuary with complex bathymetry. This leads to more mixing especially in shallow areas, and a non-uniform stratification gradient throughout NB.   Burchard and Hofmeister (2008).     : Temperature and salinity plots for (a) 1 meter below surface and (b) 0.5 meters above the seafloor at NBFSMN station locations. Buoy data is averaged from June (circles) and July (asterisk) of all available years. Standard deviation of measurements is plotted as grey error-bars. After our idealized ROMS experiment comes into steady state, station averages of temperature and salinity (triangles) are plotted.         : Time series of depth-average properties at station QP from ROMS steady state idealized experiment. (a) Tidal height (black) relative to mean sea level and northward depth-averaged velocity (grey). (b) Depth-averaged stratification (black) and rate of change (grey). (c) Components from eq. 4, including advection (blue), straining (red) and differential diffusion (yellow), convergence/divergence (purple) and error (green). (d) Advection components (eq. 4) broken up into the -direction (light blue), -direction (darkest blue). (e) Straining components (eq. 4) broken up into the -direction (light red), -direction (darkest red). and are primarily in the across-and along-bay directions.        Fig. 27. (a) Tidal height (black) relative to mean sea level and depth-averaged velocity (grey). (b) Depthaveraged stratification (black) and rate of change (grey). (c) Components from eq. 4 including advection (blue), straining (red), differential diffusion (yellow), convergence/divergence (purple) and error (green). (d) Advection components (eq. 4) broken up into the -direction (light blue) and -direction (darkest blue). (e) Straining components (eq. 4) broken up into the -direction (light red) and -direction (darkest red).