Contribution of gas concentration and transfer velocity to CO2 flux variability in northern lakes

The CO2 flux ( FCO2 ) from lakes to the atmosphere is a large component of the global carbon cycle and depends on the air–water CO2 concentration gradient (ΔCO2) and the gas transfer velocity (k). Both ΔCO2 and k can vary on multiple timescales and understanding their contributions to FCO2 is important for explaining variability in fluxes and developing optimal sampling designs. We measured FCO2 and ΔCO2 and derived k for one full ice‐free period in 18 lakes using floating chambers and estimated the contributions of ΔCO2 and k to FCO2 variability. Generally, k contributed more than ΔCO2 to short‐term (1–9 d) FCO2 variability. With increased temporal period, the contribution of k to FCO2 variability decreased, and in some lakes resulted in ΔCO2 contributing more than k to FCO2 variability over the full ice‐free period. Increased contribution of ΔCO2 to FCO2 variability over time occurred across all lakes but was most apparent in large‐volume southern‐boreal lakes and in deeper (> 2 m) parts of lakes, whereas k was linked to FCO2 variability in shallow waters. Accordingly, knowing the variability of both k and ΔCO2 over time and space is needed for accurate modeling of FCO2 from these variables. We conclude that priority in FCO2 assessments should be given to direct measurements of FCO2 at multiple sites when possible, or otherwise from spatially distributed measurements of ΔCO2 combined with k‐models that incorporate spatial variability of lake thermal structure and meteorology.

Lakes occupy only a small portion of the Earth's surface but are estimated to emit large amounts of carbon dioxide (CO 2 ) to the atmosphere with recent emission estimates ranging from 0.3 to 0.6 Pg carbon (C) per year (Holgerson and Raymond 2016;Raymond et al. 2013).The CO 2 emissions (F CO 2 ) from lakes depend on the gas transfer velocity (k) and differences in concentration of CO 2 in the surface water and its theoretical concentration when in equilibrium with the atmosphere (ΔCO 2 ): The variable k in Eq. ( 1) parameterizes the physical gas transport across the surface boundary layer and is modified by physical factors that influence turbulence, such as wind speed, direction and fetch, heating and heat loss, and currents (Banerjee and MacIntyre 2004;Read et al. 2012;Tedford et al. 2014;MacIntyre et al. 2021).In contrast, ΔCO 2 represents the balance between lake input of C, in-lake factors, and export via F CO 2 .Lake input of C comes from the surrounding catchment (Duarte and Prairie 2005;Maberly et al. 2013;Marcé et al. 2015) and is modified by catchment size and biogeohydrological properties and the magnitude and intensity of rainfall.In-lake factors constitute, for example, uptake, processing, and respiration of C by aquatic organisms (Lauster et al. 2006;Cole et al. 2007;Reis and Barbosa 2014) and are modified by light intensity, nutrient concentrations, and the standing stock of organisms.
The many factors involved makes it challenging to assess the relative influence of k and ΔCO 2 on F CO2 , especially because the effect of such factors is manifested at different temporal scales ranging from hours or days (e.g., diurnal changes in winds, temperature, primary production, and respiration) to seasons or years (e.g., seasonal and annual variability in temperature, nutrient availability, and C input).Besides, k and ΔCO 2 may interact in contrasting ways, making it difficult to interpret their individual effect on F CO 2 .For example, during periods of stratification CO 2 may accumulate in the hypolimnion and the rate at which CO 2 can be transported to the surface depends on how effective winds and cooling are at deepening the upper mixed layer.These mixing-effects can range from being transient and associated with day-night patterns of wind and temperature (Jammet et al. 2017;Czikowsky et al. 2018) to being seasonal and linked to lake mixing in spring or autumn (L opez Bellido et al. 2009;Ducharme-Riel et al. 2015;Jansen et al. 2019).
In addition, temporal patterns in ΔCO 2 and k can vary from one location to another in a lake (Natchimuthu et al. 2017).For example, at shallower depths the time for CO 2 from the sediment to reach the surface can be shorter, and the sediment is more sensitive to resuspension, which can increase near-shore ΔCO 2 (Tonetta et al. 2016).In contrast, there can be lower ΔCO 2 at nearshore locations in lakes with uptake of CO 2 by primary production in the littoral zone (Schilder et al. 2013;Peixoto et al. 2016).Winds can also produce local patterns in turbulence and k (Read et al. 2012;Schilder et al. 2013;Vachon et al. 2013) and may tilt the thermocline, causing upwelling of CO 2 resulting in local increases in F CO 2 in both small lakes (Natchimuthu et al. 2017;Mac-Intyre et al. 2021;Rudberg et al. 2021) and larger ones (Heiskanen et al. 2014;Czikowsky et al. 2018).Alongside within-lake variability in ΔCO 2 and k, variability in lake depth, size and morphology, catchment hydrology (discharge and precipitation), and temperature and wind patterns can lead to between-lake variability of k and ΔCO 2 .
Interaction effects between k and ΔCO 2 have been identified in river systems (Rocher-Ros et al. 2019), where turbulence is generally high.There, ΔCO 2 was regulated by emission rate and thus covaried with k in the sense that ΔCO 2 was low when k was high and vice versa.Although multiple studies have investigated the variability of either ΔCO 2 or k in lakes (e.g., Natchimuthu et al. 2017;Jansen et al. 2020;MacIntyre et al. 2021), studies comparing the relative influence of ΔCO 2 or k on F CO 2 variability at different temporal scales and across multiple lakes are rare.Yet, such knowledge is needed for fundamental understanding of gas fluxes in lakes and for improving process-based modeling of F CO 2 (Tan et al. 2017).It is also key for assessing when and where ΔCO 2 and k need to be measured to reliably assess F CO2 .
Here we used floating chambers and surface water samples to measure F CO 2 and ΔCO 2 , respectively, at deep and shallow locations for short-term (1-9 d) and monthly intervals in 18 lakes, in hemi-boreal, boreal, and sub-arctic regions in Sweden.We combined these measurements to derive spatially resolved values of k and identified: (i) changes in variability of F CO2 , ΔCO 2 , and k over time, (ii) the relative contributions of ΔCO 2 and k to F CO 2 variability for different time scales, and (iii) how differences in the relative contributions of ΔCO 2 and k to F CO2 variability at both short-term (1-9 d) and seasonal time scales can be related to the characteristics of lakes and other environmental variables that affect processes within lakes.Based on our findings, we also discuss how to effectively assess the variability in F CO 2 across spatiotemporal scales in lakes.

