Measurements were made on Toolik Lake (68°37.830′ N, 149°36.366′ W, WGS84 datum) at 719 m asl, and in the Toolik Field Station (TFS) environment. Toolik Lake is a relatively deep glacial lake (maximum depth ≈26 m) located on the tundra north of the Alaskan Brooks Range with a surface area of 1.5 km² (Hobbie and Kling, 2014).

Six ice-free seasons (2010–2015) were covered with eddy covariance flux measurements (Section 2.2), the results of which were presented by Eugster et al. (2020a). In addition, in 2015 an equilibrator system with a dissolved gas extraction unit was used to continuously determine the gas mixing ratios in the Toolik lake surface waters (Section 2.3). And finally, long-term observations were made as a component of the ARC LTER project since the late 1980s. The location of the TFS meteorological station (68°37.698′ N, 149°35.759′ W, 722 m asl) is ∼500 m southeast of the EC flux measurements.

Toolik lake is ice covered during 9–10 months of the year, with CO₂ and CH₄ accumulating under the ice (Kling et al., 1992), an aspect not covered here. Ice on and off dates are typically from mid-June to early July, and ice-on starts in the time period of end of August to late September. Because the Sun does not set at this northern latitude from 23 May to 19 July, there is no dark period of the day and “night” with darkness is not observed before August.

The instrumentation and methods used in this study were presented in detail in Eugster et al. (2020a). Here we briefly summarize the key aspects of the EC measurements. Flux measurements were made with a three-dimensional ultrasonic anemometer-thermometer (CSAT3, Campbell Scientific, Logan, UT, United States) in combination with a closed-path integrated off-axis cavity output spectrometer (ICOS) for CH₄ (FMA, Los Gatos Research, Inc., San Jose, CA, United States) and a nondispersive infra-red gas analyzer for CO₂ and H₂O (Li-7000, Li-Cor Inc., Lincoln, NE, United States). Instruments were mounted on a floating platform (Supplementary Figure S1A) that was moored at approximately the same location every year ≈400 m from the nearest lake shore. Depth of the lake at this location was ≈12 m. The gap-filling of missing data up to four hours (up to 8 hours for 3-hourly resolved variables) was done by linear interpolation. For longer gaps the following approach was used: 1) the median diel cycle of all available years during the time-of-day of the gap, including a 1-day margin on both sides, was used as a reference; 2) then the beginning of the gap-filling time series was connected to the endpoint of the available data using an offset correction; and 3) finally, the endpoint of the gap was connected to the first valid data point after the gap via linear scaling. In addition, the pyranometer data were offset-corrected to obtain a mean 0 W m⁻² nocturnal flux (this corrects for the offset associated with the instrument temperature of the pyranometer), after which all nocturnal values <0 W m⁻² were set to zero. The main precipitation time series only captured liquid precipitation, hence frozen precipitation in summer (which tends to be a small component) was neglected until a second, heated rain gauge was installed on 19 June 2004. There was very good agreement of the rainfall amounts measured by both gauges [major axis regression slope of 0.0998%, 95% CI 0.995–1.001; Legendre and Legendre (1998)]. Thus, for the analysis here, the gaps in the long-term precipitation time series were filled by inserting the measurements from the newer rain gauge, with no attempt to fill the remaining (mainly winter) gaps. Thus, in the case of precipitation, our analysis is biased towards the warm season with predominantly liquid precipitation. For details on flux calculation, gap-filling of missing data, and the flux footprint of the EC measurements the reader is referred to Eugster et al. (2020a).

During the 2015 ice-free season an equilibrator system with automatic dissolved gas extraction unit was deployed on the EC float (Supplementary Figure S1A). The device used was a modified version of a prototype that Los Gatos Research, Inc., produced for us in 2012–2013, from which the commercial LGR Dissolved Gas Extraction Unit (https://www.lgrinc.com/documents/LGR_Dissolved%20Gas%20Extraction%20System.pdf) evolved. The basic principle (Supplementary Figure S2) is to pull water through a cylindrical Liqui-Cel device containing a gas-permeable membrane separating the water flow from the gas flow. The under-pressure generated by pulling water through the Liqui-Cel device increases the efficiency with which dissolved gases pass the membrane and are incorporated into the air strip-flow section of the device. The strip gas containing the gases extracted from the water are then analyzed by CH₄ and CO₂ analyzers (Supplementary Figures S1,S2), and the equilibrium gas mixing ratio (mol fraction or parts per million by volume) is then calculated using the flow rates of all components of the system, the initial gas mixing ratios in the air strip-flow, and the gas extraction efficiency of the Liqui-Cel membrane.

However, we were unable to achieve the expected accuracy of the equilibration mixing ratios of CO₂ and CH₄ in this configuration because we found no simple way to accurately determine the efficiency of the Liqui-Cel membrane. The reason might have been that Toolik lake has a substantial load of particulate organic matter (POM) in the water (DelSontro, 2011) and the main water filter system was unable to remove 100% of the POM before the water stream entered the Liqui-Cel device. We thus expected that the efficiency of the Liqui-Cel device was a non-trivial function of 1) its age, 2) the elapsed time since replacement of the main water filter, and 3) the POM load in the sample water stream. Thus, we remodeled our prototype to operate in conventional closed-loop equilibration mode (Supplementary Figure S2 for the configuration used in this study).

