Seasonal increase of methane emissions linked to warming in Siberian tundra

While increasing methane emissions from thawing permafrost are anticipated to be a major climate feedback, no observational evidence for such an increase has previously been documented in the literature. Here we report a trend of increasing methane emissions for the early summer months of June and July at a permafrost site in the Lena River Delta, on the basis of the longest set of eddy covariance methane flux data in the Arctic. Along with a strong air temperature rise of 0.3 ± 0.1 °C yr−1 in June, which corresponds to an earlier warming of 11 d, the methane emissions in June and July have increased by roughly 1.9 ± 0.7% yr−1 since 2004. Although the tundra’s maximum source strength in August has not yet changed, this increase in early summer methane emissions shows that atmospheric warming has begun to considerably affect the methane flux dynamics of permafrost-affected ecosystems in the Arctic. The authors provide long-term observational evidence of an increasing trend of early summer methane emissions from a permafrost site in the Lena River Delta linked to atmospheric warming. This observed trend constitutes a major development given the thick and cold permafrost in the study area.

While increasing methane emissions from thawing permafrost are anticipated to be a major climate feedback, no observational evidence for such an increase has previously been documented in the literature. Here we report a trend of increasing methane emissions for the early summer months of June and July at a permafrost site in the Lena River Delta, on the basis of the longest set of eddy covariance methane flux data in the Arctic. Along with a strong air temperature rise of 0.3 ± 0.1 °C yr −1 in June, which corresponds to an earlier warming of 11 d, the methane emissions in June and July have increased by roughly 1.9 ± 0.7% yr −1 since 2004. Although the tundra's maximum source strength in August has not yet changed, this increase in early summer methane emissions shows that atmospheric warming has begun to considerably affect the methane flux dynamics of permafrost-affected ecosystems in the Arctic.
Methane is an important greenhouse gas contributing roughly 16% to the radiative forcing from all well-mixed greenhouse gases 1 . Current estimates 2 of the global methane source strength converge at 550-900 Tg yr −1 , with a contribution of 25 ± 14 Tg yr −1 from high-latitude wetlands in the Arctic tundra 3 . These wetlands and other inland waters form a globally relevant methane source and are the most important source of uncertainty in the global methane budget estimate 2 . Reasons for this uncertainty include the shortage of flux observation sites and short time series, the high spatial heterogeneity of permafrost-affected landscapes, the pronounced temporal dynamics of methane emissions and the complex controls on methane production, transport and oxidation [4][5][6] .
With a warming in the Arctic region that is at least twice the global average 7 , methane emissions are projected to increase in response to permafrost degradation 8 . This increase is anticipated to be driven by the remobilization of previously frozen soil organic matter or more specifically, by enhanced decomposition rates, an extension of the thawing season and a deepening of the seasonally thawed active layer into the perennially frozen layer underneath 9 . Increasing methane emissions are likely to lead to an additional rise in temperature, resulting in positive feedback with the potential to accelerate and further exacerbate climate change 10 . The resulting loss of permafrost is irreversible at centennial time scales, particularly if permafrost thaw passes a tipping point 11 . Currently, there is only moderate evidence with low agreement on whether northern permafrost regions are already releasing additional methane due to thaw 8 . Besides sparse indirect indications of increased carbon loss 12,13 , neither atmospheric measurements nor inversion or biospheric models currently provide clear trends in methane emissions from high-latitude wetlands for the recent past 14 . It also remains unclear whether rising temperatures equally affect the two counteracting processes of methane production and oxidation 6 .
Here we analysed a methane flux dataset acquired over 16 yr between 2002 and 2019 in the North Siberian Lena River Delta (72° 22 ' 24.1' N, 126° 29' 49.7' E). This dataset currently forms the longest eddy covariance methane flux record from the Arctic, allowing for the first trend analysis of directly observed methane fluxes in an Arctic ecosystem. We also examined the dynamics in the fluxes, estimated annual and seasonal flux budgets, and identified flux drivers and proxies, which explain the observed fluxes on various time scales.
Article https://doi.org/10.1038/s41558-022-01512-4 transport 22,23 . An earlier development of vascular plants also promotes a decoupling of air and soil temperatures via both stronger radiation shading and higher evapotranspiration rates.
For August, we detected a marginal decline of 0.12 ± 0.15 mmol m −2 yr −1 (0.4 ± 0.5% yr −1 ), with no statistically significant trend. Thus, the maximum source strength of the tundra has not changed. This stable state is supported by equally minor trends in the following two flux drivers: a rise in soil temperature (polygon centre, 30 cm depth) by 0.01 ± 0.02 °C yr −1 , which is similarly non-significant as the thaw depth that declined by 0.07 ± 0.13 cm yr −1 . The thaw depth measurements, however, may be biased since they were conducted by mechanical probing, which can underestimate the degradation of permafrost due to a subsiding surface resulting from ground ice thaw and drainage.
For September, we found a decrease of 0.45 ± 0.16 mmol m −2 yr −1 (2.1 ± 0.7% yr −1 ), which we did not consider a reliable trend since the two statistical tests disagreed on the significance. In addition, data availability was limited in September, leading air and soil temperature trends to change from positive to negative when their estimation was not based on all available temperature data, but only on data from those periods for which methane flux data were available. This change indicates that the selected 7 yr were not representative of the entire measurement period.

