Modeling ecosystem processes with variable freshwater inflow to the Caloosahatchee River Estuary, southwest Florida. I. Model development

Estuarine, Coastal and Shelf Science 151 (2014) 256e271

Contents lists available at ScienceDirect

Estuarine, Coastal and Shelf Science journal homepage: www.elsevier.com/locate/ecss

Modeling ecosystem processes with variable freshwater inflow to the Caloosahatchee River Estuary, southwest Florida. I. Model development Christopher Buzzelli a, *, Peter H. Doering a, Yongshan Wan a, Detong Sun a, David Fugate b a b

Coastal Ecosystems Section, South Florida Water Management District, 3301 Gun Club Rd., West Palm Beach, FL 33406, USA Department of Marine and Ecological Sciences, Florida Gulf Coast University, 10501 FGCU Blvd. South, Fort Myers, FL 33965, USA

a r t i c l e i n f o

a b s t r a c t

Article history: Received 16 April 2014 Accepted 17 September 2014 Available online 2 October 2014

Variations in freshwater inflow have ecological consequences for estuaries ranging among eutrophication, flushing and transport, and high and low salinity impacts on biota. Predicting the potential effects of the magnitude and composition of inflow on estuaries over a range of spatial and temporal scales requires reliable mathematical models. The goal of this study was to develop and test a model of ecosystem processes with variable freshwater inflow to the sub-tropical Caloosahatchee River Estuary (CRE) in southwest Florida from 2002 to 2009. The modeling framework combined empirically derived inputs of freshwater and materials from the watershed, daily predictions of salinity, a box model for physical transport, and simulation models of biogeochemical and seagrass dynamics. The CRE was split into 3 segments to estimate advective and dispersive transport of water column constituents. Each segment contained a sub-model to simulate changes in the concentrations of organic nitrogen and phosphorus  3 (ON and OP), ammonium (NHþ 4 ), nitrate-nitrite (NOx ), ortho-phosphate (PO4 ), phytoplankton chlorophyll a (CHL), and sediment microalgae (SM). The seaward segment also had sub-models for seagrasses (Halodule wrightii and Thalassia testudinum). The model provided realistic predictions of ON in the upper estuary during wet conditions since organic nitrogen is associated with freshwater inflow and low salinity. Although simulated CHL concentrations were variable, the model proved to be a reliable predictor in time and space. While predicted NO x concentrations were proportional to freshwater inflow, NHþ 4 was less predictable due to the complexity of internal cycling during times of reduced freshwater inflow. Overall, the model provided a representation of seagrass biomass changes despite the absence of epiphytes, nutrient effects, or sophisticated translocation in the formulation. The model is being used to investigate the relative importance of colored dissolved organic matter (CDOM) vs. CHL in submarine light availability throughout the CRE, assess if reductions in nutrient loads are more feasible by controlling freshwater quantity or N and P concentrations, and explore the role of inflow and flushing on the fates of externally and internally derived dissolved and particulate constituents. © 2014 Elsevier Ltd. All rights reserved.

Keywords: estuary freshwater model nutrients light seagrass

1. Introduction Many of the changes in estuarine physical, biogeochemical and biological attributes are consequences of altered patterns of freshwater inputs (Sklar and Browder, 1998; Alber, 2002). Anthropogenic manipulation of freshwater inflow alters salinity distributions that affect estuarine organisms and the structure of food webs (Livingston et al., 1997; Powell et al., 2003; Tolley et al., 2005; Gillson, 2011; Petes et al., 2012). Nutrient loading is a function of

* Corresponding author. E-mail address: [email protected] (C. Buzzelli). http://dx.doi.org/10.1016/j.ecss.2014.08.028 0272-7714/© 2014 Elsevier Ltd. All rights reserved.

freshwater discharge and the incoming concentrations of N and P (Boynton et al., 2008; Greening et al., 2011; Taylor et al., 2011). Symptoms of estuarine eutrophication are often ascribed to the increased loading of dissolved nutrients (e.g. nitrogen and phosphorus or N and P) from the adjacent watershed (Cloern, 2001; Kemp et al., 2005). Although decreasing nutrient loads can have demonstrable benefits, it is possible that an estuary may exhibit oligotrophication or other non-linear trajectories (Nixon et al., 2009; Duarte et al., 2009). While nutrient loads are often directly proportional to discharge, anthropogenic withdrawal of freshwater can confound estuary ecology in light of efforts to manage watershed nutrient inputs (Flemer and Champ, 2006).

C. Buzzelli et al. / Estuarine, Coastal and Shelf Science 151 (2014) 256e271

Variability in discharge on different scales directly affects flushing time, the input of colored dissolved organic matter (CDOM), and estuarine biogeochemistry (McPherson and Miller, 1987; Miller and McPherson, 1991; Bowers and Brett, 2008). Flushing time is one of the key modulators of allochthonous and autochthonous materials (Dettmann, 2001; Sheldon and Alber, 2006; Swaney et al., 2011; Buzzelli et al., 2013d). In fact, phytoplankton production and trophic transfer are directly influenced by freshwater inflow and flushing time (Lucas et al., 2009a; Phlips et al., 2012). In many estuaries, watershed-derived CDOM is an important component of submarine light availability for seagrass habitats (Gallegos and Kenworthy, 1996; Livingston et al., 1998; McPherson et al., 2011; Buzzelli et al., 2012; Le et al., 2013). The relationships among inflow, nutrient loading, salinity gradients, flushing time, phytoplankton, light, and seagrass survival are nonlinear and inter-dependent (Buzzelli, 2011). Salinity serves as an indicator of freshwater inflow that can be used to help define desirable conditions for biota (Chamberlain and Doering, 1998; Doering et al., 2002; SFWMD, 2012). It is possible to develop statistical models to predict salinity with changes in freshwater inflow or inlet dynamics (Morey and Dukhovskoy, 2012; Qiu and Wan, 2013). Salinity time series are essential for salinityproperty and box-model methods (Officer, 1980; Kaul and Froelich, 1984; Hagy et al., 2000; Hagy and Murrell, 2007). Box models have great utility since physical processes can be resolved at appropriate scales for application to ecological models (Kremer et al., 2010). Stand-alone ecological models for phytoplankton, sediment microalgae (SM), oysters, and seagrasses help in the evaluation and definition of environmental requirements such as salinity conditions, light availability, or nutrient sensitivity (Pinckney, 1994; Muylaert et al., 2005; Buzzelli et al., 2012, 2013a). Simulation models for water column and benthic habitats are commonly linked to more sophisticated hydrodynamic models (Cerco and Seitzinger, 1997; Cerco, 2000; Cerco and Moore, 2001; Powell et al., 2003; Cerco and Noel, 2007). However, these computerintensive models are often built upon assumptions and formulations that exceed an empirically based understanding of estuarine biogeochemical processes. Simplified modeling platforms that incorporate desirable features of various approaches are needed to help improve resource management. The quantification and evaluation of multiple stressors on estuarine ecology from days to decades requires an integrative framework linking watershed inputs, circulation, water quality, and biota (Sohma et al., 2008; Condie et al., 2012). Focused models can be used to explore potential ecological feedbacks and responses to variations in environmental drivers including altered inflow, temperature, light and nutrients, and grazing (Scavia and Liu, 2009; Lucas et al., 2009). Successful development of these models are useful to determine what can be predicted, identify missing information, and examine potential estuarine ecological responses to proposed watershed actions (Buzzelli et al., 2013c). The goal of this study was to develop and test a model of ecosystem processes with variable freshwater inflow to the subtropical Caloosahatchee River Estuary (CRE) in southwest Florida from 2002 to 2009. The modeling framework combines empirically derived inputs of freshwater and materials from the watershed, daily predictions of salinity, a box model for physical transport, and simulation models of biogeochemical and seagrass dynamics. The objectives were to 1) create a three segment model to estimate net transport between the watershed and oceanic boundaries, 2) develop biogeochemical simulation models for each segment, 3) calibrate and test model results relative to empirical data from the CRE, and (4) quantify sensitivity of model predictions to variations in key parameters.

257

2. Materials and methods 2.1. Study site The CRE is located in southwest Florida and has been altered by human activities starting in the 1880's when the river was straightened and deepened (Fig. 1; Antonini et al., 2002; Barnes, 2005). Water control structures at Lake Okeechobee (S-77) and Ortona (S-78) were completed in the 1930's with the last installed in 1966 at Olga (S-79; Antonini et al., 2002). The Franklin Lock at S79 represents the head of the CRE that extends ~48 km downstream to the Sanibel Bridge where it empties into the Gulf of Mexico. Early descriptions of the CRE characterize it as barely navigable due to extensive shoals and oyster bars near the estuarine mouth (Sackett, 1888). A navigation channel was dredged with a causeway built across the mouth of San Carlos Bay to Sanibel Island in the 1960's. 2.2. Model overview The CRE was split into 3 segments for development of box and simulation models (Fig. 1). ArcGIS was used to interpolate bathymetric data (NGVD 1988) and to delineate segment boundaries, estimate the area and average elevation of each segment, and calculate the three-dimensional volumes (Buzzelli et al., 2013b). Segment 1 extends approximately 16,098 m from S-79 to where the CRE expands. The segment has an area of 1.5  107 m2, an average sediment elevation (zseg) of 2.3 m, and a volume of 2.1  107 m3. Segment 2 extends approximately 13,405 m downstream from Segment 1 encompassing 3.0  107 m2 of mid-estuary with an average sediment elevation of 2.5 m. The volume of Segment 2 is 7.2  107 m3. Segment 3 extends 18,628 m to the Sanibel Island Bridge and is bounded to the northwest by Pine Island (Fig. 1). The average sediment elevation of Segment 3 is 1.5 m, the area is 2.7  107 m2, and the volume is 5.2  107 m3. The CRE area between S-79 and the Sanibel Bridge is 7.2  107 m2, has a total length of 48,132 m, an average sediment elevation of 2.1 m, and a volume of 1.5  108 m3. The framework includes the box model set-up (Fig. 2) and equations (Table 1), the external boundary inputs and estuarine calibration data (Figs. 3 and 4), an array of forcing functions that drive model processes, a suite of biogeochemical process equations, and ~75 mathematical coefficients (Table 2). The biogeochemical models use an integration interval (dt) of 0.03125 d ¼ 45 min over simulations spanning 2922 d from 2002 to 2009. 2.3. Salinity and transport model The box model approach was used to represent advective and dispersive exchanges between Segments 1e3 (Officer, 1980; Hagy et al., 2000; Hagy and Murrell, 2007). Each segment was assumed to be a homogeneous box of constant volume in the vertical and horizontal dimensions (Fig. 2). The approach conserves volume and mass within and among the three estuarine segments. The box model was driven by daily time series for freshwater inflow at the estuarine head (Q0) and salinity (S) for each segment (S1, S2, S3) and the downstream boundary (S4). Physical transport of a water column constituent (TRSc) was the sum of advection (Xc), lateral inputs from tributaries and ground water (TGc), and dispersion (Yc; Fig. 2; Eqn. (1)). The time series for Q0 (m3 d1) from 2002 to 2009 (2922 days) at S-79 was derived from the hydrologic data-base (DBHydro) at the South Florida Water Management District (Fig. 3A). The loadings of water column constituents at the upstream boundary (Q0c; g d1) were calculated as the product of Q0 and average monthly concentrations (g m3; Fig. 4). A watershed model was used to estimate

258

C. Buzzelli et al. / Estuarine, Coastal and Shelf Science 151 (2014) 256e271

Fig. 1. Map of the Caloosahatchee River Estuary (CRE) in southwest Florida (inset map). The CRE spans from a water control structure (S-79; filled star) to the Sanibel Bridge (open star). The estuary was divided into three segments (yellow, green, blue). Inputs of watershed materials were from monitoring stations (open circles). Concentrations of water column constituents for model calibration were from monitoring stations within the CRE (open squares). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

