*2.2. Data Analyses*

Results of instantaneous discharge measurements during July–October 2009 are reported elsewhere [12]. Most spring sites were measured under more than one tidal condition, resulting in multiple discharge samples for these sites. At 11 sites, spring discharges were measured only one time because of their relatively low flow rates, though multiple readings were recorded for each of them to either obtain an average flow rate or get the sum for the spring site if several vents are involved. Similar to the finding reported in previous measurements, the magnitude of the spring flow in Kings Bay varies greatly from one vent to another. About 15 spring sites are second order magnitude springs with mean discharges ranging between 10 and 100 cfs, while two are fifth order magnitude (1–100 gal/min, or 0.0223–0.223 cfs) or less. For the same spring site, the discharge also varies significantly from time to time. For example, a spring vent named H24 in the north portion of Kings Bay was measured on two different days. One was on 23 September 2009 and the other was on 7 October 2009. The first measurement of the discharge was 8.35 cfs, but the second flow measurement was 49.5 cfs, almost six times of the previous measurement. The total of measured mean flows from all the identified vents is about 467 cfs.

Generally, salinity is higher in southern springs than in northern springs. Except for Sites No. 1 and H24 (Figure 2), most northern springs discharge fresh water, with salinity normally less than 0.5 psu. Site No. 1 is located at the headwater of Miller Creek, which is a short waterway connecting the spring with the Crystal River (Figure 2), and has an average salinity of 1.75 psu. H24 is connected to Kings Bay through a spring run and has an average salinity of 1.14 psu. Southern springs are brackish and salinities in these springs can be 6 psu or higher (e.g., Spring Sites 38–40). Overall, the flow-weighted salinity in spring flows out of all spring vents is about 1.58 psu.

Spring temperature is generally much more stable than spring salinity in Kings Bay, both in terms of special variation and temporal variation. Most spring flows have a temperature around 23.5. The highest spring temperature was measured at H24 Spring Site, with a value of 24.94 °C. The lowest spring temperature was measured at Spring Site No. 1 located at the headwater of Miller Creek, where the average spring temperature was 22.93 °C. The flow-weighted spring temperature was 23.51 °C for all the springs in Kings Bay during the measurement period. During the coldest days in winter, spring temperature can be about 1 °C lower than in summer.

As mentioned above, real-time cross-sectional fluxes at G1 and G2 shown in Figure 2 were measured with multi-beam ADCPs during 27 July to 20 August 2009. Because cross-sectional fluxes measured at the two sites also include tidal prisms upstream of the cross sections, net spring discharges from the two groups of springs need to be adjusted as follows

$$q\_{\mathcal{R}} = q\_{\mathcal{m}} + A(\hat{\mathcal{O}}\eta/\hat{\mathcal{O}}t) \tag{1}$$

where *qg* is the net spring flow of the group, *qm* is the cross-sectional flux measured by ADCP, *A* is the total water surface area upstream of the cross section, *t* is time, and *Ș* is the water surface elevation measured at the mouth of Kings Bay station. *qg* and *qm* are positive leaving the spring group. The second term on the right hand side in the above equation is the flux due to the tidal prism. Using a geographic information system, the total water surface areas upstream of Groups 1 and 2 are found to be 613,111.20 and 1,785,544.49 square feet at the mean sea level, respectively. Because shorelines in both areas are mostly man-made vertical seawalls, both areas vary around their mean sea level values within a very small range and thus can be treated as constants. The time derivative of water surface elevation can be calculated from measured data at the USGS mouth of Kings Bay station.

Results of Equation (1) for both Groups are presented in Figure 3, along with measured water level at the USGS mouth of Kings Bay station during the same period. There were some problems with the ADCP measurement at the G1 cross section after about 9.5 days, and thus only the first 9.5 days of Group 1 spring flow data are plotted in Figure 3.

From Figure 3, one can see that spring flows for both Groups 1 and 2 exhibit strong tidal signals. Mean discharge out of Group 2 spring vents are about twice of that out of Group 1 spring vents; however, the range of discharge variation for Group 2 springs is more than five times of that for Group 1 springs. As shown in Figure 3, net spring flows from both spring groups are negatively proportional to the surface water elevation. As water level increases, spring flows decrease, and *vice versa*. When surface water level in Kings Bay increases to a certain elevation, the net spring flow from Group 2 vents becomes negative. In other words, instead of ground water being discharged out of the springs, estuarine water in Kings Bay flows into these spring vents.

