Statistical Assessment of Saltwater Intrusion in a Subtropical Estuary

Statistical assessment of saltwater intrusion in a subtropical estuary.

1. Study Site

The Suwannee River is the second largest river system in Florida by its mean annual flow of 298 m3/s and by its drainage area of 25,770 m2. It is 392 km long, and features a low-gradient stream with an average gradient of 7.5 cm per kilometer (Valle-Levinson, 2012; Light, et al., 2002). The drainage basin encompasses the Northern Highlands and the Gulf Coastal Lowlands physiographic regions. The Highlands are generally around 30-60 m above Mean Sea Level (MSL) and include streams, lakes, and ponds. The Gulf Coastal Lowlands are 30 m above MSL and are characterized by low relief, karstic topography and a high connectivity between surface water and groundwater systems (Suwannee River Water Management District, 2005). The drainage of the Suwanee River ends at the Suwanee River estuary, which extends up to 16 km upstream from the mouth (Figure 1). The estuary has a mean water depth of 2.2 m below Mean Sea Level, and it is dominated by semidiurnal tides 0f 0.7m amplitude (Tuckey & Dehaven, 2006). Observations collected during low river discharge proof that tides can reach up to 43 km upstream from the river mouth (McPherson & Hammett, 1991). But no studies have documented saltwater intrusion during strong river pulses nor its seasonality. These varying conclusions, state the requirement of further studies of the salinity gradient and landward extent of salt water at the Suwannee River.

The climate of the region is a mixture of warm temperate and subtropical conditions, with long warm summers and short mild winters (Crandall, et al., 1999). Mean annual temperature in the Florida portion of the watershed is 20°C, with minima and maxima range of 12 ºC to 27 ºC during the year. Precipitation has a bimodal distribution with peaks in March and September, varying accordingly to the latitude. The precipitation ranges from 117 – 152 cm, with the lowest in the upper basin (higher latitudes) and the highest closer to the Gulf coast—at lower latitudes (Suwannee River Water Management District, 2005). In addition to the temporal and spatial variability of rainfall across the basin, the river is also affected by the El Niño Southern Oscillation (ENSO) phenomenon. For example, evidence shows that the larger flood events of 1973, 1984 and 1998 in the Suwannee River were related to “moderate to strong” El Niño events. This is different from La Niña phase, when the river discharge is below average (Tootle & Piechota, 2004).

Figure 1. Regional location of the (a) Suwannee River estuary and (b) location of the CT stations.

Annual discharge from the Suwannee River into the Gulf of Mexico varies between the warm (May - October) and cold (November - April) seasons. Data collected at the Wilcox United States Geological Survey (USGS) station, over the period of 1942-2003, gives a mean discharge of 320 m3/s for the warm season, 255 m3/s for the cold season, and an overall annual mean of 288 m3/s. During low-flow conditions, groundwater from the Upper Floridian aquifer flows into the river and sustains the base flow. In contrast, during high-flow conditions and heavy rainfall when water from the river flows into the aquifers through small cavities and conduits (Crandall, et al., 1999; Katz, et al., 1997).

This study addresses the salinity distribution along the Suwannee River, which until now has only been described qualitatively. {Light, 2002, Hydrology`, vegetation`, and soils of riverine and tidal floodplain forests of the lower Suwannee River`, Florida`, and potential impacts of flow reductions}{Light, 2002, Hydrology`, Vegetation`, and Soils of Riverine and Tidal Floodplain Forests of the Lower Suwannee River`, Florida`, and Potential Impacts of Flow Reductions}{Light, 2002, Hydrology`, vegetation`, and soils of riverine and tidal floodplain forests of the lower Suwannee River`, Florida`, and potential impacts of flow reductions}Some reports (e.g. Light et al. 2002; Tillis, 1999) indicate that the salinity distribution varies temporally and spatially due to tides, river discharge, and winds. Higher upstream salinities occur during different conditions: low river discharge, high tides and during storm surges if the river discharge is low. Regularly salinity decreases upstream and during low tides; and vertically from bottom to surface (Light, et al., 2002). Vertical differences in salinities from bottom to surface implies some vertical stratification (Tillis, 1999), but further conclusions need a higher sampling resolution. This study reduces that gap of knowledge gap and contributes to understand the salt ocean water distribution of in the Suwannee River estuary in over temporal and spatial scales.

1. Methods

