Modeled Microbial Dynamics Explain the Apparent Temperature Sensitivity of Wetland Methane Emissions

Methane emissions from natural wetlands tend to increase with temperature and therefore may lead to a positive feedback under future climate change. However, their temperature response includes confounding factors and appears to differ on different time scales. Observed methane emissions depend strongly on temperature on a seasonal basis, but if the annual mean emissions are compared between sites, there is only a small temperature effect. We hypothesize that microbial dynamics are a major driver of the seasonal cycle and that they can explain this apparent discrepancy. We introduce a relatively simple model of methanogenic growth and dormancy into a wetland methane scheme that is used in an Earth system model. We show that this addition is sufficient to reproduce the observed seasonal dynamics of methane emissions in fully saturated wetland sites, at the same time as reproducing the annual mean emissions. We find that a more complex scheme used in recent Earth system models does not add predictive power. The sites used span a range of climatic conditions, with the majority in high latitudes. The difference in apparent temperature sensitivity seasonally versus spatially cannot be recreated by the non‐microbial schemes tested. We therefore conclude that microbial dynamics are a strong candidate to be driving the seasonal cycle of wetland methane emissions. We quantify longer‐term temperature sensitivity using this scheme and show that it gives approximately a 12% increase in emissions per degree of warming globally. This is in addition to any hydrological changes, which could also impact future methane emissions.

Methane emissions from natural wetlands tend to increase with temperature and therefore may lead to a positive feedback under future climate change. However, their temperature response includes confounding factors and appears to differ on different time scales. Observed methane emissions depend strongly on temperature on a seasonal basis, but if the annual mean emissions are compared between sites, there is only a small temperature effect. We hypothesize that microbial dynamics are a major driver of the seasonal cycle and that they can explain this apparent discrepancy. We introduce a relatively simple model of methanogenic growth and dormancy into a wetland methane scheme that is used in an Earth system model. We show that this addition is sufficient to reproduce the observed seasonal dynamics of methane emissions in fully saturated wetland sites, at the same time as reproducing the annual mean emissions. We find that a more complex scheme used in recent Earth system models does not add predictive power. The sites used span a range of climatic conditions, with the majority in high latitudes. The difference in apparent temperature sensitivity seasonally versus spatially cannot be recreated by the non-microbial schemes tested. We therefore conclude that microbial dynamics are a strong candidate to be driving the seasonal cycle of wetland methane emissions. We quantify longer-term temperature sensitivity using this scheme and show that it gives approximately a 12% increase in emissions per degree of warming globally. This is in addition to any hydrological changes, which could also impact future methane emissions.

Introduction
Methane is the second most important anthropogenic greenhouse gas after carbon dioxide (Stocker et al., 2013). It is emitted from natural wetlands as well as from some agricultural and industrial activities (Saunois et al., 2020), with natural wetlands representing ∼33% of annual global methane emissions from all sources for the year 2017 (194 vs. 596 Tg CH 4 yr −1 ). Methane emissions from natural wetlands tend to increase with temperature and are predicted to have a positive feedback under climate change (Gedney et al., 2019). Inclusion of this climate feedback from natural wetland methane emissions reduces the anthropogenic carbon budget to achieve the Paris climate agreement by up to 10% (Comyn-Platt et al., 2018).
Methane emissions from wetlands are difficult to quantify reliably using large-scale atmospheric data, since methane is advected, mixes with methane from other sources, and decays in the atmosphere (Saunois et al., 2020). Therefore, ground-based measurements at individual sites, including eddy covariance and plot-scale data, currently give us our best insights into wetland methane emissions. These studies show that methane emissions often depend strongly on both soil temperature and water table depth (Knox et al., 2019;Roulet et al., 1992;Turetsky et al., 2014;Yvon-Durocher et al., 2014).
Due to the potential feedback to the global climate, modelers make efforts to include wetland methane emissions in land surface schemes that form part of Earth system models (ESMs). These emission models vary in complexity from simple (e.g., the Joint UK Land Environment Simulator [JULES]; Best et al., 2011;Clark et al., 2011) to more complex schemes (e.g., the Community Land Model [CLM]; Lawrence, Fisher, Koven, Oleson, Swenson, Bonan, et al., 2019), with a wide variety of different schemes in use (Wania et al., 2013;Xu et al., 2016).
Due to the comparatively low spatial resolution of ESMs, whose grid cells can be hundreds of kilometers wide, the simpler wetland methane schemes often approximate the effect of water table by assuming that methane is produced only from the saturated fraction of the grid cell. This is because observations suggest that the majority of methane emissions occur when the water table is close to the surface ( Figure 2 in Turetsky et al., 2014) and methane fluxes to/from drier areas are generally very small, although recent studies show that the sink capacity of dry soils could be significant (Oh et al., 2020). More complex models such as CLM attempt to simulate the transport of oxygen into the soil, and oxidation and transport of methane through and out of the soil, which impacts how much of the produced methane reaches the atmosphere (Riley et al., 2011). Despite this variation in complexity, the model functional form for methane production is very similar across the whole current generation of ESMs, generally increasing exponentially with temperature, or following an Arrhenius function (Xu et al., 2016).
The temperature sensitivity of methane emissions at individual locations can be very high (Yvon-Durocher et al., 2014). However, the spatial temperature sensitivity (i.e., difference in emissions between on-average warmer and colder sites) is significantly smaller (Yvon-Durocher et al., 2014). This makes it challenging to parameterize the temperature sensitivity of methane production in ESMs. In Comyn-Platt et al. (2018), a range of different temperature sensitivities was used to span the uncertainty due to these apparently differing sensitivities.
In reality, methane is produced by microorganisms known as methanogenic archea or, less formally, methanogens (Bridgham et al., 2013) whose functionality is not explicitly modeled in current ESMs. Specifically, models do not account for the quantity and activity of methanogens, which can vary over time. The need for including biological processes in global soil models has been gaining more attention recently (Sulman et al., 2014;Walker et al., 2018;Wieder et al., 2013;Xenakis & Williams, 2014) and is a major focus of this study. Microbial dynamics have been included in several models for carbon dioxide production Xu et al., 2018) and also in recent models of methane emissions (Oh et al., 2020;Wang et al., 2019;Xu et al., 2015).
Our study focuses on the northern high latitudes, a region where around 40% of the world's wetlands are found (north of 50°latitude, based on Global Lakes and Wetlands Database [GLWD]; Lehner & Döll, 2004) and where future warming is likely to be higher than in any other part of the globe (Stocker et al., 2013). We also included a warmer site to test the wider applicability of our models and use independent evaluation data from around the globe.
We evaluate whether including more complex methane transport and oxidation processes, as implemented in CLM, allows us to reproduce the observed emissions better than simple functions of temperature and saturated grid cell fraction. We include a new module of microbial dynamics in an ESM land surface scheme (JULES) and demonstrate the ability of this new scheme not only to reproduce observed methane emissions but also to reconcile the difference between large-scale and local-scale temperature sensitivities. We make an initial assessment of the implications of this for the sensitivity of wetland methane emissions to global warming.

Temperature Sensitivity Terminology
In this study, we refer to temperature sensitivity through two alternative metrics: activation energy and "effective Q 10 ," which are used in the Arrhenius function (A, Equation 1).
AðT; QÞ ¼ Q 0:1T=ð1 − T=T0Þ (1) Here, T is temperature in°C, T 0 is the temperature at absolute zero (°C), and Q is approximately equivalent to the Q 10 parameter which is commonly used in ecological modeling. The original formulation of the Arrhenius function is written in terms of an activation energy (eV) instead of Q 10 (see Appendix A1 for how to convert between these two parameters).
In fact, the Arrhenius function specifically applies to single chemical reactions (Eyring, 1935) and therefore should not be applied to a biological soil system on time scales in which microbial growth and/or adaptation can take place. However, either this function or a (similar) exponential function is commonly applied on longer time scales to quantify the temperature sensitivity of biological systems (e.g., Turetsky et al., 2014;Yvon-Durocher et al., 2014), with the caveat that on these time scales, the system will not necessarily follow reaction kinetics. In fact, this discrepancy is exactly what we quantify in this study. So, for example, when we refer to "long-term" activation energy, we refer to the emergent relationship between temperature and methane emissions, not the literal activation energy of the underlying reactions. We nonetheless use the effective activation energy as a consistent metric across time scales to facilitate comparison.

