A case study of positive seasonal CH₄ emission hysteresis

As an example of seasonal hysteresis, we examined daily estimates obtained from measurements taken at the Bibai Mire in Northern Japan (JP-BBY) where $$ F_{{\mathrm{CH}}_4}$$ is insensitive to the relatively shallow WTD from 2015 to 2017²⁷ (Fig. 1b, c). Although the seasonality shown in $$ F_{{\mathrm{CH}}_4}$$ appears to follow T_air (Fig. 1a, b), a time-dependent $$ F_{{\mathrm{CH}}_4}$$-T_air relationship varies from earlier to later parts of the frost-free season (Fig. 1d–f). Specifically, plotting daily $$ F_{{\mathrm{CH}}_4}$$ as a function of T_air results in a counterclockwise loop from beginning to end of the frost-free season. Similar hysteretic patterns were found using T_soil (Supplemental Fig. 4) and gap-filled CH₄ emissions⁴⁵ (Supplemental Fig. 5), indicating that the hysteresis is not caused by time lags between T_soil and T_air resulting from heat transfer into the soil⁴⁶, and is not driven by biases caused by missing data. These hysteretic patterns suggest that $$ F_{{\mathrm{CH}}_4}$$ should not be represented as a single static function of T_air.

Seasonal CH₄ emission hysteresis among site-years

Overall, we detect positive seasonal $$ F_{{\mathrm{CH}}_4}$$ hysteresis in most site-years recorded in the FLUXNET-CH₄ database, both in terms of $$ H_A$$ and $$ H_{\mu }$$ (75–77% of site-years; Fig. 2). Consistent hysteresis patterns and magnitudes were found with monthly $$ F_{{\mathrm{CH}}_4}$$ and T_air estimates (72–74%, Supplemental Fig. 6), indicating the observed seasonal $$ F_{{\mathrm{CH}}_4}$$ hysteresis is not sensitive to temporal resolution. The non-zero $$ H_A$$ and $$ H_{\mu }$$ values demonstrate intra-seasonal changes in emergent $$ F_{{\mathrm{CH}}_4}$$–T_air dependencies among wetland and rice paddy sites across the globe, and their negatively skewed distribution indicates that the hysteretic responses are not likely to be random. Ignoring seasonal $$ F_{{\mathrm{CH}}_4}$$ hysteresis leads to overestimated (28 ± 46%) and underestimated (−9 ± 35%) $$ F_{{\mathrm{CH}}_4}$$ predictions earlier and later in the frost-free season across wetland and rice paddy sites, and such prediction bias is overlooked by using seasonally invariant T_air dependence models (−4 ± 7%, Supplemental Fig. 7). For example, $$ F_{{\mathrm{CH}}_4}$$ predictions made by a seasonally invariant emergent $$ F_{{\mathrm{CH}}_4}$$–T_air dependence at JP-BBY (i.e., black lines in Fig. 1d–f) are generally biased high and low in the earlier and later parts of the frost-free season, respectively.

Fig. 2

Predominantly positive seasonal CH₄ emission hysteresis inferred from ecosystem-scale measurements across the globe, i.e., CH₄ emissions are generally higher later in the frost-free season at the same temperature.

The distribution of normalized area of seasonal CH₄ emission hysteresis (H_A; a, b) and mean seasonal CH₄ emission hysteresis (H_μ; c, d) to air temperature among site-years derived from the FLUXNET-CH₄ database. Positive seasonal CH₄ emission hysteresis indicates higher CH₄ emissions later in the frost-free season at the same temperature (e.g., Fig. 1d–f). Red dashed lines represent no hysteresis. The corresponding boxplot of site-year specific H_A (b) and H_μ (d) derived from the FLUXNET-CH₄ database. The red central mark, and the bottom and top edges of the blue box indicate the median, and the 25^th and 75^th percentiles, respectively. The black whiskers extend to the most extreme data points not considered outliers denoted in red plus symbol.

