Evaluating alternative ebullition models for predicting peatland methane emission and its pathways via data–model fusion

Understanding the dynamics of peatland methane (CH₄) emissions and quantifying sources of uncertainty in estimating peatland CH₄ emissions are critical for mitigating climate change. The relative contributions of CH₄ emission pathways through ebullition, plant-mediated transport, and diffusion, together with their different transport rates and vulnerability to oxidation, determine the quantity of CH₄ to be oxidized before leaving the soil. Notwithstanding their importance, the relative contributions of the emission pathways are highly uncertain. In particular, the ebullition process is more uncertain and can lead to large uncertainties in modeled CH₄ emissions. To improve model simulations of CH₄ emission and its pathways, we evaluated two model structures: (1) the ebullition bubble growth volume threshold approach (EBG) and (2) the modified ebullition concentration threshold approach (ECT) using CH₄ flux and concentration data collected in a peatland in northern Minnesota, USA. When model parameters were constrained using observed CH₄ fluxes, the CH₄ emissions simulated by the EBG approach (RMSE = 0.53) had a better agreement with observations than the ECT approach (RMSE = 0.61). Further, the EBG approach simulated a smaller contribution from ebullition but more frequent ebullition events than the ECT approach. The EBG approach yielded greatly improved simulations of pore water CH₄ concentrations, especially in the deep soil layers, compared to the ECT approach. When constraining the EBG model with both CH₄ flux and concentration data in model–data fusion, uncertainty of the modeled CH₄ concentration profiles was reduced by 78 % to 86 % in comparison to constraints based on CH₄ flux data alone. The improved model capability was attributed to the well-constrained parameters regulating the CH₄ production and emission pathways. Our results suggest that the EBG modeling approach better characterizes CH₄ emission and underlying mechanisms. Moreover, to achieve the best model results both CH₄ flux and concentration data are required to constrain model parameterization.

1 Introduction

Methane (CH₄) emissions from wetlands constitute roughly one third of the global CH₄ budget (Denman et al., 2007; Saunois et al., 2020). Methane, following production under anoxic conditions, is stored belowground, oxidized into CO₂ by methanotrophs, or emitted into the atmosphere as CH₄. The emission of CH₄ is a major concern given its sustained-flux global warming potential (SGWP) of 45 (Neubauer and Megonigal, 2015). Methane emissions from wetlands cannot be simply estimated from production rates as more than 50 % of methane can be oxidized during transport to the atmosphere in various ecosystems (Conrad and Rothfuss, 1991; Teh et al., 2005; Segarra et al., 2015). The global wetland CH₄ oxidation sink has been estimated to be 40 %–70 % of CH₄ production, and it can dominate wetland CH₄ cycling (Megonigal et al., 2004). The oxidation sink depends substantially on the CH₄ emission pathways due to their different oxidation rates (Blodau, 2002).

Methane is produced at depth in the soil and then transported to the atmosphere via three primary pathways: ebullition, plant-mediated transport, and diffusion. Ebullition is the least vulnerable to oxidation as it allows CH₄ to quickly ascend in a bubble that bypasses the aerobic and anaerobic zones (Epstein and Plesset, 1950). Gases such as CH₄ can be released into the atmosphere by vascular plants (particularly sedges) after being transported through intercellular spaces (molecular diffusion) or aerenchymous tissues. Although plant-transported CH₄ bypasses aerobic zones of the soil, 20 %–90 % of plant-transported CH₄ can be oxidized in the rhizosphere or within the aerenchymous tissues where gaseous oxygen is present (Schipper and Reddy, 1996; Ström et al., 2005; Laanbroek 2010). Diffusive transport through the peat column is the slowest transport method, and therefore, CH₄ is most susceptible to oxidation as it spends the longest time transiting the aerobic and anaerobic zones (Chanton and Dacey, 1991; Megonigal et al., 2004). The relative importance of each pathway determines how much CH₄ is oxidized before it leaves the soil. Uncertainties in the relative contributions of these pathways to CH₄ emission can lead to large errors in the predictions of total CH₄ emissions (Bridgham et al., 2013). Despite their importance, the relative contributions of the CH₄ emission pathways have not been well quantified by either experimental or modeling approaches until recently (Ricciuto et al., 2021; Yuan et al., 2021).

Experimental data on the relative importance of CH₄ emission pathways are limited due to spatiotemporal heterogeneity and the difficulty in directly measuring the different pathways (Klapstein et al., 2014; Iwata et al., 2018). While most state-of-the-art land surface models (LSMs) incorporate CH₄ emission and differentiate the three transport pathways, information on the relative contribution of each pathway from modeling studies is still limited, and none of such studies has estimated the uncertainty or accuracy of the relative contributions of the emission pathways to net CH₄ emission (Bridgham et al., 2013). Comparisons between modeling approaches and empirical CH₄ data suggest that emission pathways may not be well captured by LSMs. For example, plant-mediated CH₄ transport by vascular species measured at northern peatlands accounted for 30 %–98 % of the total CH₄ emission (Shannon et al., 1996; Waddington et al., 1996), whereas model-estimated proportions in the similar ecosystems were all above 65 % (Tang et al., 2010; Wania et al., 2010). Empirical estimates also suggested that diffusion could range from 9 % to 60 % of the total CH₄ flux (Barber et al., 1988; Shea et al., 2010; Iwata et al., 2018). In contrast, modeled contributions from diffusion were always below 40 % (Tang et al., 2010; Wania et al., 2010; Peltola et al., 2018). More dramatically, modeling approaches estimated that ebullition constituted only 0 %–10 % of net CH₄ flux in natural vegetated wetlands (Tang et al., 2010; Wania et al., 2010; Peltola et al., 2018), much lower than the 10 %–64 % that was measured in experimental studies (Glaser et al., 2004; Tokida et al., 2007a, b). The uncertainties in simulated relative contributions of the pathways to net CH₄ emission in LSMs are mainly due to the lack of in situ information, inadequate representation of CH₄ processes, and unconstrained parameters used to describe emission pathways (Bridgham et al., 2013; Melton et al., 2013).

