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

Abstract

Understanding the dynamics of peatland methane (${\mathrm{CH}}_{\mathrm{4}}$) emissions and quantifying sources of uncertainty in estimating peatland ${\mathrm{CH}}_{\mathrm{4}}$ emissions are critical for mitigating climate change. The relative contributions of ${\mathrm{CH}}_{\mathrm{4}}$ emission pathways through ebullition, plant-mediated transport, and diffusion, together with their different transport rates and vulnerability to oxidation, determine the quantity of ${\mathrm{CH}}_{\mathrm{4}}$ 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 ${\mathrm{CH}}_{\mathrm{4}}$ emissions. To improve model simulations of ${\mathrm{CH}}_{\mathrm{4}}$ 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 ${\mathrm{CH}}_{\mathrm{4}}$ flux and concentration data collected in a peatland in northern Minnesota, USA. When model parameters were constrained using observed ${\mathrm{CH}}_{\mathrm{4}}$ fluxes, the ${\mathrm{CH}}_{\mathrm{4}}$ 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 ${\mathrm{CH}}_{\mathrm{4}}$ concentrations, especially in the deep soil layers, compared to the ECT approach. When constraining the EBG model with both ${\mathrm{CH}}_{\mathrm{4}}$ flux and concentration data in model–data fusion, uncertainty of the modeled ${\mathrm{CH}}_{\mathrm{4}}$ concentration profiles was reduced by 78 % to 86 % in comparison to constraints based on ${\mathrm{CH}}_{\mathrm{4}}$ flux data alone. The improved model capability was attributed to the well-constrained parameters regulating the ${\mathrm{CH}}_{\mathrm{4}}$ production and emission pathways. Our results suggest that the EBG modeling approach better characterizes ${\mathrm{CH}}_{\mathrm{4}}$ emission and underlying mechanisms. Moreover, to achieve the best model results both ${\mathrm{CH}}_{\mathrm{4}}$ flux and concentration data are required to constrain model parameterization.

1 Introduction

Methane (${\mathrm{CH}}_{\mathrm{4}}$) emissions from wetlands constitute roughly one third of the global ${\mathrm{CH}}_{\mathrm{4}}$ budget (Denman et al., 2007; Saunois et al., 2020). Methane, following production under anoxic conditions, is stored belowground, oxidized into ${\mathrm{CO}}_{\mathrm{2}}$ by methanotrophs, or emitted into the atmosphere as ${\mathrm{CH}}_{\mathrm{4}}$. The emission of ${\mathrm{CH}}_{\mathrm{4}}$ 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 ${\mathrm{CH}}_{\mathrm{4}}$ oxidation sink has been estimated to be 40 %–70 % of ${\mathrm{CH}}_{\mathrm{4}}$ production, and it can dominate wetland ${\mathrm{CH}}_{\mathrm{4}}$ cycling (Megonigal et al., 2004). The oxidation sink depends substantially on the ${\mathrm{CH}}_{\mathrm{4}}$ 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 ${\mathrm{CH}}_{\mathrm{4}}$ to quickly ascend in a bubble that bypasses the aerobic and anaerobic zones (Epstein and Plesset, 1950). Gases such as ${\mathrm{CH}}_{\mathrm{4}}$ 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 ${\mathrm{CH}}_{\mathrm{4}}$ bypasses aerobic zones of the soil, 20 %–90 % of plant-transported ${\mathrm{CH}}_{\mathrm{4}}$ 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, ${\mathrm{CH}}_{\mathrm{4}}$ 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 ${\mathrm{CH}}_{\mathrm{4}}$ is oxidized before it leaves the soil. Uncertainties in the relative contributions of these pathways to ${\mathrm{CH}}_{\mathrm{4}}$ emission can lead to large errors in the predictions of total ${\mathrm{CH}}_{\mathrm{4}}$ emissions (Bridgham et al., 2013). Despite their importance, the relative contributions of the ${\mathrm{CH}}_{\mathrm{4}}$ 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 ${\mathrm{CH}}_{\mathrm{4}}$ 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 ${\mathrm{CH}}_{\mathrm{4}}$ 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 ${\mathrm{CH}}_{\mathrm{4}}$ emission (Bridgham et al., 2013). Comparisons between modeling approaches and empirical ${\mathrm{CH}}_{\mathrm{4}}$ data suggest that emission pathways may not be well captured by LSMs. For example, plant-mediated ${\mathrm{CH}}_{\mathrm{4}}$ transport by vascular species measured at northern peatlands accounted for 30 %–98 % of the total ${\mathrm{CH}}_{\mathrm{4}}$ 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 ${\mathrm{CH}}_{\mathrm{4}}$ 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 ${\mathrm{CH}}_{\mathrm{4}}$ 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 ${\mathrm{CH}}_{\mathrm{4}}$ emission in LSMs are mainly due to the lack of in situ information, inadequate representation of ${\mathrm{CH}}_{\mathrm{4}}$ processes, and unconstrained parameters used to describe emission pathways (Bridgham et al., 2013; Melton et al., 2013).