To examine how potential controls are related to the observed seasonal $$ F_{{\mathrm{CH}}_4}$$ hysteresis, we analyzed the distribution pattern of $$ H_A$$ under different site classifications and microclimatic conditions. The majority of site-years show positive seasonal $$ F_{{\mathrm{CH}}_4}$$ hysteresis when $$ H_A$$ values are categorized into (1) different ranges of mean T_air measured in the frost-free season (Supplemental Figs. 8), (2) different wetness conditions indicated by higher and lower mean WTD later in the frost-free season (Supplemental Fig. 9), and (3) different ecosystem types (Supplemental Fig. 10). Intra-seasonal changes in emergent GPP–T_air dependencies show about equal site-year proportions of positive and negative $$ H_A$$ values (48% and 52%, respectively; Supplemental Fig. 11a), suggesting that GPP does not directly contribute to the observed seasonal $$ F_{{\mathrm{CH}}_4}$$ hysteresis. Further, predominantly positive seasonal $$ F_{{\mathrm{CH}}_4}$$ hysteresis is detected using T_soil measured at the shallowest (Supplemental Fig. 12) and deepest (Supplemental Fig. 13) site-specific soil layers, indicating substantial intra-seasonal variability in the $$ F_{{\mathrm{CH}}_4}$$-T_soil relationship. Overall, the wetland and rice paddy observations in the current FLUXNET-CH₄ database suggest that $$ F_{{\mathrm{CH}}_4}$$ are generally higher later (i.e., after reaching maximum seasonal T_air or T_soil) in the frost-free season at a given T_air and T_soil. These hysteretic responses emerged across climate zones with various GPP and frost-free season lengths, and were not directly attributable to intra-seasonal changes in T_air and T_soil (Supplemental Fig. 14).

Divergent temperature responses among sites and years

In terms of the magnitude of seasonal $$ F_{{\mathrm{CH}}_4}$$ hysteresis, intra-seasonal changes in emergent $$ F_{{\mathrm{CH}}_4}$$-T_air dependence vary substantially among site-years within each ecosystem type (Fig. 3), despite being predominantly positive (Fig. 2). For each ecosystem type, the large inter-annual (i.e., different years within the same site) and inter-site (i.e., different site-years within the same ecosystem type) variability highlights the challenge of quantifying a universal and robust emergent $$ F_{{\mathrm{CH}}_4}$$–T_air dependence across wetland and rice paddy sites. For example, using the Boltzmann–Arrhenius function (Methods) to represent the emergent $$ F_{{\mathrm{CH}}_4}$$-T_air dependence of an ecosystem type cannot accurately reflect the site- and time-specific emergent relationships between $$ F_{{\mathrm{CH}}_4}$$ and T_air (Fig. 3). A single static function of T_air thus cannot provide accurate estimates of $$ F_{{\mathrm{CH}}_4}$$, even though meta-analyses using the same functional form suggested that such a representation would lead to consistent emergent $$ F_{{\mathrm{CH}}_4}$$–T_air dependencies among aquatic, wetland, and rice paddy ecosystems³⁰. Considering intra-seasonal variability in emergent $$ F_{{\mathrm{CH}}_4}$$–T_air dependence leads to higher and lower apparent activation energies for $$ F_{{\mathrm{CH}}_4}$$ during earlier and later parts of the frost-free season, respectively (Supplemental Fig. 15a). Our findings indicate that the $$ F_{{\mathrm{CH}}_4}$$ temperature sensitivity is an emergent property that varies substantially with space and time and thus cannot be sufficiently generalized for formulating mechanistic CH₄ models, regardless of its functional form.

Fig. 3

Large differences in intra-seasonal, inter-annual, and inter-site $$ F_{{CH}_4}$$ emergent temperature dependencies are found for all examined ecosystem types.

Thin lines represent the site- and time-specific emergent dependencies of CH₄ emissions on air temperature inferred from daily measurements collected at bog (a), fen (b), marsh (c), peat plateau (d), rice paddy (e), salt marsh (f), swamp (g), and wet tundra (h) sites. Thick black lines represent ecosystem-type specific emergent dependencies of CH₄ emission on air temperature inferred from the Boltzmann–Arrhenius function that do not recognize spatial heterogeneity and temporal variability. The results inferred from earlier and later parts of the frost-free season, and full frost-free season are colored in red, blue, and black, respectively.

Factors other than temperature modulate CH₄ emissions