Flux budgets and dynamics
The methane fluxes were subject to an annual cycle for which we estimated a mean annual budget with a standard deviation of 171.5 ± 12.3 mmol m −2 yr −1 (Fig. 2). This balance ranges between the lower quartile and the median of 236.9 mmol m −2 yr −1 in a right-skewed distribution, which was obtained from 168 annual budgets at various northern tundra sites 24 . The reasons for the comparatively low fluxes at our site include the areal dominance of elevated polygon rims with non-water-saturated soils, the shallow active layer, the relatively low organic matter content 25 and one of the lowest mean ground temperatures in the northern hemisphere 26 .
Roughly 61% of the mean annual budget can be ascribed to the thawing season, which begins in late May when air temperatures start For June and July, we found statistically significant increases in methane emission of 0.41 ± 0.16 mmol m −2 yr −1 (2.1 ± 0.8% yr −1 ) and 0.44 ± 0.15 mmol m −2 yr −1 (1.7 ± 0.6% yr −1 ), respectively. The numbers in brackets denote the change relative to the monthly mean flux. These changes in early summer are equivalent to an increase of the mean annual methane budget by 0.5 ± 0.1% and are probably attributable to the significant air temperature rise of 0.3 ± 0.1 °C yr −1 in June. This warming apparently had a larger effect on methanogenesis than on methanotrophy 6,15 and corresponds to a significant increase in growing degree days (GDD) from 118 °C to 214 °C in June. Using GDD = 60 °C as a benchmark for a representative cumulative heat input while maintaining high data coverage, the warming advanced towards the frozen season by 10.8 ± 5.2 d during the period from 2002 to 2019. Contrary to air temperature, we did not detect a significant soil temperature increase in the top soil layers (1 cm and 5 cm depth) and only a small and barely significant one at 10 cm depth in June. Similarly, there were no detectable trends in near-surface hydrology. We therefore assume the atmospheric warming to have caused an earlier onset of the vegetation period 16 , which in turn was probably a main driver of increased early summer methane emissions. The following mechanisms are considered to explain the relationship between air temperature (in addition to and independently of soil temperature changes) and methane emission in early summer. Higher air temperatures lead to both earlier snowmelt and onset of the development of the vegetation. Remote sensing also shows a greening in the Lena River Delta, with the strongest trends in the southern part along the main river channels where our site is located 17 . Graminoid vascular plants such as sedges which, together with mosses, dominate the local vegetation, have been described as benefitting most from earlier snowmelt and higher air temperature in spring 18 . Earlier development of vascular plants probably promotes earlier and stronger release of labile organic substances such as root exudates and fresh root litter into the rhizosphere 19 , providing substrates for the formation of methane 20,21 . Its emission is promoted by the earlier foliating sedges that are known to foster upward methane transport through their extensive air spaces (aerenchyma), bypassing aerobic soil horizons and thus lowering the fraction of methane oxidized during   35,36 . While the source strength increased in June and July, the peak emissions during August remained unchanged. These trends are consistent in sign and similar in significance when using non-gap-filled data only.
In June, graminoid plant species begin to grow, thus facilitating upward methane transport through plant-mediated pathways 22 . The increasing fluxes reach their highest rates in early August, when the ground is thawed to a depth of approximately 50 cm. The rate of change in methane emissions during the subsequent drop in September is more pronounced than during the previous rise, leading to an asymmetric course around the flux peak. During the 4-month thawing season, the methane fluxes show a diurnal variation, with nocturnal fluxes being 7% lower on average than daytime fluxes (Fig. 3, top panels). Around early September, the vegetation period incrementally ceases, and in late September, a snow and ice cover develops, causing methane to start accumulating in a sub-surface stock between sporadic releases in the form of burst events (Fig. 4, top panels). With air temperatures already far below freezing, soil temperatures remain near 0 °C during the zero-curtain period, allowing cold-adapted methanogenic archaea to remain active 27,28 towards late November. This refreezing season accounts for about 14% of the mean annual budget. After the soil is completely frozen in December, methane fluxes continue to decrease at a slow pace by gradually depleting the soil methane pool, predominantly through diffusional processes. This frozen season, in which methanogenic microorganisms remain dormant, contributes around 25% to the mean annual budget.

