Spatiotemporal variability of light attenuation and net ecosystem metabolism in a back-barrier estuary

The light climate in back-barrier estuaries is a strong control on phytoplankton and submerged aquatic vegetation (SAV) growth, and ultimately net ecosystem metabolism. However, quantifying the spatiotemporal variability of light attenuation and net ecosystem metabolism over seasonal timescales is difficult due to sampling limitations and dynamic physical and biogeochemical processes. Differences in the dominant primary producer at a 10 given location (e.g., phytoplankton versus SAV) can also determine diel variations in dissolved oxygen and associated ecosystem metabolism. Over a one year period we measured hydrodynamic properties, biogeochemical variables (fDOM, turbidity, chlorophyll-a fluorescence, dissolved oxygen), and photosynthetically active radiation (PAR) at multiple locations in Chincoteague Bay, Maryland/Virginia, USA, a shallow back-barrier estuary. We quantified light attenuation, net ecosystem metabolism, and timescales of variability for several water properties at 15 paired channel-shoal sites along the longitudinal axis of the bay. The channelized sites, which were dominated by fine bed sediment, exhibited slightly higher light attenuation due to increased wind-wave sediment resuspension. Light attenuation due to fDOM was slightly higher in the northern portion of the bay, while attenuation due to chlorophyll-a was only relevant at one channelized site, proximal to nutrient and freshwater loading. Gross primary production and respiration were highest at the vegetated shoal sites, though enhanced production and respiration 20 were also observed at one channelized, nutrient-enriched site. Production and respiration were nearly balanced throughout the year at all sites, but there was a tendency for net autotrophy at shoal sites, especially during periods of high SAV biomass. Shoal sites, where SAV was present, demonstrated a reduction in gross primary production (GPP) when light attenuation was highest, but GPP at adjacent shoal sites where phytoplankton were dominant was less sensitive to light attenuation. This study demonstrates how extensive continuous physical and biological 25 Biogeosciences Discuss., https://doi.org/10.5194/bg-2018-335 Manuscript under review for journal Biogeosciences Discussion started: 23 July 2018 c © Author(s) 2018. CC BY 4.0 License.

individual parameter time series using either a recursive filter or median filter technique, in which values that changed from one time point to the next by more than a set threshold were flagged and replaced with a fill value (e.g. -9999). Pre-deployment calibrations and postdeployment check were performed on all EXO2 sensors following YSI procedures outlined in the EXO User Manual (Revision F). Linear corrections were obtained from either post-5 deployment calibration checks or differences between in situ values from a fouled sensor and from a cleaned, calibrated sensor at the beginning of the next deployment. If data were obviously fouled and corrections were not possible, the data were replaced with fill values. WET Labs ECO-PARSB instruments were checked to verify counts were above their "dark" count (low-count) thresholds; therefore legitimate data collected during night were also 10 discarded. All subsurface pressure data were corrected for changes in atmospheric pressure by using local barometric pressure data, from the meteorological station at CBWS when available, to give a more accurate representation of pressure caused by the overlying water." Zero values in certain parameters are common given the multi-point sensor calibrations over large ranges of variability. It is nearly impossible to accurately capture extremely high values and have 15 sensitivity at the low range. There were ice-cover events during our campaign that nearly shut down advection, resuspension, and production. Values of many parameters were zero during this time.
Time series of chlorophyll, with high concentrations during resuspension events suggest domination of suspended microphytobenthos/dead microalgae.
Agreed, we failed to note that where we discussed peaks in chlorophyll during winter resuspension 20 events in the Results section. We have added this on P11, 3-5.
Describe/discus if microalgae primary production is dominated by pelagic and/or benthic microalgae.
We do not have direct estimates of microalgal production partitioned between sediments and the water column at any of our sites, and have added the following on P14, 17-21 to constrain this uncertainty: 25 "Though we do not have direct estimates of microalgal production partitioned between sediments and the water column at any of our sites, the sediment surface is nearly covered by SAV and macroalgae at the vegetated sites, therefore we expect that benthic microalgal production is low. At site CB11, KdPAR values indicate that very little light reaches the bottom, so we presume that water-column microalgae are dominant. For site CB06, we cannot 30 constrain the partitioning between the two environments." It is concluded that the presence of SAV controls the shear stress-resuspension relation (P20, L 3). However, as the physical characteristics of especially the shallow sites suggest that microphytobenthos could be abundant, this would significantly influence the shear stress-resuspension relationship as well. Hence, SAV may not be the only/main explanation for the observed seasonal trend in shear stress-35 resuspension relation.
Agreed, the most general explanation is that depth controls the availability of light for both SAV and benthic algae; and the presence of both to some degree stabilize the bed in summer more than winter. At the deep sites, no such relationship exists, indicating no change in benthic characteristics. However, it is likely that in a vegetated sandy substrate, the seagrass would dominate bed 40 stabilization over benthic microalgae given shading and coarser sediment. In fact, the data of Ellis et increasing the light attenuation and causing a larger difference in light attenuation between the sandy, vegetated sites and the siltier, unvegetated sites. Modifying the peak value as described above from 0.017 to 0.025 (47% increase; 0.025 value was used for a different muddominated estuary by Ganju et al., 2014) increased the median light attenuation at CB06 by 17% (from 1.35 to 1.58 m-1) and at CB11 by 11% (from 1.67 to 1.86). Therefore the spatial 5 variability in light attenuation is likely robust and not confounded by variability in particle size." Determining the optimal sampling frequency is a science in itself and the results regarding the influence of sampling frequency (table 3) is interesting, but as presented it seems as an unfinished story that is not properly treated in the result and discussion section. This part should be either up-or down 10 scaled/skipped.
We believe that the concept that coarse temporal sampling masks spatial variability is an important one given the common daily sampling approach for dissolved oxygen in many impaired estuaries. We have expanded on this concept with examples from the literature in our revision, on P16, 19-23: "Many estuary sampling programs, including "citizen science" efforts, often collect one daily 15 sample of water quality parameters, including dissolved oxygen (Rheuban et al., 2016). Summers et al. (1997) demonstrated that month-long records of dissolved oxygen in a number of estuarine systems, resampled to replicate various non-continuous sampling programs, were not able to identify oxygen minima. This study expands on that finding by mimicking sampling programs within one large estuary, over an entire year." 20 The present study do not examine gradients, so this term should be avoided (e.g. P1, L 13). Use "patchiness" or "spatial variability" instead.
Agreed, we have revised throughout to "variability".
Technical correction: Table 3: State which time period (yearly, season, other?), min, max and mean is covering. Include SD for the mean values. Figure 6: add x-axis titles 25 We have revised to add time period to the caption, as well as standard deviation for the original data. We have added the x-axis (1/frequency) to Figure 6. Thank you for catching it.
Mario Hoppema (Referee) mario.hoppema@awi.de Received and published: 28 February 2020 The manuscript presents an impressive data set for four locations in an estuary, including several useful variables/parameters. The data interpretation reveals some interesting results and conclusions, which 30 are definitely worth publishing in Ocean Science.
We appreciate the support and the willingness of the Editor to provide a comprehensive review.
However, in some cases there seems to be over-interpretation, for example the difference in metabolic rates between vegetated and unvegetated sites.
In addition to our response to Reviewer #1 above, we would indicate that the daily rates are markedly 35 higher between CB3 and CB10 compare with CB6 -which is the main distinction we discussed between vegetated and unvegetated. We indicated in the revision that for mean period values, only the gross production (Pg) values between site CB03 and CB06 are statistically different.
Also the data processing should be explained in somewhat more detail.
See response above, we have added detail on the data processing on P7, 9-20.
For the net ecosystem metabolism, sometimes NEM is used, for example in Figure 11, and sometimes Pn, for example in Table 1. Please use only one single symbol or abbreviation throughout the 5 manuscript.
This is a fair comment, but the mean KdPAR is influenced strongly by the more frequent background conditions than the episodic values, which is why the mean ends up being larger at CB11, mostly caused by a higher background fDOM (Fig. 5). The partitioning of light attenuation between constituents is subject to error of the light model, but given that fDOM and chlorophyll-a are highest at that site (in terms of mean and peak values), it is expected that light attenuation would be highest 5 there (despite lower turbidity than CB06). We have revised to indicate that this is only true during periods of overlap.
P12. L13-14 "between May and October" My view of the figures says that this should be "between May and September".

