For fully formatted methods, please see attached protocol.pdf document
14C Production Methods
The approaches for estimating primary production in the study lakes using 14C
incubations differed slightly between the three research programs, but all resulted in
a similar estimate of daily epilimnetic pelagic production (mmol C m-3 d-1). In NTL
lakes, integrated samples of water from the surface of the lake to the bottom of the
epilimnion were collected between 2007 and 2013 using a 1.5 inch PVC tube
approximately every two weeks during the open water season (first described in these
lakes by Adams et al. 1993). Samples were labeled with inorganic 14C in the form of
NaHCO3 and then incubated in the lab for 3-hr across a range of light intensities with
additional dark bottles to correct for non-uptake sorption of 14C at ambient
epilimnetic water temperature. The resultant photosynthesis-irradiance (P-I) data was
used to derive P-I curves by fitting a 3-parameter photosynthesis light-inhibition
model (Platt et al. 1980) to these data. The P-I curves were coupled with concurrent,
high-frequency photosynthetically active radiation (micromol m-2 s-1; PAR)
measurements and water column light extinction data (m-1) to estimate daily primary
production (mmol C m-3 d-1) in both Sparkling and Trout Lake. Over this time period,
the availability of data for 14C production varied due to sporadic sample
contamination and equipment failures.
Methods for 14C incubations in Acton Lake were similar to those in NTL lakes.
Integrated samples were collected from the euphotic zone (usually equal to the
epilimnion) and incubated in the lab for 1-2 hr with NaH14CO3 at a range of light
intensities (including dark bottles; Fee 1990). Incubations were usually done every
two weeks (23 of 55 experiments over the four years) or more frequently (24
experiments); only 8 experiments were done at intervals greater than 2 weeks. As in
NTL lakes, P-I curves were coupled with high-frequency PAR measurements, and water
column light extinction data collected at weekly intervals. Detailed methods 14C are
described in Knoll et al. (2003).
At Castle Lake, vertical water collections were made from 13 depths between the
surface and 30 m; duplicate light and 1 dark bottle samples from each each depth were
labeled with inorganic 14C in the form of NaHCO3 and then incubated in situ at the
depth of collection for 4 hours. Detailed methods are described elsewhere (Goldman et
al. 1963, Goldman 1968). Total daily incident solar radiation was measured throughout
the summer with a LI-COR Li-200 pyrheliometer. Light profiles at the height of the
solar day are measured using a Biospherical Instruments 2104P radiometer. Daily
phytoplankton productivity rates were calculated by dividing productivity measured
during the incubation period by the fraction of the total daily PAR received during
the incubation.
Free-water O2 Metabolism Methods
The same approach was used to estimate pelagic primary production (mmol C m-3 d-1;
GPP) in all lakes using in situ time series of dissolved oxygen data (O2). Free-water
O2 production estimates were based on high frequency measurements of dissolved oxygen
(mg L-1), water temperature ( ), PAR (micromol m-2 s-1), wind speed (m s-1), and
barometric pressure (mbar). Data frequencies varied from 1 to 15 minutes based on the
research program and the year data was collected. The raw, high-frequency time series
of dissolved oxygen and water temperature were filtered to remove outliers by
excluding values that were greater than 3 and 5 standard deviations respectively from
a 7-day running average (Appendix 1A,B; sensu Phillips 2020). The choice of sampling
frequencies has implications for the processes influencing dissolved oxygen patterns
and the amount of data needed to characterize those processes (Staehr et al. 2010). In
general, frequencies between 30 minutes and 3 hours are optimal for capturing changes
driven by biological processes (Staehr et al. 2010). Thus, we extracted hourly time
series for all high frequency data by averaging observations (mean value) on the hour
of observation (n = 4-60 depending on frequency of raw data) centered on the hour
(Phillips 2020) for use in metabolism models.
Epilimnetic depth (m) was quantified from either high-frequency thermistor string
data (Trout, Sparkling, and Acton Lakes) or discrete temperature profiles (Castle
Lake). The high frequency data was filtered for outliers as outlined above and
epilimnetic depth determined using the rLakeAnalyzer package (Read et al. 2011,
Winslow et al. 2019) at the temporal frequency of the raw data. Hourly aggregate data
was then extracted based on a 1-day running average to reduce the significant amount
of noise that existed in these estimates (Appendix 1C). rLakeAnalyzer was also used to
quantify epilimnetic depth from bi-monthly water temperature profile data in Castle
Lake and linearly interpolated at hourly time steps between observations.
Exchange of dissolved gas with the atmosphere is a critical component of
metabolism models, and, while there are a number of different models for estimating
piston velocities in lentic ecosystems (Dugan et al. 2016), the model proposed by
Vachon and Prairie 2013 is robust across multiple different types of lakes (Dugan et
al. 2016) and the metabolism model leveraged in this study (see below) is quite robust
to choice in parameterization of piston velocities (Phillips 2020). Piston velocities
(m hr-1) were calculated using the LakeMetabolizer R package (Winslow et al. 2016) and
the parameterization proposed by Vachon and Prairie 2013. Light extinction
coefficients (m-1), which were typically quantified bimonthly in all lakes, were
linearly interpolated at hourly time steps between observations, and combined with
epilimnetic depth and PAR to estimate the average light levels within the epilimnion
of each lake (Staehr et al. 2012a, Phillips 2020)
The data described above was used to generate daily estimates (mmol O2 m-3 d-1) of
gross primary production (GPP), respiration (R), and net ecosystem production (NEP)
using a time-varying ecosystem metabolism model (Phillips 2020). This model differs
from many of the more commonly used metabolism models (e.g., Winslow et al. 2016) in
that the model is not fit iteratively over a daily time scale, but rather
characterizes changes across all time periods (hourly measurements across 4-7 years of
data) for a given lake in a single model fit, as well as constraining GPP and R to
positive and negative values respectively (i.e., ecologically feasible ranges;
Phillips 2020). This takes advantage of the fact that the physical and biological
processes governing ecosystem metabolism and other aspects of DO dynamics are
autocorrelated through time, which means that this shared information can be used to
inform the parameter estimates across all time points. Furthermore, this method is
statistically unified because it uses all data to fit a single model, which
facilitates characterizing the uncertainty in the ecosystem metabolism estimates
(Phillips 2020).
The model used here differs slightly from that presented in Phillips 2020 in that
we used a photoinhibition P-I curve (Steele 1962) to describe GPP (sensu Staehr et al.
2016) instead of a light saturating curve: where PI is the production rate at light
intensity I, Pmax is the maximum production rate, and Iopt is the optimal light
intensity. This photoinhibition model was chosen because recent work by Staehr et al.
(2016) found that photoinhibition in lakes was common. The model by Steele (1962) is
one of the simplest photoinhibition models (two-parameter), and, regardless of the P-I
curve formulation chosen, it is often difficult to distinguish significant differences
in model fits between different models (Aalderink and Jovin 1997). Both Pmax and Iopt,
along with the model coefficient associated with R, (see Phillips 2020) were allowed
to vary through time at a daily time scale. The degree of auto-correlation in the
parameters through time was constrained by hierarchical variance parameters in the
random walk components of the model. Attempts to fit these parameters were
unsuccessful, which is unsurprising as hierarchical variances often have poor
identifiability. Thus, the random walk variances were treated as a tuning parameters
and were selected manually such that the model converged while producing meaningful
temporal smoothing in the parameters of the photoinhibition curve.
Observed dissolved oxygen time series were fit to all years (Trout: 2007-2010,
2012; Sparkling: 2007-2013; Castle: 2014-2017; Acton: 2010-2012, 2014) simultaneously
for each lake individually (i.e., lake-specific metabolism model fitting). Missing
values in the model input data time series lead to some days having fewer than 24
observations. Although the metabolism model is robust to missing data because it fits
the entire time series simultaneously instead of in discrete daily time steps, we did
not estimate metabolism parameters for an individual day if more than two hours of
data was missing for that day (Phillips 2020). The model was fit via Stan (GET
VERSION) run in R (GET VERSION) using the rstan package (STAN Citation) as described
in (Phillips 2020). Posterior median values were used for daily production values
along with the 0.025 and 0.975 quantiles of the posterior values to characterize the
95% credible intervals. Model fits were validated by checking effective sample size,
R, tree depth, energy Bayesian Fraction of Missing Information, and divergence (see
Betancourt 2007). Metabolism parameters were not estimated when the epilimnetic depth
was shallower than the dissolved oxygen sensor ( less than 0.5m in Trout, Sparkling,
Acton; less than 3m Castle; 4% of all observations). Gross primary production values
(mmol O2 m-3 d-1) were converted to units of carbon (mmol C m-3 d-1) assuming a
photosynthetic quotient (O2:CO2) of 1.25 (Bott 1996, Hanson et al. 2003,
Wielgat-Rychert et al. 2017).
