An assessment of natural methane fluxes simulated by the CLASS-CTEM model

Natural methane emissions from wetlands and fire, and soil uptake of methane, simulated using the Canadian Land Surface Scheme and Canadian Terrestrial Ecosystem (CLASS-CTEM) modelling framework, over the historical 1850–2008 period, are assessed by using a one-box model of atmospheric methane burden. This one-box model also requires anthropogenic emissions and the methane sink in the atmosphere to simulate the historical evolution of global methane burden. For this purpose, global anthropogenic methane emissions for the period 1850–2008 were reconstructed based on the harmonized representative concentration pathway (RCP) and Emission Database for Global Atmospheric Research (EDGAR) data sets. The methane sink in the atmosphere is represented using biascorrected methane lifetimes from the Canadian Middle Atmosphere Model (CMAM). The resulting evolution of atmospheric methane concentration over the historical period compares reasonably well with observation-based estimates (correlation = 0.99, root mean square error= 35 ppb). The modelled natural emissions are also assessed using an inverse procedure where the methane lifetimes required to reproduce the observed year-to-year increase in atmospheric methane burden are calculated based upon the specified global anthropogenic and modelled natural emissions that we have used here. These calculated methane lifetimes over the historical period fall within the uncertainty range of observation-based estimates. The present-day (2000– 2008) values of modelled methane emissions from wetlands (169 Tg CH4 yr−1) and fire (27 Tg CH4 yr−1), methane uptake by soil (29 Tg CH4 yr−1), and the budget terms associated with overall anthropogenic and natural emissions are consistent with estimates reported in a recent global methane budget that is based on top-down approaches constrained by observed atmospheric methane burden. The modelled wetland emissions increase over the historical period in response to both increases in precipitation and in atmospheric CO2 concentration. This increase in wetland emissions over the historical period yields evolution of the atmospheric methane concentration that compares better with observation-based values than the case when wetland emissions are held constant over the historical period.