**Figure 3.** Measured spring flows for Group 1 (red short dashed line) and Group 2 (blue long dashed line) during 27 July–20 August 2009. The green solid line is measured water level at the USGS mouth of Kings Bay station during the same period.

As mentioned above, the inverse relationship between tides and SGD has been reported in many previous studies all over the world, including those for Kings Bay springs [2,3]. Based on flow measurement for a single spring (Spring Site 6 in Figure 2), Hammett *et al.* [3] obtained the following linear regression equation

$$q\_{\!\!\!\!-1} = 16.97 - 4.25\eta \tag{2}$$

where *q*1, in cfs, is the estimated spring flow for Spring Site 6 in Figure 2 and *Ș*, in feet NGVD 29, is water level measured at the mouth of Kings Bay station.

Equation (2) suggests that the spring discharge is only a function of water level. This is obviously not the whole story for spring discharges in Kings Bay, because one of the main driving forces that cause springs to discharge flows to Kings Bay, namely the groundwater level, also varies with time. Figure 4a shows measured daily and monthly groundwater level data during 1 January 2005 through 1 September 2010 in a nearby well called ROMP TR21-3 (Figure 1). A comparison of the water level measured at the USGS mouth of Kings Bay station with measured hourly groundwater level is shown in Figure 4b. Although the groundwater level is relatively stable in comparison with the tides in the system, it does have tidal signals in it. From Figure 4, it is clear that groundwater level contains not only high frequency variations but also low frequency variations with a time scale of a year. Hence, the consideration of groundwater level in predicting spring flow is necessary.

**Figure 4.** (**a**) Measured monthly and daily groundwater levels in ROMP TR21-3 during 1 January 2005–1 September 2010. (**b**) Comparison of measured Kings Bay tides and hourly groundwater level in ROMP TR21-3 during the 25-day continuous recordings of cross-sectional fluxes at G1 and G2.

Based on available data, including continuous cross-sectional fluxes at G1 and G2, tides at the USGS mouth of Kings Bay station, and groundwater level data at ROMP TR21-3, it was found that the following linear equation can describe the effects of tides and groundwater level on the spring flow in Kings Bay very well.

$$q = q\_0[1 + C\_1(G - \Delta G - \eta) + C\_2(\partial \eta / \partial t)]\tag{3}$$

where *q* denotes the estimated spring flow, *q*0 is the long-term mean spring flow, *G* represents the groundwater level in ROMP TR21-3, ǻ*G* is the long-term mean head difference between the groundwater level in ROMP TR21-3 and the surface water level in Kings Bay, and *C*1 and *C*2 are two parameters which are time-independent and can be determined from measured field data. The time derivative of surface water elevation (*Ș/t*) in the above equation not only allows the phase mismatch predicted by the head difference to be eliminated, but also allows higher mode oscillations shown in the measured net spring flows to be correctly matched.

Figure 5 compares time series of measured spring flows with those estimated using Equation (3). As can be seen from the figure, Equation (3) predicts spring flows very well, especially for Group 2 springs. The R<sup>2</sup> values for the match of estimated and measured SGDs are 0.72 and 0.94 for Groups 1 and 2 springs, respectively.

**Figure 5.** Time series of estimated and measured spring flows at G1 for Group 1 springs (**a**) and G2 for Group 2 springs (**b**) during 27 July–20 August 2009.

The three parameters (ǻ*G*, *C*1, and *C*2) in Equation (3) were determined through a trial and error process in obtaining the best match between estimated and measured spring flows. This trial and error process yielded two sets of (ǻ*G*, *C*1, and *C*2): (70.3 cm, í0.0088 cm<sup>í</sup><sup>1</sup> , 3.9 s cm<sup>í</sup><sup>1</sup> ) for Groups 1 and (67.06 cm, í0.0166 cm<sup>í</sup><sup>1</sup> , 79.78 s cm<sup>í</sup><sup>1</sup> ) for Group 2. Clearly, for different spring groups, these parameters are quite different, especially *C*1, and *C*2.

#### **3. Verification of SGD Estimates via Hydrodynamic Simulations**

