Drivers of diffusive CH4 emissions from shallow subarctic lakes on daily to multi-year timescales

Lakes and reservoirs contribute to regional carbon budgets via significant emissions of climate forcing trace gases. Here, for improved modelling, we use 8 years of floating chamber measurements from three small, shallow subarctic lakes (2010–2017, n= 1306) to separate the contribution of physical and biogeochemical processes to the turbulencedriven, diffusion-limited flux of methane (CH4) on daily to multi-year timescales. Correlative data include surface water concentration measurements (2009–2017, n= 606), total water column storage (2010–2017, n= 237), and in situ meteorological observations. We used the last to compute nearsurface turbulence based on similarity scaling and then applied the surface renewal model to compute gas transfer velocities. Chamber fluxes averaged 6.9±0.3 mg CH4 m−2 d−1 and gas transfer velocities (k600) averaged 4.0± 0.1 cm h−1. Chamber-derived gas transfer velocities tracked the powerlaw wind speed relation of the model. Coefficients for the model and dissipation rates depended on shear production of turbulence, atmospheric stability, and exposure to wind. Fluxes increased with wind speed until daily average values exceeded 6.5 m s−1, at which point emissions were suppressed due to rapid water column degassing reducing the water–air concentration gradient. Arrhenius-type temperature functions of the CH4 flux (E a = 0.90± 0.14 eV) were robust (R2 ≥ 0.93, p<0.01) and also applied to the surface CH4 concentration (E a = 0.88± 0.09 eV). These results imply that emissions were strongly coupled to production and supply to the water column. Spectral analysis indicated that on timescales shorter than a month, emissions were driven by wind shear whereas on longer timescales variations in water temperature governed the flux. Long-term monitoring efforts are essential to identify distinct functional relations that govern flux variability on timescales of weather and climate change.