Introduction
Earth system models (ESMs) represent physical climate system processes and their interactions with biogeochemical processes focusing primarily on the carbon cycle in the context of carbon dioxide (CO 2 ).These models are able to project how the atmospheric concentration of carbon dioxide ([CO 2 ]) will change in response to changes in anthropogenic CO 2 emissions or alternatively diagnose anthropogenic CO 2 emissions required to achieve a specific CO 2 concentration pathway (Jones et al., 2013).This capability is achieved by modelling [CO 2 ] as a prognostic variable, which itself requires modelling of the surface-atmosphere exchange of CO 2 and hence the need for land and oceanic carbon cycle components in ESMs (Friedlingstein et al., 2006(Friedlingstein et al., , 2014;;Arora et al., 2013).While most ESMs include the capability of modelling [CO 2 ] as a prognostic variable, there are only a handful of ESMs that are beginning to treat the atmospheric concen-Published by Copernicus Publications on behalf of the European Geosciences Union.tration of methane, ([CH 4 ]), as a fully prognostic variable (Collins et al., 2011;Shindell et al., 2013).
The [CH 4 ] has increased from 700 ± 25 ppb in 1750 to 1795 ± 18 ppb in 2010 (Prather et al., 2012).The [CO 2 ] has increased globally from 278 (276-280) ppm in 1750 to 390.5 (390.3-390.7)ppm in 2011.The greater global warming potential of CH 4 compared to CO 2 (84 and 28 for 20-year and 100-year time scales, respectively), has made methane the second most radiatively important greenhouse gas (GHG) after CO 2 (Myhre et al., 2013).The CO 2 radiative forcing for the period 1750-2011 is 1.82 W m −2 (associated with ∼ 112 ppm increase), while the radiative forcing for CH 4 over the same period is 0.48 W m −2 (associated with ∼ 1081 ppb = 1.08 ppm increase) (Myhre et al., 2013).As methane is a short-lived GHG with an atmospheric lifetime of around 9 years (compared to CO 2 which has an atmospheric lifetime of around 100-200 years), mitigation of anthropogenic CH 4 emissions can lead to a decrease in its atmospheric concentration within a timeframe of 10-20 years.As a result, methane is considered a short-lived climate forcer (SLCF) and thus reduction in its anthropogenic emissions offers an attractive and potentially viable target for short-term climate change mitigation policies (Shindell et al., 2012).However, to be able to address climate benefits of reduction in anthropogenic CH 4 emissions within the framework of comprehensive ESMs it is necessary to model [CH 4 ] as a prognostic variable in these models.
Treatment of [CH 4 ] as a fully prognostic variable in ESMs is hindered by at least two factors.First, the global CH 4 budget is not as well understood as for CO 2 .Our lack of ability to close the present-day global CH 4 budget is illustrated in Saunois et al. (2016) who present a recent synthesis of several studies and summarize the present-day global CH 4 budget.Saunois et al. (2016) show a large discrepancy between total CH 4 emissions, from both anthropogenic and natural sources, for the 2003-2012 period, as inferred from the top-down atmospheric inversion-based approaches (558 Tg CH 4 yr −1 ) and those based on bottom-up modelling and other approaches (736 Tg CH 4 yr −1 ).The primary reason for this discrepancy is that there are multiple sources of both natural and anthropogenic CH 4 emissions, so the bottom-up approaches that add up all the individual sources inevitably give larger total emissions than top-down approaches that are constrained by the atmospheric CH 4 burden and its loss in the atmosphere.Second, unlike CO 2 , CH 4 has a sink in the atmosphere, which requires representation of atmospheric chemistry in ESMs to properly account for the removal of CH 4 and feedbacks of methane on chemistry.CH 4 is destroyed in the troposphere and stratosphere due to its reaction with OH radicals and chlorine.This is typically very computationally expensive to represent.As an example, the model years per wall clock day simulated by the atmospheric component of the second generation Canadian Earth System Model (CanESM2; Arora et al., 2011) are reduced by a factor of around 6 when atmospheric chemistry is turned on.Despite these two challenges there are ways forward to model [CH 4 ] as a fully prognostic variable and be able to use comprehensive ESMs to ask questions that the climate modelling community has asked so far in the context of CO 2 .For example, how would future [CH 4 ] change in response to changes in anthropogenic and natural CH 4 emissions, or alternatively what should anthropogenic future CH 4 emissions be to achieve a given CH 4 concentration pathway, while anthropogenic CO 2 emissions continue to increase?In terms of emissions, as the top-down estimates of CH 4 emissions from natural and anthropogenic sources are better constrained than the bottom-up estimates they are likely to provide more robust estimates for evaluating ESMs and their CH 4 related components.The expensive atmospheric chemistry modules can be replaced with simple first-order representations of chemical losses or, ignoring the spatial variations in CH 4 concentration, the global average concentration of methane can be simulated with a box model using specified methane life times which are calculated a priori using full 3-D chemistry-climate models.Although, of course, using specified CH 4 losses implies that feedbacks of methane on methane loss rates and interactions between atmospheric chemistry and climate can be neglected.
The CLASS-CTEM modelling framework serves as the land surface component in the family of Canadian ESMs (CanESMs) (Arora et al., 2009(Arora et al., , 2011;;Arora and Scinocca, 2016) developed by the Department of Environment and Climate Change, Government of Canada, and models the landatmosphere fluxes of water, energy, and CO 2 .It consists of the Canadian Land Surface Scheme (CLASS) and the Canadian Terrestrial Ecosystem Model (CTEM).In preparation for modelling [CH 4 ] as a prognostic variable in future versions of CanESMs we have included several CH 4 related processes in the CLASS-CTEM modelling framework.These include representations of dynamic natural wetlands and their CH 4 emissions, CH 4 emissions from fires, and uptake of CH 4 by soils.This paper evaluates the simulated spatial distribution of wetlands as well as the magnitude of CH 4 emissions from wetlands and fires, and CH 4 uptake by soils against their respective present-day observationbased estimates.We also evaluate the simulated time evolution of the global sums of these fluxes for the 1850-2008 period by using a one-box model of atmospheric CH 4 burden.This one-box model requires anthropogenic CH 4 emissions, emissions from other natural sources that are not modelled in the CLASS-CTEM framework, and a representation of atmospheric sinks.The anthropogenic CH 4 emissions for the period 1850-2008 are obtained by harmonizing the RCP and EDGAR data sets, and natural emissions from sources that are not modelled are specified.Finally, the atmospheric sink of CH 4 is based on bias-corrected global atmospheric lifetime of CH 4 as computed by the Canadian Middle Atmosphere Model (CMAM).The one-box model of atmospheric CH 4 burden is used to evaluate CLASS-CTEM simulated natural CH 4 fluxes by comparing simulated evolution of global [CH 4 ] with their observation-based estimates as well as by comparing the CH 4 lifetime required to reproduce the observed evolution of global [CH 4 ] over the historical period with their observation-based estimates.
The rest of this paper is organized as follows.A brief description of the CLASS-CTEM modelling framework is presented in Sect. 2 along with the details of methane related processes that are implemented, the data sets used and the experimental protocol.Results are presented in Sect. 3 and finally, discussion and conclusions are presented in Sect. 4.  et al., 2009), CanESM2 (Arora et al., 2011), andCanESM4.2 (Arora andScinocca, 2016).CLASS simulates atmosphere-land fluxes of energy and water and it prognostically calculates the liquid and frozen soil moisture contents, and soil temperature for its soil layers, the liquid and frozen moisture contents and temperature of the single vegetation canopy layer (if present), and the snow water equivalent and temperature of a single snow layer (if present).CLASS is described in detail in Verseghy (1991), Verseghy et al. (1993), andVerseghy (2000).In the version 3.6 of CLASS used here, the thicknesses of the three permeable soil layers are specified as 0.1, 0.25 and 3.75 m, although the model can be configured to use any number of layers with specified thicknesses.The thicknesses of the permeable layers also depend on the depth to the bedrock which is specified on the basis of the global data set of Zobler (1986).For example, if the depth to bedrock is only 2 m, then the thicknesses of the permeable soil layers are taken to be 0.1, 0.25 and 1.65 m.The energy and water balance calculations are performed for four plant functional types (PFTs) (needleleaf trees, broadleaf trees, crops and grasses).CLASS operates at a sub-daily time step and a time step of 30 min is used here.
CTEM simulates the fluxes of CO 2 at the land-atmosphere boundary and in doing so models vegetation as a dynamic component of the climate system.It models photosynthesis, autotrophic respiratory fluxes from its three living vegetation components (leaves, stem, and roots, denoted by L, S, and R, respectively) and heterotrophic respiratory fluxes from its two dead carbon components (litter and soil carbon, denoted by D and H , respectively).The flow of carbon through these five carbon pools is explicitly tracked, which allows the calculation of the amount of carbon in these pools as prognostic variables.Disturbance through fire and land use change are also modelled.CTEM cannot operate without coupling to CLASS.Its photosynthesis module operates at the same time step as CLASS and requires estimates of net radiation and soil moisture from CLASS.In return, CTEM provides CLASS with dynamically simulated structural attributes of vegetation which are functions of the driving meteorological data.The amount of carbon in the leaves, stem and root components is used to estimate structural attributes of vegetation.The leaf area index (LAI) is calculated from leaf biomass using PFT-dependent specific leaf area (SLA, m 2 (Kg C) −1 ), which determines the area of leaves that can be constructed per unit leaf carbon biomass (Arora and Boer, 2005a); vegetation height is calculated based on stem biomass for tree PFTs and LAI for grass PFTs (Arora and Boer, 2005a); and rooting depth is calculated based on root biomass (Arora and Boer, 2003).Other than photosynthesis, all terrestrial ecosystem processes in CTEM are modelled at a daily time step.CTEM models its terrestrial ecosystem processes for nine PFTs that map directly to the PFTs used by CLASS.Needleleaf trees in CTEM are divided into deciduous and evergreen for which terrestrial ecosystem processes are modelled separately, broadleaf trees are divided into cold and drought deciduous and evergreen types, and crops and grasses are divided into C 3 and C 4 versions based on their photosynthetic pathways.Version 2.0 of the model is explained in detail in Melton and Arora (2016), while version 2.1 is used here which amongst other minor changes includes all methane related processes discussed below.
The methane related processes implemented in CLASS-CTEM build on the model's existing capabilities.Processes are implemented to be able to dynamically model the geographical distribution of wetlands and their methane emissions, methane emissions from fire and methane uptake by upland soils.The fractional coverage of wetlands in a grid cell is based on flat fraction within a grid cell with slope less than 0.2 % and grid-averaged soil moisture.The methodology is explained in Appendix A.

Wetland methane emissions
The dominant controls on methane emissions in nature are considered to be (1) the position of the water table below which methane is produced due to anoxic decomposition of available organic matter and above which methane is oxidized, (2) the soil temperature which determines the rate of decomposition of organic matter, (3) the availability of organic matter itself, and (4) the pathway through which methane is transferred to the atmosphere (through soil via molecular diffusion, through stems of the vascular plants and through ebullition if the water table is above the soil surface).These factors are not completely independent and their relative importance changes as environmental conditions change www.biogeosciences.net/15/4683/2018/Biogeosciences, 15, 4683-4709, 2018 (Walter and Heimann, 2000).The explicit consideration of these factors becomes more important as the spatial scale at which CH 4 emissions are being modelled reduces.For example, X. Zhu et al. (2014) show that as the modelling spatial scale reduces from 100 to 5 km, the dominant control on simulated wetland CH 4 emissions switches from soil temperature to water table depth.At the current operational resolution of around 2.81 • of CanESM (equal to about 310 km at the equator) we expect dominant controls on methane emissions to be soil moisture, soil temperature, and the availability of organic matter and this allows us to use the simple approach that we have used here.This approach also allows us to estimate CH 4 emissions from wetlands without the use of wetland specific PFTs, as at large spatial scales net primary productivity (NPP) and heterotrophic respiration are governed by climate (Chengjin et al., 2016).CH 4 emissions are simulated to occur over the wetland fraction of the grid cell.The simulated CH 4 emissions from wetlands are calculated by scaling the heterotrophic respiratory flux (R h ) from model's litter (D) and soil (H ) carbon pools, which itself depends on soil temperature, soil moisture, and the available organic matter.Heterotrophic respiration from the litter and soil carbon pools takes the following basic form, with the formulation explained in detail in Melton and Arora (2016).
where ς i represents the base respiration rate (kg C (kg C) −1 yr −1 ) at 15 • C, C i is the amount of carbon in model's litter or soil carbon pool (kg C m −2 ), f i (Q 10 ) = Q 0.1(T i −15) 10 is a Q 10 function that models the effect of temperature, T i is the temperature of litter or soil carbon pool ( • C), and f i ( ) is the function that reduces heterotrophic respiration when soils are too dry and too wet using the soil matric potential ( ) (Melton et al., 2015).The constant 2.64 × 10 −6 converts units from kg C m −2 yr −1 to mol CO 2 m −2 s −1 .
Modelled CH 4 emissions from wetlands, per unit area of a grid cell (mol CH 4 m −2 s −1 ), are calculated as where f w is the wetland fraction in a grid cell (see Appendix A), α w is the ratio of wetland to upland heterotrophic respiratory flux, and δ s converts flux from CO 2 to CH 4 units but also takes into account that some of the CH 4 flux is oxidized in the soil column before reaching the atmosphere.A value of 0.45 is used for α w , as heterotrophic CO 2 respiratory flux over lowlands is typically lower than over uplands, due to limitation by increased soil moisture including a high water table level.While the f i ( ) function in Eq. ( 1) does reduce heterotrophic respiration when soils are wet it does so using only the grid averaged soil moisture content.Wania et al. (2010) use a preferred value of δ s equal to 0.1 and Q. Zhu et al. (2014) found δ s varies between 0.1 and 0.7 with a mean value of 0.23 when calibrating their model against data from 19 sites.A value of 0.135 is used for δ s in this study.The product α w δ s thus equals 0.061, which implies that for each mol CO 2 m −2 s −1 of heterotrophic respiratory flux 0.061 mol CH 4 m −2 s −1 is generated over unit area that is deemed wetland.At large spatial scales CH 4 and CO 2 heterotrophic respiratory fluxes are expected to be highly correlated, as to the first order they are both governed by temperature and the amount of organic matter available for decomposition (Dalva et al., 2001;Q. Zhu et al., 2014).In addition, as the spatial scale increases it is possible to ignore the effect of water table depth as X.Zhu et al. (2014) illustrate.

Fire methane emissions
Fire in CLASS-CTEM is modelled using an intermediate complexity scheme, which represents both natural and human-caused fires, and accounts for all elements of the fire triangle: fuel load, combustibility of fuel, and availability of ignition sources.The fire module accounts for both natural fires caused by lightning and anthropogenic fires which are the result of ignitions caused by humans expressed as a function of population density.Increasing population density increases human-caused fire ignitions but also increases suppression of fire.The suppression of fire represents fire-fighting efforts, landscape fragmentation, and other processes which leads to a reduction in area burned and is also modelled as a function of population density.The original fire parametrization is described in Arora and Boer (2005b), which has since been adapted and used in several other DGVMs (Kloster et al., 2010;Li et al., 2012;Migliavacca et al., 2013).The fire module in CTEM v. 2.1 incorporates changes suggested in these studies as well as several new improvements which are summarized in detail in Melton and Arora (2016).The approach has been evaluated at the global scale in Arora and Melton (2018) who assess how reduction in global wildfire emissions since the 1930s leads to an enhanced land carbon sink.The two primary outputs from the fire module are fraction of area burned per grid cell and dry organic biomass burned per unit area (gC m −2 ).The dry organic matter burned is then multiplied by corresponding emissions factors to obtained emissions (g species m −2 ) for several species of trace gases and aerosols including methane (CO 2 , CO, CH 4 , H 2 , NHMC, NO x , N 2 O, total particulate matter, particulate matter less than 2.5 µm in diameter, and black and organic carbon).These emissions factors are based on an updated set by Andreae and Merlet (2001) listed in Tables 3 and 4 of Li et al. (2012).

Soil uptake of methane
The methane uptake over soil occurs over the unsaturated (upland) fraction of a grid cell that is not deemed wetland.The parameterization is based on an exact solu-tion of the one-dimensional diffusion-reaction equation in the near-surface (top) soil layer and described in detail in Curry (2007).Briefly, the methane uptake by soil is a function of diffusion of methane into soil (which depends on atmospheric methane concentration) and its subsequent oxidation by microbes.The diffusion of methane into the soil depends primarily on air filled porosity of the soil and increases as the pore volume filled by liquid and frozen moisture decreases.The oxidation of methane by microbes is a function of both soil moisture and temperature.Oxidation preferably occurs when soils are neither too dry (when microbial activity is limited by low soil moisture) nor too wet (when microbes are deprived of oxygen).Warmer temperatures favour oxidation of methane in soil and oxidation increases by about four times as soil temperature increases from 0 to 27.5 • C. Finally, the inhibition of methane uptake in cultivated soils is accounted for by a linear factor that reduces oxidation as crop fraction in a grid cell increases.

Forcing data for the CLASS-CTEM model
The  personal communication, 2012).The meteorological variables (surface temperature, pressure, precipitation, wind, specific humidity, and incident short-wave and long-wave radiation fluxes) are available at a spatial resolution of 0.5 • × 0.5 • and at a 6-hourly time interval for the period 1901-2015.These data are regridded to a spatial resolution of 2.81 • and temporally to a 30 min time step to drive the CLASS-CTEM model.Temperature, pressure, wind, specific humidity, and long-wave radiation are linearly interpolated in time while short-wave radiation is assumed to change with the solar zenith angle with maximum radiation occurring at solar noon.Following Arora (1997) the 6-hourly precipitation amount (P , mm (6 h) −1 ) is used to estimate the number of wet half-hours in a given 6-hour period and the 6-hourly precipitation amount is randomly distributed over these wet half-hours.Figure B1 in Appendix B shows the annual landaveraged temperature and precipitation (excluding Antarctica) as derived from the CRU-NCEP data.Both temperature and precipitation show an overall increase over the 20th century that continues into the 21st century, associated with the changing climate.
The land cover data are used by the model to specify the fractional coverage of CTEM's nine PFTs in each grid cell.These data are based on a geographical reconstruction of the historical land cover driven by the increase in crop area (Arora and Boer, 2010) but using the crop area data based on the LUH2 v1h version of the Hurtt et al. (2006) land cover product.The final data set consists of the fractional coverage of CTEM's nine PFTs for the period 1850-2015 at the global scale and at 2.81 • spatial resolution.The increase in crop area over the historical period leads to decrease in area of natural vegetation thus leading to deforestation.A fraction of deforested vegetation is burned but deforested biomass is also converted to paper and wood products which decompose over time leading to land use change emissions.These processes are described in detail in Arora and Boer (2010).In context of terrestrial methane budget, an increase in crop area leads to lower methane uptake by soil over the cultivated fraction of a grid cell.
The globally averaged atmospheric CO 2 and CH 4 concentrations used to drive the model are obtained from the data sets put together for the sixth phase of the Coupled Model Intercomparison Project (CMIP6) and available from input4MIPs web site (https://esgf-node.llnl.gov/projects/input4mips/, last access: April 2017).These data are shown in Fig. B2.

Observation and model-based data for CLASS-CTEM evaluation
In addition to evaluating the CLASS-CTEM simulated methane emissions from wetlands, fire and methane uptake by soils in the context of the one-box atmospheric CH 4 model as mentioned in Sect.2.2 below, we also evaluate simulated present-day wetland extent and all modelled methane fluxes directly against other model and observation-based estimates.
The CLASS-CTEM simulated wetland extent is compared against two data sets: the wetland data from the Global Lakes and Wetlands Database (GLWD; Lehner and Döll, 2004) and a new product that is formed by merging remote sensing based observations of daily surface inundation from the Surface Water Microwave Product Series (SWAMPS;Schroeder et al., 2015) with the static inventory of wetland area from the GLWD.The derivation of the second product is explained in detail in Poulter et al. (2017).SWAMPS provides estimates of fractional surface water based on data from multiple passive and active microwave satellite missions.While open water (e.g.rivers, lakes, and ocean) and inundated wetlands comprising of open plant canopies are mapped by satellites, inundation beneath closed forest canopies, and exposed wetlands with water table below the surface cannot be mapped.However, satellite data are able to provide the seasonal cycle which static data sets like GLWD cannot.The merged SWAMPS-GLWD product attempts to overcome limitations of both individual data sets.
The simulated present-day methane emissions from wetlands and fire, and methane uptake by soil, are compared to top-down estimates compiled by Saunois et al. (2016).We also compare the anthropogenic emissions we have used within the framework of one-box atmospheric methane model (Sect.2.2) with estimates from Saunois et al. (2016).
Finally, we also evaluate the model regionally over the West Siberian lowlands (WSL).This region is chosen because inversion-based methane fluxes are readily accessible over the region, which were compiled and documented for the WETCHIMP-WSL intercomparison project (Bohn et al., 2015).We compare simulated wetland extent and wetland methane emissions with observation-and inversionbased results from Bousquet et al. (2011), Kim et al. (2011), andWinderlich (2012) and participating models in the Wetland and Wetland CH 4 Intercomparison of Models Project (WETCHIMP, Melton et al., 2013) focused on the West Siberian lowlands region (WETCHIMP-WSL) (Bohn et al., 2015).Of these the Kim et al. (2011) and Winderlich (2012) are regional inversions.Kim et al. (2011) used wetland methane emissions from Glagolev et al. (2010) at 1 • spatial resolution as their prior and used the NIES-TM atmospheric transport model for the period 2002-2007.They derived climatological monthly wetland emissions optimized to match atmospheric methane concentrations obtained by aircraft sampling.Winderlich (2012) used the Kaplan (2002) wetland inventory for prior wetland emissions with the TM3-STILT global inversion system for year 2009.Their posterior monthly wetland emissions were uniquely determined for each grid cell within their domain at 1 • spatial resolution and optimized to match atmospheric methane concentrations measured at four tower observation sites located between 58 and 63 • N. The Bousquet et al. (2011) is a global inversion but uses two priors, the first based on the Matthews and Fung (1987) emissions inventory and the second based on Kaplan (2002).Bousquet et al. (2011) inversion used the Laboratoire de Météorologie Dynamique general circulation model (LMDZ) atmospheric transport model at a 3.75 • × 2.5 • grid and estimated monthly methane emissions at a 1 • spatial resolution for the period 1993-2009.Being a global inversion, they optimized atmospheric concentrations relative to global surface observations at several flask stations for methane but also other trace gases.For wetland extent, we compare the CLASS-CTEM simulated wetland extent over the WSL region with models participating in the WETCHIMP-WSL intercomparison, and GLWD and SWAMPS + GLWD products mentioned above but also the Global Inundation Extent from Multi-Satellites (GIEMS; Prigent et al., 2007;Papa et al., 2010) derived from visible and near-infrared and active and passive microwave sensors for the period 1993-2004.In addition, we also use the estimate from Peregon et al. (2009) who used a regional wetland typology map further refined by satellite image classifications to calculate the wetland extent in the WSL region.

One-box model of atmospheric methane
A one-box model of atmospheric CH 4 is used to evaluate the time evolution of simulated methane emissions from wetlands (E w ), methane emissions from fire (E f ) and the soil uptake of methane (S soil ) over the period 1850-2008.The model describes the changes in burden of atmospheric CH 4 (B) as a balance of surface emissions (consisting of natural, E N , and anthropogenic emissions, E A ) and the atmospheric (S atmos ) and surface soil sinks (S soil ).
where t is the time and Eq. ( 3) is applied at an annual time step.The atmospheric CH 4 burden (B, Tg CH 4 ) equals 2.78 times [CH 4 ] (represented in units of parts per billion, ppb) (Denman et al., 2007).The distinction between natural and anthropogenic emissions is not straightforward for fire, which contains emissions due to both natural and humancaused fires.3) can be rewritten as Biogeosciences, 15, 4683-4709, 2018 www.biogeosciences.net/15/4683/2018/where t = 1 year.Equation ( 4) can be used to evaluate simulated natural methane emissions E N in two ways.First, when all the terms on the right hand side of Eq. ( 4), including an initial value of B(t), are known, then the time evolution of B can be calculated and compared to its observation-based estimate.Second, if the time evolution of B is specified based on observations of methane concentration in the atmosphere, then the value of τ chem required to satisfy Eq. ( 4) can be calculated (see Eq. 5) and compared its observation-based estimate e.g. from Prather et al. (2012).
In Sect.3, we have used both these methodologies to assess CLASS-CTEM simulated methane emissions from wetlands (E w ), methane emissions from fire (E f ), and the soil uptake of methane (S soil ).Note that Eq. ( 4) does not include any term for oceanic methane emissions.Saunois et al. (2016) report a range of 0-5 Tg CH 4 yr −1 , with a mean value of 2 Tg CH 4 yr −1 for oceanic methane emissions.Given the large uncertainty in other components of the global methane budget, and the small magnitude of oceanic methane emissions, we have ignored this term.Equation (3) may also be used to calculate pre-industrial methane concentration assuming atmospheric methane concentration was in equilibrium with pre-industrial emissions and sinks, i.e. by setting dB dt = 0, which yields Substituting S atmos (t) = B (t) 1 − exp (−1/τ chem (t)) similar to as was done in Eq. ( 4) and solving for B (t) gives Using values of sources and sinks and τ chem corresponding to year 1850, for example, yields B (t) (Tg CH 4 ), which when divided by 2.78 yields [CH 4 ] in ppb for 1850.