the daily lateral input from tributaries and ground water to each segment (Qtgw1, Qtgw2, Qtgw3; Y. Wan, unpublished data; Fig. 3B). Total inflow for each segment (Q1, Q2, Q3) was calculated as the sum of the advective component and the lateral inflows (Eqn. 2aec). These total segment inflows were used to predict downstream transport of water column constituents across the segment boundaries (X01, X12; X23, X34; Eqn. 3aed). The loadings of water column constituents from the tributaries and ground water (TG1c, TG2c, TG3c) were the product of the lateral inflows to each segment (Qtgw123) and the corresponding average monthly input concentration (Eqn. 4aec). Monthly constituent concentrations in the tributaries and ground water were derived from Lee County, Florida monitoring stations averaged by segment (Fig. 4). Time series for the average daily salinity of each segment (S1, S2, S3) were derived using a predictive statistical model developed for the CRE (Fig. 3CeD; Qiu and Wan, 2013). The salinity time series were used to calculate the dispersive exchanges (E1, E2, E3; Eqn.5aec) and influence the growth of seagrasses. Dispersion represents the bi-directional, non-advective transport of water column constituents across segment boundaries (Y12 and Y21; Y23 and Y32; Y34 and Y43, Eqn. 6aec; Hagy and Murrell, 2007). 2.4. Biogeochemical models Each of the three segments included a water column sub-model to simulate the concentrations of phytoplankton carbon (Cphyto;

g C m3), organic nitrogen and phosphorus (ON and OP; g N or  3 P m3), ammonium (NHþ 4 or NH; g N m ), nitrate-nitrite (NOx or 3 3 3 NO; g N m ), ortho-phosphate (PO4 or PO; g P m ) and sediment microalgae (SM; g C m2). Downstream Segment 3 also had submodels for the seagrasses Halodule wrightii (Hw; g C m2) and Thalassia testudinum (Tt; g C m2). Biogeochemical processes (gross primary production, nutrient uptake, respiration, remineralization, nitrification, denitrification) were modulated by variations in temperature, depth, and submarine light (Table 1). Water temperature (T) and surface light were modeled using trigonometric functions established for 26 north latitude (Eqn. (7); Buzzelli et al., 2012). Irradiance at the water surface (I0) and photoperiod (Pphoto) were necessary to simulate variations in submarine light (Table 1; Eqns. (8)e(9)). Variable water level (h) was used to calculate depth to affect submarine light penetration. Water level was calculated hourly based on the amplitude, period, and phase of the M2 tide determined for Ft. Myers, FL. Depth (hseg) was calculated as the difference between h and zseg (Table 1; Eqns. (10) and (11)). The total attenuation coefficient for submarine light contained contributions from pure water (kw), color, turbidity (NTU), and chlorophyll a (CHL; Eqn. (12); Christian and Sheng, 2003). Attenuation due to color (kcolor) was estimated using a negative exponential relationship with average salinity of the segment (Sseg, Eqn. (13); McPherson and Miller, 1987; Bowers and Brett, 2008; Buzzelli et al., 2012). Time series for the average turbidity of each segment were derived from monitoring data

0 =S 79

C. Buzzelli et al. / Estuarine, Coastal and Shelf Science 151 (2014) 256e271

Do wn b ou nd 4

1

C

Y1

2

C

Y2

3

C

Y4

3

X

34

Y3

4

Se g

me S3 nt 3

Y3

2

TG

3

X

23

3

Se g

me S2 nt 2

Y2

1

TG

2

X

12

2

Se g

me S1 nt 1

TG

1

X

01

Q

Up b ou nd 0

Fig. 2. Schematic diagram for three-segment model of physical transport in the CRE. The numeric designations are for the upstream boundary (0), the three segments (1, 2, 3), and the downstream boundary (4). Average salinity of the segments (S) and lateral inputs of materials from tributaries and ground water (TG), downstream advective transport across segment boundaries (X), and bi-directional dispersion (Y) are shown. See text for details.

available through DBHydro (http://www.sfwmd.gov/dbhydroplsql/ show_dbkey_info.main_menu). There were specific coefficients for each of the attenuation components (kw, aNTU, aCHL, acolor, bcolor; Table 2). Submarine light (I) was calculated at mid-depth for phytoplankton (Im; Eqn. (14); Robson, 2005), the sediment surface for seagrasses (Iz; Eqn. (15)), and below the seagrass canopy for sediment microalgae (ISM; Eqn. (16); Pinckney, 1994). ISM required the predicted seagrass shoot biomass (Cshoot) and the solar angle and declinations (b and d; Eqns. (17) and (18)). The masses of the water column constituents resulted from the sum of all sources and sinks (Table 1). The mass of each constituent was divided by the segment volume to derive the concentration at each time step. Phytoplankton was a key variable since it receives external inputs of CHL from the watershed, is the primary sink for inorganic nitrogen (N) and phosphorus (P), is the primary source of autochthonous organic N and P, is a central agent in submarine light extinction, and serves as an ecological indicator (Doering et al., 2006; Buzzelli, 2011). In addition to transport, phytoplankton biomass was affected by gross production (Gphyto), respiration (Rphyto), mortality (Mphyto), exudation (Exphyto), sedimentation (Sedphyto), and the resuspension of SM (RSphyto; Eqn. (19)). Gphyto was calculated as a function of the maximum photosynthetic rate (Pm), an exponential increase with T (Tfx), and hyperbolic relationships for light and nutrients dependent upon half-saturation coefficients (Ik, kDIN, kDIP; Eqn. (19a)). Rphyto and Mphyto included constant rate coefficients (kRphyto and

259

kMphyto) modified by Tfx (Eqn. 19bec; Kuo and Park, 1994). Mphyto is assumed to include both grazing and non-grazing death of phytoplankton cells. Exphyto accounted for an estimated 10% of gross production that is lost through exudation (Eqn. (19d); Larsson and Hagstrom, 1979; Bjornsen, 1988). Phytoplankton sinking was a function of a constant rate (ksed ¼ 0.35 m d1; Kuo and Park, 1994; Lucas et al., 2009b) and hseg (Eqn. (19e); Table 2). Conversely, the source of algal biomass from the sediments resulted from the resuspension of a constant fraction of SM (kresus ¼ 0.025) multiplied by Aseg (Eqn. (19f); Table 1). The ON and OP masses in the water column of each segment resulted from transport, production (ONp and OPp), remineralization (ONrem and OPrem), and sinking from the water column to the sediment (ONws and OPws; Eqns. (20) and (21)). The in situ production of ON and OP was calculated as a fraction of the sum of the phytoplankton loss terms (Rphyto, Mphyto, Exphyto, Sedphyto) converted to N and P units using fixed molar C:N (6.625) and C:P (106) ratios, respectively (fONp ¼ fOPp ¼ 0.7; Eqns. (20a) and (21a)). ONrem and OPrem were functions of the ON and OP masses, Tfx, and the remineralization coefficient (kONrem ¼ kOPrem ¼ 0.08 d1; Table 2; Eqns. (20b) and (21b)). The particulate fractions of the ON and OP masses (fPON ¼ fPOP ¼ 0.1) sank at a rate of 0.35 m d1 (kNPsed, Table 2; Eqns. (20d) and (21c)). The masses of NH, NO, and PO in the water column of each segment were influenced by a similar series of source and sink processes (Eqn. (22)e(24), Table 1). ONrem was split in half (fNOrem ¼ 0.5) to account for the production of both NH and NO (NHp and NOp; Eqns. (22a) and (23a)). All of OPrem provided the local source of PO (Eqn. (24a)). Temperature-dependent nitrification (NHnit or NOnit) was a function of the NH mass with a basal rate of 0.01 d1 and provided a sink for NH but a source of internal NO (Eqn. (22b); Table 1; Walsh et al., 2006). Loss of NO through denitrification (NOden) was a function of the NO mass in the segment, the denitrification coefficient (kden ¼ 0.01 d1; Cornwell et al., 1999), and Tfx (Eqn. (23b)). Uptake by phytoplankton represented a significant loss term for water column NH, NO, and PO. Nutrient uptake affects water column inorganic concentrations which influence rates of gross primary production by phytoplankton (Eqn. (19a)). NHup and NOup were calculated using Gphyto converted to N units using the fixed molar CN ratio and scaled by a preference fraction based on NO x (fNOup ¼ 0.45; Table 2) and a feedback function (fbNH and fbNO; Eqns. (22c) and (23c)). POup had a similar feedback function (Eqn. (24b)). The feedback functions for NH, NO x , and PO (Eqns. (22d), (23d) and (24c)) varied from 0 to 1 and were derived via donor control where the mass of the particular constituent within each segment was modulated by limitation (NHlim, NOlim, POlim) and threshold (NHthresh, NOthresh, POthresh) values (Wiegert, 1975; Sin and Wetzel, 2002). The exchanges of NH, NO, and PO between the sediment and water (NHsw, NOsw, POsw) were estimated using empirically derived values averaged by segment and modified by Tfx (Eqns. (22e), (23e) and (24d); Howes et al., 2008; Buzzelli et al., 2013b). The negative values for kNOsw3 and kPOsw3 denoted the removal of water column NO and PO by the sediment boundary in Segment 3. The mass of SM changed with gross production (GSM), phytoplankton sinking (PHYSM), respiration, (RSM), mortality (MSM), and resuspension (RSSM; Table 1). GSM was modeled similarly to phytoplankton with Pmax, Tfx, and biomass (CSM) but did not incorporate nutrient limitation effects (Eqn. (25a)). An assumed 50% of Sedphyto was connected to the SM pool (fphySM ¼ 0.5; Eqn. (25b)). RSM resulted from the product of the respiration coefficient (kRSM ¼ 0.1 d1), Tfx, and CSM (Eqn. (25c)). MSM represented loss to grazing of SM by benthic fauna in the sediment surface layer and was estimated using C2SM and a mortality coefficient (kMSM ¼ 0.5 m2 gC1; Eqn. (25d); Buzzelli, 2008).

260

C. Buzzelli et al. / Estuarine, Coastal and Shelf Science 151 (2014) 256e271

Table 1 Equations for three segment model of the Caloosahatchee River Estuary. See Fig. 2 for diagrammatic representation and Figs. 3e5 for time series of freshwater inflows (Q0, Qtgw1, Qtgw2, Qtgw3), boundary concentrations (C0, Ctgw1, Ctgw2, Ctgw3), and segment salinity (S1, S2, S3, S4). Description

Equations

(1) Summary equation for physical transport (TRSc; g d1) of constituent (C; g m3) in segments (1,2,3) (2abc) Total inflow (Q; m3 d1) to segments (1,2,3) (3abcd) Downstream advective transport (X; g d1) of constituent (C; g m3) between segments (1,2,3) including boundary (4) (4abc) Lateral transport (TG; g d1) of constituent from tributaries and ground water (TGC; g m3) (5abc) Non-advective exchanges (E; m3 d1) between segments (1,2,3) (6abc) Dispersive transport (Y123; g d1) of constituent (C; g m3) between segments (7) Temperature (T;  C)

TRSC ¼ advection þ lateral ± nonadvection ¼ XC þ TGC ± YC

(8) Photoperiod (Pphoto; hrs) (9) Irradiance at water surface (I0; mmole m2 s1) (10) M2 Water Level (h; m) (11) Water Depth of segment (hseg; m) (12) Light extinction coefficient (kt; m1) (13) Light extinction color (kcolor; m1) (14) Light at mid-depth (Im; mmole m2 s1) (15) Light at bottom (Iz; mmole m2 s1) (16) Light at bottom with seagrass (ISM; mmole m2 s1) (17) Solar angle (sinb; l ¼ latitude in radians; t ¼ model hour in radians) (18) Solar declination (d; 4 ¼ day of year in radians) (19) Phytoplankton mass change (Cphyto; gC d1) (19a) Gross primary production (G)

Q1 ¼ Q0 þ Qtgw1; Q2 ¼ Q1 þ Qtgw2; Q3 ¼ Q2 þ Qtgw3 X01 ¼ Q0*C0; X12 ¼ Q1*C1; X23 ¼ Q2*C2; X34 ¼ Q3*C3

TG1c ¼ Qtgw1* Ctgw1; TG2c ¼ Qtgw2* Ctgw2; TG3c ¼ Qtgw3* Ctgw3 *S1 E1 ¼ ðSQ20S 1Þ

0 S1 Q0 S2 Þ E2 ¼ ðE1 S1 þQ ðS2 S3 Þ

0 S2 Q0 S3 E3 ¼ E2 S2 þQ ðS3 S4 Þ

Y12 ¼ E1C1; Y21 ¼ E1C2; Y23 ¼ E2C2; Y32 ¼ E2C3; Y34 ¼ E3C3; Y43 ¼ E3C4
T ¼ 25  5*cos



2*p*day32 365


 Pphoto ¼ 12  2*cos 2*p*day 365 " I0 ¼ MAX

Iamp *cos

!!

2*p*ðhour12Þ 2*Pphoto

#

; 0:0




 h ¼ MSL þ AM2*cos 2*p* hour1:43 TM2 hseg ¼ h  zseg kt ¼ kw þ ½kcolor  þ ½aNTU *NTU þ ½aCHL *CHL kcolor ¼ acolor *eðbcolor *sseg Þ " # Im ¼ I0 