The empirical formula, Equation (3), for estimating SGDs from the spring vents was used in an unstructured Cartesian grid hydrodynamic model named UnLESS3D, which was developed to simulate circulations, salinity transport, and thermodynamics in the Crystal River/Kings system. Figure 6 shows the unstructured grid mesh used by the UnLESS3D model for the estuarine system, which was discretized with 3030 horizontal grids and 14 vertical layers. The model was driven by measured water elevations, salinities, and temperatures at open boundaries (USGS stations near Shell Island and in Salt River) and wind shear stresses and heat flux at the water surface, which were calculated based on measured wind, solar radiation, air temperature, and relative air humidity at a weather station about 10 miles north of Kings Bay. The UnLESS3D model was also driven by the spring flows at the bottom of Kings Bay. The following describes how the spring flows from all the spring vents are estimated based on data measured during 2008–2009 and how they are incorporated in the model simulation. Comparisons of model results and measured real-time data are presented to demonstrate that the spring estimates are reasonable. Details on the theory of the hydrodynamic model can be found in [13] and its application to Crystal River/Kings can be found in [14].

**Figure 6.** Unstructured Cartesian grid mesh used in the model application to the Crystal River/Kings Bay estuary.

Using the mean spring discharges during the 2008–2009 field investigation as the long-term mean flow rates, Equation (3) was used to estimate the spring flow from each individual vent at each time step of the simulation during the model run. For those vents in Groups 1 and 2, the two sets of parameters (ǻ*G*, *C*1, and *C*2), determined through the trial and error process mentioned in the last section, were used. For other spring sites, it is further assumed that their spring flow characteristics are similar to those of either Group 1 springs or Group 2 springs, depending on the distances from the spring site to the two spring groups. For example, if the distance from a spring site to Group 1 springs is shorter than that to Group 2 springs, the set of parameters for Group 1 is used for that spring site; otherwise, the parameters (ǻ*G*, *C*1, and *C*2) for Group 2 will be used.

Another SGD to Kings Bay is hairline fractures which are spread all over the bottom of Kings Bay. In order to include this spring flow source, 40 small spring vents (one can use more) were assumed to be randomly distributed on the bottom of Kings Bay according to a uniform distribution. The long-term mean discharges of the 40 small vents were also randomly assigned with values between 0 and an assumed maximum value (e.g., 1 cfs). Because there is no way to measure flows out of all the hairline fractures in Kings Bay, a factor (Į) was used in the model to adjust the magnitudes of spring discharges from these assumed spring vents which represents the hairline fractures. The determination of this adjustment factor is a part of model calibration.

The UnLESS3D model was calibrated and verified against measured real-time water level, salinity, and temperature data at the USGS Bagley Cove and the mouth of Kings Bay stations within the simulation domain (Figure 6). The total simulation period was a 34-month period (1037 days), from 24 April 2007 to 23 February 2010. The model was calibrated against real-time data for a 150-day period during 28 December 2007–26 May 2008 after a spin-up run of 25 days. It was then verified for the remaining days before and after the 150-day calibration period.

During the model calibration process, four model parameters, including the bottom roughness, the background eddy viscosity/diffusivity, the attenuation coefficient of short wave radiation, and the flow adjustment factor (Į) for the hairline fracture, were tuned to achieve the best agreement of model results with field data of water level, salinity, and temperature.

**Figure 7.** Comparison of measured and simulated water levels at the Bagley Cove station (top panel) and the mouth of King Bay station (bottom panel) during 2 August–1 October 2007.

Comparisons of model results with measured field data at the two measurement stations inside the simulation domain (Bagley Cove and the mouth of Kings Bay) are presented in Figures 7–9. For simplicity and clarity, only an arbitrarily chosen, two-month period between Hour 66,480 and Hour 67,920 (2 August–1 October 2007) are shown here to demonstrate how model results are compared with measured field data. Figure 7 shows and compares simulated and measured water levels at the two measurement stations. The top panel of the figure is for the Bagley Cove station, and the bottom panel is for the mouth of Kings Bay station. Dashed lines are measured water level data, while solid lines are simulated water level results by the UnLESS3D model. Figure 8 shows and compares simulated and measured salinities. The top panel is for the bottom layer at the Bagley Cove station (salinity and temperature were only measured near the bottom at this station), while the middle and bottom panels are for the top and bottom layers, respectively at the mouth of Kings Bay station. Again, dashed lines are measured field data and solid lines are model results. In Figure 9, simulated and measured temperatures at Bagley Cove and the mouth of Kings Bay are shown and compared in the same manner as that in Figure 8, with dashed lines being field data and solid lines being model results.