Study lakes and sampling design
Eighteen lakes were studied between 2018 and 2020 (six lakes per year) (Supporting Information Table S1; Fig. 1).The lakes are distributed within six regions spanning sub-arctic to hemi-boreal environments in Sweden.Lakes were chosen to cover landscape characteristics and carbon and nutrient conditions representative of ranges observed in Swedish lakes (Supporting Information Fig. S1).Ranges of dissolved organic carbon (DOC), absorbance at 420 nm and total phosphorous (TP) were: 2-22 mg L À1 , 0.006-0.11cm À1 , and 4-280 μg L À1 .
All lakes were visited monthly for week-long sampling campaigns during the full ice-free period for 1 yr.Sampling was conducted on 2-5 occasions during each week-long visit, with 24-48 h between each sampling.Samples of pCO 2aq and pCO 2air were collected to derive ΔCO 2 and combined with F CO2 measurements to estimate k along three horizontal transects from shore toward the lake's center.Each transect consisted of 2-4 sampling locations that represented both shallow ($ 0-2 m) and deep (> 2 m) parts of the respective lakes (Fig. 1).Typically, shallow and deep locations constituted near-shore and offshore locations, respectively, but for logistical reasons the most central part of the lakes were not always covered.In three of the shallowest lakes, Dammsjön (DAM), Lammen (LAM), and Grinnsjön (GRI), and in BD6, where sampling was limited to one of two bays, deep sampling was partly made at depths less than 2 m and reflects longer distance to shore rather than a great depth difference.
Floating chambers were deployed at each sampling location and were used to measure both F CO 2 and pCO 2aq .The chamber design was similar to that used in Natchimuthu et al. (2017) and each chamber consisted of an inverted plastic bucket (volume: 5.9-8.6 liters, diameter: 0.3 m), covered with aluminum tape to minimize internal heating by sunlight.The chamber had floats attached to its sides allowing the chamber walls to immerse approximately 3 cm into the water.The chambers were attached to another small float with a $ 1-m line, which in turn was connected to weights; this mooring allowed the chambers to move up and down with waves and water level changes.This type of chamber design yield similar flux estimates as other methods that do not interfere with the water surface (Cole et al. 2010;Gålfalk et al. 2013).
A 0.2-to 0.3-m polyurethane tube (inner and outer diameter of 3 and 5 mm, respectively) connected to a three-way Luer-lock valve was attached through the top of the chamber for manual collection of CO 2 from inside the chamber.
We retrieved air temperature and atmospheric pressure at 2 m height, wind speed at 10 m height, and precipitation from the meteorological reanalysis model MESAN developed by the Swedish Meteorological and Hydrological Institute, which interpolates measurements from nearby weather stations combined with a meteorological model to determine hourly means of meteorological variables on a 2.5 Â 2.5 km grid (Häggmark et al. 2000).Although the MESAN model provides estimates for wind speed at 10 m height, it gives consistent (R 2 = 0.78-0.87)values of wind speed and wind gust speed compared to values measured at lake level (Rudberg et al. 2021).

Field measurements
Measurements of pCO 2aq and pCO 2air pCO 2aq was measured with two separate methods; a primary method combining headspace-and dissolved inorganic carbon (DIC) samples (pCO 2HS ) and dissociation constants from Wallin et al. (2010), used in 95% of cases and a secondary method involving measurements from equilibrated floating chambers (pCO 2FC ), used in the remaining cases.Sampling of pCO 2aq and pCO 2air is only briefly described below and details are found in the Supporting Information (Supporting Information Text S2).Sampling and measuring DIC are described in Rudberg et al. (2021).
The pCO 2HS -method involved collecting water samples at $ 0.1 m depth and then extracting them into the headspace of a syringe.We used a 140-mL syringe with 35 mL of air and equilibrated the sample for 2-3 min.Values of pCO 2aq were converted from ppm (obtained from the gas analysis) to μatm, by accounting for changes in the carbonate system equilibrium in the vial headspace during sampling (Koschorreck et al. 2021), as described in Rudberg et al. (2021).In cases where headspace equilibration samples were not collected at the specific sampling location but elsewhere in the sampling transect (21% of pCO 2HS values), we instead calculated pCO 2aq using headspace samples from the closest sampling location.
In cases where no headspace measurements were collected on the sampling transect, we estimated pCO 2aq using the pCO 2FC -method.The pCO 2FC -method involved collecting gas samples from the headspace of floating chambers after deployment times of $ 24 or 48 h, which is sufficient time for the chambers to fully equilibrate with CO 2aq (Bastviken et al. 2015).The chamber headspace was sampled with three 60-mL syringes and the gas was used to flush a 20-mL vial followed by injecting sufficient gas in the vial to create an overpressure.The sample integrity was confirmed when the overpressure in the vial was released prior to analysis in the lab (see below).Values of pCO 2aq were converted from ppm to μatm by multiplying with the pressure during sampling.There is a k-dependent delay in capturing the real pCO 2aq with floating chambers, which has been estimated to 1 and 3 h for wind speeds > 6 and < 1 m s À1 , respectively (Bastviken et al. 2015).Such delay times imply a smoothing of short-term temporal variability of pCO 2aq .The bias in pCO 2FC caused by this delay was investigated by comparing pCO 2HS and pCO 2FC -samples that were collected within a few minutes of each other.A standardized major axis (SMA) linear regression between the two methods showed no consistent bias (slope: 1.04, R 2 = 0.8; Supporting Information Fig. S2), after removing outliers (7 of 790 values).
Sampling of pCO 2air was made the same way as for the pCO 2FC -method by collecting air directly above the water surface.Assuming that variability of pCO 2air was small within a transect, it was only sampled from the location nearest to the shore at each transect.
Values of pCO 2air and pCO 2aq were converted from ppm to μatm by multiplying with the pressure during sampling.ΔCO 2 (mmol m À3 ) was derived by first subtracting transect location pCO 2air from sampling location pCO 2aq and then multiplying with Henry's constant (M atm À1 ), calculated from surface water temperatures at the time of sampling according to Weiss (1974).

