NERRS banner

System-wide Monitoring Program

Synthesis of the Water Quality Data
METHODS

Data Collection and Management
Collecting long-term trend data that captures natural variability requires adequate temporal and spatial coverage. To address the temporal coverage, the National Estuarine Research Reserve System-Wide Monitoring Project employs Yellow Springs Instrument Co. (YSITM) model 6000 or 6600 UPG data sondes to collect water quality data. These data loggers record at 30-minute intervals, relay measurements to internal memory, can run unattended for weeks at a time, and operate in depths from a few cm to greater than 150 m. The attached dissolved oxygen, conductivity, turbidity, temperature and pH sensors can be quickly replaced in the field, if necessary. The dissolved oxygen sensor provides accurate readings without the use of an auxiliary stirrer and with little drift for extended periods, although deployment duration may differ among sites based on differences in fouling and water flow. The turbidity probe features a mechanical cleaner for the optical face, thereby preventing fouling in long-term deployments. The data logger interfaces with a PC or laptop, as well as to real-time data collection platforms using telemetry systems. Built-in software exports data to spreadsheet programs and generates basic statistics and plots.

For spatial coverage, each Reserve deploys a minimum of two data loggers. At half of the Reserves, one data sonde is deployed at a reference location and serves as a long-term control and additional loggers are deployed at other locations to monitor conditions related to specific non-point source concerns (Table 1) within the Reserve (Trueblood et al. 1996). With the exception of two Reserves (North Carolina and Waquoit Bay) where the data sondes are deployed at sites with minor anthropogenic disturbances, data sondes at the remaining Reserves are deployed at impacted sites only. In addition to the data sondes, many NERR sites deploy additional water sampling instruments and measure other water quality variables, such as nutrients. Standardized protocols developed by the Reserves assure that sampling, processing and data management techniques are comparable among sites. A Centralized Data Management Office (CDMO) at the Belle W. Baruch Institute for Marine Biology and Coastal Research of the University of South Carolina provides additional quality control for data and metadata.

Table 1. Types of non-point source pollution categories addressed at different National Estuarine Research Reserves (NERR). Blank spaces indicate non-issues for the respective Reserve.

NERR
Urban

 Agricultural

 Boating
Wetland
Restoration
ACE Basin
X
Apalachicola
X
Chesap. Bay MD
X
Chesap. Bay VA
X
Delaware Bay
X
X
Elkhorn Slough
X
Great Bay
X
X
Hudson Bay
X
Jacques Cousteau
Jobos Bay
X
Narragansett Bay
X
North Inlet
North Carolina
Old Woman Creek
X
Padilla Bay
X
Rookery Bay
X
Sapelo Island
X
X
South Slough
X
Tijuana River
X
Waquoit Bay
X
Weeks Bay
X
Wells
X

Site Summaries and ancillary data
In an effort to characterize sites at each Reserve where water quality data were collected, a standardized form requesting information on sites was sent to each Research Coordinator. To supplement written descriptions of each site, maps showing location of water quality sampling sites were obtained from the NOAA Estuarine Reserves Division. Information requested included:
  • Latitude and longitude
  • Salinity range and mean for site
  • Tidal frequency, range and mean for site area
  • Creek or water body dimensions (length, mean depth, mean width)
  • Distance meter is above the bottom
  • Creek bottom habitat at site
  • Dominant marsh plants and submerged aquatic vegetation near site
  • Dominant upland vegetation
  • Upland land use
  • Activities potentially impacting the site

Data Review and Protocols for Deletion
Data from each of the 22 NERR sites were downloaded from the CDMO FTP site. Data were then imported into MS Access to obtain a continuous time series for all three years. Queries were run to detect duplicate data and verify the number of records for each day. There were 52,608 possible records for each site within a reserve. As data were reviewed, it became evident that a complete QA/QC would be necessary on all data received. Scatter plots of yearly data were graphed for each of the five key variables at each site. If graphs indicated that data were suspect or unusual from the graphs, then the Anomalous Data section of the metadata was consulted to determine if the Reserve also felt the data were suspect. If the Reserve considered the data suspect, then those data were deleted. If data were suspect or atypical and the Reserve did not make note of suspect data in the Anomalous Data section of the metadata, a second opinion was sought to determine if the data were unusual enough to delete. In both instances, the whole deployment was deleted unless an event occurred that subsequently altered records after the event, in which case data were deleted from the onset of this event until the end of the deployment.

