Recent increases in annual, seasonal, and extreme methane fluxes driven by changes in climate and vegetation in boreal and temperate wetland ecosystems

Climate warming is expected to increase global methane (CH4) emissions from wetland ecosystems. Although in situ eddy covariance (EC) measurements at ecosystem scales can potentially detect CH4 flux changes, most EC systems have only a few years of data collected, so temporal trends in CH4 remain uncertain. Here, we use established drivers to hindcast changes in CH4 fluxes (FCH4) since the early 1980s. We trained a machine learning (ML) model on CH4 flux measurements from 22 [methane‐producing sites] in wetland, upland, and lake sites of the FLUXNET‐CH4 database with at least two full years of measurements across temperate and boreal biomes. The gradient boosting decision tree ML model then hindcasted daily FCH4 over 1981–2018 using meteorological reanalysis data. We found that, mainly driven by rising temperature, half of the sites (n = 11) showed significant increases in annual, seasonal, and extreme FCH4, with increases in FCH4 of ca. 10% or higher found in the fall from 1981–1989 to 2010–2018. The annual trends were driven by increases during summer and fall, particularly at high‐CH4‐emitting fen sites dominated by aerenchymatous plants. We also found that the distribution of days of extremely high FCH4 (defined according to the 95th percentile of the daily FCH4 values over a reference period) have become more frequent during the last four decades and currently account for 10–40% of the total seasonal fluxes. The share of extreme FCH4 days in the total seasonal fluxes was greatest in winter for boreal/taiga sites and in spring for temperate sites, which highlights the increasing importance of the non‐growing seasons in annual budgets. Our results shed light on the effects of climate warming on wetlands, which appears to be extending the CH4 emission seasons and boosting extreme emissions.


| INTRODUC TI ON
Methane (CH 4 ) is the second-most important anthropogenic greenhouse gas after carbon dioxide (CO 2 ; Myhre et al., 2014), with a concerning increasing trend as atmospheric CH 4 concentration has more than doubled since 1800 (Etheridge et al., 1998).
The global emissions approached 600 Tg CH 4 year −1 in 2017 (Saunois et al., 2020), driven by human activities, including agriculture, gas from energy production and use, and waste disposal, but also by natural CH 4 sources and sinks (Saunois, Bousquet, et al., 2016;Saunois, Jackson, et al., 2016;Turner et al., 2019;Zhang et al., 2022).Global CH 4 concentrations plateaued during 1999-2006, possibly due to equilibrium in the hydroxyl CH 4 sink (Schaefer et al., 2016).However, between 2007 and 2013, a rapid rise in the global average CH 4 mole fraction of 5.7 ± 1.2 ppb per year was observed (Nisbet et al., 2016).Wetlands are the largest natural CH 4 sources in the atmosphere (Jackson et al., 2020), yet their annual emissions remain among the most uncertain (Dean et al., 2018), and their interannual changes are even greater (Nisbet et al., 2019;Saunois et al., 2020).
What makes it difficult to detect trends in wetland CH 4 emission is the fact that CH 4 emissions from natural ecosystems vary daily and seasonally.Yet, responding to abiotic factors as key drivers of CH 4 fluxes (Armstrong et al., 2015;Bansal et al., 2016Bansal et al., , 2018)), it is expected that these CH 4 fluxes (FCH 4 ) will increase following soil and air temperatures (Yvon-Durocher et al., 2014).A recent study based on seven Alaskan sites found a positive CH 4 response to the 1970-present temperature and precipitation trends across all sites (Ma et al., 2023).In 2007, atmospheric CH 4 increased in the northern high-latitude regions, which was mainly attributed to anomalously high temperatures in the Arctic (Dlugokencky et al., 2009).Similarly, in northern mid-latitudes (30-50° N), strong increases in CH 4 were registered in 2003, 2014, and 2016, likely resulting from record high temperatures during these years (including an extreme heatwave in Europe in 2003) that enhanced wetland CH 4 emissions (Nisbet et al., 2019).Wetland CH 4 emissions also exhibited abrupt increases in northern Europe and Siberia between 2016 and 2018 (Zhang et al., 2021).At boreal sites (around 50° N), emission peaks in CH 4 were found to lag the global peaks from 2015 and 2017 by 1 year; these emissions partially stemmed from wetlands in summer (Nisbet et al., 2019).Conversely, in a global scale study, the 2007-2013 increases in atmospheric CH 4 were attributed to growth from agriculture and waste (Zhang et al., 2021), while a recent study indicated strong positive wetland CH 4 feedback under current warming and changes in precipitation due to climate change (Zhang et al., 2023).
An important knowledge gap is around the contribution to annual CH 4 emissions of non-growing, or shoulder season, emissions (Bao et al., 2021;Wang et al., 2022;Zona et al., 2016).Non-growing seasons can contribute a significant fraction of annual CH 4 emissions in northern ecosystems (up to 50%; Treat et al., 2018), but their contributions to annual budgets vary significantly depending on local conditions (Delwiche et al., 2021;Treat et al., 2018;Zona et al., 2016).The freeze-up of the soil in the fall can lead to a pulse of methane emissions as soil pore spaces are filled with frozen water.
And large spring emission pulses of CO 2 and CH 4 may occur during rapid warming cycles with rain on snow, resulting in soil or lake surface thaw, soil cracking, and the release of gasses stored overwinter, as, for example, reported for a polygonal tundra in Alaska (Raz-Yaseef et al., 2017).A full understanding of non-growing season CH 4 emissions and their drivers remains lacking and would be greatly enhanced by high-frequency measurements coupled with measurements of soil temperature and other abiotic factors (Arndt et al., 2019) to understand temporal variability and its drivers, and ultimately, CH 4 trends over the past few decades.
Ground-based eddy covariance (EC) towers provide highfrequency measurements of CH 4 emissions from natural sources at ecosystem scales (Delwiche et al., 2021;Knox et al., 2019).Syntheses of EC measurements have provided powerful insights into the main drivers of methane fluxes (FCH 4 ): air temperature (Pugh et al., 2018), soil temperature (Turetsky et al., 2014;Zhu et al., 2021), atmospheric pressure (Crawford et al., 2014;Ueyama, Hirano, & Kominami, 2020;Zhang et al., 2018), and water table depth (Chen et al., 2021;Davidson et al., 2019;Turetsky et al., 2014) by aerenchymatous plants.We also found that the distribution of days of extremely high FCH 4 (defined according to the 95th percentile of the daily FCH 4 values over a reference period) have become more frequent during the last four decades and currently account for 10-40% of the total seasonal fluxes.The share of extreme FCH 4 days in the total seasonal fluxes was greatest in winter for boreal/taiga sites and in spring for temperate sites, which highlights the increasing importance of the non-growing seasons in annual budgets.
Our results shed light on the effects of climate warming on wetlands, which appears to be extending the CH 4 emission seasons and boosting extreme emissions.
climate change, climate feedbacks, extreme fluxes, greenhouse gases, hindcasting, methane fluxes, wetlands Machine learning (ML) models can be used to develop models that can extend site-specific observations into longer time periods to assess changes over time.To this end, we reconstruct FCH 4 using ML at 22 wetland, upland, and lake FLUXNET-CH 4 sites across temperate and boreal biomes.We trained our site-specific gradient boosting decision tree (GBDT) ML model using available CH 4 flux measurements from each site, with atmospheric conditions from climate reanalysis data used as predictors.Our goals were to (1) assess annual and seasonal changes in FCH 4 at each site between 1981 and 2018, (2) investigate drivers of these changes, and (3) evaluate changes in the frequency of extreme FCH 4 days (i.e., days with unusually high FCH 4 , see Methods) over this time period.We hypothesized that high emission sites would show the largest increases over time, driven by warming and wet precipitation during their non-growing seasons, as well as by more frequent days with extreme FCH 4 .