Observations of salinity and water level were obtained at six stations along the Suwannee River; recording measurements every 3 minutes for approximately one year. These measurements were used to calculate the upstream excursion of the isohaline of 2g/kg (2) and the water level. River discharge and precipitation values were obtained at the Wilcox Station from the United States Geological Survey (USGS). Wind speeds and directions were available at the Cedar Key Station from the National Data Buoy Center-NOAA, 23 km southwest from the mouth of the Suwannee River.

The sSaltwater intrusion was represented by the position of 2 . It was explained described spatially and temporally with two approaches. The A temporal approach consists of an autoregressive numerical model that predicted 2saltwater intrusion(dependent variable) based on the variability of the forcings (independent variables (forcings): river discharge, water level, precipitation and winds; and the dependent variable (response): isohaline position. For theThe second approach, a spatial approach, fitted the tidally averaged salinity values along the 6 sampling sites to a prescribed function. , iIt was assumed that the longitudinal distribution of salinity followed the shape of a hyperbolic tangent function. The longitudinal distribution of salinity was defined at each tidal cycle by the maximum salinities from the six stations. [V1]Each longitudinal distribution is associated with the coefficient β from the hyperbolic tangent fit. The coefficient β from the hyperbolic tangent fit represents the length scale of the saltwater[V2], which at the same time can be associated to forcings that modify the salinity distribution. This variability of forcings can be explained with a linear regression model. The numerical [V3]model related the salinity distribution (response) to changes in river discharge, water level, precipitation and winds (forcings). Results from the regression model where used to recalculate the hyperbolic tangent fit and compare it with the observed results of the salinity distribution.

1.1Data collection

Two sets of data were used in our analysis. The first data set consisted of salinity measurements and water levels obtained from different types of instruments: Schlumberger Diver recorders (CTD divers) and Sea-Bird SBE26 and SBE37 recorders (CT). The instruments were deployed from Nov 23, 2014 to Nov 22, 2015, a nearly continuous time series with 3-minute sampling intervals. The CTD divers, and the CTs were deployed on the bed near the edge of the river. The second data set consisted of time series of river discharge, precipitation and wind (Figure 2). River discharge and precipitation time series consisted of long-term data from the Wilcox Station 023323500 (29° 35’ 22’’ N, 82° 56’ 12’’ W). River discharge values had a sampling interval of 15 minutes and available from USGS at ( Precipitation data had daily measurements from a rain gage administered by the Suwannee River Water Management District (SRWMD). Raw data are available at ( /realtime/rain-levels.php). In addition, 10-minute wind speeds and directions were obtained from the meteorological station #CDRF1 at Cedar Key, FL, available at (

1.2Data processing

Tidal effects were removed from the river discharge using a low-pass Lanczos filter centered at 30 hours. Wind vectors were decomposed into north and east components. These composites were smoothed with a low-pass Lanczos filter centered at 70 hours to remove inertial effects and periods longer than 2 days. The water level was demeaned before applying a complex demodulation. A complex demodulation extracts the amplitude of any given frequency from a time series. In this case, the water level corresponds to the amplitude of the semidiurnal tide (12.42 hours) obtained with a complex demodulation. More information about complex demodulation is found in Thomson & Emery(2014). Raw data for precipitation was interpolated to the time range of the river discharge, wind velocities, and water level.

Absolute salinity was estimated using conductivity, temperature and pressure data from the CTD divers and CTs using the thermodynamic equation of seawater – TEOS 2010 (McDougall, et al., 2012). Longitudinal salinity distributions were calculated from six stations along the Suwannee River estuary (Figure 1b). The distance from the mouth of the river to the location where the salinity is 2 g/kg would be the isohaline position of X2. The X2 position was based on bottom measurements of salinity per each tidal cycle e.g. Monismith, et al., 2002.

Saltwater intrusion analysis – X2

A first order Markov chain regression model was applied to time series of the isohaline position X2, river discharge, water level, precipitation and wind velocity components to determine the saltwater intrusion along the Suwannee River. Estimates of saltwater intrusion were obtained for the period of Nov 23, 2014 to July 21, 2015. The maximum location (km) of the isohaline of 2 g/kg (X2) was located from a semidiurnal complex demodulation. Time series of river discharge—, water level—, precipitation—, and wind velocity components were averaged for each tidal cycle. The time series of the independent variables (, , , ,) were interpolated to the same time range of the dependent variable—X2. Discontinuities in the time series in between deployments of the CT recorders (Figure 2).

Figure 2. Data set from Nov 23, 2014 to July 21, 2015 used for the Markov chain regression model. a) Smoothed river discharge (m3/s) with a low-pass Lanczos filter centered at 30hrs; b) Water level (cm) from a semidiurnal complex demodulation; c) Precipitation (mm); d) Wind velocity components in north (blue) and east (red) following the oceanographic convention, smoothed with a low-pass Lanczos filter centered at 70hrs; e) Position of the isohaline of 2 g/kg (X2) from a semidiurnal complex demodulation along the Suwannee River.