Since net CH₄ emissions depend on transport mode, all the emission pathways must first be represented correctly in ecosystem models in order to simulate CH₄ emission accurately (Blodau, 2002; Tang et al., 2010). Compared to diffusion and plant-mediated CH₄ transport, ebullition is less certain and could be the main reason for the mismatch between simulated and observed CH₄ concentrations in deep soil layers (Peltola et al., 2018). This is because diffusion is described with Fick's law and Henry's law, which have been widely used and well tested, and plant-mediated pathway happens only within the rooting depth, which is typically shallow in wetlands with a high water table level (Iversen et al., 2018). Ebullition makes a significant contribution to the total CH₄ emissions in some wetlands (Christensen et al., 2003; Yu et al., 2014). However, this process has not been incorporated well into most state-of-the-art LSMs. Mechanistically, CH₄ ebullition occurs when the buoyancy force of a bubble exceeds the retention force. During ascent, the bubbles exchange gas with the surrounding pore water, and some of the bubbles become trapped, allowing CH₄ to re-dissolve or be oxidized within the confining layer. In modeling studies, ebullition is commonly estimated using the ebullition concentration threshold (ECT) approach. In ECT, when the pore water CH₄ concentration is larger than a defined threshold, the excess CH₄ is directly released into the atmosphere (Walter and Heimann, 2000; Zhuang et al., 2004; Wania et al., 2010; Riley et al., 2011; Xu et al., 2016). However, this approach ignores the possibility of a CH₄ bubble moving into a less saturated layer where it can subsequently dissolve and possibly be oxidized, potentially overestimating ebullition. Other methods for modeling ebullition include the ebullition pressure threshold (EPT) or the ebullition bubble growth volume threshold (EBG) to trigger ebullition (Tang et al., 2010; Zhang et al., 2012). For the EPT method, bubbles form when the CH₄ concentration exceeds a certain threshold. The EBG method describes how temperature, pressure, and gas exchange alter the bubble volume and uses maximum bubble volume as a threshold to trigger ebullition events (Fechner-Levy and Hemond, 1996; Kellner et al., 2006; Zhang et al., 2012). Peltola et al. (2018) compared these modeling approaches and concluded that EBG should be incorporated into LSMs instead of ECT or EPT given its most realistic representation of the observed temporal variability in CH₄ emissions. However, the ability of the EBG approach to represent the relative importance of CH₄ emission pathways has not been evaluated against observations.

A more realistic projection of the emission pathways requires not only an improved model structure but also more appropriate parameter values (Wania et al., 2010; Riley et al., 2011; Shi et al., 2018). Data–model fusion directly informs process-based models by synthesizing multisource data streams and thus can help determine parameter values that lie within biophysically realistic ranges and reduce model uncertainty (Williams et al., 2009; Keenan et al., 2013; Shi et al., 2015a; Liang et al., 2018). Previous studies have found that sporadic measurements of net CH₄ emissions were only useful to constrain a few model parameters, and data assimilation with only CH₄ emission (flux-based) data did not help reduce the uncertainties in emission pathways (Bridgham et al., 2013; Ma et al., 2017). In our previous study, we found that monthly CH₄ emission data could only constrain CH₄ production-related parameters such as temperature sensitivity (Q₁₀) and basal production rate of CH₄ production (Ma et al., 2017). While direct measures of CH₄ emission pathways are rare, depth-specific pore water CH₄ concentration profiles can help elucidate the relative importance of CH₄ emission pathways. Indeed, measured CH₄ concentration profiles are critical for constraining the responsive parameters associated with CH₄ emission pathways because in process-based CH₄ models, all the three emission pathways are calculated based on the CH₄ concentration in each soil layer.

To date, few modeling studies have considered CH₄ concentration data for structural improvement or parameter optimization (Zhuang et al., 2004; Wania et al., 2010; Riley et al., 2011; Zhu et al., 2014). In those studies that compared simulation results to observed pore water CH₄ concentrations, the simulated concentration profiles did not agree well with observations despite good agreements between simulated and observed CH₄ emission data (Walter and Heimann, 2000; Tang et al., 2010). Thus, when CH₄ emission pathway parameters are calibrated using only net CH₄ flux data, models may not realistically represent CH₄ production, oxidation, and emission pathways. The exclusion of concentration profile data results in poorly constrained model parameters due to equifinality, in which multiple combinations of parameters result in similar flux predictions. This can cause misunderstanding of the mechanisms of CH₄ processes. It will be problematic to use these not-yet-well-calibrated parameter sets for climate change predictions or extrapolating CH₄ fluxes from the site level to larger spatial and temporal scales as these intermediate processes may have different responses to perturbations in climate.

To address these uncertainties, we evaluated the performance of two state-of-the-art methods for modeling ebullition, EBG and ECT, against the observed net CH₄ fluxes and pore water CH₄ concentration profiles in a northern Minnesota peatland. We also compared the strength of the flux-based data and pool-based data in constraining the parameters using data–model fusion. We hypothesized that (1) the EBG approach can reproduce the observed pore water CH₄ profiles better than the ECT approach given its more mechanical representations of bubble formation, gas exchange, and release and (2) pore water CH₄ concentration data offer more information for model parameters to reduce the uncertainties in simulated CH₄ emission and its pathways.

Table 1 The SPRUCE site data used in this study.

^∗ Total CH4 emissions were measured in August, September, and October 2011, May through November 2012, July, September, and October 2013, June through December 2014, April through November 2015, and March through December 2016.

2 Methods

2.1 Site and measurements

The data we used to calibrate our model were collected from the Spruce and Peatland Responses Under Changing Environments (SPRUCE) experiment, which is conducted in the 8.1 ha S1 bog in northern Minnesota in the USDA Forest Service Marcell Experimental Forest (47^∘30.476^′ N, 93^∘27.162^′ W) to study the responses of northern peatlands to climate warming and elevated atmospheric CO₂ concentration (Hanson et al., 2017a). The mean annual temperature from 1961 to 2009 at the SPRUCE site was 3.4 ^∘C, and the mean annual precipitation was 780 mm (Sebestyen et al., 2011). The mean peat depth is 2–3 m (Parsekian et al., 2012). The dominant plant species include Picea mariana, Larix laricina, a variety of ericaceous shrubs, and Sphagnum sp. moss. The graminoids Carex trisperma and Eriophorum vaginatum, as well as the forb Maianthemum trifolium, have seasonal dieback of their aboveground tissues in this peatland. Whole-ecosystem warming levels of +0, +2.25, +4.5, +6.75, and +9 ^∘C are paired with two CO₂ treatments (ambient or ∼ 400 ppm and 900 ppm) in open-top infrastructures (12 m × 8 m). Deep peat warming began in June 2014, aboveground warming began in August 2015, and elevated CO₂ treatments began in June 2016. In this study, however, all observed data we used were only from ambient plots (no infrastructures and no warming treatment) for our research goals, and we did not explore the warming effects on CH₄ processes. Modeling CH₄ emissions in response to warming and elevated CO₂ at this experiment site can be found in Yuan et al. (2021). A complete list of data streams used in this study is included in Table 1.