Anthropogenic methane emissions
The time evolution of global E A (used in Eq. 4) for the 1850-2008 period is based on two data sets.The first data set is the decadal representative concentration pathway (RCP) anthropogenic methane emissions data set (version 2.0.5)available at 0.5  1 of Saunois et al., 2016).These two data sets are blended (or harmonized) to obtain a consistent time series of annual global anthropogenic methane emissions for the period 1850-2008.First the RCP and EDGAR data are regridded from their 0.5 and 0.1 • spatial resolutions, respectively, to the 2.81 • resolution at which the model is applied.Next, all non-biomass burning emission categories are added separately in both data sets to obtain total anthropogenic emissions for each data set.The emissions categories in both data sets are somewhat different as shown in Table 1.Emissions from fire and biomass burning are excluded because CLASS-CTEM simulates CH 4 emissions from fire explicitly.As our framework requires only total anthropogenic methane emissions (excluding biomass burning) the different emissions categories in the two data sets do not matter.The wetland emissions from rice paddies, which we do not explicitly model, are included in the specified anthropogenic emissions (in category 11 of the EDGAR data set and category 1 of the RCP data set as shown in Table 1).Equation ( 8) summarizes this harmonization methodology for a given grid cell.Fugitive from solid 5. Industrial processes and combustion 6.
Oil production and refineries 6.
Gas production and distribution 7.
Waste treatment and disposal 8.
Industrial  and annual EDGAR emissions (excluding biomass burning).The RCP and EDGAR emissions are fairly similar for the period 1970-1990 but are different after 1990.Figure 1b and c illustrate how the harmonization works for two selected grid cells based on Eq. ( 8).In Fig. 1b anthropogenic methane emissions are shown for a land grid cell where emissions first increase and then decrease.In Fig. 1c anthropogenic methane emissions are shown for an ocean grid cell, with six orders of magnitude lower emissions than the land grid cell, where emissions more or less continuously increase.In both case, the harmonization ensures that by 1970 the adjusted RCP emissions are same as the EDGAR emissions.Although for the purpose of using E A in Eq. ( 4) only its global values are required, the methodology described here yields a continuous consistent gridded data set of anthropogenic methane emissions for the period of our analysis with no abrupt jumps.