Graphs were particularly useful for detecting erroneous transposing of data. Erroneous transposition was obvious when data were lower or higher than what is usually found at that site. By adding deployment and retrieval times to the data and plots, we were able to detect data that were not removed at the beginning and end of a deployment, even when the meter was not in the water. Several criteria were used in review of the data. If data for at least two variables (usually salinity and depth) indicated the probes were exposed or partially exposed due to a meter change or low tide, then all data were deleted for that time period. If the temperature data were determined to be suspect and were deleted, then corresponding records of specific conductivity, salinity, and DO (mg/l) were also deleted (McDonald 1996). Similarly, if data for specific conductivity and salinity were determined to be suspect and were deleted, then records for DO (mg/l) were also deleted.

The data were also reviewed for occurrence of negative numbers. If depth and temperature data were negative, they were retained. Because salinity and DO data cannot be negative, those records were changed to zero. In some cases, the DO data were negative for long periods of time, suggesting that the meter was not functioning properly. Those data were evaluated to determine if deletion was appropriate. Notation was made in the Access database for all cases in which data were deleted. All deletions were annotated and compiled in electronic files and sent to the Reserves and the CDMO.

Data Analysis and Synthesis

Sampling interval
A two-sample t-test was used to compare daily mean, minimum and maximum DO (% saturation) at 30-min and 4-h sub-sampled intervals during deployments in July-August 1997 and 1998 at Big Bay Creek and St. Pierre Creek in the ACE Basin NERR to determine the effects of sampling interval on accuracy of estimates of mean, minimum, and maximum values.

Effect of deployment period on dissolved oxygen (% saturation)
Potential drift in dissolved oxygen (% sat) due to fouling necessitated determination of whether noticeable differences in hypoxia and supersaturation occurred over an entire deployment. Hypoxia and supersaturation, expressed as percent of time per deployment, were evaluated at 1, 2, 4, 7 and 14-day intervals for July-August 1997 and 1998 data. Line graphs were created in MS Excel to evaluate decay/increase in DO (% saturation) readings with respect to time post-deployment and these relationships were subsequently subjected to regression analysis. One-way Analysis of Variance (ANOVA) was used to test for significant differences in mean hypoxia and mean supersaturation for 1, 2, 4, 7, and 14-day treatment intervals. Linear regression analysis was used to determine whether there was a significant relationship between supersaturation and hypoxia (% of deployment with condition) and deployment duration. Data for analyses were not available for the following sites: Chesapeake Bay MD (Jug Bay-1997; Patuxent River-1997, 1998), Chesapeake Bay VA (Goodwin Island-1997), Jobos Bay (Stations 9 & 10-1997), Narragansett Bay (T-wharf-1998), Rookery Bay (Blackwater River-1997; Upper Henderson-1998), South Slough (both sites/years), Week’s Bay (Fish River-1998), and Waquoit Bay (both sites/years).

Descriptive Analyses by Reserve Site
Water depth (m), temperature (°C), salinity (ppt), dissolved oxygen (mg/L and % saturation), and deployment duration for each site at each Reserve were summarized independently using a suite of graphical data analysis techniques, which are summarized in Table 2.

  • Scatter plots for depth, temperature, salinity, and DO were created for each year of data at each site for detection of erroneous data and outlying data points.
  • Histogram plots for depth, temperature, salinity, and DO were created for each year of data at each site to examine the frequency distribution of data at each site and inter-annual variability.
  • Box plots of central tendency (mean, minimum, maximum, 90th percentile, and 10th percentile) were created for depth, temperature, salinity, and DO for each deployment at each site to examine seasonal and inter-annual variability.
  • Bar graphs for hypoxia (<28% sat) and supersaturation (>120% sat) events were created for each deployment at each site to examine seasonal and inter-annual variability.
  • Bar graphs for deployment duration were created for each site to examine the extent and duration of seasonal monitoring at each NERR site.

Reserve Comparisons
Frequency of hypoxia (<28% saturation) and supersaturation (>120% saturation) events were compared within and among reserves, within and between seasons, and between years. Duration of hypoxia among reserves was also examined. Mean daily salinity and daily salinity range were compared among reserves within and between seasons. SPSS was used for statistical testing.

* A two-sample t-test was used to determine whether hypoxia, supersaturation, mean daily salinity, and daily salinity range between June 21st and September 21st (1997, 1998) differed among sites within Reserves. If variances were heterogeneous by Levene’s test, then a two-sample t-test for unequal variances was used. Incomplete data for Chesapeake Bay MD (Patuxent River), South Slough (both sites), and Waquoit Bay (Metoxit Point) resulted in partial exclusion from analysis.