Since net ${\mathrm{CH}}_{\mathrm{4}}$ emissions depend on transport mode, all the emission pathways must first be represented correctly in ecosystem models in order to simulate ${\mathrm{CH}}_{\mathrm{4}}$ emission accurately (Blodau, 2002; Tang et al., 2010). Compared to diffusion and plant-mediated ${\mathrm{CH}}_{\mathrm{4}}$ transport, ebullition is less certain and could be the main reason for the mismatch between simulated and observed ${\mathrm{CH}}_{\mathrm{4}}$ 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 ${\mathrm{CH}}_{\mathrm{4}}$ 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, ${\mathrm{CH}}_{\mathrm{4}}$ 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 ${\mathrm{CH}}_{\mathrm{4}}$ 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 ${\mathrm{CH}}_{\mathrm{4}}$ concentration is larger than a defined threshold, the excess ${\mathrm{CH}}_{\mathrm{4}}$ 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 ${\mathrm{CH}}_{\mathrm{4}}$ 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 ${\mathrm{CH}}_{\mathrm{4}}$ 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 ${\mathrm{CH}}_{\mathrm{4}}$ emissions. However, the ability of the EBG approach to represent the relative importance of ${\mathrm{CH}}_{\mathrm{4}}$ 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 ${\mathrm{CH}}_{\mathrm{4}}$ emissions were only useful to constrain a few model parameters, and data assimilation with only ${\mathrm{CH}}_{\mathrm{4}}$ 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 ${\mathrm{CH}}_{\mathrm{4}}$ emission data could only constrain ${\mathrm{CH}}_{\mathrm{4}}$ production-related parameters such as temperature sensitivity ($Q_{\mathrm{10}}$) and basal production rate of ${\mathrm{CH}}_{\mathrm{4}}$ production (Ma et al., 2017). While direct measures of ${\mathrm{CH}}_{\mathrm{4}}$ emission pathways are rare, depth-specific pore water ${\mathrm{CH}}_{\mathrm{4}}$ concentration profiles can help elucidate the relative importance of ${\mathrm{CH}}_{\mathrm{4}}$ emission pathways. Indeed, measured ${\mathrm{CH}}_{\mathrm{4}}$ concentration profiles are critical for constraining the responsive parameters associated with ${\mathrm{CH}}_{\mathrm{4}}$ emission pathways because in process-based ${\mathrm{CH}}_{\mathrm{4}}$ models, all the three emission pathways are calculated based on the ${\mathrm{CH}}_{\mathrm{4}}$ concentration in each soil layer.

To date, few modeling studies have considered ${\mathrm{CH}}_{\mathrm{4}}$ 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 ${\mathrm{CH}}_{\mathrm{4}}$ concentrations, the simulated concentration profiles did not agree well with observations despite good agreements between simulated and observed ${\mathrm{CH}}_{\mathrm{4}}$ emission data (Walter and Heimann, 2000; Tang et al., 2010). Thus, when ${\mathrm{CH}}_{\mathrm{4}}$ emission pathway parameters are calibrated using only net ${\mathrm{CH}}_{\mathrm{4}}$ flux data, models may not realistically represent ${\mathrm{CH}}_{\mathrm{4}}$ 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 ${\mathrm{CH}}_{\mathrm{4}}$ processes. It will be problematic to use these not-yet-well-calibrated parameter sets for climate change predictions or extrapolating ${\mathrm{CH}}_{\mathrm{4}}$ 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 ${\mathrm{CH}}_{\mathrm{4}}$ fluxes and pore water ${\mathrm{CH}}_{\mathrm{4}}$ 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 ${\mathrm{CH}}_{\mathrm{4}}$ profiles better than the ECT approach given its more mechanical representations of bubble formation, gas exchange, and release and (2) pore water ${\mathrm{CH}}_{\mathrm{4}}$ concentration data offer more information for model parameters to reduce the uncertainties in simulated ${\mathrm{CH}}_{\mathrm{4}}$ emission and its pathways.

Table 1

The SPRUCE site data used in this study.

Purpose	Data name	Year	Period	Time step	References
Environmental variables(input) to drive theTECO model	Soil temperature at 0, 5, 10, 20, 30, 40, 50, 100, 200 $\mathrm{cm}$ depth	2011–2016	Whole year	Hourly	Hanson et al.(2015a, b, 2016b)
	Air temperature at 2 $\mathrm{m}$
	Relative Humidity at 2 $\mathrm{m}$
	Wind speed at 10 $\mathrm{m}$
	Precipitation
	Photosynthetically activeradiation (PAR) at 2 $\mathrm{m}$
Water–heat balance andcarbon cycle data tocalibrate the model	Soil moisture at 0, 20 $\mathrm{cm}$	2011–2016	Whole year	Hourly	Same as above
	Water table depth	2011–2016	Whole year	Hourly	Same as above
	Leaf, wood, root biomass	2011–2016	End of growingseason	Once a year	Hanson et al.(2018a, b, 2018)
	Soil $\mathrm{C}$ content	2012	13–15 August	Yearly	Iversen et al. (2014)
	NEE, GPP, ER fluxes	2011–2016	All year aroundexcept mid-winter${}^{\ast }$	1–2 timesa month	Hanson et al. (2014, 2016a)
Data streams used indata–model fusion	Total ${\mathrm{CH}}_{\mathrm{4}}$ emission	2011–2016	All year around except mid-winter${}^{\ast }$	1–2 timesa month	Hanson et al. (2014, 2016a, 2017b)
	Pore water ${\mathrm{CH}}_{\mathrm{4}}$ concentration at 25, 50, 75, 100, 150, 200 $\mathrm{cm}$ depth	2014–2016	Growingseason	Once a month	Wilson et al. (2016)
Generate vertical profile ofheterotrophic respirationand use in calculatingplant-mediated ${\mathrm{CH}}_{\mathrm{4}}$transport	Fine root biomass verticaldistribution	2011–2012	Growingseason	Estimated fromminirhizotron images collectedweekly	Iversen et al.(2018); Malhotraet al. (2020a, b)

${}^{\ast }$ 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 $\mathrm{ha}$ S1 bog in northern Minnesota in the USDA Forest Service Marcell Experimental Forest (47${}^{\circ }$30.476${}^{\prime }$ N, 93${}^{\circ }$27.162${}^{\prime }$ W) to study the responses of northern peatlands to climate warming and elevated atmospheric ${\mathrm{CO}}_{\mathrm{2}}$ concentration (Hanson et al., 2017a). The mean annual temperature from 1961 to 2009 at the SPRUCE site was 3.4 ${}^{\circ }\mathrm{C}$, and the mean annual precipitation was 780 $\mathrm{mm}$ (Sebestyen et al., 2011). The mean peat depth is 2–3 $\mathrm{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 ${}^{\circ }\mathrm{C}$ are paired with two ${\mathrm{CO}}_{\mathrm{2}}$ treatments (ambient or $\sim$ 400 $\mathrm{ppm}$ and 900 $\mathrm{ppm}$) in open-top infrastructures (12 $\mathrm{m}$ $\times$ 8 $\mathrm{m}$). Deep peat warming began in June 2014, aboveground warming began in August 2015, and elevated ${\mathrm{CO}}_{\mathrm{2}}$ 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 ${\mathrm{CH}}_{\mathrm{4}}$ processes. Modeling ${\mathrm{CH}}_{\mathrm{4}}$ emissions in response to warming and elevated ${\mathrm{CO}}_{\mathrm{2}}$ 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 ${\mathrm{CH}}_{\mathrm{4}}$ 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).