Measurements of F CO 2
Each floating chamber was equipped with a CO 2 sensor (K33 ELG, SenseAir AB) to measure CO 2 concentration in the chamber headspace at 5-to 6-min intervals.The sensors also logged temperature and relative humidity.These CO 2 sensors have proved reliable under non-condensing conditions (Bastviken et al. 2015).The CO 2 sensor was placed inside the chamber within a plastic box attached to the chamber top, to protect from splashes.The box was drilled with multiple 0.5to 1-cm-diameter holes to facilitate gas transport.To protect the sensors against humidity, critical areas on the circuit boards were coated with polyurethane resin.Before each field measurement the sensors were cross calibrated for approximately 3 h by comparison to reference CO 2 measurements made with an Ultraportable Greenhouse Gas Analyzer (Los Gatos Research).Sensors whose measurements deviated from the reference were not used in field measurements.Following retrieval of data from CO 2 sensors, the data were filtered to account for humidity peaks, and the CO 2 concentration was corrected for temperature and humidity effects, as described in Bastviken et al. (2015).In nine cases, we identified nonlinear patterns in our sensor measurements associated with shortterm instabilities and that sensor CO 2 increase at the time near deployment was not reflecting the long-term trend.Those cases were excluded from further analysis.
F CO 2 was calculated through an iterative procedure, accounting for measurements within a time window of 6-45 min from the chamber deployment.Measurements made within the first 6 min after deployment were discarded to minimize bias from incomplete mixing inside the chamber or possible turbulence caused by deploying the chamber.A slope representing CO 2 change over time was estimated for the first three measurements within the time window, and additional slopes were calculated after including one measurement point at a time for the whole time window.Among these slopes, the one with highest R 2 was used for calculation of F CO2 according to Eq. (3).Cases when the highest R 2 of the slopes was < 0.9 were only considered if the root mean square error was below 10 ppm to not discriminate against small fluxes, as R 2 generally decreases when the surface water CO 2 concentration is close to equilibrium with the atmosphere.For the calculation of F CO2 in Eq. ( 2), Δf =Δt represents the change in mole fraction over time (ppm s À1 ), P is atmospheric pressure (atm) at the time of sampling, V and A correspond to the volume (m 3 ) and area (m 2 ) of the chamber headspace, T is the air temperature (K) during sampling, derived from the chamber CO 2 sensor, R is the ideal gas constant (L atm K À1 mol À1 ) and 86.4 Á 10 3 is the number of seconds per day, used to convert units from s À1 to d À1 .
Values of pCO 2aq , pCO 2air and associated chamber CO 2 sensor values allowed further assessment on data quality.We checked whether the concentration change measured by the CO 2 sensor (i.e., increase or decrease) over the first hour after deployment fully equilibrated with water (pCO 2aq minus pCO 2air ), which was seen as unrealistic as it normally takes > 1 h for that to occur (Bastviken et al. 2015).To limit the effect of outliers, we discarded measurements where the above criterion was violated for at least three consecutive sensor measurements.Six hundred twenty of 1810 (34%) measurements did not pass this quality assessment and were excluded from analysis.Our motivation for this filtering approach and potential impacts on results are discussed in detail in Supporting Information Text S3.Briefly, 32% of the discarded measurements were made at such near-saturation conditions (jpCO 2aq À pCO 2air j < 100 μatm), when the uncertainty of k estimates is high (Supporting Information Fig. S3).An additional 14% of the discarded CO 2 sensor measurements used for estimating fluxes were made during reducing air temperatures when relative humidity neared 100%, conditions where condensation inside the chamber headspace can potentially bias measurements (but such conditions did not always lead to data discards).One possible explanation for the remaining ($ 50%) discards is that pCO 2aq derived from syringe sampling at $ 0.1 m depth (pCO 2HS method) was lower than the pCO 2aq in the surface water boundary layer, potentially due to microbial processes or photochemical reactions in the surface microlayer, enhancing CO 2 production (Conrad and Seiler 1988;Zhou and Mopper 1997;Cunliffe et al. 2011), which would contribute to underestimating pCO 2aq and result in overestimating k (Pajala et al. 2023), and in turn more rejected data.

Estimation of k from F CO 2 and ΔCO 2
Values of k were derived from sampling location and timespecific F CO 2 and ΔCO 2 , by solving Eq. (1) for k.The direction of F CO 2 is out of the lake (i.e., positive) when ΔCO 2 is positive, and likewise F CO2 is into the lake (i.e., negative) when ΔCO 2 is negative.Hence, k can only be positive.Nevertheless, 7% of our k-estimates were negative, which can happen near equilibrium (i.e., when ΔCO 2 is low) due to uncertainties in the measurements.These values were excluded from the analysis.