Flux controls on various time scales
The methane fluxes and the explanatory power of associated flux drivers varied strongly on time scales from diurnal to inter-annual. On the daily scale, the variations in methane fluxes were driven by friction velocity, which was also subject to a diurnal cycle (Fig. 3). This variable also influenced fluxes on the multi-day scale by driving methane bursts during stormy events that were associated with both friction velocity peaks and/or air pressure drops (Fig. 4). Friction velocity in particular affects soil-atmosphere coupling through pressure pumping that enables advective trace gas transport rates through soil pores and the snowpack to increase far beyond molecular diffusion 29,30 . This turbulence-induced process also facilitates methane transport across the air-water interface through bubble ebullition, that is, the intermittent release of undissolved methane trapped by adhesive forces, via agitation of plants and water 31 .
On the intra-annual scale, carbon dioxide fluxes determined with the eddy covariance method served as a proxy for methane fluxes between October and May (Fig. 5). During this refreezing and frozen season, the proportionality of methane and carbon dioxide fluxes (Fig. 5,top panel,left) suggests that their ratio in the pore space of soil and snow is rather invariant and that both trace gases were subject to the same physical transport processes. From June to August, thaw depth held explanatory power as an indicator of the maximum soil volume available for methanogenesis. However, this relationship was absent in September, when flux rates decreased despite near-maximum thaw depths. In September and October, surface albedo had a regulative influence, when it formed a proxy for enhancing snow coverage, and thus the decreasing coupling between ground and atmosphere. During the snowmelt in May, however, this coupling was less important due to both the ground being still largely frozen and its methane stock probably being depleted.
On the annual scale, a large part of the variability could be explained by air and soil temperatures as the fundamental variables driving metabolic processes such as methane production (Fig. 5). The soil temperature in a polygon centre at 20 cm depth was the best predictor compared with soil temperatures at other depths and locations within the polygon. This soil temperature reflects the water-saturated conditions in the centre, the timing of the annual thaw/freeze cycle and higher organic matter content in the top soil 27,32 . The annual course of the air temperature preceded that of the methane flux, implicating hysteresis effects 33 . Thus, the air temperature was most predictive, when described by two fits. The resulting apparent temperature sensitivities  (Q 10 ) were 1.9 (March to July) and 1.5 (August to February). This difference in Q 10 values indicates that the methane fluxes appeared more responsive to changes in air temperatures during the early summer than thereafter, when soil temperatures were more influential. Besides air temperatures derived from observations, the annual variation was also well represented by air temperatures based on reanalysis data 34 . Their widespread availability thus offers the opportunity to estimate annual budgets or support regional flux upscaling in remote and sparsely inhabited tundra landscapes.
On the inter-annual scale, the period from 1 July to 15 September, which centres around the mean annual flux peak on 8 August, was taken into consideration for each year. Predictive power was provided by    variables that best describe this maximum development stage of the ecosystem: soil temperature in a polygon centre at 30 cm depth, thaw depth and growing degree days (Fig. 5). The explanatory power of growing degree days was not diminished when observed air temperatures were replaced with air temperatures from reanalysis data 34 .

