Dye, B.; Jose, F., and Allahdadi, M.N., 2020. Circulation dynamics and seasonal variability for the Charlotte Harbor Estuary, Southwest Florida coast. Journal of Coastal Research, 36(2), 276–288. Coconut Creek (Florida), ISSN 0749-0208.
A hydrodynamic model was developed and validated for the Charlotte Harbor estuarine system, located in SW Florida, to elucidate freshwater fluxes within the system's various inlets during diverse hydrologic conditions. Fresh water entering the system not only varies seasonally but also, because of regulatory fresh water, releases controlling water levels within an upstream lake. The unnatural freshwater releases have been found to negatively affect the system's ecology, in particular within the Caloosahatchee River portion of the system. Neither the flood nor ebb phase exhibits uniform dominance in flushing the system's four major passes. Boca Grande Pass and Big Carlos Pass were mostly ebb dominant, whereas San Carlos Bay was largely flood dominant; neither phase dominated at Captiva Pass. The similarities and/or contradictions of these results in comparison to former field and modeling results are mainly attributed to the differences between the freshwater sources and environmental forces corresponding to each study that forces a different mass-balance condition over the estuary-bay system and, thereby, at each individual inlet. A Lagrangian particle tracking study revealed particles released within the Peace River during different hydrological conditions were comparably transported regardless of freshwater inputs and predominate wind direction. In contrast, particles released within the Caloosahatchee River were flushed into the Gulf of Mexico within 10 days during a usually wet El Niño, dry (November–April) season period whereas during the summer wet (May–October) season released particles remained in the estuary for a longer period (13 days), ultimately resulting in their further transport into Pine Island Sound and Matlacha Pass. The results also demonstrate the effect of freshwater river inputs and wind on the travel time of the neutrally buoyant particles within the estuarine system. The hydrodynamic and coupled particle tracking model serve as the first step in a forthcoming larval transport modeling study.
INTRODUCTION
The Charlotte Harbor estuarine system is located in SW Florida (Figure 1) and is considered microtidal with a semidiurnal mixed tidal range of 0.6 m (Scarlatos, 1988). Additionally, the region is regarded as being subtropical because of distinct wet (May–October) and dry (November–April) seasons (Sun, Wan, and Qui, 2016) with average (2004–16) total seasonal precipitation of 780.9 and 320.7 (mm), respectively (Dye, 2018). The estuarine system has a mean depth of 2–3 m and covers a surface area of ∼700–800 km2 (Harris et al., 1983; Poulakis et al., 2004; Zheng and Weisberg, 2004). Marine water from the Gulf of Mexico enters the estuarine system through various inlets (e.g., Boca Grande Pass, Captiva Pass, Redfish Pass, and Big Carlos Pass) between the barrier islands and San Carlos Bay (Figure 1). Generally, the estuarine system is considered to be well mixed; however, stratification can develop during times of high freshwater inflow commonly associated with the wet season (Goodwin, 1996; Weisberg and Zheng, 2003). River inputs, along with precipitation and the associated watershed runoff, serve as the estuarine system's primary sources of fresh water; Myakka and Peace Rivers supply the northern portion, the Caloosahatchee River supplies the middle region, and the Estero and Imperial Rivers supply the southern region (Figure 1). The discharge rates of the Peace, Myakka, Estero, and Imperial Rivers are directly related to drainage basin rainfall, whereas, in contrast, the freshwater inflow into the Caloosahatchee River, which feeds the San Carlos Bay portion of the estuarine system, results from drainage basin rainfall in addition to controlled regulatory releases from Lake Okeechobee (Doering and Chamberlain, 1999; Steinman, Havens, and Hornung, 2002).
Figure 1
The Charlotte Harbor estuarine system study area showing locations of water level, currents and wind measurement stations (1–10, W), particle release points (R), and surrounding geographic area, including Lake Okeechobee and S-79 Franklin Lock & Dam.
A series of lock and dam structures located within the Caloosahatchee River regulate the amount of fresh water released into San Carlos Bay from Lake Okeechobee (Figure 1). The lock and dam structures, S-77 and S-78, control the water level of Lake Okeechobee located farther upstream, whereas the final lock and dam structure, S-79, acts as a salinity barrier separating the marine and fresh water (Chamberlain and Doering, 1998).
Freshwater flow at S-79 (i.e. fresh water released from the dam) is highly seasonal with generally greater flow during the wet season as opposed to the dry season (Qui and Wan, 2013). Seasonal freshwater flow patterns play an important role in the salinity gradients, which develop below the S-79 lock and dam within the Caloosahatchee River Estuary (CRE) and San Carlos Bay (Chamberlain and Doering, 1998). During the dry season, fresh water released from S-79 is usually minimal, resulting in salt water traveling farther into the CRE. In contrast, during the wet season, San Carlos Bay can experience persistently low salinities (<10 Practical Salinity Units) because of the high volume of fresh water discharged from S-79. Water management practices related to the wide fluctuations in freshwater releases have been found to adversely affect the diverse ecosystems (e.g., oyster reefs, seagrass beds) within the CRE and San Carlos Bay (Barnes, 2005; Barnes et al., 2007; Chamberlain and Doering, 1998; Tolley, Volety, and Savarese, 2005; Volety et al., 2009).
Several hydrodynamic models have been developed for the Charlotte Harbor estuarine system (Goodwin, 1996; Weisberg and Zheng, 2003; Zheng and Weisberg, 2004, 2012) and the CRE (Qui and Wan, 2013; Scarlatos, 1988; Sun, Wan, and Qui, 2016; Wan et al., 2013). Those hydrodynamic model studies primarily focused on circulation and flushing dynamics (Goodwin, 1996; Scarlatos, 1988; Weisberg and Zheng, 2003; Zheng and Weisberg, 2004, 2012), residence time (Wan et al., 2013), and salinity modeling (Sun, Wan, and Qui, 2016). The prior model studies served as important references and provided guidance for the development of the hydrodynamic model presented in this study. The focus of this study is to quantitatively evaluate the relative dominance of four tidal inlets in flushing the Charlotte Harbor estuarine system during different hydrologic conditions. Moreover, an effort has been made to evaluate the transport of neutrally buoyant particles—proxy for invertebrate larvae in a forthcoming study—within the estuary under the influence of tides and seasonally varying freshwater discharge conditions, which will ultimately help to understand the transport of invertebrate larvae within the estuarine environment. Several new features about the study area occur, which the present study addresses. This study demonstrates the effect of different sources/forces on the circulation of the estuary-bay system by including three different simulation scenarios and including freshwater seasonality, transient wind field, and different tidal cycles. It was clearly shown that the flood/ebb dominance of the fluxes for different inlets highly depends on two major participants, i.e. river discharge and wind. Showing the location of particles at different travel times under different forces and seasons is another compelling part of the present study.
METHODS
The MIKE 21 FLOW MODEL–FM, developed by DHI Water and Environment (DHI, 2016), was employed to simulate the hydrodynamics (e.g., currents, surface elevations, volume flux) of the study area when forced with realistic wind, tide, and freshwater flow. The two-dimensional model solves shallow-water, the depth-integrated incompressible Reynolds averaged Navier-Stokes equations (DHI, 2016). The MIKE model has been successfully implemented for various coastal setting(s), particularly for the northern Gulf of Mexico (Allahdadi et al., 2017; Allahdadi et al., 2011; Freeman et al., 2015). These studies successfully analyzed coastal and nearshore circulation dynamics along the Louisiana coast. The MIKE software suit enables the hydrodynamic model to be coupled with Lagrangian-based particle tracking and agent-based modules to simulate the transport of passive and active particles, respectively (DHI, 2016).
Computational Mesh
MIKE's flexible mesh (FM) technology (DHI, 2016) was used to construct the computational mesh, comprising 20328 nodes and 33875 triangular elements (Figure 2A,B). The model domain encompasses Gasparilla Sound, Charlotte Harbor, Pine Island Sound, Matlacha Pass, the CRE, San Carlos Bay, Estero Bay, and all of the major tributaries (see Figure 1) and extends 80 km offshore into the Gulf of Mexico to allow tides to rhythmically propagate from the outer shelf into the shallow inlets and bays (Figure 2A,B; Zheng and Weisberg, 2004). The FM is categorized into five distinct regions, depending on the preferred level of resolution. Finer mesh resolutions provide more precise simulation results at the expense of increased computational time (Davids et al., 2010). The finest resolution (∼250 m) was implemented in the focus area of the forthcoming agent-based oyster larval transport studies within the CRE and San Carlos Bay. The resolution is reduced away from the focus region, with values of 700 m, 1.8 km, 2.3 km, and 2.7 km for the inner estuary, offshore, middle offshore, and outermost region, respectively.
Model Inputs
Bathymetric data were obtained from the South Florida Water Management District (SFWMD) and the General Bathymetric Chart of the Oceans (GEBCO). The SFWMD provided high horizontal–resolution (30 m) bathymetry measurements throughout the estuary and around 13 km offshore. The provided data are a compilation of LIDAR data sets from multiple missions flown by the United States Army Corps of Engineers and the United States Geological Survey (USGS) as well as high-resolution bathymetry-topography surveys conducted by the Florida Department of Environmental Protection and other supporting agencies. To supplement the model domain extending ∼80 km offshore, additional bathymetric data, spatial resolution of 30 arc-second (∼900 m), were obtained from the GEBCO online database (GEBCO, 2019).
The model was also forced by wind. Wind data were collected from the National Centers for Environmental Predictions North American Regional Reanalysis online database (NCAR, 2019). Reanalyzed 32-km resolution archived wind vectors and pressure data (three hourly) at 10 m above ground were extracted from the National Center for Atmospheric Research servers (NCAR, 2019). The grib-formatted data were extracted and reformatted into MIKE model–compatible ASCII files using in-house Matlab routines (Jose, Kobashi, and Stone, 2007).
The model domain comprises two boundaries, land and open boundary (water). The offshore water boundary extends in a rectangular shape from the northern to southern coastline of the domain (Figure 2). Offshore water level boundary conditions (tide) were extracted from a global tide model using the MIKE 21 Toolbox. The tide prediction of heights tool generated a time series of predicted tidal elevations for the model's offshore boundaries (e.g., Figure 3; DHI, 2016).
Figure 3
Times series of offshore water level (m) boundary conditions generated using MIKE 21 tide prediction Toolbox.
River discharge data were collected from the United States Geological Survey (USGS) National Water Information System (USGS Water Information System, 2019) at seven sites throughout the estuarine system (Figure 1; Figure 4; Appendix Table A1) and were included as inputs into the model. As a result of the high freshwater flow, the Caloosahatchee River discharge was introduced as two adjacent point sources, with each location releasing half of the total measured discharge.
Bottom Friction Coefficients
Georeferenced characteristics of the estuarine and nearshore bottom sediments and substrate features were collected from the USGS usSEABED online database (Buczkowski et al., 2006; USGS, 2019; Figure 5). The usSEABED database is a collection of both published and unpublished data from a variety of sources and provides point source data regarding various benthic characteristics (e.g., sediment texture and composition, seafloor hardness) (Buczkowski et al., 2006).
Figure 5
Spatial distribution of seabed roughness parameter data points available within model domain. Data compiled from usSEABED database maintained by USGS.
Particle grain size (phi converted to mm) data obtained from usSEABED enabled the approximate calculation of seabed roughness parameter in terms of Manning's roughness coefficient by using the Strickler equation 1 (Chow, 1959):
n = Manning's roughness coefficient
d50 = median grain size (mm).
Numerical Model Implementation
For the model calibration, simulations were performed for three time periods (seasons): 1 May–15 June 2016 (wet season), 15 November–30 December 2016 (dry season), 1 February–20 March 2016 (El Niño—an unusually wet, dry season, winter period). The wet and dry season time periods (Sun, Wan, and Qui, 2016), having distinct hydrologic conditions, were chosen to ensure that the model could accurately simulate contrasting seasonal hydrodynamics. Additionally, an atypical dry season winter period comprising unusually high freshwater discharges (Figure 4; Appendix Table A1) provided an opportunity to model a period commonly associated with low discharges. The atypical period resulted from the strong 2016 El Niño-Southern Oscillation (ENSO) event (Santoso, McPhaden, and Cai, 2017). The ENSO events create strong responses in dry season precipitation and related river discharge in south Florida as compared to neutral ENSO conditions (Hanson and Maul, 1991; Schmidt et al., 2001).
Model Calibration Data
Predicted tidal amplitudes and measured water-surface elevations obtained at multiple stations throughout the model domain were used to calibrate and validate the model.
Water-level data within the study area were collected from two online databases: the University of South Florida's Coastal Ocean Monitoring and Prediction System (COMPS, 2019) and the National Oceanic and Atmospheric Administration (NOAA) Tides and Currents (NOAA 2019; Table 1). The various stations provide sufficient spatial coverage of the model domain (Figure 1; Table 1), thus ensuring the model's accuracy throughout the entire domain.
Table 1
Stations with measured or predicted tides, water levels, and currents.
Additionally, data collected from an upward looking Acoustic Doppler Current Profiler (ADCP) deployed at Big Carlos Pass (station 1A; Figure 1) from 5 May–9 June 2016 enabled model validation of u and v current velocities. Measured currents averaged across the water column were used for comparison with the depth-averaged model results.
Model Calibration and Sensitivity Analysis
Bottom friction coefficients were tuned to determine the best agreement between the simulated and measured currents, water-surface elevations, and tidal amplitudes. A range of Manning's spatially uniform roughness coefficients (15–60 m1/3/s) and the spatially varying Manning's roughness coefficients calculated from the usSEABED measurements were tested to optimize model performance. Although usSEABED provided in situ seafloor measurements (Figure 5), overall model performance was greatest with a spatially constant Manning's roughness coefficient of 60 m1/3/s.
Model Performance Evaluation
Model performance was skill assessed by comparing simulation outputs with in situ current velocities from Big Carlos Pass (station 1A; Figure 1) and water-surface elevations at two stations (stations 1B and 2 in Figure 1) within the estuarine system. Current velocity measurements were collected from 5 May–9 June 2016; therefore, simulated current velocities were validated for only the wet season simulation. Additionally, predicted tidal amplitudes were obtained for the remaining stations (Figure 1). The predicted tidal amplitudes are based on the entire range of tidal harmonic constituents (37 constituents) from the St. Petersburg, Florida (Station ID 8726520) and Naples, Florida (Station ID 8725110) NOAA tidal stations. The UTide MATLAB Function (Codiga, 2011) was employed to filter tidal amplitudes from the simulated water-surface elevations to enable direct comparison between predicted and simulated tidal amplitudes.
Four statistical parameters were employed to evaluate model performance in simulating current velocities, water-surface elevations, and tidal amplitudes: root mean squared error (RMSE), coefficient of efficiency (Nash-Sutcliffe coefficient; E; Nash and Sutcliffe, 1970), coefficient of determination (R2), and index of agreement (d; Willmott, 1981):
In which OBSi are the observed values, 
are the mean observed values, SIMiare the simulated equivalent values, and 
are the mean simulated values.
The RMSE is a measure of the differences between the simulated (SIMi) and observed (OBSi) values, with lower values indicating a better relationship between the observed and simulated values. Coefficient of efficiency (E) values range from minus infinity to 1, whereas the coefficient of determination (R2) and index of agreement (d) values range from 0–1, with values closer to 1 indicating a better fit in each test.
RESULTS
Model performance was evaluated by testing various spatially constant Manning's roughness coefficients (15–60 m1/3/s) as well as the spatially varying Manning's roughness coefficients calculated from the usSEABED data. Figure 6 displays the model sensitivity at three stations (9, 7, and 5; Figure 1) spatially located throughout the model domain. At station 9, the Manning's roughness coefficient of 60 m1/3/s produced the best agreement between the modelled and predicted tidal elevation during both wet and dry season simulations. Both Manning's roughness coefficients of 45 m1/3/s and usSEABED based coefficients underpredicted the larger ebbing phase as well as generally underpredicting the flooding phase, but to a lesser extent than the ebbing phase. At station 7, all three Manning's roughness coefficients closely reproduced tidal elevations during both modelled seasons. Additionally, during the wet season at station 5, all three Manning's roughness coefficients accurately reproduced the predicted tidal heights with the Manning value of 60 generally slightly overpredicting the peak flood phase, whereas the other two coefficients underpredicted the larger neap phase. The Manning's roughness coefficient of 60 m1/3/s most accurately recreated tidal amplitudes at station 5 during the dry season as the two other coefficients generally underpredicted tidal amplitudes. Although not provided, sensitivity testing occurred at all other sites throughout the domain (Figure 1), and a spatially constant Manning's roughness coefficient of 60 m1/3/s was found most effective at simulating the hydrodynamic conditions in this shallow estuarine system.
Figure 6
Comparison of predicted and simulated tidal elevations (m) for varying Manning's roughness coefficients for the wet and dry seasons.
Model performance results for the three model simulation periods (seasons), each associated with distinct hydrologic conditions, with a spatially constant Manning's roughness coefficient of 60 m1/3/s are shown in Table 2. The wet season simulation slightly outperformed the dry season and winter period simulations; however, good overall agreement occurred between simulated and predicted tidal elevations at most locations throughout the system. In particular, the model is capable of reproducing tidal elevations at passes (e.g., stations 6 and 7), mid-estuary (e.g., stations 5 and 8), and upper reaches of the system (e.g., station 9).
Table 2
Model run performance statistics for tidal heights.
To further test model performance, measured and simulated water-surface elevations were compared at two stations: station 2 located in the upper reaches of the CRE and station 1B located within Big Carlos Pass (Figure 1). The model was capable of simulating water-surface elevations at station 1B during the wet season and showed a lowered, yet still adequate performance during the winter period. The model generally overpredicted water levels at Big Carlos Pass (station 1B), especially during spring tides. Water levels at the station located in the upper reaches of the Caloosahatchee River were less successfully modelled throughout all simulations (station 2; Table 3) as water levels were predominantly underpredicted by the model during all simulations.
Table 3
Model performance statistics for water levels for the wet and dry seasons and El Niño–winter period model run.
Model skill was also assessed by comparing measured and modelled u and v current velocities at Big Carlos Pass (Figure 1; station 1A). The model acceptably simulated u and v current velocities at station 1A throughout the wet season simulation (Table 4; Figure 7). Both simulated u and v velocities are generally of lesser magnitude than measured velocities. The model was more capable of reproducing the ebbing phase magnitudes and less successful at reproducing the peak flooding phase; however, both tidal phases were still under-simulated by the model (Figure 7). The model's sufficient performance in reproducing tidal heights, water levels, and currents allowed for the investigation of different aspects of estuarine dynamics, in particular the amount of flow through the system's various passes.
Table 4
Model performance in simulating current velocities at station 1A–Big Carlos Pass from 5 May 5–9 June 2016.
Figure 7
Comparison of measured and simulated currents (m/s) at Big Carlos Pass (station 1A, Figure 1).
Flow through Inlets and Passes
The total flow entering and exiting the system during ebbing and flooding phases was calculated to determine which, if either, phase dominates tidal exchange as well as any possible inequality between the magnitudes of spring and neap tidal flow. Each model simulation covered three consecutive spring and neap tidal cycles (El Niño—winter period: 8 February–15 March; wet season: 7 May–12 June; dry season: 29 November–20 December), thus allowing for the estimation of the total flow (m3/s) during three spring/neap tidal cycles at various passes/inlets throughout the estuarine system. The system's mixed, semidiurnal tidal cycle results in two high and low tides, with one high and low tide having greater magnitude than the other two phases. Flows through the system's various inlets and passes (Appendix Figure A1) were calculated by summing the total discharge during both flooding phases and ebbing phases separately throughout each of the three spring/neap tidal cycles. Total flows during both flooding and ebbing phases were then averaged over the three tidal cycles (Figure 8).
Figure 8
Average simulated spring and neap tide water fluxes (m3/s) through the main inlets of the Charlotte Harbor estuarine system, with error bars indicating standard deviation.
Results indicate that average flows were similar for the three model periods (seasons) at Boca Grande Pass, Captiva Pass, San Carlos Bay, and Big Carlos Pass (Figure 1) during both spring and neap tides but with less water flowing through the passes during neap tides. Boca Grande Pass and San Carlos Bay flush ∼20,000–30,000 m3/s more during the spring tides than the neap tides, whereas Captiva Pass and Big Carlos Pass flush around 5000 m3/s more during spring tides as opposed to neap tides.
During spring tides at Boca Grande Pass, the ebb tidal phase, on average, flushed more water than the flooding phase: ∼7000, 3000, and 9000 m3/s during the wet season, winter period, and dry season simulations, respectively. However, ebb phase dominance did not exist at Boca Grande Pass during neap tides for each simulation, as the flood phase flushed ∼2000 m3/s more during the wet season simulation, whereas the ebb phase flushed ∼5000 and 1400 m3/s more for the winter period and dry season simulations, respectively. Flood phase fluxes almost entirely dominated flushing at San Carlos Bay throughout both spring (∼3500, 5000, and 2500 m3/s wet season, winter period, and dry season) and neap (∼3000 and 1000 m3/s winter period, dry season) tides. However, the ebb phase flushed ∼1800 m3/s more than the flood phase during the wet season, neap tide at San Carlos Bay.
Fluxes at Captiva Pass were quite similar between spring/neap tides and flood/ebb tidal phases with the greatest average difference between tidal phases (∼1000 m3/s) taking place during the neap tide, winter period simulation. The differences between the remaining two simulations were less than ∼50 m3/s. At Big Carlos Pass, the ebb phase mostly dominated flushing (∼100–200 m3/s) for all simulations apart from the neap tide, wet and dry season, simulations where the flood phase flushed (∼150 m3/s) more.
Overall, neither the flood nor ebb phase exhibits uniform dominance in flushing at the four passes/bay. Boca Grande Pass and Big Carlos Pass were mostly ebb dominate, San Carlos Bay was largely flood dominate, and neither phases substantially dominated at Captiva Pass.
DISCUSSION
The nature and extent of substrate plays a critical role in controlling flow dynamics, especially as a friction factor. Soft sediments and materials provide less frictional drag on flow, whereas the roughness and drag of the bottom is greater for coarser materials and for benthic habitats (Hart and Finelli, 1999; Leonardi, Orlandi, and Antonia, 2007). Although usSEABED measurements were available throughout the model domain (Figure 5), simulation performance was best with a spatially uniform bottom friction coefficient. These results suggest that usSEABED measurements may be too spatially sparse in regions of less adequate model performance, viz., upper reaches of the Caloosahatchee River. The sparse data points in the Caloosahatchee River region likely introduced inaccuracies during the interpolation of the usSEABED measurements into the model domain, thereby possibly resulting in the weaker performance.
In addition, it is important to acknowledge the challenges associated with accurately modeling the upper reaches of the CRE. Zheng and Weisberg (2004, 2012) attributed difficulties in modeling the upper portion of estuary to complexities in the estuary's long, narrow shape and shallow depths. The model is capable of simulating the water-surface elevation phases (i.e. tidal phases) in the upper reaches of the estuary (station 2; Figure 1); however, the model generally underpredicted water-surface elevations. Simulation inaccuracies may have resulted from use of a spatially uniform bottom friction coefficient throughout the domain, which may not be representative of the CRE bottom substrate and the spatial distribution of various benthic habitats. Also, channelized flow within the estuarine system may not be fully captured in the simulation because of the spatial resolution of the estuarine bathymetry data and model resolution. This could be tested by implementing finer grid resolution, ∼100 m instead of the current ∼250 m within the estuary. However, the increased computational time required for the finer resolution may not be conducive to this current study, but likely necessary for the future larval transport study.
Previous modeling studies by Goodwin (1996) and Zheng and Weisberg (2004) found the estuarine system to be slightly ebb dominate (i.e. discharge is greater during ebb tidal phase) but with less ebb dominance throughout neap tides. The results are consistent with the previous findings of ebb dominance at Boca Grande Pass; however, the current model simulation results found San Carlos Bay to be flood dominant in respect to flushing as compared to Goodwin's (1996) measured (moving boat technique) discharge rates indicating ebb dominance. The main reason for these inconsistencies could be the circumstances under which the input and output fluxes to the inlets were calculated. The circulation, and thereby water exchange between different parts of the estuary-bay system, is produced as a result of different sources and forces, including freshwater discharge from rivers, tidal force, and wind stress. Hence, depending on the strength of each of these forces, mass balance over the system can be altered, and a different behavior in terms of tidal force dominancy can occur. The inconsistency could also be the result, or combined results, of inaccurate bottom roughness values (particularly within the Caloosahatchee River) and differing bathymetry (i.e. 1996 study bathymetry vs. 2007 LIDAR survey model bathymetry) resulting from the shifting sand shoals and dredging activities. The transient nature of the sand shoals was discovered as a factor in the amplification of tidal currents (de Swart and Zimmerman, 2009) at the inlet where the ADCP was deployed for current velocity measurements; therefore, it is possible that longer term, large-scale sediment transport has taken place throughout the estuarine system, in particular within San Carlos Bay.
The Charlotte Harbor estuarine system is a tidally driven microtidal estuary; however, shallow estuaries can be strongly influenced by wind forces (Cho, 2007), freshwater releases, and other meteorological-related events. This is especially important because freshwater inflows within the estuarine system are seasonally variable. Particularly, within the Caloosahatchee River portion of the system where regulated, unnaturally high and low freshwater flows have been identified as a key ecological stressor. High magnitude Caloosahatchee River freshwater flows, exceeding 130 m3/s, have been found to physically flush planktonic organisms entirely out of the system whereas extended freshwater releases, resulting in reduced salinities, can cause the mortality of benthic organisms (Chamberlain and Doering, 1998).
To better understand the affects freshwater flows impose on planktonic organisms within the system, a particle tracking study was performed.
Relative Dominance of Tidal Inlets
Freshwater inputs into the estuarine system were predominately greater during the El Niño–winter period simulation than the wet season and dry season, apart from the Myakka and Peace Rivers, which experienced similar El Niño–winter period and wet season discharges (Figure 4; Appendix Table A1). Although fluxes through the inlets and passes were similar during the three model simulations, high freshwater inputs into the estuarine system have the ability to modify the system's exchange flow (i.e. flux of fresh and salt water) and the residence time of various estuarine components (e.g., nutrients, sediment, larvae). To further study the effect of freshwater inputs on the system's residence time, a Lagrangian particle tracking study was performed.
One thousand neutrally buoyant particles were simultaneously released on the water surface (depth 0.0 m) at two point sources located in the upper reaches of the CRE (26.69° N, 81.83° W) and the Peace River (26.96° N, 82.04° W) during a spring tide's peak flood phase in each simulation (El Niño–winter period (21 February 16), wet (20 May 16), and dry (20 November 2016) seasons) (R, Figure 1). The two release locations were chosen to compare and contrast regions with natural (Peace River) and unnatural, regulated freshwater releases (CRE). Particle locations 5, 10, 15, and 20 days postrelease are provided in Figure 9.
Figure 9
Spatiotemporal distribution of neutrally buoyant particles released into the estuarine system from Peace River and Caloosahatchee River for the three seasons representing three contrasting hydrologic conditions. Panels on the left represent Peace River simulations; panels on the right represent Caloosahatchee River simulations.
Particles released in the Peace River during the El Niño–winter period simulation largely remained in the northern reaches of the estuarine system at the conclusion of the 20-day simulation (Figure 9). A portion of particles were transported westward into the Myakka River region (Figure 9, e.g., Day 10) and then continued to travel southward along the western shoreline, whereas a moderate number of particles travelled southward along the eastern shoreline. This division of particles (Figure 9, Day 15) is likely attributable to a deeper channel (∼5–6.5 m) located within the center of the northern estuarine system, as compared to the shallower depths (∼2–3.5 m) away from the channel, near the shorelines. Although three particles were flushed out of the system through Boca Grande Pass after 20 days, a majority of particles remained within the northern region of the system without being dispersed southward into Matlacha Pass or Pine Island Sound. In contrast to the Peace River particles transport, particles released at the CRE location were quickly flushed out into the Gulf of Mexico through San Carlos Bay, 10 days postrelease (Figure 9, Day 10). This is likely the result of the Caloosahatchee River's high freshwater inputs (∼187 m3/s; Appendix Table A1) as well as winds predominately blowing offshore during the simulation period (Figure 10A; Appendix Table A2). Additionally, particles were transported northward into Matlacha Pass and Pine Island Sound and ultimately flushed into the Gulf of Mexico through various passes (e.g., Redfish Pass, Boca Grande Pass).
Figure 10
Wind rose diagrams for the three simulation periods, El Niño–winter period (A), wet season (B), and dry season (C) at W—Figure 1 (26.50° N, –82.35° W). Values in the circles represent the percentage of calm corresponding to wind speeds smaller than 3 m/s.
During the wet season simulation, particles released at the Peace River location exhibited comparable, yet contrasting transport characteristics to the El Niño–winter period simulation. The particles were once again moved into the Myakka River region; however, particles were not divided between the west and east shorelines but were solely transported along the western shoreline (Figure 9, Days 15, 20). Also, all particles remained in the northern region with none being flushed into the Gulf of Mexico. In comparison to the CRE El Niño–winter period simulation, the relatively similar high Caloosahatchee River discharges during the wet season simulation also flushed particles out through San Carlos Bay into the Gulf of Mexico, but with a greater number of particles traveling northward into Matlacha Pass and continuing further into the northern portion of the system. This may have also resulted from the relaxed wind stress during the summer wet season period (Figure 10B; Appendix Table A2), which allowed the tidal forces to more effectively transport particles into the northern region.
Once again, particles released in the Peace River during the dry season simulation largely remained within the northern region after 20 days (Figure 9, Day 20). However, a higher abundance of particles persisted within the uppermost northeast portion of the system, likely the result of the low Peace River inputs (∼5 m3/s; Appendix Table A1) during the simulation. Similarly, CRE-released particles stayed within the system, mainly within the CRE section; however, a minor portion of released particles were transported into Matlacha pass. These results can likely be attributed to the low Caloosahatchee River inputs during the simulation even with the predominately offshore winds (Figure 10C; Appendix Table A2).
The Lagrangian particle tracking study provides valuable insight into the effects that freshwater inputs and wind forces impose on different regions of the estuarine system. Freshwater inputs play an important role in the circulation and mixing within the Charlotte Harbor estuarine system and, as provided in Figure 9, in how the two regions (Peace River and CRE) respond differently to the seasonal hydrologic forcing. Particles released in the Peace River were comparably transported to similar regions of the estuarine system regardless of the amount of freshwater input and predominate wind direction. This occurs apart from less particles being transported away from the release region during the dry season's lower freshwater inputs, which is to be expected. On the other hand, CRE-released particles were transported differently depending on the freshwater inputs and predominate wind direction. During the high freshwater inputs (El Niño–winter period and wet season), particles were flushed out into the Gulf of Mexico within 10–13 days depending on the simulation; however, the predominant wind direction either equally dispersed particles throughout the entire system (El Niño–winter period) or facilitated more particle transport into Matlacha Pass and furthermore into the northern section of the system (wet season). Similar to the dry season Peace River particle transport, particles largely remained in the CRE section during the dry season's low freshwater input, which is to be expected.
CONCLUSION
The developed hydrodynamic model is capable of simulating water levels, tidal heights, and currents with good accuracy for the Charlotte Harbor estuarine system, apart from the weaker performance in the upper reaches of the Caloosahatchee River. Model performance was tested with varying Manning's roughness coefficients with the spatially uniform Manning's roughness coefficient of 60 m1/3/s providing the strongest match to the predicted and measured data. If additional seabed measurements, allowing for the calculation of bottom friction, could be obtained within the CRE, those measurements, in conjunction with the spatially varying Manning's roughness coefficients provided by the SEABED data, may produce more accurate model results within the CRE. Model performance was also consistent during varying hydrologic conditions (i.e. El Niño–wet period, wet and dry seasons), which supports the model's ability to simulate the seasonally varying conditions within the system. The model's ability to simulate the system's estuarine dynamics allowed for the calculation of the estuarine system's flushing rates at various passes.
Flushing rates calculated at various passes indicated that the estuarine system is neither flood nor ebb phase dominate, but rather a combination of the two phases. Boca Grande Pass and Big Carlos Pass were found to be mainly ebb dominant, San Carlos Bay was largely flood dominant, whereas Captiva Pass was found to be neither flood nor ebb dominate. The varying tidal dominance in flushing the system is important as the system is also impacted by local wind directions. For example, predominately offshore-ward winds and high freshwater inflows during the El Niño–winter period reduced the CRE particle travel time as compared to wet season's (Figure 9) shoreward winds, yet with both still experiencing high magnitude freshwater inflows. Because relatively similar tidal forcing was applied for each of the three different modeling scenarios, the difference in distribution of particles and, thereby, their travel time for each case could mostly be attributed to the other two forcing, i.e. freshwater discharge and wind.
The particle transport model outputs from the CRE for high freshwater input scenarios provides interesting results and valuable guidance in the pursuit of the oyster larval transport modeling within Charlotte Harbor Estuarine System. Because oysters spawn throughout the year in SW Florida (Volety et al., 2009), during the wet season, oyster larvae could be transported into Pine Island Sound and Matlacha Pass and may find refuge in these sheltered water bodies during high freshwater discharge. Whereas a greater abundance may be flushed out of the system during El Niño years with high freshwater input conditions (e.g., 2016 El Niño event). Incorporating previous field-based studies (e.g., Chamberlain and Doering, 1998; Volety et al., 2009) within the system along with particle tracking will help to gain a better understanding of the impacts freshwater inflows impose on the system's planktonic organisms, in particular oyster larvae.
ACKNOWLEDGMENTS
The authors acknowledge Dr. David Fugate, Florida Gulf Coast University, for providing the ADCP data from Big Carlos Pass and the South Florida Water Management District for providing the bathymetry data. Additional appreciation and gratitude to DHI for supporting our research through the use of their MIKE software as well as DHI scientists Jonas Brandi Mortensen and Mikkel Anderson for their countless support in the development and implementation of the hydrodynamic and particle tracking model.