Revised. 10
P12, L14 "during November to April" I was say "during November to March"

Revised.
P12, L1920 "but instances of net autotrophy (Pg > Rt) occurred nearly 70% of the time at the vegetated sites" When I view Figure 11, neither autotrophy nor heterotrophy are statistically significant. This should be mentioned here. 15 It is a fair comment to suggest that we did not include statistical analyses to "prove" the existence of net autotrophy, but rather expressed the frequency of daily net autotrophic conditions. We have mentioned this on P13, L21-22: "Gross primary production and respiration were largely balanced across all sites, but the frequency of daily net autotrophy (Pg > Rt) was higher at the vegetated sites." 20 P13, L10 at 1-7 day

Revised.
P13, L23-24 "The peak in spectral density was 30% higher at CB03 than CB10" This is really hard to see, if at all, in Fig. 6. Possibly the authors could add a comment to Figure 6 (Diss Oxygen) that the curve of CB3 lies under the one of CB10. 25 We have added a comment in the caption.
P17, L26 where the canonical C:N (I suggest to add canonical, since this is the well known Redfield ratio based on data from many locations)

Revised.
P19, L1-2 "In fact, modest net autotrophy prevailed during the summer season at vegetated sites but 30 not at un-vegetated sites." This does not follow from the data in Figure 11, where NEM is around zero all through the year. See also comment above. Actually I am surprised by the uncertainty interval of NEM, which should have been formed frrom subtracting two larger terms. Because of this, I would expect it to be clearly larger and not smaller than the uncertainty of both of the terms.
If you compare NEM to Pg and Rt on an absolute basis, we agree that NEM is near zero. But our 35 statement simply refers to the mean monthly value of NEM, which although much smaller than Pg and Rt, is still greater than zero at the vegetated sites during some seasons, which we define as autotrophic. The error bars around NEM represent the standard deviation of the monthly mean estimates of NEM, which the reviewer is correct in that they are derived as the difference between Pg and Rt. But these difference calculations are performed on a daily basis, so the error bars represent the total variation of all ~30 daily Pg, Rt, and NEM rates and thus do not represent uncertainty carried over from each individual NEM calculation. 5

Revised.
Figure 6 I think CB10 in red under the figure panels must be CB11 (unvegetated) Please indicate panels (a) -(d), and then also in the caption. Please explain the x-axis; in the main text the authors talk about "peak at 12.4 h", but this is not reflected in the figure.
Revised, the 12.4 peak, at ~ 0.5 days, is now noted in the caption. 5 Figure 7 It is not clear to me how the oxygen saturation could be 130% in the middle of winter (e.g. at CB11).
The periods where oxygen saturation exceed 130% do not truly occur in winter, which we would define as mid-December to mid-March. The exception to this is during a brief period in February and early March at CB 11, where chlorophyll-a exceeded 25 ug/L indicating that phytoplankton biomass 10 was high and presumably, primary production rates. Correspondence to: Neil K. Ganju (nganju@usgs.gov) Abstract. Quantifying system-wide biogeochemical dynamics and ecosystem metabolism in estuaries is often attempted using a long-term continuous record at a single site, or short-term records at multiple sites due to sampling limitations that preclude long-term monitoring at multiple sites. However, differences in the dominant primary producer at a given location (e.g., phytoplankton versus benthic producers submerged aquatic vegetation; SAV) 10 control diel variations in dissolved oxygen and associated ecosystem metabolism, and may confound metabolicsm estimates that do not account for this variability. We hypothesize that even in shallow, well-mixed estuaries there isare strong spatiotemporal variability gradients in ecosystem metabolism due to benthic and water column properties, the influence of submerged aquatic vegetation (SAV), and ensuing feedbacks to sediment resuspension, light attenuation, and primary production. We tested this hypothesis by measuring hydrodynamic properties, 15 biogeochemical variables (fluorescent dissolved organic matter (fDOM), turbidity, chlorophyll-a fluorescence, dissolved oxygen), and photosynthetically active radiation (PAR) over one year at 15 min intervals at paired channel (unvegetated) and shoal (vegetated by eelgrass) sites in Chincoteague Bay, Maryland/Virginia, USA, a shallow back-barrier estuary. Light attenuation (KdPAR) at all sites was dominated by turbidity from suspended sediment, with lower contributions from fDOM and chlorophyll-a. However, there was significant seasonal variability in the 20 resuspension-shear stress relationship on the vegetated shoals, but not in adjacent unvegetated channels. This indicated that KdPAR on the shoals was mediated by submerged aquatic vegetation (SAV)SAV and possibly microphytobenthos presence in the summer, which reduced resuspension and therefore KdPAR. We also found that gross primary production (Pg) and KdPAR were significantly negatively correlated on the shoals and uncorrelated in the channels, indicating that Pg over the vegetated shoals is controlled by a feedback loop between benthic 25 stabilization by SAV and/orSAV microphytobenthospresence, sediment resuspension, and light availability.
Metabolic estimates indicated substantial differences in net ecosystem metabolism between vegetated and unvegetated sites, with the former tending towards net autotrophy in the summer. Ongoing trends of SAV loss in this and other back-barrier estuaries suggests that these systems may also shift towards net heterotrophy, reducing their effectiveness as long-term carbon sinks. With regards to temporal variability, we found that varying sampling frequency between 15 min and 1 d resulted in comparable mean values of biogeochemical variables, but extreme 5 values were missed by daily sampling. In fact, daily re-sampling minimized the variability between sites and falsely suggested spatial homogeneity in biogeochemistry, emphasizing the need for high-frequency sampling. This study confirms that properly quantifying ecosystem metabolism and associated biogeochemical variability requires characterization of the diverse estuarine environments, even in well-mixed systems, and demonstrates the deficiencies introduced by infrequent sampling on the interpretation of spatial variabilitygradients.