${\mathrm{CH}}_{\mathrm{4}}$ fluxes and pore water ${\mathrm{CH}}_{\mathrm{4}}$ concentrations were used for data assimilation. We averaged the data from all ambient plots measured on the same dates to represent the site-level ${\mathrm{CH}}_{\mathrm{4}}$ emissions and concentration profiles. Variations among different ambient plots were not considered in this study. In total, 45 daily ${\mathrm{CH}}_{\mathrm{4}}$ emission measurements were obtained from ambient plots from 2011 to 2016. In situ pore water ${\mathrm{CH}}_{\mathrm{4}}$ 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 $\mathrm{cm}$ depths, respectively (Wilson et al., 2016). Piezometers consisted of a $<$ 1 $\mathrm{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 $\mathrm{cm}$ depth. Samples were immediately filtered in the field through 0.7 $\mathrm{\mu }\mathrm{M}$ Whatman glass-fiber filters and stored in pre-evacuated, septum-sealed glass vials. Phosphoric acid (1 $\mathrm{mL}$, 20 %) was added to preserve each sample during shipment to Florida State University for analyzing ${\mathrm{CH}}_{\mathrm{4}}$ 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 ($\mathrm{N}$) transfer, and soil ${\mathrm{CH}}_{\mathrm{4}}$ 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 ${\mathrm{CH}}_{\mathrm{4}}$ 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 $\mathrm{C}$–$\mathrm{N}$ transfer module simulates the movement of $\mathrm{C}$ and $\mathrm{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 $\mathrm{C}$ pool (Luo and Reynolds, 1999).

Figure 1

Conceptual structure of the ${\mathrm{CH}}_{\mathrm{4}}$ emission module in TECO_SPRUCE.

[Figure omitted. See PDF]

The ${\mathrm{CH}}_{\mathrm{4}}$ module simulates the transient, vertical dynamics of ${\mathrm{CH}}_{\mathrm{4}}$ production, oxidation, and belowground transport (via ebullition, plant-mediated transport, and diffusion), as well as ${\mathrm{CH}}_{\mathrm{4}}$ 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 $\mathrm{cm}$ thick, whereas each of the rest of the layers are 20 $\mathrm{cm}$ thick. Within each soil layer, ${\mathrm{CH}}_{\mathrm{4}}$ concentration dynamics are calculated by a transient reaction equation with ${\mathrm{CH}}_{\mathrm{4}}$ production, ${\mathrm{CH}}_{\mathrm{4}}$ oxidation, released bubbles, plant-mediated transport, and the diffusion of ${\mathrm{CH}}_{\mathrm{4}}$ 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 ${\mathrm{CH}}_{\mathrm{4}}$ production ($P_{\text{ro}}$) within the catotelm is directly related to heterotrophic respiration from soil and litter ($R_{\mathrm{h}}$, $\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{2}}{\mathrm{h}}^{-\mathrm{1}}$) via the following equation: 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}},$$ where $f_{{\mathrm{CH}}_{\mathrm{4}}}$ is an ecosystem-specific conversion scaler describing the fraction of anaerobically mineralized $\mathrm{C}$ atoms becoming ${\mathrm{CH}}_{\mathrm{4}}$. The parameters $f_{\text{stp}}$, $f_{\text{pH}}$, and $f_{\text{red}}$ are environmental scalers, representing the effects of soil temperature, pH, and redox potential, respectively, on ${\mathrm{CH}}_{\mathrm{4}}$ production. Total emission of ${\mathrm{CH}}_{\mathrm{4}}$ from the soil to the atmosphere is calculated as the sum of ${\mathrm{CH}}_{\mathrm{4}}$ ebullition from saturated soil layers, plant-mediated ${\mathrm{CH}}_{\mathrm{4}}$ 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 ${\mathrm{CH}}_{\mathrm{4}}$ 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 ${\mathrm{CH}}_{\mathrm{4}}$ 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, ${\mathrm{CH}}_{\mathrm{4}}$ 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 ${\mathrm{CH}}_{\mathrm{4}}$ 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 ${\mathrm{CH}}_{\mathrm{4}}$ at a given temperature was calculated using a statistical model used by Yamamoto et al. (1976) based on the empirical data: 2$$V=\mathrm{0.05708}-\mathrm{0.001545}T+\mathrm{0.00002069}{T}^{\mathrm{2}},$$ 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 ${\mathrm{CH}}_{\mathrm{4}}$ dissolved per volume of water was converted into grams using the ideal gas law: 3$${\left[{\mathrm{CH}}_{\mathrm{4}}\right]}_{\text{thre}}=PV\cdot C/RT,$$ where ${\left[{\mathrm{CH}}_{\mathrm{4}}\right]}_{\text{thre}}$ is the maximum concentration threshold ($\mathrm{g}\mathrm{C}{\mathrm{m}}^{-\mathrm{3}}$), $P$ is the sum of the atmospheric and hydrostatic pressures ($\mathrm{Pa}$), $V$ is the Bunsen solubility coefficient as in Eq. (), the constant $C$ is the atomic weight of carbon (12 $\mathrm{g}{\mathrm{mol}}^{-\mathrm{1}}$), 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 ($\mathrm{K}$). Then the ${\mathrm{CH}}_{\mathrm{4}}$ ebullition flux can be calculated using the following equation: 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.,$$ where $E_{\text{bu}}(z,t)$ is the ebullition flux of ${\mathrm{CH}}_{\mathrm{4}}$ ($\mathrm{g}\mathrm{C}{\mathrm{h}}^{-\mathrm{1}}$) to the lowest air layer, $K_{\text{ebu}}$ is a rate constant of 1.0 ${\mathrm{h}}^{-\mathrm{1}}$ (Walter and Heimann, 2000; Zhuang et al., 2004, 2006), and $\left[{\mathrm{CH}}_{\mathrm{4}}\right](z,t)$ is the pore water ${\mathrm{CH}}_{\mathrm{4}}$ 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 ${\mathrm{CH}}_{\mathrm{4}}$ 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 ${\mathrm{CH}}_{\mathrm{4}}$ concentration exceeds the limit that the water can withhold based on Henry's law, then excess ${\mathrm{CH}}_{\mathrm{4}}$ is converted to a gaseous volume calculated using the predefined bubble ${\mathrm{CH}}_{\mathrm{4}}$ mixing ratio ($f$). This gaseous volume is divided evenly into a certain number of bubbles ($N_{\text{bub}}$). $N_{\text{bub}}$ is a unitless tuning parameter ranging between 5 and 500 in each 10 $\mathrm{cm}$ thick soil layer and 10 and 1000 in each 20 $\mathrm{cm}$ thick soil layer, which essentially controls the mass exchange rate between the gas volume and the pore water. The ${\mathrm{CH}}_{\mathrm{4}}$ exchange between the stationary bubbles and the pore water ($Q_{\text{ebu}}$) is calculated using the equation proposed by Epstein and Plesset (1950): 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),$$ where $r$ is the radius of a bubble ($\mathrm{m}$), $D_{\mathrm{w}}$ is the ${\mathrm{CH}}_{\mathrm{4}}$ diffusion coefficient in water (${\mathrm{m}}^{\mathrm{2}}{\mathrm{s}}^{-\mathrm{1}}$) calculated using the quadratic curve of observed diffusivities against temperatures (Broecker and Peng, 1974), $V_{\mathrm{w}}$ is the amount of water in this layer (${\mathrm{m}}^{\mathrm{3}}$), $c_{\mathrm{w}}$ is dissolved ${\mathrm{CH}}_{\mathrm{4}}$ concentration in the pore water, and ${\mathrm{H}}^{\mathrm{cc}}$ is the dimensionless Henry solubility of ${\mathrm{CH}}_{\mathrm{4}}$ calculated following Sander (1999). $P$, $R$, and $T$ are the same as in Eq. (). A negative value of $Q_{\text{ebu}}$ indicates ${\mathrm{CH}}_{\mathrm{4}}$ 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_{\text{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_{\text{maxfraction}}$ is the free-phase gas-filled fraction of the pore space in the soil layer above which ebullition occurs, $c_{\mathrm{b}}$ is the ${\mathrm{CH}}_{\mathrm{4}}$ concentration in a bubble ($\mathrm{mol}{\mathrm{m}}^{-\mathrm{3}}$), $V_{\mathrm{B}}$ is the total volume of all bubbles, and $\mathrm{\Delta }{V}_{\mathrm{B}}$ is the change in the total volume due to the diffusive gas exchange in Eq. (). The amount of ${\mathrm{CH}}_{\mathrm{4}}$ in all bubbles after each time step is 8$$n_{\mathrm{b}}={\displaystyle \frac{fP{V}_{\mathrm{B}}^{\prime }}{RT}},$$ where $f$ is the predefined bubble ${\mathrm{CH}}_{\mathrm{4}}$ 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.