* A two-way ANOVA was used to compare hypoxia, supersaturation, mean daily salinity, and daily salinity range at Reserve sites between seasons for 1997 and 1998 data. Seasons were defined as winter (Jan-Mar), spring (Apr-Jun), summer (Jul-Sep), and fall (Oct-Dec). If variances were heterogeneous by Levene’s test and were not homogeneous following log10 transformation, a two-sample t-test for unequal variances was used to determine whether variables were significantly different between sites. A one-way ANOVA was used to determine if seasonal differences in variables occurred at each site. If variances were heterogeneous for the one-way ANOVA, then a Kruskal-Wallis test was used to test the hypothesis of no seasonal differences. Incomplete supersaturation data resulted in exclusion of the Winchester site (South Slough) and both Chesapeake Bay-MD NERR sites from analyses. Incomplete hypoxia data resulted in exclusion of the following sites from analyses: Jobos Bay (Station 10); Sapelo Island (Marsh Landing); Week’s Bay (Fish River); and both sites at the Chesapeake Bay-MD, Great Bay, Waquoit Bay, and Jacques Cousteau-Mullica River NERRs.

* Pearson’s correlation analysis was used to determine if hypoxia and supersaturation events in July-August 1997 were correlated with hypoxia and supersaturation events July-August 1998. Incomplete data for the following sites resulted in exclusion from analysis: Narragansett Bay (T-wharf); Week’s Bay (Fish River); Waquoit Bay (Metoxit Point); Chesapeake Bay VA (Goodwin Islands); and both sites from Chesapeake Bay MD, Jobos Bay, Rookery Bay, and South Slough.

* Hypoxic events for each site were sorted into one of six time bins (<4 h, 4-8 h, 8-12 h, 12-16 h, 16-20 h, 20-24 h, and >24 h) to examine the duration of hypoxia at sites in the SWMP. Deployments with missing data were excluded from analysis. Because not all sites collected similar quantities of data during the same time of year, calculations based on the frequency and duration of actual hypoxic events were made to predict annual frequency and duration of hypoxia at sites in order to compare sites using the same scale.

System-wide analysis
Representative values for depth, salinity, temperature, and dissolved oxygen (% saturation) were calculated in order to compare all 44 sites in the NERR system. Mean depth and mean salinity were calculated as the average of monthly mean values between 1996 and 1998 for all months with data. Mean annual frequency of occurrence that sites experienced warm (³ 25°C) and cold (£ 10°C) water temperatures was used to compare water temperature among sites, rather than a single mean value. Similarly, mean occurrence of summer (Jul-Aug 1997 and 1998) hypoxia and supersaturation during the first 48 hours post-deployment was used to compare DO among sites. Hypoxia and supersaturation were correlated (Pearson’s R) with each other as well as with water depth, water temperature, salinity, and latitude using SPSS®.

Hierarchical cluster analysis was used to detect groupings in the data based on a survey of 14 site attributes designed to characterize reserve sites with respect to watershed input attributes and water body (where sites were located) filtering capabilities (Appendix A). Additional attributes were included in the data based on analysis of water quality at each of the sites. Squared Euclidean distance was used as the measure for clustering sites based on attributes, while attributes were clustered using Pearson correlation. A resulting dendrogram indicated how the clusters were formed and provided a measure of the linkage distance for clustering.

Table 2. Summary of analytical techniques for characterizing water quality data.

Analysis
Technique
Parameters
Sites
Sampling Interval
 t-test

 Daily Min, Max, and Mean
for DO (% sat)		 2
DO vs. Deployment
 Line graph

 DO Extremes, Deployment
Duration		 42*
ANOVA

 DO Extremes, Deployment
Duration		 44
Regression

 DO Extremes, Deployment
Duration		 44
Descriptive

 Scatter

 Depth, Temperature, Salinity, DO (both) 44
Histogram

 Depth, Temperature, Salinity, DO (both) 44
Box plot

 Depth, Temperature, Salinity, DO (both) 44
Bar graph
 Hypoxia, Supersaturation
 44
Bar graph Deployment Duration 40
Reserve
comparisons
(Summer)

 t-test

 Salinity (range, mean), DO
Extremes		 44#
(Seasonal)

 ANOVA/
KW		 Salinity (range, mean), DO
Extremes		 44#
Bar graph

 Salinity (range, mean), DO Extremes 44#
(Inter-annual)
 Rank testing
 Hypoxia, Supersaturation 44#