Environmental variables, including soil temperature, air temperature, relative humidity, wind speed, precipitation, and photosynthetically active radiation plots, were used as model input data. Measurements of environmental variables in ambient plots started in 2011. Soil temperature and moisture in different layers, water table depth (Hanson et al., 2015a, b, 2016b), carbon pools (leaf, wood, root, and peat soil; Hanson et al., 2018a, b; Norby and Childs, 2018; Norby et al., 2019), and community-scale fluxes, including gross primary production (GPP), net ecosystem exchange (NEE), ecosystem respiration (ER), and CH₄ flux data (Hanson et al., 2014; Hanson et al., 2016a, 2017b), were used to calibrate the modeled water–heat balance and carbon cycle similarly as in our earlier studies (Huang et al., 2017; Ma et al., 2017; Jiang et al., 2018).

CH₄ fluxes and pore water CH₄ concentrations were used for data assimilation. We averaged the data from all ambient plots measured on the same dates to represent the site-level CH₄ emissions and concentration profiles. Variations among different ambient plots were not considered in this study. In total, 45 daily CH₄ emission measurements were obtained from ambient plots from 2011 to 2016. In situ pore water CH₄ concentrations were measured monthly during the growing seasons in 2014–2016 (11 profiles in total) with the pore water samples collected from a series of piezometers permanently installed in the plots at 25, 50, 75, 100, 150, and 200 cm depths, respectively (Wilson et al., 2016). Piezometers consisted of a < 1 cm diameter pipe that limited diffusion. A total of 24 h prior to sampling, piezometers were pumped dry and allowed to recharge naturally so that the sampled water would not have been in prolonged contact with the atmosphere prior to collection. A perforated stainless-steel tube was inserted into the peat to collect samples within 0–25 cm depth. Samples were immediately filtered in the field through 0.7 µM Whatman glass-fiber filters and stored in pre-evacuated, septum-sealed glass vials. Phosphoric acid (1 mL, 20 %) was added to preserve each sample during shipment to Florida State University for analyzing CH₄ concentrations.

2.2 Model description

2.2.1 Overview of TECO_SPRUCE

For this study, we used the process-based biogeochemistry model TECO_SPRUCE (Terrestrial ECOsystem model at the SPRUCE site). The model was built with six major modules running at an hourly time step: canopy photosynthesis, soil water dynamics, plant growth, soil thermal dynamics, soil carbon–nitrogen (N) transfer, and soil CH₄ dynamics. A detailed description of these modules can be found in Weng and Luo (2008), Shi et al. (2015b), Huang et al. (2017), and Ma et al. (2017). Here we give a brief description of these modules but describe in detail how we calculated CH₄ ebullition with the EBG and ECT approaches.

The canopy photosynthesis module was mainly derived from a two-leaf model. It couples surface energy, water, and carbon fluxes. Leaf photosynthesis is estimated based on the Farquhar photosynthesis model (Farquhar et al., 1980) and the Ball and Berry stomatal conductance model (Ball et al., 1987). The soil water dynamic module has 10 soil layers and simulates water table level and soil moisture dynamics using rainfall, snowmelt, evapotranspiration, and runoff. Evaporative losses of water and associated latent heat are regulated by soil moisture in the first layer and atmospheric demand. Transpiration is determined by stomatal conductance and soil water content of the layers with roots present. When precipitation exceeds water recharge to soil water holding capacity, runoff occurs. The water table level is estimated using a simple bucket model as described by Granberg et al. (1999). The plant growth module calculates the allocation of photosynthesis carbon to individual plant pools (foliage, stem, and root), plant growth, plant respiration, phenology, and carbon transfer to the litter and soil carbon pools. Leaf onset is regulated by growing degree days (GDDs), and leaf senescence is determined by low-temperature and/or dry soil conditions. Phenology is represented by the seasonal variations in leaf area index (LAI) with LAI < 0.1 indicating the end of the growing season. The soil thermal dynamics module simulates snow cover, freezing depth, and soil temperature in 10 layers. The soil C–N transfer module simulates the movement of C and N from plants to two litter pools and three soil pools through litterfall, litter decomposition, and soil organic matter mineralization. Carbon fluxes from the litter and soil carbon pools are based on residence time and pool size of each C pool (Luo and Reynolds, 1999).

Figure 1 Conceptual structure of the CH₄ emission module in TECO_SPRUCE.

The CH₄ module simulates the transient, vertical dynamics of CH₄ production, oxidation, and belowground transport (via ebullition, plant-mediated transport, and diffusion), as well as CH₄ emissions at the soil surface–atmosphere interface (Fig. 1). The soil column is divided into 10 layers with each of the first five layers being 10 cm thick, whereas each of the rest of the layers are 20 cm thick. Within each soil layer, CH₄ concentration dynamics are calculated by a transient reaction equation with CH₄ production, CH₄ oxidation, released bubbles, plant-mediated transport, and the diffusion of CH₄ into and out of this soil layer from the lower and upper soil layer or the atmosphere for the first layer. Similar to CLM4Me (Riley et al., 2011), LPJ-WHyMe (Spahni et al., 2011; Wania et al., 2010), and TRIPLEX-GHG (Zhu et al., 2014) models, we assume that CH₄ production (P_ro) within the catotelm is directly related to heterotrophic respiration from soil and litter (R_h, $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$) via the following equation:

$$\begin{array}{}\text{(1)}& P_{\text{ro}}(z,t)=R_{\mathrm{h}}(z,t)f_{{\mathrm{CH}}_{\mathrm{4}}}{f}_{\text{stp}}\left(zt\right)f_{\text{pH}}{f}_{\text{red}},\end{array}$$

where $f_{{\mathrm{CH}}_{\mathrm{4}}}$ is an ecosystem-specific conversion scaler describing the fraction of anaerobically mineralized C atoms becoming CH₄. The parameters f_stp, f_pH, and f_red are environmental scalers, representing the effects of soil temperature, pH, and redox potential, respectively, on CH₄ production. Total emission of CH₄ from the soil to the atmosphere is calculated as the sum of CH₄ ebullition from saturated soil layers, plant-mediated CH₄ emissions from all the soil layers, and the diffused flux from the first soil layer into the atmosphere. A sensitivity test was done to decide which parameters need to be optimized by data–model fusion (Ma et al., 2017), and more detailed descriptions on the CH₄ module can be found in Ma et al. (2017).

The original method used in TECO_SPRUCE for the ebullition process was the constant concentration threshold method (Walter and Heimann, 2000; Ma et al., 2017). However, a number of factors such as atmospheric pressure, water table level, and temperature have been shown to affect ebullition (Beckmann et al., 2004; Kellner et al., 2006; Tokida et al., 2007a). Here, we used two new methods, i.e., the modified concentration threshold method (ECT) and the bubble growth volume threshold method (EBG), to describe CH₄ ebullition. In both methods, direct ebullition into the atmosphere can take place only when the water table level is at or above the soil surface; otherwise, CH₄ in bubbles is added to the soil layer just above the water table and then continues to diffuse through the soil layers to the atmosphere. Below we describe these two methods in detail.