| ME THODS
Our methodology is divided into two steps.First, we used data from the FLUXNET-CH 4 V1.0 database (Delwiche et al., 2021;Knox et al., 2019) to train per-site ML models that estimated daily FCH 4 from atmospheric drivers for 22 FLUXNET sites.Second, we applied the ML models to hindcast daily FCH 4 over 1981-2018 at the 22 sites using the Modern-Era Retrospective analysis for Research and Applications, Version 2 (MERRA-2; GMAO, 2015) climate reanalysis dataset as predictors.

| Training data
The model definition process involved ML training, validation, and prediction.For ML model training, we resampled the hourly FCH 4 from the FLUXNET-CH 4 gap-filled database to daily fluxes as the dependent variable and atmospheric data from the MERRA-2 reanalysis dataset over the period 1981-2018 as independent variables.
For each site, we selected the pixels from MERRA-2 (resolution of 0.625° × 0.5° longitude by latitude, or about 3355 km 2 ) within which the site was located.MERRA-2 reanalysis data were aggregated from 1-hourly to daily resolution.Using daily data, as opposed to monthly or annual, as input to our model allowed us to train the ML model on a wider range of values (including extremes) for the independent variables.Based on conceptual knowledge and previous evaluations of FCH 4 predictors at these sites (Knox et al., 2021;Ueyama et al., 2023), we selected the six following predictors from the MERRA-2 dataset: Mean daily surface air temperature (Tair; °C), daily precipitation (prec; mm/day), longwave (LW) radiation fluxes (lwrd; W/m 2 ), incoming shortwave (SW) radiation fluxes (swrd; W/ m 2 ), wind speed (computed from 2-meter eastward wind (Ugrd) and 2-meter northward wind (Vgrd); m/s), and atmospheric surface pressure (pres; PA).Wind speed (WS) was computed as: In addition to temperature and precipitation selected in previous studies (Delwiche et al., 2021;Knox et al., 2021), we also included LW radiation following Morin et al. (2014) and Chen et al. (2015), since LW radiation approximates near-surface air temperature and is commonly used in land-surface models (Shu et al., 2020;Zhang et al., 2018;Zhang et al., 2023).We also included the 3-and 30day rolling mean of the mean daily surface air temperature, because long-term average temperatures (especially prior to the growing season) can influence the summer FCH 4 values (Pugh et al., 2018).
The 3-and 30-day rolling means of surface temperature are also proxies for soil temperature.Furthermore, we calculated a simple representation of cumulative days with optimal photosynthetic activity using an optimal air temperature threshold that we applied to our diverse sites and vegetation types.These cumulative optimal growing degree days (hereafter, GDD) were estimated by computing cumulative temperature days with a surface temperature above 10°C (based on the global range provided by Huang et al., 2019).For sites where more than 360 days per year were typically above 10°C, we set the GDD to Not a Number (NaN) to remove its effect from the ML model training.Yet, this did not apply to any of the final 22 sites used in the hindcasting analysis (see Section 2.3).Since water table depth data are not available from MERRA-2, we instead use the cumulative 3-and 30-day rolling precipitation sum as an approximation of soil water content anomalies (Koffi et al., 2020).

| Machine learning model
We applied a GBDT algorithm for ML (Friedman, 2002).GBDTs have been widely applied to both classification and regression problems (Lu et al., 2018;Yang et al., 2018).GBDT are based on ensembles of decision trees that fit a base learner to the residuals (the gradient of the loss functional that is minimized) at each iteration (Friedman, 2002).By learning iteratively from each weak learner and combining them, it allows for a strong model.The Gradient Boosting Regressor from the Python package scikit-learn (Pedregosa et al., 2011) was used (Irvin et al., 2021).
The ML model was trained for each site separately, dividing each site's measurements into a training set (80%) and a test set (20%).The test set consisted of the first 20% of daily observations to enable us to test if the algorithm was able to predict temporal trends and reconstruct seasonal patterns (as opposed to randomly sampling the daily data).For model tuning, we used a parameter grid search, which conducted an exhaustive search over the specified hyperparameter ranges at each site.The following hyperparameter values and ranges were considered: (1) number of estimators (sequential trees to be modeled); (2) learning rate (impact of each tree on the final outcome); (3) maximum tree depth (to control overfitting); (4) minimum samples in a leaf (to control overfitting); and (5) maximum features (to consider when searching for the best split).
During training, we used a threefold cross-validation (CV) on randomly sampled and equal-sized folds to further counteract (1)