I0 *eðkt *hseg Þ kt *hseg

Iz ¼ I0 *eðkt *hseg Þ




aSH *Cshoot sin b



 ISM ¼ Iz *e d*p  cos l*cos d*p *cos t sin b ¼ sin l*sin 180 180 d ¼ 0:39637  22:9133*cosð4Þ þ 4:02543*sinð4Þ  0:3872*cosð2*4Þ þ 0:052*sinð2*4Þ Dphyto ¼ Gphyto  Rphyto  Mphyto  Exphyto  Sedphyto þ RSSM ±TRSphyto 

 
 Im c c Gphyto ¼ Pm * ðIk þI *ðeKT ðTTopt Þ Þ* ðk DIN * ðk DIP *Cphyto þDIN Þ þDIP Þ mÞ DIN

(19b) Respiration (R)

Rphyto ¼ ½kRphyto *eKT ðTTopt Þ *Cphyto

(19c) Mortality (M)

Mphyto ¼ ½kMphyto *eKT ðTTopt Þ *Cphyto

(19d) Exudation (Ex)

Exphyto ¼ ½GPPphyto *kexu  
 Sedphyto ¼ Cphyto * hksed

(19e) Sedimentation (Sed)

c

DIP

seg

(19f) Resuspension (RS) (19g) Transport (TRS) (20) Organic nitrogen mass change (ON; gN d1) (20a) Production (P) (20b) Reminerlization (Rem) (20c) Temperature effect (Tfx) (20d) Sedimentation (Sed)

RSSM ¼ ½Aseg *CSM *kresus  TRSphyto ¼ ½Xphyto þ TGphyto ±Yphyto  DON ¼ ONp  ONrem  ONws ±TRSON ONp ¼ ½fONp *CN*ðRphyto þ Mphyto þ Exphyto þ Sedphyto Þ ONrem ¼ ½ON*kONrem *Tfx  Tfx ¼ eKT ðTTopt Þ 
 sed ONsed ¼ CON *fPON * kNP h seg

(20e) Transport (TRS) (21) Organic phosphorus mass change (OP; gP d1) (21a) Production (P) (21b) Reminerlization (Rem) (21c) Sedimentation (Sed)

TRSON ¼ ½XON þ TGON ±YON  DOP ¼ OPp  OPrem  OPws ±TRSOP OPp ¼ ½fopp *CP*ðRphyto þ Mphyto þ Exphyto þ Sedphyto Þ OPrem ¼ ½OP*kOPrem *Tfx  
 sed OPsed ¼ COP *fPOP * kNP h seg

(21d) Transport (TRS) 1 (22) Ammonium mass change (NHþ 4 ; gN d ) (22a) Production (22b) Nitrification (22c) Uptake (22d) Donor-Controlled Feedback

TRSOP ¼ ½XOP þ TGOP ±YOP  DNH ¼ NHp  NHnit  NHsw  NHup ±TRSNH NHp ¼ ½ON*kONrem *Tfx *ð1  fNOrem Þ NHnit ¼ ½NH*knit *Tfx  
  Gphyto NHup ¼ *ð1  fNOup Þ*ðfbNH Þ CN


   ðNHNHlim Þ fbNH ¼ max min ðNH NH Þ ; 1:0 ; 0:0 thresh

(22e) Sediment-water (22f) Transport 1 (23) Nitrate-nitrite mass change (NO x ; gN d ) (23a) Production (23b) Denitrification (23c) Uptake (23d) Donor-Controlled Feedback (23e) Sediment-water

lim

NHsw ¼ ½kNHsw *Tfx  TRSNH ¼ ½XNH þ TGNH ±YNH  DNO ¼ NOp þ NOnit  NOden  NOsw  NOup ±TRSNO NOp ¼ ½ON*kONrem *Tfx *fNOrem  NOden ¼ ½NO*kden *Tfx  
  Gphyto NOup ¼ *ðfNOup ÞðfbNO Þ CN


   ðNONOlim Þ min ðNO ; 1:0 ; 0:0 fbNO ¼ max thresh NOlim Þ NOsw ¼ ½kNOsw *Tfx 

c

C. Buzzelli et al. / Estuarine, Coastal and Shelf Science 151 (2014) 256e271

261

Table 1 (continued ) Description

Equations

(23f) Transport 1 (24) Phosphate mass change (PO3 4 ; gP d ) (24a) Production (24b) Uptake

TRSNO ¼ ½XNO þ TGNO ±YNO  DPO ¼ POp  POsw  POup ±TRSPO POp ¼ ½OP*kOP rem *Tfx    G POup ¼ phyto CP *fbPO


   lim Þ min ðPOðPOPO fbPO ¼ max ; 1:0 ; 0:0 thresh POlim Þ

(24c) Donor-Controlled Feedback (24d) Sediment-water (24e) Transport (25) Sediment microalgae (SM; gC m2) (25a) Gross primary production (G) (25b) Phytoplankton input (PHYSM) (25c) Respiration (R) (25d) Mortality (M) (25e) Resuspension (RS) (26) Halodule wrightii (Hw; gC m2) (26a) Gross primary production (GHw)

POsw ¼ ½kPOsw *Tfx  TRSPO ¼ ½XPO þ TGPO ±YPO  DSM ¼ GSM þ PHYSM  RSM  MSM  RSSM 
  ISM *ðeKT ðTTopt Þ Þ *CSM GSM ¼ Pmax * ðIk þI SM Þ 
  Sedphyto PHYSM ¼ *fphySM Aseg RSM ¼ ½kRSM *eKT ðTTopt Þ *CSM 2 MSM ¼ kMSM *CSM RSSM ¼ ½Aseg *CSM *kresus  DHw ¼ GHw  RHw  ðMHw þ MSMw Þ  TRHw 



CHw GHw ¼ PmHw * ðI IzþIz Þ * 1  CHwmax * S

seg þbSHw

RHw ¼ ½kRHw *eKT ðTTopt Þ *CHw 


kMHw *eKT ðTTopt Þ þ MHw ¼ S

 *kSHw *CHw

kHw

(26b) Respiration (R) (26c) Mortality (M) (26d) Translocation (TR) (27) Thalassia testudinum (Tt; gC m2) (27a) Gross primary production (GHw)



 aSHw

seg þbSHw

TRHw ¼ ½ðGHw  RHw Þ*fTRHw  DTt ¼ GTt  RTt  ðMTt þ MSTt Þ  TRTt 2 0


 B 6 6 B Iz CTt * 1  CTtmax *B
GTt ¼ 6PmTt * ðIz þI kTt Þ 4 @

aSHw

1

(27c) Mortality (M)

RTt ¼ ½kRTt *eKT ðTTopt Þ *CTt 20 MTt

1

0

6B B 6B B ¼ 6BkMTt *eKT ðTTopt Þ Þ þ B
4@ @

Tt

3

C 7 7 aSTt C C*kSTt 7*CTt ðSseg cSTt Þ A 5 bS

1þe

(27d) Translocation (TR)

3

C 7 7 aSTt C C*eKT ðTTopt Þ 7*CTt ðSseg cSTt Þ A 5 bS

1þe

(27b) Respiration (R)

 *eKT ðTTopt Þ *CHw

Tt

TRTt ¼ ½ðGTt  RTt Þ*fTRTt 

Fig. 3. (A) Time series of average daily inflow (Q0; m3 d1) from S-79 from 2002 to 2009. (B) Time series of average daily Q0 (shaded) and lateral inflow from tributaries and ground water (QTGW; m3 d1) for segments 1e3. (C) Time series of average daily Q0 (shaded) and salinity for segments 1 and 2 (S1 and S2). (D) Time series of average daily Q0 (shaded) and salinity (S) for segment 3 and the downstream boundary (S3 and S4).

C. Buzzelli et al. / Estuarine, Coastal and Shelf Science 151 (2014) 256e271

NH4+ (mg L-1)

262 0.30 0.25 0.20 0.15 0.10 0.05 0.00

02 NOx- (mg L-1)

0.6 0.5 0.4 0.3 0.2 0.1 0.0

PO4-3 (mg L-1)

02 0.30 0.25 0.20 0.15 0.10 0.05 0.00

02 ON (mg L-1)

1.4 1.2 1.0 0.8 0.6 0.4 0.2 0.0

02 OP (mg L-1)

0.25 0.20

(A) Q0 NH4+

03

04

05

06

07

08

09

10

05

06

07

08

09

10

(B) Q0 NOx-

03

04

(C) Q0 PO4-3

Q0 03

04

05

06

07

08

09

10

(D) Q0 ON

03

04

Qtgw1 Qtgw2 Qtgw3

05

06

07

08

09

10

05

06

07

08

09

10

05

06

07

08

09

10

(E) Q0 OP

0.15 0.10 0.05 0.00

02 CHL (ug L-1)

40 30

03

04

(F) Q0 Chl

20 10 0

02

03

04

 3 Fig. 4. Time series of watershed concentrations of (A) NHþ 4 , (B) NOx , (C) PO4 , (D) ON, (F) OP, and (G) CHL used as inputs to the CRE model from 2002 to 2009. Each constituent contributes loading at the estuarine head (Q0) and combined tributary and ground-water inflows for segments 1e3 (Qtgw1, Qtgw2, Qtgw3).

Changes in the above-ground carbon biomasses of Hw and Tt (CHw and CTt) resulted from gross production (GHW and GTt), respiration (RHw and RTt), mortality through leaf senescence (MHw and MTt), mortality due to salinity fluctuations (MSHw and MSTt), and translocation to below-ground components (TRHw and TRTt; Eqns. (26) and (27)). GHW and GTt each included terms for the maximum rate of photosynthesis (PmHw and PmTt), light limitation (IkHw and IkTt), a biomass maximum feedback function (CHwmax and CTtmax), salinity-response functions, Tfx, and the respective shoot biomasses (Eqns. (26a) and (27a); Herzka and Dunton, 1997; Burd and Dunton, 2001; Eldridge et al., 2004). The salinity-response functions for both Hw and Tt in the CRE were derived experimentally (Doering and Chamberlain, 2000; Doering et al., 2002; Buzzelli et al., 2012). Daily values for the specific rate of Hw blade turnover, and, the length of Tt blades over a range of experimental salinities (3.5e35) were scaled from 0 to 1 and fit to hyperbolic and sigmoidal curves, respectively (aSHw bSHw; aSTt, bSTt, cSTt; Table 2). RHw, RTt, MHw, and MTt were formulated similarly as respiration for phytoplankton (Eqns. (26b) and (27b); Eqns. (26c) and (27c)). The calculation of MSHw and MSTt included

the salinity-response function and a coefficient for salinity-based mortality (kSHw ¼ kSTt ¼ 0.001 d1; Table 2; Eqns. (26c) and (27c)). Constant fractions of net production (G-R) were translocated out of the shoots to the root-rhizomes (fTRHw ¼ 0.3; fTRTt ¼ 0.2; Table 2; Eqns. (26d) and (27d)). 2.5. Calibration and sensitivity testing Calibration involved iterative adjustments of coefficients to improve the correspondence between model predictions and  3 observed concentrations of CHL, NHþ 4 , NOx , PO4 , ON, and OP among the three segments. Data were available for a series of stations within the CRE through DBHydro of the SFWMD (CES01-09). Monthly average concentrations for each variable were aggregated by segment from the field data available from 2002 to 2009. The averages and standard deviations for each variable and segment were calculated for the dry (NovembereApril) and wet (MayeOctober) seasons for both the model predictions and the field data. The agreement between the seasonal-segmented averages of the model and field concentrations were compared using simple

C. Buzzelli et al. / Estuarine, Coastal and Shelf Science 151 (2014) 256e271

263

Table 2 List of model parameters. Parameter

Value

Unit

Description

Source

Iamp MSL TM2 AM2 kw aNTU aCHL aSH acolor bcolor Pm Ik Topt KtB kDIN kDIP kRphyto kMphyto ksed kexu CN CP CCHL kresus fONp fPON kONrem kOPrem fOPp fPOP kNPsed fNOrem kden fNOup knit kRSM kMSM fphySM NHlim NHthresh NOlim NOthresh POlim POthresh kNHsw1 kNHsw2 kNHsw3 kNOsw1 kNOsw2 kNOsw3 kPOsw1 kPOsw2 kPOsw3 PmHw IkHw kRHw kMHw kSHw aSHw bSHw CmaxHw fTRHw PmTt IkTt kRTt kMTt kSTt aSTt bSTt cSTt CmaxTt fTRTt Aseagrass