Introduction 34
Inland waters are an important source of the radiatively active trace gas methane (CH4) to the atmosphere 35 Cole et al., 2007). On regional to global scales, an estimated 21-46% of ice-free 36 season CH4 emissions from lakes, ponds and reservoirs occur via turbulence-driven diffusion-limited gas 37 exchange ( gas transfer models are increasingly used, for example in regional emission budgets (Holgerson and 40 Raymond, 2016; Weyhenmeyer et al., 2015). Fluxes computed with modelled gas transfer velocities agree 41 to a certain extent with floating chambers and the eddy covariance technique in short-term 42 intercomparison campaigns (Bartosiewicz et al., 2015;Crill et al., 1988;Erkkilä et al., 2018). However, long-43 term comparisons are needed to identify weather-and climate related controls on the flux that are 44 appropriate for seasonal assessments. Considering the increased use of process-based approaches in 45 regional emission estimates (Tan and Zhuang, 2015), understanding the mechanisms that drive the 46 components of the diffusive flux is imperative for improving emission estimates. 47 48 1

.1 Drivers of diffusive CH4 emissions 49
Diffusive fluxes at the air-water interface are estimated with a two-layer model (Liss and Slater, 1974): 50 The flux [mg m −2 d −1 ] depends on the concentration difference across a thin layer immediately below 51 the air-water interface (Δ[CH4] in mg m −3 ), of which the upper boundary is in equilibrium with the 52 atmosphere ( , ) and the base represents the bulk liquid ( ), and is limited by the gas transfer 53 velocity k [m d −1 ]. k has been conceptualized as characterizing transfer across the diffusive boundary layer. 54 Other models envision exchange as driven by parcels of water intermittently in contact with the 55 atmosphere. In these surface renewal models, k depends on the frequency of the renewal events 56 (Csanady, 2001;Lamont and Scott, 1970). The resulting calculation for k is based on the Kolmogorov 57 velocity scale, uη = (εν) 1/4 where ε is dissipation rate of turbulent kinetic energy (TKE) and ν is kinematic 58 viscosity (Tennekes and Lumley, 1972). Progress has been made in understanding how to compute ε and 59 gas transfer rates as a function of wind speed and the heating and cooling at the lake's surface (Tedford 60 et al., 2014). Comparisons between models and other flux estimation methods, such as eddy covariance, 61 illustrate the improved accuracy when computing gas transfer velocities using a turbulence-based as 62 opposed to wind based models (Czikowsky et al., 2018;Heiskanen et al., 2014;Mammarella et al., 2015).

64
A key control on emissions is the periodicity at which dissolved gases are brought to the air-water 65 interface. During stratification, the density gradient makes it difficult for wind driven mixing to bring gases 66 to the surface, and they may accumulate in the stratified regions. Conversely, thermal convection 67 associated with surface cooling can deepen the mixed layer and transfer stored gas to the surface (Crill et  timescales and through the rate of sediment methanogenesis on longer timescales; the diurnal cycle of 96 insolation may have a limited effect on production because the heat capacity of the water buffers the 97 temperature signal (Fang and Stefan, 1996). Similar phase lags and amplifications may lead to hysteretic 98 flux patterns, such as cold season emission peaks due to release of gases from the hypolimnion in dimictic 99 lakes (Encinas Fernández et al., 2014; López Bellido et al., 2009) or thermal inertia of lake sediments (Zimov 100 et al., 1997). Spectral analysis of the flux and its components can improve our understanding of the flux 101 variability by quantifying how much power is associated with key periodicities (Baldocchi et al., 2001).

103
Here we present a high-resolution, long-term dataset (2010-2017) of diffusive CH4 fluxes from three 104 subarctic lakes estimated with floating chambers (n = 1306), and fluxes obtained by modelling using in situ 105 meteorological observations and surface water concentrations (n = 535). The surface renewal model is 106 used to compute gas transfer velocities. Arrhenius relationships of Δ[CH4] and fluxes of CH4 are also 107 calculated. Using spectral analysis of our time series data, we distinguish the temporal dependency of 108 abiotic and biotic controls on the flux. The effects of lake size and wind exposure are illustrated by 109 comparing results from the 3 different lakes. 110

Floating chambers 126
We used floating chambers to directly measure the turbulence-driven diffusive CH4 flux across the air-127 water interface (Fig. 1). They consisted of plastic tubs covered with aluminium tape to reflect incoming  128  radiation and were equipped with polyurethane floats and flexible sampling tubes capped at one end with  129 3-way stopcocks (Bastviken et al., 2004). Depending on flotation depth, each chamber covered an area 130 between 610 and 660 cm 2 and contained a headspace of 4 to 5 litres. Chambers were deployed in pairs 131 with a plastic shield mounted 30 cm below one chamber of each pair to deflect methane bubbles rising 132 from the sediment. Every 1-2 weeks during the ice-free seasons of 2010 to 2017, 2-4 chamber pairs were 133 deployed in Villasjön and 4-7 chamber pairs in Inre and Mellersta Harrsjön in different depth zones (Fig.  134 1).  concentrations computed with Eq. 2. The rate increase associated with the mean 24h flux corrected for 249 headspace accumulation is shown as a dashed red line (Eq. 1 with kch from Eq. 2, or Eq. 3 with c1 = 1.21). 250 Labels denote fluxes calculated from the linear regression slopes (Eq. 3, black) and from Eq. 2 (red). 251

Computing gas transfer velocities with the surface renewal model 252
We used the surface renewal model (Lamont and Scott, 1970) where the hydrodynamic and thermodynamic forces driving gas transfer are expressed, respectively, as  the TKE dissipation rate at half-hourly time intervals: 293 where * is the water friction velocity [m s −1 ], is the von Kármán constant, is depth below the water 294 surface (0.15 m, the depth for which Eq. 5 was calibrated). We determined * from the air friction velocity 295 * assuming equal shear stresses ( ) on both sides of the air-water interface; = * 2 = * 2 , and 296 taking into account atmospheric stability (MacIntyre et al., 2014; Tedford et al., 2014). is the buoyancy 297 flux [m 2 s −3 ], which accounts for turbulence generated by convection (Imberger, 1985): 298 Here, is the thermal expansion coefficient [m 3 K −1 ] (Kell, 1975 We binned data to assess correlations between the flux and environmental covariates. Half-hourly values 320 of water temperature and wind speed were averaged over the deployment period of each chamber 321 (fluxes), and over 24 hours prior to the collection of each water sample (concentrations), reflecting the 322 mean residence time of CH4 in the water column. Fluxes, concentrations and k-values were then binned in 323 10 day, 1 °C and 0.5 m s −1 bins to obtain relationships with time, water temperature and wind speed, 324 respectively. The 10 day bins typically contained at least one sampling day for each overlapping year, and 325 enabled representative averaging across years. Lake-dependent variables (e.g. flux) were normalized by 326 lake to obtain a single timeseries (divided by the lake mean, multiplied by the overall mean). 327 328

Calculation of the empirical activation energy 329
Chamber and modelled fluxes as well as concentrations were fitted to an Arrhenius-type temperature 330 function ( where is the Boltzmann constant (8.62 × 10 −5 eV K −1 ) and T is the water temperature in K. The empirical 332 activation energy ( ′ , in electron volts (eV), 1 eV = 96 kJ mol −1 ) was computed with a linear regression of 333 natural logarithm of the fluxes and concentrations onto the inverse temperature (1/K), of which b is the 334 intercept. 335 336 2.12 Timescale analysis: power spectra and climacogram 337 We computed power spectra for near-continuous timeseries of the surface sediment, water-and air 338 temperature and the wind speed according to Welch's method (pwelch in MATLAB 2018a), which splits 339 the signal into overlapping sections and applies a cosine tapering window to each section (Hamming, 340 1989). Data gaps were filled by linear interpolation. We removed the linear trend from original timeseries 341 to reduce red noise, and block-averaged spectra (8 segments with 50% overlap) to suppress aliasing at 342 higher frequencies. We normalized the spectral densities by multiplying by the frequency and dividing by 343 the variance of the original timeseries (Baldocchi et al., 2001).

345
We evaluated our discontinuous (fluxes, concentrations) and continuous (meteorology) timeseries with a 346 climacogram, an intuitive way to visualize a continuum of variability (Dimitriadis and Koutsoyiannis, 2015). 347 It displays the change of the standard deviation (σ) with averaging timescale (tavg). Variables were 348 normalized by lake to create a single timeseries at half-hourly resolution (e.g. 48 entries for each 24-hour 349 chamber flux). To compute each standard deviation (σ(tavg)) data were binned according to averaging 350 timescale, which ranged from 30 minutes to 1 year. Because of the discontinuous nature of the datasets, 351 n bins were distributed randomly across the time series. We chose n = 100000 to ensure that the 95% 352 confidence interval of the standard deviation at the smallest bin size was less than 1% of the value of σ 353 (Sheskin, 2007). To allow for comparison between variables we normalized each σ-series by its initial, 354 smallest-bin value: σnorm = σ/σinit. Chamber fluxes averaged 6.9 mg m -2 d -1 (range 0.2-32.2, n = 1306) and closely tracked the temporal 371 evolution of the surface water concentrations (mean 11.9 mg m -3 , range 0.3-120.8, n = 606), with the 372 higher values in each lake measured in the warmest months (July and August, Fig 4a,e). Diffusive fluxes 373 increased with wind speed and water temperature (Fig 4b,c). Reduced emissions were measured in the 374 shoulder months (June and September) and were associated with lower water temperatures. We also 375 observed abrupt reductions of the flux at wind speeds lower than 2 m s -1 and higher than 6.5 m s -1 . Surface 376 water concentrations generally increased with temperature and peaked in the summer months, but unlike 377 the chamber fluxes they decreased with increasing wind speed (Fig. 4f,g). Relationships with wind speed 378 were approximately linear, while relationships with temperature fitted an Arrhenius-type exponential 379 function (Eq. 7). Activation energies were not significantly different when using either surface water or 380 sediment temperature (Ea′ = 0.90 ± 0.14 eV, R 2 = 0.93, Ea′ = 1.00 ± 0.17, R 2 = 0.93, respectively, mean ± 95% 381 CI). The fluxes, concentrations and the wind speed were non-normally distributed (Fig. 4d,h,o). Surface 382 water temperatures (0.1-0.5 m) were normally distributed around the mean of each individual month of 383 the ice-free season (Fig. 4n), but the composite distribution was bimodal.

385
Fluxes computed with the surface renewal model (Eq. 1 using kmod) closely resembled the chamber fluxes 386 (Eq. 3) in terms of temporal evolution (Fig. 4a) and correlation with environmental drivers (Fig. 4b,c). Mean 387 model fluxes were slightly higher than the chamber fluxes in Villasjön and Inre Harrsjön, and slightly lower 388 in Mellersta Harrsjön (Table 2). Model fluxes were significantly different between littoral and pelagic zones 389 in Inre and Mellersta Harrsjön (paired t-tests, p ≤ 0.02), reflecting spatial differences in the surface water 390 concentration (Table 2). Similar to the chamber fluxes, the air-water concentration difference (Δ[CH4]) 391 explained most of the temporal variability of the modelled emissions; both kmod (Eq. 4) and kch (Eq. 2) were 392 functions of U10 (Fig. 4k) and did not display a distinctive seasonal pattern (Fig. 4i). Modelled fluxes 393 decreased at higher wind speeds when surface concentrations decreased, and displayed a cut-off at daily 394 mean U10 ≥ 6.5 m s −1 , similar to the chamber fluxes, but not at U10 < 2.0 m s −1 . The temperature sensitivity 395 of the modelled fluxes (Ea′ = 0.97 ± 0.12 eV, mean ± 95% CI, R 2 = 0.94) did not differ significantly from that 396 of the chamber fluxes. 397 398

400
Model fluxes for each lake were computed with distinct scaling parameter values (Supplementary Fig. 1

Meteorology and mixing regime 451
Throughout the ice-free season the lakes were weakly stratified (Table 3). Figure 5 shows a timeseries of 452 the mixed layer depth and water temperature in the deeper lakes, along with wind speed, air temperature 453 and precipitation for the ice-free period of 2017. The ice-free period consisted of two phases. In the first, 454 air and surface water temperatures were higher and the two deeper lakes stratified. Wind speeds 455 increased to mean values approaching 5 m s -1 for a few days at a time and then decreased for a day or 456 two. Deep mixing events followed surface cooling and heavy rainfall. Water level maxima and surface 457 temperature minima were observed 2-3 days after rainfall events, for example between 15 and 18 July 458 2017 (Fig. 5e). In the second phase, wind speeds were persistently higher (U10 > 5 m s −1 ), air and surface 459 water lakes indicative of residual CH4 that had not evaded immediately after ice-off, around 1 June 2017 (Fig.  470 5c,d). As residual CH4 was emitted, near surface concentrations declined, and then in the first half of the 471 stratified period (July 2017, Fig. 5d), particularly in Mellersta Harrsjön, increased with increased rainfall 472 and with temperature. During this period, kch and kmod were similar. Decreases in kch occurred when air 473 temperatures increased above surface water temperatures in the day leading to a stable atmosphere and 474 when near surface temperatures were warmer, and depending upon the lake, stratified to the surface. 475 Thus, lower fluxes occurred during the second part of the stratified period (August 2017, Fig. 5c) when 476 surface concentrations increased during warming periods when winds were light, the atmosphere was 477 stable during the day, and the upper water column was strongly stratified. Fluxes and concentrations were 478 lower in the autumn mixed periods, by which time the lakes had degassed, and with the colder surface 479 sediment temperatures, rates of production had decreased. 480 The modelled gas transfer velocity generally followed the temporal pattern of the wind speed (Fig. 4b). 481 Due to model calibration, the modelled gas transfer velocities (Fig. 4b, blue line) tracked those derived 482 from chamber observations (Fig. 4b, orange rhombuses). Discrepancies pointed to a mismatch between 483 24-hour integrated chamber fluxes and surface concentrations measured at a single point in time. For 484 example, measuring a low surface concentration in the de-gassed water column after a windy period 485 during which the surface flux was high led to an overestimated kch on 21 September 2017. Contrastingly, 486 kch was lower than kmod on 3 August 2017 due to elevated surface concentrations and a low chamber flux 487 associated with a warm and stratified period preceding water sampling.

489
The temperature of the surface mixed layer exceeded the air temperature by 1.6 °C on average (Fig. 5a), 490 such that the atmospheric boundary layer over the lakes was often unstable, particularly at night during 491 warm periods as well as during the many cold fronts. We computed an unstable atmosphere over the lakes 492 (z/LMO,a < 0, where z is the measurement height and LMO,a is the air-side Monin-Obukhov length; Foken 493 2006) ∼76% of the time during ice-free seasons. Atmospheric instability increases sensible and latent heat 494 fluxes (Brutsaert, 1982), enhancing the cooling rate. Thus, buoyancy fluxes were positive at night and 495 during cold fronts throughout the ice-free season (Fig 5b, Fig. 4i-k). The magnitude of buoyancy flux during 496 cooling periods tended to range from 10 −8 to 10 −7 m 2 s −3 in the stratified period and decreased as water 497 temperatures cooled in autumn (Fig. 4i,j). TKE dissipation rates at 0.15 m were high, with values often 498 between 10 −6 and 10 −5 m 2 s −3 , although values did fall as low as 10 −8 m 2 s −3 when winds were light. 499 Comparison of these two terms indicated that buoyancy flux during cooling was typically two orders of 500 magnitude less than ε and was only equal to it during the lightest winds (Fig. 4k). Consequently, its 501 contribution to the gas transfer coefficient was minor (Fig. 7). Averaged over all ice-free seasons (2009-502 2017) the buoyancy flux contributed only 8% to the TKE dissipation rate, but up to 90% during rare, very 503 calm periods (U10 ≤ 0.5 m s −1 , Fig. 4k) and up to 25% on during the warmest periods (Tsurf ≥ 18 °C, Fig. 4j). 504  505  506  507  508  509  510  511  512  513  514  515  516  517  518  519  520  521  522  523  524  525  526  527  528  529  530  531  532  533  534

CH4 storage and residence times 549
Residence times of stored CH4 varied between 12 hours and 7 days and were inversely correlated with 550 wind speed in all three lakes (OLS: R 2 ≥ 0.57, Fig. 6). The mean residence time was shortest in the shallowest 551 lake, and was not significantly different between the two deeper lakes (paired t-test, p < 0.01, Table 3). 552 We did not find a statistically significant linear correlation between the residence time and day of year or 553 the water temperature. CH4 storage was greatest in the deeper lakes and displayed patterns similar to the 554 surface concentrations, increasing in the warmest months with water temperature and decreasing with 555 wind speed.

Variability 585
Chamber fluxes and surface water concentrations differed significantly between lakes (ANOVA, p < 0.001, 586 n = 287, n = 365) ( Table 2). Both quantities were inversely correlated with lake surface area. CH4 587 concentrations in the stream feeding the Mire (22.2 ± 5.1 mg m −3 , n = 29, mean ± 95% CI), were significantly 588 higher than those in the lakes (Table 2). Surface water concentrations over the deep parts of the deeper 589 lakes (≥ 2 m water depth) were lower than those in the shallows (< 2 m) by 21 to 26% for Inre and Mellersta 590 Harrsjön, respectively. However, the diffusive CH4 flux did not differ significantly between depth zones in Relations between the flux and its drivers -temperature, wind speed and the surface concentration -597 manifested on different timescales (Fig. 7). Over the ice-free season both the CH4 fluxes and surface water 598 concentrations tracked changes in the water temperature. The wind speed (U10) showed less variability 599 over seasonal (CV = 7%, n = 17) than over diel timescales (CV = 12%, n = 24) and displayed a clear diurnal 600 maximum. The surface water/sediment temperature varied primarily on a seasonal timescale (CV = 601 52%/45%, n = 17), and less on diel timescales (CV = 3%/2%, n = 24). Similar to the wind speed the gas 602 transfer velocity varied primarily on diel timescales (Fig. 7), albeit with a lower amplitude. This was in part 603 because ∝ 3 4 ⁄ (Eq. 4), and because the drag coefficient, used to compute the water-side friction 604 velocity in Equation 5, increases at lower wind speeds and under an unstable atmosphere, which was 605 typically the case. The surface concentration correlated with wind speed and temperature (Fig. 4f,g), and 606 showed both seasonal and diel variability. On diel timescales Δ[CH4] and kmod were out of phase; Δ[CH4] 607 peaked just before noon, when the gas transfer velocity reached its maximum value (Fig. 7b,d). However, 608 binned means of the 1-hour chamber fluxes (Fch (1h)) were not significantly different at the 95% confidence 609 level (error bars) and did not show a clear diel pattern (Fig. 7b). Temporal patterns of fluxes and 610 concentrations were very similar between the lakes (Supplementary Fig. 2 and 3). shown in Supplementary Figures 2 and 3.

Timescale analysis 632
The spectral density plot (Fig. 8a) disentangles dominant timescales of variability of the drivers of the flux. 633 The power spectra of wind speed and temperature peaked at periods of 1 day and 1 year, following well-634 known diel and annual cycles of insolation and seasonal variations in climate (Baldocchi et al., 2001). The 635 diel spectral peak was subdued for the surface sediment temperature. For U10, the overall spectral density 636 maximum between 1 day and 1 week, and somewhat longer in spectra for the ice-free period only 637 ( Supplementary Fig. 4), corresponds to synoptic-scale weather variability, such as the passage of fronts 638 (MacIntyre et al., 2009). U10 and Tair also exhibit spectral density peaks at 1-3 weeks, which could be 639 associated with persistent atmospheric blocking typical of the Scandinavian region (Tyrlis and Hoskins, 640 2008). While the temperature variability was concentrated at annual timescales, the wind speed varied 641 primarily on timescales shorter than about a month and often shorter than a week. 642 643 The climacogram (Fig. 8b) reveals that the variability of the chamber flux and the gas transfer velocity was 644 enveloped by that of the water temperature and the wind speed, as was the surface concentration 645 difference for timescales < 5 months. The distribution of variability over the different timescales is similar 646 to that shown in the spectral density plot (Fig. 8a). The standard deviation of the water temperature did 647 not change from its initial value (σ/σinit = 1) until timescales of about 1 month, following the 1 year 648 harmonic. In contrast, most of the variability of the wind speed was concentrated at time scales shorter 649 than 1 month. The variability of the chamber and modelled fluxes first tracked that of the wind speed, but 650 for timescales longer than about 1 month the decrease in variability resembled that of water temperature. 651 The variability of the modelled fluxes followed that of the surface concentration difference rather than the 652 gas transfer velocity. However, the coarse sampling resolution of the fluxes and concentrations may have 653 led to an underestimation of both the variability at <1-week timescales (Fig. 7b) and the value of σinit. 654 Finally, the climacogram shows that kmod retains about 72% of its variability at 24-hour timescales, which 655 justifies our averaging over chamber deployment periods for comparison with kch and the computation of 656 the model scaling parameter α' (Fig. 3). 657  Mire lakes do not appear to be limited by substrate quality or quantity (Wik et al., 2018), but strongly 675 depend on temperature (Fig. 4b) The gas transfer velocity in the Stordalen Mire lakes was similar to that predicted from wind-based models 681 of Cole and Caraco (1998) and Crusius and Wanninkhof (2003) at low wind speeds (Fig. 9). Both were based 682 on tracer experiments with sampling over several days, and thus, like our approach, are integrative 683 measures. The slope of the linear wind-kch relation (OLS: 0.81 ± 0.21, slope ± 95% CI, R 2 = 0.20 and p < 0.01 684 for the individual kch estimates (small orange rhombuses in Fig. 9)) was similar to that reported by Soumis spectrum of water bodies in which the gas transfer models in Fig. 9 were developed (Table S1)

Drivers of flux 715
Methane emitted from lakes in wetland environments can be produced in situ, or be transported in from 716 the surrounding landscape (Paytan et al., 2015). The distinction is important because some controls on 717 terrestrial methane production, such as water table depth (Brown et al., 2014), are irrelevant in lakes. In 718 the Stordalen Mire lakes, the Arrhenius-type relation of CH4 fluxes and concentrations (Fig. 4b,f) together 719 with short CH4 residence times (Fig. 6) suggest that efficient redistribution of dissolved CH4 strongly 720 coupled emissions to sediment methane production. High CH4 concentrations in the stream (section 3.4) 721 further suggest that external inputs of CH4 -produced in the fens and transported into the stream with The distinct spectral peaks of temperature and U10 (Fig. 8) indicate that flux dependencies on these 747 parameters (Fig. 4b,c)  would resolve most of the variability of the ice-free diffusive CH4 flux at timescales longer than a month. 751 Advanced gas transfer models that account for atmospheric stability and rapid variations in wind shear, 752 such as we have used here, allowed us to resolve variability in flux at timescales shorter than about a 753 month. Our results are representative of small, wind-exposed lakes in cold environments, where, as a 754 result of considerable wind driven mixing, fluxes are lower than would be predicted in lakes where 755 buoyancy fluxes during heating and cooling are higher. 756

Storage and stability 757
The robust temperature-sensitivity of lake methane emissions (Fig. 4b,f)  by periodically decoupling production from emission rates (Engle and Melack, 2000). Here, enhanced CH4 760 accumulation during periods of stratification may have contributed to concentration and storage maxima 761 in July and August (Fig. 4e, 6d). However, as the CH4 residence time was invariant over the season and with 762 temperature (Fig. 6a,b), the storage-temperature relation (Fig. 6e) likely reflects rate changes in sediment 763 methanogenesis rather than inhibited mixing. For example, the highest CH4 concentrations in our dataset 764 (59.1 ± 26.4 mg m −3 , n = 37) were measured during a period with exceptionally high surface water 765 temperatures (Twater = 18.5 ± 3.6 °C) that lasted from 23 June to 30 July 2014. Emissions during this period 766 comprised 29%-56% (depending on lake) of the 2014 ice-free diffusive flux, while the peak quantity of 767 accumulated CH4 comprised <5%. Two mechanisms may explain the lack of CH4 accumulation. First, 768 stratification was frequently disrupted by vertical mixing (Fig. 5g-h)  Thus, in weakly stratified lakes with strong wind mixing, the temperature sensitivity of diffusive CH4 775 emissions may be observed without significant modulation by stratification. 776 777 Degassing (Fig. 4c,g) prevented an unlimited increase of the emission rate with the gas transfer velocity. 778 In this way, Δ[CH4] acted as a negative feedback that maintained a quasi-steady state between CH4 779 production and removal processes throughout the ice-free season. In all three lakes CH4 residence times 780 were inversely proportional to the wind speed (Fig. 6c), indicating an imbalance between production and 781 removal processes. We hypothesize that the imbalance exists because the variability of wind speed peaked 782 on shorter timescales than that of the water temperature (Fig. 8a). Changes in wind shear periodically 783 pushed the system out of production-emission equilibrium, allowing for transient degassing and 784 accumulation of dissolved CH4. The temporal variability of dissolved gas concentrations is likely higher in 785 shallow wind-exposed systems with limited buffer capacity (Natchimuthu et al., 2016(Natchimuthu et al., , 2017, and should 786 be taken into account when applying gas transfer models to small lakes and ponds.

788
Rapid degassing occurred at U10 ≥ 6.5 m s −1 (Fig. 4c). Gas fluxes at high wind speeds may have been 789 enhanced by the kinetic action of breaking waves (Terray et al., 1996) or through microbubble-mediated 790 transfer. Wave breaking was observed on the Stordalen lakes at wind speeds ≥ 7 m s −1 . Microbubbles of 791 atmospheric gas (diameter < 1 mm) can form due to photosynthesis, rain or wave breaking (Woolf and 792 Thorpe, 1991) and remain entrained for several days (Turner, 1961). Due to their relatively large surface 793 area they quickly equilibrate with sparingly soluble gases in the water column, providing an efficient 794 emission pathway to the atmosphere when the bubbles rise to the surface (Merlivat and Memery, 1983).

795
In inland waters microbubble emissions of CH4 have only been indirectly inferred from differences in CO2 796 and CH4 gas transfer velocities (McGinnis et al., 2015; Prairie and del Giorgio, 2013), and more work is 797 needed to evaluate their significance in relatively sheltered systems. 798

Timescales of variability 799
Overall, the short-term variability of the flux due to wind speed (1.1-13.2 mg m −2 d −1 ) was similar to the 800 long-term variability due to temperature (0.7-12.2 mg m −2 d −1 ) (ranges of the binned means, Fig. 4b-c). 801 The diel patterns in the mixed layer depth (Fig. 5) and the gas transfer velocity (Fig. 7d) and daytime 802 variation of the surface concentration (Fig. 7b) were indicative of daily storage-and-release cycles, 803 resulting in a model flux difference of about 5 mg m −2 d −1 between morning and afternoon; about half the 804 mean seasonal range (Fig. 7a). Diel variability of lake methane fluxes has been observed at Villasjön (eddy in model studies (Erkkilä et al., 2018). Apparent offsets between the diurnal peaks of the flux, surface 809 concentrations and drivers (Fig 7b,d) have been noted previously (Koebsch et al., 2015), but have yet to 810 be explained. Continuous eddy covariance measurements in lakes where the dominant emission pathway 811 is turbulence-driven diffusion could help characterize flux variability on short timescales (e.g. Bartosiewicz 812 et al., 2015).

814
The CH4 residence times (1-3 days) were not much longer than the diel timescale of vertical mixing (Fig.  815 5g,h). As a result, horizontal concentration gradients developed in the deeper lakes ( Table 2). The 23 ± 816 11% concentration difference between depth zones in the deeper lakes (mean ± 95%) fits transport model suggest that weekly sampling did not capture the full temporal variability of the surface concentrations. 830 Especially after episodes of high wind speeds and lake degassing (Fig. 4c,g), concentrations may not have 831 been representative of the 24-hour chamber deployment period. 832

Model-chamber comparison 833
It is fundamental to our understanding of controls on fluxes to determine why empirically derived values 834 of the model scaling parameter α' are relatively low in this study (0.17-0.31) compared to the theoretical 835 value of √2 15 ⁄ ≅ 0.37 (Katul et al., 2018), and why they were different in the three lakes. Differences 836 in α' resulted from kch, with mean (± 95% CI) values estimated at 3.5 ± 0.7 (n = 74), 3.1 ± 0.4 (n = 131) and 837 2.5 ± 0.6 (n = 142) cm h −1 in Villasjön, Inre Harrsjön and Mellersta Harrsjön, respectively, while kmod did not 838 differ significantly between lakes (ANOVA, p < 0.001). Synthesis studies show that scaling parameter values 839 can vary between 0.1 and 0.7 over the range of moderate to high dissipation rates computed for the wind sheltering, atmospheric stability, and within lake stratification and mixing. Here, the low α' value may 850 imply an underestimation of k derived from chamber observations or an overestimation of dissipation 851 rates used in the modelling of gas transfer velocities.

853
An underestimation of chamber-derived gas transfer velocities may have resulted from an overestimation 854 of Caq in Equation 1. In most freshwater systems a significant fraction of CH4 is removed through microbial 855 oxidation (Bastviken et al., 2002). This additional removal process invalidates the implicit assumption in 856 Eq. 1 and 2 that all dissolved CH4 that we measure in the surface water is emitted to the atmosphere. 857 Omitting oxidation would bias Δ[CH4] high, and kch low. The Stordalen Mire lakes remained oxygenated 858 throughout the ice-free season and CH4 stable isotopes indicate that between 24% (Villasjön) and 60% 859 (Inre and Mellersta Harrsjön) of CH4 in the water column was continually oxidized (Jansen et al., 2019). 860 This may explain not only the low scaling parameter value compared to those found with other tracers, 861 but also why ′ was higher in Villasjön (0.31, n = 67) than in the deeper lakes (0.17-0.25, n = 267) 862 ( Supplementary Fig. 1). However, more work is needed to establish how the oxidation effect partitioned 863 between CH4 reservoirs in the water column, where it would affect surface emissions, and the sediment. 864 An increase in surface concentrations which typically occurs at night would not have been manifest ( intermittent in our study (Fig. 5f-h). 879 880 Reduced gas transfer velocities and between-lake differences in kch could also be due to differences in 881 atmospheric forcing. First, the wind speed may have been lower over the lakes than on the Mire due to  884 1) was readily noticed during sample collection, particularly in Mellersta Harrsjön. Second, atmospheric 885 stability was different over the three lakes. The atmosphere was stable (z/LMO,a > 0) over Mellersta 886 Harrsjön, Inre Harrsjön, and Villasjön during 29%, 21% and 22% of the ice free periods (2009-2017), 887 respectively, with drag coefficients ∼16% lower than their neutral value during these times. The effect was 888 more pronounced when winds were light during daytime heating, with somewhat higher frequency during 889 autumn. Colder incoming stream water flowing into Mellersta Harrsjön may have contributed to lower 890 surface water temperatures in this lake (Table 3), with the discrepancy more noticeable as lake level rose 891 (Fig. 5e-h). More frequent periods with a stable atmosphere above Mellersta Harrsjön reduced sensible 892 and latent heat fluxes and are a likely cause of the increased stratification of the surface layer: water at 893 0.1 m was sometimes 0.5 °C to 2 °C warmer than at 0.3 m in Mellersta Harrsjön (5% of the time during ice-894 free seasons) when temperatures were isothermal in the upper 0.5 m in Villasjön and Inre Harrsjön. 895 Greater near-surface stratification coupled with lower winds than measured on the Mire would have led 896 to the lower values of k and α' obtained in this lake. While this analysis points to the challenges in modelling 897 fluxes when meteorological instrumentation is not situated on the lakes, it also suggests that a solution is 898 to use lower values of α' when modelling k for sheltered water bodies. 899 900 In summary, the model scaling parameter ′ computed in this study are lower than the theoretical value 901 of 0.37 and the 0.5 recently obtained in eddy covariance studies in which fluxes were measured with CO2 902 and modelled. The discrepancy may be explained by surface CH4 concentrations decreasing due to 903 microbial oxidation over the same timescale as our chamber measurements. Alternate explanations take 904 into account the magnitude of wind shear and degree of sheltering. Differences in ′ between lakes 905 indicate the care required in modelling emissions from sheltered lakes; the overall cooler surface water 906 temperatures in the lake with greater stream inflows points to a new control on emissions. That is, when 907 stream inflows lead to surface water temperatures cooler than air temperature in sheltered lakes, a stable 908 atmosphere results which leads to a reduced momentum flux, lower emissions, and a longer time over 909 which methane oxidation can occur. The cooling effect may be especially pronounced in northern 910 landscapes underlain by permafrost, where the temperature of meltwater streams and subsurface flow in 911 the active layer remains low throughout the year. Thus, these comparisons of modelled and measured 912 fluxes point to new areas of research. 913

Summary and conclusions 914
In this study we combined a unique, multi-year dataset with a modelling approach to better understand 915 environmental controls on turbulence-driven diffusion-limited CH4 emissions from small, shallow lakes. 916 Floating chambers estimated the seasonal mean flux at 6.9 mg m −2 d −1 and illustrated how the flux 917 depended on temperature and wind speed. Wind shear controlled the gas transfer velocity while thermal 918 convection and release from storage were minor drivers of the flux. CH4 fluxes and surface concentrations 919 fitted an Arrhenius-type temperature function (Ea' = 0.88-0.97 eV), suggesting that emissions were 920 strongly coupled to rates of methanogenesis in the sediment. However, temperature was only an accurate 921 proxy of the flux on averaging timescales longer than a month. On shorter timescales wind-induced 922 variability in the gas transfer velocity, mixing layer depth, and storage decoupled production from emission 923 rates. Transient changes in the lake mixing regime allowed for periodic CH4 accumulation and resulted in 924 an inverse relationship between wind speed and surface concentrations. In this way, the air-water 925 concentration difference acted as a negative feedback to emissions and prevented complete degassing of 926 the lakes, except at high wind speeds (U10 ≥ 6.5 m s −1 ).

928
Freshwater flux studies are increasingly focused on understanding mechanisms and developing proxies for 929 use in upscaling efforts and process-based models. Simple temperature-or wind-based proxies can yield 930 accurate flux estimates, but model parameters, such as Ea' and ′, must be calibrated to local conditions 931 to reflect relevant biotic and abiotic processes at appropriate timescales. Our study highlights the 932 importance of non-linear feedbacks, such as shallow lake degassing at high wind speeds, as well as 933 microbial removal processes and the need to consider the timescale over which fluxes occur relative to 934 the timescale over which CH4 can be oxidized. Such biological removal processes may violate the 935 fundamental assumption of gas transfer models that all gas measured in the surface mixing layer is emitted 936 to the atmosphere. Advanced gas transfer models can only improve the accuracy of flux estimates if they 937 are paired with observations that capture the meteorological conditions over the lake and the 938 spatiotemporal variability of dissolved gas concentrations. Therefore, field measurements remain 939 necessary to inform, calibrate and validate models. Our results indicate that the timescale of driver 940 variability can inform the frequency of field measurements necessary to yield representative datasets for 941 novel proxy development.