| Filtering FLUXNET-CH 4 sites
The gap-filled target variable FCH 4 was retrieved from the FLUXNET-CH 4 V1.0 database (Delwiche et al., 2021).For gap filling, we used the artificial neural network (ANN) algorithm described in Knox et al. (2019) and Delwiche et al. (2021).The percentage of gap-filled data per season and per year varied considerably, from 0% to about 56% (Table S1).In particular, in the winter season (DJF), gap-filled data could make up to 56%.Although they have been used in previous studies (e.g., Chang et al., 2023;McNicol et al., 2023;Ouyang et al., 2023), the use of gap-filled data is an additional source of uncertainty in our predictions.
Out of 79 available sites, 22 were selected according to three cri-  The final list of 22 sites includes different ecosystems: Fen (6), Marsh (5), Bog (3), Wet Tundra (2), Upland (2), Rice Paddy (managed; 1), Lake (1), Swamp (1), and Salt Marsh (1).Our sites span temperate as well as boreal/taiga (hereafter, boreal) zones.It should be noted that the footprint of the two upland sites is a combination of uplands and wetlands, which may explain why they are net sources of CH 4 .
To facilitate comparisons, we group sites into 16 temperate (latitude range 35.79°-56.45°and one site in the southern hemisphere at −37.39°) and six boreal sites (latitude range 61.83°-68.62°)throughout the paper.An overview of the sites analyzed in this study can be found in Figure 1a and Table 1.(3001 days).They show a considerable variability of mean annual FCH 4 across sites at temperate sites, ranging from 10 nmol CH 4 m 2 s −1 (at RU-Fy2, US-Srr, and US-PFa) to more than 100 nmol CH 4 m 2 s −1 (at US-Tw1, US-WPT, and DE-the; Figure 1b).Annual FCH 4 averages at boreal sites generally range from about 10 to 40 nmol CH 4 m 2 s −1 (Figure 1b).Seasonal changes differ between sites as well, where some sites have stronger seasonal cycles (e.g., KR-Crk, DE-Zrk) than others (e.g., US-Srr, US-Tw4).