We applied two approaches to evaluate factors regulating the emergent $$ F_{{\mathrm{CH}}_4}$$-T_air dependence and examine the degree of complexity needed in $$ F_{{\mathrm{CH}}_4}$$ parameterizations in biogeochemical models. In the first approach, we examined the effects of T_air, ecosystem-type variability (i.e., differences between ecosystem types), inter-site variability, inter-annual variability, and intra-seasonal variability on $$ F_{{\mathrm{CH}}_4}$$ predictions. Specifically, $$ F_{{\mathrm{CH}}_4}$$ estimates obtained from six sets of regression models selectively representing the above-mentioned variability (Methods; Supplemental Table 2) were evaluated to investigate how spatial and temporal complexity influences model performance. In the second approach, we trained a random-forest model (Methods) with the FLUXNET-CH₄ database to identify factors controlling the hysteresis parameter $$ a_{hys}$$ (Methods) that quantifies the functional relationship between $$ F_{{\mathrm{CH}}_4}$$ and T_air. To assess whether an observationally inferred model can be constructed for $$ F_{{\mathrm{CH}}_4}$$ estimates, we evaluated the predictive power of a hybrid model that uses the random-forest predicted $$ a_{\mathrm{hys}}$$ to describe the emergent $$ F_{{\mathrm{CH}}_4}$$–T_air dependence (Methods; Eq. 1) in each part of the frost-free season.

The seven T_air dependence models (six regression and one hybrid) can be broadly categorized into three tiers based on the absolute bias relative to the measured $$ F_{{\mathrm{CH}}_4}$$: (1) employing a universal emergent $$ F_{{\mathrm{CH}}_4}$$–T_air dependence inferred from measurements across the globe without representing spatial and temporal variability (76.2% biased); (2) including ecosystem-type variability (i.e., the emergent $$ F_{{\mathrm{CH}}_4}$$–T_air dependence is inferred from measurements collected at the same ecosystem type, so sites within an ecosystem type are uniformly represented; 63.5–63.9% biased); and (3) including ecosystem-site variability (i.e., the emergent $$ F_{{\mathrm{CH}}_4}$$–T_air dependence is inferred from measurements collected at each site; 38.1–45.9% biased) (Fig. 4). Our results suggest that representing ecosystem-type variability does not necessarily improve $$ F_{{\mathrm{CH}}_4}$$ estimates, because the absolute bias of modeled $$ F_{{\mathrm{CH}}_4}$$ is comparable with that estimated by using a universal emergent $$ F_{{\mathrm{CH}}_4}$$–T_air dependence, except for bog, peat plateau, and wet tundra sites (Fig. 4a). For each ecosystem type, the absolute bias of modeled $$ F_{{\mathrm{CH}}_4}$$ is reduced when ecosystem-site variability is represented, demonstrating the need to recognize inter-annual and inter-site variability (e.g., Fig. 3). For each T_air dependence model, the absolute bias of modeled $$ F_{{\mathrm{CH}}_4}$$ is generally higher in rice paddies and salt marshes than in other ecosystem types, suggesting that $$ F_{{\mathrm{CH}}_4}$$ in these systems are sensitive to factors other than T_air. For example, timing of irrigation, drainage, planting, and harvesting can all affect $$ F_{{\mathrm{CH}}_4}$$ dynamics in rice paddies⁴⁷.

Fig. 4

The accuracy of CH₄ emission estimates improves with better representation of the large wetland-site variability caused by varying environmental conditions.

The absolute bias relative to measured CH₄ emissions estimated by each model class for each ecosystem type (a). Blue and red bars denote the number of sites and quality-controlled daily data points within each ecosystem type, respectively (b). The abbreviations used in each model group represent air temperature (T), ecosystem-type variability (type), intra-seasonal variability (ISV), hybrid model based on random-forest estimated hysteresis parameter (hybrid), inter-site variability (site), and inter-annual variability (IAV).

Results derived from our random-forest model confirm the importance of ecosystem-site variability in regulating $$ a_{hys}$$ and thereby $$ F_{{\mathrm{CH}}_4}$$ predicted by the hybrid model in each part of the frost-free season (Supplemental Fig. 16). Our random-forest predictor importance analysis indicates that site-year specific $$ F_{{\mathrm{CH}}_4}$$ and T_air values are more important for $$ a_{hys}$$ estimates than other predictors such as latitude, GPP, and ecosystem type. The weak relationships found between seasonal $$ F_{{\mathrm{CH}}_4}$$ hysteresis and latitude (Supplemental Fig. 14h) and GPP (Supplemental Fig. 14d) are consistent with the relatively low predictor importance for $$ a_{hys}$$ found in our random-forest model. Collectively, our results demonstrate the importance of recognizing inter-site, inter-annual, and intra-seasonal variability for the interpretation of emergent $$ F_{{\mathrm{CH}}_4}$$–T_air dependence inferred from measurements across distinct site-years.