Process	Symbol	Range	Units	Definition	References
${\mathrm{CH}}_{\mathrm{4}}$ production	$f_{{\mathrm{CH}}_{\mathrm{4}}}$	[0.0, 0.7]	–	Fraction of anaerobically mineralized $\mathrm{C}$ atoms becoming ${\mathrm{CH}}_{\mathrm{4}}$	Zhuang et al. (2004); Segers (1998); Zhu et al. (2014)
	$Q_{\mathrm{10}\text{pro}}$	[0.0, 10]	–	$Q_{\mathrm{10}}$ for ${\mathrm{CH}}_{\mathrm{4}}$ production	Walter and Heimann (2000)
${\mathrm{CH}}_{\mathrm{4}}$ oxidation	${\mathrm{O}}_{max}$	[3.0, 45.0]	$\mathrm{\mu }\mathrm{mol}{\mathrm{L}}^{-\mathrm{1}}{\mathrm{h}}^{-\mathrm{1}}$	Maximum oxidation rate	Zhuang et al. (2004)
${\mathrm{CH}}_{\mathrm{4}}$ ebullition	$f$	[0.01, 0.5]	$\mathrm{mol}{\mathrm{mol}}^{-\mathrm{1}}$	${\mathrm{CH}}_{\mathrm{4}}$ mixing ratio in bubbles	Tang et al. (2010); Peltola et al. (2018)
	bubprob	[0.01, 0.5]	–	Probability that a bubble will get trapped at one layer	Tang et al. (2010); Peltola et al. (2018)
	$V_{\text{maxfraction}}$	[0.01, 0.2]	–	Maximum fraction of volume occupied by bubbles	Peltola et al. (2018)
Plant-mediated transportation	$T_{\text{veg}}$	[0.01, 15.0]	–	Capability of conducting ${\mathrm{CH}}_{\mathrm{4}}$ gas at plant community level	Walter and Heimann (2000); Zhuang et al. (2004)