Bar graph Hypoxia, Supersaturation 44#
(Duration) Bar graph Hypoxia 44
Bar graph Hypoxia (by region) 44
System Comparisons
Bar graph Depth (mean) 44
Bar graph Salinity (mean, min, max) 44
Bar graph Temperature Extremes 44
Scatter

 DO Extremes vs.
temperature, salinity, depth		 44
Scatter

 DO Extremes vs. each
other vs. latitude		 44
Correlation

 DO Extremes vs. other
parameters		 44
Cluster

 Attributes survey, all parameters 44
* See text (Effect of Deployment on Dissolved Oxygen) for listing of sites with unavailable data.
# See text (Reserve Comparisons) for listing of sites partially excluded from analyses.

Periodicity (Harmonic Regression Analysis)
The first challenge of this aspect of the data analysis was to find effective ways to look at the data with fitted models. Because data were collected every half-hour for three years, each of the 220 series (22 reserves, 2 sites each, 5 parameters) could have length on the order of 50,000 observations. Simple time plots of such data sets were too dense to reveal patterns, with the possible exception of annual periodicity. To allow for more effective visual evaluation, two interactive graphical functions (microscope and plot.week) were written in Splus®.

In microscope, the data were initially plotted in their entirety over time. The user may then select a sub-segment of the full data record to plot in greater detail. In plot.week, the data for the first week were plotted over time. The user may advance the plot to the next week's data, and so on until the series ends. Both functions can control vertical axes to allow comparability of multiple plots, and both can superimpose fitted curves on the raw data.

Harmonic regression (Bliss 1970) was used to describe the cyclic phenomena in the data and as a basis for comparing aspects of known periodicity between sites and reserves. These are methods of approximating a periodic function of known period p (i.e., 12.42 h) with weighted sums of sine and cosine waves for periods p, p/2, p/3, etc. The fractions p/2, p/3, etc. (and also their sine and cosine terms) are referred to as the harmonics of the principal period (p). The weights are usually determined by multiple regression (ordinary least squares). The technique works surprisingly well if the periodic function of interest is fairly smooth. Effects of multiple principal periods can also be modeled simply by adding terms and harmonics for each principal period.

Initial analyses attempted to model the entire data series with a large single multiple regression model using multiple harmonics of the 12.42 hour (T), 24 hour (D), 29.5 day (L), and 365.24 day (A) periods. These models also allowed for shifts in mean level of the series whenever new meters were deployed, and for deployment-specific trends. These models did not fit entire series well, even with many harmonic terms for each principal period. Simpler models did, however, fit individual deployment series well (typically on the order of 10-25 days in duration, about 500-1500 observations). The fit using similar models suggested that deployment effects on the response variable are much more complicated than simply a shift in mean level with a trend.

Lunar terms (29.5 d and harmonics) effectively accounted for spring-neap cycles, as well as having the useful property of correcting non-stationarity (i.e., breakpoints, trends) in the data; however, lunar terms created a co-linearity problem for very short deployment series. Because short series resulted in unstable coefficients, only deployments > 10 d in length were included in the harmonic regression analyses. In addition, some deployments were excluded from analyses because the data were essentially constant (usually zero) and there was no way to fit any model to a constant series.

A two-tiered approach to the analysis was ultimately selected. In the first tier, each deployment's data were modeled with its own harmonic regression. From each such fit, 35 descriptive measures were computed (Table 3). This resulted, for each variable and site, in a three-year series of approximately biweekly values for the given descriptive (one value for each deployment). After experimenting with different numbers of harmonic terms on a large number of representative series and variables, the following 37-term multiple regression model was settled upon for each of the "first-tier" analyses, approximately 10,000 deployment-level regressions:
  • Intercept (1 term)
  • Principal terms and three harmonics for 12.42 hour effects (8 terms)
  • Principal terms and three harmonics for 24 hour effects (8 terms)
  • Principal terms and three harmonics for 29.5 day effects (8 terms)
  • Pair-wise cross- products of principal harmonics between T, D, and L terms (12 terms)

The cross-product terms (at least, between T and L terms) were suggested physically because tidal cycle amplitudes varied according to a lunar cycle. Perhaps more significantly, the visual fit of the models was noticeably enhanced, especially at data extremes, by the use of cross-product terms in the regressions. Because the quality of the data was expected to be best early in the deployment, the first complete 12.42 hour cycle and 24 hour cycle were generated and their amplitude (max - min) calculated. Later profiles changed during the deployment due to cross products with the lunar terms.