1000 0.0 12.42 0.111 0.15 0.062 0.058 0.002 2.89 0.096 2.5 100 25 0.069 0.05 0.005 0.15 0.15 0.35 0.1 6.625 106 50 0.025 0.7 0.1 0.08 0.08 0.7 0.1 0.35 0.5 0.01 0.45 0.01 0.1 0.5 0.5 200,000 400,000 200,000 400,000 500,000 1,500,000 277,572 273,507 113,524 277,583 67,273 15,273 74,299 12,270 33,803 0.2 319 0.01 0.005 0.001 1.257 7.256 10 0.3 0.1 150 0.005 0.0025 0.001 1.01 4.55 16.8 600 0.2 3.36Eþ07

mmole m2 d1

Amplitude of surface irradiance Mean sea level Period of M2 tide Amplitude of M2 tide Attenuation due to water Attenuation factor for turbidity Attenuation factor for chlorophyll a Attenuation due to canopy biomass Constant for salinity-color relationship Constant for salinity-color relationship Algal max photosynthetic rate Algal half-saturation for light Optimum temperature for rate processes Rate constant for temperature effect Algal half-saturation for nitrogen Algal half-saturation for phosphorus Algal basal respiration rate Algal basal mortality rate Algal sedimentation rate Fraction of gross production to exudation Redfield C:N molar ratio Redfield C:P molar ratio Carbon:chlorophyll a ratio Fraction of sediment microalgae resuspended Fraction of phytoplankton loss to ON production Fraction of ON that is particulate ON remineralization constant OP remineralization constant Fraction of phytoplankton loss to OP production Fraction of OP that is particulate Particulate sedimentation rate Fraction of ON remineralization to NO-x Basal rate of denitrification Fraction of phytoplankton N uptake - NO-x Basal rate of nitrification Sediment microalgae basal respiration rate Sediment microalgae basal mortality rate Fraction of phytoplankton sinking to SM pool Donor-controlled limitation value for NH mass Donor-controlled threshold value for NH mass Donor-controlled limitation value for NO mass Donor-controlled threshold value for NO mass Donor-controlled limitation value for PO mass Donor-controlled threshold value for PO mass NH exchange sediment to water 1 NH exchange sediment to water 2 NH exchange sediment to water 3 NO exchange sediment to water 1 NO exchange sediment to water 2 NO exchange sediment to water 3 PO exchange sediment to water 1 PO exchange sediment to water 2 PO exchange sediment to water 3 Halodule wrightii max photosynthetic rate Halodule wrightii light constant Halodule wrightii shoot respiration rate Halodule wrightii shoot mortality rate Halodule wrightii shoot mortality rate-salinity Halodule wrightii shoot salinity constant a Halodule wrightii shoot salinity constant b Halodule wrightii shoot max biomass Halodule wrightii fraction translocated to RR Thalassia testudinum max photosynthetic rate Thalassia testudinum light constant Thalassia testudinum shoot respiration rate Thalassia testudinum shoot mortality rate Thalassia testudinum shoot mortality rate-salinity Thalassia testudinum shoot salinity constant a Thalassia testudinum shoot salinity constant b Thalassia testudinum shoot salinity constant c Thalassia testudinum shoot max biomass Thalassia testudinum fraction translocated to RR Area of SAV in segment 3

Local data *** NOAA Ft. Myers NOAA Ft. Myers Calculated from Gallegos 2001 McPherson and Miller 1987 McPherson and Miller 1987 Buzzelli et al., 1999 McPherson and Miller 1987 McPherson and Miller 1987 Muylaert et al., 2005 Grangere et al., 2009 Buzzelli et al., 1999 Buzzelli et al., 1999 James et al., 2005 James et al., 2005 Lucas et al., 2009 Kuo and Park 1994 Lucas et al., 2009 Larsson and Hagstrom 1979 Walsh et al., 2006 Walsh et al., 2006 Walsh et al., 2006 Buzzelli et al., 1999 Calibration Calibration Calibration Calibration Calibration Calibration Kuo and Park 1994 Calibration Howes et al., 2008 Calibration Walsh et al., 2006 Buzzelli et al., 1999 Buzzelli et al., 1999 Calibration Calibration-CRE data Calibration-CRE data Calibration-CRE data Calibration-CRE data Calibration-CRE data Calibration-CRE data Buzzelli et al., 2013a Buzzelli et al., 2013a Buzzelli et al., 2013a Buzzelli et al., 2013a Buzzelli et al., 2013a Buzzelli et al., 2013a Buzzelli et al., 2013a Buzzelli et al., 2013a Buzzelli et al., 2013a Burd and Dunton 2001 Burd and Dunton 2001 Burd and Dunton 2001 Burd and Dunton 2001 Calibration Calibration - Doering et al., 2002 Calibration - Doering et al., 2002 Calibration - CRE data Calibration Eldridge et al., 2004 Herzka and Dunton 1997 Eldridge et al., 2004 Eldridge et al., 2004 Calibration Calibration e Doering et al., 2002 Calibration e Doering et al., 2002 Calibration e Doering et al., 2002 Calibration e CRE data Calibration CRE data

m hours m m1 m1 NTU1 m3 mg1 m1 m2 gC1 m1 Unitless d1 mmole m2 d1  C  1 C g N m3 gP m3 d1 d1 m d1 Unitless fraction mole C moleN1 mole C moleP1 gC gChla1 Unitless fraction Unitless fraction Unitless fraction d1 d1 Unitless fraction Unitless fraction m d1 Unitless fraction d1 Unitless fraction d1 d1 m2 gC1 Unitless fraction gN gN gN gN gP gP gN d1 gN d1 gN d1 gN d1 gN d1 gN d1 gP d1 gP d1 gP d1 d1 mmole m2 d1 d1 d1 d1 Unitless Unitless g C m2 Unitless fraction d1 mmole m2 d1 d1 d1 d1 Unitless Unitless Unitless gC m2 Unitless fraction m2

C. Buzzelli et al. / Estuarine, Coastal and Shelf Science 151 (2014) 256e271

coefficients (kRphyto, kden, kONrem, kDIN) were selected because of their importance in water column biogeochemical cycling. In separate simulations each of these coefficients was varied by 0% (base), 10%, 20%, 50%, þ10%, þ20%, and þ50%. The effects of the altered coefficients on concentrations of CHL, dissolved inorganic  nitrogen (DIN ¼ NHþ 4 þ NOx ), and the total light extinction coefficient (kt) throughout the entire CRE were quantified using the average percent difference between the base and specific test simulation (%diff ¼ (test  base)/base*100). Segment-specific model concentrations of CHL and DIN were multiplied by the respective segment volumes (V1, V2, V3), summed, and divided by the total volume of the CRE to derive volume-weighted concentrations for the entire CRE (CHLCRE and DINCRE). Daily values for kt were averaged across all three segments (ktCRE). There were a total of 84 simulations to test for model sensitivity (4 coefficients  7 values  3 response variables).

correlation (r; James et al., 2005). Calibration was assumed to be complete when correlation coefficients among all water constituents in all three segments were maximized. There were limited data available to calibrate submarine light dynamics, biogeochemical rate processes, and benthic biomass estimates. Water column gross production was measured at four locations along the CRE in February, April, May, June, and July of 2009 (Phlips and Mathews, 2009). The assessment of sedimentwater exchanges of dissolved N and P from February 2008 included the determination of denitrification and sediment microalgal biomass (Howes et al., 2008; Buzzelli et al., 2013b). Data from these studies were aggregated and averaged across the entire CRE for qualitative comparison to model results. Total light extinction in the water column of the CRE was monitored at multiple stations sporadically from 2005 to 2009 (CES01-CES09; Fig. 1). The average monthly kt predicted by the model was analyzed graphically relative to the field values for each segment for this period of record. Seagrass biomass was monitored approximately monthly at four sites in the lower CRE from 2004 to 2007 (Mazzotti et al., 2007). Data on Hw and Tt shoot biomasses (gdw m2) were converted to C units (g C m2) to derive a time series for each species that were compared to model predictions (Duarte, 1990). The model framework included 23 state variables, ~75 coefficients, and hundreds of physical and biogeochemical processes. Based on previous experience with estuarine water quality models, preliminary simulations, and many calibration tests, four

2

04

05

06

07

08

09

1.5

kntu

15

0 02

10

1.0 0.5 0.0 03

04

05

06

07

08

09

10

35

3.0

(E) Segment 2

30

model data

5

2.5

25

4

salinity

-1

kcolor

5

(B) Segment 2

ktotal (m )

20

10

model data

03

2.0

salinity kchl

3 2

2.0

20

1.5

15

1.0

10

1

0.5

5 03

04

05

06

07

08

09

0 02

10

0.0 03

04

05

06

07

08

09

10 3.0

35

6

(C) Segment 3

salinity kchl

30

model data

salinity

3 2

(F) Segment 3

2.5

kcolor kntu

25

4

2.0

20

1.5

15

1.0

10

1

5

0 02

0 02

03

04

05

06

07

08

09

10

m-1

salinity

-1

ktotal (m )

3

6

-1

2.5

25

4

1

ktotal (m )

3.0

(D) Segment 1

30

5

5

Freshwater inflow to the CRE reflected the dry-wet seasonality of south Florida (Buzzelli et al., 2013c, 2013d). Releases of freshwater through S-79 in both dry and wet seasons followed the operational schedule to maintain pool elevations in the C-43 canal and Lake Okeechobee (SFWMD, 2010; Buzzelli et al., 2014a). Both the sub-tropical seasonality and water management practices that

m-1

(A) Segment 1

0 02

3.1. Freshwater inflow, salinity, and light

35

6

0 02

3. Results

m-1

264

0.5 0.0 03

04

05

06

07

08

09

10

Fig. 5. Time series of light extinction coefficient (kt m1) from 2002 to 2009. Comparison between data (open circles) and model (lines) for Segment 1 (A), Segment 2 (B), and Segment 3 (C). Time series of average segment salinity (S; gray shade) and the components of kt. Included are light extinction due to chlorophyll a (kchl; green), color (kcolor; yellow), and turbidity (kntu; black dash) and S for Segment 1 (D), Segment 2 (E), and Segment 3 (F). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

C. Buzzelli et al. / Estuarine, Coastal and Shelf Science 151 (2014) 256e271 1.5

0.15

50

(D) Segment 1 OP 40

0.9

0.09

30

0.6

0.06

20

0.3

0.03

10

04

05

06

07

08

09

10

0.00 02

03

04

05

07

08

09

0 02

10

0.15

40

0.9

0.09

30

0.6

0.06

20

0.3

0.03

10

04

05

06

07

08

09

10

0.00 02

03

04

05

06

07

08

09

0 02

10

0.15

1.5

(C) Segment 3 ON 0.9

0.09

15

0.6

0.06

10

0.3

0.03

5

05

06

07

08

09

10

0.00 02

07

08

09

10

04

05

06

07

08

09

10

09

10

(I) Segment 3 CHL 20

04

03

(F) Segment 3 OP 0.12

03

06

25

1.2

0.0 02

05

(H) Segment 2 CHL

0.12

03

04

(E) Segment 2 OP

1.2

0.0 02

03

50

(B) Segment 2 ON

g m-3

06

mg m-3

03

mg m-3

0.12

1.5

g m-3

(G) Segment 1 CHL

1.2

03

04

05

06

07

08

09

0 02

10

mg m-3

g m-3

(A) Segment 1 ON

0.0 02

265

03

04

05

06

07

08

Fig. 6. Time series of monthly concentrations in the water column from CRE monitoring data (points) and the simulation model (line). (AeC) Organic nitrogen (ON; g m3) in Segments 1e3. (DeF) Organic phosphorus (OP; g m3). (GeI) Chlorophyll a (CHL; mg m3).

0.25

0.7

(A) Segment 1 NH4+

0.30

(D) Segment 1 NOx-

0.6

0.20

g m-3

0.5 0.15

0.4

0.10

0.3

0.20 0.15 0.10

0.2 0.05 0.00 02

0.05

0.1

03

04

05

06

07

08

09

10

0.0 02

03

04

05

06

07

08

09

10

(B) Segment 2 NH4+

0.20

(E) Segment 2 NOx-

0.6

0.4

0.10

0.3

03

04

05

06

07

08

09

10

(H) Segment 2 PO4-3

0.25

0.5 0.15

0.00 02 0.30

0.7

0.25

g m-3

(G) Segment 1 PO4-3

0.25

0.20 0.15 0.10

0.2 0.05 0.00 02