Conclusion
Earth system models project an increased carbon release from permafrost-affected landscapes due to climate change, albeit there has been little agreement on timing and magnitude due to observational constraints and an incomplete knowledge on drivers and relationships 14 . To our knowledge, we provide the first observational evidence of an increasing trend of early summer methane emissions from tundra wetlands linked to atmospheric warming. The observed trends may still be moderate, yet they constitute a notable development given the very thick and cold continuous permafrost in the study area compared with most other observational sites. We also provide further support for the importance of the cold season, in which we estimate 39% of the total annual methane to be released while continuous data coverage remains a substantial challenge in the Arctic winter.

Online content
Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41558-022-01512-4.

Site description
The Lena River Delta is located within the zone of continuous permafrost on the fringe of the Laptev Sea in northern Siberia. The delta region is situated in a continental Arctic climate with a mean annual air temperature of −12.3 °C and a mean annual precipitation of 169 mm (ref. 37 ). One of the numerous islands is Samoylov Island, which has been hosting researchers from various fields since the late 1980s.
The flux tower was established on the late Holocene river terrace, whose patterned ground has been shaped by frost action processes that created an ice-wedge polygonal tundra with thermokarst lakes and ponds. This landscape is characterized by a mosaic of depressed polygon centres, which often have a moist to wet water regime supporting hydrophytic sedges and mosses, as well as elevated polygon rims, which remain moist to moderately dry and covered by mesophytic sub-shrubs, mosses and lichen 22 .

Instrumental setup
An

Flux processing
Half-hourly fluxes were calculated with the EddyPro 38 software (version 6.2.2). As part of the initial raw data processing, we implemented a double rotation tilt correction to reduce the means of both lateral and vertical wind to zero. Turbulent fluctuations of the wind and gas concentration signals were extracted by trend removal through an exponential running mean, with a time constant of 30 s during 2002 to 2009 and a linear detrending from 2010 onwards. Time lags between wind and gas concentration signals were identified by covariance maximization and compensated for by applying setup-specific plausible time lag windows and defaults. After the actual flux calculation, we subsequently corrected these fluxes using the WPL correction 39 to account for fluctuations in heat and water vapour in the air. We then conducted a correction to account for high-pass filtering 40 due to a finite averaging interval, and low-pass filtering 41 due to both sensor separation and sensor imperfection. When methane fluxes from both open-path and closed-path gas analysers were available, only the fluxes determined from closed-path gas concentration data were used for further processing. The calculated methane fluxes underwent a comprehensive quality assessment routine with various individual steps to ensure that only reliable fluxes were interpreted. We initiated this procedure with the raw data screening, where half-hourly time intervals were discarded when the measured methane concentrations did not lie between the limits 1 and 3 ppm and/or the number of spikes 42 for the three wind vectors and the methane concentration exceeded 5. Subsequently, we identified fluxes with non-steady-state conditions using a stationarity test 43 . These fluxes were removed, as well as the fluxes that did not pass the integral turbulence characteristics test 43 due to an absence of turbulent conditions. To further recognize a deficiently developed turbulence as well as instrument malfunctions, we examined the skewness and kurtosis 44 of flux-specific scalars. Fluxes were dismissed when the skewness was either less than −2 or greater than 2 and/or kurtosis was less than 1 or greater than 8. In addition to these steps, we filtered the open-path fluxes according to the following steps. To ensure that the WPL correction 39 was not distorted by poor-quality data, we rejected methane fluxes if the corresponding sensible and latent heat fluxes had failed the preceding quality assurance. Further steps involved discarding fluxes when the signal strength (that is, the available optical power in the laser path of the methane analyser) was below 10%. Fluxes were also removed when lags determined by the time lag compensation was outside the interval of −2 to +2 s. Finally, we computed the percentiles for the leftover methane fluxes, and fluxes below the 1st percentile and above the 99th percentile were excluded.
The remaining data formed the basis for further analysis and included 70,202 half-hourly methane flux records that were observed during 1,986 d within 16 yr. The thawing, refreezing and frozen seasons account for 63%, 13% and 24% of these records, respectively. To our knowledge, this is the most comprehensive eddy covariance methane flux dataset available from an Artic permafrost-affected ecosystem.