Using the methodology of Reyes-Merlo et al. (2013) as a guideline, a first order Markov chain model was developed to relate any given X2 isohaline location to its position in a previous tidal cycle and different forcings. Reyes-Merlo et al. related the X2 isohaline location to its location in a previous tidal cycle and to the river discharge, tidal current range and easterly wind velocity averaged in the latest tidal cycle. In this case, instead of using the tidal current, the model considers the water level at the mouth, both wind velocity components, and the precipitation. Because salinity intrusion is controlled to some extent by river discharge (Tillis, 1999), the analysis was subdivided in high and low discharge regimes (Simpson, et al., 2001). The autoregressive model is described by:

(1)

where the coefficients (b, c, d, e, f, g) are calculated with a least square nonlinear regression for high and low discharge regimes, is the isohaline location in a previous tidal cycle (hysteretic term), is the water level, is the river discharge, is the precipitation and are the wind velocity components (e.g. Figure 2). The coefficient a is attributed to the response time of saltwater intrusion to different forcings. In this case, river discharge influences saltwater intrusion more than other forcings{Tillis, 1999, Flow and Salinity Characteristics of the Upper Suwannee River Estuary`, Florida} (Tillis, 1999). Inspection of Figure 2, indicates that the response of X2 depends on river discharge. The time response was obtained with a lagged cross correlation between river discharge and isohaline position, following:

(2)

The numerator has the variations of each variable around their mean (), with one of them lagged, and the denominator is the standard deviation of each variable.

The lag of highest correlation coefficient () indicates the response of saltwater intrusion to changes in river discharge. The time response was obtained with the lag and the corresponding . This response is represented by coefficient a in the autoregressive Markov-chain model. Mathematically is represented by:

(3)

Coefficient a will be constant for high and low discharge regimes. Low and high discharge regimes were separated with a 250 m3/s river discharge. This value was selected after a series of repetitions for river discharges ranging from 200-400 m3/s. Saltwater intrusion was estimated for July 22, 2015 to November 23, 2015 using the Markov-chain model with the corresponding values for the coefficients—b, c, d, e, f, g, for each discharge regime (Table 1). The CT deployments corresponding to August – September and October – November 2015 were used to verify the model projections. Observed and projected values were compared with the root-mean-square error and goodness of fit, reported in the Results section.

Longitudinal saltwater distribution

The longitudinal saltwater distribution was fitted to a hyperbolic tangent function (Pritchard, 1954; Warner, et al., 2005), with maximum salinities at the ocean boundary and minimum values at the landward boundary. Following Warner et al. (2005), the hyperbolic tangent function has the form:

(4)

where is the longitudinal salinity value, is the maximum salinity at the ocean, is the longitudinal distance in kilometers, is a nondimensional parameter that establishes the position of the origin, in this case , and is the length scale of saltwater intrusion at each tidal cycle (Warner, et al., 2005). Maximum values of salinity were obtained at all station for each tidal cycle, and then fitted to the hyperbolic tangent function. Values of were calculated with the observed values of maximum salinity at all stations together. In that sense, the length scale of saltwater intrusion will be represented by one value of per tidal cycle. The salinity at the ocean was taken as a constant of , and the salinity of freshwater at the river boundary. The hyperbolic tangent function was fitted to all deployments, from November 2014 to November 2015.

Maximum values of salinity at each station are marked with circles and the lines represent salinity values along the estuary but using one value of per tidal cycle. A finer resolution along the estuary was used for smoothing purposes (Figure 7). The length scale of saltwater intrusion) accounts for the variability of different forcings that modify the saltwater intrusion. Following this reasoning, will be fitted to a linear regression model that relates the forcings of river discharge, water level, wind velocity components and precipitation, as

(5)

Coefficients m,n,o,q,r,s are calculated with a least square nonlinear regression for high and low discharge regimes. In this case, the discharge regimes will be separated by a river discharge of 330 m3/s. A different river discharge than the one used in the Markov-chain model approach because in this case all deployments were included in the analysis. Values of from Equation (5) are replaced in Equation (4) to calculate the longitudinal distribution of salinity with the modeled values.

[V1]Inconsistent with the previous 2 sentences….

Show the expression before talking about beta

[V2]What do you mean?

[V3]Numerical?