In the second tier of analysis, these “descriptives” were modeled for annual periodicity and overall averages and compared between sites and reserves via analysis of variance with pre-planned multiple comparisons. If a sufficient number of deployments (³ 10) for a given site and water quality variable remained after the first stage analysis, these deployments were subsequently examined for annual periodicity in the second stage. A period of 365.24 days (with four harmonics) was used in a second harmonic regression model for each of a subset of 27 deployment descriptive variables, for each water quality variable, at each site (> 5000 regressions, but £ 60 observations each). This model was fit using least squares and annual periodicity tested using an F-test for overall model significance. If annual periodicity was detected, fitted values were calculated.

The stage-two harmonic regression fits included the following items for each site, water quality variable, and deployment descriptive variable: (1) p-value for overall annual harmonic fit; (2) estimated amplitude of the annual cycle after gap removal (amplitude = 0 if no annual periodicity present); (3) mean of fitted values after gap removal (ordinary sample mean if no annual periodicity present); (4) standard error of the mean; and (5) the number of deployments analyzed. These summary data were passed to the final stage of analysis, which featured inter-site graphical summaries and multiple comparisons between sites within reserves. Conspicuously absent from these summaries are descriptives based on first-cycle profiles (descriptive variables 27-33, Table 3). Large standard errors accompanied these “first cycle” descriptive means, and it was decided that they were probably not reliable for this reason. Within reserves, sites were compared with pre-planned multiple comparisons with an error rate of 0.05 using the Šidák (1967) method (see Hsu 1996), which gives exact inference in this setting. Reserves were not formally compared, as the large number of tests so performed would render these comparisons uselessly conservative if the overall error rate was controlled.

Table 3. Summary variables computed for deployment-level harmonic regressions.

No. Description of Summary Measure Acronym / var. name
1. NERR site No. site
2. Water Quality Variable No. qualvar
3. Root mean square error rse
4. Model R2 R2
5. Model F statistic Fstat
6. Model error df dfe
7. Number of observations nobs
8. Model intercept yhatint
9. Median of raw response variable data rawmed
10. Mean of raw response variable data rawmean
11. Range of raw response variable data rawrange
12.

 Std. Deviation of raw response variable
data		 rawstdev
13.

 Coded Date/time of first observation of deployment startdate
14.

 Coded Date/time of last obs. of deployment enddate
15. Maximum model predicted value yhatmax
16. Minimum model predicted value yhatmin
17. Mean model predicted value yhatbar
18. Median model predicted value yhatmed
19. 90th percentile of model predicted values q90yhat
20. 10th percentile of model predicted values q10yhat
21. Partial SS for T terms SST
22. Partial SS for D terms SSD
23. Partial SS for TD crossproduct terms SSTD
24.

 Partial SS for T,D, and crossproduct terms SSTDI=SST+
SSD+SSTD
25. Partial SS for T and crossproduct terms SSTI=SST+SSTD
26.
 Partial SS for D and crossproduct terms
 SSDI=SSD+SSTD
27. First full 12.42 hour cycle amplitude tidal.amp
28. First full 24 hour cycle amplitude day.amp
29.

 First cycle total amplitude from T & D terms
 to.amp
30.

 Percent of total amplitude from D
amplitude		 day.pct.amp
31.

 Percent of total amplitude from T
amplitude
 tidal.pct.amp
32. First cycle time of day for D maximum dtmax
33. First cycle time of day for D minimum dtmin
34. Percent of SSTDI due to T terms tidal.pct.var
35. Percent of SSTDI due to D terms day.pct.var

Production and Respiration (P/R)
Use of diel oxygen curve data to calculate primary production, respiration and net ecosystem metabolism (NEM) was first proposed in the 1950’s (Odum 1956). Since then, these calculations have been applied to a wide variety of aquatic systems, including many estuaries (Odum and Hoskins 1958, Nixon and Oviatt 1972, Kemp and Boynton 1980, Oviatt et al. 1986, D’Avanzo et al. 1996, Swaney et al. 1999). This method underestimates production and respiration for inter-tidal sites because atmospheric oxygen exchanges directly with sediments or macrophyte beds, rather than into the water column, at low tide; however, this method is one of the few ways to determine an integrated measure of metabolic rates for the entire aquatic system. Similarly, metabolic rates in emergent vegetation such as marshes and forests are not captured by this technique. In this study, dissolved oxygen data from the NERR SWMP (January 1996 to December 1998) was used to calculate primary production, respiration, and net ecosystem metabolism. Data from two sites at each of 14 Reserves were analyzed. Sites included in these calculations were those for which data on water volume nutrient and chlorophyll concentrations were available. Sites were selected to represent as many regions and habitats as possible (Table 4).