Additional lake chemistry and environmental variables
Values of pH in surface water were either derived indirectly by solving the carbonic acid equilibria equation (Stumm and Morgan 1996), using monthly measurements of DIC and pCO 2aq at $ 0.1 m depth at each chamber position as described in Rudberg et al. (2021), or measured directly at $ 0.1 m depth at the lake center with a pH electrode (HACH IntelliCAL PHC201).Measurements of surface water pCO 2aq is explained above.Surface water DIC was sampled with a syringe according to Rudberg et al. (2021).In short, 5 mL of water was injected to a 22-mL N 2 -filled vial holding 100 μL concentrated phosphoric acid.The vial was kept at 20 C before the headspace gas of the vial was analyzed with gas chromatography (see below).DIC concentration in the water phase of the vial was derived by applying Henry's law to the concentration in the vial headspace.The in situ DIC concentration was calculated as the combined molar DIC amount in the water phase and headspace of the vial divided by the volume of the vial water phase.
Near-surface water (0.5-1 m depth) measurements of chlorophyll a (Chl a), total nitrogen (TN), TP, total organic carbon (TOC), DOC, and absorbance were made by sampling water with a Ruttner sampler at the deepest location of the lake.Part of this water was used to fill 125-mL acid-washed high-density polyethylene bottles that were transferred to a freezer upon arrival to the lab and kept at À18 C before being analyzed for TP.The remaining water was stored in a 4-L cubitainer until arriving at the lab the same day.At the lab, part of this water was transferred to 50-mL polypropylene tubes and stored at 4 C, to be analyzed for duplicate measurements of TN and TOC.Another part of the water was filtered using 47-mm Whatman GF/F glass fiber filters (nominal pore size 0.7 μm) with 100 mL of the filtrate transferred to two 50-mL polypropylene tubes and stored at 4 C to be analyzed for duplicate measurements of DOC and the remaining filtrate used for measurements of absorbance.The remaining water in the cubitainer was used for filtering with GF/C filters, upon which the filters were folded and wrapped in aluminum foil and stored at À18 C to be analyzed for Chl a. TP was analyzed with a segmented flow analyzer (AutoAnalyzer II, Seal Analytical) according to Murphy and Riley (1962).TN, TOC, and DOC were analyzed with a Total Organic Carbon Analyzer (TOC-V CSH, Shimadzu with a TN analysis unit), and certified standards containing 1, 2, and 5 mg L À1 TN and 10 mg L À1 TOC were used for calibration.Samples that we assumed had TOC concentrations higher than 20 mg L À1 were diluted with deionized water prior to analysis of TOC and DOC.Light absorbance at a 420 nm wavelength was measured with a spectrophotometer (Aqualog, Horiba).Chl a was extracted with pure methanol according to Swedish standard procedures (SS 28170) with a method based on Richards and Thompson (1952) and was analyzed with spectrophotometry.
Depth gradient measurements of temperature and light intensity in the photosynthetically active wavelength range (400-700 nm) were made at monthly intervals during the icefree period at the lake center.Temperature measurements were made vertically at 0.5-to 2-m intervals with a multi-parameter probe (AP-5000 Aquaprobe or LD0101 IntelliCAL HACH probe).Light intensity was measured with an underwater quantum sensor (LI-192SA LI-COR) every 10 cm down to 1 m depth and thereafter in 50-cm intervals.The slope of the regression between the logarithm of light intensity vs. depth was used to calculate light attenuation, according to Beer-Lamberts Law.
Dissolved oxygen (DO) saturation was measured at 0.1-0.25 m depth throughout the ice-free period using one centrally located PME miniDOT logger.Measurements were made at 10-min intervals or were recalculated to 10-min averages if higher measurement frequencies were used, in accordance with the sensor response time.Water temperature was obtained from the miniDOT loggers at 10-min intervals and at second intervals from thermistors (RBRsolo 2 , RBR) located within each sampling transect at 0.1 m depth.

Gas analysis
Analysis of gas samples and of headspace gas in the DIC samples were made with an Agilent 7890A gas chromatograph (Agilent Technologies, equipped with a 1.8 m Â 3.175 mm Porapak Q 80/100 column from Supelco, and a thermal conductivity detector) through automatic injection, with a 7697A headspace sampler (Agilent Technologies) attached.After confirming a linear and stable response by analyzing serially diluted samples from a certified high concentration standard (50,000 AE 1000 ppm) upon GC system evaluation, an independent certified standard of 1985 AE 40 ppm was used for calibration for each batch of samples being analyzed.Prior to analysis, vial overpressure was released by inserting a 0.5-mm hypodermic needle holding a water droplet making outflow of excess gas visible and forming a barrier to gas inflow.This setup made it possible to verify sample integrity during storage.Stopper CO 2 permeability was characterized by independent experiments and sample analyses were corrected to sample storage time.

Data analysis
Temporal variability in F CO2 , ΔCO 2 , and k Interquartile range (IQR) was used for assessing temporal variability in F CO 2 , ΔCO 2 , and k.By calculating IQR on datasets that represented different temporal scales in the respective lakes, we derived measures on lake-wise temporal variability.The temporal scales ranged from one monthly sampling campaign (1-9 d) and consecutive monthly intervals (i.e., multiple monthly sampling campaigns) to the full icefree period.For each lake and temporal scale, IQR was calculated using all possible monthly combinations of data belonging to that time scale, to not favor certain monthly combinations over others.The combined mean IQR of these individual estimates was calculated and correspond to the temporal variability for that lake and scale.