Lifetime of atmospheric methane
To use Eq. ( 4), for evaluation of CLASS-CTEM simulated annual values of E w , E f and S soil over the historical period, time-evolving annual values of τ chem are required.We obtain values of τ chem simulated by the Canadian Middle Atmosphere Model (CMAM).The CMAM is a fully interactive chemistry-climate model (CCM) that is based on a vertically extended version of the third generation Canadian Atmospheric General Circulation Model with a model lid at 95 km (Scinocca et al., 2008).The model contains a description of the important physical and chemical processes of the stratosphere and mesosphere and has been extensively assessed against observations through participation in two phases of the Chemistry-Climate Model Validation (CCMVal) activity (Eyring et al., 2006;SPARC CCMVal, 2010).Of more importance for methane, the chemistry has been extended throughout the troposphere by including cloud corrections on clear-sky photolysis rates, emissions of ozone precursors CO and NO x (NO + NO 2 ) including emissions of NO x from lightning, hydrolysis on specified tropospheric sulfate aerosols and interactive wet and dry deposition.Note that the chemical mechanism currently used in CMAM has not yet been extended to include the chemistry of non-methane hydrocarbons important in the troposphere; only methane chemistry is considered.
Results from CMAM with tropospheric chemistry have been submitted to the Atmospheric Chemistry and Climate Model Intercomparison Project (ACCMIP) (Lamarque et al., 2013).The experimental design for ACCMIP involved timeslice simulations at various points in time between 1850 and 2100, with simulations for year 2000 conditions used for assessing the model chemical climate against available presentday observations.The ACCMIP intercomparison found that the CMAM produced estimates of tropospheric chemical quantities that fell well within the range of current generation CCMs.For example, the present-day tropospheric ozone burden from CMAM was 323 Tg, vs. a multi-model mean of 337 ± 23 Tg, where the range is given as one standard deviation across the 15 participating models (Young et al., 2013).The present-day methane lifetime to reaction with OH in the troposphere was found to be 9.4 years, again well within the range of ACCMIP models of 9.7 ± 1.5 years (Naik et al., 2013).However, like the majority of ACCMIP models, the CMAM predicts a too fast removal of CH 4 by OH as compared with our best-estimate from methyl-chloroform decay of 11.2 ± 1.3 years (Prather et al., 2012).As described further below, the calculated CH 4 lifetime to chemical loss from CMAM is thus scaled to agree with observationally based estimates before being used in the one-box atmospheric model of CH 4 .
Time-dependent values of τ chem are derived from a simulation over the 1850-2014 period that uses specified seasurface temperatures and sea-ice from one member of the five-member historical ensemble performed by CanESM2 for the CMIP5.Data for 2006-2014 was taken from a continuation of that simulation for the RCP 6.0 scenario.Specified anthropogenic and biomass burning emissions of CO and NO x were taken from the CMIP5 historical database (Lamarque et al.The observation-based estimate of τ chem is obtained from Prather et al. (2012) who calculate τ CH 4 based on Eq. ( 9) as where τ OH (present day value of 11.2 years), τ strat (120 years), τ trop−Cl (200 years), and τ soil (150 years) are the lifetimes associated with the destruction of CH 4 by tropospheric OH radicals, loss in the stratosphere, reaction with tropospheric chlorine and uptake by soils, respectively, which yields a present day (corresponding to year 2010) value of τ CH 4 as 9.1 ± 0.9 years.τ chem is the methane life time associated with all chemical processes in the atmosphere.For the pre-industrial period (corresponding to year 1750), Prather et al. ( 2012) estimate τ CH 4 as 9.5 ± 1.3 years assuming τ OH to be equal to 11.76 years (based on Atmospheric Chemistry and Climate Model Intercomparison Project (ACCMIP) results; Voulgarakis et al., 2013) and lifetimes associated with other processes are assumed to stay the same.In our study, the methane soil sink is explicitly simulated in the CLASS-CTEM framework and this corresponds to the term 1 τ soil in Eq. ( 9).The remaining terms in Eq. ( 9) all correspond to atmospheric processes.Removing the 1 τ soil term in Eq. ( 9) gives us an observation-based estimate of τ chem making it consistent with CMAM's methane lifetime corresponding only to atmospheric processes and increases the observation-based estimates of pre-industrial and present-day atmospheric CH 4 lifetimes to 10.15 ± 1.45 years and 9.7 ± 1.0 years, respectively.
Figure 2 compares the τ chem values from CMAM with its observation-based estimate from Prather et al. (2012) and shows that CMAM based estimates are biased low and outside the uncertainty range of observation-based estimates, as mentioned earlier.When used in the one-box model of atmospheric methane, the lower than observed τ chem values will inevitably lead to a higher than observed atmospheric sink (S atmos ) and thus lower than observed atmospheric methane concentration even if all the other flux terms (E N , E A , and S soil ) in Eq. ( 5) are realistic.Therefore, we adjust the CMAM derived values of τ chem upwards by 15 % so that they lie within the uncertainty range of observation-based estimates of the τ chem as derived by Prather et al. (2012).These adjusted values of τ chem are also shown in Fig. 2.

Equilibrium pre-industrial simulation
The equilibrium pre-industrial simulation was initialized from zero biomass for all PFTs.The fractional coverages of CTEM's nine PFTs for the pre-industrial simulation are based on the land cover product described in Sect.2.1.2for 1850.The model was then driven with 1901-1925 CRU-NCEP climate data cycled repeatedly until the model pools reach equilibrium.The early part of the 20th century does not show any significant trends compared to the later part of the 20th century, as seen in Fig. B1 (Appendix B), so using the 1901-1925 data to spin up the model to equilibrium for 1850 conditions is reasonable.Atmospheric CO 2 and CH 4 concentration levels were set to 285 ppm and 791 ppb, respectively, corresponding to their pre-industrial 1850 levels.This pre-industrial equilibrium simulation yields initial conditions for all CLASS-CTEM prognostic variables for the transient 1851-2015 simulation.

Transient historical simulation
The transient historical simulation is performed for the period 1851-2015 and its prognostic variables are initialized from the equilibrium pre-industrial simulation as mentioned above.In this transient simulation, (1) wetland extent and its methane emissions respond to changes in climate but also increases in atmospheric CO 2 concentration (which increases net primary productivity and thus heterotrophic respiration), (2) methane emissions from fire respond to changes in climate, atmospheric CO 2 concentration, population density, and land use change, and (3) methane uptake by soil responds to climate, changes in atmospheric CH 4 concentration and changes in crop fraction.

Results
We first show the results from the transient 1851-2015 simulation and evaluate the time-evolution of CLASS-CTEM simulated global natural methane fluxes over the historical period using the one-box model of atmospheric methane described in Sect.2.2.1.This is followed by evaluation of model fluxes for the present day against observation-based estimates and for the WSL region using observation-based estimates and results from other models.

Time evolution of simulated global natural methane fluxes
Figure 3 shows the time evolution of simulated annual maximum wetland extent, methane emissions from wetlands and fire, and soil uptake by methane from the 1851-2015 transient historical simulation.In Fig. 3 (1) the simulated wetland methane emissions increase by 30 % over the historical period from about 130 to 169 Tg CH 4 yr −1 from 1850s to the present day (2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008), partly due to increase in wetland extent which increases by 8 % from 7.5 to 8.1 million km 2 from 1850s to the present day, (2) the simulated fire methane emissions decrease from their 1850s value by 20 % from about 34 to 27 Tg CH 4 yr −1 for the present day, and (3) the soil methane uptake more than doubles from its 1850s value of about 14 to 29 Tg CH 4 yr −1 for the present day.
The increase in wetland methane emissions over the historical period is due to an increase in the wetland area, driven by increase in precipitation seen in Fig. B1 (Appendix B), but also higher methane fluxes per unit area.The higher methane fluxes per unit area are caused by increase in atmospheric CO 2 concentration which increases both net primary productivity and heterotrophic respiration over the historical period as shown in Fig. C1 in Appendix C. As wetland methane emissions are proportional to heterotrophic respiration (R h ) in Eq. ( 2), an increase in R h also increases methane emissions from wetlands.
The decrease in methane emissions from wildfires is driven by a decrease in area burned, which itself is driven by an increase in crop area and population density over the historical period (Arora and Melton, 2018).The increase in cropland area decreases area burned by wildfires in the model, as croplands are not allowed to burn.In the real world cropland area also fragments the landscape, which affects the spread of fire.Direct anthropogenic influences on wildfires are more complex as accidental, as well as intentional, human-caused ignitions enhance wildfires, while anthropogenic suppression of wildfires decreases area burned and fire related emissions.Figure C2 (Appendix C) shows that the overall effect of increase in crop area and population density in the model is that area burned increases slightly up to about 1930, and then starts decreasing thereafter and this area burned pattern compares reasonably with the decadal charcoal index from the Global Charcoal Database version 3 (Marlon et al., 2008) for the full length of the historical simulation.Wind and smoke carry charcoal from fires and deposit it onto aquatic sediments and this forms the basis of sediment charcoal indices.The caveat with the comparison with sediment charcoal records is that they only provide a proxy for fire activity and indicate if fire activity is higher or lower relative to a point in time.Figure C2 also compares area burned with estimates from version 4.1s of the Global Fire Emissions Database (Randerson et al., 2012;Giglio et al., 2013) that are based on the satellite record and available only for the short 1997-2014 period.Model and observation-based average burned area over this 1997-2014 period are 483.4 and 485.5 million ha yr −1 and their trends are −5.57± 1.25 and −3.43 ± 1.05 million ha yr −2 , respectively.The negative trends indicate burned area has been decreasing.
Finally, the increase in methane uptake by soils is primarily the response to an increase in the atmospheric concentration of methane.Diffusion of methane into the soil is directly proportional to its atmospheric concentration (Curry, 2007) which more than doubles from around 790 ppb in 1850 to around 1830 ppb in 2015 (Fig. B2b).

Evaluation of simulated global natural methane fluxes
We first evaluate the CLASS-CTEM simulated global natural methane fluxes in a forward simulation where all the right hand side terms of Eqs. ( 3) and ( 4) are specified at an annual time step and the change in atmospheric methane burden dB dt is calculated every year.Although the CRU-NCEP meteorological data are available to 2015 allowing us to perform offline CLASS-CTEM simulations up until 2015, the last year for which harmonized RCP-EDGAR emissions are available is 2008.Therefore, we simulate the time evolution of atmospheric methane concentration for the period 1851-2008.In this forward calculation of atmospheric methane burden the one box model may be initialized using the observed 1850 methane concentration or using the 1850 concentration that is in equilibrium with 1850 sinks and sources.The latter is calculated by assuming dB dt = 0 as illustrated in Eqs. ( 6) and ( 7).The numerator term in Eq. ( 7) for 1850 is 214.2Tg CH 4 yr −1 (wetland emissions are 130.1, other natural emissions are specified at 25, anthropogenic emissions are 39.5, fire emissions are 33.2, and soil uptake is 13.6).In Eq. ( 7), using an adjusted CMAM τ chem value of 8.82 years (for 1850, red line in Fig. 2) yields an 1850 equilibrium [CH 4 ] for 1850 of 719 ppb compared to its observation-based value of 791 ppb.In contrast, using the midrange value of τ chem of 9.97 years based on Prather et al. (2012) yields equilibrium [CH 4 ] for 1850 of 808 ppm which compares much better with its observation-based estimate of 791 ppm.An equilibrium methane concentration in 1850 that is lower than observed, when using adjusted CMAM τ chem , implies that the atmospheric methane sink is higher and associated with lower than observed atmospheric www.biogeosciences.net/15/4683/2018/Biogeosciences, 15, 4683-4709, 2018 methane lifetime in 1850.Indeed, the adjusted CMAM τ chem for 1850 of 8.82 years is less than Prather et al. (2012) midrange estimate of 9.97 years but still within their uncertainty range as seen in Fig. 2. Regardless, as the lifetime of methane in the atmosphere is only around 10 years, the effect of initial conditions disappears by around 1870.This is shown in Fig. 4a, which compares the simulated evolution of atmospheric methane burden with its observations using the adjusted CMAM τ chem .Figure 4a shows that overall the simulated increase in methane concentration over the historical period is reasonable compared to observation-based estimates despite the various specified and modelled sources and sinks that contribute to the time evolution of atmo-spheric methane burden.In Fig. 4a An alternative approach to evaluate the modelled natural sinks and sources is to specify the rate of increase of atmospheric methane burden according to its observations and calculate the required atmospheric lifetime of methane (excluding the soil sink), given modelled natural sinks and sources, following Eq.( 5) as discussed in Sect.2.2.1.These results are shown in Fig. 4b, which compares the methane lifetimes required to achieve the observed increase in methane concentration over the historical period (black line) to within ±3 ppb (shaded area between grey lines) with observation-based estimate of atmospheric methane lifetime based on Prather et al. (2012) and the adjusted atmospheric methane lifetime from CMAM (both of which were shown earlier in Fig. 2). Figure 4b shows that for the most part, the calculated atmospheric methane lifetime stays within the uncertainty of observation-based estimates.Moreover the temporal trend after 1900 in calculated atmospheric methane lifetime compares well to the trend in the atmospheric methane lifetime from the CMAM model.Both anthropogenic emissions and methane concentration during the early part of the 1850-2008 historical period are more uncertain than during the later part and thus the differences between simulated and observation-based estimates of atmospheric methane concentration (in Fig. 4a) and methane lifetimes (in Fig. 4b) for the period 1850-1900 are not unexpected.
While the results in Fig. 4a and b provide some confidence that the magnitude and temporal evolution of simulated global natural methane sources and sinks over the historical period are reasonable they, of course, do not allow the evaluation of all of the simulated natural fluxes individually.
We also evaluate the role of increase in wetland methane emissions over the historical period (as seen in Fig. 3) on the historical methane budget.Instead of using wetland methane emissions from the transient historical simulation in the onebox model of atmospheric methane we use wetland methane emissions from the equilibrium pre-industrial simulation in which 1901-1925 CRU-NCEP meteorological data are used repeatedly and CO 2 is held constant at its pre-industrial level of 285 ppm.As a result wetland extent and methane emissions do not respond to changes in climate and increasing CO 2 , and do not increase over the historical period (as seen in Fig. 5a).Methane emissions from fire and soil uptake of methane still respond to changes in climate and increasing CO 2 .The result of using these wetland methane emissions (shown in Fig. 5a) in the framework of the one-box model of atmospheric methane is shown in Fig. 5b.In Fig. 5b, although [CH 4 ] increases overall over the historical period in response to increase in anthropogenic emissions, the result of wetland methane emissions not increasing over the historical period is that the simulated methane concentration in year 2008 is calculated to be 1667 ppm, which is 130 ppb lower than the 1797 ppb seen in Fig. 4a.  the Global Lakes and Wetland (GLWD; Lehner and Döll, 2004) and the new product formed by merging remote sensing based observations of daily surface inundation from the Surface Water Microwave Product Series (SWAMPS; Schroeder et al., 2015) with the static inventory of wetland area from the GLWD from Poulter et al. (2017), as mentioned earlier in Sect.2.1.3.Maximum wetland fraction from the model and SWAMPS + GLWD product is calculated as the maximum of 12 mean monthly values from the 13 years spanning the 2000-2012 period.Figure 6 shows that overall the model is able to capture the broad latitudinal distribution of wetlands with higher wetland fraction at northern highlatitudes and in the tropics.The model yields higher wetland fraction in the tropics than both observation-based estimates and this is due to higher wetland fraction simulated in the Amazonian region.The Amazonian region is densely forested and the SWAMPS product is unable to map wetlands beneath closed forest canopies.Biases also likely exist in the GLWD data set as parts of the Amazonian region are fairly remote.This is shown in Fig. 7 which compares the geographical distribution of simulated annual maximum wetland fraction with that from the GLWD and SWAMPS + GLWD products.The model successfully captures wetlands in the Hudson Bay lowlands, the West Siberian lowlands, the Pantanal and the region bordering Argentina, Paraguay and Uruguay in South America, Indonesia, and the low lying region around Bangladesh.In terms of differences from these observation-based data sets the model most notably overestimates wetland extent in Europe.One possible reason for this is that wetlands in Europe have been drained for agriculwww.biogeosciences.net/15/4683/2018/Biogeosciences, 15, 4683-4709, 2018 ture and our wetland parameterization does not take this into account.About two-thirds of the European wetlands that existed 100 years ago have been lost (European Commission, 1995) leading to a substantial decrease in the number and size of large bogs and marshes, and small or shallow lakes.Methane emissions from wetlands (168.9Tg CH 4 yr −1 ) are the largest of natural fluxes, as is well known, while emissions from fire (26.8 Tg CH 4 yr −1 ) and methane uptake by soil (28.7 Tg CH 4 yr −1 ) are an order of magnitude lower.As expected, the geographical distribution of methane emissions from wetlands (Fig. 8a) corresponds well to the geographical distribution of wetlands themselves (Fig. 7a) although per unit wetland area methane emissions are higher in tropics and milder temperate regions than in high-latitude regions.This is because warmer temperatures and a longer growing season in the tropical and milder temperate regions imply that wetlands can emit more methane per unit wetland area and for a longer period of time than the colder high-latitude regions with a shorter growing season.In Fig. 8b the geographical distribution of methane fire emissions shows higher values in seasonally dry tropical regions and order of magnitudes lower values in mid-high latitude regions.These results are consistent with area burned (not shown), which shows a sim-ilar pattern (Arora and Melton, 2018).Finally, the geographical distribution of methane uptake by soils in Fig. 8c shows higher methane uptake by soils in parts of arid regions (including the Sahara and the Australian outback) where soil moisture does not get too dry (so as to not excessively limit soil microbial activity) but otherwise has fairly uniform uptake in the tropics and lower values in mid-high latitude regions where lower temperatures and higher soil moisture limit methane uptake by soils.

Regional evaluation over West Siberian lowlands
While evaluation of simulated global wetland extent and wetland methane emissions, and their geographical distribution, provides confidence in model results, we further evaluate the model at a regional scale over the West Siberia lowlands.The model results are sampled for the region lying between 50 to 75    SWAMPS and GLWD product, while both the satellite-based inundation products by themselves (SWAMPS and GIEMS) show much lower values.Satellite-based products that remotely sense inundated areas can only do so when water table is above the ground and thus wetland areas inferred from these products are expected to be lower in magnitude than products that also take into account land cover, as is the case for the merged SWAMPS and GLWD product.
In Table 2, the annual maximum wetland extent is quite similar for the merged SWAMPS and GLWD product (0.55 million km 2 ) and CLASS-CTEM simulated values (0.53 million km 2 ) but the maximum occurs in different months.In Fig. 9b, for the merged SWAMPS and GLWD product the maximum wetland extent occurs in June while CLASS-CTEM simulated values show peak both in June and September.The participating models from the WETCHIMP-WSL intercomparison show a range of values for the monthly wetland extent in the WSL region (Fig. 9a).Models range from those which specify constant values with no seasonality for the wetland ex-tent to models that dynamically model wetland extent two of which show maximum wetland extent of greater than 1 million km 2 .The average annual maximum wetland extent across the participating models in the WETCHIMP-WSL is 0.70 ± 0.15 million km 2 (mean ± standard error).Finally, the Peregon et al. (2009) estimate is 0.68 million km 2 , which is somewhat higher than the CLASS-CTEM simulated value (0.53 million km 2 ) and the merged SWAMPS and GLWD product (0.55 million km 2 ).
Figure 9b compares CLASS-CTEM simulated monthly wetland methane emissions with those from models participating in the WETCHIMP-WSL intercomparison and four inversion-based estimates mentioned in Sect.2.1.3.The two Bousquet et al. (2011) inversions shown in Fig. 9b correspond to ones using the reference Matthews and Fung (1987) emissions inventory (Bousquet 2001 R) and the emissions inventory based on Kaplan (2002) (Bousquet 2001 K).All values are reported as average for the period 1993-2004 except for the Kim et al. (2011)  Table 2 compares the simulated annual wetland methane emissions from CLASS-CTEM with those from models participating in the WETCHIMP-WSL intercomparison and the four inversion-based estimates.In Table 2, the inversionbased annual wetland methane emissions vary from 3.08 to 9.80 Tg CH 4 yr −1 .The highest annual emissions in the Winderlich (2012) inversion are due to higher emissions in the shoulder months of spring and fall compared to other inversions but also non-zero emissions during winter months (December to February) as seen in Fig. 9b.The CLASS-CTEM model calculates annual wetland methane emissions of 7.76 Tg CH 4 yr −1 and the average for models participating in the WETCHIMP-WSL intercomparison is 5.34 ± 0.54 Tg CH 4 yr −1 (mean ± standard error).Of all the models and inversions only the Winderlich (2012) inversion shows substantial methane emissions for the November to April period.In CLASS-CTEM as the liquid soil moisture in the top soil layer freezes wetland extent contracts to zero (Fig. 9a) and methane emissions are shut off during the winter months.Bohn et al. (2015) note that Winderlich (2012) inversion-based estimates may have been influenced by emissions from fossil fuel extraction and biomass burning, although the seasonality of Winderlich (2012) fluxes, with nonzero emissions even in winter, is plausible.Based on yearround eddy flux measurements of methane emissions from Alaskan Arctic tundra sites, Zona et al. (2016) find that cold season (September to May) emissions account for ≥ 50 % of the annual methane flux, with the highest emissions from non-inundated upland tundra.They find a major fraction of cold season emissions occur during the "zero curtain" period, when subsurface soil temperatures are near 0 • C. Langer et al. (2015) report winter emissions from tundra ponds in Siberia as they are freezing during early winter.They analyzed concentrations of methane in bubbles (trapped in the In Table 3 the total emissions from natural sources in our framework are 199 Tg CH 4 yr −1 which are in the lower part of the range of 194-292 Tg CH 4 yr −1 compiled by Saunois et al. (2016).This is due to lower specified emissions from nonwetland sources.While our modelled emissions from wetlands of 169 Tg CH 4 yr −1 compare well with the Saunois et al. (2016) central estimate of 166 Tg CH 4 yr −1 , our specified emissions from other natural sources (including termites, geological sources and fresh water bodies) of 25 Tg CH 4 yr −1 are near the low end of their range (21-130 Tg CH 4 yr −1 ).In contrast, our anthropogenic emissions of 344 Tg CH 4 yr −1 (which include emissions from fire for consistency with Saunois et al., 2016) are higher than Saunois et al. (2016) central estimate of 319 Tg CH 4 yr −1 and towards the higher end of their range (255-357 Tg CH 4 yr −1 ).This is due to the use of EDGAR emissions which as Saunois et al. (2016) note are towards the higher end of all data sets of anthropogenic emissions.Our modelled fire emissions of 27 Tg CH 4 yr −1 are lower than Saunois et al. (2016) central estimate of 35 Tg CH 4 yr −1 .One reason for this is that while we include natural and anthropogenic fires in our framework we do not account for biofuel burning.Overall, our emissions from natural sources are 35 Tg CH 4 yr −1 lower, and emissions from anthropogenic sources are 25 Tg CH 4 yr −1 higher, than the Saunois et al. (2016) central estimates.As a result, the sum of natural and anthropogenic emissions (543 Tg CH 4 yr −1 ) in our framework is 9 Tg CH 4 yr −1 lower than Saunois et al. (2016) 2016) but also the historical evolution of atmospheric methane burden and methane's lifetime simulated using a one-box model of atmospheric methane.While our simulated and specified present-day global methane budget components lie within the uncertainty range for top-down estimates from Saunois et al. (2016) this also implies that methane emissions from individual sectors for the present-day budget can be increased or decreased within their uncertainty ranges as long as the total emissions stay the same.However, the time evolution of the atmospheric methane burden over the historical period provides additional constraints than the present-day methane budget.For example, our specified methane emissions of 25 Tg CH 4 yr −1 from other nature sources (E o ) (including termites, geological sources, wild animals, and freshwater) are lower than Saunois et al. (2016) central estimate of 68 Tg CH 4 yr −1 , although still within their uncertainty range (21-130 Tg CH 4 yr −1 ).In the absence of any information about its time evolution we have assumed that E o remains constant over the historical period.This is a plausible assumption as we do not expect emissions from termites, geological sources, wild animals, and freshwater to show a large response to changing environmental conditions over the historical period.Certainly, not as large as we saw for wetland emissions which increased by 40 Tg CH 4 yr −1 over the 1850-2008 period (Fig. 3).However, when we use a constant E o of 68 Tg CH 4 yr −1 in our framework we obtain higher than observed methane concentration throughout the historical period and the year 2008 value is 1953 ppb compared to 1797 ppb that we obtain in Fig. 4a (observed methane concentration for 2008 is 1790 ppb).This is shown in Fig. 10.Part of the reason for this may be that present-day EDGAR emissions are higher than other estimates as Saunois et al. (2016) note and that's why we had to choose a lower E o .However, the global harmonized and RCP anthropogenic emissions are fairly similar up until 1990 (see Fig. 1a) and thus had we use RCP based emissions (up until year 2000) we still would have obtained higher than observed atmo- The second constraint provided by the time evolution of the atmospheric methane burden over the historical period is that related to wetland methane emissions.As seen in Sect.3.2, in the absence of the simulated increase in wetland methane emissions from about 130 to 170 Tg CH 4 yr −1 from 1850s to the present day the simulated year 2008 atmospheric methane concentration is about 100 ppb lower than observed (as seen in Fig. 5).Assuming our RCP and EDGAR based harmonized anthropogenic emissions and their increase over the historical period is reasonably realistic, this indicates that it is very likely that wetland methane emissions have indeed increased over the historical period in response to changes in climate and increased atmospheric CO 2 concentration.The results in Sect.3.1 showed that this increase of 30 % in wetland methane emissions in driven more by an increase in methane emissions per unit area than the increase in maximum wetland extent, which increased by about 8 % over the historical period.The implication of this is that wetland methane emissions will likely keep increasing in the future in response to the increasing atmospheric concentration of CO 2 driven by higher heterotrophic respiration.
The evaluation of simulated wetland extent against the GLWD and the merged SWAMPS and GLWD product provides confidence that the model is broadly able to reproduce the geographical distribution of wetlands although, of course, some limitations remain.Over the WSL region the simulated estimates of wetland extent and wetland methane emissions are also broadly consistent with observation-based estimates.
The version of the model used here treats its land mask as binary so each grid cell is either land or water.While the model is capable of representing inland lakes using a separate tile (and simulate the resulting impact on energy and water fluxes) this functionality was not used in this study.In addition, representation of inland lakes requires modelling methane emissions from their anoxic sediment (Bastviken et al., 2004) and care needs to be taken to avoid double counting the wetland extent of inland lakes (Thornton Brett et al., 2016).
We have not evaluated the model's wetland methane emissions at the site level against observations.Wetland methane emissions are known to be spatially highly heterogeneous and temporally intermittent (e.g.Godwin et al., 2013) and CLASS-CTEM does not represent physical processes that govern methane emissions at small spatial and temporal scales.Instead the model is designed for operation at large spatial scales (> 100 km) and implementation within an Earth system model and as such only temperature, soil moisture, and substrate availability (through heterotrophic respiration) are taken into account.Water table depth, ebullition, transport through vascular plants, and PFTs specific to wetlands are not considered in our modelling framework.The corollary of this is that our model cannot be expected to reproduce wetland methane emissions at a point scale where site specific processes and conditions including water depth, ebullition, and wetland specific PFTs become more important.Similar approaches are followed by several large scale models including a recent attempt by Bloom et al. (2017) who derive wetland methane emissions using heterotrophic respiration from eight terrestrial ecosystem models.
We have also not evaluated parameter and forcing uncertainties in relation to simulated methane emissions.However, several aspects of the CLASS-CTEM model have been evaluated before against observations including photosynthesis, autotrophic and heterotrophic respiration, allocation of carbon from leaves to stem and root components, dynamic leaf phenology, fire, and land use change.These aspects of the model have been evaluated at point (Arora and Boer, 2005a;Melton et al., 2015), regional (Peng et al., 2014;Garnaud et al., 2015;Arora et al., 2016) and global (Arora and Boer, 2010;Melton andArora, 2014, 2016) scales.In regards to processes relevant to methane emissions from wetlands (which, of course, is the largest natural source) heterotrophic respiration and wetland extent are the most important.Uncertainties in these two quantities will propagate to calculated methane emissions from wetlands.The majority of in-crease in wetland methane emissions over the historical period comes from an increase in heterotrophic respiratory flux due to increase in gross and net primary productivities in response to increase in atmospheric CO 2 concentration.The rate of increase of simulated gross primary productivity is adjusted in CLASS-CTEM to obtain realistic land carbon sink from 1960 onwards (Le Quéré et al., 2018) but also to obtain a realistic amplitude of the annual CO 2 cycle in a fully coupled Earth system model simulation (Arora and Scinocca, 2016).The response of model's heterotrophic respiration to soil moisture is expressed as a function of soil matric potential as mentioned earlier in Sect.2.1.1.This response is based on Griffin (1981) who suggests that the microbial activity is optimal at an absolute soil matric potential of 0.05 MPa and decreases as the soil becomes waterlogged near 0.00 MPa or too dry near 1.5 MPa .This parameterization has been evaluated indirectly at seasonally dry locations in the Amazonia (Melton et al., 2015).
Our next step to evaluate natural methane fluxes from CLASS-CTEM is to use these fluxes in an atmospheric transport model to simulate and compare methane concentrations at selected stations to assess seasonality of simulated wetland methane emissions at large spatial scales in a somewhat more direct manner.In addition, CLASS-CTEM simulated natural fluxes can be used as a prior in a methane inversion-based system, together with anthropogenic methane emissions, to calculate optimized posterior fluxes to which the prior fluxes can be compared.Although atmospheric inversions-based systems have their own limitations (Houweling et al., 2017), the objective is to evaluate CLASS-CTEM simulated natural methane fluxes using a range of available methodologies.
Overall the results presented here suggest that the natural fluxes of methane between the atmosphere and the land, and the geographical distribution of wetland extent, simulated by the CLASS-CTEM modelling framework are sufficiently realistic to use the model to study the changes in natural methane fluxes due to changes in environmental conditions.

Figure 1 .
Figure 1.Comparison of RCP, EDGAR and their harmonized annual global anthropogenic methane emissions (a) excluding biomass burning.(b, c) illustrate the harmonization technique for a land and an ocean grid cell, respectively.

Figure 2 .
Figure 2. Comparison of atmospheric methane lifetime obtained from the Canadian Middle Atmosphere Model (CMAM) for the period 1850-2010 with observation-based estimates from Prather et al. (2012) but excluding the soil sink as explained in Sect.2.2.3.
, 2010) up to 2000 and for the RCP 6.0 scenario to 2014.Specified concentrations of long-lived greenhouse gases were fromMeinshausen et al. (2011), following RCP 6.0 from 2006-2014.Specified stratospheric aerosol surface area density fields, used to account for the effects of large volcanic eruptions, was based on the 1960-2010 database created for the Chemistry Climate Model Initiative (CCMI) extended back to 1850 following the approach described inNeely III et al. (2016).Solar variability was included by calculating the wavelength-resolved daily variability relative to the long-term average (June 1976-January 2007) from the recommended CMIP6 database(Matthes et al., 2017).
For the period 1851 to 1900 of this simulation the model is driven with meteorological data from 1901-1925 twice, similar to what was done for spinning up the model for the pre-industrial simulation.For the period 1901-2015 the meteorological data corresponding to each year are used.Time varying concentrations of atmospheric CO 2 and CH 4 are specified for the period 1851-2015.The annual timevarying fractional coverages of C 3 and C 4 crop PFTs in each grid cell are based on LUH2 v1h version of the Hurtt et al. (2006) land cover product.

