Seasonal and Inter-annual Patterns in Primary Production, Respiration, and Net Ecosystem Metabolism in Three Estuaries in the Northeast Gulf of Mexico

Measurements of primary production and respiration provide fundamental information about the trophic status of aquatic ecosystems, yet such measurements are logistically difficult and expensive to sustain as part of long-term monitoring programs. However, ecosystem metabolism parameters can be inferred from high frequency water quality data collections using autonomous logging instruments. For this study, we analyzed such time series datasets from three Gulf of Mexico estuaries: Grand Bay, MS; Weeks Bay, AL; and Apalachicola Bay, FL. Data were acquired from NOAA's National Estuarine Research Reserve System Wide Monitoring Program and used to calculate gross primary production (GPP), ecosystem respiration (ER), and net ecosystem metabolism (NEM) using Odum's open water method. The three systems represent a diversity of estuaries typical of the Gulf of Mexico region, varying by as much as two orders of magnitude in key physical characteristics, such as estuarine area, watershed area, freshwater flow, and nutrient loading. In all three systems, GPP and ER displayed strong seasonality, peaking in summer and being lowest during winter. Peak rates of GPP and ER exceeded 200 mmol O2 m−2 day−1 in all three estuaries. To our knowledge, this is the first study examining long-term trends in rates of GPP, ER, and NEM in estuaries. Variability in metabolism tended to be small among sites within each estuary. Nitrogen loading was highest in Weeks Bay, almost two times greater than that in Apalachicola Bay and 35 times greater than to Grand Bay. These differences in nitrogen loading were reflected in average annual GPP rates, which ranged from 825 g C m−2 year−1 in Weeks Bay to 401 g C m−2 year−1 for Apalachicola Bay and 377 g C m−2 year−1 in Grand Bay. Despite the strong inter-annual patterns in freshwater flow and salinity, variability in metabolic rates was low, perhaps reflecting shifts in the relative importance of benthic and phytoplankton productivity, during different flow regimes. The advantage of the open water method is that it uses readily available and cost-effective sonde monitoring technology to estimate these fundamental estuarine processes, thus providing a potential means for examining long-term trends in net carbon balance. It also provides a historical benchmark for comparison to ongoing and future monitoring focused on documenting the effect of human activities on the coastal zone.


