WETMETH 1.0: a new wetland methane model for implementation in Earth system models

Wetlands are the single largest natural source of methane (CH₄), a powerful greenhouse gas affecting the global climate. In turn, wetland CH₄ emissions are sensitive to changes in climate conditions such as temperature and precipitation shifts. However, biogeochemical processes regulating wetland CH₄ emissions (namely microbial production and oxidation of CH₄) are not routinely included in fully coupled Earth system models that simulate feedbacks between the physical climate, the carbon cycle, and other biogeochemical cycles. This paper introduces a process-based wetland CH₄ model (WETMETH) developed for implementation in Earth system models and currently embedded in an Earth system model of intermediate complexity. Here, we (i) describe the wetland CH₄ model, (ii) evaluate the model performance against available datasets and estimates from the literature, and (iii) analyze the model sensitivity to perturbations of poorly constrained parameters. Historical simulations show that WETMETH is capable of reproducing mean annual emissions consistent with present-day estimates across spatial scales. For the 2008–2017 decade, the model simulates global mean wetland emissions of 158.6 Tg CH₄ yr⁻¹, of which 33.1 Tg CH₄ yr⁻¹ is from wetlands north of 45^∘ N. WETMETH is highly sensitive to parameters for the microbial oxidation of CH₄, which is the least constrained process in the literature.

1 Introduction

Wetlands are vegetated locations that are inundated with water on a permanent, seasonal, or recurrent basis (Wheeler, 1999). In the context of this study, wetlands are defined following the latest global methane (CH₄) budget report (Saunois et al., 2020): natural ecosystems with inundated or water-saturated soils where anoxic conditions lead to the production of CH₄. Wetlands across the globe are the single largest natural source of atmospheric CH₄, accounting for approximately one-third of total global emissions (Bridgham et al., 2013; Saunois et al., 2016). Estimates of global wetland CH₄ emissions over the past few decades vary between 140 and 210 Tg CH₄ yr⁻¹ (Kirschke et al., 2013). Although different types of wetlands exist, such as bogs, fens, swamps, marshes, and floodplains (Aselmann and Crutzen, 1989; Saunois et al., 2016), the release of CH₄ from any wetland results from the balance between two biogeochemical processes (Segers, 1998): the production of CH₄ by anaerobic microbes (namely methanogens) and the oxidation of CH₄ primarily by aerobic microbes (namely methanotrophs).