Flux budget estimation
At first, we computed a mean flux and its standard deviation for each calendar day using half-hourly fluxes. These daily means formed a mean annual course that was then approximated with a Fourier series whose integration yielded the mean annual budget. Its uncertainty was obtained by processing the standard deviations with standard error propagation techniques. We used three standard deviations around the mean (μ ± 3σ) and the few missing standard deviations for some winter days were replaced with monthly means. The cumulative flux, F CH4 , was calculated as where i is the index of summation, x depicts the annual time scale as calendar days, and a, b and w symbolize fitting parameters derived through least squares regression. The mean annual course could be well reproduced as shown by a coefficient of determination of 0.93.
The mean annual budget was allocated to three seasons, which we defined on the basis of the soil's thermal state. The thawing season began on 24 May when the daily mean surface temperature (determined via the outgoing long-wave radiation) exceeded the freezing point. When the surface temperature again dropped below zero around 27 September, the following refreezing season started. The subsequent 'frozen' season began when sub-zero temperatures prevailed in the entire soil from 23 November onwards.

Flux driver search
In search of adequate predictors on different time scales, the relationship between observed fluxes and potential flux drivers was analysed by means of simple least squares regression. In doing so, we isolated the respective effect of individual flux drivers (while neglecting the mutual effect of multiple flux drivers) to compare their explanatory power. In many cases, a linear fit was applied. The uncertainty in the functional relationships was assessed with 95% confidence intervals. For characterizing the relationship with soil temperatures on the annual scale, we used an exponential approach, whose curve flattening well depicted the low but steady emissions during the frozen season.

Nature Climate Change
Article https://doi.org/10.1038/s41558-022-01512-4 F CH4 denotes the methane flux, T soil represents the soil temperature and b 1-3 refer to the coefficients to be estimated. For the curve fitting of fluxes with thaw depth as a predictor, we used the same parameterization but excluding the coefficient 'b 3 '. The effect of air temperature on methane emissions was well described by classic Q 10 relationships 45 , which provided additional information on the temperature sensitivity of the emissions.
F base describes the basal emission at the reference temperature T ref which was set to 15 °C, and the scaling factor 46 γ was held constant at 10 °C. Q 10 indicates the temperature sensitivity and specifies a value by which the emission multiplies/divides when the temperature rises/ drops by 10 °C. These two fitting parameters, however, disregard the effect of various other hydrometeorological and biogeochemical factors on methane fluxes, and thus rather state an apparent temperature response. The correlation between flux and environmental variables was assessed with the coefficient of determination.
For the determination of growing degree days, we averaged the daily maximum and minimum air temperatures for each day. This mean was adopted if it was greater than the base temperature of 0 °C; otherwise, it was set to zero. Subsequently, we determined the cumulative sum of all daily values for each year. Using data for growing degree days, soil temperatures and thaw depths in investigating the year-to-year variability depended on the availability of methane fluxes: for each year, we calculated a mean flux that was based on data between 1 July and 15 September. The means of the three predictors were then computed using only data from time intervals when fluxes were available.