Introduction
Primary production, particularly phytoplankton production, forms a critical base to estuarine food webs (Iverson 1990, Nixon and Buckley 2002, Day et al. 2013) and strongly influences oxygen and nutrient dynamics. Nutrient loading related to human population growth and intensive agriculture in estuarine drainage basins stimulates phytoplankton and macroalgal production, causing eutrophication, which has been recognized as a significant problem in estuarine and coastal waters for decades (Nixon 1995;NRC 2000;Cloern 2001). Generally, nutrient concentrations in Gulf of Mexico estuaries tend to be low compared to more well studied, urbanized estuaries (Pennock et al. 1999), although population growth in the Gulf of Mexico region is placing increasing pressures on the coastal zone (Twilley et al. 2001). In addition to alteration of nutrient loading, development increases impervious surface area, generally causing greater and flashier runoff patterns (Holland et al. 2004). Another consequence of increasing urbanization is the increased freshwater demand, which reduces downstream water supply (Benson 1981, Petes et al. 2012. Thus, hydrologic modifications of the watersheds that alter the timing and volume of freshwater flow and nutrient loads to the estuary can impact estuarine-dependent fisheries. Light availability represents another critical factor controlling primary production (Cloern et al. 1985, Pennock and Sharp 1994, Cloern 1999. The estuaries of the northern Gulf of Mexico region, all located at ∼30°N, have high solar insolation, which can support relatively high primary production throughout the year (Pennock et al. 1999). Benthic primary production is significant (Moncreiff et al. 1992, Schreiber and Pennock 1995, Murrell et al. 2009), owing to the shallow nature of these systems.
Ecosystem respiration is largely microbial activity, thus influenced by seasonal temperature patterns and seasonal patterns in the supply and quality of organic matter (Pomeroy and Weibe 2001;Murrell 2003;Hopkinson and Smith 2005). In a seminal review of the estuarine literature, Hopkinson and Smith (2005) found temperature explained 28 % of the variability in water column respiration rates. However, in systems with small allochthonous organic matter inputs, chlorophyll a concentrations tended to be a better predictor of respiration rates than seasonal patterns (Hopkinson and Smith 2005).
Net ecosystem metabolism (NEM), the balance between gross primary production (GPP), and ecosystem respiration (ER), provides an integrated measure of the trophic state of the estuary, indicating whether the system is accumulating or depleting organic matter within the system and whether there is net uptake or release of CO 2 Hollibaugh 1993, Cai 2011). Net autotrophy, (i.e., GPP>ER) implies organic carbon accumulation and net CO 2 uptake from the atmosphere. Conversely, net heterotrophy implies a depletion of the organic C pools and a net CO 2 release to the atmosphere. In estuarine environments, net heterotrophy has been commonly observed, and heterotrophy appears to be more prevalent in systems with large allochthonous organic matter inputs to the system (Kemp et al. 1997) particularly from adjacent marshes and mangroves (Caffrey 2004). Conversely, systems with high dissolved inorganic nitrogen (DIN) inputs tend towards autotrophy relative to systems with low DIN inputs (Kemp et al. 1997, Caffrey 2004.
The goal of this study was to examine seasonal and interannual variation in estuarine ecosystem metabolism and the interactions between metabolism and environmental variables such as freshwater flow, temperature and nutrients. Because Gulf of Mexico estuaries are relatively poorly known, this study aims to provide fundamental information about how these estuaries function, and serves as a useful benchmark for comparison to more well-studied systems.

Study Sites
We examined three estuaries along the northern Gulf of Mexico coastline: Apalachicola Bay, Weeks Bay, and Grand Bay (Fig. 1). These three estuaries are part of the National Estuarine Research Reserve System across the USA, which is based on partnership between the National Oceanic and Atmospheric Administration and a respective State agency (http://www.nerrs.noaa.gov.). All reserves participate in the System Wide Monitoring Program (SWMP), which supports monitoring of water quality and weather for all Reserves in the system. SWMP also provides support for consistent data quality assurance and archival through its Centralized Data Management Office (CDMO: http://cdmo.baruch.sc.edu).
The Apalachicola Bay NERR was designated in 1979 with the Florida Department of Environmental Protection as the state partner. It is the second largest NERR in the country, with extensive open water, submerged and emergent wetland vegetation, tidal flats, and unconsolidated bottom The upland areas of the Reserve include forested floodplain and barrier island (Little St. George Island) habitat (ANERR, 1998). Much of the area adjacent to the Reserve and Apalachicola Bay (over 85 % of Franklin County) is also publicly owned or held in the public trust, including the Apalachicola National Forest and the Tate's Hell State Forest (Edmiston, 2008).
Outside of public lands, adjacent land use includes agriculture, silviculture, commercial, and residential uses.
Apalachicola Bay is ∼450 km 2 in area and lies at the terminus of the Apalachicola River, itself formed by the confluence of the Chattahoochee and Flint Rivers at the Florida-Georgia border. The Apalachicola Bay watershed is ∼52,000 km 2 and originates in northern Georgia (Table 1). The bay is bounded by a series of barrier islands with several inlets that exchange with the Gulf of Mexico. Water residence time in Apalachicola Bay is relatively short, ranging from ∼3 to 25 days, due to high freshwater flow and the relatively open exchange with the Gulf of Mexico (Mortazavi et al. 2000a). The timing and duration of flood events strongly influences the salinity distribution and nutrient loads to Apalachicola Bay (Mortazavi et al. 2000a;Pennock et al. 1999;Twilley et al. 1999). The river and bay support commercial and recreational fisheries of finfish, shrimp, and oysters. Apalachicola Bay provides 90 % of Florida's and 10 % of the US oyster harvest (FFWCC 2011;NMFS 2011). Critical management concerns include potential river diversions for urban and agriculture use that would change the salinity distributions and residence times, thus threatening the valuable oyster fishery (Table 1).
The Weeks Bay NERR was designated in 1986 and the Alabama Department of Conservation and Natural Resources became the State partner in 2000. The watershed lies within Baldwin County, AL, which has experienced rapid human population growth over the past 20 years. Urban areas increased 92.5 % from 1990 to 2000 with a 42.9 % increase in population (Cartwright 2002). Baldwin County attracts ∼1 million tourists annually and supports and active sport fishing industry (Baldwin County Wetlands Conservation Plan 2005).
The Weeks Bay estuary is a small (7 km 2 ), shallow (mean depth 1.4 m) sub-estuary of Mobile Bay (Schroeder et al. 1990;Miller-Way et al. 1996; Table 1). The Fish River is the major freshwater source, contributing 73 % of total freshwater input, with the remainder coming from the Magnolia River. The mouth of Weeks Bay empties into Mobile Bay proper and exchange between the two bays varies with river discharge, wind and tidal forcing (Miller-Way et al. 1996). Management concerns center on increased urban and agriculture runoff,  Table 2 leading to high nutrient concentrations in the Bay, the loss of wetland habitat and alteration of the food web from high nutrient concentrations (Table 1). One notable example was an extensive harmful algal bloom that occurred during a drought in 2007, including species such as: Prorocentrum minimum , Karlodinium veneficum , Kryptoperidinium folaceum , Skeletonema sp., Cylindrotheca , Akashiwo, Nitzschia, and Oscillatoria (S. Phipps, unpublished data). The Grand Bay NERR, located in southeastern Jackson County, MS, was established in 1999 with the Mississippi Department of Environmental Quality. It is a small and relatively pristine estuary (∼60 km 2 ) with a watershed area of ∼263 km 2 (Woodrey 2007a, CDMO 2012 Table 1). The upland areas are comprised of pine savanna forest, maritime forest, coastal bayhead and cypress swamps, including rural residential and several industrial facilities along the western edge of the Reserve near Bangs Lake. Like all estuaries along this portion of the Gulf of Mexico coastline, the tides are predominately diurnal and of low amplitude (<0.6 m). Wind forcing is therefore is a relatively strong driver of short-term variability in water level. The estuary opens to Mississippi sound and has no major river inputs (Woodrey 2007b); thus, the estuary tends to be more marine in character compared to Weeks Bay and Apalachicola Bay. Freshwater input into the estuary is primarily local runoff from bayous and tidal creeks, including Bayou Cumbest, Bayou Heron, and Bangs Lake (Woodrey 2007b). Nutrient loading to Grand Bay is relatively small, with ambient nutrient concentrations often below detection. One notable event introduced high phosphate levels in April 2005 as a result of a breach in a containment levee from a gypsum stack near Bangs Lake (Dillon and Walters 2007). The open-water areas support widespread submerged aquatic vegetation, including widgeon grass (Ruppia maritima) with smaller patches of shoal grass (Halodule wrightii ). The muddy intertidal areas support scattered, unconsolidated, or fringe oyster reefs. The shoreline includes extensive freshwater, brackish and salt marsh habitat and unvegetated salt flats or pannes.

Data Sources
For this study, we acquired water quality, nutrient, chlorophyll a , and weather data from NOAA's National Estuarine Research Reserve Centralized Data Management Office (CDMO: http://cdmo.baruch.sc.edu/). The water quality data included high frequency (15 or 30 min) records from YSI water quality multimeters deployed at two to four sites per estuary (total of nine sites, see Table 2). Sondes were deployed 30 cm from the bottom. Measurements included temperature, salinity, and dissolved oxygen (DO); water depth, pH, and turbidity. For the metabolism analysis, temperature, salinity, and dissolved oxygen data were used. The time span of data collection ranged from 5.4 to 15.8 years ( Table 2). The three longest records (>15 years) included two sites at Weeks Bay (WB and FR) and one site at Apalachicola Bay (EB).
Weather data were collected at 15-min intervals since 2002 at most sites (Table 2) were processed and then archived at CDMO. Measurements included air temperature, relative humidity, barometric pressure, precipitation, photosynthetically active radiation (PAR), wind speed, and wind direction. Because the period of record for weather data did not completely overlap the water quality data, we summarized key weather parameters as climatological monthly means by hour. The wind summaries were used to estimate the gas exchange at the air-sea interface.
Dissolved inorganic nutrients and chlorophyll samples were collected at monthly intervals at most sites beginning Per estuarine area, based on Southeast regional SPARROW model (Hoos and McMahon 2009) in 2002 (Table 2). Flow data for the Apalachicola and Fish Rivers over the study period were obtained from the USGS website (http://water.usgs.gov/data/explorer/). The USGS SPARROW model for the southeastern USA (Hoos and McMahon 2009) was used with the SPARROW Decision Support System to determine total nutrient loading to each estuary (http://cida.usgs.gov/sparrow/).

Metabolism Analysis: Theory and Practice
We employed the "open-water method" pioneered by Odum (1956) using DO time series data to infer ecosystem metabolism parameters described by general mass balance equation where C is DO concentration (mmol m −3 ), t is time (h), P is the photosynthesis rate (mmol m −3 h −1 ), and R is the respiration rate (mmol O 2 m −3 h −1 ), and D is the rate of oxygen exchange across the air-water interface (mmol O 2 m −3 h −1 ). Gas exchange at the air-water interface was calculated as: where k a is the volumetric reaeration coefficient (h −1 ), C s is the DO saturation concentration (mmol O 2 m −3 h −1 ), calculated as a function of water temperature and salinity (Benson and Krause 1984). The coefficient, k a , is itself dependent on wind speed, barometric pressure, air temperature, water temperature, salinity, and depth of the water column at the sampling site. For k a , we used a modified form of unified equation developed by Ro and Hunt (2006), as implemented by Thebault et al. (2008).
where H is mean water depth (m) at the sampling site, D w is diffusivity of oxygen in seawater (m 2 s −1 ), ν w is the kinematic viscosity of seawater (m 2 s −1 ), ρ a is the density of air (kg m −3 ), ρ w is the density of seawater (kg m −3 ), and U 10 is the wind speed normalized to 10 m height above ground level (m s −1 ). Further details of the calculations are well documented in Thebault et al. (2008). Because weather data were only available for a portion of the study period (Table 2), we generated monthly mean wind speeds for each hour of the day for k a calculations (i.e., Fig 2). The diffusion corrected DO fluxes (dC /dt −D ) were averaged separately during day and night periods to compute hourly rates of apparent primary production (P) and nighttime respiration (R ), respectively (units: mmol O 2 m −3 h −1 ). Respiration rates were assumed constant during day and night, thus daily R was calculated as R ×24, and daily P was calculated as (P−R )×day length (h). These volumetric daily rates of production and respiration were multiplied by water depth at each site to yield areal ecosystem gross production (GPP) and ER. Finally, NEM was calculated as: The open water method has been applied in a wide variety of aquatic ecosystems, though it depends on several key assumptions that need to be considered in interpreting the Included for each site are the number of data records, the number of years covered in the data record and the date when monitoring was initiated results. One of the most important assumptions is that the water mass being monitored is reasonably homogenous and well mixed, essentially having common metabolic history. Thus, at sites where physical mixing dominates over biological processes, metabolic rates may be less reliable (Kemp and Boynton 1980, Caffrey 2003, Collins et al. 2013. A similar assumption is that, in shallow systems, the water column is well mixed, thus in frequent contact with both the sedimentwater and air-water interfaces. If true then all submersed primary producers contribute to the observed dissolved oxygen fluxes including phytoplankton, seagrass, and benthic macro-and microalgae. Similarly, both water column and benthic heterotrophs contribute to the observed respiration. Diurnal oscillations of oxygen concentrations are commonly observed in aquatic environments, increasing during daylight hours as photosynthesis exceeds respiration and declining during night when respiration exceeds photosynthesis. Inevitably, there are periods when this regular diurnal pattern is not evident, which may point to various other causative factors. First, the lack of strong diurnal oxygen cycles may occur during periods of low metabolism, such as relatively oligotrophic systems, or seasonally, during dormant winter periods. Second, anomalous diurnal oxygen cycles can also be caused by advection of water masses past the sensor that have a distinct metabolic history (Odum and Hoskins 1958;Kemp and Boynton 1980;Caffrey 2003) or by certain wind conditions that increase the air-sea exchange (Ziegler and Benner 1998). When the open water calculations yield anomalous results-either negative GPP or positive ERthose values are frequently discarded (Caffrey 2003, Collins et al. 2013), though few investigators specify how anomalous data were handled. One possibility not often considered is that Hour Wind Speed (m s -1 ) Fig. 2 Hourly climatological monthly average wind speed for January, April, July, and October. a Grand Bay. b Weeks Bay. c Apalachicola Bay anomalous values may simply reflect periods of low metabolism, below the effective method detection limit. Thus, discarding such values would upwardly bias estimates of metabolism.
We evaluated some of these concerns in a bit more detail using a subset of the time series from Weeks Bay, reasoning that the four Weeks Bay sites tend to span the estuarine gradient from freshwater (FR) to mesohaline (WB), so should represent a diversity of conditions. The analysis consisted of comparing the daily variability in salinity, sigma-t, and water depth with calculated metabolism values. We reasoned that daily variability in water quality parameters should provide a relative measure of the degree of advection occurring at each site. We anticipated that days with high variability in either salinity, sigma-t or water depth would be more likely to coincide with anomalous metabolic rates. We also looked at the change in PAR between days, anticipating that days with low PAR would be more likely to have anomalous metabolic rates.

Statistics
Relationships among the variables were examined with crosscorrelation analysis (CCA) using SPSS-PASW Version 18. This analysis is appropriate for time series data which tend to be autocorrelated; hence, values from adjacent time steps are not independent. CCA is analogous to Pearson product correlation analysis, allowing one to describe the relationship between two time series (r values). In addition, lags or offsets between the two time series can also be determined. The reported correlation coefficients have no time lags, unless specifically noted. The traditional P value denoting significance was not reported, however, upper and lower 95 % confidence intervals were used to evaluate the significance of the relationships between the two time series. Temporal trends in time series were evaluated using nonparametric Kendall Tau analysis, allowing one to describe interannual or seasonal trends (Jassby and Powell 1990). The Kendall Tau analysis was conducted only for the longest time series from each estuary (Apalachicola Bay EB, Weeks Bay FR, Weeks Bay WB, and Grand Bay BL).

Environmental Conditions
Wind conditions at all reserves exhibited a classic diurnal sea breeze pattern with higher winds during the day and lower winds at night (Fig. 2). Overall, wind speeds tended to be highest at Grand Bay (mean, 3.9 m s −1 ), followed by Apalachicola Bay (mean, 2.8 m s −1 ) and Weeks Bay (mean, 2.0 m s −1 ); however, these systematic differences may in part reflect differences in the immediate surroundings near each weather station. Seasonal water temperature patterns were virtually identical in the three estuaries, ranging from 10 to 12°C in winter to over 30°C during the summer (Fig. 3a). In contrast, there were distinctive patterns in riverine discharge (Fig. 3b) and salinity (Fig. 3c) across estuaries. In Grand Bay, monthly average salinities ranged from 10 to 30 and were less seasonally variable than either Apalachicola Bay or Weeks Bay (Fig. 3c). In Weeks Bay, monthly average salinity was lower than the other two estuaries, rarely exceeding 20. Weeks and Apalachicola Bays exhibited consistent seasonality with lowest salinity in spring and highest during the fall, reflecting the prevailing freshwater flow regime (Fig. 3b). Monthly average flow from the Apalachicola and Fish Rivers were negatively correlated with monthly average salinity in the respective estuaries (r =−0.71, Apalachicola Bay; r =−0.44, Weeks Bay). Salinity among the different sites within each estuary was strongly coherent, with cross-correlation coefficients exceeding 0.93 with no time lag among sites within Weeks and Grand Bays. In Apalachicola Bay, salinities at the EB and CP sites were strongly cross-correlated with no time lag (r =0.91), while salinity at the DB site lagged 1 month behind EB and correlation between the sites was less strong (r =0.57).
These three estuaries had sharply contrasting DIN and chlorophyll a concentrations, while DIP concentrations were low, often below detection limits in all estuaries (Table 3). Grand Bay had the lowest DIN concentrations, followed by Apalachicola Bay and Weeks Bay, respectively (Table 3). Surprisingly, in Apalachicola Bay, even though the EB site is fresher than the open bay sites, DIN at this site was often lower than either CP or DB. This site also had the highest chlorophyll a concentrations, about double those of CP (Table 3, Fig. 4). Within Weeks Bay, DIN was highest (typically >30 μM) at the river sites (FR and MR), and somewhat lower at the open bay sites (MB and WB) ( Table 3).
Chlorophyll a concentrations among the three estuaries reflected the overall nutrient regime, being highest in Weeks Bay, intermediate in Apalachicola Bay, and lowest in Grand Bay (Table 3, Fig. 4). Chlorophyll a was often higher during summer months at both sites in Grand Bay and Apalachicola Bay CP. Chlorophyll a at DB and EB sites in Apalachicola Bay peaked in May. Distinct seasonal patterns in chlorophyll a were not observed in Weeks Bay (Fig. 4). While median turbidities were similar at all sites (Table 4), peak turbidity was highest in Weeks Bay. Turbidity and salinity were negatively correlated at Grand Bay PC (r =−0.55), Weeks Bay MR (r =−0.23), and Weeks Bay WB (r =−0.26).

Metabolism
The time series of ecosystem GPP and ER showed distinct seasonality at most sites, with high rates during summer and low rates during winter (Fig. 5). In Grand Bay, both GPP and ER were the lowest of the three estuaries with peak summer rates ranging between 130 and 160 mmol O 2 m −2 day −1 for GPP and between −150 and −170 mmol O 2 m −2 day −1 for ER (Fig. 5a). Apalachicola was intermediate with summer peak rates ranging from 200 to 250 mmol O 2 m −2 day −1 for GPP and from −220 to −300 mmol O 2 m −2 day −1 for ER (Fig. 5c). Weeks Bay had the highest metabolic rates of the three estuaries with summer rates ranging from 200 to 400 mmol O 2 m −2 day −1 for GPP and from −300 to −450 mmol O 2 m −2 day −1 for ER (Fig 5b). Apalachicola and Grand Bays showed strong seasonal cycling, whereas Weeks Bay had relatively weak seasonal patterns, especially in GPP (Fig. 6).
Rates among individual sites within each estuary were often very similar (Figs. 5 and 6). One exception was the CP site in Apalachicola Bay, which was consistently about 30 % lower in GPP than the other sites. In Weeks Bay, two sites (FR and MR) showed no distinct seasonal pattern in GPP, whereas the other Weeks Bay sites (MB and WB) had more noticeable seasonal pattern; the cross-correlation coefficient between temperature and gross production at WB was high (r =0.77), reflecting the distinct seasonality. Similarly in Grand Bay, GPP and temperature were strongly correlated (r =0.68), but was weaker at EB in Apalachicola Bay where GPP and temperature were not significantly correlated (r =0. 10). The variation in metabolic rates across the range of temperatures is shown in Fig. 7. The curved lines represent the expected relationship, assuming a nominal Q 10 of 2, fit to the mean value for each estuary. Thus, the relative distribution of points Temperature and salinity data were averaged across the NERR sites used in this study. Flow data from Fish River near Silver Hill (site: 02378500) and from the Apalachicola River near Sumatra (site: 02359170) were acquired from USGS NWIS data portal (http://waterdata.usgs.gov/nwis). There was no gauged river flow for Grand Bay Median (5th and 95th percentile) values over the period of record around the curves provides a visual way to evaluate the fit with a nominal temperature curve. For ER, all three sites show a similar trend with temperature, with the Q 10 curves appearing to bisect the data over the entire temperature range. For GPP, there was a similar temperature trend, especially for Grand Bay and Apalachicola Bay. For Weeks Bay, GPP values fell progressively further below the Q 10 curve as temperature increased (Fig. 7b). When viewed in relation to NEM (Fig. 8), the relatively subtle temperature trends in the component measures became more pronounced. In all three estuaries, there was a significant negative relationship between NEM and temperature with systems becoming more heterotrophic as temperature increased (r =−0.5).
The relationships among salinity, DIN concentrations, and gross production were generally weak. While salinity and gross production were significantly correlated in all three estuaries (data not shown), the cross-correlation coefficient ranged between −0.3 at Grand Bay (gross production was high when salinity was low) to +0.3 at Apalachicola Bay and + 0.5 at Weeks Bay (gross production was highest when salinity was high). In addition, the gross production time series lagged the salinity time series by 2-3 months depending on the estuary, indicating a delay in the response of gross production to changes in salinity. The relationship between DIN concentrations and gross production was either non-significant (Grand Bay) or relatively weak (r =−0.24 with 1 month lag, Apalachicola Bay; r = −0.46 with no lag, Weeks Bay). Although weak, this negative relationship suggested that primary producers in Apalachicola and Weeks Bays effectively removed DIN from the water column.
Respiration rates also showed distinct seasonality (Fig. 7), with cross-correlation coefficients between temperature and respiration greater than 0.7 at Weeks Bay WB and Grand Bay BL. As with gross production, respiration rates were highest at Weeks Bay and lowest in Grand Bay. In all three estuaries, ER and GPP rates were strongly correlated (r >0.80). Within each estuary, lower salinity sites such as Weeks Bay FR and Apalachicola Bay EB had higher respiration than at the higher salinity sites. All sites were net heterotrophic for much of the year, most strongly so in the summer (Fig. 6), when temperatures were high (Fig. 8). Sites were generally net autotrophic or near balance between November and February.

Interannual and interestuary comparisons
The temperature, salinity, and metabolism time series ranged from 5 to over 15 years at these sites. The record included both dry (1999, 2000, and 2007) and wet (1998, 2003, 2005, and 2010) years (Fig. 3). In Weeks Bay, there was a significant increase in salinity over time of about 3 % year −1 (p=0.02) at both FR and WB sites, while respiration declined by 0.8 % year −1 (p =0.04) at the WB site, but not the FR site. Gross production was highest during 2007 at FR, MB, and MR sites. While gross production at Apalachicola Bay EB was generally highest when salinity was above 8, there were no significant interannual trends in salinity or gross production in Apalachicola or Grand Bays, nor were there consistent differences in metabolism between drought and wet years in these two estuaries. Primary production and respiration were strongly coupled (r >0.8, Fig. 9). Annual mean metabolism indicated that all three estuaries were net heterotrophic for all years with the sole exception of 1998 in Apalachicola Bay. The monitoring sites in Weeks Bay had the highest means and ranges in DIN and chlorophyll a concentrations (Table 3); DIN concentrations were greater in wet years compared to drought years. In contrast, Apalachicola Bay showed no clear interannual patterns in DIN and chlorophyll a concentrations between wet and dry years. Similarly, Grand Bay had minimal interannual variability in salinity, DIN, or chlorophyll a.

Assumptions of open water method
The use of high frequency dissolved oxygen data for estimating ecosystem metabolism parameters has seen wide use in a variety of freshwater, estuarine, and marine systems, but it is important to evaluate the underlying assumptions and limitations of the method. The choice of a suitable sampling site is Primary production converted to C units assuming PQ=1. Apportionment into water column and benthic fractions (in bold) assuming an average of 27 % benthic production, estimated from the literature from similarly shallow non seagrass systems a Includes seagrass production perhaps the most important factor to consider. For example, locating a site near the center of well-mixed embayments might be desirable because it is expected to be less influenced by tidal exchange from nearby wetlands or by upland runoff. While a data sonde deployment is typically at a fixed location, it integrates processes occurring over spatial and temporal scales related to tidal excursion (D'Avanzo et al. 1996, Collins et al. 2013. Thus, estuarine metabolism measurements from between two and four sites per estuary may differ from than studies with more extensive spatial, but more limited temporal coverage. In deeper systems with persistent stratification, deployment of surface and bottom sondes may be necessary to more accurately estimate total system metabolism, as well as delineation of the depth of surface and bottom layers. In this study, all sites were very shallow (<2 m) and were rarely (GND PC, GND BL, APA CP, and WKB WB) or intermittently (APA EB APA DB, WKB FR, WKB MB, and WKB MR) stratified (L. Levi, personal communication, K. Cressman, personam communication, and S. Phipps, unpublished data). A previous study examined paired sonde deployments in Apalachicola Bay (EB) finding that calculations of GPP and ER were not significantly different between surface and bottom layers (Caffrey 2003). The percentages of anomalous metabolism values (i.e. −GPP or +ER) ranged from 5 % at Grand Bay PC to 24 % at Weeks Bay FR and were similar to previous analyses at these sites (Caffrey 2004) and similar to that found in other studies (Collins et al. 2013). We expected that periods of high advection at a given site might account for most of the anomalous values, and that high intraday variability in salinity, sigma-t, or water  Fig. 5 Time series of monthly mean gross production (GPP) and total respiration (ER) in a Grand Bay, b Weeks Bay, and c Apalachicola Bay. The symbols denote different sites within each estuary. The lines represent the average across sites depth would be a good predictor of advection. However, this was not the case; values with high daily variability in these physical characteristics were not more likely to be anomalous than those with low daily variability (data not shown). We also examined several weeklong data records during which large rain events in the watershed led to rapid and extreme changes in salinity and PAR. However, these weather events did not appear to be related to instances of anomalous metabolic rates. Collins et al. (2013) used a signal processing approach to transform the dissolved oxygen that eliminated both short term (<2.5 h) and longer term (>25 h) variability. While this approach may have some advantages and holds promise for future studies, these techniques did not eliminate anomalous metabolism values, In fact, Collins et al. (2013) reported a similar percentage of anomalous values in their analysis to what we observed here. It seems likely that some proportion of the anomalous values simply reflect cases where the values fall beneath an effective detection limit for this method, thus may be more or less randomly distributed, therefore confounding a simple explanation or means of handling those data.
Wind mixing leading to air-sea exchange can be an important component of the metabolism calculations, and is especially important in these shallow systems. A variety of studies have examined different methods for parameterizing air-sea exchange of dissolved oxygen (Hartman and Hammond 1984, Marino and Howarth 1993, D'Avanzo et al. 1996, Ziegler and Benner 1998, Ro and Hunt 2006, Thebault et al. 2008, Collins et al. 2013, although a systematic comparison of how each method influences metabolic rates is yet to be done. The use of wind speed derived estimates of the reaeration coefficient (k a ) provides a more realistic model of the air-sea exchange compared to using a constant k a . However, many existing data records may not have wind data available, or the location of the wind sensor may not be optimal, For example, relocation of the Weeks Bay weather station to an interior location after damage by Hurricanes Ivan, Dennis, and Katrina, resulted in lower average wind speeds. The frequency of winds >3 m s −1 occurred 30 % of the time before relocation and<5 % of the time after relocation (S. Phipps, unpublished data). Thus, colocation of wind sensors and sondes is ideal for this analysis. Because we did not have complete wind records at these sites, we generated hourly climatological values for each month of the year using available data. However, this approach cannot accurately capture small scale variations from normal wind patterns, be they either very windy or calm, which might have contributed to some of the anomalous values. In addition, currents may also influence air-sea exchange (Hartman and Hammond 1984). Thus, more detailed physical modeling of the systems would help to constrain the relative importance of advective transport and air sea exchange and lead to more accurate and precise estimates of ecosystem metabolic parameters.

Seasonal Response
The seasonal ranges in temperature, salinity, and nutrients in these three estuaries were similar to those reported in other Gulf of Mexico estuaries (Pennock et al. 1999, Murrell et al. 2007). Despite the smaller temperature span relative to temperate systems (∼10-30°C vs 0-30°C, respectively), both gross primary production and respiration showed strong seasonality. However, the temperature effect was stronger for respiration than for primary production, leading to more heterotrophic conditions during summer (Figs. 6 and 8), similar to numerous other studies (Caffrey 2004, Yvon-Durocher et al. 2010Vaquer-Sunyer et al. 2012;Yvon-Durocher et al. 2012). We used CCA to explore how environmental factors such as salinity and nutrient inputs covary with the peaks in  Fig. 7 Relationship between temperature and monthly mean gross production (GPP) and total respiration (ER) for each site in: a Grand Bay, b Weeks Bay, and c) Apalachicola Bay. The symbols denote the different sites within each estuary. The curved lines represent the nominal expected temperature response assuming Q 10 =2 fitted through the mean of all values. The dotted line bisecting each panel is the zero line summer metabolism. While solar insolation is maximal during the summer, the light available for both planktonic and benthic primary production depends on water column light attenuation. While data are limited, long term records of Secchi disk depth from Apalachicola Bay show little seasonality (L. Levi, personal communication), implying relatively uniform light attenuation. There are far fewer light attenuation measurements for Grand or Weeks Bays, so it is difficult to make a full evaluation. However, a recent study reported that K d ranged from 0.5 to 3.2 m −1 and used these values in a resource modeling exercise to conclude that light limitation was less important than nutrient limitation in controlling phytoplankton productivity (Amacker 2013).
The timing and magnitude of nutrient loading to coastal waters is a key factor controlling primary production in systems as varied as Chesapeake Bay (Kemp et al. 2005), the Baltic Sea (Rydberg et al. 2006), Neuse River (Mallin et al. 1993, Boyer et al. 1993, Moreton Bay (Eyre and McKee 2002), Pensacola Bay (Murrell et al. 2007), Elkhorn Slough (Caffrey et al. 2007), and Long Island Sound (Collins et al. 2013). Often, there is a lag between delivery of nutrients and the response by primary producers (Boyer et al. 1993, Caffrey et al. 2007, Murrell et al. 2007). In addition, the location of the peak response of primary producers to nutrient pulses varies within and among estuaries based on the interaction between mixing characteristics and residence times. This pattern of enhanced phytoplankton production related to nitrogen inputs has been reported in Apalachicola Bay, with strong seasonal variation in the timing and mode of response (Mortazavi et al. 2000b). While we have no direct measure of monthly nutrient inputs to these systems, we examined two proxies: salinity, which reflects the freshwater inputs, and DIN concentrations. Of the two, salinity provides more coherent patterns because DIN concentrations within the estuary are influenced by  Fig. 8 Relationship between temperature and monthly mean net ecosystem metabolism (NEM) for each site in: a Grand Bay, b Weeks Bay, and c Apalachicola Bay. The symbols denote the different sites within each estuary. The zero line bisecting each panel denotes balance uptake and removal processes within the estuary. In fact, in both Apalachicola Bay and Weeks Bay, DIN was negatively correlated with gross production, consistent with uptake by primary producers influencing observed DIN distributions. However, in Apalachicola Bay and Weeks Bay, reduced salinity was followed by declines in gross production rather than increases. This may have been the result of reduction in light availability due to higher turbidity, and low residence times resulting in phytoplankton washout during low salinity events. Only in Grand Bay were periods of low salinity followed by enhanced primary production.

Interannual trends
One somewhat surprising result was the relatively small range in annual rates of primary production (Fig. 9) despite the large variation in DIN and salinity (Fig. 3). However, in Weeks Bay at the WB site, there was a significant long-term increase in salinity, lower DIN concentrations during dry years and a concomitant decrease in metabolic rates. This has been observed in other estuarine and coastal systems time series (Lindahl et al. 1998;Rydberg et al. 2006). So why did we see an interannual relationship between freshwater flow, nutrient input and primary production in Weeks Bay, but not Apalachicola Bay or Grand Bays? Previous studies in Apalachicola Bay have shown a relationship between riverine nitrogen loading, chlorophyll a and phytoplankton production (Mortazavi et al. 2000b).
One major difference between this analysis and other studies is that most long term measurements of metabolism have focused on measurements of phytoplankton production (Rydberg et al. 2006;Steinberg et al. 2001) as opposed to whole system metabolism as in this study. To the extent that the assumptions of open water method are met, production and respiration estimates include both water column and benthic processes. While benthic respiration in these shallow systems can be a significant contributor to total system respiration (Murrell et al. 2009;Mortazavi et al. 2012), we focus here on the role of benthic production.
Many northern Gulf of Mexico estuaries are so shallow that the euphotic zone extends throughout the water column to the sediment water interface over relatively large portions of the estuaries (Pennock et al. 1999;Turner 2001). Recent measurements of light attenuation from all three systems indicated that the benthos receives about 5 % of surface irradiance in Weeks Bay and Apalachicola Bay and about 15 % in Grand Bay (based on K d values from Amacker 2013), similar to results from 10 years of Secchi disk measurements in Apalachicola Bay (L. Levi, personal communication). Thus, benthic primary production is likely a significant component of total production in these systems. In Gulf of Mexico estuaries, it can represent up to 32 %, of total primary production (Table 4; Schreiber and Pennock 1995, Moncreiff et al. 1992, Murrell et al. 2009). This percentage is similar to observations from other tropical and temperate estuaries (Borum and Sand-Jensen 1996; Table 4). In addition, benthic chlorophyll a can represent a significant component of the total stock of chlorophyll a (Cahoon 1999).
In order to evaluate the potential role of benthic productivity, we assumed that benthic production was an average of 27 % of total system production, based on literature values from non seagrass systems (Table 4). Based on this back of the envelope calculation, estimated benthic production in Weeks Bay was 226 g C m −2 year −1 , about double that in Apalachicola or Grand Bays (Table 4). These values are in the same range as those measured in other estuaries (Schreiber and Pennock 1995, Borum and Sand-Jensen 1996, Cahoon 1999. It is unclear what the seasonal or interannual variability in benthic production is in these three estuaries. Benthic primary producers may be less sensitive to riverine nutrient inputs and more dependent on benthic regeneration of nutrients (McGlathery et al. 2001). Is phytoplankton production stimulated by nutrient inputs during wet years, while benthic production is stimulated during drought years due to greater water clarity? Seasonal and interannual variations in light attenuation likely drive benthic primary production, however, more focused studies are needed to address these questions.

Interestuary Comparisons
The metabolism results from this study are consistent with those measured in other systems: Weeks Bay falls in the upper end of the range and Grand and Apalachicola Bays fall more towards the mid range of the literature values (Table 4). The estimated relative contribution of water column and benthic production at our sites was similar to that found in temperate estuaries (Boynton et al. 1982;Borum and Sand-Jensen 1996). Estimated water column production in Apalachicola Bay is similar to prior measurements of phytoplankton  Fig. 9 Relationship between annual mean GPP and ER across the three estuaries studied. The symbols denote the average value of all sites within each estuary production (Mortazavi et al. 2000b), while production in Weeks Bay was somewhat higher than prior measurements (Schreiber and Pennock 1995). The strong relationship between production and respiration at all sites (r >0.8) suggests that respiration is strongly linked to autochthonous sources of organic matter (Fig. 9). However, all sites were net heterotrophic on average indicating that allochthonous sources of organic matter also contribute to system metabolism ( Fig. 9). Interestingly, the magnitude of net heterotrophy was similar among the estuaries despite differences in watershed area and land cover between these systems. In this study, the watershed nutrient inputs appeared to strongly influence the rates of autotrophic biomass and production For example, the total nitrogen loading from the Weeks Bay watershed was >140,400 kg N km −2 year −1 , twice that of Apalachicola Bay (Table 1, Fig. 10). Consequently, it is not surprising that production and chlorophyll a in Weeks Bay was the highest of all three estuaries. Literature relationships between nitrogen loading and primary production in Weeks Bay and Apalachicola Bay were very similar to predicted relationships, whereas Grand Bay had much higher production than predicted from nitrogen load ( Fig. 10; Nixon et al. 1996;Conley et al. 2000). This suggests that Grand Bay may receive other nitrogen inputs not accounted for in the SPARROW model. Potential sources include diffuse groundwater inputs or from nearby Mobile Bay via Mississippi Sound (Pennock et al. 1999).

Management Implications
The use of long-term SWMP datasets to estimate metabolic rates is a valuable tool, especially because they are transferrable across estuaries within the NERR system due to the consistent monitoring methodology. Thus, estimates of metabolism can be compared in a meaningful way to evaluate the magnitude or direction of changes over time. Reserves in the northern Gulf of Mexico are particularly susceptible to climate change impacts due to their latitude, geomorphology and in the case of Apalachicola Bay, a large watershed. These shallow systems will be impacted by increasing sea level rise, which result in increasing salinity, increasing inundation of near shore areas, and most likely, more intense storm effects. Because of the concomitant nature of these impacts, it may be difficult to identify simple cause and effect relationships between climate change measures and trophic state of estuarine and coastal systems, but system metabolic parameters should be a sensitive indicator of changes in system function. While the resiliency of an estuary is difficult to measure because of its dynamic nature, this tool provides a potentially useful metric in understanding the response of estuaries to changing climate, and other impacts, whether natural or anthropogenic.
Overarching management goals for these estuaries are to reduce alterations to hydrology, to limit nutrient inputs and to increase public understanding of the linkages between anthropogenic impacts and critical estuarine habitat functions. In Weeks Bay, non point source pollution and harmful algal blooms are the dominant management concerns. An increase in urban and agricultural development in the watershed since the 1990s has led to increased impervious surface area causing more rapid and flashy runoff, an overall increase in residence time, and higher nutrient loading in the watershed (Basnyat et al. 1999). Large and prolonged blooms of dinoflagellates occurred during the drought in 2007, perhaps a result of increased residence time in the system. This was also a period of high gross production in the upper reaches of the Bay. Management goals for Weeks Bay are to preserve and restore natural shoreline, riparian zones and wetlands with emergent plant assemblages, but implementation is politically limited. In contrast, a critical management issue for Apalachicola Bay is the reduction in freshwater flow and nutrient inputs due to upstream diversion. Maintenance of nutrient inputs is particularly important during the summer when phytoplankton production peaks (Mortazavi et al. 2000b;Putland 2005) and when phytoplankton are most likely to be nutrient limited (Mortazavi et al. 2000b;Amacker 2013). Thus, adequate freshwater flow to the estuary is essential for phytoplankton production and the oyster fishery and higher trophic levels which are dependent on phytoplankton production (Chanton and Lewis 2002). While Grand Bay is a relatively pristine system, phosphorous releases from the gypsum spill in 2005 and subsequent recent episodes (K. Cressman, personal communication) are a serious management concern. Primary production in Grand Bay  (Nixon et al. 1996, solid line), and from Danish estuaries (Conley et al. 2000, dashed line). The means from each of the three northern Gulf of Mexico estuaries in this study are shown as filled circles, which are identified by the labels at the bottom of the panel. The error bars represent the range in means among the sites. The label, Nitrogen loading, refers to total nitrogen (TN), except for Nixon et al. (1996), which is DIN following the 2005 spill was similar that of other years which is consistent with bioassay results showing that nitrogen rather than phosphorus stimulates phytoplankton growth in this estuary (Amacker 2013). Thus, these systems each have unique management concerns, but each also benefits from sustained monitoring programs that provide place-based scientific information. Such information is critical in working with private landowners and to the general public through education and coastal training programs as well as to provide information needed by state regulatory agencies.
Long-term estimates of primary production and ecosystem respiration are rare in the estuarine literature, yet they provide fundamental information about the trophic status of estuarine environments. The traditional experimental methods of measuring primary production and respiration are relatively labor intensive, thus there are few regions with sufficient resources to maintain such a measurement program for more than a few years. The appeal of the open water method is that it uses more readily available and cost-effective sonde monitoring technology to estimate these fundamental estuarine rate processes. When collected consistently over a long period of time, this approach makes it possible to look for long-term trends, but perhaps more importantly, it provides a historical benchmark for evaluating future patterns. Additionally, it may be particularly valuable in providing context for more focused short term studies.
The role of estuaries and coastal oceans in global carbon budgets continues to be a topic of active research and debate (Smith and Hollibaugh 1993;Gattuso et al. 1998;Cai 2011), yet a clear understanding of the factors that affect the carbon balance is crucial for making informed management decisions. Clearly the net metabolic balance of estuaries and the coastal zone are sensitive to temperature and eutrophication, but will also be affected by emerging environmental changes such as sea level rise and ocean acidification. Thus, maintenance of active monitoring networks is critical to evaluate the role of human activities in modulating these processes.