Both CH₄ production and oxidation in wetlands are sensitive to changes in climate conditions. For instance, soil warming accelerates the microbial activity with a higher response for methanogenic than methanotrophic activity (Bridgham et al., 2013; Dunfield et al., 1993; Segers, 1998). At the landscape or larger scale, increased wet conditions tend to enhance methanogenic activity to the detriment of methanotrophic activity (Duval and Radu, 2018; Helbig et al., 2017; Kim, 2015). In turn, wetland CH₄ emissions can affect the global climate through changes in atmospheric CH₄ levels and associated radiative forcing (Dean et al., 2018; O'Connor et al., 2010). Analyses of ice cores suggest that CH₄ emissions from tropical and northern wetlands contributed significantly to climate changes during past glacial–interglacial transitions (Loulergue et al., 2008; Rhodes et al., 2017).

The interactions between climate conditions and wetland CH₄ emissions translate into a positive feedback loop that has the potential to amplify changes in global mean surface air temperature, which is a major concern for future climates (Dean et al., 2018; O'Connor et al., 2010). Research on feedbacks between the physical climate and biogeochemical cycles is generally conducted with three-dimensional (3-D) fully coupled Earth system models (ESMs) (Arora et al., 2013). Over the past decade, these ESMs have proven very useful to investigate and inform international climate policies, such as the accounting of carbon emissions required to avoid the risk of dangerous climate change (Zickfeld et al., 2009) and achieve the goals of the Paris Agreement (Tokarska and Gillett, 2018). However, biogeochemical processes regulating CH₄ emissions in wetlands are not commonly included in fully coupled ESM simulations.

In the past, several process-based models have been developed to investigate the response of wetland CH₄ emissions to climate variability and climate change (Hodson et al., 2011; Hopcroft et al., 2011; Pandey et al., 2017; Paudel et al., 2016; Shindell et al., 2004; Zhang et al., 2018; Zhu et al., 2015). These wetland CH₄ models are generally embedded in terrestrial or land surface models and forced with observational datasets or reanalysis products (Melton et al., 2013; Wania et al., 2013; Xu et al., 2016). A second application for wetland CH₄ models has been to quantify the climate response to wetland CH₄ emissions (Gedney et al., 2004, 2019; Z. Zhang et al., 2017). In this case, results from wetland CH₄ models are used in climate–carbon cycle model emulators to assess their impact on radiative forcing (Gedney et al., 2019; Z. Zhang et al., 2017). These modelling studies have contributed to advance research on the possible evolution of wetland CH₄ emissions in the 21st century (Koven et al., 2011; Shindell et al., 2004), the magnitude of their impact on the global climate (Gedney et al., 2019; Z. Zhang et al., 2017), and their implications for international climate policy (Comyn-Platt et al., 2018). However, their quasi-coupling methods do not reflect the complete feedback loop between climate conditions and wetland CH₄ emissions as expected in the natural world. So far, only 1-D and 2-D models of the northern high-latitude regions have been applied for simulating the feedback between climate conditions (temperature changes) and wetland CH₄ emissions in a fully coupled mode (Schneider von Deimling et al., 2012, 2015).

The implementation of process-based wetland CH₄ models in fully coupled ESMs is needed in order to advance research on wetland CH₄–climate feedbacks in the context of global climate projections (Dean et al., 2018). In particular, this addition to Earth system modelling should be beneficial to ongoing research on the permafrost carbon feedback (Nzotungicimpaye and Zickfeld, 2017; Schuur et al., 2015) and the remaining carbon budget for achieving the goals of the Paris Agreement (Rogelj et al., 2019).

This paper introduces a wetland CH₄ model developed for implementation in ESMs and currently embedded in an Earth system model of intermediate complexity (EMIC). Our study aims at developing a computationally efficient process-based model for simulating large-scale wetland CH₄ emissions constrained with sparse observations. Section 2 gives an overview of processes regulating CH₄ emissions in wetlands. Section 3 provides the model description and an outline of performed model simulations. Section 4 describes the model calibration and choice of parameter values. Section 5 presents the model performance evaluation. Section 6 describes the model sensitivity to poorly constrained parameters. Sections 7 and 8 are for discussions and conclusions, respectively.

2 Overview of processes regulating methane emissions in wetlands

2.1 Microbial production of methane

Wetlands host several communities of microbes adapted to the predominant anoxic conditions of these environments (Bridgham et al., 2013). Some of these microbes are methanogens, which decompose organic matter for their metabolism and produce CH₄ as a by-product of their respiration (McCalley et al., 2014; Segers, 1998). The organic matter decomposed by methanogens in wetlands originates from litter fall, root exudates, dead plants, and dissolved organic carbon (Bridgham et al., 2013; Conrad, 2009; Girkin et al., 2018; Mitsch and Mander, 2018). In the northern permafrost region, carbon from thawed soils constitutes an additional source of organic matter to methanogens (Kwon et al., 2019; Olefeldt et al., 2013).

There are three pathways through which methanogens produce CH₄ from soil organic matter (Le Mer and Roger, 2001; Segers, 1998; Whalen, 2005). The first pathway (acetotrophic methanogenesis) is operated by methanogens that rely on acetate for their metabolism, resulting in the production of both CH₄ and carbon dioxide (CO₂) (Bridgham et al., 2013; Whalen, 2005). The second pathway (hydrogenotrophic methanogenesis) is operated by methanogens that produce CH₄ through CO₂ reduction in the presence of hydrogen (Bridgham et al., 2013). The third pathway (methylotrophic methanogenesis) is operated by methanogens that use methylated substrates (e.g. methanol, methylamines, and dimethylsulfide) for their metabolism (Zalman et al., 2018).

Rates of CH₄ production in wetlands are generally highest in upper anoxic layers due to several factors such as the quality of organic matter and the spread of active microbial populations. For instance, in comparison to soil layers at depth, where organic matter can be recalcitrant to microbial decomposition, the organic matter in near-surface soil layers is more labile due to fresh inputs from litter fall and vegetation mortality (Treat et al., 2015; Walz et al., 2017; Wild et al., 2016). Furthermore, observations at various sites show that methanogenic activity decreases as depth increases (Bridgham et al., 2013; Cadillo-Quiroz et al., 2006).

Increasing soil temperatures stimulate the dynamics and growth of methanogenic communities in wetlands, resulting in an increase in CH₄ production rates (Bridgham et al., 2013; Segers, 1998). However, several studies indicate that there is an optimal temperature for methanogenic activity between 25 and 30^∘ C (Dean et al., 2018; Dunfield et al., 1993). Other factors promoting the occurrence of CH₄ production in wetlands include the persistence of anoxic conditions as well as soil pH varying between acidic and neutral (Dunfield et al., 1993; Segers, 1998).

2.2 Microbial oxidation of methane

In wetlands, methanotrophs (CH₄-oxidizing microbes) populate oxic portions of the soil column (Bridgham et al., 2013; Conrad, 2009; Whalen, 2005). Such oxic portions are primarily soil layers close to the surface which are in contact with the atmosphere, commonly near and above the water table (Bridgham et al., 2013; Le Mer and Roger, 2001; Segers, 1998). In the presence of vascular plants, other oxic portions of the soil column can be found near the roots due to the downward transport of oxygen (O₂) through plant aerenchyma (Kwon et al., 2019; Whalen, 2005). All these oxic portions of the soil column constitute the so-called oxic zone, which is predominantly made of soil layers near and above the water table (Bridgham et al., 2013; Conrad, 2009; Segers, 1998). Methanotrophs consume CH₄ that ascends from the zones of production at depth to the overlying oxic zone for their metabolism, and they primarily produce CO₂ as part of their respiration (Bridgham et al., 2013; Segers, 1998).

While O₂ has been considered for years to be the only electron acceptor involved in the microbial oxidation of CH₄, there is a growing evidence of the occurrence of CH₄ oxidation under anoxic conditions operated by anaerobic microbes that rely on alternate electron acceptors such as nitrate and sulfate (Dean et al., 2018). However, although anaerobic CH₄ oxidation in marine environments has been well established for decades (Hoehler et al., 1994; Reeburgh, 1976), this process remains poorly investigated in wetlands despite its potential importance for the CH₄ cycle (Gauthier et al., 2015; Smemo and Yavitt, 2011).

In analogy to CH₄ production, CH₄ oxidation is influenced by changes in soil temperatures (Bridgham et al., 2013; Segers, 1998). For instance, CH₄ oxidation rates increase during the summer because of intensified microbial activity but also the availability of substantial CH₄ in response to increased soil temperatures (Segers, 1998). However, the temperature response for CH₄ oxidation is generally lower than that for CH₄ production (Bridgham et al., 2013; Dean et al., 2018; Dunfield et al., 1993; Segers, 1998).

2.3 Mechanisms transporting methane to the atmosphere

Various mechanisms exist for transporting CH₄ produced in wetlands to the atmosphere. Three transport mechanisms are well documented in the literature and generally monitored in situ (Bridgham et al., 2013; Whalen, 2005): the diffusion of CH₄, whereby molecules of CH₄ slowly ascend the overlying water column; the ebullition of CH₄, whereby bubbles of CH₄ rapidly ascend towards the soil surface; and the transport of CH₄ through the aerenchyma of vascular plants. However, other transport mechanisms for CH₄ in wetlands have been revealed: the hydrodynamic transport of CH₄ in the form of upwelling caused by temperature gradients primarily at nighttime (Poindexter et al., 2016), and the transport of CH₄ through tree stems (Bridgham et al., 2013; Conrad, 2009; Pangala et al., 2017) whose driving processes are still not well understood (Barba et al., 2019).

Methane oxidation is highly dependent on the predominant transport mechanism for CH₄. The water table position plays a crucial role in affecting what fraction of the produced CH₄ reaches the atmosphere (Blodau, 2002; Moore and Roulet, 1993; Segers, 1998). When the water table is well below the surface, methanotrophs may oxidize all of the diffusing CH₄ before the gas reaches the atmosphere (Segers, 1998). In the presence of vascular plants, a lower fraction of the produced CH₄ is oxidized because these plants allow the gas to bypass the oxic zone where methanotrophs are hosted (Blodau, 2002; Bridgham et al., 2013; Segers, 1998). In the case of ebullition, which often occurs episodically, CH₄ may escape to the atmosphere with reduced opportunities for oxidation (Bridgham et al., 2013; Whalen, 2005). How CH₄ oxidation relates to the transport of CH₄ through tree stems (Barba et al., 2019) or by hydrodynamic processes (Poindexter et al., 2016) is not well established.

2.4 A synopsis of wetland methane dynamics

Figure 1 illustrates vertical profiles of soil organic content, CH₄ concentration, and CH₄ oxidation rates in a soil column with and without inundation at the surface based on principles outlined in the literature (Blodau et al., 2004; Whiticar and Faber, 1985). In general, the water table position determines the maximum depth at which O₂ is available in the soil column (i.e. the oxic–anoxic interface). When the surface is flooded and the water is stagnant (Fig. 1a), O₂ diffuses slowly into the soil column and may only be present in a portion of the upper soil layer which is in contact with the atmosphere. Under such predominantly anoxic conditions, CH₄ production occurs throughout the soil column and the concentration of CH₄ mirrors soil organic content – eventually with a small reduction near the surface due to CH₄ oxidation. A modest amount of ascending CH₄ may be oxidized throughout the soil column, although with highest oxidation rates near the surface where some O₂ may be available as an electron acceptor. The combination of high CH₄ production and only modest CH₄ oxidation in the soil column results in large CH₄ emissions into the atmosphere.

Figure 1 Illustrated vertical profiles of soil organic content, CH₄ concentration, and oxidation rates in a soil column with inundation at the surface (a) and without inundation at the surface (b). The vertical profiles are based on principles outlined in the literature (Blodau et al., 2004; Whiticar and Faber, 1985). As an illustrative simplification, the soil organic content is assumed to be identical in (a) and (b). Furthermore, the CH₄ concentration profile is assumed to mirror CH₄ production rates at depth within the soil column (see explanatory text in Sect. 2.4). The blue horizontal line illustrates the water table position, and the dashed red horizontal line illustrates the oxic–anoxic interface or maximum depth at which O₂ is available in the soil column. The relative magnitude of CH₄ flux in the soil column is shown by the upward arrow to the right, also characterizing the relative magnitude of CH₄ emissions into the atmosphere.

When the flooding recedes, O₂ becomes more prevalent in the upper soil column where the CH₄ concentration decreases following a slowdown or shutdown of CH₄ production as aerobic microbes dominate the competition for organic matter (Fig. 1b). CH₄ production persists below the oxic–anoxic interface where the concentration of CH₄ mirrors soil organic content owing to the predominant anoxic conditions. Ascending CH₄ becomes subject to substantial oxidation in the soil column with the highest oxidation rates above the oxic–anoxic interface where O₂ is abundant. The combination of decreased CH₄ production and substantial CH₄ oxidation in the soil column results in small or no CH₄ emissions into the atmosphere.

3 Model description and simulations

3.1 The wetland methane model: WETMETH

Microbial production and oxidation of CH₄ are parameterized in WETMETH using a multilayer ground structure with information on the moisture distribution, the amount of organic matter (carbon content), and the average temperature in each soil layer. These soil variables are commonly simulated by ESMs. Figure 2 provides a schematic representation of WETMETH for a soil column with and without inundation at the surface. By configuration, it is considered that CH₄ emissions in WETMETH may occur not only from inundated locations but also from non-inundated ecosystems with a relatively high level of soil moisture content (Saunois et al., 2016, 2020).

Figure 2 Illustration of the developed wetland CH₄ model (WETMETH) and the dynamics of wetland CH₄ processes as represented in the model. This schematic representation depicts a soil column (model grid box) with inundation at the surface (a) and without inundation at the surface (b). The soil column is shown here with multiple layers of unequal thicknesses. The blue area at the surface of (a) represents the inundated surface area. The blue sections in the different soil layers of (a) and (b) represent water-saturated zones. For both (a) and (b), the dashed red horizontal line illustrates the oxic–anoxic interface, and the orange vertical arrow shows the relative thickness of the oxic zone or oxic zone depth (z_oxic). Larger CH₄ emissions are expected to occur when the soil surface is flooded than when it is not, due to relatively high CH₄ production and moderate CH₄ oxidation in the soil column.

3.1.1 Parameterization of methane production

For any land location, the rate of CH₄ production in an underlying soil layer i (P_i in kg C m⁻³ s⁻¹) is parameterized as follows:

$$\begin{array}{}\text{(1)}& P_i=S\left({\mathit{\theta }}_i\right)C_{i}r{Q}_{\mathrm{10}}^{\frac{T_i-T_{\mathrm{0}}}{\mathrm{10}}}\exp \left(-{\displaystyle \frac{z_i}{{\mathit{\tau }}_{\mathrm{prod}}}}\right),\end{array}$$

where S(θ_i) is the fraction of the soil layer that is saturated with water, and C_i is the amount of soil carbon (in kg C m⁻³) in the layer. Here we consider C_i to be the aggregate of all sources of soil carbon (i.e. organic matter) such as litter fall and root exudates. The product of S(θ_i) and C_i represents the organic matter (in kg C m⁻³) available for microbial decomposition under anoxic conditions. When the soil surface is not flooded (Fig. 2b), dry soil layers (S(θ_i)=0) are assumed to be predominantly oxic and not producing CH₄ (P_i=0), mostly due to aerobic microbes dominating the competition for organic matter which results in the starvation of methanogens (Segers, 1998).

The global factor r is the specific CH₄ production rate (in kg kg⁻¹ s⁻¹), which can be defined as the mass of CH₄-C that is produced per kilogram of available soil carbon per unit of time. A meta-analysis of incubated soil samples from various anaerobic landscapes indicates that r can vary between 0.3 and 27.2 µg of CH₄-C per gram of soil carbon per day (equivalent to the range from 3.5 × 10⁻¹² to 3.1 × 10⁻¹⁰ kg kg⁻¹ s⁻¹) depending on the landscape type, relative water table position, and soil depth (Treat et al., 2015). In this first version of WETMETH, we combine all possible pathways for CH₄ production in wetlands (see Sect. 2.1) without distinguishing fast and slow pathways. Section 4.1 discusses the choice of the value for r as part of the model calibration.

Table 1 Model parameters for methane production and oxidation.

^a This value is equivalent to 22.8 µg CH₄-C produced per gram of soil carbon per day. ^b A temperature-dependent Q₁₀, approximating 4.2 for a wide range of temperatures, is used instead (see Appendix A). ^c The reference temperature is used to define an optimal temperature for CH₄ production (see Appendix A).

The expression $Q_{\mathrm{10}}^{\frac{T_i-T_{\mathrm{0}}}{\mathrm{10}}}$, which depends on the average layer temperature T_i (in kelvin, K) and a baseline temperature T₀ (273.15 K), represents the temperature dependency of CH₄ production expressed with a Q₁₀ coefficient as commonly done to approximate the sensitivity of biological processes to a temperature change of 10 K (Hegarty, 1973). While some biological processes double in rate with a warming of 10 K, several studies report a higher temperature sensitivity for CH₄ production (i.e. Q₁₀ > 2), although with large uncertainties (Lupascu et al., 2012; Sjögersten et al., 2018; Walz et al., 2017; Whalen, 2005). Nevertheless, a meta-analysis of temperature-response studies suggests an average Q₁₀ of about 4.2 for CH₄ production in pure cultures of methanogens (Hoehler and Alperin, 2014; Yvon-Durocher et al., 2014) in agreement with previous estimates (Blodau, 2002). In order to account for uncertainties with this coefficient and define the occurrence of an optimal temperature for CH₄ production (Dunfield et al., 1993; Metje and Frenzel, 2007; Schipper et al., 2014), a temperature-dependent Q₁₀ is considered in WETMETH. Its mathematical formulation is $Q_{\mathrm{10}}\left(T_i\right)=\mathrm{1.7}+\mathrm{2.5}\tanh \left[\mathrm{0.1}\right(T_{\mathrm{ref}}-T_i\left)\right]$, where T_ref=308.15 K is a reference temperature that is used to define an optimal temperature for CH₄ production (Table 1). This formulation is defined by analogy to a mathematical expression used for soil respiration in another study (Wu et al., 2016), and it enables one to account for an optimal temperature for CH₄ production of ∼ 300.15 K (i.e. 27^∘ C), which is consistent with previous studies (Dunfield et al., 1993; Metje and Frenzel, 2007). Additional information on this formulation and its implications for the temperature dependency of CH₄ production are provided in Appendix A. Furthermore, CH₄ production in WETMETH is assumed to shut down in frozen soil layers, although research suggests that slow microbial activity can occur at temperatures below 273.15 K (Panikov and Dedysh, 2000; Rivkina et al., 2004).

The expression $\exp \left(-\frac{z_i}{{\mathit{\tau }}_{\mathrm{prod}}}\right)$, which depends on the depth of soil layer i relative to the surface (z_i in m, positive downwards), describes the declining effect of various environmental controls on CH₄ production with depth that are generally unresolved by ESMs. These environmental factors include the quality of organic matter and the spread of methanogens among other factors (Bridgham et al., 2013; Koven et al., 2015; Treat et al., 2015; Walz et al., 2017; Wild et al., 2016). Here, τ_prod (in m) is a scaling parameter for CH₄ production. The choice of the value for τ_prod is discussed later as part of the model calibration (see Sect. 4.1).

The total amount of CH₄ produced in the soil column (P in kg C m⁻² s⁻¹) is calculated as follows:

$$\begin{array}{}\text{(2)}& P=\underset{i=\mathrm{1}}{\overset{i=k}{\int }}{P}_i\mathrm{d}{z}_i,\end{array}$$

where P_i (in kg C m⁻³ s⁻¹) is the rate of CH₄ production in the soil layer i from Eq. (1), dz_i (in m) is the thickness of the soil layer i, and k represents the bottom-most soil layer. This amount of CH₄ (P) is then subject to oxidation in transit to emission into the atmosphere.

3.1.2 Parameterization of methane oxidation and net methane emissions

Methane oxidation is parameterized based on the amount of CH₄ produced in the soil column and the relative thickness of the oxic zone. Specifically, the total amount of CH₄ oxidized in the soil column (O_x in kg C m⁻² s⁻¹) and net CH₄ emissions to the atmosphere (E in kg C m⁻² s⁻¹) are calculated as

$$\begin{array}{}\text{(3)}& {\displaystyle }{O}_x=P\left(\mathrm{1}-\exp \left(-{\displaystyle \frac{z_{\mathrm{oxic}}}{{\mathit{\tau }}_{\mathrm{oxid}}}}\right)\right),\text{(4)}& {\displaystyle }E=P-O_x,\end{array}$$

which is equivalent to the following expression:

$$\begin{array}{}\text{(5)}& E=P\exp \left(-{\displaystyle \frac{z_{\mathrm{oxic}}}{{\mathit{\tau }}_{\mathrm{oxid}}}}\right),\end{array}$$

where P (in kg C m⁻² s⁻¹) is the total amount of CH₄ produced in the soil column as defined in Eq. (2), z_oxic (in m) is the relative depth (positive downwards) to the oxic–anoxic interface (Fig. 2), and τ_oxid (in m) is a scaling parameter for CH₄ oxidation. As for τ_prod, the choice of the value for τ_oxid is discussed as part of the model calibration (see Sect. 4.2).

Regarding z_oxic, we assume that O₂ may be present in soil layers unsaturated with water as well as in a shallow oxic–anoxic transition zone within the uppermost soil layer saturated with water (Fig. 2). In this first development of WETMETH, we consider a constant thickness (z_oatz) of 0.05 m for the oxic–anoxic transition zone, with its bottom defined as the oxic–anoxic interface (Frolking et al., 2002; Singleton et al., 2018). When the soil surface is inundated, z_oatz is identical to z_oxic (Fig. 2a). Otherwise, z_oatz is only a fraction of z_oxic (Fig. 2b). The penetration of O₂ into the soil and its dynamics with changing moisture conditions can be complex depending on site-specific factors such as the soil composition (Estop-Aragonés et al., 2012) and the presence of vascular plants (Brune et al., 2000). In addition, methanotrophs may be present at depth (> 0.05 m) below the water table, probably following some adaptation to low-O₂ conditions (Singleton et al., 2018). Nevertheless, the approach applied here for z_oxic is reasonable for ESMs not resolving O₂ dynamics and microbial communities in the soil.

For Eq. (3), the expression $\left(\mathrm{1}-\exp \left(-\frac{z_{\mathrm{oxic}}}{{\mathit{\tau }}_{\mathrm{oxid}}}\right)\right)$ represents the fraction of P that gets oxidized in transit to emission into the atmosphere. Various studies report estimates of CH₄ oxidation as a fraction of produced CH₄ in the soil column (Blazewicz et al., 2012; Le Mer and Roger, 2001; Roslev and King, 1996; Segers, 1998; Singleton et al., 2018). From sample-to-sample and site-to-site, however, CH₄ oxidation exhibits a broad range of values ranging from less than 20 % to more than 95 % depending on the sampled soil depth ranges, whether or not potential CH₄ oxidation under anoxic conditions is considered, the monitored transport mechanisms for CH₄, and many other factors (Blazewicz et al., 2012; Couwenberg et al., 2010; Jauhiainen et al., 2005; Kwon et al., 2019; Le Mer and Roger, 2001; Moosavi and Crill, 1998; Roslev and King, 1996; Segers, 1998; Singleton et al., 2018; Whalen, 2005). Nevertheless, the largest fractions of oxidized CH₄ are generally associated with the deepest water tables or oxic–anoxic interfaces (Bridgham et al., 2013; Couwenberg et al., 2010; Jauhiainen et al., 2005; Roslev and King, 1996; Segers, 1998; Whalen, 2005).

The parameterization described in Eq. (3) is a simple approach for characterizing CH₄ oxidation in the soil column. Such a parameterization is practical when there is little knowledge on the soil chemistry (e.g. O₂ and alternate electron acceptors), the dynamics of methanotrophs, and other environmental factors exerting a control on CH₄ oxidation (Blazewicz et al., 2012; Blodau, 2002; Dean et al., 2018; Kwon et al., 2019; Singleton et al., 2018; Smemo and Yavitt, 2011). Most importantly, this parameterization considers the net effect of all mechanisms transporting CH₄ from the anoxic soil layers where the gas is produced to the atmosphere. The oxidized CH₄ is assumed to produce CO₂ that becomes part of the soil respiration routinely simulated by ESMs.

3.2 The embedding Earth system model

WETMETH has been embedded in the University of Victoria Earth System Climate Model (UVic ESCM), an Earth system model of intermediate complexity (EMIC) (Weaver et al., 2001). A modified version of the EMIC based on UVic ESCM 2.9 (Eby et al., 2009) is used here. The UVic ESCM consists of a 3-D ocean general circulation model that is coupled to a dynamic–thermodynamic sea ice model, a 2-D (vertically integrated) energy–moisture balance model for the atmosphere, and a land surface model (Weaver et al., 2001). The land surface model is a modified version of the Met Office Surface Exchange Scheme (MOSES) with 14 ground layers of unequal thickness extending down to a depth of 250 m that can simulate permafrost processes such as freeze–thaw dynamics (Avis et al., 2011). The top eight ground layers (∼ 10 m in total depth) are soil layers and contribute to the water cycle, whereas the bottom six ground layers are bedrock layers (Avis et al., 2011). In the hydraulically active layers, porosity and permeability are determined based on the relative abundance of prescribed sand, clay, and silt-sized particles. Water phase changes are determined over a range of soil temperatures to determine the fraction of frozen and unfrozen water in the ground (Avis et al., 2011). All components of the UVic ESCM have a horizontal grid resolution of 3.6^∘ in longitude and 1.8^∘ in latitude (Eby et al., 2009; Weaver et al., 2001).

Wetlands in the UVic ESCM are identified in grid cell areas based on soil moisture content and topography. Model grid cells in which wetlands can occur are those with unfrozen soil moisture contents greater than 65 % of the saturated moisture content in the upper soil layer for at least 1 d in a year (Avis et al., 2011). Instead of using a fixed global threshold value for topography (Avis et al., 2011), the version of the UVic ESCM used here identifies wetland coverage at the sub-grid scale following a TOPMODEL approach for global models (Gedney and Cox, 2003). Appendix B describes a minor modification applied to this TOPMODEL approach. Section 5.1 presents an evaluation of wetlands simulated by the UVic ESCM.

The UVic ESCM includes a representation of the global carbon cycle. The terrestrial carbon cycle is simulated using the Top-down Representation of Interactive Foliage and Flora including Dynamics (TRIFFID), a dynamic global vegetation model that is coupled to the land surface model (Avis et al., 2011; Meissner et al., 2003). TRIFFID defines the state of the terrestrial biosphere in terms of soil carbon as well as the structure and coverage of five plant functional types (PFTs): broadleaf trees, needleleaf trees, shrubs, C₃ grasses, and C₄ grasses (Cox, 2001; Matthews et al., 2004; Meissner et al., 2003). Terrestrial carbon gain occurs through photosynthesis that is simulated as a function of atmospheric CO₂ concentration, shortwave radiation, air temperature, humidity, and soil moisture. Soil carbon gain occurs through litter fall and vegetation mortality. The present-day permafrost carbon pool is simulated by the UVic ESCM following a method that approximates the effect of long-term freeze–thaw cycles on the vertical distribution of carbon in permafrost-affected soils, a process referred to as cryoturbation (MacDougall and Knutti, 2016). Soil carbon can occur in the top six ground layers (∼ 3.35 m in total depth). Terrestrial carbon loss occurs through autotrophic respiration by plants and heterotrophic respiration by soil microbes (Matthews et al., 2004; Meissner et al., 2003). By configuration, permafrost carbon can only be lost through microbial respiration, and this heterotrophic respiration is assumed to shut down in frozen soil layers (MacDougall et al., 2012; MacDougall and Knutti, 2016). Through TRIFFID, all terrestrial carbon fluxes in the UVic ESCM are integrated with a 30 d time step (Meissner et al., 2003).

The marine carbon cycle in the UVic ESCM is represented with organic and inorganic carbon cycle models (Eby et al., 2009). The organic carbon cycle is based on marine biology simulated with a nutrient–phytoplankton–zooplankton–detritus (NPZD) ecosystem model (Schmittner et al., 2008). The inorganic carbon cycle model simulates the air–sea exchange of CO₂ and ocean carbonate chemistry following the protocols of the Ocean Carbon-Cycle Model Intercomparison Project (OCMIP) (Orr, 1999; Weaver et al., 2001). Dissolved inorganic carbon is treated as a passive tracer that is subject to ocean circulation (Weaver et al., 2001). Carbonate dissolution in ocean sediments is simulated with a model of respiration in marine sediments (Archer, 1996; Eby et al., 2009).

3.3 Model simulations

For this research, three series of model simulations are performed with the UVic ESCM in its standard fully coupled mode and including WETMETH parameterizations:

1. Firstly, the UVic ESCM is spun up for ∼ 5000 years at year 1850 conditions to allow the model to reach an equilibrium climate state representing the pre-industrial period.

2. Secondly, a transient run from 1850 to 2019 is performed in order to evaluate the model performance. This transient run is based on prescribed CO₂ concentration and other forcing data from the Phase 5 of the Coupled Model Intercomparison Project (CMIP5) (Taylor et al., 2012). The UVic ESCM is driven by historical data from 1850 to 2005 and by Representative Concentration Pathway (RCP) 8.5 data from 2006 to 2019. Figure S1 in the Supplement illustrates how the simulated historical climate conditions compare to observations in terms of global mean surface air temperature.

3. Thirdly, a set of transient runs from 2000 to 2009 is performed to analyze the model sensitivity to poorly constrained parameters. This set of model simulations (sensitivity runs) is performed by perturbing values of poorly constrained parameters associated with wetland CH₄ processes.

4 Choice of model parameter values

Here, we describe the choice of three WETMETH parameters (r and τ_prod for CH₄ production; τ_oxid for CH₄ oxidation) as part of the model calibration. These model parameters are tuned to observations from northern high-latitude regions due to the scarcity of large-scale datasets from other regions. The model calibration against northern observations is based on the assumption that tuned parameter values will be valid across the globe, which is an important limitation as will be discussed later. Nonetheless, this approach is deemed reasonable given the present state of data availability. Section 5.1 describes northern wetlands simulated by the UVic ESCM as part of the model validation.

4.1 Methane production parameters

Parameters for CH₄ production in WETMETH are calibrated against maximum CH₄ production rates measured in laboratory incubations of soil samples from several anaerobic environments across northern high-latitude regions (> 50^∘ N). These potential CH₄ production rates are obtained from a synthesis dataset, which includes information on other environmental variables such as the relative depth of the soil samples (Treat et al., 2015).

To allow a fair model–data comparison, measured CH₄ production rates with corresponding soil bulk density from the sites of origin are converted into kilograms of carbon per cubic metre per second (kg C m⁻³ s⁻¹) (see Appendix C​​​​​​​). Furthermore, measurements from landscapes identified as uplands and lakes (in the dataset) are excluded from the dataset used in this model calibration. The remaining measurements are potential CH₄ production rates in soil samples from landscapes identified (in the dataset) as wetlands, floodplains, and lowlands across Alaska.

In order to set values for r and τ_prod from Eq. (1), the depth profile of simulated CH₄ production rates across Alaska for the year 2000 is tuned to that of the measurements. By setting r to 22.8 µg CH₄–C produced per gram of soil carbon per day (equivalent to 2.6 × 10⁻¹⁰ kg kg⁻¹ s⁻¹) and τ_prod to 0.75 m, we obtain a depth profile of simulated CH₄ production rates that compares fairly well to that of potential CH₄ production rates from the laboratory incubations (Fig. 3). These default values for r and τ_prod are listed in Table 1. Section 6 presents a sensitivity analysis on these model parameters.

Figure 3 Vertical profiles of simulated and potential CH₄ production rates from wetlands across Alaska. Potential CH₄ production rates are measurements from laboratory incubations of soil samples collected from various anaerobic ecosystems (Treat et al., 2015). Both simulated and measured CH₄ production rates are shown here with a log-transformed axis (base-10 logarithmic scale).

4.2 Methane oxidation parameter

Unlike for CH₄ production, there are no published large-scale measurements of CH₄ oxidation rates that could be used in this research for the calibration of CH₄ oxidation. For that reason, CH₄ oxidation in WETMETH is indirectly calibrated via CH₄ emissions. A synthesis dataset of seasonal and annual CH₄ emissions from various terrestrial sites across temperate, boreal, and Arctic regions is used to this end (Treat et al., 2018). The model calibration focuses on annual CH₄ emissions from sites north of 50^∘ N for which many data points are available in the dataset.

While most data points are from direct measurements of CH₄ emissions, some data points are associated with different modelling methods for estimating CH₄ emissions (Treat et al., 2018). To allow a fair model–data comparison, only data points associated with direct measurements of CH₄ emissions are included in the model calibration. Furthermore, measurements from lakes, uplands, and alpine landscapes are excluded from this model calibration. In particular, the exclusion of data points from uplands and alpine landscapes sorts out measurements of terrestrial CH₄ uptake (negative CH₄ flux). The retained data points (n = 119) include measurements by chambers (85.7 %), flux towers (13.4 %), and a combination of flux towers and chambers (0.8 %).

The model calibration in this section aims at choosing a value of τ_oxid from Eq. (4) such that the range (minimum – maximum) of annual CH₄ emissions across northern wetlands (> 50^∘ N) simulated by the UVic ESCM is comparable to that of annual CH₄ emissions from the data points (0.1–60.6 g CH₄ m⁻² yr⁻¹). By setting τ_oxid to 0.0146 m, we constrain simulated CH₄ emissions from northern wetlands (specifically, grid cell CH₄ emissions divided by the inundated fraction of the grid cell) from 2000 to 2009 in the range of 0.04–65.6 g CH₄ m⁻² yr⁻¹. This tuned value for τ_oxid is listed in Table 1 and implies that ∼ 97 % of the CH₄ produced in the soil column gets oxidized in transit to emission when the soil surface is inundated. Section 6 presents a sensitivity analysis on this model parameter.

5 Evaluation of the model performance

5.1 Wetlands

Figure 4 shows the latitudinal distribution of wetland areas simulated by the UVic ESCM in comparison to two global datasets. The first dataset is Global Inundation Extent from Multi-Satellites (GIEMS), which is based on remotely sensed inundation areas (Papa et al., 2010; Prigent et al., 2001, 2007, 2012). The second dataset is Surface Water Microwave Product Series – Global Lakes and Wetlands Database (SWAMPS-GLWD), which is based on a combination of information from satellites and maps of inundated areas in order to reduce uncertainties associated with the distribution of global wetlands (Poulter et al., 2017). The comparison between the model and the datasets is done over 2000-2007, which is the overlap period for the datasets. Over this period, the UVic ESCM simulates an annual maximal extent of ∼ 12.6×10⁶ km² for global wetlands, whereas GIEMS and SWAMPS-GLWD estimate ∼ 9.3×10⁶ km² and ∼ 10.6×10⁶ km², respectively.

Figure 4 Latitudinal distribution of wetland areas simulated by the UVic ESCM over the 2000–2007 period in comparison to two global datasets: GIEMS and SWAMPS-GLWD. The comparison period corresponds to the overlap period for the two datasets. The wetland areas are summed across latitude bins of 3^∘.

The UVic ESCM agrees better with SWAMPS-GLWD in regions north of 40^∘ N, although with some underestimations around 55^∘ N, and relatively well with GIEMS between 20 and 40^∘ S (Fig. 4). However, the model simulates overly small wetland areas between 20 and 30^∘ N when compared with both GIEMS and SWAMPS-GLWD. While our model could be underestimating wetland areas in this latitude zone, inundated areas estimated by GIEMS include rice paddies which prevail in tropical and sub-tropical regions (Prigent et al., 2007, 2012). Rice paddies are likely not represented in SWAMPS-GLWD, as there were efforts to only include natural wetlands during the development of this dataset (Poulter et al., 2017). In comparison to GIEMS and SWAMPS-GLWD, our model simulates small wetland areas in Southeast Asia especially near Bangladesh (Figs. 5, 6).

Figure 5 Average wetland extents (inundated fractions of grid cells) across the globe over the 2000–2007 period as simulated by the UVic ESCM (a) in comparison to two datasets: (b) GIEMS and (c) SWAMPS-GLWD. The datasets are regridded to 3.6^∘ × 1.8^∘ for a fair comparison with the UVic ESCM. The comparison period corresponds to the overlap period for the two datasets.

Figure 6 Differences in global wetland extents (inundated fractions of grid cells) between two datasets (GIEMS and SWAMPS-GLWD) and the UVic ESCM over the 2000-02007 period: (a) SWAMPS-GLWD – GIEMS, (b) UVic ESCM – GIEMS, and (c) UVic ESCM – SWAMPS-GLWD. The comparison period corresponds to the overlap period for the two datasets.

Between 20^∘ N and 20^∘ S, the UVic ESCM simulates a bimodal distribution of the wetland extent that is consistent with the two datasets although the model simulates overly large wetland areas (Fig. 4). Unlike for GIEMS and SWAMPS-GLWD, wetlands simulated by the UVic ESCM are widespread in Amazonia as well as in West and Central Africa (Figs. 5, 6). Although the UVic ESCM could be overestimating the extent of wetlands in some of these equatorial regions, it is possible that GIEMS and SWAMPS-GLWD do not detect inundated areas in densely forested regions due to forest canopies. Recent studies suggest that tropical wetlands are commonly underestimated in large-scale datasets (Dargie et al., 2017; Gumbricht et al., 2016).

Conversely, it is possible that the UVic ESCM overestimates tropical wetland areas due to soil hydraulic properties being unrepresented in the model. A potential cause for the overestimation of tropical wetlands in our model is the standard approach for simulating global hydrology in land surface models based on the concentration of only sand, clay, and silt in the soil. A recent study suggests that the inclusion of Ferralsols (weathered soils with micro-aggregated particles that are common in the humid tropics) in a global terrestrial model can help improve the simulation of tropical wetlands (Gedney et al., 2019).

Outside of the tropics, the UVic ESCM does a better job at simulating the distribution of wetlands in sub-Arctic and Arctic regions (Fig. 7). The model simulates the occurrence of wetlands (i.e. surface inundation) across the West Siberian Lowlands (WSL) in Russia, the Hudson Bay Lowlands (HBL) in Canada, and over other parts of northern Canada in agreement with both SWAMPS-GLWD and GIEMS (Fig. 7). However, some disagreements between the UVic ESCM and the two datasets can also be identified: (i) in comparison with GIEMS, the UVic ESCM simulates more wetland area in the Hudson Bay Lowlands (HBL) as well as widespread wetlands in parts of northern Eurasia (Figs. 7b, S2b); (ii) in comparison with SWAMPS-GLWD, the model simulates less wetland area over the WSL and northern Canada, including the HBL, and more wetland area in parts of Europe (Figs. 7c, S2c).

Figure 7 Average wetland extents (inundated fractions of grid cells) north of 45^∘ N over the 2000–2007 period as simulated by the UVic ESCM (a) in comparison to two datasets: (b) GIEMS and (c) SWAMPS-GLWD. The datasets are regridded to 3.6^∘ × 1.8^∘ for a fair comparison with the UVic ESCM. The comparison period corresponds to the overlap period for the two datasets.

Statistical analyses show that (i) the UVic ESCM agrees better with SWAMPS-GLWD than with GIEMS at both the regional and global scale, and (ii) the model shows better agreement with the two datasets across northern regions than at the global scale. For details on the statistical evaluation, the reader is referred to Table S1.

5.2 Wetland methane emissions

Given the relative coarse grid resolution of the UVic ESCM, the model validation with respect to wetland CH₄ emissions focuses on large-scale emissions such as regional, zonal, and global emissions. Moreover, this model validation focuses on northern high-latitude regions because observations and estimates of wetland CH₄ emissions from other regions (e.g. the tropics) are scarce. This focus is further justified by the fact that our model better simulates the distribution of wetlands in northern high-latitude regions than in the tropics (see Sect. 5.1). Indeed, the extent of wetlands is a major control for wetland CH₄ emissions simulated by process-based models and probably the primary contributor to related uncertainties (Melton et al., 2013; Saunois et al., 2020; B. Zhang et al., 2017).

5.2.1 Northern high-latitude emissions

The UVic ESCM simulates total CH₄ emissions from northern wetlands that are in the range of recent estimates. Over the 2013–2014 period, the model simulates mean annual emissions of 33.2 Tg CH₄ yr⁻¹ for wetlands north of 45^∘ N (Table 2). These CH₄ emissions are consistent with estimates from recent upscaled flux measurements (UFMs) over the same period based on a random forest (RF) algorithm and three wetland maps (Peltola et al., 2019): 30.6 ± 9.2 Tg CH₄ yr⁻¹ (RF-DYPTOP), 31.7 ± 9.4 Tg CH₄ yr⁻¹ (RF-PEATMAP), and 37.6 ± 11.8 Tg CH₄ yr⁻¹ (RF-GLWD) (Table 2). Table S2 shows that the UVic ESCM has no preferential agreement with one of the three UFMs.

Table 2 Mean annual wetland CH₄ emissions simulated by the UVic ESCM in comparison to estimated emissions from the literature. All emissions are reported in teragrams of methane per year (Tg CH₄ yr⁻¹), and uncertainties are provided for estimates from the literature. Three periods are used to allow a fair comparison between the UVic ESCM and estimates from the literature where possible: 2008–2017, as in the latest global CH₄ budget report (Saunois et al., 2020); 2013–2014, as for recent upscaled flux measurements across the northern high latitudes (Peltola et al., 2019); and 1993–2004, as for the WETCHIMP model ensemble (Melton et al., 2013). Principal methods used in the different references for estimates are reported in the last column: top-down (TD) methods including inverse models (IM), and bottom-up (BU) methods including upscaled measurements (UM) and process-based models (PM).

^∗ These reported estimates are model ensemble means. For the West Siberian Lowlands, the range between the inverse models is 3.1–9.8 Tg CH₄ yr⁻¹ (Bohn et al., 2015). For tropical wetlands, the range between the process-based models is 85–184 Tg CH₄ yr⁻¹ (Melton et al., 2013).

Figure 8 shows the spatial distribution of simulated CH₄ emissions in comparison to the three UFMs. When compared with each other, the three UFMs exhibit substantial differences primarily attributed to the distinct wetland distributions (Peltola et al., 2019). Considering the general pattern and magnitude of wetland CH₄ emissions, the UVic ESCM agrees with either two or all three UFMs over key source regions such as the Hudson Bay Lowlands (HBL), the West Siberian Lowlands (WSL), western Europe, and south-central Canada (Fig. 8).

Figure 8 Average CH₄ emissions from wetlands north of 45^∘ N over the 2013–2014 period as simulated by the UVic ESCM (a) in comparison to three datasets (upscaled flux measurements): (b) RF-DYPTOP, (c) RF-GLWD, and (d) RF-PEATMAP. The datasets are regridded to 3.6^∘ × 1.8^∘ for a fair comparison with the UVic ESCM. The comparison period corresponds to the overlap period for the three datasets.

The UVic ESCM simulates less CH₄ emissions over parts of northeastern Canada and Fennoscandia compared with the UFMs (Fig. 8). However, the three UFMs do not necessarily agree on both the distribution and magnitude of wetland CH₄ emissions in these regions. Furthermore, the UVic ESCM does not simulate wetland CH₄ emissions in southern Eurasia (45–60^∘ N; 40–135^∘ E), whereas the three UFMs suggest that CH₄ can be emitted from sporadic wetlands in this region (Fig. 8). Overall, the mismatch between the UFMs and our model in terms of northern CH₄ emissions can be primarily attributed to differences in the wetland extent but also to the spatial distribution of soil carbon simulated by the UVic ESCM (MacDougall and Knutti, 2016).

In terms of mean annual emissions from key source regions, the UVic ESCM simulates 2.9 Tg CH₄ yr⁻¹ for the Hudson Bay Lowlands (HBL) over the 2013–2014 period (Table 2). Although these emissions are lower than estimates by the three UFMs (3.1–6.5 Tg CH₄ yr⁻¹) (Peltola et al., 2019), estimates by inverse models (2.0–3.4 Tg CH₄ yr⁻¹) over this region are comparable to our model results (Miller et al., 2014; Pickett-Heaps et al., 2011; Thompson et al., 2017). Furthermore, the UVic ESCM simulates total wetland emissions of 4.1 Tg CH₄ yr⁻¹ for the West Siberian Lowlands (WSL) over the 2013–2014 period (Table 2). Regional estimates based on the three UFMs are higher (4.9–8.5 Tg CH₄ yr⁻¹) than our model results over the same period (Peltola et al., 2019), whereas previous observation-based estimates for the WSL suggest regional wetland emissions (3.9 ± 1.3 Tg CH₄ yr⁻¹) that are similar to our model results (Glagolev et al., 2011). Estimates by inverse models over the WSL are relatively high but comparable to our model estimates (Table 2): 6.1 ± 1.2 Tg CH₄ yr⁻¹ (Bohn et al., 2015) and 6.9 ± 3.6 Tg CH₄ yr⁻¹ (Thompson et al., 2017).

The UVic ESCM is also evaluated with respect to wetland CH₄ emissions over the 2000–2009 and 2008–2017 decades, which are both reference periods for the latest global CH₄ budget report (Saunois et al., 2020). For wetlands north of 40^∘ N, the UVic ESCM simulates emissions of 37.7 Tg CH₄ yr⁻¹ over the 2000–2009 decade and 38.5 Tg CH₄ yr⁻¹ over the 2008–2017 decade. These wetland CH₄ emissions are consistent with recent estimates (37.4 ± 7.2 Tg CH₄ yr⁻¹) from data-constrained model ensembles over the same region (Treat et al., 2018). For wetlands north of 45^∘ N, the model simulates total emissions that are in the range of estimates for the 2013–2014 period discussed earlier (32.4 Tg CH₄ yr⁻¹ over 2000–2009 and 33.1 Tg CH₄ yr⁻¹ over 2008–2017). For Pan-Arctic wetlands (> 60^∘ N), the UVic ESCM simulates emissions of 17.4 Tg CH₄ yr⁻¹ over the 2000–2009 decade and a similar amount over the 2008–2017 decade (Table 2). These wetland CH₄ emissions correspond to the upper limit of bottom-up estimates (2–18 Tg CH₄ yr⁻¹) from the latest global CH₄ budget report (Saunois et al., 2020).

Figure 9 Seasonal variations of CH₄ emissions from wetlands north of 45^∘ N over the 2013–2014 period as simulated by the UVic ESCM in comparison to three upscaled flux measurements (RF-DYPTOP, RF-GLWD, and RF-PEATMAP). The dashed lines show the uncertainty range for the upscaled flux measurements.

Figure 9 shows seasonal cycles of CH₄ emissions from wetlands north of 45^∘ N over the 2013–2014 period as simulated by the UVic ESCM and estimated from the three UFMs (Peltola et al., 2019). The pattern and magnitude of simulated seasonal emissions compare well to that of the UFMs. For both the model and UFMs, minimal emissions vary between 0.2 and 0.6 Tg CH₄ per month and occur in December, whereas peak emissions are well below 10 Tg CH₄ per month and occur in July (Fig. 9). However, simulated peak emissions (∼ 8.5 Tg CH₄ per month) are relatively higher than peak emissions for the UFMs (range of best estimates from 5.6 to 7.5 Tg CH₄ per month). Moreover, in comparison to the three UFMs, the UVic ESCM simulates lower CH₄ emissions between December and May but higher CH₄ emissions between July and September (Fig. 9).

The UVic ESCM simulates the occurrence of wetland CH₄ emissions during the non-growing season. For wetlands north of 45^∘ N, our model simulates total emissions of 2.1 Tg CH₄ yr⁻¹ between November and March. The UFMs predict total emissions of 4.6–10.2 Tg CH₄ yr⁻¹ during these cold months (Peltola et al., 2019). For wetlands north of 60^∘ N, the UVic ESCM simulates emissions of 1.2 Tg CH₄ yr⁻¹ from October through May in agreement with recent estimates (1.6 ± 0.1 Tg CH₄ yr⁻¹) from data-constrained model ensembles for these months (Treat et al., 2018). Based on our calculations, the three UFMs predict about 3.5–4.5 Tg CH₄ yr⁻¹ emitted from wetlands north of 60^∘ N between October and May. Overall, this analysis shows that WETMETH is capable of simulating non-negligible CH₄ emissions from northern wetlands during cold months as emphasized by recent studies (Treat et al., 2018; Zona et al., 2016).

5.2.2 Global emissions

The UVic ESCM simulates total emissions of 155.1 and 158.6 Tg CH₄ yr⁻¹ from global wetlands over the 2000–2009 and 2008–2017 decades, respectively. According to the latest global CH₄ budget report, these wetland emissions are in the mid-range of bottom-up estimates (102–179 and 102–182 Tg CH₄ yr⁻¹) but close to the lower limit of top-down estimates (153–196 and 155–200 Tg CH₄ yr⁻¹) over the two respective decades (Saunois et al., 2020). Previous bottom-up estimates are significantly high (Melton et al., 2013; Saunois et al., 2016), which is primarily due to possible double counting of emissions from wetlands and other inland water areas (Saunois et al., 2020; Thornton et al., 2016) in addition to uncertainties associated with the extent of wetlands and model parameterizations (Melton et al., 2013). Table 2 summarizes the comparison between the model results and estimates from the latest global CH₄ budget report for the 2008–2017 decade.

Figure 10 Average methane emissions from global wetlands over the 2001–2004 period as simulated by the UVic ESCM (a) in comparison to three process-based model ensembles: (b) GCP-CH4, (c) WetCHARTs, and (d) WETCHIMP. The model ensembles are regridded to 3.6^∘ × 1.8^∘ for a fair comparison with the UVic ESCM. The comparison period corresponds to the overlap period for the three model ensembles.

Figure 10 shows the spatial distribution of simulated wetland CH₄ emissions over the 2001–2004 period in comparison to three process-based model ensembles: GCP-CH4 (Poulter et al., 2017), WetCHARTs (Bloom et al., 2017), and WETCHIMP (Melton et al., 2013). The UVic ESCM simulates few CH₄-emitting areas over Southeast Asia in comparison to the three model ensembles. The potential underestimation of wetland CH₄ emissions in that region is associated with the relatively few wetland areas simulated by the UVic ESCM (see Sect. 5.1). In tropical Africa, our model simulates too many CH₄-emitting locations in comparison with the model ensembles (Fig. 10), which is also associated with the distribution of simulated wetlands (see Sect. 5.1). Nevertheless, the UVic ESCM simulates the occurrence of wetland CH₄ emissions in key source regions such as the Amazon and Congo River basins, South Sudan (Sudd swamps), and the Indonesian islands (Fig. 10). For the Amazon and Congo River basins, however, the UVic ESCM simulates lower wetland CH₄ emissions than predicted by the model ensembles (Fig. 10). This can be due to either the consideration of an optimal temperature for CH₄ production (around 27^∘ C) in our model unlike many other process-based models or the fact that model parameters in this study are tuned to northern estimates.

Figure 11a shows the latitudinal distribution of simulated wetland CH₄ emissions compared with the model ensembles. Interestingly, although GCP-CH4 and WetCHARTs are based on the same wetland dataset (SWAMPS-GLWD) (Bloom et al., 2017; Poulter et al., 2017), their zonal wetland CH₄ emissions are very different, especially near the Equator and across northern high-latitude regions (Fig. 11a).

Figure 11 (a) Latitudinal distribution of wetland methane emissions simulated by the UVic ESCM over the 2001–2004 period in comparison to three process-based model ensembles: GCP-CH4, WetCHARTs, and WETCHIMP. The comparison period corresponds to the overlap period for the three model ensembles. (b) Latitudinal emission intensity (methane emissions per unit of wetland area) simulated by the UVic ESCM over the 2001–2004 period in comparison to the three process-based model ensembles. GCP-CH4 and WetCHARTs both use SWAMPS-GLWD as prescribed wetlands. The wetland methane emissions and emission intensities are summed across latitude bins of 3^∘.

Using the three model ensembles as reference, the UVic ESCM simulates significantly lower wetland CH₄ emissions around the Equator (Fig. 11a), despite the fact that the model simulates overly large equatorial wetland areas (Fig. 4). In fact, wetland emission intensities (emissions per unit of wetland area) by the UVic ESCM are lower than those by the model ensembles between 10^∘ S and 10^∘ N (Fig. 11b) due to relatively large wetland areas but small CH₄ emissions in equatorial regions (cf. Fig. 4 and Fig. 11a). As previously discussed, the relatively small CH₄ emissions simulated by the UVic ESCM in equatorial regions can be associated with either the optimal temperature for CH₄ production considered in WETMETH but not in most other process-based models or the fact that model parameters in this study are tuned to northern estimates.

Furthermore, the UVic ESCM simulates more wetland CH₄ emissions between 10 and 20^∘ N than the three model ensembles (Fig. 11a), and this can be attributed to the widespread wetlands in West and Central Africa simulated by our model (Figs. 5, 6). In addition, the UVic ESCM simulates significantly less wetland CH₄ emissions between 20 and 35^∘ N in comparison to the WETCHIMP ensemble (Fig. 11a), and this can be attributed to the relatively small wetland areas simulated by the UVic ESCM in Southeast Asia where some models include agricultural wetlands such as rice paddies. Moreover, wetland emission intensities by the UVic ESCM feature low variability with latitude unlike the three model ensembles (Fig. 11b). Such a relative lack of variability can be attributed to two factors: (i) both wetland areas and CH₄ emissions simulated by the UVic ESCM feature relatively low variability with latitude compared with the datasets and model ensembles (Figs. 4, 11a); and (ii) as previously discussed, our model likely simulates overly large wetland areas but overly small CH₄ emissions around the Equator, implying a lack of variability across tropical latitudes.

Despite the various discrepancies between the UVic ESCM and both model ensembles regarding the distribution of wetland CH₄ emissions in the tropics, our model simulates mean annual CH₄ emissions from tropical wetlands that are in the range of estimates from the literature (Table 2). For the 1993–2004 period, the UVic ESCM simulates tropical wetland CH₄ emissions of 105.5 Tg CH₄ yr⁻¹, whereas the WETCHIMP ensemble predicts 126 ± 31 Tg CH₄ yr⁻¹ (Melton et al., 2013). Another study suggests a lower mean value (90 ± 77 Tg CH₄ yr⁻¹) for wetland CH₄ emissions in the tropics, although with large uncertainties (Sjögersten et al., 2014). Indeed, several studies indicate that wetland CH₄ emissions in the tropics are highly uncertain due to limited ground-based measurements and poorly delimitated wetland extent (Dargie et al., 2017; Gumbricht et al., 2016; Hu et al., 2018; Pangala et al., 2017; Saunois et al., 2020).

6 Model sensitivity to poorly constrained parameters

We performed a set of 30 model runs with perturbed parameter values (sensitivity runs) over the 2000–2009 decade in order to analyze the model sensitivity to poorly constrained parameters (T_ref, r, τ_prod, z_oatz, and τ_oxid). For each parameter, we increased or decreased the default value by 10 %, 20 %, and 30 % while holding values for other parameters constant (fixed to default values). We then compared results from the sensitivity runs to the model simulation with all parameter values set to default values (control run). This comparison focuses on the total simulated global (90^∘ S–90^∘ N), northern (45–90^∘ N), and tropical (30^∘ S–30^∘ N) wetland CH₄ emissions over the 2000-2009 decade.

Our results show that the model sensitivity varies with the different parameters and across regions (Fig. 12). Among the five poorly constrained parameters, the UVic ESCM is most sensitive to perturbations of the two parameters for CH₄ oxidation (z_oatz and τ_oxid) at both the global and regional scale. For z_oatz, a decrease (increase) of the default parameter value by 10 %–30 % results in an augmentation (reduction) of default wetland CH₄ emissions by 41 %–179 % (29 %–64 %) at both the global and regional scale (Fig. 12j–l). For τ_oxid, a decrease (increase) of the default parameter value by 10 %–30 % implies a reduction (augmentation) of default wetland CH₄ emissions by 32 %–77 % (37 %–120 %) at both the global and regional scale (Fig. 12m–o).

Figure 12 Analysis of the model sensitivity to perturbations of poorly constrained parameters: T_ref, r, τ_prod, z_oatz, and τ_oxid. For each parameter, the default value is increased or decreased by 10 %, 20 %, and 30 % while values of other parameters are held constant (to default values). The model sensitivity is analyzed with respect to global (90^∘ S–90^∘ N), northern (45–90^∘ N), and tropical (30^∘ S–30^∘ N) wetland methane emissions. Vertical axes show the ratio of the resulting emissions to the default emissions.

The UVic ESCM is also very sensitive to perturbations of T_ref, but this sensitivity is more pronounced for tropical regions than northern regions (Fig. 12a–c). We recall that T_ref is used to define an optimal temperature for CH₄ production in WETMETH through the Q₁₀ formulation (see Sect. 3.1.1 and Appendix A). For northern regions, a decrease (increase) of T_ref by 10 %–30 % results in a reduction (augmentation) of default wetland CH₄ emissions by 5 %–21 % (3 %–5 %). For tropical regions, however, a decrease (increase) of T_ref by 10 %–30 % results in a reduction (augmentation) of default wetland CH₄ emissions by 34 %–82 % (33 %–75 %). Globally, a decrease (increase) of T_ref by 10 %–30 % results in a reduction (augmentation) of default wetland CH₄ emissions by 26 %–66 % (24 %–55 %). The model sensitivity to perturbations of r is linear across all regions (Fig. 12d–f). Lastly, the model is least sensitive to perturbation of τ_prod across the globe (Fig. 12g–i).

7 Discussions

7.1 WETMETH in the spectrum of wetland methane models

A recent study reviewed 40 models of CH₄ emissions in terrestrial ecosystems (predominantly rice paddies and natural wetlands) and classified them into three categories based on their level of complexity: relatively simple models, relatively mechanistic models, and mechanistic models (Xu et al., 2016). Relatively simple models are those that simulate net CH₄ emissions based on soil carbon or other environmental factors without explicit representations for the different CH₄ production and oxidation pathways as well as mechanisms transporting CH₄ to the atmosphere. Relatively mechanistic models are those that account for at least one transport mechanism for CH₄ release in addition to representing CH₄ production and oxidation with simple functions. Mechanistic models are more comprehensive and explicitly simulate different pathways for both CH₄ production and oxidation, more than two mechanisms for CH₄ release, and their environmental controls. Based on this classification, WETMETH is a relatively simple model in the sense that it does not distinguish pathways for CH₄ production and oxidation as well as the various mechanisms transporting CH₄ to the atmosphere.

Although some wetland CH₄ models are claimed to be embedded in ESMs (Xu et al., 2016), none of these models are currently run in fully coupled models with feedbacks between climate conditions and the global carbon cycle. Most of these models are rather implemented in dynamic vegetation models or uncoupled land surface components of climate models (Arora et al., 2018; Eliseev et al., 2008; Hodson et al., 2011; Riley et al., 2011; Ringeval et al., 2011; Wania et al., 2009). Nonetheless, relatively simple models present the ideal level of complexity for the current generation of ESMs. More complex models generally imply detailed soil chemistry for O₂ and alternate electron acceptors (Riley et al., 2011; Wania et al., 2010), different carbon substrates and their effects on CH₄ production (Grant, 1998; Lovley and Klug, 1986), and an explicit representation of the dynamics of different microbial communities (Grant, 1998; Xu et al., 2015), which all require comprehensive soil chemistry or model parameters that are currently not common in ESMs (Xu et al., 2016). Process parameterizations in mechanistic models generally imply too many degrees of freedom, making it difficult to constrain model parameters against sparse observations. Furthermore, mechanistic models may be too computationally demanding for fully coupled ESM runs without a proportional benefit for large-scale simulations of wetland CH₄ emissions.

The particularity of WETMETH among relatively simple models is that the model accounts for an optimum temperature for CH₄ production, a calibration of depth-dependent CH₄ production rates against potential CH₄ production rates from laboratory incubations, dynamic CH₄ oxidation based on the vertical distribution of soil moisture, and the potential for CH₄ emissions in non-inundated ecosystems with a relatively high level of soil moisture content. In conclusion, WETMETH is simple enough to be compatible with ESMs and yet complex enough to implicitly simulate biogeochemical processes regulating wetland CH₄ emissions.

7.2 Limitations for WETMETH

The developed wetland CH₄ model is associated with several limitations, which are linked to either its level of complexity or the scarcity of large-scale datasets for model calibration:

1. The present state of global wetland modelling assumes generic wetlands without distinguishing their different types (Melton et al., 2013; Poulter et al., 2017). Like many other large-scale models of the current generation, WETMETH would not be appropriate for investigating the contribution from particular wetland types to regional or global CH₄ emissions (Aselmann and Crutzen, 1989).

2. As WETMETH is not based on a comprehensive soil biochemistry module and does not include the different pathways for CH₄ production and oxidation, the model is not suited for investigating the role of specific biological and chemical controls on wetland CH₄ emissions (Bridgham et al., 2013; Kwon et al., 2019).

3. WETMETH does not simulate the contribution from wetland-specific vegetation species to CH₄ emissions, although some of these species can either lead to high emissions (e.g. sedges are vascular plants that can transport CH₄ through their aerenchyma) or low emissions (e.g. mosses are non-vascular plants that have been shown to develop a symbiotic relationship with methanotrophs) (Bridgham et al., 2013; Chen and Murrell, 2010).

4. Ebullition and aerenchyma of vascular plants allow CH₄ produced in wetlands to escape to the atmosphere with little opportunity for oxidation (Segers, 1998; Whalen, 2005). Moreover, stems of woody trees are important conduits for CH₄ emissions in Amazonia, a major source region in the world (Pangala et al., 2017). By considering the net effect of all mechanisms transporting CH₄ to the atmosphere, WETMETH presents a limitation for investigating the relative contribution of transport mechanisms to CH₄ emissions across regions and at the global scale.

5. Methane produced in northern wetlands can be stored underneath frozen soil during the winter and be released abruptly upon spring thaw (Mastepanov et al., 2013; Song et al., 2012). WETMETH does not currently feature such a storage of CH₄ in the soil column, which is probably more relevant for small-scale (site) and short-term (daily) than large-scale (regional) and long-term (seasonal) emissions (Fig. 9).

6. While the existence of an optimal temperature for CH₄ production in wetlands is relatively well established in the literature (Dean et al., 2018), there are currently no estimates of such an optimal temperature for different climate zones across the globe. Previous studies suggest a range of 25–30^∘ C for such an optimal temperature based on measurements of CH₄ production in northern wetlands (Dunfield et al., 1993; Metje and Frenzel, 2007). In WETMETH, we use a global value for this optimal temperature (∼ 27^∘ C) which is assumed to be valid for CH₄ production in both tropical and extra-tropical wetlands (see Sect. 3.1.1 and Appendix A). However, our sensitivity analysis suggests that, in the present-day climate, wetland CH₄ emissions in the tropics are much more dependent on the optimal temperature for CH₄ production than wetland CH₄ emissions in the boreal and Arctic regions (see Sect. 6). The optimal temperature for CH₄ production in WETMETH, along with other factors such as areal wetland extents, contributes to inter-model differences in simulated wetland CH₄ intensities in the tropics (see Fig. 11b).

7. As presented in this study, poorly constrained WETMETH parameters are tuned to estimates from northern high-latitude regions because large-scale datasets from other regions are scarce (see Sect. 4). A strong limitation comes with the assumption that the chosen parameter values are representative for CH₄ production and oxidation across the globe. However, the applied model calibration remains a reasonable approach given the scarcity of observations for wetland CH₄ production, oxidation, and emissions at the global scale.

Despite these limitations and the model simplicity, WETMETH is skilful when it comes to the simulation of mean seasonal, annual, and decadal wetland CH₄ emissions at the regional, hemispheric, and global scale (see Sect. 5.2). The implementation of WETMETH in a fully coupled ESM should advance research on the interactions between climate change and wetland CH₄ emissions in the context of global climate projections.

8 Conclusions

This paper introduces WETMETH – a process-based wetland CH₄ model developed for implementation in ESMs. WETMETH is currently embedded in the UVic ESCM, a fully coupled EMIC. WETMETH is a computationally efficient model, applicable globally, and of appropriate complexity with respect to the current state of wetland CH₄ modelling. Unconstrained model parameters are tuned to potential CH₄ production rates from incubated soil samples and CH₄ emissions from northern wetlands due to the scarcity of large-scale datasets from other regions. Nevertheless, WETMETH skilfully reproduces estimates of mean annual CH₄ emissions over the past few decades at the regional, hemispheric, and global scale.

Despite the importance of tropical wetlands in the global CH₄ budget (Kirschke et al., 2013; Saunois et al., 2016) and climate change (O'Connor et al., 2010; Z. Zhang et al., 2017), their areal extent and associated CH₄ emissions remain highly uncertain in both the literature and modelling work (including this study) due to a combination of limited ground-based measurements and process understanding (Pangala et al., 2017; Saunois et al., 2020; Sjögersten et al., 2014) as well as a low accuracy from remotely sensed products, especially over the dense rainforests of Indonesia, Amazonia, and the Congo River basin where new peatlands continue to be discovered to date (Dargie et al., 2017). Large-scale wetland mapping is a field of ongoing research (Tootchi et al., 2019), and further model development should focus on the improvement of wetland simulations in the tropics. In parallel, a compilation of tropical wetland CH₄ measurements from various sources into synthesis datasets would be beneficial for constraining wetland CH₄ processes in large-scale models.

The inclusion of wetland CH₄ processes in a fully coupled ESM allows one to advance the research on the feedback between climate change and wetland CH₄ emissions. The implementation of WETMETH in the UVic ESCM constitutes an ideal tool for investigating interactions between climate conditions and wetland CH₄ emissions from decadal to longer timescales. Of particular importance is the permafrost carbon feedback to climate change, in which CH₄ emissions from northern wetlands are expected to play an important role (Nzotungicimpaye and Zickfeld, 2017).

Appendix A: Temperature-dependent Q₁₀ coefficient for methane production

Figure A1 illustrates the different shapes of the temperature-dependency function for CH₄ production $\left(Q_{\mathrm{10}}^{\frac{T_i-T_{\mathrm{0}}}{\mathrm{10}}};T_{\mathrm{0}}=\mathrm{273.15}\mathrm{K}\right)$ across a range of temperatures when considering (i) a constant Q₁₀ of 4.2 and (ii) a temperature-dependent Q₁₀ coefficient given by $Q_{\mathrm{10}}\left(T_i\right)=\mathrm{1.7}+\mathrm{2.5}\tanh \left[\mathrm{0.1}\right(T_{\mathrm{ref}}-T_i\left)\right]$, where T_ref=308.15 K. The temperature-dependent Q₁₀(T_i) implies an optimal temperature for CH₄ production in WETMETH of around 300.15 K. When Q₁₀(T_i) decreases to reach negative values, its value in WETMETH is set to 10⁻³ to represent a very small methanogenic response to temperature changes (Fig. A1).

Figure A1 (a) Differences between a constant Q₁₀ coefficient and a temperature-dependent Q₁₀(T_i) coefficient, and (b) implications for the temperature-dependency expression for CH₄ production ($Q_{\mathrm{10}}\left[\right(T_i-T_{\mathrm{0}})/\mathrm{10}]$). The temperature-dependent coefficient $Q_{\mathrm{10}}\left(T_i\right)=\mathrm{1.7}+\mathrm{2.5}\left[\tanh \right(\mathrm{0.1}(\mathrm{308.15}-T_i)\left)\right]$ allows one to account for uncertainties in the Q₁₀ coefficient and to define an optimal temperature for CH₄ production of around 300.15 K (dashed vertical line). The freezing point of water is shown at 273.15 K (continuous vertical line).

Appendix B: Applied minor modification to the TOPMODEL approach

The TOPMODEL approach implemented in the UVic ESCM is based on the formulation by Gedney and Cox for global land surface models (Gedney and Cox, 2003). This approach combines the simulated hydrology with a prescribed topographic index to determine the occurrence of wetlands (surface inundation) and soil moisture heterogeneity at the sub-grid scale. The occurrence of wetlands is simulated in an area whose local topographic index (Λ) satisfies the following condition:

$$\begin{array}{}\text{(B1)}& {\mathrm{\Lambda }}_{min}\le \mathrm{\Lambda }\le {\mathrm{\Lambda }}_{max},\end{array}$$

where Λ_min is a lower threshold that can be related to undersaturated conditions, and Λ_max is an upper threshold that can be related to oversaturated conditions.

In the initial work by Gedney and Cox, Λ_min depends on the transmissivity of the entire soil column (T(0)), the transmissivity of the soil column below the mean water table depth (z_w) of the grid box (T(z_w)), and the mean topographic index (Λ_mean). It is calculated as ${\mathrm{\Lambda }}_{min}=\ln \frac{T\left(\mathrm{0}\right)}{T\left(z_{\mathrm{w}}\right)}+{\mathrm{\Lambda }}_{\mathrm{mean}}$. While Λ_mean is static and prescribed with a topographic index map, both transmissivities (T(0) and T(z_w)) are simulated and non-static for a specific grid cell. Hence, Λ_min is a non-static and grid-dependent threshold. Unlike Λ_min, Λ_max is a static and global threshold. This threshold is applied to constrain the occurrence of wetlands in areas of stagnant water based on the assumption that locations where the water table rises well above the surface would be characterized by streamflow.

For the current study, a minor modification is applied to the above TOPMODEL approach. The revision consists of using a non-static and grid-dependent Λ_max instead of a static and global threshold. Following the formulation by Comyn-Platt et al. (2018), an expression for Λ_max that depends on Λ_min is currently used in the UVic ESCM. This threshold is defined as follows:

$$\begin{array}{}\text{(B2)}& {\mathrm{\Lambda }}_{max}={\mathrm{\Lambda }}_{min}+{\mathrm{\Lambda }}_{\mathrm{range}},\end{array}$$

where Λ_range is a global tuning parameter (Λ_range=0.93 in the version of the UVic ESCM used in this study).

In summary, unlike the initial work by Gedney and Cox (2003), the modified TOPMODEL approach considers two non-static and grid-dependent thresholds (Λ_min and Λ_max) for the identification of wetlands across the globe.

Appendix C: Unit conversion for potential methane production rates

Here, we describe steps followed for converting units of maximum CH₄ production rates measured in laboratory incubations from a soil weight basis (µg C g DW⁻¹ h⁻¹) to a soil volume basis (kg C m⁻³ s⁻¹). This unit conversion relies on the soil bulk density (BD in g cm⁻³) from the site of origin. The following two steps illustrate the applied unit conversion. In the first step, the potential CH₄ production rates (P_d,0) are converted from micrograms of carbon per gram dry weight per hour (µg C g DW⁻¹ h⁻¹) to micrograms of carbon per cubic centimetre per hour (µg C cm⁻³ h⁻¹) as follows:

$$\begin{array}{}\text{(C1)}& P_{\mathrm{d},\mathrm{1}}=\left(\mathrm{BD}\right)P_{\mathrm{d},\mathrm{0}}.\end{array}$$

Then, the conversion of P_d,1 from micrograms of carbon per cubic centimetre per hour (µg C cm⁻³ h⁻¹) to kilograms of carbon per cubic metre per second (kg C m⁻³ s⁻¹) is done as follows:

$$\begin{array}{}\text{(C2)}& P_{\mathrm{d},\mathrm{2}}={\displaystyle \frac{\mathit{\delta }}{\mathit{\gamma }}}{P}_{\mathrm{d},\mathrm{1}},\end{array}$$

where δ encompasses the conversion factors from micrograms (µg) to kilograms (kg) and from cubic centimetres (cm⁻³) to cubic metres (m⁻³) ($\mathit{\delta }={\mathrm{10}}^{-\mathrm{3}}$ kg m⁻³); and γ is the number of seconds per hour (γ=3600 s).