1 Introduction
Back-barrier estuaries are biologically productive environments that provide numerous ecological, recreational, and economic benefits. Submerged aquatic vegetation (SAV) proliferates in these environments due to relatively shallow bathymetry and sufficient light availability, providing habitat for many fish and crustaceans (Heck and Orth 1980) as 15 well as enhancing wave attenuation . Primary production in back-barrier estuaries and similar shallow marine ecosystems is relatively high given the shallow bathymetry, benthic light availability, and sometimes large SAV beds (e.g., Duarte and Chiscano 1999). Benthic communities within shallow ecosystems host other primary producers where SAV are absent, including microphytobenthos (Sundbäck et al. 2000) and various forms of macroalgae, and the relative contribution of these producers is altered by nutrient enrichment (e.g., McGlathery 20 2001, Valiela et al. 1997). In deeper, unvegetated habitats, phytoplankton may also contribute significantly to primary production, where the balance between water column and benthic primary production is dependent on depth, light availability, and nutrient levels, but it is unclear if total ecosystem primary production is affected by these factors (Borum and Sand-Jensen 1996).
A fundamental control on estuarine primary production is light availability, which is affected by bathymetry for 25 benthic primary producers, but is also a function of other spatiotemporally dynamic variables (e.g., nutrient availability, sediment type). Models of light attenuation consider the role of suspended sediment, phytoplankton, and colored dissolved organic matter (CDOM) concentrations in the water column, either through empirical formulations (Xu et al. 2005) or detailed models of scattering and absorption properties (Gallegos et al. 1990).
Suspended-sediment concentrations are controlled by processes that vary on a variety of time scales (minutes, weeks, months), including bed composition, bed shear stress, resuspension, and advective inputs of sediment from external sources. In contrast, phytoplankton concentrations are a function of light and water column nutrients, while 5 CDOM is driven by the input of terrestrial material through freshwater loading. The relative contributions of these constituents to light attenuation is dependent on local conditions, can vary spatially and seasonally based on external forcings (wind-induced resuspension, external inputs), and thus requires high-frequency measurements over space and time to measure all aspects of variability.
A wealth of literature describes the relationship between light and photosynthesis for marine photoautotrophs, 10 including for SAV the role of self-shading, overall water-column conditions, and light attenuation (Sand-Jensen et al. 2006, Kemp et al. 2004, Duarte 1991. Light-photosynthesis interactions for SAV and phytoplankton can differ substantially, given that phytoplankton are vulnerable to water-column structure and mixing while SAV are rooted to sediments, SAV are vulnerable to epiphytic fouling under nutrient-enriched conditions (Neckles et al. 1993), and SAV can engineer their own light environment through attenuation of flow velocity, wave energy, and therefore bed 15 shear stress and resuspension (Hansen and Reidenbach 2013). Multiple studies have revealed self-reinforcing feedbacks within SAV beds, where SAV shoots and roots stabilize the sediment bed to reduce sediment resuspension, increasing the local net deposition of water-column particulates and improving local light conditions (e.g., Gurbisz et al. 2017). These feedbacks are complex, however, and depend on bed size, aboveground biomass, and other factors (e.g., Adams et al. 2016). Consequently, healthy SAV beds are likely to generate high rates of 20 primary production relative to adjacent unvegetated areas, but each habitat may respond differently to tidal, diurnal, and event-scale variations in physical forcing.
The spatiotemporal variability of light attenuation and primary production in back-barrier estuaries over seasonal timescales is not well-constrained due to large variabilitygradients in benthic habitats, nutrient and carbon concentrations, and circulation. Nonetheless, numerous net ecosystem metabolism (NEM) estimates have been made 25 with limited spatiotemporal information due to the inherent difficulty of continuous measurements at multiple locations. For example, Caffrey (2004) synthesized NEM across numerous estuaries, however some estuaries were represented by 1-2 sampling locations within the system. Conversely, Howarth et al. (2014) quantified NEM at multiple locations in a small, eutrophic estuary, but sampling was limited to ~ 100 d over 7 years, thereby adding uncertainty to annual rates. Indeed, Staehr et al. (2012) suggest that undersampling and variability are both large sources of uncertainty in NEM estimates that are not well-understood.
We hypothesize that the type of benthic habitat and biophysical environment control spatial variations in primary 5 production via their influence on light availability via local attenuation of light by resuspended particles and watercolumn phytoplankton biomass. Consequently, the quantification of spatial variability in ecosystem metabolism and subsequent 'up-scaling' to system-wide rates is influenced by sampling frequency. Therefore, the objective of this study is to quantify sub-hourly variations in light attenuation, water-column properties, and net ecosystem metabolism and examine the relationship among these properties across habitats in a back-barrier estuary,

10
Chincoteague Bay (Maryland/Virginia, USA) using a year-long deployment of high-frequency sensors. We first describe the observational campaign and analytical methods, followed by analysis of the light attenuation and associated forcing mechanisms. We then quantify gross primary production, respiration, and net ecosystem metabolism, which have not been studied comprehensively in this estuary, and relate it with the variability of light attenuation across different habitats. Given the large spatial variability in bed sediment type, bathymetry, and 15 dominant vegetation, we aim to quantify how temporal variations in light-attenuating substances, wave dynamics, dissolved oxygen, and metabolic rates are linked across these spatially distinct environments. Our conclusions highlight the importance of quantifying spatiotemporal variability in these processes, which indicate feedbacks between physical and ecological processes in estuarine environments that should be considered when evaluating future ecosystem response.

Site description
Chincoteague Bay, a back-barrier estuary on the Maryland/Virginia Atlantic coast ( Fig. 1), spans 60 km from Ocean 25 City Inlet at the north to Chincoteague Inlet in the south. A relatively undeveloped barrier island separates the estuary from the Atlantic Ocean. The mean depth is 1.6 m, with depths exceeding 5 m in the inlets. The central basin depths are approximately 3 m, and the eastern, back-barrier side of the bay is characterized by shallower vegetated shoals; the western side is deeper with no shoals (Fig. 1).
The coastal ocean tide range approaches 1 m, but is attenuated to less than 0.10 cm in the center of the bay, where water levels are dominated by wind setup and remote offshore forcing (Pritchard, 1960); subtidal water level 5 fluctuations throughout the bay approach 1 m (Nowacki and Ganju, 2018). River and watershed constituent inputs are minimal with the highest discharge and lowest salinities near Newport Bay in the northwest corner of the bay.
Atmospheric forcing is characterized by episodic frontal passages in winter with strong northeast winds; summer and fall exhibit gentler southwest winds. Waves within the bay are predominantly locally generated, with substantial dependence on wind direction and fetch due to the alignment of the estuary along a southwest to northeast axis.   Table 1). Instruments were recovered, downloaded, and serviced three times during that period (October 2014, January 2015, and April 2015). The sites were chosen to span a range of water depth, benthic type, and proximity to the coastal ocean; these represent a subset of the eight sites considered by Nowacki and 20 Ganju (2018) in their sediment flux study. Here we include only the sites where light climate and/or water quality were measured comprehensively. Beginning in the southern portion of the estuary, site CB03 is within a seagrass meadow (primarily Zostera marina) on the eastern edge of the southern basin, at 1 m depth. Tide range is approximately 0.3 m (Nowacki and Ganju, 2018); the seabed is mainly composed of sand and silt, with organic content < 4% (Ellis et al., 2015). Moving northward, site CB06 is within the main channel north of the nominal 25 boundary between the northern and southern basins, with a siltmud-dominated bed, organic content of 6%, and no vegetation, and a depth of 3 m, and a slightly diminished tide range of 0.2 m. Site CB10 is within a seagrass meadow on the eastern side of the northern basin, with a sand-dominated bed and < 4% organic content, at 1 m depth, and tide range of 0.15 m. Lastly, site CB11 is in the central basin of Newport Bay within the northwest portion of the estuary, at a depth of 2.5 m, and tide range of 0.2 m. The seabed at this site is silt-dominated, with the highest observed organic content (>8%). A meteorological station was deployed at site CBWS, approximately 3.8 m above the water surface at Public Landing, Maryland on the central, western shore of Chincoteague Bay. Seagrass The shallow-water platform described by Ganju et al. (2014) was deployed at sites CB03 and CB10, within sandy patches of the seagrass meadows. The platform was designed to measure hydrodynamic, biogeochemical, and light 10 parameters in the bottom half of a 1 m water column, and consists of an RBR Virtuoso D|Wave pressure recorder (+/-0.05% accuracy), a pair of WetLabs ECO-PARSB self-wiping photosynthetically active radiation (PAR, 400-700 nm; +/-5% accuracy) sensors; a YSI EXO2 multi-parameter sonde measuring temperature, salinity, turbidity, dissolved oxygen, chlorophyll-a fluorescence, fluorescing dissolved organic matter (fDOM, a proxy for CDOM), pH, and depth (parameter-dependent accuracy available at https://www.ysi.com/EXO2); and a Nortek Aquadopp 15 acoustic Doppler current profilerADCP measuring water velocity profiles (2 MHz standard at CB10 and 1 MHz high resolution at CB03; +/-1% accuracy). All instruments except for the upper PAR sensor were mounted at 0.15 meters above the bed (mab) on a weighted fiberglass grate approximately 1 m x 0.5 m. The lower PAR sensor was recessed inside a PVC tube protruding from the bottom of the frame. The upper PAR sensor was mounted at 0.45 mab; the upper and lower sensors provide an estimate of light attenuation KdPAR over the PAR spectrum (400-700 20 nm), calculated as: where dz is the distance between the two PAR sensors (0.3 m in this case). Light attenuation was calculated only between the hours of 1030 and 1530, when the angle of the sun relative to the deployment location was closest to 0 degrees. All sensors sampled at 15 min intervals, except for the wave recorders, which burst-sampled at 6 Hz every 25 3 min. Temperature, turbidity, and inner filter effects (IFE) have been shown to alter fDOM measurements (Baker, 2005;Downing et al., 2012). We corrected fDOM measurements to account for temperature, turbidity, and IFE, according to Downing et al. (2012). Measurements of fDOM at turbidities > 50 NTU were removed due to interference with the fluorescence signal. Chlorophyll-a concentration was calculated from sensor-based fluorescence measurements by regressing fluorescence against discrete measurements of chlorophyll-a made on four dates (August and October 2014, January and April 2015) at all stations (Fig. S1). In short, water was collected in the field at the time of sensor sampling, and filtered through 0.7 m GF/F filters, which were wrapped in foil, and 5 frozen until laboratory analysis using standard methods (EPA Method 445.0). Non-photochemical quenching (NPQ) was accounted for by removing fluorescence measurements during periods of daylight (upper PAR sensor> 150 W/m 2 , or ~40% of the record), and interpolating to fill gaps. Platforms at sites CB06 and CB11 were identical to platforms at CB03 and CB10 except for the omission of PAR sensors. Light attenuation at those sites was estimated using the model of Gallegos et al. (2011), discussed below; this enabled comparison between vegetated and 10 unvegetated sites using the same fundamental variables to compute light attenuation.
Quality control and quality assurance checks were performed on all data and are described in detail in the data report (Suttles et al., 2017). These included removal of obvious spikes in individual parameter time series using either a recursive filter or median filter technique, in which values that changed from one time point to the next by more than a set threshold were flagged and replaced with a fill value (e.g. -9999). Pre-deployment calibrations and post-15 deployment check were performed on all EXO2 sensors following YSI procedures outlined in the EXO User Manual (Revision F). Linear corrections were obtained from either post-deployment calibration checks or differences between in situ values from a fouled sensor and from a cleaned, calibrated sensor at the beginning of the next deployment. If data were obviously fouled and corrections were not possible, the data were replaced with fill values. WET Labs ECO-PARSB instruments were checked to verify counts were above their "dark" count (low-20 count) thresholds; therefore legitimate data collected during night were also discarded. All subsurface pressure data were corrected for changes in atmospheric pressure by using local barometric pressure data, from the meteorological station at CBWS when available, to give a more accurate representation of pressure caused by the overlying water.
We investigated the frequency response of constituents using spectral density estimates (Welch's overlapped segment averaging estimator) and their uncertainties (Bendat and Piersol 1986) for dissolved oxygen, turbidity, 25 chlorophyll-a, and fDOM. Uncertainty bounds on the spectra were used to test if differences in spectral density between sites were significant, where non-overlapping uncertainty envelops suggest significant differences. We also performed wavelet coherence analyses (Grinsted et al. 2004) between site CB03 (shoal) and CB06 (channel) for each constituent; wavelet coherence indicates co-variability of two signals at varying frequencies through time, and can identify mechanistic linkages between processes that have complex temporal coupling. The analysis was limited to these two sites due to the more complete data coverage. Statistics of the constituent time-series (mean, maxima, and minima) were computed with the original 15 min data, and with variable sampling intervals (1 h, 2 h, 24 h) to 5 investigate the influence of sampling resolution on spatial variabilitygradients in biogeochemical variables.
Further details on collection protocols and access to the time-series data are reported by Suttles et al. (2017). The hydrodynamic results of this field campaign have been described in detail by prior studies Nowacki and Ganju 2018); for the purposes of this paper we focus on the spatiotemporal variability of constituent concentrations and light attenuation.

Estimation of light attenuation contributions
We estimated light attenuation (for periods with missing PAR data or at sites with no PAR data) and relative contributions from turbidity, chlorophyll-a, and CDOM using the method of Gallegos et al. (2011). This formulation computes spectral attenuation in terms of suspended and dissolved constituents including the effects of water, CDOM, phytoplankton, and non-algal particulates (NAP, e.g., detritus, minerals, bacteria). We include absorption 15 by four components: (1) absorption by water was assumed to follow the spectral characteristics of pure water; (2) CDOM absorption was taken proportional to fDOM concentration, with a negative spectral slope ) set to sg=0.02 nm -1 (Oestreich et al., 2016); (3) phytoplankton absorption was proportional to chlorophyll-a concentration and with the spectrum shape normalized by the absorption peak at 675 nm (initial value for peak absorption was taken as a=0.0235 m 2 (mg chl a) -1 , within the range provided by Bricaud et al. 1995); and (4) 20 non-algal absorption was taken as proportional to the suspended-sediment concentration with a spectral shape (Bowers and Binding, 2006) that included a baseline of cx1=0.0024 m 2 g -1 (Biber et al. 2008), an absorption crosssection of cx2=0.04 m 2 g -1 (Bowers and Binding 2006), and a spectral slope of sx=0.009 nm -1 (Boss et al. 2001). The backscattering ratio of water was set at 0.5, while CDOM is considered non-scattering (Mobley and Stramski 1997), and the particulate effective backscattering ratio bbx was initially set at 0.017.