The values of bubprob, $f$, and $V_{\text{maxfraction}}$ are dependent on the soil texture, porosity, water content, etc. and have been found to significantly affect the modeled ${\mathrm{CH}}_{\mathrm{4}}$ 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_{\text{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_{\text{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

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\},$$ where $Z_i\left(t\right)$ is the $i$th observation stream at time $t$, $X\left(t\right)$ is the model simulated value, and ${\mathit{\sigma }}_i\left(t\right)$ 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.

Data assimilation runs	Ebullition approachesembedded with TECO	Observation data streams used for constraining the parameters
ECT_F	ECT	${\mathrm{CH}}_{\mathrm{4}}$ fluxes
ECT_FC	ECT	${\mathrm{CH}}_{\mathrm{4}}$ fluxes $+$ pore water ${\mathrm{CH}}_{\mathrm{4}}$ concentration profiles
EBG_F	EBG	${\mathrm{CH}}_{\mathrm{4}}$ fluxes
EBG_FC	EBG	${\mathrm{CH}}_{\mathrm{4}}$ fluxes $+$ pore water ${\mathrm{CH}}_{\mathrm{4}}$ concentration profiles

Parameters directly regulating ${\mathrm{CH}}_{\mathrm{4}}$ 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_{\mathrm{10}\text{pro}}$, $O_{max}$, and $T_{\text{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 ${\mathrm{CH}}_{\mathrm{4}}$ emission and pore water ${\mathrm{CH}}_{\mathrm{4}}$ 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 ${\mathrm{CH}}_{\mathrm{4}}$ 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 ${\mathrm{CH}}_{\mathrm{4}}$ flux data. Then we compare results from EBG_F and EBG_FC to show the ability of pore water ${\mathrm{CH}}_{\mathrm{4}}$ concentration data to help constrain the parameters related to the ${\mathrm{CH}}_{\mathrm{4}}$ emission pathways. Model performance was evaluated against the observed data using root mean square error (RMSE). Model uncertainties in pore water ${\mathrm{CH}}_{\mathrm{4}}$ 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.

Parameter	Ebullition method	Observation data	Posterior distribution mean $\pm$ SD	MLE${}^{\ast }$	Parameter class
$f_{{\mathrm{CH}}_{\mathrm{4}}}$	ECT	Flux	0.16 $\pm$ 0.016	0.16	Bell-shaped
	EBG	Flux	0.17 $\pm$ 0.023	0.17	Bell-shaped
	EBG	Flux $+$ concentration	0.15 $\pm$ 0.021	0.15	Bell-shaped
$Q_{\mathrm{10}\text{pro}}$	ECT	Flux	3.0 $\pm$ 0.85	3.0	bell-shaped
	EBG	Flux	2.69 $\pm$ 0.82	2.69	Bell-shaped
	EBG	Flux $+$ concentration	3.21 $\pm$ 1.07	3.21	Bell-shaped
$O_{max}$	ECT	Flux	22.8 $\pm$ 12.1	–	Flat
	EBG	Flux	22.5 $\pm$ 12.1	–	Flat
	EBG	Flux $+$ concentration	22.4 $\pm$ 12.0	–	Flat
$T_{\text{veg}}$	ECT	Flux	7.7 $\pm$ 4.0	–	Flat
	EBG	Flux	5.8 $\pm$ 4.0	–	Flat
	EBG	Flux $+$ concentration	1.43 $\pm$ 0.46	1.43	Bell-shaped
f	EBG	Flux	0.11 $\pm$ 4.0	0.11	Edge-hitting
	EBG	Flux $+$ concentration	0.29 $\pm$ 0.46	0.29	Bell-shaped
bubprob	EBG	Flux	0.22 $\pm$ 0.87	0.19	Edge-hitting; big variance
	EBG	Flux $+$ concentration	0.25 $\pm$ 0.015	0.23	Bell-shaped
$V_{\text{maxfraction}}$	EBG	Flux	0.1 $\pm$ 0.13	0.08	Bell-shaped
	EBG	Flux $+$ concentration	0.11 $\pm$ 0.12	0.1	Bell-shaped

${}^{\ast }$ 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 ${\mathrm{CH}}_{\mathrm{4}}$ emission data (ECT_F). The green lines are the PPDs from the EBG model structure trained with ${\mathrm{CH}}_{\mathrm{4}}$ emission data (EBG_F). The purple lines are the PPDs from EBG trained with both ${\mathrm{CH}}_{\mathrm{4}}$ emission and ${\mathrm{CH}}_{\mathrm{4}}$ concentration data (EBG_FC). Well-constrained parameters have a unimodal distribution, whereas poorly constrained parameter distributions tend to be flat.

[Figure omitted. See PDF]

3 Results

3.1

Parameter optimization using ${\mathrm{CH}}_{\mathrm{4}}$ flux data with different model structures

The ${\mathrm{CH}}_{\mathrm{4}}$ production-related parameters, $f_{{\mathrm{CH}}_{\mathrm{4}}}$ (fraction of anaerobically mineralized $\mathrm{C}$ atoms becoming ${\mathrm{CH}}_{\mathrm{4}}$) and $Q_{\mathrm{10}\text{pro}}$ (temperature sensitivity of ${\mathrm{CH}}_{\mathrm{4}}$ production), were constrained using the ${\mathrm{CH}}_{\mathrm{4}}$ 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_{\mathrm{10}\text{pro}}$ decreased from 3.0 to 2.69 in the EBG approach compared to the ECT approach. In contrast, $O_{max}$ and $T_{\text{veg}}$ were not constrained in either model with large uncertainties in model-estimated ${\mathrm{CH}}_{\mathrm{4}}$ 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 ${\mathrm{CH}}_{\mathrm{4}}$ emission flux data, $V_{\text{maxfraction}}$ (maximum fraction of volume occupied by bubbles) was constrained with a unimodal shaped posterior distribution (Fig. 2g), $f({\mathrm{CH}}_{\mathrm{4}}$ 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 ${\mathrm{CH}}_{\mathrm{4}}$ emissions from static chamber measurements. Error bars are the standard errors. Blue lines are ECT-model-simulated ${\mathrm{CH}}_{\mathrm{4}}$ emissions based on the parameter probability distributions (PPDs) constrained by ${\mathrm{CH}}_{\mathrm{4}}$ flux data (ECT_F, RMSE $=$ 0.61). Green lines are EBG-model-simulated ${\mathrm{CH}}_{\mathrm{4}}$ emissions based on the PPDs constrained by ${\mathrm{CH}}_{\mathrm{4}}$ flux data (EBG_F, RMSE $=$ 0.53). Purple lines are EBG-model-simulated emissions based on the PPDs constrained by both ${\mathrm{CH}}_{\mathrm{4}}$ 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.

[Figure omitted. See PDF]

3.2 Evaluations of model structures against the observed data

Using ${\mathrm{CH}}_{\mathrm{4}}$ emission data to constrain the parameters, EBG-simulated ${\mathrm{CH}}_{\mathrm{4}}$ 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 ${\mathrm{CH}}_{\mathrm{4}}$ emissions (Fig. 3c) than ECT (Fig. 3a). The simulated contributions from plant-mediated transport, diffusion, and ebullition were 40.7 $\pm$ 8.0 %, 35.7 $\pm$ 8.7 %, and 23.5 $\pm$ 9.4 %, respectively, in ECT_F (Fig. 3b) and 38.4 $\pm$ 13.9 %, 38.7 $\pm$ 9.9 %, and 22.7 $\pm$ 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 ${\mathrm{CH}}_{\mathrm{4}}$ 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 ${\mathrm{CH}}_{\mathrm{4}}$ flux data (ECT_F). Green lines are the EBG-model-simulated concentrations based on the PPDs constrained by ${\mathrm{CH}}_{\mathrm{4}}$ flux data (EBG_F). Purple lines are the EBG-model-simulated concentrations based on the PPDs constrained by both ${\mathrm{CH}}_{\mathrm{4}}$ 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.

[Figure omitted. See PDF]

The ECT model constrained by ${\mathrm{CH}}_{\mathrm{4}}$ flux data could not reproduce well the patterns of the observed pore water ${\mathrm{CH}}_{\mathrm{4}}$ concentrations, especially in the deep peat layers (RMSE $=$ 0.77, Fig. 4, ECT_F). When calibrated by ${\mathrm{CH}}_{\mathrm{4}}$ flux data alone, the EBG approach for ebullition captured deep layer ${\mathrm{CH}}_{\mathrm{4}}$ 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 ${\mathrm{CH}}_{\mathrm{4}}$ concentration profiles, especially in the deep peat layers, mainly due to the unconstrained $T_{\text{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 ${\mathrm{CH}}_{\mathrm{4}}$ processes

For the ECT approach, as described earlier, assimilating the observed ${\mathrm{CH}}_{\mathrm{4}}$ flux could constrain parameters such as $f_{{\mathrm{CH}}_{\mathrm{4}}}$ and $Q_{\mathrm{10}\text{pro}}$. However, when using both observed ${\mathrm{CH}}_{\mathrm{4}}$ 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 ${\mathrm{CH}}_{\mathrm{4}}$ emissions and pore water ${\mathrm{CH}}_{\mathrm{4}}$ concentrations.

In contrast to the ECT approach, incorporation of pore water ${\mathrm{CH}}_{\mathrm{4}}$ concentration data in the EBG approach greatly improved parameter estimations. While $T_{\text{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 ${\mathrm{CH}}_{\mathrm{4}}$ flux and pore water ${\mathrm{CH}}_{\mathrm{4}}$ concentration data streams were used in the data–model fusion (Table 4, Fig. 2, EBG_FC). Compared to EBG_F, the parameter $T_{\text{veg}}$ was well constrained to a very small range of 1.43 $\pm$ 0.46, and the parameter bubprob was also well constrained to a range of 0.25 $\pm$ 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_{\mathrm{10}\text{pro}}$ 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 $\pm$ 4.0 to 0.29 $\pm$ 0.46 when the pore water ${\mathrm{CH}}_{\mathrm{4}}$ concentration information was added into data assimilation. All the emission pathway-related parameters ($T_{\text{veg}}$, bubprob, f, and $V_{\text{maxfraction}}$) were well constrained once the pore water ${\mathrm{CH}}_{\mathrm{4}}$ concentration profile information was added to data–model fusion. However, incorporation of the pore water ${\mathrm{CH}}_{\mathrm{4}}$ 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 ${\mathrm{CH}}_{\mathrm{4}}$ 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 $\pm$ 4.9 %, 58.1 $\pm$ 5.1 %, and 9.9 $\pm$ 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). ${\mathrm{CH}}_{\mathrm{4}}$ flux and concentration data together reduced the uncertainty of the modeled ${\mathrm{CH}}_{\mathrm{4}}$ 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 ${\mathrm{CH}}_{\mathrm{4}}$ concentration profiles was decreased mainly due to the well-constrained parameters regulating the ${\mathrm{CH}}_{\mathrm{4}}$ 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., ${\mathrm{CH}}_{\mathrm{4}}$ emission and pore water ${\mathrm{CH}}_{\mathrm{4}}$ concentration data, in simulating peatland ${\mathrm{CH}}_{\mathrm{4}}$ emission and its pathways.

4.1

Better representing ${\mathrm{CH}}_{\mathrm{4}}$ emission and pore water ${\mathrm{CH}}_{\mathrm{4}}$ 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 ${\mathrm{CH}}_{\mathrm{4}}$ emissions ($R^{\mathrm{2}}=$ 0.63) when compared with the ECT ($R^{\mathrm{2}}=$ 0.56) and EPT ($R^{\mathrm{2}}=$ 0.35) methods (Peltola et al., 2018). We also found that the EBG-simulated ${\mathrm{CH}}_{\mathrm{4}}$ 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 ${\mathrm{CH}}_{\mathrm{4}}$ emissions as EBG, the mean annual ${\mathrm{CH}}_{\mathrm{4}}$ 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 ${\mathrm{CH}}_{\mathrm{4}}$ stored belowground and distinctly different distributions of ${\mathrm{CH}}_{\mathrm{4}}$ along the soil profiles. In line with their results, we found that the ECT method produced much higher pore water ${\mathrm{CH}}_{\mathrm{4}}$ concentrations than the EBG method, especially in the deep layers (Fig. 4).

Of the few modeling studies that compared results with observed belowground ${\mathrm{CH}}_{\mathrm{4}}$ concentration, Walter et al. (2001) simulated ${\mathrm{CH}}_{\mathrm{4}}$ 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 ${\mathrm{CH}}_{\mathrm{4}}$ concentrations within the first 50 $\mathrm{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 ${\mathrm{CH}}_{\mathrm{4}}$ concentration profile, although a mismatch in the concentration profile remained, especially from the model that best reproduced observed ${\mathrm{CH}}_{\mathrm{4}}$ emissions. The recently developed land surface model (LSM) with a new microbial-functional group-based ${\mathrm{CH}}_{\mathrm{4}}$ 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 ${\mathrm{CH}}_{\mathrm{4}}$ production and net fluxes, but ${\mathrm{CH}}_{\mathrm{4}}$ concentrations in soil layers deeper than 1 $\mathrm{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 $\pm$ 4.9 % for PMT, 58.1 $\pm$ 5.1 % for diffusion, and 9.9 $\pm$ 6.1 % for ebullition). In our study, when training the modified ECT model with both ${\mathrm{CH}}_{\mathrm{4}}$ 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 ${\mathrm{CH}}_{\mathrm{4}}$ fluxes and the patterns of pore water ${\mathrm{CH}}_{\mathrm{4}}$ concentrations no matter the combinations of parameters used. In contrast, the EBG method could capture observed ${\mathrm{CH}}_{\mathrm{4}}$ emissions and the patterns of pore water ${\mathrm{CH}}_{\mathrm{4}}$ 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 ${\mathrm{CH}}_{\mathrm{4}}$ emissions, the flux-constrained parameter distributions varied between the two methods. For example, $f_{{\mathrm{CH}}_{\mathrm{4}}}$ increased, but $Q_{\mathrm{10}\text{pro}}$ 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 ${\mathrm{CH}}_{\mathrm{4}}$ production and emission processes and the underlying mechanisms.

4.2

Pool-based ${\mathrm{CH}}_{\mathrm{4}}$ 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 ${\mathrm{CH}}_{\mathrm{4}}$ emissions did not necessarily reproduce a realistic pore water concentration profile (Figs. 3 and 4). By comparing the parameter posterior distributions trained by observed ${\mathrm{CH}}_{\mathrm{4}}$ emissions with and without observed pore water concentration profiles using the same model structure, TECO_SPRUCE_EBG, we revealed that ${\mathrm{CH}}_{\mathrm{4}}$ emission data could constrain the ${\mathrm{CH}}_{\mathrm{4}}$ production-related parameters $f_{{\mathrm{CH}}_{\mathrm{4}}}$ and $Q_{\mathrm{10}\text{pro}}$ and ebullition-related parameter $V_{\text{maxfraction}}$ very well. The ebullition-related parameter $f$ was edge hitting, but parameter bubprob and plant transport-related parameter $T_{\text{veg}}$ remained unconstrained, causing large uncertainty in simulated ebullition and plant-mediated transport (Table 4). However, by training the model with both ${\mathrm{CH}}_{\mathrm{4}}$ emission and pore water ${\mathrm{CH}}_{\mathrm{4}}$ concentration data, the parameters regulating ${\mathrm{CH}}_{\mathrm{4}}$ production, plant transport, and ebullition were all well constrained (Fig. 2). This is because the vertical concentration profile of ${\mathrm{CH}}_{\mathrm{4}}$ is a balance between the dynamic ${\mathrm{CH}}_{\mathrm{4}}$ production, oxidation, and three emission pathways. The constrained parameters contributed to a more accurate estimation of pore water ${\mathrm{CH}}_{\mathrm{4}}$ 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 ${\mathrm{CH}}_{\mathrm{4}}$ emission data mainly constrained parameters that represented instant responses to temperature change ($Q_{\mathrm{10}\text{pro}}$) and input rate from the source pool ($f_{{\mathrm{CH}}_{\mathrm{4}}}$). The pore water ${\mathrm{CH}}_{\mathrm{4}}$ concentration data contributed to constraining the allocation rates of ${\mathrm{CH}}_{\mathrm{4}}$ to the different emission pathways. Due to the different information contained between ${\mathrm{CH}}_{\mathrm{4}}$ 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 ${\mathrm{CH}}_{\mathrm{4}}$ models.

It needs to be noted that there is a large disagreement in simulated relative contribution by ebullition between ${\mathrm{CH}}_{\mathrm{4}}$ 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 ${\mathrm{CH}}_{\mathrm{4}}$ 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 ${\mathrm{CH}}_{\mathrm{4}}$ 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 ${\mathrm{m}}^{\mathrm{2}}$ will be deployed in each plot to measure ${\mathrm{CO}}_{\mathrm{2}}$ and ${\mathrm{CH}}_{\mathrm{4}}$ fluxes. Along with the continued ${\mathrm{CH}}_{\mathrm{4}}$ 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 ${\mathrm{CH}}_{\mathrm{4}}$ 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 ${\mathrm{CH}}_{\mathrm{4}}$ 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 ${\mathrm{CH}}_{\mathrm{4}}$ concentrations (Tang et al., 2010; Riley et al., 2011). In LSMs, plant transport exclusively dominates ${\mathrm{CH}}_{\mathrm{4}}$ 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 ${\mathrm{CH}}_{\mathrm{4}}$ flux and pore water ${\mathrm{CH}}_{\mathrm{4}}$ 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 ${\mathrm{CH}}_{\mathrm{4}}$ modeling studies and suggest future directions for both modeling and experimental efforts, namely (1) the under-described ${\mathrm{CH}}_{\mathrm{4}}$ processes in models and (2) the lack of observational data to constrain key ${\mathrm{CH}}_{\mathrm{4}}$ processes in the models. More explicit ${\mathrm{CH}}_{\mathrm{4}}$ processes are needed in modeling ${\mathrm{CH}}_{\mathrm{4}}$ 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 ${\mathrm{CH}}_{\mathrm{4}}$ dynamics in wetland soils are based on the premise that the oxidation of ${\mathrm{CH}}_{\mathrm{4}}$ occurs only in aerobic environments (Wania et al., 2010; Riley et al., 2011; Bridgham et al., 2013). However, the anaerobic oxidation of ${\mathrm{CH}}_{\mathrm{4}}$ may be an important sink for ${\mathrm{CH}}_{\mathrm{4}}$, 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 ${\mathrm{CH}}_{\mathrm{4}}$ model was developed accounting for both aerobic and anaerobic ${\mathrm{CH}}_{\mathrm{4}}$ oxidations, and this model has been validated against the concentration of ${\mathrm{CH}}_{\mathrm{4}}$ and ${\mathrm{CO}}_{\mathrm{2}}$ from incubation data (Xu et al., 2015). In Xu et al. model, 7 out of total 33 key ${\mathrm{CH}}_{\mathrm{4}}$ process parameters control ${\mathrm{CH}}_{\mathrm{4}}$ oxidation, and their values vary widely across different ecosystems and environmental conditions. The incorporation of anaerobic ${\mathrm{CH}}_{\mathrm{4}}$ oxidation into LSMs may help improve the calculations of ${\mathrm{CH}}_{\mathrm{4}}$ oxidation if the uncertainties from these ${\mathrm{CH}}_{\mathrm{4}}$-oxidation-related parameters can be reduced.

While more comprehensive and process-based models for simulating all the processes or mechanisms involved in ${\mathrm{CH}}_{\mathrm{4}}$ 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 ${\mathrm{CO}}_{\mathrm{2}}$ becoming ${\mathrm{CH}}_{\mathrm{4}}$ ($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_{\text{veg}}$ was used as a proxy of the ability of the whole plant community (e.g., biomass and abundance) to emit ${\mathrm{CH}}_{\mathrm{4}}$. 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_{\text{veg}}$ was well constrained by the observed ${\mathrm{CH}}_{\mathrm{4}}$ emission and concentration data at a range of 1.43 $\pm$ 0.46, which indicated that the ability of the plant community to emit ${\mathrm{CH}}_{\mathrm{4}}$ at this site was low (compared to its prior knowledge of 0.01–15, Table 2). Empirical measurements of plant-mediated ${\mathrm{CH}}_{\mathrm{4}}$ 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 ${\mathrm{CH}}_{\mathrm{4}}$ emission pathways is critical to mechanistically model future ${\mathrm{CH}}_{\mathrm{4}}$ dynamics. Acknowledging that pore water ${\mathrm{CH}}_{\mathrm{4}}$ concentration is the driving force for each emission pathway, we evaluated the ability of two ebullition modeling approaches to reproduce observed ${\mathrm{CH}}_{\mathrm{4}}$ 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 ${\mathrm{CH}}_{\mathrm{4}}$ emissions and ${\mathrm{CH}}_{\mathrm{4}}$ concentration profiles much better than the modified ebullition concentration threshold approach (ECT), especially for ${\mathrm{CH}}_{\mathrm{4}}$ concentrations in the deeper soil layers. By assimilating the net ${\mathrm{CH}}_{\mathrm{4}}$ emission and belowground ${\mathrm{CH}}_{\mathrm{4}}$ concentration data into the models, we substantially reduced the uncertainties of modeled ${\mathrm{CH}}_{\mathrm{4}}$ emissions from the emission pathways involved. While net ${\mathrm{CH}}_{\mathrm{4}}$ efflux data are often the only data stream for ${\mathrm{CH}}_{\mathrm{4}}$ model validations, we recommend that more attention be given to in situ measurements of the pore water ${\mathrm{CH}}_{\mathrm{4}}$ 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 ${\mathrm{CH}}_{\mathrm{4}}$ 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 ${\mathrm{CH}}_{\mathrm{4}}$ emission and its transport processes are more realistic in response to climate change scenarios. Future studies should also include anaerobic ${\mathrm{CH}}_{\mathrm{4}}$ oxidation into LSMs and constrain its parameters to better predict wetland ${\mathrm{CH}}_{\mathrm{4}}$ emissions.

Code and data availability

All data sets from this study are publicly available at project websites. Relevant measurements were obtained from the SPRUCE web page (http://mnspruce.ornl.gov/, last access: 23 April 2022), as well as the archival ftp site (ftp://sprucedata.ornl.gov, last access: 23 April 2022). The data can be acquired upon reasonable request. Results used to generate figures are available at 10.5281/zenodo.5722449 (Ma et al., 2021).

The supplement related to this article is available online at: https://doi.org/10.5194/bg-19-2245-2022-supplement.

Author contributions

SM and YL designed the project. SM, JJ, YH, and XL carried out modeling study. RMW, JPC, CMI, AM, PJH, and SB provided experimental data for model evaluations and parameter optimization. SM, LJ, and YL prepared the manuscript. All authors contributed to analyzing and interpreting the modeling results, as well as improving the manuscript.

Competing interests

The contact author has declared that neither they nor their co-authors have any competing interests.

Disclaimer

Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Acknowledgements

Part of this research was carried out at the Jet Propulsion Laboratory, California Institute of Technology, under a contract with the National Aeronautics and Space Administration. We thank Zhenggang Du's help with discussing the data–model fusion techniques.

Financial support

This work was primarily funded by subcontract 4000158404 from Oak Ridge National Laboratory to Northern Arizona University. Oak Ridge National Laboratory is managed by UT-Battelle, LLC, for the U.S. Department of Energy under contract DE-AC05-00OR22725. The SPRUCE (Spruce and Peatland Responses Under Changing Environments) project is supported by the Biological and Environmental Research program in the U.S. Department of Energy's Office of Science.

Review statement

This paper was edited by Akihiko Ito and reviewed by two anonymous referees.