Figure 3 .
Figure 3.Time evolution of simulated natural methane fluxes (shown on the primary y axis) and annual maximum wetland extent (shown on the secondary y axis) by CLASS-CTEM for the 1851-2015 period in the transient historical simulation.

Figure 4 .
Figure 4. Comparison of simulated methane concentration over the historical period with its observation-based estimates (a).The simulation may be initialized from the 1850 observed methane concentration (solid green line) or from an 1850 concentration that is in equilibrium with 1850 specified methane sources and sinks (dashed blue line) as explained in Sect.3.2.Panel (b) compares the methane lifetimes required to achieve the observed increase in methane concentration over the historical period (black line) to within ±3 ppb (shaded area between grey lines) with observation-based estimates of atmospheric methane lifetime based on Prather et al. (2012) and the adjusted atmospheric methane lifetime from CMAM.

Figure 5 .
Figure 5.Time evolution of simulated natural methane fluxes (shown on the primary y axis) and annual maximum wetland extent (shown on the secondary y axis) by CLASS-CTEM for the 1851-2008 period for the case when wetland extent and methane emissions are not allowed to respond to changing climate and increase atmospheric CO 2 over the historical period (a).Panel (b) shows the simulated methane concentration over the historical period, together with its observation-based values, when the natural fluxes shown in panel (a) are used within the framework of the one-box model of atmospheric methane.