Equilibrium gas mixing ratios were obtained by scaling the raw mixing ratio measurements obtained from the equilibrator device with the net difference between the gas flow from the Liqui-Cel membrane to the analyzers and the strip gas flow measured with a gas flow meter, c e q = c m e a s u r e d ⋅ F g a s − F s t r i p g a s ⋅ c a m b i e n t F g a s − F s t r i p g a s (1) with ambient mixing ratio c_ambient of 390 ppmv (parts per million by volume) for CO₂ and 1.85 ppmv for CH₄. The gas flow to the analyzer, Fgas, was ≈5.5·10⁻⁴ m³ s⁻¹ (a 408 cm³ sample cell refreshed every 7.4 s). The strip gas flow Fstripgas was continuously measured by the dissolved gas extraction unit (Supplementary Figures S1B,C). Due to the high variability of the raw mixing ratio measurements, a local polynomial regression smoothing (loess function in R) was applied in the computation of c_measured to obtain the final equilibrium mixing ratio c_eq. Supplementary Figure S3 shows the data used for quality checking this procedure.

In 1988 long-term meteorological measurements started, and the variables used in this study are shown in Supplementary Table S1. From the available variables, air temperature, barometric pressure, short-wave incoming radiation, lake depth, lake temperature, precipitation, relative humidity (RH), soil temperatures, and wind speed were used for trend analysis. For pre-1996 temperature and RH at 5 m height, we gap-filled from a local station when available, and to correct a degraded RH sensor an offset-correction was applied that sets the moving 7-days maximum to 100% and trims data to 100%. Vapor pressure deficit (VPD) was calculated from the air temperature Ta (in °C), relative humidity RH (in %) and barometric pressure Pa (in hPa) using the procedure by Buck (1981) to yield V P D = ( 1 − R H 100 % ) ⋅ 6.1121 ⋅ e 17.502 ⋅ T a T a + 240.97 ⋅ ( 1.0007 + 3.46 ⋅ 10 − 6 ⋅ P a ) (2)

Snow depth and wind direction were not included in the trend analysis. Snow depth is only measured since 2014 (Supplementary Table S1) and moreover during the ice-free season an extensive snow cover is missing and thus it is not expected that this variable strongly influences summertime CO₂ and CH₄ fluxes. Wind direction also was excluded from the trend analysis because small trends in wind direction are not expected to influence lake–atmosphere GHG effluxes substantially.

For the trend analyses we used the Theil-Sen median slope estimator (Sen, 1968; Akritas et al., 1995) that takes account for the fact that time-series do not meet the model assumptions of ordinary least-square regression (the x-axis is not a random variable, and values are equally spaced in the time-series and do not follow a random normal distribution). This trend slope estimator was then combined with a quantile regression approach (Easterling, 1969; Lejeune and Sarda, 1988; Patel, 2021); within each year the quantiles of the available measurements are determined for each percentile of the empirical distribution for that year. Then the Theil-Sen median trend slope fit was calculated for each new time-series of the quantiles in each percentile of each year. This approach allows identification of trends that are not uniform across the entire gradient of values observed in a given time-series. To estimate the long-term overall trend of each monitored variable we then used only the percentiles with significant trends to obtain the median change per decade and its 95% confidence interval.

Statistical analyses were carried out with R version 4.1.2 (R Core Team, 2014). To compute the Theil-Sen median trend slope the trend package (Pohlert, 2020) was used, which provides the sens.slope function (Sen, 1968). Significance of the slope was determined using the nonparametric Mann-Kendall test provided by the mk.test function. Major axis (orthogonal) regression was computed using the lmodel2 function of the lmodel2 package (Legendre and Legendre, 1998).

A principal component analysis (PCA) of CO₂ flux measured during 2015 over Toolik Lake combined the potential driving variables at the same temporal resolution as CO₂ fluxes. This analysis includes meteorological variables, eddy flux variables, and temperature and gas mixing ratios in the water, including temperature of the upper mixing layer in the lake (0–4 m depth) and temperature near or below the typical summer thermocline (5 m depth). If arrows are shown perpendicular to each other, the variables are independent (uncorrelated) in the selected PCA axes. Figure 4A shows CO₂ flux as the dominant variable in the second PCA axis. The two first PCA axes in Figure 4A explain a total of 59.5% of the variance in the selected dataset (42.4 % and 17.1%). F_CO2 increases as atmospheric CO₂ density (in mmol m⁻³) increases, but is almost unrelated to CO₂ mixing ratio corrected for variability in density of air, which includes barometric pressure, air temperature, and atmospheric water vapor content. In contrast, the CO₂ mixing ratio gradient across the air-water interface, ΔCO₂ (in ppmv), has a weak, direct negative effect on F_CO2. The sign conventions of F_CO2 and ΔCO₂ are that F_CO2 is positive when gas evades from the water to the atmosphere, and ΔCO₂ is positive when the equilibration mixing ratio of CO₂ in the surface water is higher than the corresponding atmospheric value.

The transfer of CO₂-rich waters from depths lower in the epilimnion or the metalimnion toward the lake surface will enhance flux to the atmosphere. This is seen in Figure 4A as the first temperature below the typical thermocline (water temperature at 5 m depth) correlates positively with ΔCO₂, whereas the epilimnion temperatures down to 4 m depth are negatively correlated with ΔCO₂. There is a slight but important increase in correlation between the epilimnion temperatures and F_CO2 with increasing depth, indicating that temperature effects at the thermocline where CO₂ is typically higher (Supplementary Figure S5) are more important for F_CO2 than is the surface temperature.