| Validation of ML predictions
Applying the models to predict daily FCH 4 over the test set exhibited reasonable accuracy (  An overview of the validation metrics for daily FCH 4 at each site can be found in Table 2.

| Reconstruction of CH 4 fluxes
We predictions to monthly and annual emission rates for the subsequent analyses.

| Change in annual and seasonal emissions
To assess change over the period , we started by comparing emissions over two periods, 1981-1989 and 2010-2018, to test whether significant changes were detectable using an independent two-sample t-test on the annual average site FCH 4 .In addition, we also analyzed seasonal changes for both periods.Broad categories of seasons were used based on calendar quarters (following Ito et al., 2023 andMiller et al., 2016).We then assessed the relative contributions of the seasonal FCH 4 to the total annual emissions at each site for both periods by calculating the fraction of the total annual fluxes contributed by each season.This enabled us to detect changes in the contribution from each season separately from changes in annual emissions.In addition to comparing differences between the two aforementioned periods, we also computed the linear trend over time on the anomaly in annual It should be noted that ML algorithms tend to smooth natural variability and underestimate extreme events (Steininger et al., 2023).Nevertheless, the presumed effect of systematic and random error sources on our results was further minimized by -The use of relatively large data sets.For example, we used periods with 9 years of data (in particular, we used a reference period 2010-2018), which diminish the influence of random noise.
-The fact that we are not reporting absolute changes but relative changes between two periods.It is expected that the use of relative changes will modulate the effects of comparable biases affecting estimates corresponding to the periods being compared.

| Predictor importance
To evaluate the importance of meteorological predictors, we compared variable importance using Shapley-Additive explanation values (SHAP; Lundberg & Lee, 2017, Lundberg et al., 2018).SHAP is

| FCH 4 hindcasting
According to ML hindcasted fluxes, the average annual FCH 4 was equal or greater at all sites during 2010-2018 than during the 1981-1989 reference period, with a significant increase at half the sites (11 out of 22).Across all sites, the percent increased by a median of 5% (range −1 to +24%) and 8% (range + 6 to +11%), across temperate and boreal sites, respectively (Figure 2).
The greatest increase in FCH 4 was at the temperate fen site DE-Zrk, where during the last three decades FCH 4 values have risen by about 24%, and at the fen sites US-Los and DE-Hte, where FCH 4 has surged by 14% and 13%, respectively (Figure 2; Table S2).
FCH 4 values at the United States (U.S.) marsh sites, US-Tw4 and US-Tw1, on the other hand, showed no significant change.Fen and wet tundra sites as well as the rice and lake sites have experienced higher increases in FCH 4 than marshes (apart from the marsh at US-WPT with an 10% increase from 1981-1989 to 2010-2018; Figure 2; Table S2).
Despite the strong interannual variability (Figure 3), for some sites, strong positive trends of increasing annual FCH 4 are apparent (see red dotted line for sites with a significant trend in Figure 3).At

| Seasonal flux reconstructions
Between the reference and study periods, the fall (September, October, November (SON) months) and summer (June, July, August, (JJA) months) experienced the largest FCH 4 increases, while winter (December, January, February, DJF) and spring (March, April, May (MAM)) had the least (Figure 4).In the fall and summer months, US-NC4, RU-Fy2, and NZ-Kop (Figure 4b).However, our results do suggest that fall CH 4 emissions have become progressively more important (Figure 4d).In the fall (September, October, November (SON)), half of the sites (11) show an increase in FCH 4 values of ca.
We also computed the changes from 1981-1989 to 2010-2018 as the relative contribution of each season to the annual CH 4 emission (Figure S7).Here too, the fall and summer months showed the strongest growth (Figure S7c RU-Fy2 is this time at the low end, as MAM is only contributing about 5% in both periods (1981-1989 and 2010-2018).

| Climatic drivers of FCH 4 changes
Given that the responses of FCH 4 to atmospheric drivers are complex and nonlinear, finding a direct correlation between them is unreasonable (Sturtevant et al., 2016).Nonetheless, changes in climatic drivers can be observed at many sites (see Figure S14).and RU-Che: 0.50°C).In fall (SON), at sites that have undergone the most notable rises in FCH 4 (Table S2), temperature increased between 0.67°C (US-ORv) and 1.68°C (Ru-Che/RU-Ch2).Across all sites, temperatures have risen the most in fall (0.97°C) and winter (0.99°C) relative to spring (0.77°) and summer (0.5°C); while in winter the temperature increase has likely not been sufficient to trigger higher FCH 4 emissions, the strongest relative increases in FCH 4 in fall across sites are likely due to the strong increase in surface temperatures.
We found similar trends for shortwave radiation, except for the RU-Che site, which received about 7 W/m 2 less radiation in With GDD and SW radiation being important predictors, these rises have likely further contributed to the trends predicted by our hindcast analysis.
Changes in precipitation show a more diverse picture: while in summer the six sites with significant changes in FCH 4 all show an increase in precipitation (between 0.04 and 0.64 mm/day), seven sites (US-Srr, US-WPT, US-Tw4, US-ORv, US-Myb, FI-Si2, and SE-Deg) with significant changes in FCH 4 have undergone a decrease in precipitation in fall (SON) of −0.10 to −0.77 mm/day (Figure S14).Yet, among these sites, precipitation was only an important predictor at the marsh sites US-Tw4 (2nd most important predictor; see S8) and US-Myb (4th most important).However, as discussed above, our precipitation and related indices may not have been able to fully represent anoxia controls, or the reduction in precipitation may not have been critical for triggering a decrease in FCH 4 .
To better understand the drivers of FCH 4 , we processed the feature importance for each site (Figure 5; Figures S8/S9).At all sites, the longerterm (3 or even 30-day mean) temperature was a more important driver of FCH 4 than the concurrent daily temperature, while the cumulative growing degree days (GDD) also played a vital role across sites.While soil temperature has been suggested to have higher predictive power than air temperature (Rinne et al., 2018), we found a high correlation between the two variables measured at almost all sites (Figure S10) for surface temperatures above −2°C.Other relevant predictors are shortwave radiation at temperate sites and precipitation (30-day sum) at boreal sites.
Although surface temperature is the main driver of FCH 4 across sites, a more detailed picture of the drivers at each site becomes apparent using the SHAP values (Figures S8/S9

| Ecological drivers of FCH 4 changes
In addition to climate drivers, we also evaluated how site-specific features such as site type and dominant vegetation type influenced annual and seasonal changes in FCH 4 both between the two periods (Figure 6) and for the trend over time (Figure S11).
In both analyses, we found that high-emitting sites such as fens and marshes and sites where aerenchymatous vegetation was the dominant plant type showed the most significant (p < .05)changes over time (Figure S11).Lower-emitting sites, such as wet tundra and bogs, also had significant changes when aerenchymatous vegetation was present (Figure 6; Figure S11).Upland, lakes, swamps, salt marsh, and shrubbed fen showed relatively few significant increases.

| Extreme FCH 4 values
Analyzing the changes in the frequency of extreme FCH 4 days Given their high annual FCH 4 (shown by marker size), this highlights fen vulnerability to climate change.Bogs, marshes, and wet tundra also show significant differences.Conversely, we detected only one significant trend in the uplands.As expected, aerenchymatous vegetation (typically associated with the largest annual fluxes) showed more changes than other vegetation types.
season increased significantly across all sites, on average by about 55% at temperate sites and by 80% at boreal sites.The surges were the highest in the fall (SON) and summer (JJA) months and the lowest in the winter (DJF) and spring (MAM).In the fall and summer months, extreme FCH 4 days more than doubled between the two study periods at 59% and 54% of our sites, respectively (Table S3).
On the other hand, in winter and spring, the number of sites where extreme FCH 4 days occurred increased to only 14% and 18% of our sites, respectively (Table S3).Similar to the results for seasonal fluxes, we found considerable differences between sites.For instance, during the winter season (DJF), the number of days with extreme fluxes decreased at three sites (US-Los, US-NC4, and US-Tw4), while they doubled or even tripled (to 9.8-14.2days) at US-Myb, RU-Ch2, and RU-Che (Table S3).Notably, all sites increase (though only 13 sites substantially) extreme FCH 4 days in fall (SON), which seems to imply that this is the season where across

| DISCUSS ION
We used ML models to reconstruct FCH 4 over 1981-2018 using data from 22 sites spanning a range of ecosystem types (various wetland types, upland, and lake) that include temperate and boreal biomes to evaluate trends and drivers of FCH 4 over the last few decades.From 1981-1989 to 2010-2018, our predicted FCH 4 time series indicates modest median increases of 5% (range −1 to +24%) and 8% (range + 6 to +11%) at temperate and boreal sites, respectively, driven primarily by higher surface temperatures.
Our models also predicted substantial changes in seasonal CH 4 contributions, suggesting a significant increase in CH 4 emissions in the fall season across most sites.Both the annual and seasonal changes were strongest at high-emitting sites such as fens and marshes and at sites with aerenchymatous vegetation.Lastly, the importance of extreme FCH 4 days (>95th percentile of the daily FCH 4 values) increased by 55% at temperate sites and 80% at boreal sites on average over the last four decades.In the absence of multidecadal FCH 4 measurement time series, our ML predictions can help unearth the pathways in which flux seasonality and extreme events may have contributed to annual increases over the last four decades.The daily ML predictions underpinning this analysis can investigate these detailed temporal trends and evolutions that cannot be detected with global scale models with coarser timesteps (Chang et al., 2023;McNicol et al., 2023).Days of extreme FCH 4 may play a more vital role under climate change, and their characterization allowed us to detect changes in contributions to total seasonal fluxes over a period of 40 years.

| Temperature driven increases in FCH 4 over the last four decades
Ten out of our 22 sites had significantly greater and increasing annual FCH 4 from the 1980s to the 2010s.The most important explanatory variable for these increases was surface air temperature as well as 3-and 30-day rolling averages of air temperature (Figure 5).This is in good agreement with previous findings that temperature is an important driver of FCH 4 not just in Arctic wetland systems (Ishizawa et al., 2019;Malhotra & Roulet, 2015) but across all ecosystem scales (Bansal et al., 2023;Christensen et al., 2003;Knox et al., 2019;Yvon-Durocher et al., 2014).Day of year was also an important predictor and may have been a proxy for a variety of seasonal variables such as gross primary productivity or the water table fluctuation, which strongly covary with FCH 4 (Knox et al., 2016) and follow a marked seasonal cycle (Irvin et al., 2021).
In good agreement with previous studies (Delwiche et al., 2021;Knox et al., 2021;McNicol et al., 2023), precipitation, while important, was not as important as temperature-related variables.We used cumulative precipitation metrics as proxies for anaerobic conditions at each site, given the lack of consistent water table data.
However, our precipitation and related indices may not have been able to fully represent anoxia controls.Alternatively, for sites where moisture is not limiting and water table depths are relatively stable, FCH 4 would indeed be expected to predominantly respond to temperature changes (Knox et al., 2021;Moosavi et al., 1996).It may also be possible that the pixel value of precipitation retrieved from the MERRA-2 dataset is not representative of the site conditions.
Previous studies (Bechtold et al., 2019) have shown that the coarse resolution of reanalysis datasets like MERRA-2 may contribute to discrepancies between observed and modeled groundwater table depths (Bechtold et al., 2019).given that non-growing season CH 4 emissions can be a significant fraction of annual emissions (Treat et al., 2018).In the North Siberian Lena River Delta (400 km further north and 3500 km of the Russian sites RU-Che and RU-CH2), Rößger et al. reported slightly different seasonal trends for RU-Che and RU-CH2: while since 2004 they observed a significant increase in FCH 4 of about 1.9% per year in June and July, no conclusions could be drawn for August and September emissions (where we found significant increases in FCH 4 ), although the trend seemed to be slightly negative (though not significant).The differences to our findings could stem from the lower mean annual air temperature in the Lena River Delta (the site has one of the lowest mean ground temperatures in the northern hemisphere) that only leads to notable changes in summer; indeed, atmospheric warming was suggested to be the main drivers of the early summer increase due to an earlier onset of the vegetation period (Rößger et al., 2022).
EC measurements in northern Alaska, on the other hand, are consistent with our findings; during the last two decades, terrestrial fall CH 4 enhancements were associated with the "zero curtain" condition (where soils remain unfrozen around 0°C), which have extended by 2.6 days per year into winter (Arndt et al., 2019).

| Relative contribution of non-growing seasons to annual FCH 4
The contribution of non-growing seasons to the total annual FCH 4 is substantial across our sites, amounting up to 25% in spring and ranging from 15% to 40% in fall (Figure S7).While observations in northern wetlands suggest that non-growing seasons (September-May) contribute about 45% to the total annual FCH 4, model estimates from the Global Carbon Project Model Intercomparison regularly underestimate the contribution during this period by about 40% (Ito et al., 2023).Our estimates (45%-55%; Figure S7) are well aligned with previous observations, emphasizing the importance of extending FCH 4 measurements to periods of the year outside the core growing season to further enhance modeling efforts.

| Higher signal of enhanced methanogenesis in fall
The interannual variability in the contribution of fluxes has also changed considerably between 1981-1989 and 2010-2018, due to changes during the summer and fall months (Figure S7).We expected to see more significant changes in the winter and spring as well, but given that the summer and fall have higher emissions as a baseline, it is not surprising that their changes and contributions are the most significant.Mechanistically, it is also possible that enhanced methanogenesis driven by warming-induced increases in plant productivity had a higher signal in the fall than in the spring.Relative to spring, fall temperatures may still be high enough for microbial activity, and increased summer plant productivity may mean that substrates are not limiting, leading to greater methane flux per unit temperature (Chen et al., 2021).However, these mechanisms require further testing, and our analysis can only suggest that summer and fall FCH 4 may be disproportionately important in a warmer world.

| Site and vegetation type vary in FCH 4 trends
The temporal trends that we observed in FCH 4 can, in part, be explained by site and vegetation type.High-emitting sites with aerenchymatous vegetation tended to have the highest changes over time (Figure 6; Figure S11).Fens stand out as an ecosystem that shows strong change across many sites both seasonally and annually (Figure 6; Figure S11).Fens, in addition to being high-FCH 4 systems with a higher abundance of aerenchymatous plants, are also considered to have lower resilience to climate change given their sensitivity to hydrological change (Wu & Roulet, 2014)

| Increased occurrence of extreme FCH 4 days
We found that high-emitting fen sites also showed higher incidences of extreme FCH 4 flux days (Figure S12).

| Limitations
Our analysis relies on the broad annual and seasonal pattern of  (Koh, 2023;Velthoen et al., 2023).
In terms of predicting seasonality, model predictions were capable of capturing the shape of the seasonal growth decline even if absolute FCH 4 values were incorrect.This would suggest that our conclusions about seasonal proportional changes may be more robust than absolute seasonal changes.Note that our conclusions are based on using 3-month seasons (DJF, MAM, JJA, and SON).While this is a simplified definition of seasons, it made it possible to compare our results across sites and with results from prior relevant studies (e.g., Ito et al., 2023;Miller et al., 2016;Rößger et al., 2022).
Regarding extreme FCH 4 , our GBDT model treats daily fluxes as independent, aside from cumulative and lead predictors, which could misrepresent CH 4 transport and release processes affecting the frequency of episodic extreme fluxes.We acknowledge that the patterns discussed here are based on predictions statistically linked to meteorological drivers, which limits the depth of interpretation of patterns.We further acknowledge that our ML is limited in its ability to extrapolate through time over the six sites with negative R 2 in the test dataset, denoting ML predictions worse than the average, which we nonetheless decided to include in our analysis as it adequately represents interannual changes in FCH 4 .We hope that several of the temporal and seasonal changes described here can be investigated with EC measurements once records are sufficiently long.
Finally, it should be noted that our site-specific models have been developed using the FLUXNET-CH 4 dataset, which implies that they were trained on sites with relatively large amounts of FCH 4 that are not necessarily representative of their respective ecosystems.
Caution should therefore be taken for modeling groups in the interpretation of the results.The models are aimed at providing sitespecific dynamics and characteristics for hindcasting rather than representations of the entire ecosystem.

| CON CLUS ION
We developed a ML model approach for reconstructing CH Laboratory (NNH20ZDA001N).We also acknowledge the PI's listed in Table 1 for providing the data used in this study.
. Unfortunately, most CH 4 open-path sensors on EC towers have become operational only during the last decade, precluding the establishment of annual-todecadal trends over the last decades.
teria: First, to ensure sufficient training data for our ML model, only sites with at least two full years of measurements were considered (730 days), leaving us with 30 sites.The 730 days did not have to be continuous.Second, sites were filtered according to the ML model's predictive performance over a threefold cross-validation score of daily fluxes from the training set.Only 24 sites with a coefficient of determination >0.1 over the cross-validation training dataset were selected, that is, in this step we removed sites CH-Cha, DE-Sfn, US-Bi1, US-HO1, US-Snd, and US-Twt.Third, two additional sites (US-Sne and US-Uaf) with a highly negative value (< −100) on the coefficient of determination of the testing set were not considered in the hindcast.The relatively poor performance of the site-specific models for these eight sites is shown in Figures S3 and S6 .
The measurements for training at the 22 selected FLUXNET-CH 4 V1.0 sites range from about 2 years (778 days) to about 8 years F I G U R E 1 (a) Map of the 22 FLUXNET-CH 4 wetland, upland, and lake sites used in this study.(b) Boxplot of the mean annual methane (CH 4 ) fluxes in taiga/boreal (green) and temperate (purple) sites according to the FLUXNET-CH 4 V1.0 dataset.In each box, the central mark (white stripe) indicates the median, and the edges indicate the 25th and 75th percentiles.The whiskers extend to the maximum and minimum non-outlier data.The base map in (a) is the Terrestrial Ecoregions of the World version 2 from the World Wildlife Fund (WWF; Olson et al., 2001).Map lines delineate study areas and do not necessarily depict accepted national boundaries.TA B L E 1 Site characteristics for methane flux for 22 wetland, upland and lake sites.
), except for the NZ-Kop in the test dataset and some overestimations at sites JP-BBY, US-ORv, and US-Tw4.The model also had some issues with capturing certain peak fluxes (e.g., at DE-Hte, US-Los, and RU-Fy2).
FCH 4 relative to the average annual FCH 4 over the 1981-1989 reference period.2.7 | Extreme daily fluxes2.7.1 | Days of extremely high fluxesBesides monthly fluxes, we also assessed the changes in the number of days of extremely high FCH 4 .In order to be considered "extreme," the daily FCH 4 has to be greater than a threshold.For each season and for each site, the threshold was taken as equal to the 95th percentile of the corresponding daily FCH 4 values available for each quarter (season) over the reference period 1981-1989 (i.e., we created a pool of about 92 days (per season) × 9 years at each site, which resulted in about 828 days per site for each season).From this pool of data, we computed the 95th percentile.On average, the number of days per season of extremely high FCH 4 corresponded to about 4.6 days (5% of ca.92 days) over the reference period 1981-1989.2.7.2 | Contributions of extremely high fluxes to annual fluxesWe additionally computed the 95 th percentile of the daily FCH 4 values for each season and for each site for the 2010-2018 period to identify increases in extreme FCH 4 days.This allowed us to quantify any extreme daily FCH 4 occurrences greater than the identified threshold of 4.6 days per season.Finally, we computed the contribution of these extremely high FCH 4 days to the total seasonal flux in the 2010-2018 period.
based on Shapley values, which take the average of the marginal contributions of each feature across all permutations.As opposed to traditional methods for ML model interpretation, SHAP is consistent over differently parameterized trees and therefore more reliable for assessing feature importance across models.SHAP values determine if the effect of the feature is associated with a higher or lower CH 4 flux value (FiguresS8/S9).Based on the SHAP values, we ranked the first, second, and third most important predictors of the 22 sites to summarize the rankings based on individual models.Lastly, to evaluate the ecological drivers of our observed trends, we also evaluated our response indices (significant p-value from t-test between periods, slope and significant p-value of trend over time, etc.) by site type and dominant vegetation type at each site (extracted from the FLUXNET-CH 4 database).

the
German fen sites DE-Hte and DE-Zrk, the FCH 4 anomaly (compared to the reference period 1981-1990) has been positive at least since the year 2000.Some of the Scandinavian sites (SE-Deg, FI-Lom, and FI-Si2) as well as the bog in Japan (JP-BBY) and the marsh in the U.S. (US-WPT) show a similar pattern of mostly positive anomalies since 2000.

FCH 4
increased significantly (p < .05) between the two study periods at 50% and 20% of our sites, respectively (Figure4d,c).On the other hand, in the winter and spring months, we only saw significant increases at one site each (Figure4a,b).However, it is worth noting that seasonal FCH 4 changes from 1981-1989 to 2010-2018 vary considerably within biomes and between sites.For instance, in boreal spring (MAM), mean FCH 4 values increased by more than 10% at six sites (DE-Dgw, DE-Hte, DE-Zrk, KR-Crk, US-Los, and US-WPT), while over the same period and season, FCH 4 values decreased at ,d).As expected, in the northern hemisphere, the summer months (JJA) made up for the highest share of annual fluxes, contributing 20-70%, except for RU-Fy2 (close to 0%) and the southern hemisphere site in New Zealand, NZ-Kop (about 20%).During fall (SON), boreal and most temperate sites contribute 30-70% to annual fluxes, both at boreal and at temperate sites.Although 10 sites indicate a noticeable increase in FCH 4 values in boreal fall from 1981-1989 to 2010-2018, overall, we could not find a significant increase in the contribution of the SON season to the total annual CH 4 flux (likely given the relatively smaller contribution of SON fluxes to the total annual flux, while FCH 4 also increased in summer).Hence, an increase in the relative contribution of SON was not significant.Also, as expected, the boreal winter season (DJF) contributes little to the annual FCH 4 , between 5-12% at boreal sites and less than 20% at temperate sites.Temperate upland RU-Fy2 is an exception, where DJF FCH 4 contributes about 70% to total annual fluxes (Figure S7a).The other exception is the NZ-Kop site, which is located in the southern hemisphere, and hence DJF corresponds to the summer season.The relative share of fluxes from boreal spring (MAM) ranges from 5% to about 23% both at temperate and boreal sites, where F I G U R E 2 Boxplot of reconstructed annual methane (CH 4 ) flux averaged over the period 1981-1989 and over the period 2010-2018 for 22 wetland, upland, and lake sites.In each box, the central white stripe indicates the median, and the box edges indicate the 25th and 75th percentiles.The whiskers extend to the maximum and minimum non-outlier data.Asterisks indicate degrees of significance in the annual methane flux.
For example, between 1981For example, between  -1989For example, between   and 2010For example, between  -2018, the highest increases in FCH 4 occurred at sites where temperatures have gone up the most (at DE-Hte by 1.34°C; DE-Zrk: 1.51°C; Kr-Crk: 0.83°C; summer (JJA) between 2010-2018 than in 1981-1989 despite the increase in temperature.GDD (see Section 2.1) has increased at F I G U R E 3 Anomalies of the reconstructed mean annual methane (CH 4 ) flux at temperate and boreal sites relative to the 1980-1989 mean.Red trendlines are statistically significant at an alpha of 0.05, while gray trendlines are not significant.three sites in summer (JJA): DE-Hte (5 days/season), DE-Zrk (11), and KR-Crk (8), while the other sites did not show significant changes.In fall (SON), GDD has increased considerably at DE-Hte (6 days), US-WPT (12), US-ORv (10), DE-Dgw (14), and FI-Si2 (8).
).As expected, shortwave (SW) radiation has a positive impact on the SHAP value at temperate sites (DE-Hte, DE-Zrk, NZ-Kop, US-Myb, US-Tw4), meaning it is an important predictor of high FCH 4 values.Although precipitation is listed only a few times under the top three predictors, it is still an important predictor, especially the cumulative monthly precipitation (e.g., DE-Dgw, NZ-Kop, RU-Fy2, US-Los, US-Myb, US-NC4, and US-PFa).At boreal sites, GDD (above 10°C; see methods) is also an essential predictor of FCH 4 .F I G U R E 4 Boxplot of reconstructed methane fluxes (FCH 4 ) averaged per season-year for two periods 1981-1989 and 2010-2018 for 22 wetland, upland, and lake sites.(a) Winter: December-January-February (DJF); (b) Spring: March-April-May (MAM); (c) Summer: June-July-August (JJA); and (d) Fall: September-October-November (SON).Each box summarizes the seasonal-year averaged flux rate per period (n = 9).The central mark (white stripe) indicates the median, and the edges indicate the 25th and 75th percentiles.The whiskers extend to the maximum and minimum non-outlier data.
(i.e., days with extreme methane fluxes, see Section 2.7) provide another possible mechanism by which we can explain the FCH 4 trends over time.Between the reference (1981-1989) and study (2010-2018) periods, the number of extreme FCH 4 days per F I G U R E 5 First, second, and third most important features for machine learning (ML) models at each of the 22 temperate and boreal sites.The 30-day mean temperature is the most important feature at half of the analyzed sites.DOY, day of year; GDD, growing degree days; SW, shortwave.F I G U R E 6 Distribution of significant changes (inferred using p-values of a t-test) between reference (1981-1989) and current periods (2010-2018) shown by site and dominant vegetation class.Different seasons are shown using colors.Note that the p-values are shown on a log scale for ease of visualization.Fens had the most occurrences of significant p-values (<.05; data points to the left of the red line).
biomes and wetland classes the number of days with extreme fluxes (with regards to the particular season) has increased most consistently.The number of extreme FCH 4 days are the highest in summer, winter, and fall at sites with aerenchymatous vegetation and in fall at sites with sphagnum moss (FigureS12).While extreme FCH 4 days are by definition about 4.6 days per season (5% of ca.92 days), the share of these extreme days in the total seasonal flux corresponds to 10-15% of the FCH 4 on average across sites (FigureS13).The share of extreme FCH 4 days in the total seasonal flux is the highest during winter (DJF) at boreal sites and during spring (MAM) at temperate sites, with median values of 15% and 16%, respectively.Furthermore, the share of extreme FCH 4 days in the total seasonal flux has slightly increased from 1981-1989 to 2010-2018, suggesting that extreme FCH 4 days are increasing in magnitude, especially in the non-growing seasons.