Figure 6
Figure 6 compares the zonally averaged maximum wetland fraction over land with observation-based estimates based on

Figure 6 .
Figure 6.Comparison of simulated zonally averaged maximum wetland fraction over land with observation-based estimates based on the Global Lakes and Wetland (GLWD; Lehner and Döll, 2004) and a new product that is formed by merging remote sensing based observations of daily surface inundation from the Surface Water Microwave Product Series (SWAMPS; Schroeder et al., 2015) with the static inventory of wetland area from the GLWD as explained in Poulter et al. (2017).

Figure 8
Figure8shows the geographical distribution of methane emissions from dynamic wetlands (panel a) and fire (panel b) and the soil uptake of methane (panel c) simulated by the CLASS-CTEM model.The figures also show the global total of these fluxes averaged over the 2000-2008 period for later comparison with estimates fromSaunois et al. (2016).Methane emissions from wetlands (168.9Tg CH 4 yr −1 ) are the largest of natural fluxes, as is well known, while emissions from fire (26.8 Tg CH 4 yr −1 ) and methane uptake by soil (28.7 Tg CH 4 yr −1 ) are an order of magnitude lower.As expected, the geographical distribution of methane emissions from wetlands (Fig.8a) corresponds well to the geographical distribution of wetlands themselves (Fig.7a) although per unit wetland area methane emissions are higher in tropics and milder temperate regions than in high-latitude regions.This is because warmer temperatures and a longer growing season in the tropical and milder temperate regions imply that wetlands can emit more methane per unit wetland area and for a longer period of time than the colder high-latitude regions with a shorter growing season.In Fig.8bthe geographical distribution of methane fire emissions shows higher values in seasonally dry tropical regions and order of magnitudes lower values in mid-high latitude regions.These results are consistent with area burned (not shown), which shows a sim-