When using a universal emergent $$ F_{{\mathrm{CH}}_4}$$-T_air dependence that only represents a generic T_air sensitivity of $$ F_{{\mathrm{CH}}_4}$$ (i.e., the top row in Fig. 4a), the resulting $$ F_{{\mathrm{CH}}_4}$$ predictions substantially underestimate the range of $$ F_{{\mathrm{CH}}_4}$$ measured across wetland and rice paddy sites (Fig. 5a). This generic T_air sensitivity of $$ F_{{\mathrm{CH}}_4}$$ flattens the high temporal and spatial variability^10,11 that strongly controls the timing and magnitude of $$ F_{{\mathrm{CH}}_4}$$, reinforcing the need to parameterize factors other than T_air in CH₄ models. Including factors other than a generic T_air sensitivity of $$ F_{{\mathrm{CH}}_4}$$ (i.e., the bottom row in Fig. 4a) improves $$ F_{{\mathrm{CH}}_4}$$ predictions (Fig. 5b, c), which suggests that $$ F_{{\mathrm{CH}}_4}$$ and emergent $$ F_{{\mathrm{CH}}_4}$$–T_air dependence strongly depend on site- and time-specific environmental conditions. Therefore, models should mechanistically represent CH₄ biogeochemistry, because site- and time-specific emergent $$ F_{{\mathrm{CH}}_4}$$-T_air dependence cannot be accurately parameterized everywhere and all the time. Although many CH₄ models parameterize methanogenesis, methanotrophy, and CH₄ transport for $$ F_{{\mathrm{CH}}_4}$$ modeling¹³, only three of 40 recently reviewed CH₄ models mechanistically represent CH₄ biogeochemistry based on explicit microbial dynamics²⁹. Consequently, implementing process-based representations of CH₄ biogeochemistry in CH₄ models is necessary to improve $$ F_{{\mathrm{CH}}_4}$$ predictions across ecosystem and global scales. Such efforts are imperative because the $$ F_{{\mathrm{CH}}_4}$$ prediction error can increase substantially with increased $$ F_{{\mathrm{CH}}_4}$$, especially for the relatively simple parameterization that only represents a generic T_air sensitivity of $$ F_{{\mathrm{CH}}_4}$$ (Fig. 5c).

Fig. 5

CH₄ emission prediction error increases substantially as measured CH₄ emission increases.

The performance of CH₄ emissions modeled by the regression models that only include a universal emergent CH₄ emission temperature dependence (a), and those that include site- and time-specific conditions (b). The root mean square errors associated with the regression models used in (a) and (b) (bars, left axis) and number of data points (green line, right axis) for measured CH₄ emission bins (c). Two of the 27,130 daily observations have CH₄ emission above 1600 mg C m⁻² d⁻¹, which are not shown for the ease of representation. Lighter colors in the density scatter plot represent denser data points. Solid blue and dashed black lines represent the linear best-fit and one-to-one lines, respectively. The abbreviations used in each model group represent air temperature (T), intra-seasonal variability (ISV), inter-site variability (site), and inter-annual variability (IAV).

Limitations and implications

Additional measurements and analysis of factors controlling methanogenesis, methanotrophy, and CH₄ transport will be needed to investigate the cause of the predominantly positive seasonal $$ F_{{\mathrm{CH}}_4}$$ hysteresis we observed across wetland and rice paddy sites. When anoxic conditions are prevalent and T_soil is the most important driver regulating $$ F_{{\mathrm{CH}}_4}$$^26,27 (e.g., Supplemental Fig. 4), the observed positive seasonal $$ F_{{\mathrm{CH}}_4}$$ hysteresis is consistent with the higher $$ F_{{\mathrm{CH}}_4}$$ driven by higher substrate availability later in the frost-free season²⁵. We identified some environmental drivers affecting the emergent $$ F_{{\mathrm{CH}}_4}$$–T_air dependence at sites where the necessary measurements were available: (1) When WTD drops below the critical zone of CH₄ production later in the frost-free season³¹, the reduced  FCH4 may drive negative seasonal $$ F_{{\mathrm{CH}}_4}$$ hysteresis in a given site-year (e.g., the Kopuatai bog in New Zealand (NZ-Kop), Supplemental Fig. 3). (2) $$ F_{{\mathrm{CH}}_4}$$ may become more sensitive to T_air changes under higher salinity⁴⁸, and our results indicate that seasonal $$ F_{{\mathrm{CH}}_4}$$ hysteresis shifts from positive to negative with increased salinity (e.g., the Sacramento-San Joaquin Delta of California in USA (US-Myb), Supplemental Fig. 17).