Increasing wind speed (U) and turbulence (u_*, u′², v′², w′²) all reduce—not increase—F_CO2, almost in perfect agreement with net radiation, which represents atmospheric stability over the water surface. With increases in net radiation, wind speed and turbulence also increase. Atmospheric stability is negative (z/L < 0 and Ri < 0; e.g., Stull, 1988 and Supplementary Figure S6) most times of the day, except in the afternoon when it is in the near-neutral to slightly unstable range (ca. 1200–2100 h AKDT; data not shown). This indicates the important difference of a lake surface from a terrestrial surface, where over land atmospheric stability is negative when net radiation is high, but positive (stable) at night during times where the lake surface water temperatures are warmer than the terrestrial atmosphere. These conditions where a cooler atmosphere overlies a warmer surface water, leading to convection in the water, lead to the substantial peak in F_CO2 at night between 0000 and 0600 h AKDT (Figure 2B), whereas barometric pressure has no influence on F_CO2. In summary, CO₂ effluxes from Toolik Lake are strongly governed by environmental drivers and thus are a key component in the second PCA axis, indicating a close link between F_CO2 and the observed driver variables.

In contrast to CO₂ efflux from Toolik Lake, the CH₄ fluxes are much less coupled with the same driver variables (Figure 4B) that influence CO₂ evasion (Figure 4A). The first two PC axes (explaining 59.5% of the data subset used in this analysis) do not include F_CH4, thus we only inspect PC axes 3 and 4 which explain 22.8% and 14.3% of the remaining variance, respectively (Figure 4B). For F_CH4 the wind direction has a strong influence with northerly directions (direction towards 360°) reducing F_CH4. This reflects the local daytime foreland–mountain wind system with northwesterly winds (from the coast to the Brooks range) typically prevailing from 0800 to 1400 h AKDT, and south-easterly dominance (winds from the mountains to the foreland of the North Slope) from 1500 to 0600 h AKDT.

However, in contrast to CO₂, equilibrium methane mixing ratios in the water are supersaturated at all times (Figure 3) by a median 27.3 ppmv (95% CI 13.2–51.3 ppmv), and thus both atmospheric CH₄ mixing ratio and the gradient measured acros s the water interface have almost no influence on F_CH4 (the two arrow are almost at a right angle with F_CH4 in Figure 4B, indicating statistical independence). However, not quite as expected, air temperature (represented by three redundant measurements) has a negative effect on F_CH4 with warmer air temperatures reducing CH₄ evasion from the lake. According to Figure 4B warmer air temperatures are associated with lower wind speed and atmospheric turbulence. Also barometric pressure has an influence as high pressure tends to suppress F_CH4. Higher wind speeds and turbulence only have a weak enhancing effect on F_CH4, but it is clearly in contrast to F_CO2 where the same variables seem to reduce F_CO2. Net radiation (and thus atmospheric stability) and water temperatures of the epilimnion were unrelated with F_CH4; only the upper metalimnion temperature (5 m depth) seems to positively influence F_CH4. In the biweekly sampled CH₄ mixing ratio profiles the peak CH₄ mixing ratio is often found at the depth of the thermocline (Supplementary Figure S5), whereas in contrast to the substantial CO₂ storage in the metalimnion there is no such CH₄ storage in deeper waters of Toolik lake.

The diel variation of CH₄ mixing ratio in the water is substantial (Figure 3A), with the lowest median of 18.0 ppmv at 0400 h AKDT and the highest median of 37.6 ppmv at 1900 h AKDT. This range, however, does not translate to a clear diel cycle in F_CH4 (Figure 3B) because the water is strongly supersaturated with CH₄ at all times. Median fluxes during low and high water mixing ratio periods of the day are all on the order of 0.02–0.06 µg CH₄ m⁻² s⁻¹, with no strong pattern in diel variation (Figure 3B).

The CO₂ efflux from Toolik Lake is almost insensitive to mean horizontal wind speed (U) (see also Zappa et al., 2007) (Figure 5A; all regression fits to Toolik Lake data are given in Supplementary Table S2). It even appears that the highest observed wind speeds in 2015 actually reduced the CO₂ flux slightly (Figure 5A). If the data are analyzed separately for conditions with subsaturated (Figure 5C) vs. supersaturated (Figure 5D) CO₂ mixing ratios in the lake surface waters, and U is replaced by the friction velocity u_* (representing mechanical turbulence directly driving the gas flux), then CO₂ effluxes are highest under lower turbulence conditions (u_* < 0.2 m s⁻¹) when conditions are subsaturated, and become net neutral at u_* ≥ 0.2 m s⁻¹ (Figure 5C). In contrast, if the surface waters are supersaturated with CO₂, then CO₂ fluxes are similar to the subsaturated condition at low u_*, but become higher and more variable at u_* ≥ 0.2 m s⁻¹ (Figure 5C). This is consistent with surface waters that are supersaturated because increased stirring (via increasing u_*) should more completely outgas the surplus CO₂ in the surface waters.

A more direct and clear relationship of CO₂ efflux from Toolik lake was found when the temperature difference between the water surface and the air was used as a reference (Figure 6A). In the range where T_w is colder or less than 5°C warmer than the air, the CO₂ fluxes appear weakly driven by the temperature difference. But at T_w > (T_a + 5°C) an exponential increase of CO₂ efflux with increasing temperature difference was found (Figure 6A).