Figure 7 .
Figure 7.Comparison of geographical distribution simulated annual maximum wetland fraction with observation-based estimates based on the Global Lakes and Wetland (GLWD; Lehner and Döll, 2004) and a new product that is formed by merging remote sensing based observations of daily surface inundation from the Surface Water Microwave Product Series (SWAMPS; Schroeder et al., 2015) with the static inventory of wetland area from the GLWD as explained in Poulter et al. (2017).

Figure 8 .
Figure 8. Geographical distribution of annual emissions from wetlands (a) and fire (b), and the soil sink (c) simulated by the CLASS-CTEM model.The data are averaged over the 2000-2009 period.

Figure 9 .
Figure 9.Comparison of CLASS-CTEM simulated wetland extent (a) and wetland methane emissions (b) over the West Siberian lowlands (WSL) region with those from models participating in the WETCHIMP-WSL intercomparison and observation-and inversion-based estimates as discussed in Sect.3.5.
inversion which reports fluxes for year 2005 and the Winderlich (2012) inversion which corresponds to year 2009.
central estimate (552 Tg CH 4 yr −1 ).The total sink strength is calculated to be 538 Tg CH 4 yr −1 in our framework which compares well with the Saunois et al. (2016) estimate of 546 Tg CH 4 yr −1 .Saunois et al. (2016) do not provide uncertainty ranges for the atmospheric and total sink.The modelled atmospheric (509 Tg CH 4 yr −1 ) and soil (29 Tg CH 4 yr −1 ) sinks in Table 3 also compare well with Saunois et al. (2016) estimates of 514 and 32 Tg CH 4 yr −1 , respectively.4 Discussion and conclusions The offline evaluation of natural methane fluxes simulated by the CLASS-CTEM modelling framework presented here is the first step in making atmospheric methane concentration a prognostic variable in the family of Canadian earth system models.The evaluation is based on comparison of presentday fluxes with existing observation-and model-based estimates compiled by Saunois et al. (

Figure 10 .
Figure 10.Simulated methane concentration over the historical period, together with its observation-based values, when other non-wetland natural methane emissions (E o ) are specified at 68 Tg CH 4 yr −1 over the historical period.
Prather et al. (2012)ogenic emissions excluding those from fire and biomass burning (E A excl fire ) data are explained in Sect.2.2.2.The atmospheric sink S is calculated as a first-order loss process from methane's lifetime τ chem in the atmosphere as S atmos (t) = B (t) 1 − exp (−1/τ chem (t)) .An estimate of τ chem is obtained from the Canadian Middle Atmosphere Model (CMAM) with chemistry and compared to an observation-based estimate fromPrather et al. (2012)as later shown in Sect.2.2.3.With S atmos represented in terms of τ chem Eq. ( Saunois et al., 2016)unois et al. (2016)global CH 4 budget (as shown later in Sect.3) we consider all emissions from fire as anthropogenic, although the CLASS-CTEM model calculates fire emissions due both to lightning and human-caused ignitions.Natural emissions (E N = E w + E o ) consist of modelled wetlands emissions (E w ) and emissions from other natural sources (E o ) (including termites, geological sources, wild animals, and freshwater), which we specify at 25 Tg CH 4 yr −1 (consistent with, but towards, the lower end of the range natural emissions, 21-130 Tg CH 4 yr −1 , as deduced by top-down approaches summarized inSaunois et al., 2016).The reason for specifying E o at 25 Tg CH 4 yr −1 is discussed later in Sect. 4. Anthropogenic emissions (E A = E A excl fire + E f ) consist of specified emissions from all anthropogenic sources excluding fire (E A excl fire ) and fire emissions, which we explicitly model (E f ).
data set is part of version 4.2 of the Emission Database for Global Atmospheric Research (EDGAR, http://edgar.jrc.ec.europa.eu/overview.php?v=42, last access: October 2014) available at 0.1 • spatial resolution and available for the period 1970-2008.These data sets were selected because the RCP data set provides the anthropogenic methane emissions going back to 1850 and the EDGAR data set provides the anthropogenic methane emissions for more recent years since 1970.The EDGAR data set was also chosen because amongst the recent anthropogenic data sets this is the only data set which provides gridded anthropogenic methane emissions at an annual time scale (see Table web-apps/tnt/RcpDb, last access: October 2014).The value for the first year of each decade is assumed to correspond to the rest of the decade.So the 1850 value corresponds to the 1850-1859 decade, the 1860 value corresponds to the 1860-1869 decade, and so on.The second anthropogenic methane emissions

Table 1 .
, where E A, RCP and E A, EDGAR represent annual anthropogenic methane emissions (excluding biomass burning) in the RCP and EDGAR data sets, respectively, E A, RCP adjusted represent the adjusted RCP emissions and t is the time in years from 1850 to 1970.The harmonization algorithm adjusts the annual total anthropogenic methane emissions in the RCP data, for each 2.81 • grid cell, from 1850 to 1970 such that by the time RCP emissions reach 1970 they are the same as the EDGAR's total emissions excluding biomass burning.As a result of this harmonization, the largest change is made to RCP emissions for 1970 and the smallest change is made for year 1851.The RCP emissions for 1850 are not changed.Emissions categories for EDGAR and RCP anthropogenic methane emissions.
The final harmonized time series for E A is obtained by concatenating E A, RCP adjusted for the period 1850 to 1970 and E A, EDGAR for the period from 1971 to 2008.Figure1ashows the harmonized time series of global anthropogenic methane emissions along with the decadal RCP www.biogeosciences.net/15/4683/2018/Biogeosciences, 15, 4683-4709, 2018 , the coefficient of correlation between observation-based (black line) and simulated (green and blue lines)[CH 4] is 0.99 and root mean square error (RMSE) is 35 ppb (green line) and 41 ppb (blue dashed line).The simulated values are somewhat overestimated from 1885 to 1980 and underestimated from 1980 to 2005.

Table 2 .
Bohn et al. (2015)S-CTEM simulated annual methane emissions and annual maximum wetland extent for the West Siberia lowlands (WSL) region with models participating in the WETCHIMP-WSL intercomparison and observation-and inversion-based estimates as discussed in Sect.3.5.Numbers shown are mean ± standard error fromBohn et al. (2015)for models participating in the WETCHIMP-WSL intercomparison.Standard error is not available for all inversions.All values are reported as average for the period 1993-2004 unless otherwise noted.

Table 3 .
Comparison of the components of the present day methane budget based on this study with those from Saunois et al. (2016) (based on synthesis of published studies).The values used in this study are averaged for the period 2000-2008 (as the last year of the version of EDGAR emissions used is 2008) while Saunois et al. (2016) values correspond to the 2000-2009 period.
CTEM values correspond to the 2000-2008 period since the EDGAR anthropogenic emissions were available only until 2008 at the time of this study and thus the one-box model of atmospheric methane is also run up until 2008.For clarity, Table 3 also identifies which fluxes are modelled by CLASS-CTEM, which are specified and which are based on atmospheric methane lifetimes.