As for the emergent $$ F_{{\mathrm{CH}}_4}$$–T_soil dependence, our results suggest that the functional relationship between $$ F_{{\mathrm{CH}}_4}$$ and T_soil may vary non-monotonically along the soil profile. For example, the positive seasonal $$ F_{{\mathrm{CH}}_4}$$ hysteresis inferred from T_soil measured at 16 cm depth is stronger than those at 8 and 32 cm depths at US-Myb (Supplemental Fig. 18). Such a non-monotonic relationship indicates that the magnitude of seasonal $$ F_{{\mathrm{CH}}_4}$$ hysteresis is not simply caused by time lags between T_soil and T_air, suggesting that factors other than temperature can strongly control $$ F_{{\mathrm{CH}}_4}$$. T_soil measured at depths where methanogenesis is occurring will be needed to rigorously examine the emergent dependence of $$ F_{{\mathrm{CH}}_4}$$ on T_soil across the globe, but such depth-dependent measurements are not yet available among sites in the FLUXNET-CH₄ database. To improve understanding of mechanisms leading to seasonal $$ F_{{\mathrm{CH}}_4}$$ hysteresis, we urge further long-term measurements on factors modulating CH₄ biogeochemistry (e.g., WTD, T_soil, microbial activity, and substrate availability), especially in the tropics and the Southern Hemisphere, both of which are sparsely represented in the FLUXNET-CH₄ database. Although seasonal $$ F_{{\mathrm{CH}}_4}$$ hysteresis occurs across seasonal climate and latitudinal gradients (Supplemental Fig. 14), better-representing ecosystems south of 30 °N could affect the partitioning of negative and positive seasonal $$ F_{{\mathrm{CH}}_4}$$ hysteresis inferred from existing measurements. While our synthesis in tropical and subtropical regions shows intra-seasonal changes in emergent $$ F_{{\mathrm{CH}}_4}$$–T_air dependence (Supplemental Fig. 19), future studies are needed to examine seasonal $$ F_{{\mathrm{CH}}_4}$$ hysteresis in wetlands south of 30 °N (that account for about 75% of global wetland $$ F_{{\mathrm{CH}}_4}$$⁶).

The observed seasonal $$ F_{{\mathrm{CH}}_4}$$ hysteresis provides a benchmark to evaluate modeled $$ F_{{\mathrm{CH}}_4}$$ functional responses and should inform and motivate CH₄ model development and refinement. Studies have shown that temporal variations in $$ F_{{\mathrm{CH}}_4}$$ are strongly modulated by substrate and microbial dynamics^33,49,50, which may explain the substantial seasonal $$ F_{{\mathrm{CH}}_4}$$ hysteresis identified in our wetland and rice paddy sites. For example, a model that explicitly represents substrate and microbial dynamics reproduced the observed hysteretic $$ F_{{\mathrm{CH}}_4}$$ to temperature relationships in several wetlands with different vegetation and hydrological conditions³³. Such dynamics could be parameterized in the terrestrial components of Earth system models⁴⁹. Our synthesis thus provides observational evidence for incorporating substrate and microbial dynamics into next generation CH₄ models.

Using the largest available database of ecosystem-scale CH₄ emissions measured by eddy covariance flux towers, we show that the apparent relationships between CH₄ emissions and air and soil temperatures are hysteretic and vary strongly with sampling location and measurement period. Approximately 77% of site-years recorded in the wetland and rice paddy subset of the FLUXNET-CH₄ database⁴⁰ show that CH₄ emissions become higher later in the frost-free season at the same air temperature. This predominantly positive seasonal CH₄ emission hysteresis may be driven by substrate-mediated higher CH₄ production²⁵ later in the frost-free season³³. Changes in environmental conditions also modulate seasonal CH₄ emission hysteresis and thus ecosystem-scale CH₄ emissions.