Methane fluxes do not show the same response to the temperature gradient between water and air (Figure 6B). If the water is 5°C colder than the air, the CH₄ effluxes were low and steadily increased to an optimum when the water temperature was around 8°C warmer than the air. If the difference was even larger, then CH₄ efflux decreased slightly (Figure 6B).

In this study we measured both concentrations in water and air, and fluxes of CO₂ and CH₄ fluxes (an advance since the beginning of air–sea exchange measurements, see Wanninkhof et al., 2009). Contrastingly, if only gas gradients are measured across the water–air interface, e.g., by syringe sampling, then the piston velocity normalized to a Schmidt number of 600 (k₆₀₀) must be estimated from horizontal wind speed measurements, although some authors claim that the relation to horizontal wind speed often is weak (Zappa et al., 2007). Most of the models used to predict k₆₀₀ from wind speed at 10 m height (U₁₀) are derived for open ocean environments with a large fetch without obstacles. Over an inland lake with limited surface area such as Toolik, these models of k₆₀₀ inadequately reflect the local conditions. Figure 7 shows that our direct measurements of the k₆₀₀ piston velocity for CO₂ are almost constant, the highest values occur at low U₁₀, and only at U₁₀ > 10 m s⁻¹ do our estimates fall in the range of k₆₀₀ that any model, including the MacIntyre et al. (2010) model that considers the buoyancy flux term, would predict for a stratified lake (Supplementary Table S3). This can be interpreted with the fact that it is actually u_*, not U₁₀, that drives vertical turbulent exchange of gases along the gradient between water surface and atmosphere. If measurements are made over a small lake, then the distance travelled over the lake surface is too short to bring u_* into equilibrium with U₁₀. There is a much higher level of turbulent mixing present in air masses that travelled over rough terrestrial tundra, and thus u_* remains relatively high in comparison to U₁₀.

The k₆₀₀ piston velocity for CH₄ differs from that of CO₂ and shows the expected exponential increase with increasing U₁₀ (Figure 8). However, there is a substantially higher turbulent exchange efficiency at low U₁₀ than expected from models, with k₆₀₀ around 14 cm h⁻¹ even at very calm wind conditions (Figure 8). At U₁₀ > 6 m s⁻¹ the k₆₀₀ for CH₄ then merges with the predicted values that most models provide (Figure 8).

The high turbulent exchange efficiency at low wind speeds, and the rather weak dependence of k₆₀₀ on U₁₀ in general, result in an almost constant CO₂ efflux (Figure 9A) and CH₄ efflux (Figure 9B) and mostly independent of the air-water difference in gas mixing ratio.

Figure 10A clearly shows that the u_*-to-U₁₀ relationship is substantially enhanced over Toolik Lake as compared to what would be expected over the open ocean (red dashed line in Figure 10A). At U₁₀ < 2 m s⁻¹, u_* is rather constant around 0.08 m s⁻¹, which is a factor 2–5× higher than what the Charnok relationship predicts (Figure 10A). Moreover, at low wind speeds the turbulent mixing not only depends on the mean U₁₀ during an average interval, but also on the history of U₁₀ during the previous interval. Over Toolik Lake the threshold of equal antecedent wind speeds is around U₁₀ ≈ 5 m s⁻¹ (Figure 10B); if U₁₀ is below that threshold, it is most likely that U₁₀ was higher with more turbulent mixing during the previous time interval, and if U₁₀ was above that threshold, the U₁₀ in the previous time interval was most likely lower.

Long-term trends are available for a few monitoring variables that cover the period from 1988 (or later) to 2020. For CO₂ and CH₄ gas exchange over Toolik Lake we focus on the long-term trends of lake depth (Figure 11A), lake surface water temperature (Figure 11B), horizontal wind speed (Figure 11C), rainfall (Figure 12), air temperature (Figure 13A), relative humidity (Figure 13B), short-wave incoming radiation (Figure 13C), barometric pressure (Figure 13D), and soil temperatures at the surface (in moss at the transition from green active tissue to brown peat) and at the deepest level recorded (1.5 m depth, corresponding to the maximum extension of the thawed active layer in summer) (Figure 14).

Because trends tend to be weak or statistically insignificant at annual aggregation of the long-term monitoring variables, all trends were determined via quantile regression that determines the trend slope for each quantile of the full dataset collected in each observation year. In this way, opposing trends at low, intermediate, or high values (quantiles) of the respective monitoring variable can be detected. Depending on the variable and expectation, we either tested the significance of the trend with a one-sided test (for existence of monotonically increasing or decreasing trends), or with a two-sided test to determine if any trend is significantly different from a null trend. Lake depth (Figure 11A) significantly decreased in all quantiles by −8.5 cm per decade, thereby reducing the water pressure on the lake bottom sediments which could increase the flux of gases produced at saturation in the sediments and released as bubbles to the water body (particularly the insoluble CH₄), and potentially thus increase the GHG efflux from organic sediments. However, we have observed very few ebullition events for CH₄ in Toolik Lake, which are the most expected events resulting from changing air pressure (Eugster et al., 2020a).

Lake water temperature did not show any significant trend during the warm season (Figure 11B), only winter conditions with surface water temperatures <2°C show a significant warming trend on the order of 0.1°C per decade (Figure 11B). Wind speed shows a decreasing trend at the lower 50% of the wind speed distribution, whereas wind speeds in the typical gentle breeze range (Beaufort 3, 3.6–5.1 m s⁻¹ according to WMO, Hasse and Isemer, 1986; Arguez and Vose, 2011) remained almost unchanged in the past 33 years. Only the 1% most extreme wind speeds have increased more substantially by 0.3 m s⁻¹ per decade, but this trend is not significant when we test for greater than normal wind speed extremes.