Variation in F CO 2 explained by ΔCO 2 or k
We used dominance analysis to estimate the relative influence of ΔCO 2 and k on F CO 2 variability (Budescu 1993;Azen and Budescu 2003).In dominance analysis, relative influence of predictor variables is derived by first calculating the additional contributions that each predictor variable provides to the model R 2 across all possible linear subset models (in our case two subset models, one with only ΔCO 2 included and one with only k included as the predictor variable).The additional contribution of each predictor is calculated as the increase in R 2 that results from adding that predictor to the regression model, and results in an output between 0 and 1, where higher output indicates greater relative influence of that predictor variable to the response variable.
We based our model used for dominance analysis on the logarithm of Eq. ( 1), satisfying the requirement of linear relationships between predictor variables, log(ΔCO 2 ) and log(k), and the dependent variable, log(F CO 2 ).We used absolute values of F CO 2 , ΔCO 2 , and k in dominance analysis because logarithms cannot be calculated on negative values.It was also preferrable to use absolute values to not bias results in cases of negative F CO 2 (otherwise, k will be positive when F CO 2 is negative, and the relationship between F CO2 and k will be non-linear if positive and negative F CO 2 are used in the same regression; Supporting Information Fig. S4).The model is displayed in Eq. (3).
Due to the small number of measurements with negative F CO 2 , we did not conduct separate dominance analyses for sets of values associated with negative and positive F CO 2 , respectively.Instead, all values were combined in the same analysis.By doing so, we presume that k and ΔCO 2 influence F CO 2 in similar ways when the flux is directed into (4% of cases) and out of (96% of cases) the lake, respectively.
Influence indexes were derived at various temporal scales by applying dominance analysis on separate subsets of data that represented measurements of F CO 2 , ΔCO 2 , and k from different temporal scales.The method was designed as an iterative Monte-Carlo procedure.The temporal scales were the same as those used for estimating temporal variability in our data (explained above).Similar to that method, individual estimates of influence indexes were calculated for each time scale, consisting of different monthly combinations of data to not favor certain monthly combinations over others (i.e., for a temporal scale of 2 months, one set included data from March and April, another April and May, and so forth).The method followed the general procedure: (1) The first month in the dataset was set as the starting month.(2) Data belonging to the starting month was added to data for consecutive month(s), where the total of months compiled depended on the monthly scale considered.(3) We identified the smallest number (n min ) of observations for any of the months included in the data.(4) We extracted a random subset of n min unique samples from each month of data, to ensure similar weights for months with unequal sample sizes.(5) We calculated influence indexes of ΔCO 2 and k to F CO 2 variability for the subset of samples with dominance analysis; 6. Steps 4 and 5 were repeated 100 times.(7) We calculated the respective mean relative influence of ΔCO 2 and k to F CO 2 variability from the iterations.( 8) Steps 1-7 were repeated but using different starting months as reference.(9) After separate mean relative influences had been estimated accordingly for all monthly combinations that represented the same time scale, the collective mean of these values was calculated and correspond to lake mean relative influence (here denoted I) for that specific time scale.
In the text, we refer to lake mean relative influence of ΔCO 2 and k on F CO 2 variability at a given temporal scale as I CO 2 and I k , respectively.The sum of I CO2 and I k is always 1, which means that I CO 2 >0.5 indicates greatest influence of ΔCO 2 on F CO 2 .Likewise, I CO 2 < 0.5 indicates greatest influence of k on F CO2 .We refer to short-term relative influence as the shortest temporal scale included here (one monthly sampling campaign; 1-9 sampling days) and refer to long-term relative influence as covering the full ice-free period in our lakes.
We calculated I CO 2 and I k for all sampling depths combined in each lake.We also investigated spatial within-lake variability by comparing results for shallow and deep sampling locations, according to definitions of "shallow" and "deep" mentioned earlier and displayed in Fig. 1.Due to limited sampling at deep locations in BD3 and BD6 (Supporting Information Fig. S5), those lakes were excluded from analyses concerning I CO 2 and I k at deep sampling locations.

Comparison with lake characteristics
To explore how much short-term and long-term I CO2 could be linked to lake characteristics, univariable linear, logarithmic, and exponential ordinary least squares (OLS) regressions as well as multivariable linear OLS regressions of I CO2 (one value per lake) were made with a range of variables (Table 1).
In multivariable regressions, the potential predictor variables listed in Table 1 were used in combinations of two at a time.Shallow and deep location I CO 2 were statistically different (p < 0.05, t-test) at long-term scales among the lakes, hence regressions were conducted on shallow and deep locations, respectively.At short-term scales, shallow and deep location I CO 2 were not statistically different ( p = 0.62, Kruskal-Wallis test) and were thus combined in the regression analysis.
All lake water and weather variables were compiled by first averaging values by sampling occasion and then calculating a total average for the respective lakes.Values of wind speed, averaged by both sampling occasion and the full year, were analyzed.Lake water and meteorological weather variables thus represent mean values for the whole lake in all regressions (i.e., in regressions with whole-lake, shallow, and deep I CO 2 ).Variable values for the respective lakes are listed in Supporting Information Table S2.
Model stability was evaluated by repeating the regression analysis but leaving one lake out each time.In cases where model significance was determined by the inclusion of one lake alone, the model was deemed unstable and was excluded from analysis.Furthermore, predictor variables included in multiple regressions were examined for collinearity by calculating variance inflation factors (VIF), which is a measure of the maximum correlation between the respective predictor variables.Significant regressions were only accepted if VIF was < 5 for the predictor variables (Sheather 2009).Finally, models were ranked based on adjusted R 2 (calculated according to Eq. 4).

Statistical analysis
Probability values below 0.05 have been considered statistically significant to reject null-hypotheses in statistical tests and regressions.Reported R 2 values are adjusted to the number of observations (n) and predictors (p) in the model In boxplots, boxes indicate the distribution of data, where the box itself represents the range of the middle 50% of data (IQR), the horizontal lines in each box show the median value, the whiskers above and below boxes show the first and last quartile of data within AE 1.5 Â IQR, and points represent outliers (> j1.5 Â IQRj).

Results
Magnitude and variability of F CO 2 , ΔCO 2 , and k Mean ice-free F CO2 , ΔCO 2 , and k ranged between 2 and 73 mmol m À2 d À1 , 4 and 78 mmol m À3 , and 0.8 and 3.2 m d À1 , respectively, for the different lakes (Table 2; Fig. 2).F CO 2 was predominantly positive in our lakes (i.e., flux out of the lake) but was negative (i.e., flux into the lake) for 3%, 6%, 16%, and 38% of measurements in oligotrophic BD6 and our three eutrophic lakes, GUN, VEN, and SOD.Mean F CO2 and ΔCO 2 corresponded well (r = 0.81) and were linked to light absorbance at 420 nm (r = 0.74 and 0.79) which in turn was linked to DOC (r = 0.9).No clear trend was identified between F CO2 and k (Fig. 2).Variability in ΔCO 2 increased from temporal scales of 1-9 d up to the full ice-free period (10 months or less).In contrast, variability in k was reduced or leveled out at scales > 5 months in most lakes (Fig. 3).Variability in F CO 2 differed between lakes.In lakes VEN, GUN, LAM, and GRI, F CO 2 variability increased at lower time scales but was reduced at scales > 6-8 months.In most other lakes, F CO 2 variability either increased or was leveled out over time.