Oxygen is produced as a product of photosynthesis, which causes oxygen concentrations to generally increase during the day. At night, dissolved oxygen concentration decreases due to respiration. Diffusion of oxygen across the air-water interface also affects the water column oxygen concentrations. Diffusion, or air-sea exchange, was estimated by multiplying the exchange coefficient by the percent saturation (DOsat). We assumed that the exchange coefficient was 0.5 g O2 m-2 hr-1 when the oxygen concentration was zero.

Air-sea exchange = (1-(DOsat,t2+DOsat,t1)/200)* .05^*dt

Estimation of oxygen exchange using this formula is reasonable under most conditions; however, because the rate of diffusion is dependent on wind speed (Copeland and Duffer 1964, Hartman and Hammond 1984, Marino and Howarth 1993), oxygen exchange may be underestimated during periods of high winds and overestimated during calm periods. For each time interval, oxygen flux was calculated as the change in oxygen concentrations (DO) minus air-sea exchange.

Oxygen flux = (DOt2 - DOt1) – Air-sea exchange

Net production was calculated as the sum of oxygen fluxes during daylight hours. Net respiration was calculated as the sum of oxygen fluxes at night. Net production and net respiration were used to calculate gross production and total (day + night) respiration rate. Constant respiration was assumed during the day and night; thus, night respiration divided by total number of night hours represented the hourly respiration rate. Total respiration was calculated as the hourly respiration rate multiplied by 24. Gross production was equal to net production plus daytime respiration and was calculated by adding net production to the hourly respiration multiplied by the number of daylight hours. Net ecosystem metabolism (NEM) was calculated by subtracting total respiration from gross production and was reported in terms of oxygen and carbon units. Oxygen was converted to carbon assuming a photosynthetic quotient of 1.25 and a respiration quotient of 1.0 (Kemp et al. 1997).

A major assumption of the diel oxygen curve method is that water masses passing by the sensor are laterally and vertically homogenous (i.e., they have the same metabolic history). In areas where physical processes such as advection and diffusion mask site-specific biological processes, metabolic rates will likely be underestimated (Kemp and Boynton 1980). Deployments for which gross production estimates were less than zero and/or total respiration estimates were greater than zero failed to meet this assumption and were excluded from calculation of photosynthesis and respiration. To determine how effective the diurnal curve technique was at estimating metabolic rates, the percent of observations meeting the assumptions was determined for each site.

The effect of instrument drift (i.e., due to fouling) on metabolic rates (gross production, total respiration and net metabolism) was determined for each site using a paired t-test comparing the record for the entire deployment to the first 2 days of the deployment. If instrument drift was determined to be significant, then only the first 2 days of deployment were used for the subsequent statistical analyses. Average and standard error of metabolic rate estimates were calculated for each site. Correlation analysis was used to determine relationships between temperature and salinity and metabolic rates for each site. Daily rates of gross production and total respiration were compared using a paired t-test. If production and respiration rates were significantly different from one another, net metabolism was significant. Results are reported as non-significant when p>0.05.Table 4. Summary of dominant habitat near 27 sites used in ecosystem metabolism analysis.

Reserve

 Site
Dominant Habitat
near site
ACE Big Bay Salt marsh
St Pierre Salt marsh
Apalachicola Surface Estuarine
Bottom Estuarine
Chesap Bay MD Jug Bay Fresh marsh
Patuxent River Park Fresh marsh
Chesap Bay VA Goodwin Island Eelgrass
Taskinas Creek Brackish marsh
Elkhorn Slough Azevedo Pond Uplands
South Marsh Salt marsh
Great Bay Great Bay Buoy Estuarine
Squamscott River Estuarine
Hudson River Sawkill Uplands
Tivoli South Bay Fresh marsh
Jobos Bay Station 9 Mangrove
Station 10 Mangrove
Narragansett Bay Potters Cove Estuarine
T-wharf Estuarine
North Inlet-
Winyah Bay		 Oyster Landing

 Salt marsh
Thousand Acre Creek Salt marsh
Padilla Bay Bay View Eelgrass
Joe Leary Slough Upland
Rookery Bay Upper Henderson Mangrove
Blackwater River Mangrove
Weeks Bay Fish River Estuarine
Weeks Bay Estuarine
Waquoit Bay Central Basin Macroalgae