2.2.2 Ebullition approach based on the concentration threshold (TECO_SPRUCE_ECT)

With the concentration threshold approach, we assume that bubbles form when the CH₄ concentration exceeds a certain threshold based on the equilibrium concentration defined by Henry's law. Instead of using a constant value for the threshold, in this study, we allowed the threshold to fluctuate with atmospheric pressure, water column pressure, and soil temperature, following the method proposed by Wania et al. (2010). The maximum solubility of CH₄ at a given temperature was calculated using a statistical model used by Yamamoto et al. (1976) based on the empirical data:

$$\begin{array}{}\text{(2)}& V=\mathrm{0.05708}-\mathrm{0.001545}T+\mathrm{0.00002069}{T}^{\mathrm{2}},\end{array}$$

where V is the Bunsen solubility coefficient, defined as the volume of gas dissolved per volume of water at atmospheric pressure and a given temperature. The volume of CH₄ dissolved per volume of water was converted into grams using the ideal gas law:

$$\begin{array}{}\text{(3)}& {\left[{\mathrm{CH}}_{\mathrm{4}}\right]}_{\text{thre}}=PV\cdot C/RT,\end{array}$$

where ${\left[{\mathrm{CH}}_{\mathrm{4}}\right]}_{\text{thre}}$ is the maximum concentration threshold (g C m⁻³), P is the sum of the atmospheric and hydrostatic pressures (Pa), V is the Bunsen solubility coefficient as in Eq. (2), the constant C is the atomic weight of carbon (12 g mol⁻¹), the gas constant R is 8.3145 ${\mathrm{m}}^{\mathrm{3}}\mathrm{Pa}{\mathrm{K}}^{-\mathrm{1}}{\mathrm{mol}}^{-\mathrm{1}}$, and T is the temperature (K). Then the CH₄ ebullition flux can be calculated using the following equation:

