Coupled simulation of high‐frequency dynamics of dissolved oxygen and chlorophyll widens the scope of lake metabolism studies

High frequency time series of dissolved oxygen (DO), delayed chlorophyll fluorescence (Chl) and relevant background variables were recorded on nearly 900 d during 7 yr in large, shallow, meso‐eutrophic Lake Balaton (Hungary). Novel models were developed for coupled simulation of diel dynamics of DO and Chl using sequential learning and uncertainty assessment in a Bayesian framework. Despite the generally good model fit for both variables, the uncertainty of the metabolic estimates was high, due primarily to the identification problem of individual metabolic processes. Deviations between observed and simulated DO concentrations suggested that neglect of transient stratification might be responsible for the bulk of the systematic model errors. Net ecosystem production (NEP) was uncertain. Unless air‐water oxygen exchange can be estimated from direct measurements, the free‐water DO method cannot reliably estimate NEP. Gross primary production (GPP) could satisfactorily be hindcasted assuming non‐linear multiplicative dependence on Chl, water temperature and light. Hindcast of community respiration (CR) was less successful, possibly due to the impact of local benthic respiration. Results suggested a major shift in lake metabolism at about 16°C. Below and above this temperature, 70% and 90% of net primary production could be utilized by heterotrophs within a day, respectively. Indirect evidence suggested that biomass‐specific net primary production was determined by phosphorus. The large difference between reproductive rates and net growth rates estimated from GPP and Chl and from daily change in Chl, respectively indicated that loss rates of phytoplankton were as important determinants of algal dynamics as reproductive rates.

In a carbon-based biosphere, every single individual and each ecosystem take part in the global biogeochemical cycle of carbon (C). There is a considerable interest in understanding dynamics of C cycling in various types of lakes (also called "lake metabolism"), because it reflects overall ecosystem functioning (Odum 1956;Batt et al. 2013). Studies of lake metabolism have recently gained momentum in the context of climate change that brought into question whether lakes are sinks or sources of CO 2 at regional and global scales Tranvik et al. 2009). Since oxygenic photosynthesis and aerobic respiration dominate the biological carbon cycle (Holland and Turekian 2011), cycles of carbon and O 2 are inherently associated in most aquatic ecosystems. Freewater assessment of lake metabolism on the basis of dissolved O 2 concentration (DO) has a long tradition (Odum 1956;Vollenweider et al. 1974), because it is easier to measure DO than inorganic carbon, and because dissolution of O 2 in water is a pure physical process in contrast to the complicated chemical equilibrium of CO 2 (Stumm and Morgan 1981;Hanson et al. 2003). With the spread of relatively cheap, automated, reliable DO sensors, an increasing number of lake metabolism studies are based on the free-water DO method (Staehr et al. 2010a). The relationship between C and DO metabolism, however, is not straightforward.
An approximate stoichiometric coupling between CO 2 and O 2 exists only in oxygenic photosynthesis and aerobic respiration. Anoxygenic photosynthesis, anaerobic respiration and fermentation tend to decouple C and O 2 cycles in aquatic ecosystems. Photo-oxidation of organic carbon further modifies the C : O 2 stoichiometry of metabolism in humic lakes (Gran eli et al. 1996). Besides these complications, molar quotients of both photosynthesis and respiration (O 2 : C and C : O 2 , respectively) exhibit taxon-specific and physiological variability (Falkowski and Raven 1997). While inputs and outputs of DO other than gas exchange with the atmosphere are usually negligible, catchment-borne inputs of organic carbon may profoundly influence lake metabolism (Hanson et al. 2003;Tranvik et al. 2009;Staehr et al. 2010b). Odum (1956) has developed a simple method to estimate gross primary production (GPP) and community respiration (CR) from diel DO curves measured in streams: where rates are expressed as g O 2 m 23 d 21 . X is the atmospheric exchange rate (positive when O 2 escapes from the water). Inflow and outflow of DO via advection have been neglected. Local input from the watershed, IN l;w , has also been neglected because it was small and did not show diel variability. Odum (1956) solved his equation using the "bookkeeping" method -the only option when DO data are sparse (hourly or less frequent). Although nowadays high frequency DO measurements (order of minutes) are common, most studies adopted Odum's method to estimate metabolic rates (Cole et al. 2000;Hanson et al. 2003;Lopez-Archilla et al. 2004;Staehr et al. 2010a). Dynamic models with various levels of complexity have also been developed that allow assessment of the uncertainty in metabolic estimates (Hanson et al. 2008;Solomon et al. 2013;Honti et al. 2016).
Uncertainties that showed up in modeling were most often explained by neglected physical processes, primarily vertical and horizontal transport (Hanson et al. 2008;Staehr et al. 2010a;Solomon et al. 2013;Rose et al. 2014). Indeed, a significant effect of transport on single point DO measurements has repeatedly been demonstrated (Lauster et al. 2006;Sadro et al. 2011;Staehr et al. 2012b;Antenucci et al. 2013;McNair et al. 2015). Surprisingly, the major source of uncertainty, the identification problem has been overlooked in nearly every study of lake metabolism (but see Hanson et al. 2008;Honti et al. 2016), while it has attracted much attention in other disciplines, such as economics and hydrology (Manski 1993;Beven 2001). The identification problem means that quite different model structures and/or more than one set of parameters of the same model generate observationally equivalent distributions and therefore there is no way to decide which model is the "best" one and which parameter set is the "true" one. According to Eq. 1 (assuming IN l,w 5 0), observed changes in DO should first be split into two unknown components: net ecosystem production (NEP; NEP 5 GPP-CR) and atmospheric gas exchange (X). To cope with the identification problem, metabolic studies applied empirical or semi-empirical functions to estimate X (Odum 1956;Cole and Caraco 1998;Staehr et al. 2010a). The choice of X, however, predetermined NEP and error in X propagated directly into uncertainty in NEP (Vollenweider et al. 1974;Gelda and Effler 2002;Dugan et al. 2016). A second identification problem emerges when NEP is decomposed into GPP and CR. Since night-time GPP is zero, CR (more precisely, CR-X) can be recognized with relatively high confidence when transport does not interfere seriously (Hanson et al. 2008). During the day, however, GPP and CR could only be identified at the cost of poorly supported assumptions ("bookkeeping" method) or a set of poorly recognizable parameters (mechanistic models). Although the respiration of both autotrophic and heterotrophic organisms depends on organisms' physiological status and environmental conditions (Geider and Osborne1989;Pringault et al. 2009), the "bookkeeping" approach assumed that CR was identical in the dark and in the light (Staehr et al. 2010a). Dynamic models linked GPP and CR to one or more drivers (light, temperature, and algal biomass), usually without constraining any of the parameters by direct observations (for a review, see McNair et al. 2015).
A better identification of the underlying processes in a dynamic system can be achieved by developing tighter theories, and/or collecting richer data (Manski 1993). Nowadays phytoplankton and/or dissolved organic matter fluorescence is measured at high frequency in many lakes, simultaneously with DO (http://www.gleon.org; http://www. dkit.ie/networking-lake-observatories-europe). Combining these datasets within an appropriate modeling framework may partially solve the identification problem and promote understanding of details of lake metabolism that have not been explored by existing models.
The present study aimed at coupling dynamics of DO and phytoplankton by mechanistic models of various complexity using high frequency records from nearly 900 d in a 7-yr long period from large, shallow, meso-eutrophic Lake Balaton (surface area is 594 km 2 , mean depth is 3.1 m), Hungary. Our objective was to explore whether these complex models could provide a consistent and deeper insight into lake metabolism than models aiming at capturing only DO dynamics. We adopted the modeling approach of Honti et al. (2016), which is based on sequential learning and uncertainty assessment in a Bayesian framework. The long dataset allowed us to statistically analyze metabolic rates as a function of environmental drivers at various time scales.

Study site
In shallow, elongated Lake Balaton strong, steady north, north-west winds generate relatively closed large-scale water circulations in the four basins (Shanahan et al. 1986;Fig. 1). Large longitudinal and transverse seiches are characteristic of the lake (Muszkalay 1973;Shanahan et al. 1986). Wind-induced waves cause frequent sediment resuspension (Luettich et al. 1990;Istv anovics et al. 2004). The main inflow (the Zala River) drains about half (2750 km 2 ) of the total watershed into the smallest and shallowest Basin 1 (surface area is 38 km 2 ; mean depth is 2.3 m, maximum depth is 3 m, mean water residence time is 3 months), where we operate our monitoring station (Fig. 1).
Balaton is a calcareous lake (in Basin 1, mean pH is 8.4, alkalinity is 4 mEq L 21 , calcium concentration is 46 g Ca m 23 , magnesium concentration is 57 g Mg m 23 ).
The water is always turbid due to intense calcite precipitation and sediment resuspension. Secchi depth averages 0.5 m in Basin 1 during the ice-free period. Macrophyte growth is restricted to about 5% of surface area by low light availability and high wave exposure (Entz and Sebesty en 1946;Vir ag 1998;Istv anovics et al. 2008).
Rapid eutrophication during the 1970s (Herodek 1986) has successfully been managed by measures taken from the late 1980s. Although N 2 -fixing filamentous cyanobacteria still dominate summer phytoplankton in Basin 1 in most years (Padis ak et al. 2006), recovery from eutrophication has been observed from the late 1990s (Istv anovics and Somly ody 2001). In the early 2000s, the climate sensitivity of the lake has unexpectedly come into the focus of management. After several decades of maintaining water level in a gradually narrowing range prescribed on socio-economic rather than ecological grounds, the lake lost one third of its volume during the severe drought in 2000-2004. Climate change has been predicted to increase the frequency of low water levels (Honti and Somly ody 2009).

High frequency measurements
Our monitoring station was situated in the muddy northern littoral (43845 0 14 00 N and 17814 0 59 00 E; Fig. 1), where aquatic macrophytes had low abundance during the study period. Water depth at the station varied between 0.9 m and 1.9 m as estimated from a reference depth -water level and daily mean water levels of the lake (http://www.hydroinfo. hu). Measurements were carried out between 2009 and 2015, usually from early April to late October. Sensors were inspected and cleaned manually twice a week according to a maintenance protocol. We kept strict electronic record of every intervention and event that might concern the quality of measurements.
In April 2009, we deployed an electrochemical DO sensor (WTW Oxi, Germany) in the middle of the water column. This sensor was replaced with an optical sensor (Hach HQD, U.S.A.) in May 2009. From 2010, two optical DO sensors were positioned at known depths (0.3-0.6 m and 0.8-1.2 m from the bottom depending on water level). Sensors were calibrated in water-saturated air. DO concentration and water temperature were recorded every minute, taking a 30 s average.
Wind speed and direction, as well as global radiation (400-1100 nm) were recorded at 6.1 m from the bottom every minute (Mettech Ltd., Hungary). The height of the anemometer above the water surface was calculated from mean daily water depth. Wind speed was corrected for Istv anovics and Honti Coupled simulation of DO and Chl dynamics 10 m (U 10 , m s 21 ) according to the logarithmic rule (Smith 1988). Turbidity was recorded every 6 s with lightbackscattering sensors (Wetlabs, U.S.A.) facing downwards. The number of turbidity sensors decreased from 5 to 2 during the study period due to technical reasons. The concentration of chlorophyll (Chl) and photosynthetic properties of phytoplankton were estimated every 35-45 min by the delayed fluorescence (DF) technique, which detects only living, actively photosynthesizing algae (Gerhardt and Bodemer 1998;Istv anovics et al. 2005). The fully automated DF equipment (Tett Ltd., Hungary) sampled water from the middle of the water column into a dark chamber. After 5-10 min of dark adaptation, algae were illuminated with white LEDs at saturating light intensity. DF decay kinetics was measured in the dark for 100 s. To measure light-dependence of DF, dark-adapted algae were subjected to 9-18 light levels up to 400 lmol photon m 22 s 21 before their DF signal was measured for 1 s (dimension: Nominal DF Unit, NDFU; Honti and Istv anovics 2011). The hyperbolic tangent function of Jassby and Platt (1976) was fitted to the data using Solver in Microsoft Excel to estimate a B (NDFU [mg Chl] 21 [lmol photon m 22 s 21 ] 21 ), the biomass-specific initial slope that reflects efficiency of light utilization and P B max (NDFU [mg Chl] 21 ), the biomass-specific maximal rate of photosynthesis.

Manual measurements
To convert vertically averaged turbidity (NTU) to the coefficient of diffuse light attenuation (K d , m 21 ), underwater light profiles were measured manually each week at 0.1 m depth intervals using a 4p quantum sensor (LiCor, U.S.A.). K d and mean turbidity were linearly related (K d 5 0:078 NTU 1 0:612, r 2 5 0.89, n 5 92).
The concentration of Chl was estimated from DF data using conversion factors obtained from linear regressions between DF kinetic integral and weekly manual measurements of Chl. Water samples were filtered onto Whatman GF/F glass fiber filters. Chlorophyll a was extracted for 24 h in cold acetone in the dark, and measured with a TD700 fluorimeter (Turner Designs, U.S.A.) using the non-acidification optical kit (EPA Method 445). Conversion factors were different from year to year because of major upgrades of the DF equipment. The determination coefficient of linear regressions exceeded 0.85 in each year (n 5 20-25).
Potential phosphorus (P) limitation of GPP was assessed using cell quota of surplus phosphorus (Q SP , mg P [g C] 21 ). Since surplus P has not been measured in the present study, we used weekly data obtained between 2000 and 2004 at our monitoring site (Istv anovics et al. 2004). Triplicate water samples were filtered onto pre-washed cellulose acetate membrane filters (0.2 lm pore diameter; Whatman) and extracted in boiling distilled water for 1 h (Fitzgerald and Nelsson 1966). Concentration of soluble reactive P (SRP) was determined in the extracts according to Murphy and Riley (1962). Carbon-based cell quota was estimated by assuming a C/Chl ratio of 50.
Daily mean discharge (Q, m 3 s 21 ), daily loads of SRP (molybdenum blue method), total P (measured as SRP after H 2 SO 4 -H 2 O 2 digestion), nitrate (the salicylate method) and weekly load of dissolved organic C (DOC; CHN autoanalyser) were obtained from the West Transdanubian Water Directorate at the inlet gauge of the Zala River that transports about 90% of total loads to Basin 1.

Metabolic models
The most parsimonious model for simulating coupled dynamics of DO and phytoplankton biomass requires three state variables: DO, autotrophic biomass (B a , mg Chl m 23 ) and heterotrophic biomass (B h , g C m 23 ): where the superscript i denotes instantaneous rates. Each rate is expressed as g O 2 m 23 d 21 . Community respiration was split into its autotrophic (R i a ) and heterotrophic (R i h ) components (CR i 5 R i a 1 R i h ). The quotient of photosynthesis and the reciprocal of the quotient of autotrophic respiration were assumed to be equal [q O 2 : Chl ð Þ ] and were expressed in terms of chlorophyll (g O 2 [mg Chl] 21 ). Respiratory quotient of heterotrophic organisms q C : O 2 ð Þ had the dimension g C (g O 2 ) 21 . GHP i was gross heterotrophic production. A metabolic model cannot fully reflect the diversity of trophic interactions. Equation 2 was based on two simplifying assumptions. First, the contribution of higher than secondary trophic levels to community metabolism was negligible and thus, heterotrophic organisms included only bacteria, grazers, and detritivores. Second, net sedimentation of phytoplankton was zero.
To address the lack of high frequency information on B h and GHP i , two solutions were possible. The first option was to further simplify Eq. 2 by assuming that sub-diel changes in B h were negligible (Model 1; Table 1, Supporting Information Table S1). This assumption did not preclude day-to-day and seasonal variability in B h , since the rate of heterotrophic respiration at 208C (R 20 h ) was calibrated and its changes should reflect changes in B h at scales over a day. The second option was to simulate daily change in heterotrophic biomass (Model 2; Supporting Information Table S1). While Model 1 decreased the power of the identification problem compared to Odum's model (Eq. 1), Model 2 provided a more detailed description of lake metabolism than the models published so far on free-water DO dynamics (McNair et al. 2015).
GPP i was a hyperbolic tangent function of mean photosynthetically active radiation (PAR) in the mixed water column (PAR i , lmol photon m 22 s 21 ; Jassby and Platt 1976) with a van't Hoff-Arrhenius type dependence on water temperature (T i , 8C; Talling 1966; Supporting Information Table S1). Respiration was also temperature-dependent. GHP i was proportionate to B a and independent of B h . We chose this description because in Lake Balaton the most important consumers of algae and algal detritus were chironomid larvae besides bacteria (Specziar and V€ or€ os 2001). Since maximum clearance rate of filter feeding zooplankton did not exceed 10% of water volume per day (G.-T oth et al. 1986), grazing was omitted. X i was a linear function of saturation deficit and piston velocity (Supporting Information Table S1). DO i sat was approximated from water temperature by a 5 th order polynomial function (Garcia and Gordon 1992). The Schmidt number (Sc i ) was calculated from T i according to Wanninkhof (1992). Piston velocity (k i 600 , m d 21 at Sc 5 600) should depend on both wind speed and convective currents (MacIntyre et al. 2010;Dugan et al. 2016). Since no direct measurements were available on atmospheric gas exchange in Lake Balaton to choose the proper empirical function for k i 600 , we used an iterative approach. First we ran Model 1 on the whole data set with k 600 as an adjustable parameter. Thereafter we examined the relationship between estimated k i 600 and wind speed during periods of cooling, warming and constant temperature. The k i 600 5 f(U i 10 ) relations were similar in these three periods. Therefore, the final lake-specific function used in both Models 1 and 2 was based on wind speed only:

Processing of input variables
The measured time series were preprocessed before using them as inputs to the models. A day was defined as the period between two successive sunrises. We selected days when each input variable, as well as DO and Chl data were available during the whole day. The number of such days was 887. Four additional days were also accepted when DO measurements ceased during the night, but at least half of the dark period was covered by data. Of the input variables, T i , U i 10 , incident PAR (PAR i 0 ) and K i d were averaged for 10 min. T i was also averaged vertically. DF data (P B max and a B ) were averaged for a day and converted from their original dimensions to their dimensions used in the model. Based on previous results (Honti and Istv anovics 2011), a conversion factor of 0.004 mg O 2 NDFU 21 d 21 was adopted. Considering that mean measured P B max was in the order of 10 5 NDFU [mg Chl] 21 , this corresponded to 0.4 g O 2 [mg Chl] 21 d 21 , a realistic value (Stefan and Fang 1994). PAR i was estimated according to Reynolds (1997): We assumed that PAR was 47% of global radiation (J m 22 s 21 ) and mean energy of photons was 218 kJ mol 21 in the visible range (400-700 nm). Mixing depth (Z mix ) was taken equal to the depth of the water column in spite of frequent but transient stratification that usually lasted for a few hours (V€ or€ os et al. 2010).

Parameters and model calibration
The number of adjustable parameters was reduced before fitting the models by finding constant values for certain parameters. We assumed that the molar respiratory quotient of heterotrophs was 1, that is q C : O 2 ð Þ5 0.375 g C (g O 2 ) 21 . Similarly, molar quotients of both photosynthesis and autotrophic respiration were assumed to equal 1. Since these latter quotients were expressed in terms of Chl, a C/Chl ratio was needed to find the value of q O 2 : Chl ð Þ . Although C/Chl ratios of phytoplankton might vary in a wide range depending primarily on resource availability and temperature (Geider et al. 1998), we found a good agreement between phytoplankton biomass and daily mean DO amplitude (A DO , g O 2 m 23 ) assuming a C/Chl ratio of 50 mg C (mg Chl) 21 (Supporting Information Fig. S1). A DO was defined as the difference between the 99 th and 1 st percentiles of DO concentrations measured during a day. Using this C/Chl ratio, q O 2 : Chl ð Þequaled 0.133 g O 2 (mg Chl) 21 . To determine the temperature coefficients of photosynthesis and respiration (h P and h R , respectively) we ran Model 1 on a sample of days with different parameter values. With "wrong" temperature coefficients, P B; 20 max and CR 20 strongly co-varied with T. The final choice was h P 51.06 and h R 51.09. With these values, Spearman's rank correlation coefficients between daily mean T and T-dependent parameters of both models (P B; 20 max , CR 20 and R 20 h ) were low (0 > r s > 20.1; n 5 891). Daily respiration of algae was taken as 10% of autotrophic biomass (f 20 a 50.1 d 21 ; Reynolds 1997). In the observed range of daily mean water temperatures (10.6-31.68C), 4-27% of biomass could be respired in a day. Daily R a represented 4 to over 100% of daily GPP (average was 11%). The final number of adjustable model parameters was three in Model 1 (P B; 20 max , a B , CR 20 ) and four in Model 2 (P B; 20 max , a B , R 20 h , k GHP ; cf. Supporting Information Table S1). Models were calibrated to DO data that were averaged both vertically and for 10 min time steps and to daily mean Chl using the method of Honti et al. (2016). Calibration was done in a sliding window of 3 d. In each calibration unit, Bayesian parameter inference and uncertainty analysis were performed. Measured daily mean values of P B max and a B with a prescribed coefficient of variation (CV; CV 5 standard deviation/mean) of 20% were taken as prior distributions in each calibration unit. Non-measured parameters (CR 20 , R 20 h and k GHP ) were estimated by a sequential Bayesian learning procedure, in which posterior distributions obtained in the preceding day became prior distributions in the actual day (for the philosophy of this approach, see Honti et al. 2016). Structural errors of the metabolic model were described by a first-order autoregressive error model for DO (Reichert and Schuwirth 2012). A single likelihood function incorporated fit to both the diel DO curve and daily mean algal biomass (Honti et al. 2016). Posterior distribution of parameters was sampled for each calibration window with Markov chain Monte Carlo sampling using the Metropolis algorithm (Gamerman 1997). The differential equations of the model (Eq. 2, Supporting Information Table S1) were solved with the LSODA solver (Hindmarsh 1983;Petzold 1983). Instantaneous rates were integrated over the day; in the following only daily integrals are presented and analyzed.

Statistical methods
A locally weighted regression/smoothing technique (LOWESS; Cleveland 1979) was used to visualize non-linear relationships between daily rates of metabolism and their drivers. Multivariate non-linear regression models were fitted to square root transformed area-specific daily GPP and CR (Y, g O 2 m 22 d 21 ): where Y max was the maximum transformed daily rate, j denoted driver variables and Q j f j reflected the assumption that the effect of various drivers (f j ) was independent and multiplicative. For the sake of simplicity, each f j was expressed by the same function: where a5 x crit 2x In Eq. 7, b 5 5 was used to force f j x ð Þ to attain a near-zero value at x min . Regression was calculated by using a conventional likelihood function assuming independent, normally distributed errors. Square root transformation of data before model fitting was needed to improve the normality of model residuals.
Autocorrelations of variables were measured by Pearson product-moment correlation coefficients, whereas correlations between variables were characterized by Spearman's rank correlation coefficients. Due to the large number of data (n > 730), correlations and linear regressions were highly significant even at diminutive r. Instead of using statistical significance, we considered a correlation or regression useful if the resolution power index (N class ) of Prairie (1996) exceeded a threshold value of 2.

Raw data
There was a considerable year-to-year variability in both Chl and external nutrient loads to Basin 1 (Supporting Information Fig. S2). Although N 2 -fixing cyanobacteria occasionally made up more than half of phytoplankton biomass in most years (L. V€ or€ os, unpubl.; Supporting Information Fig.  S2), large summer blooms developed only in 2010 and 2015. The highest Chl concentrations coincided with low external nutrient loads. Concentration of DOC was relatively constant in the Zala River and thus, DOC load (L DOC ; g C m 22 month 21 ) was primarily dependent on discharge.
To characterize the typical diel course of DO and its drivers, we pooled data from each day (n 5 891) and calculated the median and 95% uncertainty range (2.5-97.5% quantiles) of variables every 10 min (Fig. 2). Median DO ranged between 8.1 and 9.7 g O 2 m 23 during this "average" day. As expected, peak DO coincided with maximum T and was delayed behind the highest PAR by 4 h. DO saturation ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 600=Sc p ) had its minimum before dawn and rose up to about 1.1 d 21 by sunset (Fig. 3), following the course of median water temperature and wind speed (Figs. 2, 3). Thus, gas exchange tended to be higher during the day than during the night (Fig. 3). Consequently, an error in X might more strongly influence the estimates of GPP than CR. During the "average" day, there was a net oxygen release of 64 mg O 2 m 23 d 21 from the lake.
Of the environmental drivers of metabolism, water temperature was strongly autocorrelated (Table 2). Autocorrelations of, and correlations between other drivers were too weak (r < 0.75) to reach the resolution power of N class 2.

Simulation results
There was a good agreement between simulated and observed time series of Chl (Model 1: Chl sim 5 5:34 6 0:19 ð Þ 1 0:91 6 0:01 ð Þ 3Chl obs , r 2 5 0.94, n 5 891, N class 5 5.24; Model 2: Chl sim 5 2:74 6 0:15 ð Þ 1 0:94 6 0:01 ð Þ 3Chl obs , r 2 5 0.96, n 5 891, N class 5 6.48). In general, both models fitted reasonably to the DO time series. Measurement noise was an order of magnitude higher in the first 18 d of the study when DO was measured with an electrochemical sensor than later on when we used optical sensors. In the latter period, the standard error of fit ranged between 0.05 g O 2 m 23 and 0.62 g O 2 m 23 . Relative error was obtained by dividing standard error with the amplitude of diel DO change. Relative error remained below 10% in 59 and 66% of days fitting Models 1 and 2, respectively (Fig. 4). We exemplified a particularly good, an average and a bad fit of Model 2 by selecting days when frequency of relative error was 0.025, 0.5, and 0.975 (Supporting Information Fig. S3). Since the two models performed similarly well and yielded similar estimates for shared metabolic rates, only results obtained with Model 2 will be presented in detail.
To assess the influence of neglected transport, we determined the distribution of the DO error (DDO 5 DO obs -DO sim ) for every 10 min time step (n 5 128,172). Mean wind speed and mean temperature difference between the two DO sensors (DT, 8C m 21 ) were calculated for each bin of DDO (Fig.  4b,c). In the range, where 100 or more DDO values were available, the error tended to increase with decreasing wind speed and increasing strength of stratification. The same pattern emerged by integrating the area between observed and modeled DO curves (not shown). However, integration resulted in a relatively small sample size (n 8500) due to strong autocorrelation of deviations. Considering the "average" day, the peak in DDO coincided with maximal DT (Fig. 4d).
Despite their wide ranges, GPP and CR showed relatively smooth seasonal variations (Supporting Information Fig. S4). This resulted in high autocorrelation of GPP and CR, that partly reflected the autocorrelation of A DO and partly that of the main drivers (Table 2). Net primary production (NPP, NPP 5 GPP -R a ; [g O 2 m 23 d 21 ]) exceeded net heterotrophic production (NHP, NHP 5 GHP-R h ; g O 2 m 23 d 21 ) by an order of magnitude (Supporting Information Fig. S4). Uncertainties of metabolic estimates were compared by their mean daily CVs calculated from posterior parameter distributions. Both model structure and the success of recognizing metabolic processes influenced the CVs (Table 3). Community respiration could be identified with the highest confidence, whereas uncertainty of net ecosystem production exceeded the most likely value by up to several orders of magnitude. Although Model 2 had one more parameter than Model 1, it generally yielded less certain estimates for metabolic processes than Model 1 because of the higher level of the identification problem.
There was a strong linear relationship between A DO and GPP that slightly exceeded the N class 5 2 threshold, whereas the relationship between A DO and CR had weaker resolution power (Fig. 5a). This was reasonable, since A DO approximately equaled the sum of GPP and half of CR. As expected, GPP and CR as well as NPP and GHP were linearly related (Fig. 5). According to Model 1, 79% of GPP was respired within 1 d (N class 5 2.88). Model 2 resulted in a somewhat lower CR relative to GPP (74.5%, N class 5 2.38). Of net primary production, 92% could be channeled to gross heterotrophic production in the day of production (Fig. 4c; N class 5 4.07). Heterotrophic respiration contributed by 6.7-99% to CR, its mean share was 88%. R h peaked at an intermediate phytoplankton biomass of 20-40 mg Chl m 23 . After an initial increase with heterotrophic biomass, GHP leveled off at B h 5-10 g C m 23 . Most of the variability in NEP (88%) was explained by gas exchange (Fig. 5f).

Drivers of metabolism beyond sub-diel scales
We related GPP and CR with daily mean values of several potential drivers including wind speed, significant wave height, temperature, DT, buoyancy frequency, etc. using LOWESS smoothing for visualization. Highly scattered, nonlinear patterns emerged with algal biomass, temperature and PAR as independent variables (Supporting Information Fig.  S5). Multivariate non-linear regressions fitted with these drivers captured 80% and 58% of variability in GPP and CR, respectively. Temperature limitation caused most of the variability in both GPP and CR. GPP was limited moderately by PAR (Fig. 6).
Biomass-specific NPP declined with increasing algal biomass (Fig. 7a). A similar pattern was observed in the phosphorus limitation factor (f P ) that was calculated from cell quota of surplus P according to the Droop model (Droop   1973; f P 5 Q SP 2Q SP;0 À Á =Q SP ), assuming a maintenance cell quota, Q SP;0 of 6 mg P [g C] 21 . Daily NPP was averaged for weeks. Weekly mean NPP ranged between 0.05 g C m 22 d 21 and 4.8 g C m 22 d 21 during the study period (Fig. 7b). Its seasonal trend was remarkably similar from early May until late July (weeks 18-30) in each year, followed by rapid differentiation across years. The lowest late summer productivity was observed in 2014 when summer was relatively cold and cyanobacteria failed to dominate (cf. Supporting Information Fig. S2). The highest productivity occurred in the hot summers of 2010 and 2015 when cyanobacteria made up 80-94% of phytoplankton biomass (Supporting Information Fig.  S2; L. V€ or€ os, unpubl.). Using the exponential growth equation, the reproductive rate of phytoplankton (l, d 21 ) was approximated as l 5 ln(1 1 NPP/B a ). Mean weekly l was as high as 0.76 d 21 during the study period (Fig. 7c). At the same time, 96% of net growth rates calculated from daily change in phytoplankton biomass were between 20.5 d 21 and 10.5 d 21 (not shown).
Weekly, monthly and seasonal mean NPP were unrelated to external nutrient loads. R h was not related to DOC load. Monthly mean NPP exceeded the respective L DOC by a mean factor of 25 (the range was from 1.3 to 139). Gross heterotrophic production tended to increase with the increasing NPP/ L DOC ratio (Fig. 7d).

Discussion
Alternative metabolic models -is there a superior model?
A few studies have compared performance of various metabolic models over short periods. Assuming that transport could be neglected, both Hanson et al. (2008) and McNair et al. (2015) found that the simplest DO model (Odum 1956) fitted nearly as well to the data as models that incorporated more biological complexity when the basic assumption on transport was met and vice versa, even complex models failed when the impact of transport was significant. In line with these studies, our more complex Model 2 yielded marginally better fits than Model 1, and both models performed poorly in certain days (Fig. 4, Supporting Information Fig. S3). Deviations between observed and modeled DO tended to increase with decreasing wind speed and increasing strength of stratification (Fig. 4), indicating that the lack of full mixing was the primary source of systematic error. Strong radiative forcing could build up relatively stable thermal stratification in calm days; the temperature difference between surface and bottom waters occasionally reached 58C (V€ or€ os et al. 2010). Sea and land breezes were strong enough to result in a periodic daily wind pattern on the background of much stronger synoptic winds (Fig. 2) but they could not prevent development of transient stratification, the mean duration of which was 5.1 h per day (V€ or€ os et al. 2010). Coupling hydrodynamical and metabolic models might reduce systematic error (Staehr et al. 2012b;Antenucci et al. 2013;McNair et al. 2015). Further research, however, is needed before this coupling will be generally feasible. In Lake Balaton, hydrodynamic modeling has a long tradition (Muszkalay 1973;Shanahan et al. 1986;Kr amer et al. 2004). Lake-specific 1D-3D models reasonably capture water movements (basin-wide circulations, seiches, windinduced waves) but they lack the required accuracy to be meaningfully coupled with a local heat budget (Torma and Kr amer 2017) or DO model.
Goodness of fit to observed data is an important albeit not an absolute measure of usefulness of a model. In contrast to classical calibration methods, where the single objective is to minimize fit error, the Bayesian parameter inference makes a formal statistical compromise in fit to match a priori expectations that reflect independent knowledge on the examined system. We had three basic expectations. First, recognizing that a metabolic model without a transport term was structurally imperfect (Lauster et al. 2006;Van de Bogert et al. 2007;Hanson et al. 2008;Sadro et al. 2011;Staehr et al. 2012b;Antenucci et al. 2013;McNair et al. 2015), we assumed that model errors were significantly autocorrelated. Although 60% of modeled DO values (n 5 128,172) deviated by no more than 0.2 g O 2 m 23 from measured DO, autocorrelation of errors occasionally resulted in large departures from measured DO curves and in high parameter uncertainty (Fig. 4, Supporting Information Fig. S4; Table 3). Second, we expected that day-to-day variability in community respiration and photosynthetic parameters should be smooth in a resilient ecosystem. Applying sequential Bayesian learning and model calibration in a 3 d sliding window (Honti et al. 2016), no irreversible shifts occurred in estimated metabolic rates during the 7 yr long study period to falsify the resilience hypothesis. This could be seen from smooth seasonal changes (Supporting Information Fig. S4) and high autocorrelation of estimated metabolic rates (Table 2). Third, acknowledging that productivity and biomass of phytoplankton were closely coupled, we expected that diel DO curves and biomass dynamics of algae could simultaneously be modeled when daily biomass data were available for calibration. This could successfully be done both in this study and in other lakes (Honti et al. 2016).
We agree with Hanson et al. (2008) that multiple metabolic models can be equivalent in terms of fit quality, making the simplest good performing model the best for practical applications. At the same time, too simple models hide a series of scientifically important information. Heuristic values, such as insight into the consistency of metabolic estimates do distinguish between models that fit similarly well to the data. In published models, mass balance of DO is open, giving no chance to check consistency of model outputs. In contrast, both of our models showed that metabolic rates were consistent with observed changes in phytoplankton biomass. Additionally, the biologically more complex Model 2 allowed us to discriminate between autotrophic and heterotrophic production and respiration (Supporting Information Table S1). In this abstract sense, Model 2 was superior to Model 1 and to published metabolic models. Complexity of models should only be increased if additional information was available to keep the identification problem under control. Developing Model 2, we intended to include grazing into GHP besides bacterial decomposition and detritivory. Grazing should depend on both autotrophic and heterotrophic biomass, whereas the supply of exudates and algal detritus can be thought of as a function of only B a (Supporting Information Table S1). In Lake Balaton, where filter feeding zooplankton were insignificant (G. -T oth et al. 1986), the impact of grazing could not be recognized from DO and Chl data. In lakes where zooplankton are important, grazing may leave its signature in the DO and Chl time series and may potentially be identified in a metabolic model. In scientific applications, complexity of metabolic models should be driven by data availability and basic limnological knowledge on the study lake.

Uncertainty of metabolic estimates -can NEP reliably be estimated from diel DO curves?
In spite of the reasonable fit of our models (Fig. 3), uncertainty of metabolic estimates was high for several reasons (Table 3; Supporting Information Fig. S4). Uncertainty arising from neglect of transient stratification burdened each metabolic rate, as explained by Coloso et al. (2011). In most cases, Model 1 yielded rate estimates with lower CV than Model 2 because of better identification of processes. For shared rates, CVs followed the order CR < GPP < NEP (Table  3). Unlike NEP, CR and particularly GPP were strongly associated with, and therefore constrained by the diel amplitude of DO (A DO ; Table 2; Fig. 5).
Community respiration could be identified with the highest confidence, because nighttime decrease in DO depended primarily on CR and thus, about one third of the observations supported the identification of this single process. Uncertainty was introduced by gas exchange and transport, the latter being most evident when the DO signal increased during the night (also observed by Hanson et al. 2008). GPP was more uncertain than CR, because its magnitude was partly dependent on uncertain CR and because it was more strongly affected by gas exchange than respiration (Table 2; Fig. 2). Moreover, similar uncertainty of GPP in the two models suggested that measurement error might also contribute: time series of mean underwater PAR might occasionally be inconsistent with measured photosynthetic parameters (P B max and a B ) that were used in our models as prior parameter distributions. Uncertainty of NEP was disappointingly large considering either its CV (Table 3) or standard deviation (in Model 2, mean SD was 3.7 g O 2 m 23 d 21 , range was from 0.1 g O 2 m 23 d 21 to 11.8 g O 2 m 23 d 21 ). According to Model 1, NEP averaged at 20.03 g O 2 m 23 d 21 during the study period, whereas in Model 2, the balance was slightly positive (0.075 g O 2 m 23 d 21 ). Similarly, Cremona et al. (2014) observed that various models resulted in either net heterotrophy or in net autotrophy in Lake Võrtsj€ arv (Estonia). In agreement with our modeling results, low organic content of sediments (1-4%; M at e 1987) indicated that primary production was closely balanced by mineralization at long time scales (years, decades) in Lake Balaton. NEP is a small difference between two large, uncertain metabolic rates (GPP and CR) that does not have its own constraining drivers (Eq. 1, Supporting Information Table S1), therefore its calculated value is incidental.
Besides its high uncertainty, NEP strongly depended on the choice of the gas exchange function in a recent ensemble modeling study (Dugan et al. 2016). This is not surprising, since from a decomposition perspective, NEP equals the difference between observed changes in DO concentration and X (Eq. 1). Accordingly, in the present study X explained 88% of the variability in NEP (N class 5 3.74, n 5 891; Fig. 5). Proper description of gas exchange is crucial to reliably estimate NEP and to assess the role of lakes as sources or sinks of carbon. In the present study, we aimed at iteratively deriving a lake-specific function for atmospheric gas exchange coefficient from diel DO curves. Since atmospheric exchange rates have not been measured in Lake Balaton to verify our k f estimates, we cannot offer this approach as a well-tried recipe. Nevertheless, it seems worth developing independently verified methods to estimate gas flux from high frequency DO curves in lakes of different sizes, morphometry, and wind exposure.
In spite of the uncertain and incidental nature of NEP, comparative studies have revealed meaningful patterns in NEP along broad gradients of nutrient and DOC availability (Hanson et al. 2003;Staehr et al. 2010b;Solomon et al. 2013). The possible explanation for the emergence of these patterns is that dependence of NEP on gas flux cannot fully mask the large differences in diel changes in DO saturation and concentration among lakes that cover a spectrum of trophic types and states. This, however, does not make the absolute values of NEP estimates less incidental or more credible. Consequently, ranges of TP and DOC concentrations where lakes are predicted to be net autotrophic or heterotrophic need revision when these ranges have been derived from high frequency DO measurements.

Ecological observations beyond the autotrophicheterotrophic dichotomy
The interannual stability of weekly NPP between late spring and late July (Fig. 7c, Week 15-30) was consistent with the observation of Honti et al. (2007) that in the lack of a sufficiently long-lasting unidirectional selective pressure, phytoplankton were in a "ground state" characterized by low biomass and hectic compositional changes. In contrast, variability of weekly NPP from late July reflected year-to-year variability of growth conditions in blooming sub-spaces of phytoplankton identified by Honti et al. (2007).
Reproductive rates of algae estimated from NPP often reached 1-1.2 d 21 , the average was 0.76 d 21 during the study period (Fig. 7). These rates reasonably agreed with reproductive rates obtained under natural conditions using different methods (Schnoor and Di Toro 1980;Reynolds and Irish 1997). Net growth rates estimated from daily change in B a (20.5 d 21 to 10.5 d 21 ) also agreed with commonly reported values (Sommer 1981;Reynolds and Irish 1997;Istv anovics et al. 2005). The large difference in the rates of reproduction and net growth indicated substantial mortalities of phytoplankton implying that algal dynamics depended as much on loss processes (exudation of excess photosynthates, sedimentation, grazing, physiological mortality) than on growth (Schnoor and Di Toro 1980;Crumpton and Wetzel 1982;Reynolds and Irish 1997). Strong linear regression between NPP and GHP (Fig. 5) suggested that in Lake Balaton heterotrophs utilized on average, 90% of NPP in the same day when biomass was produced.
There was no relationship between external nutrient loads and phytoplankton biomass or metabolic rates in this study (Supporting Information Figs. S2,S4). Similar exponential decrease in the upper bounds of both biomass-specific NPP and P limitation factor with increasing algal biomass (Fig. 7) was an indirect indication that growth rates might be increasingly P-determined (f P < 0.5) when biomass exceeded 2 g C m 23 (40 mg Chl m 23 ). These observations were supported by findings that P desorption from resuspended sediments was the prime source of P supply of seriously phosphorus deficient summer phytoplankton in Lake Balaton (Lijklema et al. 1986;Istv anovics and Herodek 1995;Istv anovics et al. 2004).
When both NPP and NHP were positive (503 d out of 891), NHP averaged at 30% of NPP. This value agreed well with the estimate of Cole et al. (1988) who found that bacterial and zooplankton production made up a mean 20% and 12% of NPP, respectively in a broad selection of fresh and saltwater ecosystems. Since high concentrations of finegrained particles make Lake Balaton into a stressful habitat for filter-feeding zooplankton (G. -T oth et al. 1986), the contribution of these two groups of heterotrophs to NHP might differ from the estimate of Cole et al. (1988).
Daily area-specific GPP could reasonably be hindcasted by multivariate non-linear regression using daily mean areaspecific biomass of phytoplankton, T, and PAR as independent variables (Fig. 6, Supporting Information Fig. S5). Hindcast of area-specific CR was considerably weaker using either B a , and T (Fig. 6, Supporting Information Fig. S5) or B a 1B h , and T. The pragmatic consequence was that CR could be predicted with a roughly similar uncertainty from either B a and T or from the multivariate non-linear regression model of GPP combined with the linear regression between GPP and CR (cf. Fig. 5, Supporting Information Fig. S5). Our explanatory variables overlapped with the most powerful predictors of GPP and CR in other productive lakes (Harris 1978;Staehr and Sand-Jensen 2007;Staehr et al. 2010b).
Of the predictor variables, water temperature explained the largest variability in both GPP and CR. The temperature limitation factor increased more steeply with T for CR than for GPP (Fig. 6) as also noted by Staehr et al. (2010b). The f T values of these to processes were equal at about 168C, suggesting a major shift in lake metabolism below and above this temperature. Indeed, splitting the dataset at T 5 16.08C, parameters of the linear regression between volumetric GPP and CR in the warm season (CR5 0:38 6 0:08 ð Þ 1 0:91 6 0:02 ð Þ 3GPP, r 2 5 0.67, n 5 811, N class 5 2.29) became nearly identical to those of the NPP vs. GHP relationship (Fig. 5). Although the regression between GPP and CR did not reach the resolution power index of 2 in the cold season (r 2 5 0.39, n 5 80), the slope was significantly lower (0.71 6 0.10) whereas the intercept was significantly higher (0.59 6 0.14) than in the summer. Thus, about 90% and 70% of GPP might immediately be respired in the summer and spring/autumn, respectively, and a higher fraction of CR might be supported by allochthonous DOC when productivity was low in the cold season. Since we usually failed to capture the spring diatom bloom (Supporting Information Fig. S2) that might have developed prior to the start of our monitoring program, phytoplankton biomass was low in the 80 d colder than 168C. This left open the question, how much increased availability of autochthonous carbon sources might trigger community respiration at low temperature.
The temperature limitation factors of GPP and CR attained a value of f T 5 0.99 at 308C and 268C, respectively. The high temperature, at which GPP was fully released from T-limitation might be explained by the high optimum temperature for growth of Cylindrospermopsis raciborskii (around 318C; Shafik et al. 2001;Briand et al. 2004;Kov acs et al. 2016), the dominant summer species in years when cyanobacteria blooms have developed in Lake Balaton (Padis ak and Reynolds 1998). During the study period, 17% and 2% of days were warmer than 268C and 308C, respectively. In a changing climate, both frequencies are expected to increase. A disproportionately larger increase in the frequency of days above 268C than 308C might potentially increase gross productivity and reduce net productivity in Lake Balaton, increasing the chance of net C release. This prediction may hold true in other lakes dominated by cyanobacteria.
In contrast to the considerable increase in daily GPP with T and autotrophic biomass, the impact of PAR was modest (Fig. 6). Moreover, GPP showed an unexpected increase when underwater light exceeded about 120-130 lmol photon m 22 s 21 (Supporting Information Fig. S5). This was surprising because GPP i was a hyperbolic tangent function of PAR i in our models (Supporting Information Table S1). The lack of a strong light limitation might be due to the low light requirement of species selected by the generally poor underwater light climate that prevailed in Lake Balaton (Reynolds 1997). PAR tended to exceed 125 lmol photon m 22 s 21 for several days only in July and August, during the potential blooming period of C. raciborskii. This cyanobacterium has low light requirement (I k is 15-30 lmol photon m 22 s 21 ) and it is tolerant to light intensities up to several hundreds of lmol photon m 22 s 21 (Shafik et al. 2001;Briand et al. 2004;Kov acs et al. 2016). The coincidence of high temperature, physical stability and relatively good light conditions might enhance photosynthesis and growth of C. raciborskii and this could lead to the unexpected increase in GPP in the range of the highest PAR values (Supporting Information Fig. S5). Benthic primary production might also contribute to this upward trend in GPP (see below). Several studies have demonstrated a significant contribution of benthic/littoral metabolism in both seasonally stratified and shallow lakes, even though availability of benthic habitats was limited in some of these lakes (Lauster et al. 2006;Van de Bogert et al. 2007;Sadro et al. 2011;Idrizaj et al. 2016). A rich and abundant algal flora inhabits the sediments of Lake Balaton (Uherkovich and Lantos 1987). Typically, Fragilaria construens dominate the phytobenthos (V€ or€ os et al. 2000). Benthic primary production was comparable under clear ice to pelagic production in summer in mesotrophic Basin 4 of the lake (Herodek and Ol ah 1973). P B max of benthic algae collected at shallow (< 1 m) sites in summer varied between 10 and 20 mg C [mg Chl] 21 h 21 in laboratory incubations ( € Uveges et al. 2011). Assuming a 10 h photoperiod, these rates translated into 270-540 mg O 2 [mg Chl] 21 d 21 . For comparison, biomass-specific GPP averaged at 175 mg O 2 [mg Chl] 21 d 21 and ranged between 0 and 400 mg O 2 [mg Chl] 21 d 21 in our study. Thus, lightsaturated benthic photosynthesis had a high potential to influence measured DO time series. I k values of phytobenthos were estimated to range between 240 lmol photon m 22 s 21 and 940 lmol photon m 22 s 21 ( € Uveges et al. 2011). The maximum, 97.5 th percentile and median 10 min light intensities that reached the sediment surface were 222 lmol photon m 22 s 21 , 40 lmol photon m 22 s 21 , and 3.4 lmol photon m 22 s 21 , respectively, at our monitoring site during the study period (n 5 78,424). Thus, local benthic photosynthesis was unlikely to influence DO concentrations significantly even after taking into consideration that I k estimates of phytobenthos by € Uveges et al. (2011) might be too high. (Samples were illuminated with daylight fluorescent tubes, the emission spectrum of which is different from the underwater light spectrum in a turbid lake.) Halving the depth range observed at our site (0.9-1.9 m), the lowest laboratoryderived I k value of 240 lmol photon m 22 s 21 occurred at the bottom for a shorter or longer period nearly each day. In this way horizontal currents associated with transverse seiches might transport oxygen-rich water masses from the vicinity of the shore to the monitoring site about 100 m offshore. Conclusively, the mechanism how benthic photosynthesis might influence offshore DO concentrations might be the same in shallow, turbid lakes as in deep ones (Lauster et al. 2006;Van de Bogert et al. 2007;Sadro et al. 2011). Unlike local benthic photosynthesis, local benthic respiration might significantly influence our CR estimates. This may partly explain that the deviation between observed and modeled 10 min DO values tended to increase with transient stratification (Fig. 4), and could be one of the reasons why the multivariate hindcast of CR was less successful than that of GPP (Fig. 6).

Metabolism in Lake Balaton
Weekly area-specific net primary production averaged at 1.5 g C m 22 d 21 during the study period and typically varied between 2 and 4 g C m 22 d 21 during the summer (Fig. 7). These values were similar to those estimated at around the turn of the century Istv anovics et al. 2004) and lower by a factor of 2-3 than those estimated during the eutrophication in the 1970s (Herodek 1986) using different variants of the 14 C technique. Compared to other lakes where metabolic rates have been estimated from high frequency DO measurements, Basin 1 of Lake Balaton ranked among the most productive third of systems (Laas et al. 2012;Staehr et al. 2012a;Solomon et al. 2013).
Monthly mean NPP exceeded monthly external load of DOC by a mean factor of 25. Considering that biological availability of DOC ranged from 9% to 15% at the mouth of the Zala River irrespective of the season (T oth et al. 2007), bacterial production should have been fueled almost exclusively by NPP. Nevertheless, the tendency of increasing GHP with NPP/L DOC ratios (Fig. 7) suggested that microbial utilization of allochthonous DOC might occasionally be important in Basin 1 of Lake Balaton. The intercept of the GPP vs. CR regression (Fig. 5) suggested that a daily mean respiration of 0.73 g O 2 m 23 d 21 might be supported by DOC.

Implications for high frequency lake metabolism studies
The role of lakes as sinks or sources or CO 2 is a focal question in most contemporary studies of lake metabolism. We argue that lake metabolism studies that use the free-water DO technique cannot reliably estimate net ecosystem production unless independent lake-specific information is available on gas exchange between the air and water.
Transient stratification and full mixing may cause as much uncertainty in metabolic estimates in shallow lakes as in deep ones. Therefore, vertical DO profiles may bear important information on mixing processes in shallow lakes, too.
Complex models of lake metabolism that incorporate relevant biological, chemical, and/or physical data (e.g., algal biomass, DOC concentration, velocity field) may help to assess consistency of metabolic estimates and may reveal details of lake metabolism that are beyond the scope of simple models.