4. 2 |
Increasing summer and fall contributions to annual FCH 4Our seasonal analyses suggest that summer and fall months have increased in absolute FCH 4 and in relative importance within annual emissions.Sites that had high increases in the annual FCH 4 between 1981-1989 and 2010-2018 also tended to show significant seasonal increases (e.g., German fen and lake sites, DE-Dgw, DE-Hte, and DE-Zrk).In addition to understanding summer dynamics, the fall contribution to annual emissions is an important knowledge gap For example, fen flux increases since 1981 may have been driven by an increase in extreme FCH 4 flux days, particularly in the summer for sites with aerenchymatous vegetation and in the fall and spring for sites with sphagnum as the dominant vegetation.On the other hand, surprisingly, upland CH 4 trends, while non-significant, tended to show increases in summer extreme FCH 4 days.In the wet tundra, spring days of extreme FCH 4 values emerge as being important.Overall, these extreme FCH 4 days can contribute anywhere from 10% to 40% of the seasonal budgets of sites (FigureS13) and thus should be considered in future estimations of CH 4 budgets.The high contribution of extreme FCH 4 in winter at boreal sites may suggest that a few warm days in winter cause most of the FCH 4 .In a warmer world with higher winter temperatures, more extreme FCH 4 days are to be expected.While winter fluxes still play a minor role in annual fluxes, these extreme fluxes may increase their importance in the future.
4 fluxes over 1981-2018 from ground-based EC measurements at 22 tower sites.We found that half (11) of the sites show an increasing trend of methane fluxes (FCH 4 ), driven mainly by an increasing surface temperature.These trends in FCH 4 vary across ecosystem and vegetation type, with high-CH 4 -emitting sites such as fens with aerenchymatous vegetation showing the strongest signal of change.Our study also suggests that climate change may increase the contribution of summer and fall to annual CH 4 budgets.The contribution of extreme FCH 4 days to the total FCH 4 has hardly been addressed until now, and we found that the number of days with extreme FCH 4 values, which currently account for 10-40% of the total seasonal fluxes, has increased during the last four decades, particularly in fens with aerenchymatous vegetation in summer and in the fall and spring for sites dominated by sphagnum moss.Our approach was made possible by the growing understanding of FCH 4 drivers and their ability to predict FCH 4 .While hindcasting is currently necessary to evaluate multi-decadal changes in site-level FCH 4 , extending observation records throughout the year is key to understanding the impact of climate change on the CH 4 budget going forward.AUTH O R CO NTR I B UTI O N S Sarah Feron: Conceptualization; formal analysis; investigation; methodology; software; validation; visualization; writing -original draft; writing -review and editing.Avni Malhotra: Formal analysis; software; writing -original draft.Sheel Bansal: Writing -original draft; writing -review and editing.Etienne Fluet-Chouinard: Writing -original draft; writing -review and editing.Gavin McNicol: Conceptualization; writing -original draft; writing -review and editing.Sara Knox: Writing -review and editing.Kyle Delwiche: Writing -review and editing.Raul R. Cordero: Conceptualization; writing -original draft; writing -review and editing.Zutao Ouyang: Software; writingreview and editing.Zhen Zhang: Data curation; writing -review and editing.Benjamin Poulter: Writing -review and editing.Robert B. Jackson: Conceptualization; writing -review and editing.ACK N OWLED G M ENTS The support of ANID (ANILLO ACT210046) is gratefully acknowledged.AM and EF-C were supported by COMPASS-FME, a multiinstitutional project supported by the U.S. Department of Energy, Office of Science, Biological, and Environmental Research program as part of the Environmental System Science Program.The Pacific Northwest National Laboratory is operated for DOE by Battelle Memorial Institute under contract DE-AC05-76RL01830.SB was partially funded by the Genomic Science Program of the U.S. Department of Energy, Office of Science, Biological and Environmental Research program, grant DE-SC0023084.Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government.GM was partially supported by a National Aeronautics and Space Administration Carbon Monitoring System subaward to the University of Illinois Chicago via Lawrence Berkeley National Table2), with a median site R 2 of 0.54 (range: −7.32 to 0.9) and a median root mean square error (RMSE) of 34 μmol m-2 s-1 (range: 2.01-122.9).It should be noted that the metrics from the model may slightly change if no random seed is set in the model training process.The model is particularly limited in its ability to extrapolate through time over the six sites, with negative R 2 denoting ML predictions worse than the average.Yet, we preserved these sites in the hindcasting because the model adequately represents interannual changes in FCH 4 , which matter for the interpretation of the trend over 1981-2018.We also tested if the ML model was able to replicate seasonal patterns across the selected 22 sites by comparing the FCH 4 seasonality of the FLUXNET-CH 4 V1.0 measurements against the sea- reconstructed daily FCH 4 values for the selected 22 sites over 1981-2018 using MERRA-2 predictors and the trained ML model.
Note: Please note that the R 2 can be negative, which occurs at certain sites.Refer to https:// sciki t-learn.org/stable/ modul es/ model_ evalu ation.html#r2-score for details.Abbreviations: MAE, mean absolute error; RMSE, root mean square error.TA B L E 2 Machine learning model prediction accuracy over the training and test sets based on 22 boreal and temperate sites.
to novel meteorological conditions that had not yet been measured during the measurement period of EC towers.It may also be possible that the pixel value of precipitation retrieved from the MERRA-2 dataset is not representative of the site conditions, thus contributing to discrepancies between observed and modeled groundwater table depths.The timing and overall mag- ML-predicted daily CH 4 fluxes.Our ML models extended severalfold the length of daily flux observation, and the models were likely exposed