$$\begin{array}{}\text{(4)}& E_{\text{bu}}(z,t)=\left\{\begin{array}{l}{K}_{\text{ebu}}\left(\left[{\mathrm{CH}}_{\mathrm{4}}\right](z,t)-{\left[{\mathrm{CH}}_{\mathrm{4}}\right]}_{\text{thre}}\right)\\ \text{if}\left[{\mathrm{CH}}_{\mathrm{4}}\right]>{\left[{\mathrm{CH}}_{\mathrm{4}}\right]}_{\text{thre}}\\ \mathrm{0.0}\\ \text{if}\left[{\mathrm{CH}}_{\mathrm{4}}\right]\le {\left[{\mathrm{CH}}_{\mathrm{4}}\right]}_{\text{thre}}\end{array}\right.,\end{array}$$

where E_bu(z,t) is the ebullition flux of CH₄ (g C h⁻¹) to the lowest air layer, K_ebu is a rate constant of 1.0 h⁻¹ (Walter and Heimann, 2000; Zhuang et al., 2004, 2006), and [CH₄](z,t) is the pore water CH₄ concentration in soil depth z at model time step t.

2.2.3 Ebullition approach based on the bubble growth volume threshold (TECO_SPRUCE _EBG)

In contrast to the concentration threshold approach, the EBG approach uses bubble volume as a threshold to trigger ebullition events (Fechner-Levy and Hemond 1996), and it has been applied to model CH₄ ebullition (Kellner et al., 2006; Zhang et al., 2012; Peltola et al., 2018). The total bubble volume in each soil layer is calculated and updated continuously based on the ideal gas law and Henry's law. In detail, if CH₄ concentration exceeds the limit that the water can withhold based on Henry's law, then excess CH₄ is converted to a gaseous volume calculated using the predefined bubble CH₄ mixing ratio (f). This gaseous volume is divided evenly into a certain number of bubbles (N_bub). N_bub is a unitless tuning parameter ranging between 5 and 500 in each 10 cm thick soil layer and 10 and 1000 in each 20 cm thick soil layer, which essentially controls the mass exchange rate between the gas volume and the pore water. The CH₄ exchange between the stationary bubbles and the pore water (Q_ebu) is calculated using the equation proposed by Epstein and Plesset (1950):

$$\begin{array}{}\text{(5)}& Q_{\text{ebu}}={\displaystyle \frac{\mathrm{4}\mathit{\pi }r{D}_{\mathrm{w}}{N}_{\text{bub}}}{V_{\mathrm{w}}}}\left(c_{\mathrm{w}}-{\displaystyle \frac{{\mathrm{H}}^{\mathrm{cc}}fP}{RT}}\right),\end{array}$$

where r is the radius of a bubble (m), D_w is the CH₄ diffusion coefficient in water (m² s⁻¹) calculated using the quadratic curve of observed diffusivities against temperatures (Broecker and Peng, 1974), V_w is the amount of water in this layer (m³), c_w is dissolved CH₄ concentration in the pore water, and H^cc is the dimensionless Henry solubility of CH₄ calculated following Sander (1999). P, R, and T are the same as in Eq. (3). A negative value of Q_ebu indicates CH₄ transfer from the bubbles back to the pore water. This reverse gas exchange mechanism has not been included in other ebullition methods but has been revealed as an important process in empirical studies (McGinnis et al., 2006; Rosenberry et al., 2006). The ebullition flux E_bu(z,t) is then calculated when the bubble volume at a certain depth (z) exceeds the volume threshold (V_max) within the time step t:

$$\begin{array}{}\text{(6)}& {\displaystyle }{V}_{max}=V_{\text{maxfraction}}\cdot V_{\mathrm{w}},\text{(7)}& {\displaystyle }{E}_{\text{bu}}(z,t)=\left\{\begin{array}{l}{c}_{\mathrm{b}}\cdot \left(V_{\mathrm{B}}+\mathrm{\Delta }{V}_{\mathrm{B}}-V_{max}\right),\\ \text{if}V+\mathrm{\Delta }V>V_{max}\\ \mathrm{0.0}\\ \text{if}V+\mathrm{\Delta }V<V_{max}\end{array}\right.,\end{array}$$

where V_maxfraction is the free-phase gas-filled fraction of the pore space in the soil layer above which ebullition occurs, c_b is the CH₄ concentration in a bubble (mol m⁻³), V_B is the total volume of all bubbles, and ΔV_B is the change in the total volume due to the diffusive gas exchange in Eq. (5). The amount of CH₄ in all bubbles after each time step is

$$\begin{array}{}\text{(8)}& n_{\mathrm{b}}={\displaystyle \frac{fP{V}_{\mathrm{B}}^{\prime }}{RT}},\end{array}$$

where f is the predefined bubble CH₄ mixing ratio as mentioned earlier, and $V_{\mathrm{B}}^{\prime }$ is the updated total bubble volume after each time step. Excess bubbles will be released into the lowest air layer within one time step unless they are trapped in the soil profile. To determine if a bubble will be trapped, we adopted an approach similar to Peltola et al. (2018), assuming that the probability for a bubble to be trapped within a certain soil layer is a predefined constant number (bubprob); thus the bubbles formed in deeper layers would have a larger probability of being trapped during ascent. In contrast, all other ebullition modeling methods assume that no bubbles will get trapped.

Table 2 Parameters used for data–model fusion.

The values of bubprob, f, and V_maxfraction are dependent on the soil texture, porosity, water content, etc. and have been found to significantly affect the modeled CH₄ fluxes, the layers where bubbles were formed, and the number of ebullition events (Zhang et al., 2012; Peltola et al., 2018). The tuning parameter, N_bub, however, has a minimal effect on modeled ebullition (Peltola et al., 2018). In this study, we used the empirical values measured from other sites or the values used in other models as the prior ranges of bubprob, f, and V_maxfraction in our models (Table 2). Then we constrained these parameter values via data–model techniques so that the model estimation of the ebullition process was more accurate.

2.3 Data–model fusion

We used the Markov chain Monte Carlo (MCMC) method based on the Metropolis–Hasting algorithm (Metropolis et al., 1953) to optimize the posterior distribution of parameters and explore model uncertainty. Both the observed data and simulated results were rescaled to a daily emission unit for comparison. The prior range for each parameter was assumed to be uniformly distributed, which indicates that all values within the range have equal likelihood. We also assumed that errors between observations and model simulations independently follow a normal distribution with a zero mean. The cost function weights the mismatch between observations and the modeled corresponding variables, represented by

$$\begin{array}{}\text{(9)}& p\left(Z\mathrm{|}\mathit{\theta }\right)\propto \exp \left\{-{\sum }_{i=\mathrm{1}}^{\mathrm{2}}{\sum }_{t\in Z_i}{\displaystyle \frac{\left[Z_i\right(t)-X(t){]}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}_i^{\mathrm{2}}\left(t\right)}}\right\},\end{array}$$

where Z_i(t) is the ith observation stream at time t, X(t) is the model simulated value, and σ_i(t) is the standard deviation of observation error estimates.

The parameter space was explored for 50 000 iterations during the optimization process. The new parameter value at the current step was based on the accepted parameter in the previous step by a proposed distribution. The current value was accepted if the observation–model difference was reduced or otherwise passed with a random probability. We ran five chains of 50 000 simulations and used the Gelman–Rubin statistic (Gelman and Rubin, 1992) to check the convergence of sampling chains. The first half of the accepted parameters were discarded as the burn-in period, and the second half were used for posterior analysis. More details on sampling and the cost function can be found in Xu et al. (2006).

Table 3 Details for data assimilation runs.

Parameters directly regulating CH₄ emission pathway and belowground dynamics and their prior ranges used for data assimilation are listed in Table 2. Specifically, we selected four parameters (i.e., $f_{{\mathrm{CH}}_{\mathrm{4}}}$, Q_10pro, O_max, and T_veg) from the TECO_SPRUCE_ECT and seven parameters (all the seven parameters in Table 2) from the TECO_SPRUCE_EBG during data assimilation. The prior ranges were determined by combining information from empirical measurements or modeling studies from peatland ecosystems. The in situ CH₄ emission and pore water CH₄ concentration data from ambient plots (Table 1) were used as observations to constrain model parameters. In order to evaluate how a proper model structure and constrained parameter values help improve model-simulated CH₄ emission pathways, we conducted four data assimilation runs with the TECO_SRUCE model, as shown in Table 3.

We illustrate the improvement from model structure by comparing ECT_F and EBG_F, which were calibrated using the observed CH₄ flux data. Then we compare results from EBG_F and EBG_FC to show the ability of pore water CH₄ concentration data to help constrain the parameters related to the CH₄ emission pathways. Model performance was evaluated against the observed data using root mean square error (RMSE). Model uncertainties in pore water CH₄ concentrations were quantified as the standard deviation across all soil layers in each of the model runs listed in Table 3.

Table 4 Parameter values for the posterior distribution of parameters.

^∗ MLE: maximum likelihood estimation.

Figure 2 Posterior distribution of parameters that govern methane emission processes from data–model fusion. Parameters are defined in Table 2. Panels (a) and (b) are parameters related to methane production, (c) is a parameter related to methane oxidation, (d) is a parameter related to plant transport, and (e–g) are parameters related to ebullition. X axes indicate the prior ranges of parameters used for data–model fusion. The blue lines are the parameter posterior distributions (PPDs) from the ECT model structure trained with CH₄ emission data (ECT_F). The green lines are the PPDs from the EBG model structure trained with CH₄ emission data (EBG_F). The purple lines are the PPDs from EBG trained with both CH₄ emission and CH₄ concentration data (EBG_FC). Well-constrained parameters have a unimodal distribution, whereas poorly constrained parameter distributions tend to be flat.

3 Results

3.1 Parameter optimization using CH₄ flux data with different model structures

The CH₄ production-related parameters, $f_{{\mathrm{CH}}_{\mathrm{4}}}$ (fraction of anaerobically mineralized C atoms becoming CH₄) and Q_10pro (temperature sensitivity of CH₄ production), were constrained using the CH₄ emission flux data for both TECO_SPRUCE_ECT and TECO_SPRUCE_EBG (Table 4, Fig. 2a and b). However, the maximum likelihood estimates (MLEs) for parameters varied between the two models (Table 4), with $f_{{\mathrm{CH}}_{\mathrm{4}}}$ slightly increased from 0.16 to 0.17 and Q_10pro decreased from 3.0 to 2.69 in the EBG approach compared to the ECT approach. In contrast, O_max and T_veg were not constrained in either model with large uncertainties in model-estimated CH₄ oxidation and plant transport parameters (Table 4, Fig. 2c and d).

Of the three ebullition-related parameters used only in the EBG approach, when assimilating only the CH₄ emission flux data, V_maxfraction (maximum fraction of volume occupied by bubbles) was constrained with a unimodal shaped posterior distribution (Fig. 2g), f(CH₄ mixing ratio in bubbles) was edge hitting with a marginal distribution downward (Fig. 2e), and bubprob (probability that a bubble will get trapped at a certain soil layer) was constrained with a wide, slightly domed distribution (Fig. 2f).

Figure 3 Observed versus modeled methane emissions (a, c, and e) and simulated relative contributions (%) of plant-mediated transportation (PMT), diffusion, and ebullition under ambient conditions (b, d, and f). Black dots are observed CH₄ emissions from static chamber measurements. Error bars are the standard errors. Blue lines are ECT-model-simulated CH₄ emissions based on the parameter probability distributions (PPDs) constrained by CH₄ flux data (ECT_F, RMSE = 0.61). Green lines are EBG-model-simulated CH₄ emissions based on the PPDs constrained by CH₄ flux data (EBG_F, RMSE = 0.53). Purple lines are EBG-model-simulated emissions based on the PPDs constrained by both CH₄ flux and concentration data (EBG_FC, RMSE = 0.52). The mid-lines and shaded areas correspond to the means and standard deviations, respectively, from 500 model simulations with parameters randomly drawn from the posterior distributions. Relative contributions (%) are the daily mean values calculated from the simulations.

3.2 Evaluations of model structures against the observed data

Using CH₄ emission data to constrain the parameters, EBG-simulated CH₄ emissions (RMSE = 0.53, Fig. 3c) had a better agreement with observations than ECT (RMSE = 0.61, Fig. 3a). In addition, EBG simulated a smaller seasonal variability in CH₄ emissions (Fig. 3c) than ECT (Fig. 3a). The simulated contributions from plant-mediated transport, diffusion, and ebullition were 40.7 ± 8.0 %, 35.7 ± 8.7 %, and 23.5 ± 9.4 %, respectively, in ECT_F (Fig. 3b) and 38.4 ± 13.9 %, 38.7 ± 9.9 %, and 22.7 ± 9.4 %, respectively, in EBG_F (Fig. 3d). Compared to ECT (Fig. 3b), EBG simulated a smaller contribution from ebullition but more frequent ebullition events (Fig. 3d).

Figure 4 Observed versus simulated pore water CH₄ concentration profiles. Black dots are observed concentrations measured from piezometer samples. Blue lines are the ECT-model-simulated concentrations based on the parameter probability distributions (PPDs) constrained by CH₄ flux data (ECT_F). Green lines are the EBG-model-simulated concentrations based on the PPDs constrained by CH₄ flux data (EBG_F). Purple lines are the EBG-model-simulated concentrations based on the PPDs constrained by both CH₄ flux and concentration data (EBG_FC). All mid-lines and shaded areas correspond to the means and standard deviations, respectively, from 500 model simulations with parameters randomly drawn from the posterior distributions.

The ECT model constrained by CH₄ flux data could not reproduce well the patterns of the observed pore water CH₄ concentrations, especially in the deep peat layers (RMSE = 0.77, Fig. 4, ECT_F). When calibrated by CH₄ flux data alone, the EBG approach for ebullition captured deep layer CH₄ concentrations much better than the ECT approach (RMSE = 0.33, Fig. 4, EBG_F). The observed concentration profiles lay within the 95 % probability intervals and the means were comparable to observed patterns. However, the EBG model structure simulated a relatively large range of CH₄ concentration profiles, especially in the deep peat layers, mainly due to the unconstrained T_veg and bubprob controlling the plant transport and ebullition pathways, respectively (Fig. 2, EBG_F).

3.3 Comparison of the flux- and pool-based data in constraining the parameters for simulating CH₄ processes

For the ECT approach, as described earlier, assimilating the observed CH₄ flux could constrain parameters such as $f_{{\mathrm{CH}}_{\mathrm{4}}}$ and Q_10pro. However, when using both observed CH₄ flux and concentration data to constrain the parameters of TECO_SPRUCE_ECT (i.e., the ECT_FC run), no parameter set was accepted within the observational uncertainty range, indicating that the ECT model structure failed to simultaneously simulate the dynamics of CH₄ emissions and pore water CH₄ concentrations.

In contrast to the ECT approach, incorporation of pore water CH₄ concentration data in the EBG approach greatly improved parameter estimations. While T_veg and bubprob were not constrained by flux-based observation data alone (Table 4, Fig. 2, EBG_F), they were well constrained to a unimodal distribution when both CH₄ flux and pore water CH₄ concentration data streams were used in the data–model fusion (Table 4, Fig. 2, EBG_FC). Compared to EBG_F, the parameter T_veg was well constrained to a very small range of 1.43 ± 0.46, and the parameter bubprob was also well constrained to a range of 0.25 ± 0.015 with less uncertainty under EBG_FC (Table 4, Fig. 2). The parameter $f_{{\mathrm{CH}}_{\mathrm{4}}}$ decreased from 0.17 in EBG_F to 0.15 in EBG_FC, whereas Q_10pro increased from 2.69 in EBG_F to 3.21 in EBG_FC (Table 4, Fig. 2). Moreover, the formerly constrained range of parameter f under EBG_F shifted from 0.11 ± 4.0 to 0.29 ± 0.46 when the pore water CH₄ concentration information was added into data assimilation. All the emission pathway-related parameters (T_veg, bubprob, f, and V_maxfraction) were well constrained once the pore water CH₄ concentration profile information was added to data–model fusion. However, incorporation of the pore water CH₄ concentration data in data assimilation with the TECO_SPRUCE_EBG did not improve the constraint of O_max.

In terms of the model's performance in fitting observed CH₄ emission patterns, the two parameterization methods for the EBG approach were comparable, with RMSE of 0.53 under EBG_F and RMSE of 0.52 under EBG_FC (Fig. 3c and e). However, the simulated contributions from plant-mediated transport, diffusion, and ebullition by EBG_FC, which were 31.8 ± 4.9 %, 58.1 ± 5.1 %, and 9.9 ± 6.1 %, respectively (Fig. 3f), varied greatly from those simulated by EBG_F (Fig. 3d). The contribution from ebullition under EBG_FC was much less than that under EBG_F (Fig. 3d). CH₄ flux and concentration data together reduced the uncertainty of the modeled CH₄ concentration profiles by 78 %–86 % compared to the flux data alone for data–model fusion, with RMSE reducing from 0.33 in EBG_F to 0.12 in EBG_FC (Fig. 4). The uncertainty in modeled CH₄ concentration profiles was decreased mainly due to the well-constrained parameters regulating the CH₄ production and emission pathways (Fig. 2a, b, and d–g).

4 Discussion

In this study, we evaluated two alternative model structures with two data streams, i.e., CH₄ emission and pore water CH₄ concentration data, in simulating peatland CH₄ emission and its pathways.

4.1 Better representing CH₄ emission and pore water CH₄ concentrations by the ebullition bubble growth (EBG) model

Previous studies suggested that the EBG method of modeling ebullition agreed better with the observed temporal variability in CH₄ emissions (R²= 0.63) when compared with the ECT (R²= 0.56) and EPT (R²= 0.35) methods (Peltola et al., 2018). We also found that the EBG-simulated CH₄ emissions (RMSE = 0.53) had a better agreement with observations than the ECT method (RMSE = 0.61). Ebullition events simulated by EBG had a higher frequency but a smaller magnitude than those obtained from ECT, which is consistent with on-site minirhizotron observations of small bubbles around fine roots (Fig. S1 in the Supplement). Although the ECT method was able to simulate a similar seasonal pattern of CH₄ emissions as EBG, the mean annual CH₄ emission was 17.8 % lower compared with the EBG method. Peltola et al. (2018) reported that the different ebullition modeling approaches simulated significantly different amount of CH₄ stored belowground and distinctly different distributions of CH₄ along the soil profiles. In line with their results, we found that the ECT method produced much higher pore water CH₄ concentrations than the EBG method, especially in the deep layers (Fig. 4).

Of the few modeling studies that compared results with observed belowground CH₄ concentration, Walter et al. (2001) simulated CH₄ concentration with an early generation ECT method. This method used a constant concentration threshold that was tuned to match the observed concentration data, but they also found discrepancies with observed data (only CH₄ concentrations within the first 50 cm soil were compared). Tang et al. (2010) compared the EPT method with the early generation ECT method and found that EPT had an improved CH₄ concentration profile, although a mismatch in the concentration profile remained, especially from the model that best reproduced observed CH₄ emissions. The recently developed land surface model (LSM) with a new microbial-functional group-based CH₄ module incorporated, i.e., the ELM_SRUCE model, used the modified ECT method for the ebullition process but incorporated the acetoclastic and hydrogenotrophic pathways for methanogenesis, as well as anaerobic and aerobic oxidations (Ricciuto et al., 2021). This model could accurately predict the seasonal cycle of CH₄ production and net fluxes, but CH₄ concentrations in soil layers deeper than 1 m were still not well simulated (Ricciuto et al., 2021) and led to different estimates of emission pathways (23.5 % for plant-mediated transportation (PMT), 15.0 % for diffusion, and 61.5 % for ebullition) from our study (31.8 ± 4.9 % for PMT, 58.1 ± 5.1 % for diffusion, and 9.9 ± 6.1 % for ebullition). In our study, when training the modified ECT model with both CH₄ emission and pore water concentration data, no parameter set was accepted, which suggested that the ECT method was not able to simultaneously reproduce both the magnitude of observed CH₄ fluxes and the patterns of pore water CH₄ concentrations no matter the combinations of parameters used. In contrast, the EBG method could capture observed CH₄ emissions and the patterns of pore water CH₄ concentration profiles simultaneously (Figs. 3 and 4).

Moreover, we found that although both the ECT and EBG methods were able to represent the general patterns of observed CH₄ emissions, the flux-constrained parameter distributions varied between the two methods. For example, $f_{{\mathrm{CH}}_{\mathrm{4}}}$ increased, but Q_10pro decreased in EBG compared to ECT (Table 4, Fig. 2), which might be attributed to the confounding effects of missing/inappropriate model structures on parameter estimation because different combinations of model parameter values or structures can give similar model outputs (Williams et al., 2009). More studies are needed to further explore model structures and parameter optimization methods to best simulate CH₄ production and emission processes and the underlying mechanisms.

4.2 Pool-based CH₄ concentration data reduced the uncertainty of the emission pathways

Our study suggests that even using a more reasonable model structure, i.e., EBG, parameter sets that resulted in good simulations of CH₄ emissions did not necessarily reproduce a realistic pore water concentration profile (Figs. 3 and 4). By comparing the parameter posterior distributions trained by observed CH₄ emissions with and without observed pore water concentration profiles using the same model structure, TECO_SPRUCE_EBG, we revealed that CH₄ emission data could constrain the CH₄ production-related parameters $f_{{\mathrm{CH}}_{\mathrm{4}}}$ and Q_10pro and ebullition-related parameter V_maxfraction very well. The ebullition-related parameter f was edge hitting, but parameter bubprob and plant transport-related parameter T_veg remained unconstrained, causing large uncertainty in simulated ebullition and plant-mediated transport (Table 4). However, by training the model with both CH₄ emission and pore water CH₄ concentration data, the parameters regulating CH₄ production, plant transport, and ebullition were all well constrained (Fig. 2). This is because the vertical concentration profile of CH₄ is a balance between the dynamic CH₄ production, oxidation, and three emission pathways. The constrained parameters contributed to a more accurate estimation of pore water CH₄ concentration (RMSE = 0.12) and better-constrained emission pathways (Table 4, Fig. 4).

Previous studies have emphasized the importance of combining carbon-pool data with carbon-flux data to improve estimated ecosystem carbon exchange. For example, Richardson et al. (2010) reported the initial leaf pool size could not be constrained until biomass information was combined with flux data. Du et al. (2015) also found that carbon flux data could constrain parameters reflecting instant responses to environmental changes such as temperature sensitivity, while pool-based data mainly contained information that could help constrain transfer coefficients. GPP and ER data could effectively constrain parameters that were directly related to flux data, such as the temperature sensitivity of heterotrophic respiration, the carbon allocation to leaves, and leaf turnover rate (Fox et al., 2009). In our study, the CH₄ emission data mainly constrained parameters that represented instant responses to temperature change (Q_10pro) and input rate from the source pool ($f_{{\mathrm{CH}}_{\mathrm{4}}}$). The pore water CH₄ concentration data contributed to constraining the allocation rates of CH₄ to the different emission pathways. Due to the different information contained between CH₄ flux and concentration data, we highly recommend that both types of measurements should be made when possible and that both data streams should be used when constraining CH₄ models.

It needs to be noted that there is a large disagreement in simulated relative contribution by ebullition between CH₄ flux-data-constrained models (i.e., 0.13 % by the ECT approach with a constant concentration threshold (Ma et al., 2017), 23.5 % by the modified ECT approach with varied concentration thresholds in our study, and 22.7 % by the EBG approach in our study, as well as the EBG approach constrained with both CH₄ flux and concentration data (9.9 %) (Figure 3). This suggests the urgent need of observed data for separating these relative contributions in field experiments, possibly through (1) having continuous total emission flux measurement (Susiluoto et al., 2018) despite being hard to deploy and calibrate in the field, (2) separately measuring diffusive, plant-mediated-transport, and ebullition fluxes despite being technically challenging, and (3) measuring belowground CH₄ concentration profile as suggested in our study. At the SPRUCE experiment site, starting in the summer of 2022, two auto-chambers with a footprint of 0.2 m² will be deployed in each plot to measure CO₂ and CH₄ fluxes. Along with the continued CH₄ profile measurement, these whole set of observations will provide the opportunity to further evaluate these discussed approaches for improving model simulations.

4.3 Broader impacts and implications

Large uncertainties exist in understanding future wetland CH₄ emissions in response to projected climate change, which result from inappropriate model structure and insufficient parameterizations even after the uncertainties in wetland areas are considered (Melton et al., 2013; Luo and Schuur, 2020). Decades of modeling research on CH₄ has evolved to a stage that the emission pathways are explicitly calculated with various complexities, but determining the accuracy and uncertainty of individual pathways still requires more research (Xu et al., 2016). Currently, models fail to reproduce diffusive fluxes by more than 40 % mainly due to the lack of validations by the pore water CH₄ concentrations (Tang et al., 2010; Riley et al., 2011). In LSMs, plant transport exclusively dominates CH₄ emission in all wetland types tested (Tang et al., 2010; Wania et al., 2010; Peltola et al., 2018). However, according to the experimental studies, each of the three emission pathways can dominate, depending on the wetland type, vascular plant coverage, and the height of the water table (Whiting and Chanton, 1992, Shannon et al., 1996). By assimilating empirical data of both CH₄ flux and pore water CH₄ concentration data, our data–model fusion study proposes a more reasonable model structure and more robust parameter estimates with greatly reduced uncertainties.

Our results also implicate barriers of current CH₄ modeling studies and suggest future directions for both modeling and experimental efforts, namely (1) the under-described CH₄ processes in models and (2) the lack of observational data to constrain key CH₄ processes in the models. More explicit CH₄ processes are needed in modeling CH₄ emission and its pathways. For example, in this study, the maximum aerobic oxidation rate (O_max) was always poorly constrained with wide, slightly domed distributions (Fig. 2c) regardless of what observation data were being assimilated into the models. This poor constraint might partly result from the missing anaerobic oxidation process in the models. In current process-based models, much of the descriptions of CH₄ dynamics in wetland soils are based on the premise that the oxidation of CH₄ occurs only in aerobic environments (Wania et al., 2010; Riley et al., 2011; Bridgham et al., 2013). However, the anaerobic oxidation of CH₄ may be an important sink for CH₄, sometimes reducing emissions by over 50 % in experimental studies (Smemo and Yavitt, 2011; Gupta et al., 2013; Segarra et al., 2015). Recently, a microbial-functional-group-based CH₄ model was developed accounting for both aerobic and anaerobic CH₄ oxidations, and this model has been validated against the concentration of CH₄ and CO₂ from incubation data (Xu et al., 2015). In Xu et al. model, 7 out of total 33 key CH₄ process parameters control CH₄ oxidation, and their values vary widely across different ecosystems and environmental conditions. The incorporation of anaerobic CH₄ oxidation into LSMs may help improve the calculations of CH₄ oxidation if the uncertainties from these CH₄-oxidation-related parameters can be reduced.

While more comprehensive and process-based models for simulating all the processes or mechanisms involved in CH₄ emissions are laudable, observations on such specific processes are critical to constrain parameters and reduce model uncertainty. Without sufficient data to evaluate such processes or to calibrate models, developing such complex models to explicitly simulate these processes could also introduce large uncertainties. Increased model complexity only contributes to the improved forecasting if parameters can be calibrated adequately by observed data (Famiglietti et al., 2021). If there were not enough observational data for model calibrations, increased complexity can lead to even worse forecast skills than the intermediate-complexity models (Shi et al., 2018; Famiglietti et al., 2021). Currently, similar to our model, many process-based biogeochemistry models (e.g., CLM, LPJ, TRPLEX, JULES, TEM) also use a parameter that varies with soil conditions to describe the potential ratio of CO₂ becoming CH₄ ($f_{{\mathrm{CH}}_{\mathrm{4}}}$), which is due partly to the limitation of data availability. Another example of observational data hindering model development is the unconstrained parameters to calculate plant-mediated transport. Although new algorithms and parameters to calculate plant aerenchyma transport have been added to LSMs to represent this mechanism more realistically, the parameters such as tiller radius, number of tillers, cross-section area of tillers, and the tiller porosity are highly idealized and poorly constrained (Wania et al., 2010; Riley et al., 2011). In the TECO_SPRUCE model used in this study, the parameter T_veg was used as a proxy of the ability of the whole plant community (e.g., biomass and abundance) to emit CH₄. Root growth was simulated by a phenological process using LAI and temperature, and in situ fine root profile measurements were used as a proxy for vertical rooting distributions (Iversen et al., 2018). T_veg was well constrained by the observed CH₄ emission and concentration data at a range of 1.43 ± 0.46, which indicated that the ability of the plant community to emit CH₄ at this site was low (compared to its prior knowledge of 0.01–15, Table 2). Empirical measurements of plant-mediated CH₄ transport at the same study site supported our model results (Scott Bridgham, personal communications). This finding can also be explained given that the diversity and abundance of aerenchymous plants at our study site were low, similar to many other northern ombrotrophic bogs.

5 Conclusions

Understanding relative contributions of CH₄ emission pathways is critical to mechanistically model future CH₄ dynamics. Acknowledging that pore water CH₄ concentration is the driving force for each emission pathway, we evaluated the ability of two ebullition modeling approaches to reproduce observed CH₄ emissions and pore water concentration profiles at a large-scale manipulated experimental site in a northern Minnesota, USA, peatland. The ebullition bubble growth volume threshold approach (EBG) fits the observed CH₄ emissions and CH₄ concentration profiles much better than the modified ebullition concentration threshold approach (ECT), especially for CH₄ concentrations in the deeper soil layers. By assimilating the net CH₄ emission and belowground CH₄ concentration data into the models, we substantially reduced the uncertainties of modeled CH₄ emissions from the emission pathways involved. While net CH₄ efflux data are often the only data stream for CH₄ model validations, we recommend that more attention be given to in situ measurements of the pore water CH₄ concentrations and assimilations of the concentration data for model parameterization. Since the relative ratio of the emission pathways (ebullition, plant-mediated transport, and diffusion) determines how much CH₄ is oxidized before leaving the soil due to their different transport rates and vulnerability to oxidation, we also suggest that the EBG approach should be incorporated into land surface models (LSMs) so that the projections of both CH₄ emission and its transport processes are more realistic in response to climate change scenarios. Future studies should also include anaerobic CH₄ oxidation into LSMs and constrain its parameters to better predict wetland CH₄ emissions.