0.05

0.1

03

04

05

06

07

08

09

10

0.25

0.0 02

03

04

05

06

07

08

09

10

0.35

(C) Segment 2 NH4+

0.20

0.00 02

03

04

05

06

07

08

09

10

0.15

(F) Segment 3 NOx-

0.30

(I) Segment 3 PO4-3

0.12

g m-3

0.25 0.15

0.20

0.09

0.10

0.15

0.06

0.10 0.05 0.00 02

0.03

0.05

03

04

05

06

07

08

09

10

0.00 02

03

04

05

06

07

08

09

10

0.00 02

03

04

05

06

07

Fig. 7. Time series of monthly concentrations in the water column from CRE monitoring data (points) and the simulation model (line). (AeC) Ammonium 3 3 3 1e3. (DeF) Nitrate þ nitrite (NO x ; g m ). (GeI) Ortho-phosphate (PO4 ; g m ).

08

(NHþ 4;

09

10 3

gm

) in Segments

0.28 0.08 0.90 0.38 0.79 0.03 0.07 0.42 0.45 0.02 0.65 0.37 0.88 0.98 0.65 0.94 0.51 0.69 0.63 0.75 0.76 0.07 0.08 0.51 0.02 0.02 0.01 0.02 ± ± ± ± 0.06 0.05 0.04 0.05 0.24 0.18 0.11 0.19 ± ± ± ± 0.57 0.49 0.40 0.48 0.01 0.01 0.01 0.02 ± ± ± ± 0.07 0.05 0.03 0.05 0.05 0.05 0.02 0.06 ± ± ± ± 0.15 0.15 0.07 0.12 0.01 0.01 0.01 0.02 ± ± ± ± 0.05 0.03 0.02 0.03 4.0 3.8 2.6 4.1 ± ± ± ± 14.2 11.6 8.5 11.4 0.01 0.02 0.01 0.02 ± ± ± ± 0.06 0.06 0.03 0.05 0.06 0.12 0.10 0.03 ± ± ± ± 1.20 0.98 0.52 0.90 0.04 0.03 0.01 0.05 ± ± ± ± 0.12 0.10 0.03 0.09 0.06 0.07 0.03 0.09 ± ± ± ± 0.20 0.16 0.04 0.13 0.02 0.01 0.01 0.02 ± ± ± ± 3.9 5.2 1.1 4.6 ± ± ± ±

0.06 0.05 0.03 0.05

0.21 0.02 0.02 0.44 0.14 0.34 0.43 0.51 0.67 0.76 0.02 0.75 0.91 0.84 0.54 0.68 0.003 0.26 0.05 0.04 0.04 0.03 0.03 0.03 0.43 0.31 0.25 0.33 6.8 5.1 3.6 5.2 0.05 0.05 0.04 0.05 0.99 0.81 0.59 0.79 4.5 2.5 1.5 3.4 ± ± ± ±

Dry Seg 1 7.4 2 7.1 3 3.9 CRE 6.2 Wet Seg 1 9.4 2 8.8 3 3.2 CRE 7.1

CHL

0.05 0.04 0.04 0.04

± ± ± ±

0.02 0.02 0.01 0.02

0.28 0.11 0.07 0.15

± ± ± ±

0.07 0.09 0.03 0.11

0.07 0.05 0.04 0.06

± ± ± ±

0.02 0.01 0.02 0.02

ON

± ± ± ±

0.01 0.10 0.07 0.19

OP

± ± ± ±

0.01 0.02 0.02 0.01

CHL

± ± ± ±

2.3 1.5 0.7 2.1

0.14 0.13 0.09 0.12

± ± ± ±

0.04 0.04 0.03 0.04

0.29 0.29 0.15 0.25

± ± ± ±

0.06 0.06 0.05 0.09

0.07 0.05 0.03 0.05

± ± ± ±

0.01 0.01 0.01 0.02

ON

± ± ± ±

0.13 0.08 0.06 0.12

OP

± ± ± ±

0.02 0.01 0.01 0.01

0.61 0.50 0.09 0.63

ON PO3 4 NO-x NHþ 4 CHL

r

PO3 4 NO x NHþ 4 Model

PO3 4 NO x

The agreement between field observations and model predictions varied by constituent (CHL, ON, OP, NH, NO, PO), segment (1, 2, 3), year (2002e2009), season (dry vs. wet), and month (Figs. 6 and 7). The model provided better predictions of ON concentrations in Segments 1 and 2 during the 2002e2006 wet period compared to drier conditions from 2007 to 2009 (Fig. 6AeB). Agreement was reversed for ON in Segment 3 as model concentrations better approximated monthly average values during the last three years of simulation (Fig. 6C). The model generally underestimated average ON concentrations within the CRE (0.33 and 0.48 g m3 vs. 0.79 and 0.90 g m3) with correlation coefficients between the data and model of 0.51 (dry) and 0.42 (wet; Table 3). Average Model OP was similar in magnitude to field concentrations in both the dry and wet seasons (0.03 and 0.05 g m3 vs. 0.05 and 0.05 g m3). The model overestimated OP concentrations in Segment 2 from 2002 to 2007 (Table 3; Fig. 6E). CHL was the most variable of the constituents in the calibration data (Table 3). The variance in average monthly CHL observed in the CRE increased with the onset of drier conditions in 2007. CHL concentrations from the model were within range of observations in Segments 1 and 2 from 2002 to 2006 (Fig. 6GeH). The correlations between the CHL model and data were 0.63 and 0.51 in the dry and wet seasons, respectively (Table 3). While the model best approximated data observations in Segment 1 (rdry ¼ 0.61 and rwet ¼ 0.76), the model over-estimated CHL in the wet season across all three segments (14.2 vs. 9.4 mg m3, 11.6 vs. 8.8, 8.5 vs. 3.2 mg m3). Average monthly concentrations of NH in the model and field observations varied greatly from 0.01 to 0.24 g m3 in Segments 1 and 2 (Fig. 7AeB). NH concentrations were reduced and less variable in Segment 3 (Fig. 7C). In situ NH concentrations were slightly greater in wet (0.05 g m3) vs. dry (0.04 g m3) seasons (Table 3). The model overestimated NH concentrations for all three segments in the dry season (0.14 vs. 0.05 g m3, 0.13 vs. 0.04 g m3, 0.09 vs. 0.04 g m3) with comparatively low correlation coefficients of 0.003, 0.26, and 0.05 (Table 3). The reduced NH concentrations

NHþ 4

3.2. Model calibration

Data

support regional flood protection are influenced by climatic variations in rainfall (SFWMD, 2012). Climatic variability accounted for comparatively wet conditions from 2002 to 2006 followed by drier conditions and drought across south Florida from 2007 to 2009. Q0 ranged from 0.0 to 5.4  107 m3 d1 during 2002e2009 (Fig. 3). Inflow was generally greater from 2002 to 2006 before the onset of drought. Lateral inflows from tributaries and ground water (Qtgw123) followed patterns of Q0 with values approximately an order of magnitude less than discharge at S-79 (Fig. 3B). The largest segment (Segment 2) had the greatest amount of lateral freshwater inflow. S ranged from 0 to 20 in Segment 1, 0 to 30 in Segment 2, <5 to 35 in Segment 3, and <5 to 40 in San Carlos Bay (Fig. 3CeD). The simulated light extinction coefficient (kt) in Segments 1e3 fluctuated seasonally and inter-annually from 2002 to 2009 (Fig. 5AeC). Both total submarine light and the correspondence of average kt between the model and data was improved in the downstream Segment 3 relative to Segments 1 and 2 from 2005 to 2009. Reduced availability of submarine light (e.g. higher kt) was associated with increased Q0 (Fig. 5DeF). Attenuation due to color ranged from ~0.5 to >2.5 m1 and was the dominant component of kt in the upper estuary. Increased values of kchl and kNTU ranged from 0.1 to 1.0 m1 with greater values in the wet seasons. All three components of light attenuation were <0.5 m1 as S approached 30 in Segment 3 during the drought conditions from 2007 to 2009 (Fig. 5F).

OP

C. Buzzelli et al. / Estuarine, Coastal and Shelf Science 151 (2014) 256e271 Table 3 3 3 3  3 Statistical summary of model calibration from 2002 to 2009. Concentrations of water column chlorophyll a (Chl-a; mg m3), ammonium (NHþ 4 ; g m ), nitrate þ nitrite (NOx ; g m ), ortho-phospate (PO4 ; g m ), organic nitrogen (ON; g m3), and organic phosphorus (OP; g m3) from the field data and model were analyzed by segment (seg ¼ 1, 2, 3) and season (dry vs. wet). Dry season results are provided in the upper panel and wet season results are provided in the lower panel. Given are the seasonal average and standard deviation for the data (left), model (shaded; center), and the correlation coefficient (r) between data and model seasonal averages (right). The last row for each season provides the average and standard deviations across the entire estuary.

266

C. Buzzelli et al. / Estuarine, Coastal and Shelf Science 151 (2014) 256e271

predicted for the wet season agreed better with field observations (r ¼ 0.75). Model representation of NO concentration was very good in Segments 1 and 2 from 2002 to 2009 as NO in the upper CRE tracked intra- and inter-annual changes in Q0 (Fig. 7DeE). In the dry and wet seasons, correlation coefficients with the monitoring data were 0.91 and 0.88 (Segment 1) and 0.84 and 0.98 (Segment 2; Table 3). The influence of Q0 diminished in Segment 3 where concentrations were approximately 50% of those upstream and the correlations were not as strong (rdry ¼ 0.54 and rwet ¼ 0.65, respectively). While the average concentration of NO in the CRE predicted by the model (0.12 g m3) was similar to the field average (0.13 g m3; r ¼ 0.94) in the wet season, the model overestimated NO in the dry season (0.25 vs. 0.15 g m3; r ¼ 0.68; Table 3). Maximum PO concentrations approached 0.25 g m3 in Segments 1 and 2 but only 0.15 g m3 in Segment 3 (Fig. 7GeI). PO observations exhibited a climatic shift from reduced concentrations during the wet period (2002e2006) to increased values from 2007 to 2009. Model predictions of water column PO concentrations in Segments 1 and 2 approximated those observed from 2002 to 2006 (Fig. 7GeH). By contrast, agreement between the model and field observations of PO in Segment 3 was minimal from 2004 to 2005 and 2008e2009 (Fig. 7I). While the Segment 3 correlation was only 0.02 in dry season, the value increased to 0.65 in the wet season (Table 3). Overall, the model represented average PO concentrations in the CRE in the dry season (0.05 vs. 0.06 g m3; r ¼ 0.75) although wet season values were underestimated (0.05 vs. 0.09 g m3; r ¼ 0.37). Model SM biomass (mg CHL m2) followed intra- and interannual patterns of submarine light in all three segments (Fig. 8A). Given a more favorable light environment in Segment 3, SM concentrations were 2-3X greater than in Segments 1 and 2 despite the presence of the seagrass canopy. Average monthly concentrations in Segments 1 and 2 ranged from 2 to 10 mg CHL m2 from 2002 to 2006. While SM biomass remained low in Segment 1, values increased to >20 mg CHL m2 in Segment 2 from 2007 to 2008 before declining to <10 mg CHL m2 in 2009. Sediment chlorophyll averaged 21.7 ± 26.2 mg m2 throughout the CRE in February 2008 (Howes et al., 2008; Buzzelli et al., 2013b). The average and standard deviation from all three segments over the 2922 d simulation was 10.1 ± 7.3 mg CHL m2. The responses of Hw and Tt to changes in salinity and submarine light were evident in both the model and field observations (Fig. 8BeC). Shoot biomass declined for both seagrasses with decreased salinity (and associated decreased light) in 2004e2005. By contrast, increased salinity in 2007 contributed to seagrass recovery. Variability in shoot biomass was great at the field sites. While the maximum shoot biomass of Hw approached 5 g C m2 in 2007, the biomass of Tt was >400 g C m2. The model reflected the observed increase (2004) and subsequent decline (2005) in Hw shoot biomass (Fig. 8B). However, Hw above-ground biomass recovered faster in 2006 than predicted. In contrast, the model underestimated the shoot biomass of Tt in 2004e2005 (Fig. 8C). Model Tt biomass corresponded to the observed decrease (2006) and subsequent increase (2007). The model overestimated the shoot biomass of Tt in 2007. 3.3. Sensitivity testing The basal rate of phytoplankton respiration (kRphyto ¼ 0.15 d1) was varied from 0.135 d1 (10%) to 0.225 d1 (þ50%) in sensitivity tests (Table 4). Decreasing this value by 10%, 20%, and 50% led to relative increases in CHLCRE of 6.9%, 14.46%, and 42.12%, respectively. The effects on the concentration of DIN in the CRE were compelling as the 50% value for kRphyto (0.075 d1) led to 23.73%

