Abstract. Bubbling is an important pathway of methane emissions from wetland ecosystems. However the concentration-based threshold function approach in current biogeochemistry models of methane is not sufficient to represent the complex ebullition process. Here we revise an extant process-based biogeochemistry model, the Terrestrial Ecosystem Model into a multi-substance model (CH₄, O₂, CO₂ and N₂) to simulate methane production, oxidation, and transport (particularly ebullition) with different model complexities. When ebullition is modeled with a concentration-based threshold function and if the inhibition effect of oxygen on methane production and the competition for oxygen between methanotrophy and heterotrophic respiration are retained, the model becomes a two-substance system. Ignoring the role of oxygen, while still modeling ebullition with a concentration-based threshold function, reduces the model to a one-substance system. These models were tested through a group of sensitivity analyses using data from two temperate peatland sites in Michigan. We demonstrate that only the four-substance model with a pressure-based ebullition algorithm is able to capture the episodic emissions induced by a sudden decrease in atmospheric pressure or by a sudden drop in water table. All models captured the retardation effect on methane efflux from an increase in surface standing water which results from the inhibition of diffusion and the increase in rhizospheric oxidation. We conclude that to more accurately account for the effects of atmospheric pressure dynamics and standing water on methane effluxes, the multi-substance model with a pressure-based ebullition algorithm should be used in the future to quantify global wetland CH₄ emissions. Further, to more accurately simulate the pore water gas concentrations and different pathways of methane transport, an exponential root distribution function should be used and the phase-related parameters should be treated as temperature dependent.