25
Given the high variability in turbidity and suspended-sediment concentrations due to wind-wave resuspension, we varied the value of bbx as a function of turbidity to achieve the best agreement between the model and observations.
At turbidity below 50 NTU, bbx is constant at 0.017, and then linearly declines to 0.0024 at turbidities above 220 NTU. This modification implies that backscattering by mineral particles dominates at low turbidities, while large organic aggregates dominate at the highest turbidities (Gallegos et al. 2011), and has the effect of enhancing attenuation by particulate-induced absorption at the highest turbidities. The lowest values of backscattering to scattering ratio bbx in the literature range from 0.005 (Snyder et al. 2008) to about 0.002 (Chang et al. 2004). The 5 chosen value of 0.0024 was obtained from Loisel et al. (2007). Morel and Bricaud (1981) described the bbx ratio as decreasing with increasing absorption, which would be consistent with high turbidity situations. McKee et al. (2009) described a decrease of the bbx ratio with increasing concentrations in a mineral-rich environment. The maximum sediment concentrations (15 mg/l) in that study were higher than previously mentioned studies, but smaller than concentrations observed in this environment (Nowacki and Ganju 2018). Determining the appropriate backscattering 10 ratio at high turbidities is still poorly constrained and the minimum value used in this study might even be an overestimation (a sensitivity analysis is described in section 3.1). The calibration of the light model at shoal sites (CB03 and CB10) and subsequent application at channel sites assumes that the vertical variability in the bottom 0.3 m of the water column is similar among these sites.