267

Fig. 8. (A) Time series of predicted average daily biomass of sediment microalgae (SM; g C m2) for Segment 1 (black dash), Segment 2 (green), and Segment 3 (blue) from 2002 to 2009. (B) Time series of Halodule wrightii average monthly shoot biomass data (filled circles; g C m2) and predicted daily biomass (black line). The average daily salinity of Segment 3 is shown (gray area). (C) Time series of Thalassia testudinum average monthly shoot biomass data (filled circles; g C m2) and predicted daily biomass (black line). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

less DIN. Decreasing this coefficient resulted in linear decreases in light penetration (kt ¼ 1.87, 3.92, and 11.49 m1) while increased values decreased light extinction (kt ¼ 1.71, 3.27, and 7.32 m1). The denitrification rate (kden ¼ 0.1 d1) had comparatively less influence on CHLCRE, DINCRE, and ktCRE than the other coefficients tested (Table 4). Decreasing this coefficient by 10%, 20%, and 50% led to increased concentrations of CHLCRE (0.38%, 0.78%, 2.03%) and DINCRE (0.47%, 0.96%, 2.54%). Increasing kden by 10%, 20%, and 50% led to decreased concentrations of CHLCRE (0.37%, 0.73%, 1.75%) and DINCRE (0.46%, 0.90%, 2.15%). There were implications for submarine

268

C. Buzzelli et al. / Estuarine, Coastal and Shelf Science 151 (2014) 256e271

Table 4 Results from sensitivity tests of Caloosahatchee River Estuary (CRE) simulation model. Through a series of model runs four individual parameters (kRphyto, kden, kONrem, kDIN) were varied by 50%, 20%, 10%, þ10%, þ20%, and þ50% to test effects on the average chlorophyll a in the CRE (CHLCRE), the average concentration of DIN (DINCRE), and the average total light attenuation coefficient (ktCRE). Results indicate the percent difference between the base and test values (% difference ¼ (Test  Base)/Base*100) averaged over 2922 days of simulation from 2002 to 2009. Model test

Value

CHLCRE

DINCRE

ktCRE

kRphyto Base 10% 20% 50% þ10% þ20% þ50%

0.15 0.135 0.120 0.075 0.165 0.180 0.225

6.90 14.46 42.12 6.30 12.07 26.82

5.46 10.55 23.73 5.85 12.10 33.13

1.87 3.92 11.49 1.71 3.27 7.32

kden Base 10% 20% 50% þ10% þ20% þ50%

0.1 0.009 0.008 0.005 0.011 0.012 0.015

0.38 0.78 2.03 0.37 0.73 1.75

0.47 0.96 2.54 0.46 0.90 2.15

0.13 0.27 0.71 0.13 0.25 0.59

kDIN Base 10% 20% 50% þ10% þ20% þ50%

0.05 0.045 0.040 0.025 0.055 0.060 0.075

2.62 5.77 21.34 2.26 4.28 9.52

4.87 9.75 23.99 4.84 9.61 23.44

0.75 1.69 6.93 0.61 1.15 2.51

kONrem Base 10% 20% 50% þ10% þ20% þ50%

0.08 0.072 0.064 0.040 0.088 0.096 0.120

3.17 6.57 18.46 2.97 5.75 13.14

3.37 6.95 19.44 3.19 6.22 14.45

1.09 2.26 6.28 1.03 1.99 4.56

light as ktCRE increased with reduced kden (0.13%, 0.27%, 0.71%) and decreased with increased kden (0.13%, 0.25%, 0.59%). Decreasing the nitrogen half-saturation coefficient (kDIN ¼ 0.05 g m3) by 10% and 20% increased CHLCRE (2.62% and 5.77%) and decreased DINCRE (4.87% and 9.75%; Table 4). The 50% decrease in kDIN promoted a 21.34% increase for CHLCRE, but led to a 23.99% decrease in DINCRE. Light extinction responded as ktCRE increased by 0.75%,1.69%, and 6.93% with lower levels of kDIN. The 50% increase in kDIN increased DINCRE by 23.44% but decreased CHLCRE by only 9.52% (Table 4). Submarine light was degraded by 2.51%. CHLCRE and DINCRE concentrations exhibited a similar series of decreases (3.17% and 3.37%; 6.57% and 6.95%; 18.46% and 19.44%) with decreased values (10%, 20%, 50%) for the ON remineralization rate (kONrem ¼ 0.08 d1). ktCRE decreased by 1.09%, 2.26%, and 6.28% as kONrem was reduced from 0.08 g m3 to 0.04 d1 (Table 4). Increasing kONrem by 10% and 20% triggered proportional increases in CHLCRE (2.97% and 5.75%), DINCRE (3.19% and 6.22%), and ktCRE (1.03% and 1.99%). The three response variables in the CRE increased by 13.14%, 14.45%, and 4.56% when kONrem ¼ 0.120d1. 4. Discussion Variations in freshwater inflow have ecological consequences for estuaries ranging among nutrient loading and eutrophication, flushing and transport, and high and low salinity impacts on biota

(Dettmann, 2001; Lucas et al., 2009a; Alber, 2002; Gillson, 2011; Gonzalez and Drake, 2012). The potential effects of altered discharge are heightened in sub-tropical estuaries with heavily modified watersheds and managed freshwater supplies (Childers et al., 2006; Petes et al., 2012; Buzzelli et al., 2013d). The Caloosahatchee River Estuary (CRE) in southwest Florida has a watershed characterized by extensive agriculture and urbanization, is influenced by regulated freshwater outflow from Lake Okeechobee, and contains valuable biological resources (Chamberlain and Doering, 1998; Doering et al., 2002, 2006; SFWMD, 2010). Effective management of watershed and estuary resources requires the ability to link freshwater inflows to downstream ecological responses. Mathematical models accommodate this linkage by integrating many sources of information while providing a platform to conduct exploratory experiments (Sohma et al., 2008; Swaney et al., 2008; Buzzelli, 2011; Condie et al., 2012). This study resulted in a unique framework that combines empirically-derived inputs of freshwater and materials from the watershed, predicted salinity distributions, box modeling of estuarine transport, biogeochemical simulation modeling, and in situ time series for calibration (Hagy and Murrell, 2007; Buzzelli, 2008; Buzzelli et al., 2012; Wan et al., 2013). This approach was highly effective because it 1) utilized local monitoring data, 2) implemented a straightforward and transferrable physical model, and 3) afforded the capability to execute simulations in seconds on a desk top computer. These attributes expedite model parameterization, sensitivity analyses, and the application to local and regional freshwater inflow scenarios. However, model conceptualization, equation structure, parameterization, and uncertainty must be explored before they can be incorporated into management schemes (James et al., 2005). While the assessment of model uncertainty can be a specialized topic (Borsuk et al., 2004; Lehrter and Cebrian, 2010), general testing of the present model demonstrated that it has the potential to be an effective tool for the management of watershed and estuary resources. The performance and reliability of simulation models are dependent upon the availability of data on biogeochemical rate processes and concentrations. Unfortunately, insufficient and highly variable data exist for the CRE (Buzzelli et al., 2013c). While monthly watershed concentration data were adequate for use as external drivers, data on submarine light, water column rate processes, and biogeochemical cycling at the sediment-water interface within the CRE were lacking. Despite its regional importance and suggested sensitivity to nutrient loading, it is unfortunate that phytoplankton and nutrient turnover in the CRE has not been quantified more extensively (Buzzelli et al., 2013d). The few available measured rates of gross primary production and denitrification had great variability at the estuary scale (Howes et al., 2008; Phlips and Mathews, 2009). As opposed to other locations, the absence of information at the sediment-water interface necessitated rudimentary formulations for the vertical exchanges of particulate and dissolved materials (Smits and van Beek, 2013; Brady et al., 2013). There were only ~3.5 years of light penetration data and ~4 years of seagrass biomass estimates at four locations. Recent measurements of light quantity and quality are being examined in greater detail (Z. Chen and P. Doering, unpublished data). Monitoring of seagrass habitats since 2008 has focused on community composition and not biomass (SFWMD, 2012). Thus, model equations and coefficients were largely derived from literature sources. The CRE responds to external DIN loading by increasing rates of primary production, remineralization, and DIN concentrations (Buzzelli et al., 2013d). However, DIN accounts for only ~30% of total N loading to the estuary with the dynamics of dissolved ON remaining largely unknown (Smith and Hollibaugh, 2006; Eyre

C. Buzzelli et al. / Estuarine, Coastal and Shelf Science 151 (2014) 256e271

et al., 2011). Based on the correlation coefficients, the model was a better predictor of ON in the upper estuary during wet conditions because dissolved and particulate organic nitrogen are more associated with low salinity (Loh, 2008). The reactivity of allochthonous vs. autochthonous ON sources has implications for phytoplankton dependence upon external vs. internal supplies of DIN. Whereas þ NO x concentrations were directly related to freshwater inflow, NH4 was less predictable due to the prevalence of internal cycling in the dry season (Buzzelli et al., 2013b). Although CHL concentrations were variable, the comparatively simple model proved to be a reasonable predictor in time and space. It is possible that model over-estimation of chlorphyll a concentrations in the downstream segment 3 during the wet season affects predicted partitioning of submarine light extinction (Buzzelli et al., in press). The overall agreement between model chlorphyll a and the monitoring data was similar to that from other studies (Cerco, 2000;Muylaert et al., 2005; Sohma et al., 2008; Condie et al., 2012; Smits and van Beek, 2013). Differences with other studies were primarily because the CRE model was derived for only three segments, a single vertical layer, and relied on comparatively fewer rate-process coefficients. The model was sensitive to variations in both kRphyto and kDIN. The basal rate of respiration was 6% of the maximum growth rate of phytoplankton (Lucas et al., 2009b). Decreasing the respiration rate greatly increased phytoplankton biomass, gross production, and the uptake of inorganic N from the water column. Reducing kDIN relieved nitrogen limitation leading to increased gross production and DIN removal. Adjusting these coefficients offers a potential pathway to improve model fitness. However, seasonal differences in predicted CHL concentrations (dry ¼ underestimate; wet ¼ overestimate) likely resulted because phytoplankton were assumed to be equally  þ influenced by NHþ 4 and NOx . Differential preference for NH4 over NO would emphasize internal sources of DIN thereby augmenting x CHL in the dry season and mitigating it in the wet season. This offers a potentially realistic adjustment that could be made to the CHL and DIN formulations. However, it is difficult to justify in the absence of empirical or laboratory information on water column carbon and nitrogen transformations in the CRE. The model accounted for the magnitude but not necessarily the timing of Halodule wrightii maximum biomass. Conversely, the timing of Thalassia testudinum growth was captured although the maximum biomass did not exactly match the variable field observations. Overall, the model reflected variations in seagrass biomass despite the simplicity of the model formulations. Intra- and interannual variations in seagrass biomass are functions of complex physiological responses to submarine light, salinity, nutrient availability, and epiphyte characteristics (Doering et al., 2002; Lee et al., 2007; Kahn and Durako, 2008; Dawes and Avery, 2010; Herbert et al., 2011). Seagrasses, especially those in oligotrophic environments, have the ability to translocate reserves of C, N, and P internally among shoots and root-rhizomes and between individual shoots at the patch scale (Burd and Dunton, 2001; Marba et al., 2002). Additionally, seagrasses have the capacity to access sediment DON supplies to satisfy nitrogen requirements (Vonk et al., 2008). Future improvements to the seagrass sub-model include development of below-ground components, nutrient uptake, bidirectional translocation, and epiphyte effects (Buzzelli et al., 1999; Buzzelli et al., in press). Presently, the model is being used to quantify effects of light vs. salinity on the growth and survival of seagrasses in the lower CRE (Buzzelli et al., in press). This modeling framework was designed to link watershed inputs with ecological indicators for better management of both freshwater and estuarine resources. Such a model is required because the inter-relationships among water temperature, freshwater inflow, salinity, flushing, nutrient loading, primary

269