Influence of ΔCO 2 and k on F CO 2 variability
The lake mean relative influence of ΔCO 2 on F CO2 variability was generally below 0.5 at temporal scales between 1-9 d and 6 months, meaning that k had most influence on F CO 2 variability at those temporal scales.However, I CO2 in the lakes studied here increased with greater temporal scale and at scales longer than 7 months average lake I CO 2 was instead slightly above 0.5 (i.e., ΔCO 2 had most influence on F CO2 ).Despite this general trend there was large between-lake Table 2. General information including duration of sampling campaigns, total number of sampling campaigns (n campaigns ) and observations (n obs ) and mean, standard deviation (SD), and minimum-maximum range of F CO 2 , ΔCO 2 , and k for the respective lakes.Statistics are calculated as collective sampling month means (i.e., by first averaging for each monthly sampling campaign and then calculating a combined average of those).For calculation of SD, only sampling months with > 3 measurements have been considered.Lake Sampling period n campaigns n obs   variability, with differences in I ranging between 0.3 and 0.5 depending on temporal scale.The change in I with temporal scale (increase for I CO 2 and decrease for I k ) was generally greatest at deeper sampling locations (Fig. 4) and when only these deep sampling locations were considered, I CO 2 had a clear positive trend with increased timescale in all lakes except in BD3, BD6, SOD, and GRA (Supporting Information Fig. S5).Accordingly, I CO2 at deep sampling locations was significantly greater than at shallow locations when data for the full open water period in the lakes were used (p < 0.05).

Short-term regulation of ΔCO 2 and k on F CO2 variability
For short-time periods (1-9 d) catchment slope, catchment fraction of open field, and lake area were positively related to relative influence of ΔCO 2 on F CO2 variability and combinations of those variables produced the best models (R 2 = 0.54-0.73;Fig. 5a,b; all equations and statistics are shown in Supporting Information Table S3).Of these three variables, catchment slope and fraction of open field were correlated to wind speed (r = 0.72 and 0.84).
Other significant models (R 2 = 0.37-0.48)combined either catchment slope, catchment fraction of open field, or lake area, with lake physiochemical variables (TOC, absorbance, ΔCO 2 , Chl a, DIC, DOC, TP, and light attenuation).Similar degrees of explanation were identified when catchment fraction of open field was combined with either latitude or air temperature (R 2 = 0.4-0.42),and when lake area was combined with wind speed (R 2 = 0.42 and 0.41 using mean values from sampling occasions and the full year, respectively; Fig. 5c).
Ice-free period regulation of ΔCO 2 and k on F CO 2 variability No significant relations between long-term I CO2 and driver variables were found for shallow sampling locations.When deep sampling locations were considered, yearly precipitation, lake volume, and ice-free period duration had positive relationships with long-term I CO 2 on their own (R 2 = 0.39-0.5;all equations and statistics are shown in Supporting Information Table S4).The degree of explanation increased when two of these variables were combined together (R 2 = 0.61-0.72;Fig. 6a,c), when precipitation was combined with catchment area or lake-to-catchment area ratio (R 2 = 0.7-0.72)or when ice-free period duration was combined with mean depth (R 2 = 0.64, Fig. 6b).Lake volume was also significant in regressions with DIC, TP or Chl a, but those multiple regressions mostly explained deviation in lake SOD, which had much higher values of DIC, TP, and Chl a than other lakes.

Discussion
Our study of variability in F CO 2 shows that variable influence changes depending on time scale.Over shorter time scales, k generally had greater influence than ΔCO 2 on F CO2 variability.On longer time scales, the relative influence of ΔCO 2 for F CO 2 variability increased and became the dominant factor for F CO2 variability over deep locations, which was in line with the seasonal pattern of ΔCO 2 observed in most of our lakes (Supporting Information Fig. S7) and the greater variability observed over longer time scales for ΔCO 2 compared to k (Fig. 3).Nevertheless, k was still the dominant factor for F CO 2 variability at shallow near-shore locations even on longer time scales.These time-scale differences can be generalized to ΔCO 2 contributing more than k to F CO 2 variability on longer time scales over deep locations, and k contributing more than ΔCO 2 to F CO 2 variability at shallow near-shore locations regardless of temporal scale (Fig. 7).Natchimuthu et al. (2017) found that ΔCO 2 contributed more than k to F CO 2 variability over weekly scales, whereas k contributed more at hourly scales.However, such findings were reported for one lake only and the influence of ΔCO 2 on F CO 2 variability exceeded influence of k on F CO 2 variability in only some of the 18 lakes we studied.Thus, the relative importance of k and ΔCO 2 for F CO 2 differ over time and space in lakes and among lakes.