Model Descriptions
We used two current state-of-the-art land surface models which are used in major ESMs: CLM and JULES. As well as running the standard versions of these models, modified versions of both JULES and CLM were also used (denoted JULES-microbe and CLM-prod.only, respectively) resulting in four different model versions in total. In this section, we first describe the methane production schemes in CLM and JULES models (sections 2.2.1 and 2.2.2); we then describe the new microbial model that we developed in JULES (section 2.2.3); and finally, we describe the modified CLM scheme that was used to test the impact of transport and oxidation processes (section 2.2.4).

Joint UK Land Environment Simulator
JULES is the land surface model used in the UK Earth System Model (UKESM) (Sellar et al., 2019). It is a community model that represents the surface energy balance, heat and water fluxes, snowpack dynamics, vegetation dynamics, soil biogeochemistry, and carbon and nitrogen fluxes Clark et al., 2011;Harper et al., 2016;Wiltshire et al., 2020). As well as being used in UKESM, JULES takes part in multimodel analyses such as the Inter-Sectoral Model Intercomparison Project (Rosenzweig et al., 2017) and the Global Carbon Project (Friedlingstein et al., 2019;Saunois et al., 2020) and has been used to make global projections, for example, of future hydrology, permafrost thaw, and carbon and methane emissions/feedbacks Chadburn et al., 2015a;Comyn-Platt et al., 2018;Gedney et al., 2019). JULES has been described in detail in the above references; here we include only a brief description of the methane emissions scheme.
The wetland methane scheme in JULES calculates methane emission from soil temperature (T) and substrate availability (C) (Equation 2), and this is then multiplied by grid box saturated fraction (calculated using a topographic index approach) to give the grid box methane emissions (Gedney et al., 2004). Recently, the scheme was updated to calculate methane production on multiple vertical soil layers and sum the flux from each layer to give total emissions (Comyn-Platt et al., 2018). The sum of methane production from each layer (i) is weighted by an exponential decay factor (τ), which empirically represents oxidation (see Equation 2 below). There are several options for the substrate to be used for methane production (Gedney et al., 2019), but for this study, as in Comyn-Platt et al. (2018), we assume that methane comes from soil carbon, thus neglecting root exudates. The overall scale factor for methane production (k 1 , Equation 2) has traditionally been tuned for each new JULES simulation in order to give global total wetland methane emissions in line with other estimates (e.g., from Saunois et al., 2016). However, in this study, we fit k 1 using site-level data and examine the resulting global emissions (see section 3.4).
The following is an equation representing methane emitted from the ith soil layer (saturated area only): Here, z i are depths at centers of layers, dz i are layer thicknesses, and A is the Arrhenius function (Equation 1). See Appendix A1 for more details.