production, and submarine light are non-linear in sub-tropical estuaries such as the CRE (Buzzelli, 2011). The goal was to quantify the independent and cumulative effects of multiple environmental drivers on water quality and seagrass survival in order to define optimal inflows, nutrient loads, and salinity distributions. This integration of watershed inputs, salinity gradients, and biogeochemical simulations provides the platform to conduct mathematical experiments to help reach this goal. In part II of this study (Buzzelli et al., in press) the calibrated model was used to explore three inflow-related questions for the CRE. First, how are variations in freshwater inflow and phytoplankton biomass related to patterns of submarine light extinction in time and space (McPherson et al., 2011; Buzzelli et al., 2012; Le et al., 2013)? The model was used to investigate the relative importance of CDOM vs. CHL in submarine light availability throughout the CRE (McPherson and Miller, 1987; Gallegos and Kenworthy, 1996; Livingston et al., 1998; Christian and Sheng, 2003). Second, if reductions in external nutrient loads are desirable, should management aim to regulate N and P inputs by controlling discharge or inflow concentrations? This question is directly relevant for planned water reservoirs vs. treatment wetlands in the upper watershed of the CRE (SFWMD, 2012). Third, what is the range of freshwater inflows that benefit submarine light conditions and seagrass survival in the lower CRE? The model provides an expeditious and efficient approach to address these important questions in a quantitative capacity. 5. Conclusions  This study resulted in a framework that combines empiricallyderived inputs of freshwater and materials from the watershed, predicted salinity distributions, box modeling of estuarine transport, biogeochemical simulation modeling, and in situ time series for calibration.  The correlations between the model representation and empirically derived concentrations of organic nitrogen and phosphorus, inorganic nitrogen and phosphorus, and chlorophyll a varied for each constituent with estuary segment (1,2,3), season (dry vs. wet), and annual freshwater inflow (2002e2009).  The capability to connect watershed inputs with ecological indicators helps with the management of both freshwater and estuarine resources because temperature and seasonality, freshwater inflow, salinity, flushing, nutrient loading, primary production, and submarine light are inter-related in sub-tropical estuaries such as those in south Florida.  The modeling framework is being used to generate an experimental series of simulations to quantify potential effects of variable inflow and nutrient loads on submarine light penetration and seagrass survival in the Caloosahatchee River Estuary. Acknowledgments We would like to acknowledge R. Thomas James for input during model development and constructive comments on the manuscript. Karen Bickford at the Lee County Natural Resource Office provided the concentration of water column constituents in the tributaries. The manuscript greatly benefitted from comments provided by Susan Gray, Linda Lindstrom, Zhiqiang Chen, Cassondra Thomas, Patricia Gorman, and two anonymous reviewers. References Alber, M., 2002. A conceptual model of estuarine freshwater inflow management. Estuaries 25, 1246e1261.

270

C. Buzzelli et al. / Estuarine, Coastal and Shelf Science 151 (2014) 256e271

Antonini, G.A., Fann, D.A., Roat, P., 2002. A Historical Geography of Southwest Florida Waterways. In: Placida Harbor to Marco Island, vol. 2. National Seagrant College Program, Silver Spring, Maryland. Barnes, T., 2005. Caloosahatchee estuary conceptual ecological model. Wetlands 25, 884e897. Bjornsen, P.K., 1988. Phytoplankton exudation of organic matter: why do healthy cells do it? Limnol. Oceanogr. 33 (1), 151e154. Borsuk, M.E., Stow, C.A., Reckhow, K.H., 2004. A Bayesian network of eutrophication models for synthesis, prediction, and uncertainty analysis. Ecol. Model. 173, 219e239. Boynton, W.R., Hagy, J.D., Cornwell, J.C., Kemp, W.M., Greene, S.M., Owens, M.S., Baker, J.E., Larsen, R.K., 2008. Nutrient budgets and management actions in the Patuxent River Estuary, Maryland. Estuaries Coasts 31, 623e651. Bowers, D.G., Brett, H.L., 2008. The relationship between CDOM and salinity in estuaries: an analytical and graphical solution. J. Mar. Syst. 73, 1e7. Brady, D.C., Testa, J.M., DiToro, D.M., Boynton, W.R., Kemp, W.M., 2013. Sediment flux modeling: calibration and application for coastal systems. Estuar. Coast. Shelf Sci. 117, 107e124. Burd, A.B., Dunton, K.H., 2001. Field verification of a light-driven model of biomass changes in the seagrass Halodule wrightii. Mar. Ecol. Prog. Ser. 209, 85e98. Buzzelli, C., Wetzel, R.L., Meyers, M.B., 1999. A linked physical and biological framework to assess biogeochemical dynamics in a shallow estuarine ecosystem. Estuar. Coast. Shelf Sci. 49, 829e851. Buzzelli, C., 2008. Development and application of tidal creek ecosystem models. Ecol. Model. 210, 127e143. Buzzelli, C., 2011. Ecosystem modeling in small sub-tropical estuaries and embayments. In: Wolanski, E., McLusky, D.S. (Eds.), Treatise on Estuarine and Coastal Science, vol. 9. Academic Press, Waltham, Massachusetts, pp. 331e353. Buzzelli, C., Robbins, R., Doering, P., Chen, Z., Sun, D., Wan, Y., Welch, B., Schwarzchild, A., 2012. Monitoring and modeling of Syringodium filiforme (Manatee Grass) in the southern Indian River Lagoon. Estuaries Coasts 35, 1401e1415. Buzzelli, C., Parker, M., Geiger, S., Wan, Y., Doering, P., Haunert, D., 2013a. Predicting system-scale impacts of oyster clearance on phytoplankton production in a small sub-tropical estuary. Environ. Model. Assess. 18, 185e198. Buzzelli, C., Chen, Z., Coley, T., Doering, P., Samimy, R., Schlezinger, D., Howes, B., 2013b. Dry season sediment-water exchanges of nutrients and oxygen in two Florida estuaries: patterns, comparisons, and internal loading. Fla. Sci. 76 (1), 54e79. Buzzelli, C., Doering, P.D., Wan, Y., Gorman, P., Volety, A., 2013c. Simulation of potential oyster density with variable freshwater inflow (1965e2001) to the Caloosahatchee River Estuary, southwest Florida, USA. Environ. Manag. 52, 981e994. Buzzelli, C., Wan, Y., Doering, P., Boyer, J.N., 2013d. Seasonal dissolved inorganic nitrogen and phosphorus budgets for two sub-tropical estuaries in south Florida, USA. Biogeosciences 10, 6721e6736. Buzzelli, C., Boutin, B., Ashton, M., Welch, B., Gorman, P., Wan, Y., Doering, P., 2014a. Fine-scale detection of estuarine water quality with managed freshwater releases. Estuaries Coasts 37, 1134e1144. Buzzelli, C., Doering, P., Wan, Y., Sun, D., Chen, Z., 2014b. Modeling ecosystem processes with variable freshwater inflow to the Caloosahatchee River Estuary, southwest Florida. II. Inflow and nutrient loading, submarine light, and seagrasses. Estuar. Coast. Shelf Sci. (in press). Chamberlain, R.H., Doering, P.H., 1998. Freshwater inflow to the Caloosahatchee Estuary and the resource-based method for evaluation. In: Treat, S.F. (Ed.), Proceedings of the Charlotte Harbor Public Conference and Technical Symposium. South Florida Water Management District, West Palm Beach, FL, pp. 81e90. Childers, D.L., Boyer, J.N., Davis, S.E., Madden, C.J., Rudnick, D.T., Sklar, F., 2006. Relating precipitation and water management to nutrient concentrations in the oligotrophic “upside-down” estuaries of the Florida Everglades. Limnol. Oceanogr. 51, 602e616. Christian, D., Sheng, Y.P., 2003. Relative influence of various water quality parameters on light attenuation in Indian River Lagoon. Estuar. Coast. Shelf Sci. 57, 961e971. Cerco, C.F., Seitzinger, S.P., 1997. Measured and modeled effects of benthic algae on eutrophication in Indian River-Rehoboth Bay, Delaware. Estuaries 20 (1), 231e248. Cerco, C.F., 2000. Phytoplankton kinetics in the Chesapeake Bay eutrophication model. Water Qual. Ecosyst. Model. 1, 5e49. Cerco, C.F., Moore, K.A., 2001. System-wide submerged aquatic vegetation model for Chesapeake Bay. Estuaries 24 (4), 522e534. Cerco, C.F., Noel, M.R., 2007. Can oyster restoration reverse cultural eutrophication in Chesapeake Bay? Estuaries Coasts 30 (2), 331e343. Cloern, J.E., 2001. Our evolving conceptual model of the coastal eutrophication problem. Mar. Ecol. Prog. Ser. 210, 223e253. Condie, S.A., Hayes, D., Fulton, E.A., Savina, M., 2012. Modelling ecological change over a half a century in a subtropical estuary: impacts of climate change, land-use, urbanization, and freshwater extraction. Mar. Ecol. Prog. Ser. 457, 43e66. Cornwell, J.C., Kemp, W.M., Kana, T.M., 1999. Denitrification in coastal ecosystems: methods, environmental controls, and ecosystem level controls, a review. Aquat. Ecol. 33, 41e54. Dawes, C., Avery, W., 2010. Epiphytes of the seagrass Halodule wrightii in Hillsborough Bay, Florida: a 14 year study in an estuary recovering from eutrophication. Fla. Sci. 73 (3/4), 185e195.

Dettmann, E.H., 2001. Effect of water residence time on annual export and denitrification of nitrogen in estuaries: a model analysis. Estuaries 24 (4), 481e490. Doering, P.H., Chamberlain, R.H., 2000. Experimental studies on the salinity tolerance of turtle grass, Thalassia testudinum. In: Bortone, S.A. (Ed.), Seagrasses: Monitoring Ecology, Physiology, and Management. CRC Press, Boca Raton, FL, pp. 81e98. Doering, P.H., Chamberlain, R.H., Haunert, D., 2002. Using submersed aquatic vegetation to establish minimum and maximum freshwater inflows to the Caloosahatchee Estuary, Florida. Estuaries 25, 1343e1354. Doering, P.H., Chamberlain, R., Haunert, K.M., 2006. Chlorophyll a and its use as an indicator of eutrophication in the Caloosahatchee Estuary, Florida. Fla. Sci. 69, 51e72. Duarte, C.M., 1990. Seagrass nutrient content. Mar. Ecol. Prog. Ser. 67, 201e207. Duarte, C.M., Conley, D.J., Carstensen, J., Sanchez-Camacho, M., 2009. Return to Neverland: shifting baselines affect eutrophication restoration targets. Estuaries Coasts 32, 29e36. Eldridge, P.M., Kaldy, J.E., Burd, A.B., 2004. Stress response model for the tropical seagrass Thalassia testudinum: the interactions of light, temperature, sedimentation, and geochemistry. Estuaries 27, 923e937. Eyre, B., Ferguson, A.J.P., Webb, A., Maher, D., Oakes, J.M., 2011. Denitrification, Nfixation, and nitrogen and phosphorus fluxes in different benthic habitats and their contribution to the nitrogen and phosphorus budgets of a shallow oligotrophic sub-tropical coastal system (southern Moreton Bay, Australia). Biogeochemistry 102, 111e133. Flemer, D.A., Champ, M.A., 2006. What is the future fate of estuaries given nutrient over-enrichment, freshwater diversion, and low flows? Mar. Pollut. Bull. 5, 247e258. Gallegos, C.L., 2001. Calculating optical water quality targets to restore and protect submersed aquatic vegetation: overcoming problems in partitioning the diffuse attenuation coefficient for photosynthetically active radiation. Estuaries 24, 381e397. Gallegos, C.L., Kenworthy, W.J., 1996. Seagrass depth limits in the Indian River Lagoon (Florida, U.S.A.): application of an optical water quality model. Estuar. Coast. Shelf Sci. 42, 267e288. Gillson, J., 2011. Freshwater flow and fisheries production in estuarine and coastal systems: where a drop of rain is not lost. Rev. Fish. Sci. 19 (3), 168e186. Gonzalez-Ortegon, E., Drake, P., 2012. Effects of freshwater inputs on the lower trophic levels of a temperate estuary: physical, physiological or trophic forcing? Aq. Sci. 74, 455e469. Greening, H.S., Cross, L.M., Sherwood, E.T., 2011. A multi-scale approach to seagrass recovery in Tampa Bay, Florida. Ecol. Restor. 29 (1e2), 82e93. Hagy, J.D., Boynton, W.R., Sanford, L.P., 2000. Estimation of net physical transport and hydraulic residence times for a coastal plain estuary using box models. Estuaries 23 (3), 328e340. Hagy, J.D., Murrell, M.C., 2007. Susceptibility of a northern Gulf of Mexico estuary to hypoxia: an analysis using box models. Estuar. Coast. Shelf Sci. 74, 239e253. Herbert, D.A., Perry, W.B., Cosby, B.J., Fourqurean, J.W., 2011. Projected reorganization of Florida Bay seagrass communities to the increased freshwater inflow of Everglades restoration. Estuaries Coasts 34, 973e992. Herzka, S.Z., Dunton, K.H., 1997. Seasonal photosynthetic patterns of the seagrass Thalassia testudinum in the western Gulf of Mexico. Mar. Ecol. Prog. Ser. 152, 103e117. Howes, B.H., Schlezinger, D., Samimy, R., 2008. The Characterization and Quantification of Benthic Nutrient Fluxes in the Caloosahatchee River Estuary. Project Summary Report. Coastal Systems Program, School for Marine Science and Technology, University of Massachusetts Dartmouth, p. 56. James, R.T., Bierman, V.J., Erickson, M.J., Hinz, S.C., 2005. The Lake Okeechobee water quality model (LOWQM) enhancements, calibration, validation, and analysis. Lake Reserv. Manag. 21 (3), 231e260. Kahn, A., Durako, M., 2008. Photophysiological responses of Halophila johnsonii to experimental hyposaline and hyper-CDOM conditions. J. Exp. Mar. Biol. Ecol. 376, 230e235. Kaul, L.W., Froelich, P.N., 1984. Modeling estuarine nutrient geochemistry in a simple system. Geochimica Cosmochimica Acta 48, 1417e1433 dr. Kemp, W.M., Boynton, W.R., Adolf, J.E., Boesch, D.F., Boicourt, W.C., Brush, G., Cornwell, J.C., Fisher, T.R., Glibert, P.M., Hagy, J.D., Harding, L.W., Houde, E.D., Kimmel, D.G., Miller, W.D., Newell, R.I.E., Roman, M.R., Smith, E.M., Stevenson, J.C., 2005. Eutrophication of Chesapeake Bay: historical trends and ecological interactions. Mar. Ecol. Prog. Ser. 303, 1e29. Kremer, J.N., Vaudrey, J.M.P., Ullman, D.S., Bergondo, D.L., LaSota, N., Kincaid, C., Codiga, D.L., Brush, M.J., 2010. Simulating property exchange in estuarine ecosystem models at ecologically appropriate scales. Ecol. Model. 221, 1080e1088. Kuo, A.Y., Park, K., 1994. A PC-based Tidal Prism Water Quality Model for Small Coastal Basins and Tidal Creeks. Spec. Report No. 324 in Applied Marine Science and Ocean Engineering, College of William and Mary, Gloucester Point, VA, p. 119. Larsson, U., Hagstrom, A., 1979. Phytoplankton exudate release as an energy source for the growth of pelagic bacteria. Mar. Biol. 52, 199e206. Le, C., Hu, C., English, D., Cannizzaro, J., Chen, Z., Kovach, C., Anastasiou, C.J., Zhao, J., Carder, K., 2013. Inherent and apparent optical properties of the complex estuarine waters of Tampa Bay: what controls light? Estuar. Coast. Shelf Sci. 117, 54e69. Lee, K.-S., Park, S.R., Kim, Y.K., 2007. Effects of irradiance, temperature, and nutrients on growth dynamics of seagrasses: a review. J. Exp. Mar. Biol. Ecol. 350, 144e175.