Short-term influence on F CO 2 variability across lakes
For short-time periods (1-9 d) catchment slope, catchment fraction of open field, and lake area were most relevant for relative influence of ΔCO 2 on F CO 2 variability.Of these three variables, catchment slope and fraction of open field were correlated to wind speed.The third variable, lake area, is in turn important for wind fetch (Read et al. 2012;Vachon et al. 2013;Klaus and Vachon 2020).The above suggests high influence of ΔCO 2 on F CO 2 variability on a short-term basis in lakes with characteristics favoring high winds (Fig. 5c).A lower relative influence of k on F CO 2 variability in such lakes, which usually have high values of k (Schilder et al. 2013), may seem counterintuitive.However, the magnitude of k does not reflect variability in k, and in lakes where k is commonly stable from constant wind exposure, k also varies less and is thereby less important for F CO 2 variability.Similarly, the short-term relative influence of k on F CO 2 variability was high in small lakes sheltered by forests (LJR, NOR, NBJ, and NAS; 1-3 ha and 69-89% forest coverage), where wind speeds were below the lake average and ΔCO 2 was generally high (Supporting Information Fig. S7).There, variability in k between calm (most of the time) and higher wind speed (occasionally) can be greater than in lakes with consistently high wind speeds.In a shortterm study of a small boreal lake located directly beside NBJ, MacIntyre et al. (2021) found that k was primarily wind driven yet highly variable even at low wind speeds, which supports the possibility of variable k in small wind sheltered lakes.
Short-term relative influence of ΔCO 2 on F CO 2 variability was higher in lakes with low ΔCO 2 .In our data, the Fig. 6.Multivariate regressions for differences in long-term (ice-free period) lake mean relative influence of ΔCO 2 on F CO2 variability (I CO2 ) at deep sampling locations, including the variables lake volume, duration of ice-free period (t ice-free ), yearly precipitation and mean lake depth.In each regression, one of the variables, (a) precipitation, (b) mean depth, and (c) t ice-free , is displayed in colors ranging from low (white) to high (blue) values.BD3 and BD6 are excluded in the regressions due to limited sampling at deep locations.Dashed horizontal lines in the panels represent the divide between higher influence of ΔCO2 on F CO2 variability (above line) and higher influence of k on F CO2 variability (below line).
unproductive clearwater lakes (BD3, BD4, BD6, SGA) and the most eutrophic lake (SOD) have low levels of ΔCO 2 .In these low-ΔCO 2 lakes, even small alterations in ΔCO 2 can be of high relative importance for F CO 2 .In such situations, the interactions between ΔCO 2 and k can be extra important for F CO2 variability.For example, if minor relative changes in k are still large enough to induce sufficient transport of CO 2 from, for example, epilimnetic sediments to surface waters to generate greater relative changes in ΔCO 2 , this would imply a greater influence of ΔCO 2 than k on F CO 2 variability.
Lower levels of ΔCO 2 in more eutrophic lakes may instead be caused by higher rates of primary productivity than respiration during daytime (when most of our samples were collected), where primary production produces O 2 and consumes CO 2 and respiration produces CO 2 and consumes O 2 .
The effect of primary productivity (and thus its consumption of CO 2 ) depends on factors that vary on daily and diel timescales, such as temperature and light availability (Gomez-Gener et al. 2021).Thus, higher primary productivity can positively affect the relative influence of ΔCO 2 on short-term variability of F CO2 .In lakes where respiration and primary productivity is important for surface water CO 2 , alterations in DO saturation can reflect surface water CO 2 concentrations (Cole et al. 2000;Hanson et al. 2003;Hanson et al. 2006).We found strong correlations between DO saturation and ΔCO 2 (r = À0.88) and they also had similar importance in regressions (Supporting Information Table S3).Therefore, we believe that primary productivity can control short-term relative influence of ΔCO 2 on F CO 2 variability in eutrophic lakes.

Long-term influence on F CO 2 variability across lakes
Duration of the ice-free period, lake volume and yearly precipitation dominated among variables being most relevant for long-term relative influence of ΔCO 2 on F CO 2 variability at deep lake locations and explained up to 74% of variability when used in multiple regressions (Supporting Information Table S4).Deep large-volume lakes are likely to have more gradual and prolonged seasonal mixing periods.The effect of such mixing periods on F CO 2 variability may be accentuated in lower-latitude lakes with longer ice-free periods, where a warmer climate may enhance hypolimnetic CO 2 production and accumulation.Prolonged seasonal mixing periods and enhanced accumulation can in turn generate greater seasonal variability in F CO2 and ΔCO 2 when deep-water mixes to the surface.The precipitation effect on seasonal F CO 2 variability can be related to transport of CO 2 from the catchment to the lake with surface water or via soil and ground water.If greater mean yearly precipitation implies greater variability in precipitation and CO 2 transport to the lake, this will also trigger greater variability in F CO2 and ΔCO 2 .
Factors controlling the influence of ΔCO 2 vs. k on F CO 2 within lakes Short-term relative influence of ΔCO 2 on F CO 2 variability at shallow and deep locations was similar among our lakes, with differences in I CO2 between these two depth categories averaging at 10% for our lakes.For longer time periods, differences were greater (19% on average) and whereas we found variables that could explain a substantial portion (R 2 ≤ 0.71) of betweenlake variability in I CO 2 at deep sampling locations (Fig. 6), we found none that could explain such patterns at shallow locations.This could be due to complex interactions between ΔCO 2 and k.For example, at shallow locations a greater k will enhance not only gas exchange rates to the atmosphere, but also the gas exchange with shallow sediments, sustaining a steeper gas concentration gradient (del Giorgio and Williams. 2005;Bergström et al. 2010;Bartosiewicz et al. 2015).Hence, k variability can have a dual effect on F CO2 and thereby emerge as the most important factor for F CO 2 variability from shallow waters.In effect, ΔCO 2 variability at shallow locations may, similarly to k, be more linked to short-term patterns in wind speed and wind direction and less linked to seasonal patterns.
In addition, the use of lake mean values in regressions also means that some variables, such as wind speed, may not adequately represent near-shore conditions or may represent it differently in different lakes.This may have resulted in a lack of models for relative influence on F CO2 at shallow locations.Furthermore, other drivers for long-term relative influence on F CO 2 at shallow locations than those considered here may be important, such as groundwater or riverine input (Striegl and Michmerhuizen 1998;Humborg et al. 2010;Maberly et al. 2013), catchment geology and productivity (Maberly et al. 2013), sediment characteristics (Bergström et al. 2010;Chmiel et al. 2016), internal waves and seiches (Heiskanen et al. 2014;MacIntyre et al. 2021), and macrophyte density (Peixoto et al. 2016).These factors need to be addressed in other studies.Overall, we here show that the consideration of temporal and spatial scales to disentangle processes important for spatial and temporal variability is important when comparing patterns in F CO 2 within and among lakes.