15
The basic concept and method for computing community production and respiration (and ecosystem metabolism) was developed by Odum and Hoskin (1958) and, with numerous modifications, has been used since for estimating these rate processes in streams, rivers, lakes, estuaries and the open ocean (Caffrey 2004;Howarth et al. 2014). The technique is based on quantifying increases in oxygen concentrations during daylight hours and declines during nighttime hours as ecosystem rates of net primary production and respiration, respectively. The sum of these two 20 processes over 24 h, after correcting for air-sea exchange, provides an estimate of net ecosystem metabolism. We The oxygen data used to make metabolic computations were obtained from sensors deployed near-bottom in relatively shallow waters. Our metabolic computations assume that the water-column is well-mixed, which is 10 necessary for the air-water flux correction to be valid and for the oxygen time-series to be representative of the combined water-column and sediments (e.g., Murrell et al. 2018). We tested the validity of this assumption at our two deeper sites (which were 2-3 m deep) using monthly vertical profile data for dissolved oxygen, temperature, and salinity from the Maryland Department of Natural Resources from 1999-2014. These data indicate instantaneous vertical dissolved oxygen differences of >1 mg/L on occasion over the 15 y record, generally during the productive 15 summer months. The long-term mean vertical oxygen difference, however, was less than 0.5 mg/L during September to May, but between 0.5 and 0.9 mg/L during June-August, indicating that although we did not measure vertical gradients during our study, they do occasionally develop.

Time-series of turbidity, chlorophyll-a, fDOM, dissolved oxygen, and light attenuation
Turbidity ranged from near zero to a maximum of over 400 NTU at site CB06 during a winter storm (December 2014) that induced waves exceeding 0.7 m (Figs. 2-5). The highest peak turbidities were observed at CB06, while sites CB03, CB10, and CB11 had similar statistical distributions of turbidity (Table 3), though the differences 25 between mean values over the entire period were not statistically significant (Table 3). Turbidity at shoal sites was highest in the winter, between November-April, while turbidity at site CB06 did not display a similar seasonal signal. In general, tidal resuspension appeared to be minimal although tidal advection after large wind-wave resuspension events was observed . Spectrally, most of the energy in the turbidity signal was found at subtidal frequencies (i.e., > 1 day), with little correspondence between sites at tidal frequencies (Fig. 6).
Chlorophyll-a concentrations peaked just below 50 g/L at site CB11, and below 30 g/L at all other sites (Figs. 2- 5). Mean Cchlorophyll-a concentrations were not statistically different comparable across sites CB03, CB10, and CB06, but were nonetheless on average twice as high at CB11 (Table 1)

10
Despite removing effects of non-photochemical quenchingNPQ, an attenuated diel signal remains in the spectra that may be a partial artifact or consistent with day-night variations in PAR (Fig. 6). The diurnal signal was comparable across all stations and was stronger than the tidal signal (Fig. 6).
Concentrations of fluorescing dissolved organic matter (fDOM) varied between 0 and 70 QSU, with the highest peak values consistently over the overlapping period of record observed at CB11 in Newport Bay (Figs. 2-5).

15
Concentrations were similar between the other sites, with lowest fDOM in the winter, possibly due to reduced dissolved organic material production in the surrounding watershed, and ultimately in freshwater runoff to the estuarybiological activity. Periodic large decreases in fDOM coincided with increases in turbidity due to wind-wave resuspension; despite correcting the fDOM time-series for turbidity interference, the fluorescence measurement was likely attenuated beyond correction. Most of the energy in the fDOM signal was at subtidal frequencies ( Fig. 6),

20
there was a distinct peak at the M2 tidal frequency (12.42 h or ~ 0.5 days) corresponding to advection of fresher, fDOM-elevated water on ebb tides.
Dissolved oxygen percent saturation ranged from a minimum of 19% at site CB11 in the summer, to a maximum of over 200% at site CB10 (Fig. 7), with maximum diel variability in the summer. The channel sites (CB06 and CB11) showed substantially attenuated diel variations as compared to the shoal sites (CB03 and CB10). Diel fluctuations in 25 the winter were typically less than the subtidal changes. Spectral analysis clearly shows the dominance of 24 h diel fluctuations at all sites (Fig. 6), with higher energy at the vegetated shoal sites (CB03, CB10).
Direct measurements of KdPAR at sites CB03 and CB10 were partially confounded by instrument fouling and malfunction of one or both sensors on each platform. Nonetheless, we successfully captured a wide range of conditions, with peak light attenuation of approximately 10 m -1 occurring at sites CB03 and CB10 in the winter during a sediment resuspension event . Median KdPAR was approximately 1 m -1 at both shoal sites, but event-driven magnitudes exceeded 2 m -1 multiple times during the deployment. The KdPAR data were used to 5 calibrate the light model (Fig. S2), which we implemented to reconstruct missing data at sites CB03 and CB10, and to estimate light attenuation at sites CB06 and CB11 where constituent concentrations were measured but no light data were collected. The full time-series of estimated light attenuation at the four sites demonstrates the strong seasonal variability of light attenuation at all sites (Figs. 2-5), with peak attenuation occurring during winter storms.
Wave-induced sediment resuspension and advection were responsible for increased turbidity, which accounted for 10 approximately 40% of the light attenuation at sites CB03, CB10, and CB11, and 61% at site CB06. Light attenuation at site CB11, with its proximity to freshwater and nutrient sources, was highest overall during the overlapping period of record and more highly influenced by chlorophyll-a and fDOM than at other sites. Median KdPAR was highest at the two unvegetated channel sites, and lowest at the two vegetated shoal sites (Table 2), though the period-means are not statistically different. Turbidity, the strongest driver of KdPAR , was generally lower during the warm season at the 15 vegetated shoal sites and the turbidity generated for a given wind-wave induced bed shear stress was significantly reduced in summer at the vegetated sites, while there was no significant seasonal variation in the turbidity-shear stress relationship at unvegetated channel sites (Fig. 8).
Sediments at site CB06 and CB11 tend to be finer (Ellis et al., 2015) and may be more susceptible to flocculation, which would induce error in the light model. Therefore we ran the light model with variation in the maximum 20 particle backscatter ratio (bbx) which is particle size dependent. An increase in flocculated suspended-sediment would increase this ratio, thereby increasing the light attenuation and causing a larger difference in light attenuation between the sandy, vegetated sites and the siltier, unvegetated sites. Modifying the peak value as described above from 0.017 to 0.025 (47% increase; 0.025 value was used for a different mud-dominated estuary by Ganju et al., 2014) increased the median light attenuation at CB06 by 17% (from 1.35 to 1.58 m-1) and at CB11 by 11% (from 25 1.67 to 1.86). Therefore the spatial variability in light attenuation is likely robust and not confounded by variability in particle size. Wavelet coherence of dissolved oxygen between site CB03 and site CB06 was maximized at diel timescales (Fig. 9), corresponding to co-varying oscillations due to daytime production and night-time respiration. Coherence was also high at lower frequencies, especially during winter when oscillations at both sites were small due to reduced production and respiration. Turbidity was coherent between channel and shoal sites primarily at subtidal timescales, corresponding to episodic multi-day resuspension events and generally high turbidity during most of the winter.