There is a highly significant increasing trend of rainfall events. On a daily basis the number of hours with precipitation has increased since 1988, both at the annual scale (Figure 12A) and during the terrestrial growing season months (May–September, Figure 12B). For days with up to 4 h of precipitation the increase in rainfall hours closely follows the + 5% increase curve during the growing season and follows the + 4% increase in the range of 4–7 h day⁻¹. The trend on days with >7 h of precipitation is less dramatic, although still statistically significantly compared to a null trend. Interestingly, the increasing duration of rain events is associated with significantly decreasing amounts of rainfall (−0.15 mm h⁻¹ per decade; Figure 12C). This combination of prolonged rainfall events with lower precipitation amounts may reduce the risk of terrestrial flushing or erosion due to strong storms, and reduce the terrestrial export of DOC, thus supplying the lake with less carbon to be respired to CO₂ and CH₄.

During the warm season no significant trends in air temperature were found (Figure 13A). Significantly increasing air temperatures only affect winter temperatures < −8°C (Figure 13A), but the trend is rather substantial with a 1.2°C increase per decade.

Trends of relative humidity measured at 5 m a.g.l (above ground level) were separately determined for temperatures above and below freezing (Figure 13B). Above freezing, the high relative humidity (>85%) showed no significant trends. This is in line with the observed rainfall trends with prolonged duration of rainfall events, which tend to dominate conditions with high relative humidity. Although statistically significant, the decreasing trend in relative humidity at humidity values < ∼80% is only a −1.6% reduction per decade, a value that most likely does not strongly affect GHG fluxes from Toolik lake.

Short-wave incoming radiation (Figure 13C) showed a roughly +0.5% increase over 33 years of measurements, which corresponds to a short-wave radiation increase of 9.0 [0.7–21.6] W m⁻² per decade. This decadal trend is 3.3 times the estimate of the human-caused radiative forcing in 2019 as compared to 1750 [2.72 (1.96–3.48) W m⁻²; IPCC AR6 WGI (2021; p.13: A.4.1)], but is in agreement with the global brightening that replaced the trend towards global dimming when the Toolik Lake long-term measurements were initiated in 1988 (Wild, 2005; Wild et al., 2012).

Barometric pressure has significantly increased at all quantiles by an average 1.4 [1.1–2.2] hPa per decade (Figure 13D), against a background range of ∼900–1,000 hPa. Converted to a water column pressure this corresponds to a 1.43 cm increase in water level per decade. This increase in barometric pressure counteracts the effect of the decreasing lake depth trend, but its magnitude is only one sixth of the concurrent decreasing trend in lake depth (Figure 11A).

Soil temperatures from the air-moss interface down to 1.5 m depth show a very consistent warming trend ranging from the lowest values at the soil surface (i.e., in the moss layer) of 0.7–1.0°C per decade (Figures 14A, B) to 1.1°C per decade at 1.5 m depth (Figures 14C, D), which is the transition zone from the maximum extent of the summer active-layer thaw and the permafrost. The warming trend is significant and consistent across all soil depths for temperatures < −1°C and > 1−5°C, but around the freezing point there is no significant trend seen at any depth (Figure 14). This lack of trend at the freezing point is regulated by the physics of the phase transition from solid ice to liquid water; the phase transition remains near 0.0°C irrespective of climate warming trends, but the response is how much heat content is stored in either the deeper permafrost soil or in the seasonally-thawed active layer.

Perhaps the most intriguing finding of this study is that measured piston velocities (k) for both CO₂ and CH₄ at wind speeds <10 m s⁻¹ exceed what flux models would predict (Figures 7, 8), and the friction velocities (u_*) at all wind speeds exceed the expected values given the Charnock relationship (Charnock, 1955; Figure 10A). Even updated Charnock relationships (Edson et al., 2013; Jimenez and Dudhia, 2018) that predict increased wind stress at a given U₁₀ still underestimate our measured u_*, especially at lower wind speeds (Figure 10A). Thus the higher friction velocity imparted to the water surface in our study should generate greater TKE in the surface water to enhance gas exchanges (k in Figures 7, 8).

There are several possible explanations for the high gas transfer velocities we measured at medium to low wind speeds (<10 m s⁻¹, Figures 7, 8). First, Munk (1947) proposed that flow over the ocean is laminar up to the Kelvin-Helmholtz instability around 6.5 m s⁻¹ and only becomes turbulent at higher wind speeds (Supplementary Figures S6, S7). Fetch is unlimited over the open ocean, and hence even at low wind speeds an equilibrium can establish between the atmospheric flow and the ocean depending on the previous status of the atmospheric and ocean turbulence (Supplementary Figure S6). Contrastingly, smaller lakes have a limited fetch and are thus more influenced by the surrounding, relatively rough landscape that generates more atmospheric turbulence, and the previous status is always “turbulent” under such conditions (Supplementary Figures S6, S7). There is no wind sheltering at Toolik Lake due to the low canopy structure of the surrounding tundra, but topographical rises of 20 m within 100 m of the shoreline are common and hills can be 30 m above the lake surface within 300 m of the lake shore. This complex terrain, coupled with a lack of wind sheltering to reduce wind speed and turbulence in the air mass reaching the lake (e.g., Markfort et al., 2010), act to generate more TKE in the atmosphere over the lake than is typically generated at the same wind speeds in smooth landscapes or over large lakes or oceans.