Implications for assessing F CO2
Clearly, both k and ΔCO 2 control F CO 2 variability, but with a relative importance that differs between temporal scales.For time periods ranging from individual monthly sampling campaigns (1-9 d) to several months, k is a key variable although ΔCO 2 sets base-line conditions for fluxes.Over seasonal time periods ΔCO 2 is a main driver in many lakes, especially in lakes with large volumes and long ice-free periods, allowing for prolonged seasonal mixing periods and greater storage of DIC, or in lakes within large and precipitation-prone catchments, allowing for greater input of carbon through the catchment (Fig. 7).
The contribution from shallow and deep areas, respectively, to total lake F CO 2 corresponded to the lake area that they encompassed (shallow areas encompassed 10%-70% of lake area; Supporting Information Fig. S6), showing that both shallow and deep areas are equally important for whole lake F CO 2 .Sampling designs that combine spatially distributed measurements of monthly pCO 2aq with single-location estimates from k models that account for turbulence both within and above the lake (Tedford et al. 2014;Czikowsky et al. 2018;MacIntyre et al. 2021), may provide reliable estimates of ice-free F CO2 at deep central parts of lakes but are likely to bias F CO 2 at shallow near-shore parts of lakes.To limit bias in such sampling designs, additional measurements at shallow locations as input to long-term modeling using k would be beneficial.However, and importantly, the large within-and among-lake differences in dominance of ΔCO 2 vs. k on F CO2 combined with previous concerns regarding lack of reliability in commonly used wind-based k models (Schilder et al. 2013;Klaus and Vachon 2020) lead to the general advice to prioritize direct spatiotemporally distributed measurements of F CO 2 whenever possible or alternatively to assess k from models and with sampling designs that incorporate spatial variability of lake thermal structure and meteorology (e.g., MacIntyre et al. 2021).
Constraining the current uncertainty in F CO 2 is critical for accurate regional and global budgets.Our work here shows that variability over the short term depends on the gas transfer velocity but over the long term depends on variability in the concentration gradient of CO 2 at the air-water interface.We have shown that the latter depends on lake volume, precipitation, and duration of the ice-free period in boreal and sub-arctic lakes, implying that duration of seasonal mixing periods, hypolimnetic storage of DIC and input of carbon from the catchment are key factors.While modeling of emissions from lakes and reservoirs on global scales uses typical values of k and of ΔCO 2 over broad spatial areas (Raymond et al. 2013), our method based on lake-attributes and landscape characteristics provides an approach for further refinement.Further improvement will come from studies such as this one which identify key drivers and with development of more effective low-cost sensor measurements for CO 2 assessments, facilitating estimates of CO 2 emissions from diverse environments.

Fig. 1 .
Fig. 1.Location of lakes and sampling locations of F CO2 and pCO 2aq within the study lakes, with red and blue markers denoting shallow and deep sampling locations, respectively.Blue color toning represents depth intervals of 2 m, calculated in ArcGIS (3D Analyst, Natural Neighbor).The colors in the main map correspond to lakes within frames of similar color.Horizontal white scale bars have length of 100 m.Orthophotos were retrieved from Lantmäteriet, Sweden.

Fig. 2 .
Fig.2.Magnitude and variability of F CO2 , ΔCO 2 , and k for all study lakes.Lakes are sorted by median F CO2 , displayed with red lines.Darker shading of boxes in the boxplot represents greater light absorbance measured at 420 nm.Colors of lake names in the horizontal axis are connected to similarly colored lake regions in Fig.1.Colors above lake names denote trophic state (blue: oligotrophic, yellow: mesotrophic, green: eutrophic) based on TP, according toCarlson (1977).

Fig. 3 .
Fig. 3. Temporal variability of F CO2 , ΔCO 2 , and k, represented by the IQR of measured values aggregated by periods of 1-10 monthly sampling campaigns.The black lines show the average across lakes with at least eight consecutive months of data.

Fig. 4 .
Fig. 4. Distribution of lake mean relative influence of ΔCO 2 on F CO2 variability (I CO2 ) at different temporal scales-(a) 1-10 sampling campaigns and (b)full ice-free periods-for shallow (red) and deep (blue) sampling locations in our study lakes.Lakes BD3 and BD6 are excluded in the analysis for deep locations due to limited sampling there.Boxes show the distribution of I CO2 values from the different lakes (each box with its horizontal lines, whiskers, and points represents maximum one I CO2 value per lake) and colored areas cover values within one standard deviation from the mean.Colored numbers below each boxplot pair denote the number of lakes included in the assessments.The dashed horizontal line indicates the separation between greater influence of ΔCO 2 on F CO2 variability (above the line; influence of ΔCO 2 > 0.5) and greater influence of k on F CO2 variability (below the line; influence of ΔCO 2 < 0.5).Individual estimates of relative influence used for shallow and deep sampling locations in the lakes are displayed in Supporting Information Fig.S5.Asterisk in (b) imply significant differences for the relative influence of ΔCO 2 on F CO2 variability (and likewise significant differences for the relative influence of k on F CO2 ) at shallow and deep sampling locations.

Fig. 5 .
Fig. 5. Multivariate regressions for short-term (1-9 d) lake mean relative influence of ΔCO 2 on F CO2 variability (I CO2 ), using data from all sampling depths, together with lake area and (a) mean catchment slope, (b) catchment fraction of open field, and (c) sampling mean wind speed.In each regression, the variable lake area is displayed in colors ranging from low (white) to high (blue) values on a logarithmic scale.

Fig. 7 .
Fig. 7. Conceptual flow-chart for lake mean relative influence of ΔCO 2 and k on F CO2 variability.Red and blue lines and boxes denote shallow and deep locations, respectively.The gray box denotes shallow and deep locations combined.Text in colored boxes separate between results (normal font) and implications (bold font).The bottom box summarizes potential explanations for results.