5
Phase lag was minimal during times of high coherence for both of these parameters. Coherence for both chlorophylla and fDOM was minimal throughout the year, though the former demonstrated sporadic coherence at diel and multi-day timescales. The December-April period of generally high chlorophyll-a was similar between these sites and manifested as increased coherence at ~512-1024 h (21-42 d).
Statistics of the constituent time-series were computed with temporal sampling intervals of 15 min (original 10 sampling interval), 1 h, 2 h, and 24 h ( Figure 10;  (Table 3). The effect of these sampling intervals is discussed later in the manuscript.

Net ecosystem metabolism
Estimates of ecosystem metabolism displayed strong seasonal variability, with elevated rates of Pg and Rt during warm months across all stations (Figs. 11, 12; Table S1). High rates of Pg and Rt persisted between May and 20 OctoberSeptember, when temperature and PAR peaked seasonally, and were consistently lower during November to April March (Figs. 11, 12), and Rt was exponentially related to temperature across sites (although less so at CB06).
Daily Mmetabolic rates were clearly higher at vegetated shoal sites (CB03, CB10), consistent with the strong diurnal signal in oxygen at these sites (Figs. 6, 7) and the presence of seagrass (Table 1). Temperature was strongly associated with rates of respiration across all sites, and respiration reached peak values under conditions of high 25 temperatures and high rates of Pg (Fig. 12). Gross primary production and respiration were largely balanced across all sites, but the frequency of instances of daily net autotrophy (Pg > Rt) wasoccurred highernearly 70% of the time at the vegetated sites, and were also persistent at CB06 (Fig. 13). Pg and Rt were more balanced at CB11, but we did not have enough data to capture the complete seasonal cycle at this site (Fig. 11).
The relationship between rates of Pg and light availability was site specific. At the vegetated sites, CB03 and CB10, Pg and KdPAR were significantly negatively correlated, with high rates of Pg occurring during the periods of lowest KdPAR and highest surface PAR (Fig. 14). At CB06 and CB11, variations in Pg and KdPAR were not significantly 5 correlated and KdPAR was higher and Pg was lower at these sites than at CB03 and CB10 (Tables 1, 2; Figs. 13, 14).
Over the entire record, mean values were not statistically different except for mean Pg between sites CB03 and CB06 (Table 1). In summary, the highest daily metabolic rates we measured occurred during warm periods in vegetated shoals, where wave attenuation by seagrass reduced turbidity and KdPAR.

4 Discussion
Net ecosystem metabolism and other water-quality parameters are typically measured over limited spatial (e.g. single point) and temporal (e.g. days-to-weeks) scales. An analysis of a comprehensive suite of high-frequency biological and physical measurements in Chincoteague Bay over an annual cycle revealed the primary drivers of light attenuation, the role of light attenuation in driving variations in gross primary production, the primary 15 timescales of biogeochemical variability, and the effect of habitat type (i.e., vegetated versus un-vegetated; nutrient enrichment) on oxygen variability and net ecosystem metabolism. Turbidity dominated light attenuation variability and varied considerably at 1-7 day time scales, consistent with the frequency of storm passage. Storm-associated wind waves were the specific driver of resuspension and turbidity, and reduction of bed shear stress and turbidity in SAV-dominated shoal environments during summer increased light availability in these habitats. As a consequence,

20
Pg was substantially higher in SAV-dominated shoals compared with adjacent plankton-dominated sites and Pg was negatively correlated with KdPAR in the shoals, highlighting the role of short-term variability in light availability driving Pg. High rates of Pg and Rt in shoal environments led to much higher diurnal and seasonal-scale variability in dissolved oxygen in these habitats. Though we do not have direct estimates of microalgal production partitioned between sediments and the water column at any of our sites, the sediment surface is nearly covered by SAV and 25 macroalgae at the vegetated sites, therefore we expect that benthic microalgal production is low. At site CB11, KdPAR For site CB06, we cannot constrain the partitioning between the two environments.

Interpretation of spectral signals
Spectral analysis of the constituent time-series elucidates significant differences between channel and shoal sites, as 5 well as longitudinally within the estuary. With regards to dissolved oxygen, the peak in diurnal energy at shoal sites is markedly stronger relative to channel sites with no overlap in uncertainty bounds between channel and shoal sites, indicating higher local production/respiration due to SAV presence. Between sites CB03 and CB06, the diurnal peak in spectral density is reduced by 94%, and reduced by 83% between site CB10 and CB11. The peak in spectral density was 30% higher at CB03 than CB10, perhaps due to higher SAV biomass at that site. This indicates the 10 dominant role of seagrass in controlling oxygen concentrations locally, and the relatively weaker oxygen dynamics at the channelized, unvegetated sites. From a sampling perspective, this also suggests that interpreting oxygen dynamics with limited spatial resolution is confounded by heterogeneous benthic coverage.
The highest spectral densities for turbidity were observed at site CB06, across all frequencies. At the lowest frequencies, corresponding to multi-day storms, there was no overlap in the uncertainty bounds. Nowacki and Ganju

15
(2018) showed that sediment fluxes at this site were an order of magnitude larger than the other three sites, due to a strong response to wind events, which are typically multi-day events. The central location of the site suggests that local sediment transport would respond consistently to the seasonally variable winds, which tend to act along the axis of the estuary. Resuspension over shoals on either end of the estuary, and subsequent advection cause the spatial integration of resuspension processes throughout the estuary at this main channel site.

20
Spectral density peaks of chlorophyll-a were similar between sites with significant overlap across frequencies, but site CB11 demonstrated the highest total spectral density, primarily at subtidal frequencies. The reduced exchange with other portions of the estuary, local nutrient inputs, and longer residence times increase the potential for phytoplankton growth at this station, which is known to be locally eutrophic (Fertig et al. 2009). This coincides with the lowest oxygen minima across the four sites (Table 3), indicating eutrophication. Again, in regards to sampling 25 strategies, characterizing spatial variabilitygradients in eutrophication is confounded by the competing timescales of transport and biogeochemical processes.
While the tidal signal in fDOM was pronounced (peak at 12.42 h or 0.5 days), the bulk of the spectral density was in the lower frequency band at sites CB10 and CB11, in the northern half of the estuary. This likely represents a longitudinal gradient in freshwater input, with decreased salinity in the north where inputs are largest. The inverse 5 correlation between salinity and fDOM in Chincoteague Bay was identified in prior work (Oestreich et al. 2016), and highlights the roles of tidal circulation, residence time, and freshwater input in biogeochemical and light attenuation patterns.

Channel-shoal coherence of water quality parameters
10 Wavelet coherence analyses demonstrate that though certain parameters are tightly coupled between channel and shoal, other parameters are not. The coherence of dissolved oxygen at diel and longer timescales, along with nearly no phase lag, is congruent with diel and seasonal patterns of production and respiration. The "breathing" of the estuary (Odum and Hoskin, 1958) appears as a near-simultaneous process regardless of location, and demonstrates the importance of local biogeochemical processes throughout the system. With regards to turbidity and sediment 15 resuspension, the coherence at subtidal timescales indicates that channels integrate resuspension processes over the entire estuary; minimal phase lag suggests that a rapidly evolving local wave climate in response to episodic winds controls both resuspension and rapid advection through the channels (Nowacki and Ganju, 2018). Conversely, chlorophyll-a demonstrates coherence only over the seasonal timescale, and therefore advection and other tidal processes do not appear to control dynamic exchange of phytoplankton or benthic microalgae. Lastly, the lack of 20 coherence in fDOM shows that freshwater input is relatively minor in the southern half of the estuary, and despite a significant peak in spectral density at the tidal timescale, the advection of freshwater (and fDOM) between channel and shoal is not significant. Despite a qualitatively coherent annual signal in fDOM, the analysis cannot detect coherence at this timescale with a single year of data. These complex relationships between parameters highlight both the intricate response to spatiotemporal variability in biophysical environments, and the necessity of high-25 frequency, spatially comprehensive measurements to understand estuarine ecosystem behavior.

Influence of temporal sampling resolution on spatial variabilitygradients
Many estuary sampling programs, including "citizen science" efforts, often collect one daily sample of water quality parameters, including dissolved oxygen (Rheuban et al., 2016). Summers et al. (1997) demonstrated that month-long records of dissolved oxygen in a number of estuarine systems, resampled to replicate various non-continuous sampling programs, were not able to identify oxygen minima. This study expands on that finding by mimicking 5 sampling programs within one large estuary, over an entire year. Though mean biogeochemical values were relatively insensitive to sampling interval, maximum values were significantly modulated for all parameters, while minimum values for dissolved oxygen were also affected. With regards to dissolved oxygen specifically, we find that daily sampling dampens the spatial variability in maxima and minima between sites. For example, a daily sampling program would not detect hypoxic conditions at site CB11, and would only observe a 13% difference in 10 dissolved oxygen as compared to nearby site CB10, while the actual difference was over 30%. This is relevant given the ubiquitous daily sampling programs in many estuaries, which cover numerous locations with infrequent sampling. This result suggests that characterizing differences in water-column conditions across space requires sampling at timescales finer than 1 day, especially in highly metabolic environments where diel oscillations in dissolved oxygen are large. Continuous monitoring at multiple locations in an estuary is difficult due to fouling and 15 instrument limitations (e.g. power and memory), especially across large systems over an entire year. However, it may be more informative to conduct shorter term, continuous monitoring deployments at select locations rather than daily monitoring at many locations. In the case of Chincoteague Bay, characterizing the oxygen environment would be possible with seasonal, month-long deployments at four sites. This type of data would yield both an accurate assessment of minimum and maximum values, and allow for spectral analysis to capture tidal timescale variability.

Spatiotemporal variability of light attenuation
Light attenuation can appear to be controlled by uncoupled biological (i.e., nutrient loading and phytoplankton blooms) and physical (sediment resuspension) processes, but in reality the feedbacks between physics and biology are consistently present in estuaries. Wave-induced suspended-sediment resuspension is primarily responsible for light attenuation in shallow lagoons like Chincoteague Bay, and seagrass clearly modulates the magnitude and 25 spatiotemporal variability in sediment resuspension through the attenuation of wave energy (e.g., Hansen and Reidenbach 2013). This wave-induced resuspension and turbidity occurs on the time scale of periodic wind events observed in this system. During summer months, when vegetation densities are ostensibly the highest, the dependence (i.e. slope) of turbidity on bed shear stress decreases significantly at the vegetated shoal sites, while it is not significantly different at the unvegetated channel sites (Fig. 8). In the winter, when seagrass is largely absent, the slope increases at the shoal sites, indicating the influence of seagrass on bed stabilization in the summer. The role of seagrass in modulating their physical environment is demonstrated by this improvement in light availability when 5 vegetation re-establishes in the warmer months. Furthermore, the dependence of Pg on KdPAR at the vegetated shoal sites (Fig. 14) completes the positive feedback loop between benthic stabilizationseagrass, sediment resuspension, light attenuation, and primary production. It is important to note that the seasonal change in resuspension response at the shallow sites may also be partially attributed to increased microphythobenthos coverage in summer. However it is likely that in these vegetated sandy substrates, seagrass would dominate bed stabilization over benthic microalgae 10 given shading and coarser sediment. In fact, the data of Ellis et al. (2015) indicate organic matter percentages of less than 2% on the sandy shoals.
Our measurements underscore the ability of comprehensive continuous measurements to capture multiple scales of variation in light attenuation. For example, subsampling our light attenuation measurements under fair weather conditions (wave height < median) leads to underestimation of median KdPAR by as much as 27% (sites CB06 and 15 CB10). Subsampling to periods with wave heights less than the 84 th percentile leads to underestimation by 13% (site CB06). The effect of these changes can be large if sparse data are used to assess trends in biogeochemical variables, or to drive ecological models. In situations where continuous monitoring of light attenuation is difficult, light modeling using continuous measurements of turbidity, chlorophyll-a, and fDOM is a suitable proxy. With a reasonable number of discrete light attenuation samples, it is possible to calibrate a light model and estimate 20 seasonal changes in contributions to light attenuation using the constituent measurements. This becomes useful when attempting to determine the relative influence of physical and biogeochemical processes on spatiotemporal variations in light attenuation, and the potential benefits of restoration activities.
Apart from Newport Bay in the northwest corner of the estuary, Chincoteague Bay is relatively unimpacted by external nutrient and organic material inputs, as compared to other lagoonal systems on the U.S. East Coast (e.g.

25
Indian River Lagoon, Phlips et al. 2002;Great South Bay, Kinney and Valiela 2011). For example, the relatively urbanized watershed surrounding Barnegat Bay, New Jersey, contributes larger nutrient loads to the northern portion of the estuary (Kennish et al. 2007), however the highest light attenuation is observed in the southern portion, where sediment concentration are highest due to the availability of fine bed sediment and wind-wave resuspension ). In fact, in Barnegat Bay light attenuation is lowest in the more eutrophic, poorly flushed (Defne and Ganju 2015) northern portion due to coarser bed sediments and smaller wave heights (and less resuspension). In the more nutrient-enriched regions of the Maryland coastal bays north of Chincoteague Bay, chlorophyll-a is also much 5 higher and likely makes a larger contribution to light attenuation (Boynton et al. 1996). In contrast, spatial variabilitygradients in the light field within Chincoteague Bay isare driven by geomorphology (depth), the presence of submerged aquatic vegetation (CB03 and CB10), and the higher contribution of chlorophyll-a at a nutrientenriched site (CB11; Boynton et al. 1996). Thus, generalizing the light attenuation in back-barrier estuaries is hampered by these subtle, habitat-specific controls on attenuation, but improvements in continuous monitoring will 10 lead to an increased understanding of these controls.