Flux modelling
For filling the gaps in the methane flux time series, we implemented a random forest 47,48 , that is, a supervised algorithm which is robust with respect to noise and is able to deal with nonlinear relationships, complex interactions as well as missing data. Its error rates compare favourably with other machine learning techniques. The forest consists of a large collection of individual regression trees that operate as an ensemble. Each tree is based on a bootstrapped dataset of observational features, with duplicate samples being allowed. Growing the roots and splits of each tree occurs through a random subset of features in each bootstrapped dataset. Therefore, each leaf corresponds to the average flux rate from a distinctly different cluster of observations. The randomly varying data input (in terms of both samples and features) in the course of the growing process results in a wide variety of de-correlated trees, which constitutes a favourable way to resolve the bias-variance trade-off in the random forest. Its performance was further enhanced by manually tuning inherent hyperparameters (minimum leaf size, maximum number of splits, number of learning cycles), while its complexity for optimal predictive power was assessed through 10-fold cross validation.
On the basis of the outcome of the previous regression analysis, we constrained the model with daily means of (1) air temperatures and soil temperatures from a polygon centre in 20 cm depth to reflect methane production, (2) air pressure and friction velocity to account for methane transport and (3) surface albedo and thaw depth to provide additional explanatory power during the shoulder seasons and the summer flux peak period when the temporal variability was greater than in the winter. To avoid multicollinearity among the six predictor variables, we carried out two tests 49 . The largest condition index and variance inflation factor were 7.9 and 3.3, respectively, and thus well below 30 and 10, above which these scores are indicative of problematic multicollinearity. The model was able to reproduce the observed fluxes well and yielded, in the absence of heteroscedasticity, a coefficient of determination of 0.75. The mean absolute error was 138.2 µmol m −2 d −1 , which translates to a mean relative error of 18.7%.

Flux trend evaluation
To detect changes in emission rates of individual months over several years, we only used months with a minimum fraction of 30% of observed fluxes and full data coverage after gap-filling. Importantly, using non-gap-filled data for trend analyses resulted in trends consistent in sign and similar in significance. Thus, using gap-filled data compensates for uneven data distribution across the years and allows for robust trend estimation without introducing artefacts. Following these criteria resulted in June, July, August and September being eligible for analysis, for which data from 9, 11, 13 and 7 yr were respectively available.
We assessed upward or downward changes in methane fluxes over time by fitting a trend line following the Theil-Sen approach. The major advantage of this method is its insensitivity to outliers, whose impact would naturally be larger in our 16 yr dataset compared with the 30 yr period that is traditionally used for the characterization of climates. The magnitude of the trend corresponded to the slope of the obtained trend line. We further divided this absolute trend by the monthly flux mean to obtain the relative trend. The associated uncertainty was given using the 95% confidence interval of the trend line. The trend estimation for other environmental variables was conducted in an equivalent manner.
The statistical significance of the trend was examined using two modifications of the two-sided Mann-Kendall test 35,36 that both take into account the effect of serial correlation in time series on the probability of trend detection. The objective of using two tests was to find a trade-off between the power of both tests and their capability to preserve the commonly adopted significance level of 0.05 against the background of limited data availability 50 . Within the given uncertainty range, we considered a trend present or absent when both tests led to the same outcome, whereas opposed test statistics indicate minor confidence in the decision on trend detection.

Data availability
The dataset analysed in this study is available from the corresponding author upon request and at the European Fluxes Database Cluster under site ID RU-Sam.