**Figure 8.** Comparison of measured and simulated salinities near the bottom at the Bagley Cove station (top panel), in the top layer at the mouth of King Bay station (middle panel), and in the bottom layer at the mouth of Kings Bay station (bottom panel) during 2 August–1 October 2007.

**Figure 9.** Comparison of measured and simulated temperatures near the bottom at the Bagley Cove station (top panel), in the top layer at the mouth of King Bay station (middle panel), and in the bottom layer at the mouth of Kings Bay station (bottom panel) during 2 August–1 October 2007.

From Figures 7–9, it can be seen that the modeled water levels, salinities, and temperatures agree well with measured real-time data at the two USGS stations within the simulation domain. For the entire 34 months of the simulation period, the overall mean errors for water level, salinity, and temperature are 0.90 cm, í0.06 psu, and 0.03 °C, respectively, while overall mean *R*<sup>2</sup> values for water level, salinity, and temperature are 0.98, 0.75, and 0.90, respectively. The mean skill assessment parameters of Willmott [15] are 0.99 for water level, 0.91 for salinity, and 0.98 for temperature. Detailed skill metrics of the model performance can be found in [14].

Because freshwater loading is one of the most important factors controlling circulations, salinity distributions, and thermodynamics in an estuary, a good match of model results with field data is impossible without a good quality of freshwater loading data that are used in the model. As mentioned above, 99% of fresh water loading to Kings Bay comes from the spring flows. As such, a good match between model results and measured data shown in Figures 7–9 not only suggests that the UnLESS3D model works well in simulating circulations, salinity transport processes and thermodynamics in the Crystal River/Kings Bay estuary, but also indicates that the above procedure in estimating SGD from the spring vents at the bottom of Kings Bay is adequate. In other words, the 3D modeling using the UnLESS3D model implies that the empirical formula expressed in the form of Equation (3), which links SGD in Kings Bay with tides and groundwater level, is applicable to all the spring vents in Kings Bay.

#### **4. Summary and Conclusions**

In an effort to quantify spring flows entering Kings Bay and to study how these submarine groundwater discharges are influenced by tides and groundwater level, a two-phase field investigation was conducted during 2008–2009. The first phase was to carry out an inventory study to search, identify, label, and measure the dimensions of all the detectable spring vents in Kings Bay, while the second phase included measuring instantaneous flow rates out of all the detectable spring vents and gauging cross-sectional fluxes at two spring runs continuously for about 25 days. The inventory study identified about 40 additional springs which had never been documented in any previous studies of the system.

Field data of instantaneous measurements of spring flows showed that the total mean flow from all of the detectable spring vents was about 467 cfs during the measurement period. This number is lower than those of previous USGS studies. Effects of tides on spring flows can be clearly seen in these instantaneous measurements.

From real-time field data measured at two spring runs and the available water level and groundwater level data, it is found that the spring flow rate is a linear function of the head difference between groundwater level and water surface elevation in the bay. It is also found that the spring flow increases with the increase of the time derivative of the surface water elevation. Based on the analysis of the real-time field data, an empirical formula that relates spring flow rate with water surface elevation and groundwater level is obtained. This spring flow rate formula was used in an unstructured three-dimensional hydrodynamic model to simulate tidal circulation, salinity transport processes, and thermodynamics in the spring-fed estuary. A good match between model results with measured real-time data of water level, salinity, and temperature is achieved and suggests that the empirical formula is appropriate for estimating SGDs out of the spring vents in Kings Bay.

#### **Acknowledgments**

I would like to thank my colleagues at the SWFWMD for their support and encouragements during this study. Michael Sidney Flannery co-managed the field investigation project in Kings Bay with me, which was contracted to VHB, Inc., for whom Gary Serviss was the project manager. Many discussions and meetings with Sid and Gary led to the successful completion of the project in the field.

#### **Conflicts of Interest**

The author declares no conflict of interest.