C. Buzzelli et al. / Estuarine, Coastal and Shelf Science 151 (2014) 256e271 Lehrter, J.C., Cebrian, J., 2010. Uncertainty propagation in an ecosystem nutrient budget. Ecol. Appl. 20 (2), 508e524. Livingston, R.J., Niu, X., Lewis, F.G., Woodsum, G.C., 1997. Freshwater input to a Gulf estuary: long term control of trophic organization. Ecol. Appl. 7, 277e299. Livingston, R.J., McGlynn, S.E., Niu, X., 1998. Factors controlling seagrass growth in a gulf coastal system: water and sediment quality and light. Aquat. Bot. 60, 135e159. Loh, A., 2008. Mixing and Degradation of Riverine Dissolved Organic Nitrogen in the Caloosahatchee Estuary. South Florida Water Management District, West Palm Beach, FL, p. 35. Lucas, L.V., Thompson, J.K., Brown, L.R., 2009a. Why are diverse relationships observed between phytoplankton biomass and transport time? Limnol. Oceanogr. 54, 381e390. Lucas, L.V., Koseff, J.R., Monismith, S.G., Thompson, J.K., 2009b. Shallow water processes govern system-wide phytoplankton bloom dynamics: a modeling study. J. Mar. Syst. 75, 70e86. Marba, N., Hemminga, M.A., Mateo, M.A., Duarte, C.M., Mass, Y.E.M., Terrados, J., Gacia, E., 2002. Carbon and nitrogen translocation between seagrass ramets. Mar. Ecol. Prog. Ser. 226, 287e300. Mazzotti, F.J., Pearlstine, L.G., Chamberlain, R., Barnes, T., Chartier, K., DeAngelis, D., 2007. Stressor-response Models for Seagrasses, Halodule Wrightii and Thalassia Testudinum. Final report to the South Florida Water Management District and the U.S. Geological Survey. University of Florida, Florida Lauderdale Research and Education Center, Fort Lauderdale, Florida, p. 19. McPherson, B.F., Miller, R.L., 1987. The vertical attenuation of light in Charlotte Harbor, a shallow, subtropical estuary, southwestern Florida. Estuar. Coast. Shelf Sci. 25, 721e737. McPherson, M.L., Hill, V.J., Zimmerman, R.C., Dierssen, H.M., 2011. The optical properties of greater Florida Bay: implications for seagrass abundance. Estuaries Coasts 34, 1150e1160. Miller, R.L., McPherson, B.F., 1991. Estimating estuarine flushing and residence time in Charlotte Harbor, Florida, via salt balance and box model. Limnol. Oceanogr. 36 (3), 602e612. Morey, S.L., Dukhovskoy, D.S., 2012. Analysis methods for characterizing salinity variability from multivariate time series applied to the Apalachicola Bay estuary. J. Atmos. Ocean. Technol. 29, 613e628. Muylaert, K., Tackx, M., Vyverman, W., 2005. Phytoplankton growth rates in the freshwater tidal reaches of the Schelde Estuary (Belgium) estimated using a simple light-limited primary production model. Hydrobiologia 540, 127e140. Nixon, S.W., Fulweiler, R.W., Buckley, B.A., Granger, S.L., Nowicki, B.L., Henry, K.M., 2009. The impact of changing climate on phenology, productivity, and benthicpelagic coupling in Narragansett Bay. Estuar. Coast. Shelf Sci. 82, 1e18. Officer, C., 1980. Box models revisited. In: Hamilton, P., Macdonald, R.B. (Eds.), Estuarine and Wetland Processes, Marine Sciences Series, vol. 11. Plenum Press, NY, pp. 65e114. Petes, L.E., Brown, A.J., Knight, C.R., 2012. Impacts of upstream drought and water withdrawals on the health and survival of downstream estuarine oyster populations. Ecol. Evol. 2 (7), 1712e1724. Phlips, E.J., Mathews, L., 2009. Primary Production in the Caloosahatchee River/San Carlos Bay Ecosystem (Final report to the South Florida Water Management District). Phlips, E.J., Badylak, S., Hart, J., Haunert, D., Lockwood, J., O'Donnell, K., Sun, D., Viveros, P., Yilmaz, M., 2012. Climatic influences on autochthonous and allochthonous phytoplankton blooms in a subtropical estuary, St. Lucie Estuary, Florida, USA. Estuaries Coasts 35 (1), 335e352. Pinckney, J., 1994. Development of an irradiance-based ecophysiological model intertidal benthic microalgal production. In: Krumbein, W.E., Paterson, D.M., Stal, L.J. (Eds.), Biostabilization of Sediments. Oldenburg, pp. 55e83.

271

Powell, E.N., Klinck, J.M., Hofmann, E.E., McManus, M.A., 2003. Influence of water allocation and freshwater inflow on oyster production: a hydrodynamic-oyster population model for Galveston Bay, Texas, USA. Environ. Manag. 31 (1), 100e121. Qiu, C., Wan, Y., 2013. Time series modeling and prediction of salinity in the Caloosahatchee River Estuary. Water Resour. Res. 49, 1e13. Robson, B.J., 2005. Representing the effects of diurnal variations in light on primary production on a seasonal time-scale. Ecol. Model. 186, 358e365. Sackett, J.W., 1888. Survey of the Caloosahatchee River, Florida. Report to the Captain of the United States Engineering Office, St. Augustine, FL. Scavia, D., Liu, Y., 2009. Exploring estuarine nutrient susceptibility. Environ. Sci. Technol. 43 (10), 3474e3479. SFWMD, September 16, 2010. Final Protocols for Lake Okeechobee Operations. South Florida Water Management District. SFWMD, 2012. South Florida Environmental Report. Caloosahatchee River Watershed Protection Plan 2012 Update, vol. 1. Appendix 10e2. Available on-line at: www.sfwmd.gov/SFER. Sheldon, J.E., Alber, M., 2006. The calculation of estuarine turnover times using freshwater fraction and tidal prism methods: a critical evaluation. Estuaries Coasts 29, 133e146. Sin, Y., Wetzel, R.L., 2002. Ecosystem modeling analysis of size-structured phytoplankton dynamics in the York River Estuary, Virginia (USA). I. Development of a plankton ecosystem model with explicit feedback controls and hydrodynamics. Mar. Ecol. Prog. Ser. 228, 75e90. Sklar, F.H., Browder, J.A., 1998. Coastal environmental impacts brought about by alterations to freshwater flow in the Gulf of Mexico. Environ. Manag. 22 (4), 547e562. Smith, S.V., Hollibaugh, J.T., 2006. Water, salt, and nutrient exchanges in San Francisco Bay. Limnol. Oceanogr. 51 (1/2), 504e517. Smits, J.G.C., van Beek, J.K.L., 2013. ECO: a generic eutrophication model including comprehensive sediment-water interaction. PLoS ONE 8 (7), e68104, 1e25. Sohma, A., Sekiguchi, Y., Kuwae, T., Nakamura, Y., 2008. A benthic-pelagic coupled ecosystem model to estimate the hypoxic estuary including tidal flat e model description and validation of seasonal/daily dynamics. Ecol. Model. 215, 10e39. Swaney, D., Scavia, D., Howarth, R.W., Marino, R.M., 2008. Estuarine classification and response to nitrogen loading: Insights from simple ecological models. Estuar. Coast. Shelf Sci. 77, 253e263. Swaney, D., Smith, S.V., Wulff, F., 2011. The LOICZ biogeochemical modeling protocol and its application to estuarine ecosystems. In: Wolanski, E., McLusky, D.S. (Eds.), Treatise on Estuarine and Coastal Science, vol. 9. Academic Press, Waltham, Massachusetts, pp. 135e160. Taylor, D.I., Oviatt, C.A., Borkman, D.G., 2011. Non-linear responses of a coastal aquatic ecosystem to large decreases in nutrient and organic loadings. Estuaries Coasts 34, 745e757. Tolley, S.G., Volety, A.K., Savarese, M., 2005. Influence of salinity on the habitat use of oyster reefs in three southwest Florida estuaries. J. Shellfish Res. 24 (1), 127e137. Vonk, J.A., Middleburg, J.J., Stapel, J., Bouma, T.J., 2008. Dissolved organic nitrogen uptake by seagrasses. Limnol. Oceanogr. 53 (2), 542e548. Walsh, J.J., Jolliff, J.K., Darrow, B.P., Lenes, J.M., Milroy, S.P., Remsen, A., Dieterle, D.A., Carder, K.L., Chen, F.R., Vargo, G.A., Weisberg, R.H., Fanning, K.A., MullerKarger, F.E., Shinn, E., Steidinger, K.A., Heil, C.A., Tomas, C.R., Prospero, J.S., Lee, T.N., Kirkpatrick, G.J., Whitledge, T.E., Stockwell, D.A., Vilareal, T.A., Jochens, A.E., Bontempi, P.S., 2006. Red tides in the Gulf of Mexico: where, when, and why? J. Geophys. Res. 111 (C11003), 1e46. Wan, Y., Qiu, C., Doering, P., Ashton, M., Sun, D., Coley, T., 2013. Modeling residence time with a three-dimensional hydrodynamic model: linkage with chlorophyll a in a sub-tropical estuary. Ecol. Model. 268, 93e102. Wiegert, R.G., 1975. Simulation models of ecosystems. Annu. Rev. Ecol. Syst. 6, 311e338.