Community Land Model
CLM (version 5) is a community model that accounts for the advanced full carbon and nitrogen cycle, described in Lawrence, Fisher, Koven, Oleson, Swenson, Bonan, et al. (2019). It forms the land surface component for the Community Earth System Model version 2 (CESM2; http://www.cesm.ucar.edu/models/cesm2/). CLM represents several aspects of the land surface, including surface heterogeneity, and consists of components or submodels related to land biogeophysics, the hydrologic cycle, biogeochemistry, human dimensions, crops, and ecosystem dynamics. The specific processes related to ecosystem productivity include vegetation composition, structure, and phenology; radiation; heat fluxes; soil and snow; hydrology; carbon-nitrogen cycling; dynamic land cover change; and crops with interactive management. See Lawrence, Fisher, Koven, Oleson, Swenson, Bonan, et al. (2019) for details.
The methane emission scheme in CLM includes processes for transport and oxidation in the soil and is described and evaluated in detail in Meng et al. (2012) and Riley et al. (2011). Methane production in each layer is calculated as a fraction of the aerobic respiration in the grid cell, multiplied with an additional exponential temperature dependence as well as some other factors, for example, a function of pH (Riley et al., 2011). The baseline aerobic respiration depends on soil temperature and carbon content in a similar way to the JULES scheme (Equation 2). However, aerobic respiration also depends on soil water contents with an optimum at a certain level of soil moisture, and this response is passed directly to the methane production, despite that in reality the "optimum" methane production occurs under full saturation (anaerobic conditions).
Methane emission calculations are performed twice, once assuming the average grid cell water contents and once assuming a 100% saturated soil column. The total emissions are then weighted according to the fractions of the grid cell that are considered saturated/unsaturated according to a topographic index approach. Note that the baseline respiration used in both calculations is based on the water contents of the unsaturated part, which may introduce spurious effects in the saturated fraction of the grid cell.
Oxygen content and soil methane concentration are vertically resolved and simulated separately for saturated and unsaturated grid box fractions. The methane oxidation rate is colimited by oxygen concentration and methane concentration. Methane can leave the soil via ebullition, diffusion, and aerenchyma transport. Similarly, oxygen can enter the soil via diffusion and aerenchyma transport. Oxidation rate also increases with temperature according to a Q 10 function. This means that the fraction of produced methane that gets emitted to the atmosphere can vary dynamically depending on all of these processes (Riley et al., 2011).

JULES-Microbe
A microbial model of methane production was developed and subsequently implemented in JULES. JULES coupled with the microbial model is referred to as JULES-microbe. Methane is produced under anoxic conditions by a group of microorganisms known as methanogenic archea, or methanogens (Garcia et al., 2000), which are the microbes to which our microbial model refers. This section briefly describes the model and more details are given in Appendix A.
The basic structure of the microbial model is illustrated in Figure 1. Organic material in the soil is gradually decomposed to dissolved substrate. Dissolved substrate is then converted to microbial (methanogenic) biomass, methane, and carbon dioxide. Methanogenic biomass can also decrease due to mortality (i.e., cell lysis), and dead biomass is returned back to the dissolved substrate pool. We focus on soil organic carbon as a substrate source for proof of concept and intend to include root exudates in the future (see section 4).
Decomposition of soil organic matter to dissolved substrate happens in several stages (Christy et al., 2014), but hydrolysis has been shown to be the rate-limiting step (Mata-Alvarez et al., 2000). It is not performed by the methanogens but rather by hydrolytic microorganisms (Christy et al., 2014). We choose to model it as purely a chemical process based on the evidence from Veeken and Hamelers (1999) (see Appendix A2).
Since the default JULES methane emissions scheme is also modeled as a chemical process, we are able to set the parameters to give the same total soil carbon decomposition in both JULES and JULES-microbe. The difference is that in JULES, the decomposed organic matter is assumed to be transformed immediately to methane without the intermediate dissolved phase.
The rate of consumption of dissolved substrate depends on the microbial biomass and the rate of microbial respiration, which itself increases with temperature, substrate availability/concentration, and the level of microbial activity (see Equation A2). The consumption of the substrate is then partitioned between growth of methanogenic biomass (with a carbon use efficiency parameter) and methane and carbon dioxide emissions (50% to each).
The mortality rate increases with temperature, and the level of microbial activity increases to a maximum in good growth conditions (i.e., high growth rate, which occurs under warm conditions with high substrate) and reduces to a minimum level in poor growth conditions, representing dormancy. Mortality rate reduces with activity level as well as respiration/growth rate, which means that dormancy provides an advantage for survival (see Appendix A1). Additionally, for frozen soil, when the liquid water content falls below a certain level, water is no longer able to move between pores (Olesen et al., 2001). Therefore, below this threshold water content, we make a reduction to substrate production (see Appendix A1).
All parameters in the model are estimated from observational studies except for the overall decomposition scale factor (k 1 ) and the oxidation parameter (τ) , both of which are also used in the original JULES scheme (Equation 2) and we optimize to fit site-level observations (section 2.3.4). We later use CLM-prod.only (section 2.2.4) to compare the fitted oxidation parameter (τ) with the process-based oxidation scheme used in CLM (notably explicit oxygen transport and supply and explicit transport of methane out of the soil; section 2.2.2), and we use independent, globally distributed methane emission observations to validate the overall scale factor (k 1 ) (see sections 3.1 and 3.4).
The full equations of the microbial model and the full set of parameter values are given in Appendix A1. The source code of the microbial scheme is available in a JULES branch (see Data Availability Statement).

CLM-Prod.Only
In order to test the impact of the detailed representation of methane transport and oxidation processes in CLM, a modified model version was also analyzed, where only the production of methane in the saturated part of the grid cell was used from the original scheme, and then instead of applying transport and oxidation processes, the simple oxidation function from JULES (expð−τzÞ, Equation 2) was applied. Emissions from the unsaturated part of the grid cell were set to zero.
This bridges the gap between JULES and CLM by making a "JULES-like" version of the CLM methane emissions scheme but which has methane production provided by CLM and is therefore directly comparable with the CLM model results.

Site-Level Data and Simulations 2.3.1. Site-Level Data
We used detailed site-level observations of methane emissions from seven different sites ( Figure S1), all of which include at least some wetland area. The sites span high-latitude climatic conditions as well as including one midlatitude site (a 22-year-old restored tule/catail freshwater wetland on Twitchell Island, CA). The sites are described in detail in Appendix B. Only in relatively recent years has it been practical to run eddy covariance to measure methane emissions (Hendriks et al., 2007;Peltola et al., 2019). This is extremely valuable since methane is produced by soil microbial activity, and therefore, these data can give an insight into what is going on below the ground. Where sites were not 100% wetland, we filtered the eddy covariance data according to wind direction when wetlands were in certain directions, or divided it by the average wetland fraction when the footprint was mixed. Note that by "wetland," here, we would ideally refer to areas where the water table is close to the surface of the soil, to be as consistent as possible with how the models define the saturated fraction (see section 2.3.4 below). In practice, we were only able to select relevant land cover types, for example, bogs/fens. A time series of more than 1 year of methane emission data was available for most sites.
Soil carbon profiles were measured at all sites, with varying resolution (see Table B1). Each site also measured a continuous time series of soil temperatures (hourly or half-hourly) at multiple soil depths (thus covering every variable in Equation 2). Soil temperatures were taken from wet areas in order to be representative of the conditions where methane is produced. Note that soil temperatures from a dry spot in the same site can be very different (see, e.g., Göckede et al., 2017).
Meteorological data were available from local met stations or weather stations, and continuous snow depth data were used to estimate the snowfall, based on an estimated "fresh snow density," as snowfall is difficult to measure directly (see Chadburn et al., 2017, for more details). The data available from each site are listed in Table B1.

Site Forcing and Ancillary Files
Meteorological forcing data to run CLM and JULES were prepared using observations from the respective sites combined with reanalysis data for the grid cell containing the site. The WATCH Forcing Data and WATCH Forcing Data ERA-Interim (Weedon et al., 2014) were corrected to match the local meteorological stations at the sites and combined to provide a continuous time series of forcing from 1901 until the end of 2016. The method is described in Chadburn et al. (2017).
Soil thermal and hydraulic properties for JULES were derived from the measured profiles of soil carbon, by first using carbon content to estimate the volumetric fraction of organic matter in the soil and second by combining organic and mineral soil properties as described in Chadburn et al. (2015b). Mineral soil parameters were estimated based on the soil types described in site publications (see Appendix B for references). CLM similarly used the observed soil carbon profiles including organic matter density and soil texture (sand and clay ratios) at 10 soil levels down to 3.8 m to determine soil thermal and hydraulic properties, following Lawrence and Slater (2008) and Letts et al. (2000).

Model Configuration and Spin-Up for Sites
The four land surface model versions (JULES, JULES-microbe, CLM, and CLM-prod.only) were run with the site-specific forcing and ancillary data described in section 2.3.2. Details of the model configurations are given in Table 1. JULES was set up following the UKESM configuration (Sellar et al., 2019;Walters et al., 2019), but with some alterations (see Table 1). Both JULES and JULES-microbe were spun up to equilibrium by repeating the forcing data from 1901 to 1920 and then run in forward mode until the end of 2016. Additionally, the simulations were re-run forcing the soil moisture to maintain saturated conditions, to better represent wetlandsat the expense of correctly simulating the non-wetland fraction of the grid cell. For the future, soil tiling approaches are being considered such that both could be simulated at once (e.g., Swenson et al., 2019 Riley et al., 2011). Additionally, the model implemented a simplified inversion approach to dynamic representation of spatial inundation fraction for methane production based on Prigent et al. (2007). Inversion parameters were prescribed at 0.9°× 1.25°resolution for calculating inundation fraction based on total water storage. The meteorological forcing data with height set to 2 m were cycled over the years 1901-1920 to spin up the model to equilibrium in the "accelerated decomposition" mode (Thornton & Rosenbloom, 2005). The model then continued in normal mode until the end of 2016.

Analysis and Parameter Fitting
For JULES, the k 1 and τ parameters in Equation 2 were fitted by running the JULES methane scheme (essentially, Equation 2) using observed soil carbon and soil temperature, for the sites Lompolojänkkä, Abisko, and Twitchell. A full plausible range of k 1 and τ values were tested and chosen to optimize the 3.8 3.8 3.0 3.0 Soil hydrology levels 10 10 14 14 Vertically resolved oxygen concentration Koven et al. (2013). b Koven et al. (2017).
annual mean methane emissions, by minimizing root mean squared error (RMSE). These sites were chosen due to having 100% wetland fraction at least in certain wind directions-a mixed footprint would introduce additional uncertainty-and methane observations covering close to a full year continuously (or more), so that an annual total could be calculated. Parameter values are listed in Table A1 and optimization results shown in Figure S4. Note that the same value of τ is also used in CLM-prod.only, and the same values of both k 1 and τ are used in JULES-microbe where they represent the same processes.
We validated the microbial scheme by similarly using observed soil carbon and soil temperature to run the JULES-microbe scheme. All seven sites were used. To combine observed soil carbon and temperatures, soil temperature profiles and carbon profiles were matched up using an interpolated profile of soil carbon and dividing it into layers based on the vertical depths at which the soil temperatures were measured. This removed model errors in simulating soil carbon and temperature and allowed a more direct validation of the methane emissions.
The models represent wetlands as a "saturated fraction" of the grid cell which updates every time the grid cell gets wetter or dryer. This differs from the real-world concept of a wetland ecosystem, which is much more spatially static but whose water table level may fluctuate. We attempted to reconcile these differences and use comparable quantities from models and observations. In particular, we used the methane emissions per square meter of "saturated" area in the models and similarly from observations-by estimating as far as possible the area in the flux tower footprint for which the water table would be close to the surface, which is what the modeled saturated fraction represents. JULES, JULES-microbe, and CLM-prod.only did not simulate any emissions from the non-saturated fraction of the grid cell. For CLM, the emissions per square meter were typically calculated by dividing the total grid box emissions by the simulated saturated fraction. For some sites, CLM simulated soil moisture to be essentially saturated for the whole grid cell and thus produced methane almost equally from the whole grid cell but (erroneously) did not calculate a saturated fraction close to 100%. In that case, emissions from the whole grid cell were considered. Furthermore, the saturated fraction at Lompolojänkkä was simulated as zero. In this case, the full methane dynamics are still simulated for the saturated part of the grid cell (even though its areal fraction is zero), but they are not included in the grid cell output. The method of dividing grid cell emissions by saturated fraction clearly could not be applied here, so for this site, methane emissions were estimated as net production (production − oxidation) in the saturated fraction.
To analyze temperature sensitivity of model outputs and observations for the sites, we estimated activation energy on a seasonal time scale by fitting the Arrhenius equation (Equation 1) to the model or data in question. This was done using monthly mean temperature and log of monthly mean methane emissions. Monthly means were used because the observations would be expected to have a larger day-to-day variability than the models, for example, due to atmospheric turbulence (Wille et al., 2008) which is not modeled, or due to measurement errors, and that would lead to regression dilution in the observations relative to the models when the seasonal activation energy was calculated.

Large-Scale Data and Global Simulations
Global simulations were performed with JULES and JULES-microbe following the same configurations as in Table 1. Global atmospheric forcing data were from CRUNCEP. CRUNCEPv7 is a combination of 6-hourly NCEP data from 1901 corrected to CRU climatology and observations updated to 2014 (Harris et al., 2014). The models were spun up for 300 years with CRUNCEP data from 1901. The soil was set to be permanently saturated and the soil thermal/hydraulic properties were prescribed to resemble peat (see Table A2), thus giving a realistic temperature simulation for wetlands similar to what was done for the JULES site simulations (see section 2.3.3). Not all wetlands will have high organic matter contents, but without spatially resolved information on wetland soil properties, we decided that using this soil type would be the most reasonable choice, as it was the case at our seven test sites.
In order to estimate global total methane emissions, the emissions per square meter of saturated fraction were output from the model and post-processed with the Surface WAter Microwave Product Series (SWAMPS)-GLWD data set (Poulter et al., 2017), in order to reduce bias caused by errors in the model simulation of saturated area. SWAMPS-GLWD is an integration of the SWAMPS (Schroeder et al., 2015) and the GLWD (Lehner & Döll, 2004), thus incorporating both ground-based and satellite-based information (Poulter et al., 2017). While such an approach will not remove all errors, SWAMPS-GLWD is the data set used in the Global Carbon Project methane budget assessments (e.g., Saunois et al., 2020), so any analysis we make using this data set should be comparable to their results. The mean inundated area in SWAMPS-GLWD is around 4.6 million km 2 . It provides gridded monthly inundated areas from 2000 to 2012, so we calculated the global total emissions from JULES and JULES-microbe for the same time period, assuming that this inundated area corresponds to the saturated grid cell fraction from the models. In order to calculate the long-term trend in global methane emissions, we applied the mean seasonal cycle of inundated fraction from SWAMPS-GLWD outside its measured time period, effectively assuming that there was no long-term trend in wetland area (see section 3).
For other analyses of global simulations, emissions per square meter of saturated fraction were considered and SWAMPS-GLWD was therefore not required. For independent model validation, site-level but globally distributed annual mean wetland methane emission/air temperature data were combined from Yvon-Durocher et al.   transmission of seasonal variability to deeper soil ( Figure 2, top right panel), giving a more realistic temperature simulation. Therefore, for the rest of the analysis, we use the saturated JULES version. The seasonal cycle of soil temperature in CLM is too large, but overall, we conclude that the soil temperature simulations in both models are reasonable (see also Figures S2 and S3). The Twitchell site is cooled by input of both ocean water and meltwater from nearby mountains, so overestimation of soil temperature by the models is expected (note that soil temperatures in the JULES-microbe and CLM-prod.only simulations are identical to their parent models; only the methane scheme differs). Figure 3 shows the methane emissions per square meter of saturated area from all four model versions compared directly with the observations for all seven study sites. There are clear differences in the magnitude of the emissions both between model versions and between models and observations.

Comparing Land Surface Schemes With Detailed Site-Level Observations
We first compare the green and orange lines in Figure 3, which represent CLM and CLM-prod.only. CLM includes the detailed transport and oxidation schemes and CLM-prod.only does not (see section 2.2.4), so this comparison indicates the impact of the transport and oxidation processes. There are significant differences in magnitude between these two model versions at a few of the sites, and CLM generally has a higher short-term variability than CLM-prod.only (more "spiky" emissions). However, there are no major differences in the seasonal dynamics. In fact, the correlation coefficient of the weekly mean methane emissions in CLM and CLM-prod.only is between 0.94 and 0.98 (R 2 0.89-0.97). Thus, we conclude that including transport and oxidation processes and/or emissions from the unsaturated part of the grid cell does not significantly improve the predictive power of the model on a seasonal time scale, which is consistent with the findings of McNorton et al. (2016). This is also consistent with studies such as that by King et al. (1998) which shows that adding or removing aerenchyma for methane transport impacts the magnitude of the emissions, but that the temporal dynamics are similar.
The oxidation rate given by the fitted τ parameter (τ ¼ 6:5, Table A1) is apparently reasonably consistent with the more sophisticated oxidation scheme in CLM, since the magnitudes of methane emissions in CLM and CLM-prod.only are similar at the majority of sites. This parameter value results in oxidation of around 70% of the produced methane at the sites, which is consistent with observational estimates that 40-70% (Bridgham et al., 2013) or 60-90% (Popp et al., 2000) of produced methane is oxidized. 10.1029/2020GB006678

Global Biogeochemical Cycles
The land surface model simulations of methane production depend strongly on both the underlying carbon source and the soil temperature (e.g., Equation 2). We have seen that soil temperatures are reasonable ( Figure 2). However, the model simulations of soil carbon show anything from between less than half the amount, to more than twice as much soil carbon as the site observations. Figure 4 shows the error in the annual mean methane emissions compared with the error in the simulated soil carbon. These fall close to the 1-1 line which shows that for the sites where soil carbon is correctly simulated, the models are able to simulate approximately correct methane emissions (correlation coefficients of 0.91, 0.75, 0.88, and 0.84 for CLM, CLM-prod.only, JULES, and JULES-microbe, respectively). Thus, the soil carbon appears to be a major source of error (SD of percentage error ∼74%), but the annual mean methane emission rate per unit of soil carbon is better simulated (SD of percentage error ∼37%).
While the simulation of annual mean emissions is reasonable in the models, there are clearly major differences in the seasonal dynamics ( Figure 3). Both versions of CLM display a peak in methane emissions much earlier in the season than the observations, except at Seida, where the observed methane emissions also have an early peak. In addition, the winter emissions in CLM are too low. On the other hand, the winter emissions in JULES are too high, especially compared to the height of the summer peak (i.e., the seasonal cycle is too small).
These seasonal dynamics are seen clearly in Figure 5, which shows the relationships between methane emissions and soil temperature. The simplicity of the JULES scheme is apparent since the normalized emissions are related directly to temperature. In the observations, there is a clear tendency for lower methane emissions in the first half of the summer season than the second half (a "lag" in emissions). This is not recreated by JULES or CLM (in fact, CLM shows the opposite). However, by including the dynamics of the methanogens, JULES-microbe is the only model that is able recreate the observed dynamics.
The methanogens in JULES-microbe are in the long term limited by the total production of dissolved substrate, whose calculation is essentially the same as the calculation of methane emissions in JULES, which explains why the annual mean emissions from the two models are similar ( Figure 4). However, the microbes tend to go dormant over winter, so substrate accumulates and the methanogens then have extra substrate to use during summer, which is why the summer emissions are then greater in JULES-microbe than in JULES. Since the methanogens have gone dormant over winter, it takes time for methanogens to reactivate in the spring, which introduces a lag in methane emissions and can explain the lag pattern seen in both the observations and JULES-microbe in Figure 5.
In CLM, the methane emissions are extremely small below the freezing point, essentially dropping down to zero. This is not seen in the observations, which exhibit gradually declining emissions down to at least -15°C ( Figure 5). On the other hand, JULES-microbe follows the observed behavior reasonably well ( Figure 5).
Freezing is known to cause significant reduction in microbial biomass (Kumar et al., 2013), but this impact is highly variable and it tends to regrow after the initial reduction (Segura et al., 2017). Therefore, we did not include this in the JULES-microbe model and the results nonetheless seem reasonable. Only one site (Chersky) has significant data during very cold conditions, so more data are needed to confirm this.

Temperature Sensitivity on a Site Level
JULES-microbe appears to recreate the observed relationship between temperature and methane emissions in Figure 5. In order to quantify this, we applied the Arrhenius function to these data (notwithstanding that there are lags so it won't be exact) and calculated the effective activation energy on a seasonal time scale, which can also be written as an "effective Q 10 " (Equation 1). This is shown in Figure 6 and supports the conclusion that the temperature sensitivity in JULES-microbe is closest to the observed sensitivity. It is also consistent with the estimates from Yvon-Durocher et al. (2014). The activation energy for CLM comes out extremely high due to the drop-off of emissions to zero during winter, but the activation energy of summertime emissions alone shows that the temperature sensitivity is in fact too low during summer in CLM.
We conclude that the JULES-microbe model is best able to recreate the observed dynamics, both of effective activation energy and the lagged response between the temperature changes and methane emissions. Thus, adding methanogens into a very simple methane scheme is sufficient to reasonably reproduce observations. The underlying errors in soil carbon and soil temperature in JULES impact the magnitude of methane emis-

10.1029/2020GB006678
Global Biogeochemical Cycles sions, so that while the temperature response is recreated well by JULES-microbe ( Figures 5 and 6), the magnitude of emissions is too high or low at some sites ( Figure 3). Therefore, in Figure 7, we apply the JULES and JULES-microbe methane schemes using the observed soil temperatures and soil carbon profiles to drive the model. Here it becomes clear that the original JULES scheme can simulate the correct magnitude of methane emissions at each site (i.e., recreating the spatial variability), but its seasonal cycle is too small, and that adding the microbial scheme results in a more accurate recreation of seasonal dynamics at almost every site. This is also clear in Figure S6 which shows the same data as a plot of observed versus modeled values. Table 2 shows that RMSE is reduced for every site except Twitchell (which we discuss in section 3.3).
In JULES-microbe, methanogenic dynamics recreate the hysteresis in the observed temperature response. However, the pool of labile carbon also changes through the seasons, so could this instead explain the observed response? Carex spp. and Eriophorum spp. are a major source of root exudates which methanogenic archaea convert into methane. However, these perennial sedges grow quickest after snowmelt, when the amount of incoming solar radiation is high. Substrate availability has been measured to be highest at the beginning of the growing season, diminishing toward the autumn (Ström et al., 2012). Moreover, freeze-thaw processes during winter may degrade organic matter, adding to the pool of labile carbon early in the growing season (Mastepanov et al., 2013). It is therefore more likely that substrate availability would cause the opposite pattern, with a higher temperature sensitivity in early summer. The same is true for the dynamics of methanotrophs, which may increase oxidation rate-and thus reduce emissions-in the later part of the summer (Wagner et al., 2003). To further evaluate model performance and to ask the question of whether the JULES-microbe model gives the right answers for the right reasons, we compare the model against available microbial-related data in Appendix C, which provides additional support for our conclusions.

Water Table Dynamics
It is apparent in Figure 7 that at Twitchell, the microbial dynamics improve winter emissions but lead to a large discrepancy in spring. At this site, a time series of water table depth has been recorded. While the water table depth is generally maintained at a level a little above the soil surface,  Global Biogeochemical Cycles the observations show that it occasionally drops lower, which could impact the methane emissions. We tested a simple scheme where methane emissions from any layer were set to zero when the water table reached within 6 cm above that layer, and microbial growth rate was similarly set to zero (which induces dormancy). The 6-cm threshold was chosen because oxygen supply is greater than 50% of its unlimited value down to around 6 cm below the water table (Frei et al., 2012). It is known that substrate availability is reduced and methane production suppressed during drying and for some time after rewetting (Segers & Kengen, 1998). This "lag" period ranges from around 10 to 20 days in Segers and Kengen (1998). We set the recovery period as 12 days because this gave the lowest RMSE between model and observations. In addition, we dilute the substrate when the water level increases since it is known to come from inflow from nearby mountains and fresh estuary water. These changes removed a large part of the discrepancy in spring (Figure 7), and looking at the time series, we see that during the first few years when major water table fluctuations occurred, the dynamics are noticeably improved ( Figure S5) and overall RMSE is reduced ( Table 2). While this is a promising indication that where water table dynamics are known, their impact could be modeled, the question of how to apply this in a global model requires substantial further investigation. Changes in water table in global models are currently represented by a change of saturated area. This means that while methane emissions can be scaled by saturated fraction, the microbial dormancy and biomass changes that occur in our Twitchell scheme and their longer-term effects on the dynamics after rewetting cannot yet be captured in global runs.
Since the inclusion of water table fluctuations lowers methane fluxes at Twitchell in the early growing season, this raises the question of whether the hydrology could be an alternate explanation for the lag seen at the rest of the sites. However, we don't expect the same at the Arctic sites included in this study, since snowmelt leads to high moisture levels early on in the growing season and this suppresses methane oxidation. Following this initial wet state, the water table generally lowers due to evapotranspiration and improved drainage through active layer deepening, only to rise quickly following rain storms (see, e.g., Parmentier et al., 2011;Sachs et al., 2010). Inclusion of water table dynamics would arguably improve some dynamics within the season in response to a long drought or intense rain but does not explain the lagged response in emissions.
Since Twitchell is somewhat of a "special case" among the test sites, with stronger hydrological dynamics than the other sites, it is more difficult to confirm the applicability of the microbial model outside the high latitudes. In fact, larger seasonal water table fluctuations are more likely to occur in lower latitudes where both precipitation and evaporation rates are commonly higher. Therefore, including these dynamics may be necessary for accurately modeling warmer sites. Through a parameter sensitivity test (Figures S9 and S10), we find that the model is sensitive to the parameter choices at Twitchell, although the simulation is generally improved in winter compared with the original JULES model, and often improved in summer, even without considering the water table impacts ( Figure S10). Keeping this in mind, we perform an initial investigation of global methane dynamics with the JULES-microbe model in the following section.

Temperature Sensitivity on a Global Scale
When running the JULES and JULES-microbe models globally, we find that the simulated methane emissions are somewhat higher than observation-based estimates, which is probably to be expected since the model is missing hydrological impacts that act to reduce emissions (as shown in section 3.3 above). As discussed, the global simulation was run with forced saturated soils and peatland thermal/hydraulic properties for the whole globe. The resulting emissions for each grid cell were multiplied by the inundated fraction in the grid cell according to SWAMPS-GLWD. This resulted in an annual wetland emission of around 260 TgC/yr in the period 2000-2012. This is above the upper end of the current estimates for wetland methane emissions in this time period, approximately 130-230 TgC/yr (Saunois et al., 2016). However, we showed in section 3.3 that including water table dynamics was necessary to accurately simulate methane emissions at Twitchell. In this scheme, and in general (Zhong et al., 2020), a drop in water table below the soil surface leads to reduced methane emissions. Therefore, our global simulations, which assume that the water table is always at the surface, are necessarily an upper bound on the possible wetland methane emissions. Furthermore, there are types of wetlands that emit little methane, for example, due to pH effects which can inhibit the acetate fermentation pathway (i.e., reducing the potential sources of methane production), or due to the inhibitory effect of sphagnum mosses forming a symbiosis with methanotrophs (Bridgham et al., 2013;Kip et al., 2010). Such wetlands would be included in SWAMPS-GLWD, again explaining why our model may overestimate wetland methane emissions. Nonetheless, the fact that extrapolating parameters that were fitted at a site level to the whole globe results in a global total emission of the right order of magnitude suggests that these parameter values are reasonable. Figure 8 shows a large number of annual mean methane emissions from a wide range of observational studies and compares this to every single grid cell from the JULES and JULES-microbe global simulations, in terms of emissions per square meter of saturated fraction. The simulations are within the range of these observations, while they are toward the upper quartile, which is consistent with the overestimation of the global total discussed above. Importantly, the temperature sensitivity of methane emissions spatially-that is, across the globally distributed site observations (activation energy 0.53 eV)-matches well with the temperature sensitivity simulated by JULES-microbe (activation energy 0.49 eV). JULES exhibits a somewhat lower spatial temperature sensitivity (activation energy 0.39 eV). These activation energies are smaller even that the underlying sensitivity of the decomposition of organic matter (to form dissolved substrate) in the model (0.64 eV), which is to be expected because soil carbon is generally lower in warmer soils (Koven et al., 2017), which partly counteracts the increase of decomposition with temperature.
In contrast, Yvon-Durocher et al. (2014) showed that the activation energy of methane emissions on a seasonal basis has a much higher value of almost 1 eV (0.86-1.07 eV). Figure 9 shows that the effective seasonal activation energy in JULES-microbe is close to 1 eV in many areas. Furthermore, the vast majority of sites used in Yvon-Durocher et al. (2014) (marked on Figure 9b) coincide with areas that have simulated activation energies of 0.8-1.1 eV in JULES-microbe, and thus, the model appears to be consistent with the analysis of Yvon-Durocher et al. On the other hand, in JULES, the seasonal activation energy is generally less than 0.6 eV ( Figure 9a). Thus, microbial dynamics resolve the apparent discrepancy between local temperature sensitivity and large-scale spatial patterns of methane emissions. Figure 10 shows the long-term (60-year) trend in the globally averaged methane emissions and soil temperatures. This shows a clear increase of methane emissions with temperature which is equivalent to

10.1029/2020GB006678
Global Biogeochemical Cycles an activation energy of 0.57 eV for JULES and somewhat stronger response equivalent to 0.71 eV from JULES-microbe. The response in JULES-microbe is small compared to its seasonal-scale activation energy (0.7-1.1 eV, Figure 9b) but implies approximately a 12% increase in methane emissions for every degree of warming of the land surface. This means that one degree of warming from present day would lead in an order of 20 TgC of additional methane emissions per year or around 0.6 GtC a year in CO 2 equivalent on 20-year time scale.
While this does not account for hydrological changes, we refer to Gedney et al. (2019), who analyzed in detail the possible changes in wetland extent under future climate change scenarios. There was not a consensus as to whether the change in wetland extent would act to increase or decrease methane emissions in the future (Gedney et al., 2019; Figure S7). Under all scenarios, the impact of wetland extent was much smaller than that of temperature change, by around an order of magnitude. This indicates that large-scale hydrological dynamics are unlikely to substantially alter our conclusion.

Conclusion and Outlook
Our study shows that two current state-of-the-art land surface models (CLM and JULES) can recreate the magnitude but not the observed seasonal cycle of methane emissions at a range of wetland sites. The seasonal cycle in JULES is too small, and the seasonal cycle in CLM is too large (due to emissions dropping to zero in winter), and in fact, the temperature sensitivity of methane emissions in summer is too small in both models, as shown by their effective activation energies ( Figure 6). We have shown that JULES-microbe, a relatively simple model of microbial dynamics that represents the growth and activity of methanogens, is better able to simulate the seasonal cycle of methane emissions at the study sites (see, e.g., Figure 5). The methanogens use less substrate than is produced during winter and thus have extra substrate for growth during the summer, which leads to an amplification of the seasonal cycle in line with observations. This results in a high apparent activation energy of methane emissions on a local, seasonal basis ( Figure 6).
The microbial model accounts for the discrepancy between observed seasonal temperature response and the spatial temperature dependence of methane emissions (Yvon-Durocher et al., 2014). The long-term global temperature response of methane emissions ( Figure 10) is equivalent to an effective activation energy of around 0.7 eV, which falls somewhere between the seasonal and spatial values, and suggests that JULES-microbe would predict future methane emissions toward the lower end of recent estimates in Comyn-Platt et al. (2018).
While the ESM methane schemes are able to produce reasonable results for annual total emissions at each site when the soil carbon quantities are correctly simulated, the models do not simulate soil carbon well (ranging between less than half to more than double the observed values, Figure 4). This highlights a priority for future model development.
In addition, we showed that including a simple approximation of water table impacts led to a better simulation of the methane emission dynamics at Twitchell Island (Figures 7 and S5). Note that a drop in water table is currently represented in global models as a change in saturated area, and the impact of changing water table on microbial dynamics that we simulated at Twitchell could not be captured by this scheme. One possible solution would be to represent wetlands as separate "tiles" within a grid cell (which can receive water from surrounding tiles; see, e.g., Swenson et al., 2019). Such a scheme would also benefit the simulation of soil carbon since non-linear processes (e.g., aerobic vs. anaerobic respiration) mean that soil carbon in wetlands can be very different from the grid cell average. The models are also missing potentially important controls on substrate dynamics such as the input of dissolved substrate from root exudates (see, e.g., Ström et al., 2003) and leaching of substrate from the soil. Root exudates are simulated as part of the nitrogen scheme in JULES, and leaching is simulated in the separate DOC model. Therefore, linking nitrogen, root

10.1029/2020GB006678
Global Biogeochemical Cycles exudates, and leaching to the microbial methane scheme will be possible in the near future. Note that decomposition of fresh plant litter has been shown to be a large component of methane emissions (Chanton et al., 1995) and this is included in JULES and CLM. Evaluation of all of these processes is, however, still limited by the scarcity of observational studies that partition these different carbon fluxes.
There are also many other microbial functional groups that may be important for global biogeochemical cycling. A priority is to incorporate microbial processes in the aerobic part of the soil model (i.e., for carbon dioxide production) (Walker et al., 2018). Methanotroph dynamics were not investigated in this study, although we note that Yang et al. (2017) considered both processes and found that methanogenesis drove the landscape-level methane flux. Methanotroph community dynamics would likely oppose the effect of methanogens on the seasonal cycle of methane emissions, since methanotrophs would also be more active during the summer (e.g., due to increased methane and oxygen supply). The interaction between methanogenesis and methanotrophy is complex (Yang et al., 2017), and few studies measure their dynamics separately. However, in principle, we could extend our work to explicitly model both communities in the future. Further investigation of these and other microbial functional groups is undeniably important.
There are few direct measurements of methanogenic biomass, substrate, and so forth in conjunction with methane emissions that can be used to directly validate the microbial model. Our study provides evidence that microbial dynamics can produce qualitatively different model behavior and therefore provides more incentive to better measure and understand their dynamics and to include them in ESMs.

A1 Model Equations
In this section, the equations are laid out and their terms defined. In the following section (section A2), the rationale behind the choice of functions and parameters is explained.
Methane production in JULES (without the microbial scheme) is calculated as follows: the flux (F CH 4, i , kgC/m 2 /s) from the ith soil layer is where z is vertical depth in the soil (m) (z ¼ 0 at soil surface and takes positive values below), dz i is the thickness of layer i, and τ (m −1 ) is a decay factor determining depth dependence of oxidation. C i (kg/m 2 ) is the total grid box soil carbon in layer i after weighting the individual carbon pools (DPM …HUM, Figure 1) by their relative decomposition rates, such that the weighting factors add up to 1. Q is a parameter determin- Note. The first three rows are used in both JULES and JULES-microbe. The remaining seven parameters are introduced to determine microbial dynamics. Since only the total soil carbon was available from observations, but the JULES scheme uses four soil carbon pools weighted by their decomposition rate, k 1 had to be scaled based on the average relative pool sizes in JULES.

10.1029/2020GB006678
Global Biogeochemical Cycles ing temperature sensitivity, T i (°C) is the soil temperature in layer i, and A is the Arrhenius function describing dependence of reaction rates on temperature (Equation 1). k 1 is an overall scale factor which is traditionally tuned to give realistic global emissions but in this study was fitted to site-level data (section 2.3.4). The flux is summed over all soil layers and in a large-scale simulation is then multiplied by the grid box saturated fraction to give total emissions.
The microbial model used in JULES-microbe introduces three additional variables, dissolved substrate (S), methanogenic biomass (B) (both in kgC/m 3 ), and methanogenic activity (m act ), which represents a fraction of the maximum possible activity level at any time. These variables are calculated for each soil layer. The equations for substrate and methanogenic biomass dynamics are where the first term represents substrate production via hydrolysis of soil organic matter and the second term is the substrate consumption by methanogens, and where the first term is growth, the middle term is mortality, and the last term is maintenance respiration. Mortality is set to zero when the growth term is greater than the maintenance term. C is soil carbon density in kgCm −3 . k 2 and k d are baseline rates of respiration and mortality (respectively) at 0°C (hr −1 ). f is the carbon use efficiency. ϕ(S, T) is the substrate limitation term (discussed below). α is the maintenance respiration as a fraction of the maximum respiration. A 1 and A 2 are the Arrhenius function (Equation 1) with different Qs (given here, two alternative forms): where, in the first version of the equation, T is soil temperature in Celsius, T 0 ≃ −273.15°C is the temperature at absolute zero, and the Q j parameters (Q 1 and Q 2 ) are approximately equivalent to the often-used "Q 10 ." On the right-hand side of the equation, k is the Boltzmann constant and E M, j is the activation energy in eV. In the analysis, both activation energy and "effective Q 10 " (¼ Q j ) are often shown, so, for reference, the conversion between E M and Q j is E M; j ¼ 0:1T 2 0 k lnQ j . The methanogenic activity level, m act is downregulated or upregulated based on growth conditions. The minimum activity level is such that the potential rate of carbon use for growth/maintenance matches the maintenance level (f m act; min ¼ α, see Equation A3), the maximum is m act ¼ 1, and the dynamics of m act are where "growth" is the methanogenic growth rate (first term in Equation A3) and the increase/decrease of m act stops if it reaches its maximum/minimum value. Note that the maintenance respiration (αk 2 A 2 B in Equation A3) does not depend on m act and is therefore never downregulated (since it is the minimum respiration rate for survival). When growth <μ, net growth of methanogenic biomass is also set to zero JULES-microbe also includes acclimation of both respiration and mortality based on temperature, which acclimates slowly in response to multiannual temperature change. The acclimation factor is updated each time step (typically around 60 minutes), allowing it to diffuse slowly toward its "instantaneous" value (d a mean =dt ¼ ða instant − a mean Þ=ev, where ev ¼ 5 years). The acclimation term modifies all of the microbial activities, so namely, k 2 ¼ k 2;0 a, k d ¼ k d; 0 a, and where Q ev ¼ 0:55Q 2 (Q 2 being the Q parameter in the Arrhenius function (A 2 ) used in Equation A3).
When the soil is frozen below a certain volumetric liquid water content (θ) given by θ ¼ 0:02b (Olesen et al., 2001) (where b is the exponent in the hydraulic suction curve), substrate production (k 1 ) is reduced by a constant factor of f rz ¼ 0:5.
Finally, the methane production for each soil layer, i, is set to half of the substrate consumption (second term in Equation A2), based on the theoretical 50-50 split between carbon dioxide and methane produced from such a system (Conrad, 1999): This results in a total saturated-area methane emission flux as the sum of methane production from all soil layers, multiplied by the oxidation factor: Note that after carbon from the soil carbon pool (C) enters the substrate pool (S), the carbon balance is maintained by assuming that substrate consumption (second term in Equation A2) is partitioned between growth of methanogenic biomass (first term in Equation A3) and methane emission (Equation A7) and that the remaining balance is emitted at CO 2 . Carbon balance is maintained within the microbial pools by returning necromass back to the substrate pool.

A2 Rationale and Parameter Choices
The new parameters and functions introduced into JULES-microbe are all estimated from observations, as we will now discuss. A list of parameter values and references is given in Table A1. We also performed a sensitivity study to investigate the impact of these parameter choices ( Figures S9 and S10). Varying the parameters by ±25% showed that the parameter controlling the maximum growth rate of the methanogens had the most impact on the simulation, namely, f, k 2 , and Q 2 (see the first term in Equation A3). Parameters controlling mortality, turnover, adaptation, dormancy, and substrate limitation had a smaller effect on the dynamics ( Figure S9). This increases our confidence in our parameter choices because there are more data available on maximum respiration and growth rates than the other processes.
As discussed in section 2.2.3, hydrolysis has been shown to be the rate-limiting step in organic matter decomposition (Mata-Alvarez et al., 2000), and we choose to model it as purely a chemical process based on the evidence from Veeken and Hamelers (1999) who conclude that enzymes are not limiting in this process. They showed that hydrolysis has an Arrhenius-type behavior with an activation energy of around 0.66 eV (which is equivalent to a Q 10 of approximately 2.8; see Equation A4); thus, the parameter Q 1 (Equation A2) is set to 2.8. We also changed the dependence on soil carbon in Equation A2 to be C 2/3 since the hydrolyzing bacteria are only able to access the surfaces of the organic matter (a 2/3 power is generically true for 2 vs. 3 dimensions and does not make an assumption about shape), and the C dependence in JULES (Equation 2) was therefore also changed to C 2/3 for consistency for the simulations in this study. This means that k 1 in our study has different units from previous studies using JULES. For k 1 and τ (Equations A2 and A8), we use the values that were fitted for JULES based on the site-level data (see section 2.3.4 and Table A1).
When substrate is not limiting, ϕðS; TÞ ¼ 1 and substrate consumption increases with temperature according to the Arrhenius function, as we would expect of a biochemical process (see right-hand term in Equation A2). Yvon-Durocher et al. (2014) tell us that the activation Note. These were chosen as the mean of "layer 1" and "layer 2" in Dankers et al. (2011) (Table 2), since the upper few soil layers are the most important for methane emissions.

10.1029/2020GB006678
Global Biogeochemical Cycles energy should be around 0.93 eV at the community level, which is equivalent to a Q 10 value of around 4.3, so Q 2 is set to 4.3 (Table A1). The maximum respiration rate at a given temperature (i.e., a value for k 2 ) can be estimated from Figure 1 in Price and Sowers (2004). We estimate the maximum value at 0°C to be around 0.01 hr −1 . However, a more recent study shows that this parameter in fact consistently acclimates to the mean annual temperature of the environment (Bradford et al., 2019). Observing that the instantaneous exponential temperature response in Bradford et al. (2019) was reduced by a factor of approximately 0.55 after acclimation, we modified k 2 and k d by dividing by an Arrhenius function with a sensitivity of 0.55Q 2 (Equation A6). Time scales for adaptation to any long-term temperature changes were set to approximately 5 years, so that adaptation would only follow changes in the mean annual temperature and not seasonal changes. Exactly what this time scale should be and the mechanisms behind the acclimation will be considered in future work.
The microbes go dormant in poor growth conditions (which are generally considered to be the driver of dormancy; Rittershaus et al., 2013), which are defined as the growth rate falling below a threshold, μ. The threshold μ was estimated from Ström et al. (2012). Since this study measures a time series of substrate concentration (S), temperature, and methane emissions, it is possible to estimate both A 2 ϕ and the increase and decrease of active biomass through the season (since the methane emissions increase relative to A 2 ϕ in the middle of the season and subsequently decline). The chosen value of μ (Table A1) is also consistent with the minimum growth rate at which microbes will grow in incubation studies such as Rivkina et al. (2000), and in line with this, we also do not allow the methanogens to grow when below this threshold. The minimum activity level was chosen to match the maintenance respiration rate which we estimated from Price and Sowers (2004), and we consider the rate of reactivation or dormancy to follow the maximum rate of biomass growth, assuming that some similar cellular-level processes are required. This is doubtless an approximation, but it is also consistent with the reactivation/dormancy rate that we estimated from Ström et al. (2012). When combined with acclimation, the threshold temperature for reactivation (and therefore for growth) will vary with the temperature of the ecosystem, which is consistent with observations, which show a variety of minimum temperature thresholds for growth of methanogens (Nozhevnikova et al., 2001). Minimum growth temperatures can even be as high as 30°C (Zinder et al., 1984). In our simulations, the minimum temperature for growth ranges from approximately -1°C at Samoylov up to +10°C at Twitchell and changes by no more than 1°C over the observed time period.
There is little information as to what determines the mortality and turnover of methanogenic biomass, particularly in field conditions. We approximate the mortality rate by an Arrhenius function similar to the growth rate (Equation A3) and only activate this process when maintenance respiration cannot be balanced by carbon assimilation (i.e., growth)-in other words, in very poor conditions. Otherwise, we assume that turnover is due to maintenance respiration. Pavlostathis et al. (1990) estimated a turnover rate under good growth conditions which was 0.025 times their maximum growth rate, which is consistent with our maintenance respiration parameter α. König et al. (1985), on the other hand, starved their methanogens and saw a mortality rate approximately equal to maximum growth rate; thus, we set the mortality rate (k d ) equal to the growth rate (fk 2 ). While there is a lot of uncertainty here, the model is not sensitive to these parameters (see Figure S9).
Hydrolysis has been shown to continue in frozen soil (Segura et al., 2017) with a reduced rate, hence the reduction of the substrate production term (k 1 ) when the liquid water content falls below the threshold at which it is able to flow between pores (Olesen et al., 2001). We chose a constant factor assuming that the production becomes limited by where microbes are in contact with the organic matter surfaces, which we estimated from a crude estimation of organic matter surface area and microbial population. The sensitivity test showed that this parameter had only a small impact on the simulations ( Figure S9). Since the methanogens' activity level is modeled, they will downregulate their activity in frozen conditions anyway, so we do not apply such a factor to methanogenesis.
Finally, the functional form of ϕ is derived from first principles, as described below. We first consider one substrate molecule moving randomly around one microbe, the probability that it will reach the microbe in a given time period, δt, is proportional to δt (longer time step: more chance to reach it) and the rate at which it is moving, v (faster moving: more chance to reach the microbe), as well as factors such as size of microbe which we consider to be constant. Assuming that substrate molecules move independently of each other, the number of molecules that a microbe receives in δt will be the number of substrate molecules in the vicinity of the microbe multiplied by the probability of each one reaching the microbe, that is, number of molecules per microbe (mpm), ∝ δt × v × S.
If on average every microbe receives one substrate molecule in each time step δt, but they are distributed randomly, then, for example, some microbes will receive two and some will receive none. Estimating this fraction using a random number generator suggests that around 75% of them will receive one or more substrate molecule in this case. Using the random number generator, we can produce a curve of fractions for different "average number of substrate molecules per microbe" (mpm) which asymptotes to 1 (i.e., all microbes receive a substrate molecule) at high concentrations of substrate. This curve is fitted well by the function tanhððmpmÞ 0:8 Þ.
The relevant time period, δt, is the time period in which the microbes will process a substrate molecule, which becomes smaller at higher T as the respiration rate increases, following the inverse of the maximum respiration rate, that is, ∝ 1/A. The rate of diffusion of substrate molecules, v, increases with temperature according to Ouattara et al. (2000) (note that this value is for acetic acid [acetate], a commonly used methanogenic substrate; however, a range of different molecules have similar temperature sensitivities of diffusion in water; Maharajh & Walkley, 1973). Therefore, the final form of ϕ is where ρ is an undetermined parameter. We estimated ρ based on Ström et al. (2012), by comparing the variability in methane emissions against the variability in substrate concentrations at each point in time (where temperature and mean S are known). Note that this form matches observed behavior, in which more substrate-limited systems are less temperature sensitive.
depressions covering the majority of the mire. The flux tower was located in a wet part of the mire next to a shallow lake. Such wet areas are dominated by sedges (Carex spp., Eriophorum angustifolium), often with Salix spp. growing along the margins. Peat thickness in the mire ranges between 0 and 3 m. Mean annual temperature has been increasing and has fluctuated around 0°C since the 1990s (Callaghan et al., 2010). A detailed description of the site can be found in Jammet et al. (2017). Typical profile of soil organic carbon content was as in Chadburn et al. (2017). The snow depth time series from the Abisko research station was used for snowfall correction, which although it was not from the mire led to a realistic soil temperature simulation, so it was considered suitable for this study.

B2 Chersky
The Chersky observation site is situated within the floodplain of the Kolyma River in Northeast Siberia. The terrain is flat, but with a fine-scale mosaic of microsite conditions. The largest fraction of the terrain remains waterlogged for most of the growing season, supporting a wet tussock tundra ecosystem dominated by cotton-grass meadows, tussock-forming carex species, and low shrubs. Soils consist of an organic layer of 15-20 cm overlaying alluvial mineral soils (silty clay), with some organic material present in deeper layers due to cryoturbation. The mean annual temperature is -11°C. For more information, please refer to Göckede et al. (2017), Kittler et al. (2017), and Kwon et al. (2017). Soil carbon density was estimated from separate profiles: bulk density and soil carbon percentage (Corradi et al., 2005;Kwon et al., 2017). Since the values were from different profiles, the depth of the organic layer was different so the carbon density was estimated from the average bulk density above and below the organic layer, applied to carbon percentages (where the end of the organic layer was clear based on a sharp decrease in carbon percentage and respective increase in bulk density).

B3 Kytalyk
The Kytalyk site is located in an oligotrophic tundra area in the Indigirka lowlands, Russia. The area is underlain by continuous ice rich permafrost. Microrelief at the surface consists of high-and low-centered ice wedge polygons and low palsas (icerich mounds). The typical active layer depth during summer ranges from around 10 cm in dry areas to 50 cm in wet areas. The station is situated on a river terrace, where the vegetation consists primarily of Betula nana, Eriophorum angustifolium, Eriophorum vaginatum, Salix pulchra, and Sphagnum sp., with a typical height of 10-50 cm. Typical profile of soil organic carbon content was as in Chadburn et al. (2017). Mean annual air temperature is around -13°C. More details can be found in Parmentier et al. (2011) and van der Molen et al. (2007).

B4 Lompolojänkkä
Lompolojänkkä is an open mesotrophic sedge fen located in northern Finland. A small stream flows through the fen feeding water to the site and keeping the peat profile continuously water saturated. There is a relatively dense vegetation layer dominated by Betula nana, Menyanthes trifoliata, Salix lapponum, and Carex spp. with mean vegetation height of 40 cm. The peat depth is up to 3 m at the center of the fen; an average pH value of 5.5 was measured for the top peat layer. The mean annual temperature of -1.4°C has been measured at the nearest long-term weather station during the period 1971-2000 (Pirinen et al., 2012). For a more detailed description of Lompolojänkkä, see Aurela et al. (2009) and Lohila et al. (2010). Soil carbon and bulk density for upper soil layers was from Metzger et al. (2015) (Table A1).

B5 Samoylov
Samoylov Island is situated in the Lena River delta in Northern Siberia. Samoylov is a polygonal tundra site dominated by the Drepanocladus revolvens-Meesia triquetra-Carex chordorrhiza community in the wet polygon centers and trenches or collapsed ridges and the Hylocomium splendens-Dryas punctata community in the well-drained plateaus, polygon rims, and elevated centers. The mean annual air temperature is -12.3°C. Typical profile of soil organic carbon content was as in Chadburn et al. (2017). Soil temperatures were taken from a polygon center to represent the saturated part of the landscape. For details of the site, see Boike et al. (2019) and references therein.

B6 Seida
Seida is a sedge-dwarf shrub tundra site located in North West Russia. The landscape is dominated by upland tundra with Betula nana, Salix spp., Empetrum nigrum, Ledum palustre, and Vaccinium spp. being the most important plants. A significant portion of the landscape is also composed of peat plateaus with low water table, and thermokarst lakes, willow stands, and Carex fens with high water tables surrounding the peat plateaus. Mean annual temperature in the region is -5.7°C. More detailed description can be found in Biasi et al. (2014) and Marushchak et al. (2013). Typical landscape units soil carbon profiles were available for Seida (Hugelius et al., 2011), and the Carex fen land cover type was used, since this represents the saturated area where the majority of methane emissions originate.

B7 Twitchell
The Twitchell wetland west pond is a 22-year-old restored tule/catail freshwater wetland on Twitchell Island, CA (USA). The site is managed with the aim of maintaining the water table approximately 25 cm above the soil surface at all times. The mean annual air temperature is around 15.5°C. The water is a mix of cold ocean water from an estuary and snow melt from nearby mountains, so it is cooler than may be expected based on the air temperature. Further details can be found in, for example, Miller et al. (2008) and Miller and Fujii (2010). Upper soil layers carbon and bulk density values were taken from Miller et al. (2008) (Table 2).

Appendix C: Evaluation of the JULES-Microbe Model: Does It Give the Right Answers for the Right Reasons?
Modeling microbial dynamics appears to reproduce the seasonal methane emissions at wetland sites well (Figures 7 and S6). However, in order to further evaluate the internal mechanisms of the JULES-microbe model, we compared it against the limited observational studies that are available to evaluate dissolved substrate, biomass, and biomass activity/dormancy.
In Liu et al. (2011), methanogen quantities at three sites in China were measured and related to the soil organic carbon fraction. This relationship is compared against simulated methanogenic biomass carbon by JULES-microbe in Figure C1. Note that this comparison is not exact because converting from cells g-dry-weight −1 to gC m −3 requires assumptions about the soil density (which would vary between sites and may therefore alter the relationship to some extent) and the carbon mass of the cells, neither of which were stated in Liu et al. (2011). We show two alternative conversions in Figure C1, assuming a high, mineral soil density (1 g/cm 3 ), a low, peat soil density (0.1 g/cm 3 ), and a more typical organic soil density of 0.2 g/cm 3 . Given the non-exact nature of this relationship, the overall methanogenic biomass quantities in the model are generally within the uncertainty range of the observations, and the relationship with soil carbon is similar to that observed, so overall, this comparison is reasonable.
We also compared the microbial model to two incubation studies. In the first study (van Hulzen et al., 1999), both the dissolved substrate and methane production were tracked over time, at a range of different temperatures ( Figure S7). The model recreates the observed dynamics of both acetate and methane from that study fairly well, showing an initial increase in substrate availability until the methanogens have activated sufficiently to start consuming it, at which point it reduces back to a smaller, equilibrium value and methane production continues at a constant rate. The second study (Juottonen et al., 2008) sampled soil from a site at four different times of year, incubated it, and measured the methane production ( Figure S8). They found the lowest methane production  (Liu et al., 2011). The relationship in Liu et al. (2011) is based on number of cells rather than biomass carbon and on soil carbon per gram of soil; therefore, there may be a range of soil densities or cell carbon masses. Hence, several lines are shown for different soil densities and the cell biomass carbon is assumed to be 20 fg, which is commonly used and corresponds to a reasonable cell diameter of 0.6 μm (Jabłoński et al., 2015) and a carbon content of 20% (Bratbak &amp;Dundas, 1984). from soils sampled in October and August and highest in February and May, which JULES-microbe recreates, along with producing the same order of magnitude of methane emissions. The model does not recreate the results exactly (in particular, October production from JULES-microbe is too high and February too low). However, since we did not model the site in the study and simply compared the most similar site that was modeled (Lompolojänkkä), the correspondence is unlikely to be exact. The high-temperature dynamics are also not recreated, but we deliberately did not include high-temperature dynamics in JULES-microbe for simplicity, since soils are unlikely to reach +40°C. These caveats aside, the fact that the approximate magnitude of methane production in soils sampled at different times of year can be captured by JULES-microbe suggests that methanogen activity levels and substrate concentrations vary in a realistic way (also noting that the original JULES scheme would not simulate any differences at all).
In situ substrate concentrations were measured at a site in Greenland in Ström et al. (2012), and these compare well with the modeled values. Specifically, their range of substrate concentration values during summer at 10-cm depth is 0.5-3 mgCL −1 and JULES-microbe simulates values in summer of around 2 mgCL −1 at similar temperature sites. Ström et al. (2012) also show a clear increase of methane emissions relative to the temperature and substrate concentration toward the peak of the growing season and a decline later in the season, which JULES-microbe recreates due to an increase and decrease in methanogen activity levels.
Altogether, these studies provide validation that the internal mechanisms of the model structure are realistic.