Spatiotemporal variability in metabolism and relationship with benthic ecosystem and light attenuation
Clear differences in the magnitude and timescales of oxygen dynamics were present across habitats in Chincoteague Bay. Spectral density at the diurnal time-scale for dissolved oxygen was significantly higher at the two vegetated 15 sites, where high rates of primary production during the day and associated respiration at night led to large changes in dissolved oxygen. Although this pattern is unsurprising given the high metabolic rates in these dense seagrass beds, this study is one of few studies that have clearly documented this pattern inside and outside of SAV beds within an estuary over annual timescales. Elevated metabolic rates in macrophyte dominated habitats relative to phytoplankton dominated habitats have been documented in other systems (D'Avanzo et al. 1996), where the 20 canonical C:N molar ratio of phytoplankton (C:N = 6.6) is much less than that measured for Chincoteague Bay Z. marina (C:N = 25.8±7.2), indicating higher carbon and oxygen metabolism for a given amount of nitrogen (e.g., Atkinson and Smith 1983). Elevated magnitudes of daily oxygen change have also been associated with elevated phytoplankton biomass in Chincoteague Bay (e.g., Boynton et al. 1996) and other systems (e.g., D' Avanzo et al. 1996). This is consistent with the fact that the site at Newport Bay (CB11), which had the highest plankton 25 chlorophyll-a concentrations across all sites, had a larger spectral density for dissolved oxygen at the diel frequency than at CB06. Newport Bay is known to have elevated nutrient inputs and concentrations associated with land-use in its watershed (Boynton et al. 1996), and spectral density for chlorophyll-a was highest at this location (Fig. 6), suggesting that this site is responding to eutrophication.
Seasonal variations in metabolic rate estimates were typical of temperate regions, but were variable across space. At all stations, rates of Pg and Rt were highest in warmer months when incident PAR and temperature are highest in this region (e.g., Fisher et al. 2003), but rates were especially high in the later summer (August and September) at the 5 vegetated sites. This late-summer period is typically when SAV biomass is highest and where the self-reinforcing effect of wave attenuation, reduced resuspension, and reduced KdPAR allow for high rates of primary production (Figs. 8,14). Indeed, the relationship between KdPAR and Pg was strongest at these vegetated sites (Fig. 14), displaying the sensitivity to benthic primary production to light attenuation, which was dominated by turbidity.
Elevated temperatures during this period also stimulate respiration (e.g., Marsh et al. 1986), which has been 10 documented at the ecosystem scale in other SAV-dominated coastal systems (Howarth et al. 2014, Boynton et al. 2014). Seasonal variations in metabolic rates were comparatively lower at the unvegetated, non-eutrophic site (CB06), where low phytoplankton biomass and the absence of macrophytes leads to limited primary production and respiration rates.
The metabolic balance (ratio) between primary production and respiration is a metric of interest in coastal aquatic 15 ecosystems as it provides an indication of whether the system is a relative source or sink of carbon with respect to the atmosphere (Staehr et al. 2006). Many river-dominated estuaries and estuaries flanked by extensive tidal marshes are net heterotrophic, where respiration exceeds primary production (e.g., Kemp and Testa 2011). Pg and Rt are often balanced or indicate modest autotrophy in shallow lagoons where SAV is present (e.g., Ferguson andEyre 2010, D'Avanzo et al. 1996). Autotrophy tends to persist in these habitats as submerged macrophytes like SAV generate 20 high biomass where light availability is relatively high given low nutrient concentrations, low phytoplankton biomass, and sediment trapping within SAV beds. In fact, modest net autotrophy (on a monthly basis) prevailed during the summer season at vegetated sites but not at un-vegetated sites. The rates of Pg and Rt we estimated, which were typically between 200 and 400 mmol O2 m -2 d -1 , are comparable to those measured in other temperate ecosystems dominated by Zostera marina (Howarth et al. 2014) and net autotrophy in these low-nutrient 25 environments is consistent with other coastal lagoons in the region (Giordano et al. 2012, Stutes et al. 2007). While the overall metabolic balance (Pg/Rt) of Chincoteague Bay is difficult to assess, given large differences in metabolic balance across a range ofgradients of depth and nutrient enrichment, we made a simple calculation where mean rates of net ecosystem metabolism at CB10 were multiplied by the area of SAV in 2015 (30 km 2 ) and rates from CB06 are multiplied by the remaining area without SAV (276 km 2 ). The calculation assumes that rates measured at these stations are representative of similar habitats across Chincoteague Bay, and does not account for the potential of benthic photosynthesis to occur in shallow regions not occupied by SAV. This approach generates an estimate of net 5 metabolism for Chincoteague Bay of 317±461 Mmol O2 yr -1 , which indicates autotrophy basin-wide, but the uncertainty is high. Given ongoing loss of SAV in this and other back-barrier systems, it is likely that the metabolic balance may shift towards heterotrophy in the future.
Pg and Rt were tightly coupled in the Chincoteague Bay stations during 2014-2015, suggesting that primary producers were the dominant sources of respiration in the ecosystem. At the vegetated stations, low plankton 10 chlorophyll-a levels indicate that SAV were the dominant contributors to ecosystem metabolism, though macroalgae and epiphytic algae were also present. Alternatively, high chlorophyll-a levels at CB11 (Newport Bay) suggest that phytoplankton were dominant contributors to metabolism, as this site is known to be nutrient enriched. Alternate sources of primary production and respiration that could drive metabolism include benthic micro-and macro-algae and heterotrophic bacteria. Previously measured sediment oxygen uptake (SOU) rates during summer in sediment 15 incubations without light suggest a mean SOU of 56.2 (±24.7) mmol O2 m -2 d -1 (Bailey et al. 2005) which is ~50% of the mean summer respiration rates at CB06, but only 15% of rates measured at the phytoplankton-dominated CB11. Rates of benthic photosynthesis were not available at this site, but with a mean KdPAR of 1.35 at CB06 and a water-column depth of 3 m, only a small fraction of surface PAR would be expected to reach the sediments. Thus, respiration at CB06 might be evenly split between the water-column and sediments, which fits well with previous 20 cross-system comparisons of coastal marine ecosystems (Kemp et al. 1992, Boynton et al. 2018).

Conclusions
We tested the hypothesis that spatiotemporal differences in benthic habitat would exert a strong control on biogeochemical variables and primary production in an otherwise well-mixed estuary. We found a clear linkage 25 between SAV presence/absence, dependence of turbidity on shear stress, the ensuing light attenuation, and gross primary production, reinforcing the positive feedback loop between seagrass presence and primary production.
Vegetated shoal sites exhibited higher metabolic rates, reduced sediment resuspension, and reduced light attenuation as compared to unvegetated channel sites. Light attenuation was dominated by wind-wave induced sediment resuspension, which was minimized during summer months in vegetated shoal sites when SAV aboveground biomass was at its peak. Furthermore, we found a clear causal linkage between SAV presence, the reduction of turbidity at a given shear stress, lowered light attenuation, and elevated gross primary production, reinforcing the 5 positive feedback loop between seagrass presence and primary production. As a consequence, despite relatively balanced ecosystem production and respiration across all sites, high rates of primary production led to modest net autotrophy at sites vegetated by SAV. Future shifts towards phytoplankton-dominant habitat, given ongoing SAV loss, will likely shift these systems towards net heterotrophy and reduce their effectiveness as carbon sinks.
Interpreting spatial differences in biogeochemistry was strongly influenced by temporal sampling resolution, where 10 daily sampling reduced spatial variability between sites and attenuated peak values significantly. These results demonstrate the value of high-temporal resolution measurements at multiple locations within a given estuary, due to the strong interplay between geomorphology, light attenuation, SAV dynamics, and ecosystem metabolism. The mechanisms identified here demonstrate a need to consider feedbacks between biological and physical processes in estuaries, especially when constructing deterministic models or evaluating future ecosystem responses to 15 eutrophication and climate change.

Code availability
Freely available MATLAB toolboxes, as cited in the text, were used for analyses.

Competing interests
The authors declare that they have no conflict of interest.

Disclaimer
Use of brand names is for identification purposes only and does not constitute endorsement by the U.S. Government.        Shaded areas indicate 90% uncertainty bounds of spectral density estimates.

20
Shaded areas near bottom of each panel indicate areas that are not statistically robust and should be interpreted with caution. Figure 9. Wavelet coherence between time-series at sites CB03 (vegetated shoal) and CB06 (unvegetated channel), for four water quality parameters. Increased coherence at a given period indicates co-variability of the time-series at that time; for example, increased coherence at ~24 h for dissolved oxygen during most of the record indicates co-

25
varying diel oscillations at both sites during most of the year. Direction of arrows indicates phase; arrows pointing to the right indicate that signals are in phase. Figure 10. Mean (dots), maxima and minima (bars) for each parameter using different temporal sampling intervals (15 min, 1 h, 2 h, 24h). Spatial variability in dissolved oxygen is most impacted by coarse temporal resolution, with spatial differences in minima and maxima largely eliminated at resolution of 1 d.
5 Figure 10. Mean (dots), maxima and minima (bars) for each parameter using different temporal sampling intervals (15 min, 1 h, 2 h, 24h). Spatial gradients in dissolved oxygen are most impacted by coarse temporal resolution, with spatial differences in minima and maxima largely eliminated at resolution of 1 d.       Table 1. Water-column properties at the four study sites, SAV presence/absence, and metabolic rate estimates.
Temperature data are minimum and maximum values (min, max), salinity is the annual mean, chlorophyll-a is the annual mean, and metabolic rates estimates are means (±standard deviationSD) for the August-September (peak SAV growth) period. Standard deviations of all parameters are sufficiently large such that differences between means are not statistically significant except for Pg at site CB06 relative to site CB03.

CB06 CB10 CB11
Depth (       for four water quality parameters; A) dissolved oxygen, B) turbidity, C) chlorophyll-a, and D) fDOM. Increased coherence at a given period indicates co-variability of the time-series at that time; for example, increased coherence at ~24 h for dissolved oxygen during most of the record indicates co-varying diel oscillations at both sites during 5 most of the year. Direction of arrows indicates phase; arrows pointing to the right indicate that signals are in phase.
Shaded areas near bottom of each panel indicate areas that are not statistically robust and should be interpreted with caution. Figure 10. Mean (dots), maxima and minima (bars) for each parameter using different temporal sampling intervals (15 min, 1 h, 2 h, 24h). Spatial variabilitygradients in dissolved oxygen isare most impacted by coarse temporal resolution, with spatial differences in minima and maxima largely eliminated at resolution of 1 d.

5
Figure 12. Relationship between respiration (Rt) and water temperature at four sites, with data coloration scaled to gross primary production (Pg).