Second, the memory effect of turbulent mixing in the atmosphere brings well-mixed conditions from the rough terrestrial surface to the smooth lake water body, but in this small lake even a kilometer of fetch is unlikely sufficient to bring the turbulence state of the near-surface atmosphere into equilibrium with the water surface. However, at the same time the “young” waves generated when the fetch is small are typically shorter and steeper than waves closer to equilibrium with wind speeds in large lakes or oceans, and this steepness can increase the gas exchange (Edson et al., 2013) and increase the form stress on the water surface (Sullivan et al., 2018).

Third, the smaller, shorter waves that develop over shorter-fetch waters for a given wind speed are less likely to break than are larger waves and thus less likely to inject bubbles into the surface water. Bubble-mediated gas exchange is biased toward invasion (Woolf and Thorpe, 1991), and bubbles generated during wave breaking reintroduce atmosphere into the near surface and that reduces gas exchange (e.g., Emerson and Bushinsky, 2016). Note that as wind speed and wave height increase the bubble-mediated gas invasion will also increase, thereby reducing our measured k values compared to model predictions; this behavior is seen in Figures 7, 8.

Fourth, open-water surface roughness tends to increase at shallower depths (e.g., Taylor and Yelland, 2001), which leads to increased surface drag over shallow waters typically unaccounted for in models (Jiménez and Dudhia, 2018). Taken together, these processes tend to increase the measured k values compared to values predicted from models, including those developed for use in lakes and that incorporate buoyancy flux (e.g., MacIntyre et al., 2010). At present, we have no clear means of separating the individual contributions of these potential explanations for the higher than predicted gas exchange velocities we measured.

It is also unclear at present why there is a difference between the flux behavior of CO₂ and CH₄ (e.g., Figures 5–9). Although both of these gases are low enough in solubility to be generally controlled by water-side dynamics (Wanninkhof et al., 2009; Garbe et al., 2014), recent work has highlighted that bubble-dynamics vary between gases (Goddijn-Murphy et al., 2016; Rosentreter et al., 2017), and even gas solubility and transfer is differentially affected in response to films on the water surface (Mesarchaki et al., 2015).

A PCA indicated a strong coupling of CO₂ fluxes with environmental variables (meteorological variables, eddy flux variables, and temperature and gas mixing ratios in the water, including temperature of the upper mixing layer in the lake, 0–4 m depth, and temperature near or below the typical summer thermocline at 5 m depth). An increased dissolved CO₂ mixing ratio does not automatically increase the CO₂ efflux from the lake. It seems to represent the important transfer of CO₂-rich waters at the bottom of the epilimnion or in the metalimnion (Supplementary Figure S5) being transferred into the upper epilimnion, thereby increasing CO₂ mixing ratio, but not directly F_CO2 (Figure 4).

The CH₄ fluxes are much less coupled with the same driver variables (Figure 4B) that influence CO₂ evasion (Figure 4A). We only inspected PC axes 3 and 4 where F_CH4 plays a role and which explain 22.8% and 14.3% of the remaining variance (after PC1 and PC2), respectively (Figure 3B). It is probably the constantly high supersaturation of CH₄ (Figure 3A) that leads to a constant CH₄ efflux with a minor diel cycle (Figure 3B), and this makes F_CH4 much less dependent on environmental variables than is F_CO2 (Figure 4A).