Our results demonstrate that the relationship between CH₄ emissions and temperature is an emergent property that varies substantially across space and time. A direct integration of measurements across the globe (e.g., inferring a generic temperature sensitivity of CH₄ emissions) may not improve CH₄ model parameterization because such an approach oversimplifies factors controlling CH₄ emissions. Therefore, meta-analyses of CH₄ biogeochemistry should recognize the large intra-seasonal, inter-annual, and inter-site variability of biotic and abiotic conditions that regulate ecosystem-scale CH₄ emissions. Collectively, our analyses highlight the importance of observing and modeling spatial heterogeneity and temporal variability for the modeling of CH₄ biogeochemistry. Since most existing CH₄ models are developed using empirically based CH₄ production or emission temperature dependencies²⁹, our study motivates models to mechanistically represent methanogenesis, methanotrophy, and CH₄ transport to refine estimates of global CH₄ emissions and climate feedbacks⁵¹.

FLUXNET-CH₄ database

The FLUXNET-CH₄ initiative is led by the Global Carbon Project (https://www.globalcarbonproject.org) in coordination with regional flux networks (in particular AmeriFlux and the European Fluxes Database) to compile a global CH₄ flux database of eddy covariance and supporting measurements encompassing freshwater, coastal, natural and managed wetlands, and uplands⁴⁰. Database descriptions, including existing sites, data standardization, gap-filling, and partitioning, have been detailed previously in Knox et al.⁴⁰. We used daily mean temperature (air and soil), gross primary productivity as partitioned from net CO₂ exchange measurements, precipitation, WTD, wind speed, atmospheric pressure, and CH₄ emissions compiled at the 48 wetland and rice paddy sites (Supplemental Table 1) currently recorded in the FLUXNET-CH₄ database. Soil temperature is often measured at different depths among different sites, and only about half of the wetland sites report WTD in the current FLUXNET-CH₄ database⁴⁰. We analyzed the soil temperature reported at the shallowest and deepest measured soil layers at each site to investigate their effects on regulating CH₄ emissions. The wetland and rice paddy data (207 site-years with 62,384 site-days as of this publication) were categorized into eight CH₄ emitting ecosystem types: bog, fen, marsh, peat plateau, rice paddy, salt marsh, swamp, and wet tundra, based on previous classification^52,53. While gap-filled data are examined, they are not included in our discussion to eliminate potential biases caused by the gap-filling procedure⁴⁵.

Frost-free season

We define the frost-free season as the period when the observed temperature (air or soil) is >0 °C to investigate the emergent temperature responses to CH₄ emissions ($$ F_{{\mathrm{CH}}_4}$$) during the biologically active season across distinct climatic zones. Other data sampling thresholds, such as above-zero GPP and above 5% of annual GPP maximum, were examined, and positive seasonal $$ F_{{\mathrm{CH}}_4}$$ hysteresis is identified in 68–81% of site-years (Supplemental Figs. 20, 21), consistent with those inferred from frost-free season. We chose to present the frost-free season results because substantial GPP (e.g., above 5% of annual GPP maximum) is detected when air temperature is well below 0 °C (Supplemental Fig. 22) that may complicate our discussion of varying $$ F_{{\mathrm{CH}}_4}$$ led by temperature changes.

Emergent temperature dependence calculation and the hysteresis parameter $$ a_{\mathrm{hys}}$$

Emergent dependence of CH₄ emission ($$ F_{{\mathrm{CH}}_4}$$) on temperature (air or soil) is determined by fitting frost-free-season daily measurements of $$ F_{{\mathrm{CH}}_4}$$ and air and soil temperatures with a quadratic equation (Eq. 1), the Boltzmann–Arrhenius equation (Eq. 2), and first, second, third, and fifth order polynomials. Daily $$ F_{{\mathrm{CH}}_4}$$ estimates made by site- and time-specific emergent $$ F_{{\mathrm{CH}}_4}$$ temperature (air or soil) dependence models based on the above-mentioned functional forms show comparable root mean square errors (Supplemental Fig. 2). Results inferred from the quadratic equation (Eq. 1) are selected because (1) its functional form is mathematically consistent with the second-order polynomial equation of temperature for methanogenesis inferred from the MacroMolecular Rate Theory^54,55; and (2) it can prescribe seasonal $$ F_{{\mathrm{CH}}_4}$$ hysteresis with a single site- and time- specific parameter ($$ a_{\mathrm{hys}}$$, defined below).

The fits based on the quadratic equation were forced to pass through the origin (assuming zero $$ F_{{\mathrm{CH}}_4}$$ at 0 °C, discussed below) and $$ F_{{\mathrm{CH}}_4}$$ measured at maximum seasonal temperature in each site-year using the Matlab (MathWorks Inc., 2019, version 9.7.0) polyfix function (downloaded from https://www.mathworks.com/matlabcentral/fileexchange/54207-polyfix-x-y-n-xfix-yfix-xder-dydx). The resulting emergent dependence of $$ F_{{\mathrm{CH}}_4}$$ on temperature at any given time period can thus be represented as:

The symbols used in Eq. 1 denote CH₄ emission ($$ F_{{\mathrm{CH}}_4}\left(T\right)$$, mg C m⁻² d⁻¹), hysteresis parameter ($$ a_{\mathrm{hys}}$$, mg C m⁻² d⁻¹ °C⁻²), daily mean temperature ($$ T$$, °C; air or soil), maximum seasonal temperature ($$ T_{\max }$$, °C), and CH₄ emission measured at maximum seasonal temperature ($$ F_{{\mathrm{CH}}_4,T_{\max }}$$, mg C m⁻² d⁻¹). Therefore, the functional relationship between and temperature, described by a quadratic equation (Eq. 1), is only determined by the value of hysteresis parameter ($$ a_{\mathrm{hys}}$$) and site-year variables ($$ F_{{\mathrm{CH}}_4,T_{\max }}\mathrm{and}{T}_{\max }$$).

The two constraints (passing through the origin and $$ F_{{\mathrm{CH}}_4}$$ measured at maximum seasonal temperature) imposed in Eq. 1 are intended to force the two (earlier and later part of the frost-free season) emergent $$ F_{{\mathrm{CH}}_4}$$ temperature (air or soil) dependencies to form a closed apparent hysteresis loop for each frost-free season. By doing so, seasonal $$ F_{{\mathrm{CH}}_4}$$ hysteresis can be quantified as the normalized area enclosed by the two fits, and intra-seasonal changes can be consistently compared among site-years across distinct climate zones. Ignoring $$ F_{{\mathrm{CH}}_4}$$ around 0 °C has small effects on the magnitude and distribution of seasonal $$ F_{{\mathrm{CH}}_4}$$ hysteresis inferred from the current FLUXNET-CH₄ database, although substantial $$ F_{{\mathrm{CH}}_4}$$ may continue when air temperature is around or below 0 °C^32,42. To quantify the effect of ignoring $$ F_{{\mathrm{CH}}_4}$$ around 0 °C, we replaced the constraint of zero $$ F_{{\mathrm{CH}}_4}$$ at 0 °C by the mean $$ F_{{\mathrm{CH}}_4}$$ measured between −0.5 and 0.5 °C at 0 °C for each site-year, and found that the resulting patterns of seasonal $$ F_{{\mathrm{CH}}_4}$$ hysteresis (Supplemental Fig. 23) are consistent with those assuming zero $$ F_{{\mathrm{CH}}_4}$$ at 0 °C (Fig. 2).

Seasonal CH₄ emission hysteresis

We apply a quadratic equation (Eq. 1) to calculate the emergent dependence of CH₄ emission ($$ F_{{\mathrm{CH}}_4}$$) on temperature at the earlier ($$ F_{{\mathrm{CH}}_4,\mathrm{earlier}}\left(T\right)$$) and later ($$ F_{{\mathrm{CH}}_4,\mathrm{later}}\left(T\right)$$) part of the frost-free season separated by maximum seasonal temperature ($$ T_{\max }$$). Two metrics are used to quantify the observed seasonal $$ F_{{\mathrm{CH}}_4}$$ hysteresis: (1) Normalized area of seasonal $$ F_{{\mathrm{CH}}_4}$$ hysteresis ($$ H_A$$), defined as the area enclosed by emergent dependencies of $$ F_{{\mathrm{CH}}_4}$$ on temperature inferred from earlier and later parts of the frost-free season (i.e., $$ H_A=\frac{{\int }_0^{T_{\max }}\left(F_{{\mathrm{CH}}_4,\mathrm{later}}\left(T\right)-F_{{\mathrm{CH}}_4,\mathrm{earlier}}\left(T\right)\right)dT}{\max \left(\mathrm{abs}\left(F_{{\mathrm{CH}}_4,\mathrm{earlier}}\left(T\right),F_{{\mathrm{CH}}_4,\mathrm{later}}\left(T\right)\right)\right)\cdot T_{\max }}$$); and (2) mean seasonal $$ F_{{\mathrm{CH}}_4}$$ hysteresis ($$ H_{\mu }$$), defined as the difference between mean daily $$ F_{{\mathrm{CH}}_4}$$ inferred from measurements taken between later and earlier parts of the frost-free season. In each site-year, positive seasonal $$ F_{{\mathrm{CH}}_4}$$ hysteresis occurs when higher $$ F_{{\mathrm{CH}}_4}$$ are measured later in the frost-free season at a given air or soil temperature. Hysteretic patterns are similar when using either air temperatures (Fig. 1) or soil temperatures (Supplemental Fig. 4), and with either gap-filled (Supplemental Fig. 5) or non-gap-filled (Fig. 1) $$ F_{{\mathrm{CH}}_4}$$⁴⁵. Results derived from air temperature (Fig. 2), soil temperature measured at the shallowest soil layer (Supplemental Fig. 12), and soil temperature measured at the deepest soil layer (Supplemental Fig. 13) all indicate predominantly positive seasonal $$ F_{{\mathrm{CH}}_4}$$ hysteresis across the wetland and rice paddy sites. We chose to present results derived from air temperature for its longer and more continuous record in the wetland and rice paddy subset of FLUXNET-CH₄ database, although soil temperature has been shown to be a better predictor for $$ F_{{\mathrm{CH}}_4}$$^33,42. Specifically, there are 207, 112, and 97 site-years of measurements of air temperature, soil temperature measured at the sallowest soil layer (0–18.3 cm), and soil temperature measured at the sallowest soil layer (32–50 cm), respectively.

Temperature dependence model groups

The measurements extracted from the FLUXNET-CH₄ database were analyzed by seven air temperature (T_air) dependence model groups (six regression models and a hybrid model) to evaluate factors modulating CH₄ emission predictions. We design the six regression models to selectively represent the effects of ecosystem-site variability and ecosystem-type variability on CH₄ emission prediction by labeling data points into different groups. The relationship between CH₄ emission and T_air is analyzed at each part of the frost-free season, each site-year, each site, and each ecosystem type to quantify intra-seasonal, inter-annual, inter-site, and ecosystem-type variability, respectively (Supplemental Table 2). For the hybrid model, we use the hysteresis parameter predicted by our random-forest model to inform the quadratic equation (Eq. 1) for CH₄ emission estimates. The performance of each T_air dependence model group was evaluated to determine the most important model components required for accurate CH₄ emission estimates.

Random-forest model selection

We used random-forest model selection to identify the most important predictors of the hysteresis parameter $$ a_{\mathrm{hys}}$$ (Eq. 1) that determines the functional form of emergent CH₄ emission air temperature dependence and thereby wetland CH₄ emissions ($$ F_{{\mathrm{CH}}_4}$$). Instead of $$ F_{{\mathrm{CH}}_4}$$, the hysteresis parameter $$ a_{\mathrm{hys}}$$ was analyzed, so the results can provide useful information on the source of observed $$ F_{{\mathrm{CH}}_4}$$ hysteresis with an understandable functional form (Eq. 1). Moreover, the most important predictors identified by the machine-learning approach can be compared with the results derived from the other approach using a range of temperature dependence model groups (Supplemental Table 2).

Ten potential predictors were selected for their relatively high predictor importance to $$ a_{\mathrm{hys}}$$: seasonal branch (i.e., earlier or later part in the frost-free season), GPP cumulated in a seasonal branch, precipitation cumulated in a seasonal branch, maximum seasonal temperature, mean temperature in a seasonal branch, ecosystem type, latitude, site, site-year, and $$ F_{{\mathrm{CH}}_4}$$ measured at maximum seasonal temperature. Other potential predictors, including observational year, mean WTD in a seasonal branch, mean wind speed in a seasonal branch, and mean atmospheric pressure in a seasonal branch were examined and showed limited predictive power on $$ a_{\mathrm{hys}}$$. Four potential predictors (seasonal branch, ecosystem type, site, and site-year) were labeled as categorical data and the rest were labeled as numerical data in our random-forest model. The random-forest model selection was performed by the Statistics and Machine-Learning Toolbox in Matlab (MathWorks Inc., 2019, version 9.7.0).

Apparent activation energy for CH₄ emissions

We quantify the apparent activation energy for CH₄ emissions by fitting frost-free-season daily measurements of CH₄ emission and air temperature with the Boltzmann–Arrhenius equation of the form: