Diel, seasonal, and inter-annual variation in carbon dioxide effluxes from lakes and reservoirs

Accounting for temporal changes in carbon dioxide (CO2) effluxes from freshwaters remains a challenge for global and regional carbon budgets. Here, we synthesize 171 site-months of flux measurements of CO2 based on the eddy covariance method from 13 lakes and reservoirs in the Northern Hemisphere, and quantify dynamics at multiple temporal scales. We found pronounced sub-annual variability in CO2 flux at all sites. By accounting for diel variation, only 11% of site-months were net daily sinks of CO2. Annual CO2 emissions had an average of 25% (range 3%–58%) interannual variation. Similar to studies on streams, nighttime emissions regularly exceeded daytime emissions. Biophysical regulations of CO2 flux variability were delineated through mutual information analysis. Sample analysis of CO2 fluxes indicate the importance of continuous measurements. Better characterization of short- and long-term variability is necessary to understand and improve detection of temporal changes of CO2 fluxes in response to natural and anthropogenic drivers. Our results indicate that existing global lake carbon budgets relying primarily on daytime measurements yield underestimates of net emissions.


Introduction
The global carbon budget is rapidly changing in response to anthropogenic emissions (Hanson et al 2006, Friedlingstein et al 2020. Prior studies have estimated that 0.14-0.64 Pg C-CO 2 is annually released to the atmosphere through lakes and reservoirs (Cole et al 1994, 2007, Aufdenkampe et al 2011, Ciais et al 2013, Raymond et al 2013, Holgerson and Raymond 2016, DelSontro et al 2018, Drake et al 2018. However, most of these estimates are made with relatively limited sampling, generally constrained to the open-water or summer season during the daytime, and with limited consideration of interannual and shorter-scale variation (Butman et al 2018, Ran et al 2021. Underrepresentation of temporal changes and variability of carbon dioxide (CO 2 ) flux in existing CO 2 flux inventories may bias estimates of lake CO 2 emissions (Deemer et al 2016, Klaus et al 2019. Recent studies have found nighttime emissions exceeding daytime emissions because of potential uptake in lakes (Shao et al 2015), reservoirs (Liu et al 2016), and rivers (Gómez-Gener et al 2021). A lack of frequent and long-term CO 2 observations also limits our ability to differentiate natural CO 2 flux variations from the consequences of anthropogenic perturbations (Hasler et al 2016). Multiyearscale time series that capture sub-annual variability of the aquatic CO 2 flux remain rare (Huotari et al 2011, Shao et al 2015, Finlay et al 2019. Traditional in-situ aquatic sampling methods for CO 2 concentrations and fluxes in natural and artificial freshwaters also come with high uncertainty (Golub et al 2017, Baldocchi et al 2020, with one source being the heterogeneity of littoral and pelagic lake CO 2 fluxes (Erkkilä et al 2018, Spafford andRisk 2018).
Advances in the past several decades, however, have enabled more long-term, continuous highfrequency (hourly) measurements of CO 2 flux in freshwater ecosystems, which are capable of capturing the dynamics of air-water fluxes at time scales of hours to years (Eugster et al 2003, Huotari et  At these time scales, CO 2 fluxes have been shown to respond to variations in photosynthesis and respiration rates (Cole et al 2007), wind speed and direction (Podgrajsek et al 2015), carbonate equilibria (Atilla et al 2011), ecosystem metabolism (Provenzale et al 2018), convective mixing (Eugster et al 2003, Mammarella et al 2015, internal waves (Heiskanen et al 2014), ice phenology (Reed et al 2018), chlorophyll a concentration (Shao et al 2015), atmospheric turbulence generated by surrounding topography (Eugster et al 2022), and hydrological and carbon inflows/outflows Kortelainen 2005, Weyhenmeyer et al 2015). These sources of variation may be overlooked by low-frequency and seasonrestricted sampling that dominate freshwater science (Desai et al 2015).
A growing number of previous studies were conducted using eddy covariance (EC) flux towers, which measure ecosystem-scale air-water CO 2 fluxes (Vesala et al 2006) and has gained prominence for use in freshwaters (Vesala et al 2012). While its application over lakes has mostly covered short periods of time (e.g. Eugster et al 2003, Vesala et al 2006, Podgrajsek et al 2015, an increasing number of sites are now measuring lake-atmosphere fluxes continuously over multiple years (Huotari et al 2011, Mammarella et al 2015, Shao et al 2015, Franz et al 2016, Reed et al 2018, Eugster et al 2020. Other methods for high frequency sampling have also included the use of forced diffusion automated chambers (Spafford andRisk 2018, Rudberg et al 2021). Here, to identify modes of CO 2 flux variability missed by infrequent sampling from lakes and reservoirs, we quantify diel to inter-annual dynamics of CO 2 fluxes directly measured by EC towers in 13 lakes and reservoirs in the Northern Hemisphere, the first synthesis of this kind. Mutual information analysis is utilized to determine sources of temporal CO 2 flux variability at each site. Additionally, this dataset contains sites with a wide range of nutrient status (i.e. eutrophic, oligotrophic, and mesotrophic), therefore allowing for comparisons between temporal changes of CO 2 fluxes among different freshwater systems.

Study sites
Data on air-water CO 2 exchange and meteorological drivers were acquired from 19 study sites across the Northern Hemisphere, with at least one season of observations between 2005-2015, of which 13 were retained here for analysis (tables 1 and S1). The six remaining submitted sites were withheld for challenges in meeting uncertainty and gap filling criteria (see supplemental methods). These sites were collected based on organization of a workshop (Desai et al 2015) and an open call through listservs. Selected sites included 9 lakes and 4 reservoirs, mostly located between 40 • N and 68 • N latitude, coinciding with the largest area of Earth covered with lakes. Eight sites had data available over multiple seasons, but only a few also had measurements during winter ice cover. Lake area ranged from 0.036 km 2 to 623 km 2 (median: 15.2 km 2 ), with median mean depth of 6 m (range: 0.6-11 m); most developed a seasonal thermocline and were dimictic or monomictic (table S1). Two water bodies had a significant fraction of submerged and emergent macrophytes (SE-Tam and DE-Zrk) within the footprint of the flux tower.

Measurements
The EC technique directly measures the exchange of momentum, heat and matter (water vapor, CO 2 , or other trace gasses) at the air-water interface and is a reliable method for measuring surface exchanges with the atmosphere (Vesala et al 2006). The flux towers were located on floating platforms, lake shoals or islands, or on shore depending on the site (Supplemental Material Text; table S1). The towers were additionally equipped with instruments providing half-hourly to hourly measurements of biophysical variables (e.g. net radiation (Rnet), air temperature (T air ) and humidity, photosynthetically-active radiation (PAR), 3D wind direction and speed, water surface temperature (T water ), aquatic CO 2 or O 2 concentration, and water level), although data availability and frequency varied among the sites. Data were harmonized to uniform formats and units, screened for fetch, de-spiked, and gap-filled using a common flux post-processing standard prior to calculation of diel and monthly averages (Pastorello et al 2020;Supporting Material Text). Flux footprint and wind directional screening were applied by each investigator to retain only those half-hourly observations representative of the water body. Note that a negative CO 2 flux indicates uptake by the ecosystem from the atmosphere and a positive flux means the reverse. All data are published in the Environmental Data Initiative repository (Golub et al 2022).

Data analysis
We analyzed the half-hourly CO 2 fluxes and three major groups of biophysical covariates. The first group included variables related to wind forcing acting on the water surface (i.e. friction velocity, wind speed, momentum flux). The second group encompassed variables related to temperature cycles and proxies of energy in the system (i.e. T air , T water , ∆T (T water − T air ), sensible (H) and latent heat (LE) fluxes). The last group included the variables associated with solar radiation-proxies for primary productivity (i.e. ∆pCO 2 (pCO 2water − pCO 2air ), PAR). To determine the standardized difference between two means with repeated unpaired measurements and imbalanced population sizes, we used the Cohen's d test where the mean difference between the mean daily CO 2 fluxes is divided by the pooled variance. A coefficient d of 0.20, 0.50, 0.80 indicates small, medium, and large standardized differences between the two means, respectively.
To determine the degree of the net ecosystem exchange (NEE) predictability by biophysical drivers (i.e. T air , T water , H and LE, friction velocity (Ustar), and Rnet), we performed mutual information analysis (MI). MI measures the statistical dependence of a variable 'Y' on another variable 'X' . Using the marginal and joint probability distributions of 'X' and 'Y' , it expresses the proportion of bits needed to represent 'Y' that is redundant given the knowledge of 'X' . There are also no parametric assumptions regarding the relationship between the two variables, making MI capable of discerning linear and non-linear associations (Knox et al 2021). More information on MI may be found in Fraser and Swinney (1986). Ultimately, this method can reveal dependencies between two variables with co-varying factors, making it a useful approach for ascertaining NEE dependencies on ecosystem variables (Knox et al 2021). We use mutual information scores (MIS) to determine the relative strength of each bivariate interaction. High and low MIS corresponds to a strong and weak link between NEE and biophysical variables, respectively.
To account for driver impacts on different temporal scales, we used the maximal-overlap discrete wavelet transform (MODWT) approach to decompose half-hourly data into four temporal scales (hourly, diel, multiday, and seasonal). MODWT involves applying a high-pass wavelet filter and lowpass scaling filter to the time series. This allows a decomposition of the time series into multiple fine and coarse scales, enabling analysis at varying timescales (Percival 1995, Percival andWalden 2000). We used the 'Wavelet Methods for Time Series Analysis' in MATLAB's 'Wavelet Toolkit' to decompose the data (Cornish et al 2003). Further details on the application of this method may be found in the supplement and Sturtevant et al (2016).
Finally, to quantify potential sample bias from infrequent sampling and identify optimal sampling approaches, a random sample analysis was conducted on the gap-filled oligotrophic (US-Too; 2012), mesotrophic (FI-Van; 2016), and eutrophic (DE-Zrk; Table 1. Comparison of ice-free CO2 flux at temporal (i.e. annual, seasonal, diurnal and nocturnal) scales derived from high-frequency eddy covariance measurements over lakes and reservoirs. CO2 fluxes are presented as the mean ± one standard deviation. The numbers in brackets represent the number of observations integrated at a given time scale.  2014) lake/reservoir data with the smallest data gaps. One thousand random samples without replacement were taken for each of the following times: daytimeonly (DT), daytime/nighttime-only (DT/NT), summer mid-day (SMD), growing season (GS), and annual. DT and NT were defined as 10:00-15:30 h and 22:00-03:30 h (local times), respectively. Hours between 11:00-13:30 h were considered mid-day while the GS counted fluxes between March 1st and September 30th. Flux counts of 1, 5, and 10 were taken at each temporal resolution to test dependencies of percent error (PE) on flux count within and between sites. To obtain a single flux value, the counts containing 5 and 10 fluxes were averaged. This sampling algorithm was created using Python version 3.8.3.

Magnitude of CO 2 fluxes
Study sites represented a wide range of nutrientcolor status and physical characteristics of water bodies, and as a result spanned a range of daily CO 2 fluxes, although with some common elements (figure 1). The mean daily CO 2 flux across all sites was 0.43 ± 0.34 µmol CO 2 m −2 s −1 (mean ± standard deviation, (SD); range: −0.075-1.25 µmol CO 2 m −2 s −1 ) with only 6% of observations indicating neutral fluxes or net CO 2 uptake. The spread of time-resolved fluxes varied between 102% and 798% of the site-specific daily mean (figure 1). Reservoirs had smaller but more variable fluxes relative to the lakes (0.32 ± 0.71 vs. 0.41 ± 0.31 µmol CO 2 m −2 s −1 , respectively), although the reservoir sample size is smaller and more geographically restricted. Two thirds of sites had at least 66% of daily fluxes within the cross-site flux mean ±1 SD (Cohen's d: 0.02 < d < 0.76).
Annually, all sites were CO 2 sources to the atmosphere, except for DE-Zrk and LA-NT2, with large variability across sites (figure 2). This was also the case when comparing the same lake or reservoir type. On a single day (on average), the mesotrophic and eutrophic lakes and reservoirs were the largest and smallest carbon sources, respectively. While most sites were a greater carbon source during the nighttime relative to the daytime, the difference in hourly fluxes was small (0.5 µmol C m −2 s −1 ), except for DE-Zrk.

Temporal variability
Averaged diel CO 2 changes had regular patterns of daytime minima and nighttime maxima across all sites in most months ( figure 2(a)). Daytime hourly fluxes were on average 35% (range 7%-60%) lower than nighttime fluxes, though in 94% of site-days, those were still net positive emissions. Despite the commonly observed daytime CO 2 flux dip, the flux decrease was large enough to convert our sites to daily net sinks of CO 2 in only 11% of site-months. The mean uncertainty of diel CO 2 was strongly influenced by extreme observations, with 192% mean uncertainty, but only 79% median uncertainty.
Maximum diel flux amplitudes occurred typically in July and August and ranged 0.24-1.09 µmol CO 2 m −2 s −1 . Relative to the summer amplitudes, shoulder season CO 2 flux amplitudes were on average 44%-49% smaller in May and September and 26%-37% in April and October. Monthly to sub-annual CO 2 flux variability was nearly twofold higher compared to the diel flux variation. Surprisingly, we found frequent sub-monthly (20-30 d) variability across all water bodies, regardless of the system's physical or biogeochemical conditions. While most site-level variability fluctuated around the CO 2 flux averages, for some, amplitudes scaled with flux minima and maxima (figure S1). Sites with multi-year data had relatively consistent sub-annual patterns across years, although the timing and amplitudes of submonthly variability varied among lake-years. When integrated over time-resolved daily CO 2 fluxes, both sub-monthly and sub-annual modes of variability accounted for two thirds of the site-level daily variability (range 10%-190%). Mean and median uncertainty were 167% and 67% of mean daily CO 2 flux, respectively.
Once scaled to ice-free season annual emissions, and assuming zero fluxes during ice cover, we found all water bodies were net sources of CO 2 , despite missing any ice off/on related fluxes (table 1). The cross-site mean and SD of 23 site-years was 95 ± 49 gC m −2 yr −1 (range: 14-224 gC m −2 yr −1 ). Inter-annual variability (IAV) was calculated as a SD of annual CO 2 flux for each site with multiyear data (figure S2). The mean cross-site IAV was 22 gC m −2 yr −1 (25%) and ranged between 4 and 44 gC m −2 yr −1 (3%-58%).

Drivers of CO 2 fluxes
While the continuous data allowed capturing CO 2 flux variability at multiple temporal scales, we still had a limited capacity to attribute which factors and processes governed the CO 2 flux. We found small standardized differences between CO 2 fluxes among site groups belonging to the three humic states (d < 0.01), medium differences between oligotrophic and eutrophic states (d = 0.24), and large CO 2 differences between mesotrophic and oligotrophic states (d = 0.66), and between mesotrophic and eutrophic states (d = 0.72). Commonly observed biophysical covariates explained an average of 32% of variance in half-hourly CO 2 fluxes. Wind-related variables were identified as a key to explaining CO 2 flux variability in eight out of 13 sites. Biophysical variables related to exchanges of heat at the airwater interface, particularly ∆T and turbulent energy exchange (H and LE), both correlated with CO 2 flux. The fitted regressions were non-linear and highly variable across sites, owing to ecosystem differences and presence of confounding factors (e.g. differential responses to co-dependent covariates).
Mutual information analysis revealed different drivers to be responsible for CO 2 fluxes on different temporal scales (figure 4). On hourly scales, CO 2 flux at all sites was predicted mostly by T air and T water . The strongest links were found to occur at LA-NT2 and DE-Zrk, both being eutrophic systems. Analysis on diel scales yielded a similar result. On multi-day scales, however, more linkage between CO 2 flux and drivers was found at CA-Dar, SE-Mer, and FI-Pal (all oligotrophic). While the seasonal scale MI analysis was subject to many gaps, it did show a more uniform CO 2 flux prediction magnitude across all sites and drivers relative to other timescales.

Temporal analysis
Random sampling among different temporal resolutions resulted in large differences between mean sampled NEE and mean continuous annual NEE ( figure S3). For DE-Zrk and FI-Van, the greatest PE was for samples taken during SMD, calculated to be 868% ± 26% and 38% ± 2% (mean ± range), respectively. US-Too experienced the largest error during NT sampling, with a PE of 87% ± 31%. Increasing the number of NEE values per sample (i.e. going from 1 to 5 to 10 samples with the latter two NEE values calculated as the average) gave sporadic results, in that, agreement sometimes improved (FI-Van during GS) and sometimes worsened (US-Too during nighttime). DT/NT and annual sampling were the most representative of continuous annual NEE among all sites regardless of lake/reservoir type. GS sampling showed PE that was well within the typical uncertainty for EC flux measurements (∼20%) for FI-Van and US-Too. Sampling on an annual scale further constrained PE, including even DE-Zrk in addition to FI-Van and US-Too.
Bias was also calculated using the same random sampling method. This was defined as subtracting the mean continuous annual NEE from the mean sampled NEE at each site. Overall, DE-Zrk showed high bias across all temporal resolutions (−0.8 ± 5.2 µmol CO 2 m −2 s −1 ) while US-Too and FI-Van had bias that was an order of magnitude smaller (figure 3).

Forcing differences by type
Through this synthesis, T water is shown to be a major predictor of lake and reservoir NEE, which is consistent with past work of Zwart et al (2019) and Eugster et al (2020). There is a high degree of spatiotemporal variability between these two variables. For example, NEE at LA-NT2 and DE-Zrk (eutrophic reservoir and eutrophic shallow lake respectively) was most highly predicted by T water on short timescales (hourly and diel), suggesting these ecosystems may be most susceptible to evading CO 2 following climate warming. This large link may also be explicable through lake type. Eutrophic lakes are defined as being nutrient rich, meaning they contain larger phosphorus, nitrogen, or dissolved organic carbon concentrations (DOC) than their oligotrophic counterparts (Reed et al 2018). On multiday timescales, however, the distinguishability of the NEE/T water linkage in the eutrophic sites is less prominent when compared to the oligotrophic and mesotrophic sites. Another variable with high NEE predictability was T air , though it is possible that in some cases fluxes have an indirect relationship with T air via its impact on DOC input from land (Sobek et al 2005). Due to lack of significant findings, we excluded asynchronous driver analysis in this work. We acknowledge that these drivers are present in other studies covering freshwater ecosystems such as in Sturtevant et al (2016).
Additionally, we acknowledge that there may have been CO 2 flux variabilities due to lithological or geographical variables not discussed here. For example, López et al (2011) found CO 2 dynamics to be weakly coupled to calcareous systems. León-Palmero et al (2020) noted reservoir CO 2 fluxes to be dependent on . Sample analysis for the mesotrophic (blue), eutrophic (green), and oligotrophic (red) lakes and reservoirs with the least data-gaps. Each bar shows the bias between mean sampled NEE and mean continuous annual NEE (mean continuous annual NEE subtracted from mean sampled NEE). Samples were taken without replacement. A zoomed in version of the plot is shown to better distinguish differences between FI-Van and US-Too. lithology, with a CO 2 source and sink being reported for calcareous and siliceous watersheds, respectively. The work also found an impact of reservoir mean depth on methane fluxes, which highlights that these influences are not limited to CO 2 . Water chemistry and carbon isotope data would also provide great insights on the contributions of dissolved inorganic carbon (DIC) transformation and organic carbon mineralization (Zhong et al 2018), especially in cases of damming where DIC transport is affected (Wang et al 2020). These processes have been documented to play an important role in CO 2 outgassing, but the limitations of our dataset prevented us from conducting further analysis on these processes. As a result, future studies must take these characteristics and processes into account due to their prominent role in influencing lake and reservoir greenhouse gas flux magnitudes across multiple temporal scales.

Unresolved temporal variation in CO 2 fluxes
CO 2 fluxes from lakes and reservoirs exhibited large variability at diel to inter-annual scales, which could comprise unresolved sources of uncertainty or bias in current estimates of annual CO 2 fluxes from infrequent and season-restricted sampling. Although our study lakes were not randomly selected and cannot be directly used to upscale (Stanley et al 2019), they were broadly reflective of common mid-latitude freshwater systems spanning a broad range of humic-status and mixing regimes. Additional considerations for measurements across lake size and catchment area (Hanson et al 2007, Holgerson andRaymond 2016) and hydrological setting  would be required to design a representative estimate for global upscaling.
We were able to investigate, however, the role of temporal variation on a range of systems that broadly reflect many lakes and reservoirs. Our reported continuous daily fluxes corresponded to the upper end (88th percentile) of previously published flux magnitudes (table S2). The observed temporal variation suggests that temporal restrictions in sampling may add a significant source of underestimation bias in existing inventories of CO 2 fluxes from lakes and reservoirs of similar type and size (Klaus et al 2019).
In particular, we found large diel variation in all study sites, with routinely higher emissions at night, consistent with a recent study over rivers (Gómez-Gener et al 2021). Diel reduction of dissolved CO 2 concentrations and fluxes is often associated with ecosystem metabolism (Hanson et al 2003) and supported by negative correlations with The vertical axis lists each site colored by nutrient level (red, blue, and green indicate oligotrophic, mesotrophic, and eutrophic lakes and reservoirs respectively), while the x-axis shows drivers listed from left to right as air temperature, water temperature, latent heating, sensible heating, friction velocity, and net radiation. High mutual information score (MIS) is represented by dark blue (large NEE predictability) and low MIS is represented by white and light blue (low NEE predictability). PAR, or even algal blooms (Shao et al 2015, Ouyang et al 2017. Water temperature (Provenzale et al 2018), carbonate equilibria fluctuations (Atilla et al 2011), water-side convection (Eugster et al 2003, Mammarella et al 2015, Podgrajsek et al 2015, and internal waves (Heiskanen et al 2014) can additionally govern diel CO 2 dynamics. Our observed diel amplitudes were within 21%-43% of sub-hourly flux amplitudes derived from dissolved CO 2 concentrations (Hanson et al 2003, Morales-Pineda et al 2014 or previously published EC-measured fluxes (Vesala et al 2006, Liu et al 2016. Our results support the conclusion that existing global lake carbon budgets that rely primarily on daytime measurements are underestimates of net emissions. We also found common sub-monthly modes of CO 2 flux variability across all of our sites. Similar variability in the continuous observations have been reported for dissolved CO 2 (Huotari et al, 2009 (Franz et al 2016, Eugster et al 2020, indicating the prevalence of oscillatory patterns in CO 2 time series at both sides of the air-water interface. Variability has been previously attributed to the interplay of wind forcing (Liu et al 2016), upwellings of CO 2rich waters (Morales-Pineda et al 2014), biologicallydriven (metabolic and trophic) changes in carbonate equilibria (Atilla et al 2011, Ouyang et al 2017, convective mixing (Eugster et al 2003, Huotari et al 2009, non-local processes (Esters et al 2020), and T water (Atilla et al 2011. However, this is the first study to find a consistent pattern in a wide range of systems, regardless of size. We also observed changes to the prevalence of underlying submonthly CO 2 flux variability through the year at several sites, likely reflecting seasonal ecosystem changes, such as spring/fall turnover (Baehr and DeGrandpre 2004), radiative and heat exchanges (Heiskanen et al 2014), and hydrological inflows (Vachon et al 2017).

Implications for the global carbon budget
After our daily fluxes were scaled to site-specific annual CO 2 emission fluxes, estimates were in the upper end of previously reported estimates for lakes and reservoirs (table S2). All systems were sources of CO 2 in most years, although there have been sites that reported significant carbon sinks (e.g. Shao et al 2015, Reed et al 2018, and additional propagation of uncertainty from data gap filling and filtering (e.g. of nighttime uptake) can weakly push some of our study sites toward sinks. While our lakes are not fully representative of all lakes on Earth, we postulate that improved temporal resolution of site-level CO 2 fluxes is one of the sources of differences between this study and published annual fluxes (table S2). The results also imply that a proposed recommended number of samples per year (4-8) (Natchimuthu et al 2017, Klaus et al 2019 is likely insufficient to constrain annual CO 2 fluxes from lakes and reservoirs. Rather, our sampling analysis suggests to increase nighttime and open-water season observations, preferably weekly, which would reliably increase the accuracy of annual estimates, given our observed diel and submonthly variations.
Additionally, sites with multiple years of data all showed non-trivial interannual variation. The estimate of average IAV of CO 2 fluxes (25%) is modest compared to that (88%) observed in terrestrial ecosystems (Baldocchi et al 2018), partially reflecting the lower number and diversity of ecosystems with multiyear measurements or more buffering against climate extremes by large water bodies. However, given that CO 2 flux from freshwaters positively scales with the productivity of terrestrial ecosystems at shorter timescales (Walter et al 2021, Hastie et al 2018, it is possible that the interannual variation of carbon input from land will propagate onto CO 2 evaded through freshwaters (McDonald et al 2013, Drake et al 2018, providing a possible pathway to better characterize freshwater IAV. Neglecting this variation is an additional source of bias in our current view on global CO 2 emissions from lakes and reservoirs. Given that the CO 2 fluxes are affected at both sides of the air-water interface (Wanninkhof et al 2009), a better measure of the contribution of lakes to the global carbon cycle will also require reporting and synthesis of additional continuous waterside data (e.g. temperature, dissolved CO 2 and O 2 ), site-level ecosystem characteristics (e.g. nutrientcolor legacies, ecosystem metabolism, and aquatic vegetation such as algae), surrounding topography (Eugster et al 2022), and sampling an increased site diversity within climatic zones (Lehner and Döll 2004). With more frequent air and aquatic observations, we will better constrain CO 2 fluxes at different time scales, assess the prevalence of temporal patterns in CO 2 fluxes, and reduce uncertainty in eddy flux measurements over freshwaters (e.g. Ejarque et al 2021) and therefore improve model estimates of responses of these ecosystems to climate change. These estimates will also allow for a robust representation of different climate variable effects on these fluxes (Sobek et al 2005). Such work will be needed to quantify and evaluate landscape (Buffam et al 2011 to global (DelSontro et al 2018) carbon budget components from lakes and reservoirs.

Conclusions
Across the 13 study sites with EC flux observations, all lakes and reservoirs were, on average, net annual sources of CO 2 to the atmosphere. However, the time series revealed large diel to (sub)-monthly CO 2 flux variability among the sites that represent a broad range of biogeochemical and physical site characteristics. These modes of variability accounted for two thirds of daily and a quarter of annual CO 2 flux variation, with sub-annual variability dominating over diel and inter-annual flux variabilities. After integrating these modes of variability into time-resolved fluxes, the CO 2 flux estimates were at the upper end of published CO 2 emissions for lakes and reservoirs. Our results support the idea that long-term, frequent measurements during both day and night of carbon dynamics in freshwater aquatic systems are critical to resolve lake NEE magnitudes and detect long-term trends of lake carbon fluxes. Omitting these temporal scales will not only limit our knowledge of lake NEE, but also restrict our understanding of biophysical driver impacts. We advocate for establishing and maintaining a long-term observation network that combines EC flux measurements with highly detailed site-specific carbon budget studies over key lake and reservoir ecosystems representing broader geographical gradients.

Data availability statement
The data that support the findings of this study are openly available at the following URL/DOI: 10.6073/ pasta/87a35ca843d8739d75882520c724e99e. the support from the ACCC Flagship funded by the Academy of Finland (337549 and 337552) and ICOS-Finland by University of Helsinki and the Ministry of Transport and Communication. GAW was financially supported by the Swedish Research Council (VR: Grant Nos. 2016-04153 and 2020-03222). The deployment of the EC at the Nam Theun 2 Reservoir (Lao PDR) was funded by Electricité de France (EDF) and Nam Theun Power Company (NTPC). TS and FK were supported by the Helmholtz Association of German Research Centres through grants to TS (Grant VH-NG-821) and FK (Grant PD-129), the Helmholtz Climate Initiative REKLIM (Regional Climate Change), and infrastructure funding through the Terrestrial Environmental Observatories Network (TER-ENO). TS and FK thank the (staff of the) Department Chemical Analytics and Biogeochemistry at the Leibniz-Institute of Freshwater Ecology and Inland Fisheries (Berlin) for providing water chemistry data for DE-Zrk. WE and GK acknowledge support from NSF-DEB 1637459 and OPP 1936769. The deployment of the flux tower at CA-Eastmain was supported by Hydro Quebec. Flux observations at US-OWC were funded by the Ohio Department of Natural Resources, by NOAA's National Estuarine Research Reserves' Davidson Fellowship, and by US Department of Energy awards DE-SC0021067 and DE-SC0022191. The co-authors express gratitude for the kindness and contributions of posthumous co-author Werner Eugster.

Open research
We have deposited all EC lake observations and gapfilled values in the Environmental Data Initiative repository Golub et al (2022). Several sites are also accessible from FLUXNET affiliated archives as noted in table S2.

Author contribution statement
M G designed experimental protocol and conducted the data syntheses. N K -A. conducted additional analyses and revisions. A R D, N K -A., and M G wrote the manuscript. T V, I M, G B, and G W supervised research, contributed observations, and edited the manuscript. All other authors contributed with flux observations and commented on the manuscript.