The process-level relationships between GHG mixing ratios in the surface waters and the eddy-covariance flux measurements of CO₂ and CH₄ were only available from summer 2015, but seasonal flux measurements were carried out at Toolik lake during the ice-free seasons 2010–2015. In Eugster et al. (2020a) we provided a detailed comparison of interannual variability of the seasonal fluxes. In summary, the 2012 season provided the highest CO₂ and CH₄ effluxes with 1.3 g CO₂ m⁻² day⁻¹ (30 mmol m⁻² day⁻¹ and 3.2 mg CH₄ m⁻² day⁻¹ (0.20 mmol m⁻² day⁻¹), respectively. During the 2010 season with the lowest CO₂ efflux only 30% of the 2015 efflux was observed, and in the case of CH₄ the following 2011 season had the lowest efflux which was 35% of the 2015 maximum flux. This high interannual variability is not easily explained by the long-term trends of the monitored variables that influence gas exchange across the water surface in the short term and at the process-level. However, given that carbon input (CO₂, CH₄, DOC, and POC) to the lake from streams can be quite variable and can affect microbial processing of DOC to CO₂ or CH₄ (e.g., Kling et al., 2000; Crump et al., 2003, Crump et al., 2007), it is possible that variable inputs of carbon to the lake affect the interannual variability more than meteorological variables that govern gas exchange. Wu et al. (2013) have also reported that diffusive gas fluxes from two lakes were primarily attributable (30%–45%) to inputs and respiration of terrestrial DOC.

A linear, mixed-effect model combining the Eugster et al. (2020a) data with the monitoring variables used in this paper, suggests different processes governing the variability of the CO₂ fluxes (Table 1) versus the CH₄ fluxes (Table 2). The mixed-effect model shows that warm soil temperatures at 10 cm depth increase CO₂ fluxes by 26.5 ± 3.3% (Table 1), but rainfall 12 h ahead reduces CO₂ fluxes substantially by −37.7% ± 15.6%. It is possible that at warm times or in warm years that soil temperatures tend to increase organic matter decomposition in soil and thus increase the transport of CO₂ and CH₄ to the lake, or to increase the production of DOC and POC that feed into Toolik Lake. But rainfall, especially heavier rainfall intensities, are expected to increase effluxes of CO₂ and CH₄ from the lake (Guérin et al., 2007). Stronger winds, and stronger winds in combination with warmer 10-cm soil temperatures, also exert a negative effect on CO₂ fluxes on the order of −10% (Table 1). From other ecosystems it is also known that the near-surface soils are the most important substrate for respiration, in combination with soil temperature but with a dominance of availability of organic matter (e.g., Robinson et al., 2022). In arctic tundra there is substantial organic peat available across the upper soil profile, and thus the importance of the 10-cm conditions (and not the moss or 5 cm temperature) for respiration and decomposition may be most relevant as soils warm and eventually thaw in the Arctic.

The driving forces for CH₄ effluxes include a set of variables that each contribute only a minor percentage to the flux: warmer moss temperatures (i.e., temperature measured at the transition from green living moss matter to brown moss peat), a dry atmosphere (higher vapor pressure deficits, VPD), and warm temperatures in the topsoil (5 cm) in combination with warmer temperatures near the bottom of the active layer (100 cm). A strong negative effect is quite prominent when we observe turbulent mixing above the water surface (u_*) (−23.0 ± 0.5%, Table 2) or high u_* in combination with rainfall intensity 12 h ahead (−11.3% ± 3.1%, Table 2). That the processes affecting CH₄ fluxes differ from those affecting CO₂ is not surprising given the differences in relationships with physical forcing between the two gases (e.g., Figures 4–8). In addition, CH₄ is only produced under anaerobic conditions and will be oxidized to CO₂ when it is exposed to an aerobic environment.

To translate the past and present functional relationships between CO₂ and CH₄ fluxes and their environmental driving forces into the future, we analyze changes in the observed long-term monitoring variables and how they correlate with the GHG fluxes in Toolik Lake. Because the trends of most variables are looking roughly 30 years into the past, our projection is an estimate for conditions 30 years into the future (i.e., a time horizon of 2050). This extrapolation of trends over the last 30 years assumes that, in general, human activities leading to GHG emissions today will continue to increase at the same pace as in the last 30 years (the business as usual or RCP8.5 scenarios). This assumption is made in part because downscaled, regional predictions of future climate for the North Slope of Alaska are unavailable (see Hobbie and Kling, 2014). But projecting the trend observed during the past 30 + years to the next 30 years is not the same as using an IPCC scenario for predictions. Thus, with continued warming the IPCC projections expect that both CO₂ and CH₄ emissions to the atmosphere will increase in the years to come, despite the low confidence (or too little data) in the potential role of Arctic warming (IPCC AR6 WGI, 2021).

In summary, the relevant long-term trends in monitored potential driver variables of CO₂ and CH₄ fluxes are 1) the decreasing trend in lake depth, 2) the increasing trend in barometric pressure, 3) the decreasing trend of the lower 50% of wind speeds, 4) the 4%–5% increase in daily hours with rainfall in combination with the decreasing trend of rainfall intensity (precipitation amount per hour), 5) the increasing global radiation with global brightening since the early 1990s, and 6) the increasing soil temperatures above and below freezing (but not near the freezing point) across the entire 1.5 m depth profile equipped with temperature sensors. Not included in this discussion are 1) photosynthesis, 2) air temperatures that only show significant warming during the winter and thus outside the ice-free period covered with GHG flux measurements in this study, and 3) relative humidity that does not appear to have a strong influence on GHG effluxes from Toolik Lake. Wind direction trends were not included in this analysis (see Eugster et al., 2020a for wind direction analyses).

Starting with substrate availability for CO₂ and CH₄ production in Toolik Lake, the warming trend of the soil and increasing permafrost thaw (Turetsky et al., 2020) is probably an important positive feedback that will increase carbon inputs from land via stream inflow, groundwater, and overland flow to lakes in permafrost terrain (Hobbie and Kling, 2014; Vonk et al., 2015). If part of this increased input of C is in the form of particulate C that settles to the lake bottom, decreasing trends of lake depth might then amplify the outgassing of carbon gases produced from this particulate organic matter (less pressure to keep insoluble gases in the sediments). However, our earlier study (Eugster et al., 2020a) found no relevant role of CO₂ and CH₄ production in the bottom sediments of Toolik Lake, as deduced from the lack of evidence for ebullition. Moreover, the cold hypolimnion temperatures around 5–6°C even in peak summer limit biological activity and thus decomposition of C-rich substrates deposited in the lake bottom sediments. A study by Bayer et al. (2019) suggests that future gas fluxes from permafrost lakes are dependent on carbon inputs from the catchment.

The increasing trend in barometric pressure counteracts the effect of the decreasing trend in lake depth; that is, as the former reduces outgassing of GHGs from the bottom sediments, the latter promotes outgassing when the hydrostatic pressure on GHGs produced in lake sediments decreases. The net result would be toward enhanced outgassing because the magnitude of the barometric pressure trend is only one sixth of the opposing pressure trend exerted by lake depth. Hence we expect that a long-term trend of increasing barometric pressure may influence the short-term timing of GHG efflux from the lake, but the more important lake depth change most likely influences any long-term change in GHG fluxes derived from gases lost from lake bottom sediments.

Because of the weak relationship with CO₂ flux and either U₁₀ or u_*, a reduction of wind speeds below 1 m s⁻¹ or above 5 m s⁻¹ only slightly reduces the average F_CO2 level of 2.5 µg CO₂ m⁻² s⁻¹; at the highest wind speeds there is one data point showing that the CO₂ flux is reduced to near zero (Figure 5A). An assessment of the trend of increasing wind speeds on global ocean gas exchange with the atmosphere (Wanninkhof and Trenanes, 2017) suggested a slight increase in flux to the atmosphere, and more pronounced effects at lower wind speeds. That study did not consider gas fluxes from inland waters, which would generally experience lower wind speeds than do oceans. However, our results indicate that wind speed may be a poor predictor of and underestimate gas effluxes, which could be considered in future predictions of gas fluxes from all surface waters.

The unexpected increasing trend of daily hours with precipitation but decreasing trend in rainfall intensity (mm h⁻¹) indicates that the risk of storm erosion and thus transport of particulate organic matter to the inlet stream and further to Toolik Lake should decrease compared to the past. This could impose an important negative feedback on substrate availability for heterotrophic respiration in Toolik Lake that could even reduce the GHG fluxes from the lake. However, this potential effect must be balanced against the possibility of erosion and carbon inputs to the lake due to thermokarst failures, which are increasing in many arctic areas most likely due to permafrost thaw (Lewkowicz and Way, 2019). In addition, Cherry et al. (2014) predicted a wetter future for the Toolik region, which could increase the amount of DOC exported from land and respired within the lake (Kling et al., 2000, Kling et al., 2014). Finally, rainfall can enhance gas flux when raindrops add mixing energy to the water surface (Ho et al., 2004; Harrison et al., 2010), and longer periods of rainfall, independent of rain intensity, may thus increase gas fluxes. Overall, the effect of changes in precipitation on lake gas fluxes are complex and difficult to predict.

Increasing short-wave incoming radiation of course means more energy available at the Earth’s surface, which reduces atmospheric stability of the near-surface boundary layer and increases turbulent mixing in the atmosphere and thus the transfer of energy to the water surface that affects gas exchange (Erickson, 1993). At the same time, light-energy inputs warm surface waters, and increase buoyancy and water column stability, thus reducing TKE that can drive gas exchange across the lake surface (Imberger, 1985). On land, some of the energy will be used to increase soil heat flux; this is not in the set of long-term monitoring variables, but is key for the heat required to sustain the observed warming trends both below and above freezing in soils as discussed above. At present it is unclear how light energy increases will affect the balance between higher turbulence on land versus lower turbulence in surface waters with respect to controlling gas exchange.

The brightening of the atmosphere (Wild, 2005) reflected by the significantly increasing trend of short-wave incoming radiation has the potential to increase CO₂ fluxes, but most likely has no substantial effect on CH₄ fluxes. In case of the CO₂ fluxes we observed a positive influence of large T_w to T_a differences (Figure 6A), which most likely would be further amplified by increasing radiation inputs. In addition, increased photon flux to surface waters will increase photomineralization of DOC to CO₂ (Cory et al., 2014; Cory and Kling, 2018). Although greater solar radiation has the potential to increase algal photosynthesis, it is unlikely to be important in ultra-oligotrophic arctic lakes (Supplementary Figure S4). However, the warming trend in lake surface temperature (Figure 11B) was not significant at summer temperatures, but a significant warming trend of 0.1°C per decade under the ice cover (water temperature <2°C; Figure 11B) has the potential to increase under-ice microbial respiration, leading to a greater release of GHG accumulating under the ice during winter. However, this warming trend is quite weak and therefore unlikely to strongly impact microbial respiration under ice in the next 30 years.

Given the combined and often complex and counterbalancing effects of changing environmental drivers of gas flux, it is unlikely that Toolik Lake GHG fluxes will reach a tipping point where the system changes state from a source of gas to the atmosphere to a gas sink. The trends in drivers appear to be statistically significant but of low magnitude, most likely modifying present-day CO₂ and CH₄ effluxes in the range of an estimated factor 0.9–1.3. The lower value (0.9) is derived from an average reduction of –0.14 m s⁻¹ in mean wind speed of 3.0 m s⁻¹ per decade (–0.42 m s⁻¹ in 30 years), which translates to almost unchanged F_CO2 and a −10% change in F_CH4 (based on Supplementary Table S2, Figures 5A,5B). The upper estimate (1.3) is derived from changes in the terrestrial environment that may stimulate organic matter respiration and transport to surface waters; for example, moss temperatures (Figure 14A) increasing by 0.7°C per decade (2.1°C in 30 years) above the average temperature of 7.0°C for conditions >0°C. If this temperature increase is linearly related to the inflow of increased CO₂ and CH₄ to the lake as permafrost thaws (and not including DOC transport and subsequent oxidation), we expect a −30% increase in fluxes from the lake.

It is unclear how the identified trends in environmental drivers and the estimated impacts on GHG fluxes apply to other lakes. However, there may be some similarities in response because many of the processes that influence gas fluxes described here (e.g., rainfall and carbon loading, temperature and respiration, photomineralization of DOC) occur in all surface waters. We suggest that lakes deep enough to stratify seasonally may behave similarly to Toolik Lake, while shallower lakes or wetlands may respond differently to environmental change (e.g., greater ebullition than observed in Toolik). In addition, catchment differences in vegetation and surface topography may contribute to highly site-specific responses of lakes to changing environmental conditions.