Quantifying Wetland Methane Emissions with Process-based Models of Different Complexities

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 4 , O 2 , CO 2 and N 2) 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 4 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.


Introduction
Methane (CH 4 ) emitted from natural wetlandsis a significant component of its atmospheric budget.Biogeochemistry and atmospheric inversion models estimate the total wetland emissions to be 100-230 Tg CH 4 y −1 , around 25% of the global emissions into the atmosphere under the current climate condition (Denman et al., 2007).Inverse modeling estimates the strengths of various CH 4 sources and sinks by comparing the model simulated CH 4 concentrations to spatially discrete and temporally continuous observations of the atmospheric CH 4 concentrations (e.g.Houweling et al., 1999).Since all sources/sinks are treated simultaneously in the inversion, the total CH 4 emissions into the atmosphere can be well constrained.However, there are various limitations including the sparse in-situ observation networks of atmospheric CH 4 and unclear sources and sinks due to insufficient understanding of the biogeochemical processes.As a result, the estimates for different sources/sinks from inverse Published by Copernicus Publications on behalf of the European Geosciences Union.
modeling are usually subject to great uncertainties.Processbased models integrate and extrapolate the knowledge from field studies at limited sites to regional and global scales.Because of sparse site-level information and inadequate representation of CH 4 processes in these models, the uncertainties in the quantification from biogeochemical modeling are also substantial (e.g., Walter et al., 2001;Zhuang et al., 2004Zhuang et al., , 2009;;Denman et al., 2007).
To date, a group of process-based models with different complexities have been developed to quantify the spatial and temporal patterns of wetland CH 4 emissions.Among them, the one-substance models are widely used (e.g., Walter and Heimann, 2000;Zhuang et al., 2004;van Huissteden et al., 2006).These models focus on CH 4 only, and assume that methanogenesis and methanotrophy occur in anoxic and oxic zones, respectively, which are spatially separated by the position of water table.In contrast, the two-substance model considers CH 4 and O 2 simultaneously, and the methanogenesis and methanotrophy occur according to the status of both gases in soils (e.g., Arah and Kirk, 2000).This is accomplished by introducing the inhibition effect of O 2 on CH 4 production and the competition for O 2 between heterotrophic respiration and methanotrophy.As such, CH 4 oxidation and heterotrophic respiration dominate in the oxic zone while CH 4 production dominates in the anoxic zone.The twosubstance models have been used in modeling CH 4 emissions from rice paddies (Matthews et al., 2000), and showed reasonable results compared with field measurements.Other existing models are conceptually of either one-substance or two-substance model structure (e.g., Potter, 1997;Zhang et al., 2002).
In biogeochemistry models, three pathways for gas transport are considered: (1) molecular diffusion, (2) plant-aided transport and (3) ebullition, though some models lump the three pathways together (e.g., Cao et al., 1995;Sass et al., 2000;Zhang et al., 2002).Ebullition, if considered explicitly, is often modeled as a threshold phenomenon using the Heaviside function with some universally prescribed threshold concentration of the dissolved gases (e.g., Walter and Heimann, 2000;Matthews et al., 2000).Field and analytical studies suggest such a simple algorithm does not fully represent the physical processes of ebullition (Bazhin, 2001(Bazhin, , 2004;;Baird et al., 2004;Tokida et al., 2005Tokida et al., , 2007)).Specifically, several factors have not been considered in the concentrationbased threshold function algorithms: (1) the composition of the bubbles affected by multiple substances such as CO 2 and N 2 ; (2) the effects of the hydrostacy affected by water table dynamics and atmospheric pressure variation (Bazhin, 2001;Tokida et al., 2005Tokida et al., , 2007) ) and (3) the ebullition threshold defined in terms of gas volumes is fuzzy rather than deterministically predictable because of possible re-dissolution and gas entrapping, during the course of ebullition (Martens and Klump, 1980;Kellner et al., 2006;Coulthard et al., 2009).
In this study, we revise the CH 4 module in a biogeochemistry model, the Terrestrial Ecosystem Model (TEM) (Zhuang et al., 2004) by incorporating the effects of multiple substances in a soil profile and a probabilistic pressurebased algorithm for ebullition.We apply the revised model to two temperate peatland ecosystems to demonstrate the importance of considering the effects of multiple substances in soils on episodic emissions during atmospheric pressure changes (Mattson and Likens, 1990).We also demonstrate the retardation effects of increases in standing water depth on CH 4 effluxes when different model complexities are assumed (Jauhiainen et al., 2005;Zona et al., 2009).

Overview
We developed a four-substance CH 4 module within a biogeochemistry model, the Terrestrial Ecosystem Model (Zhuang et al., 2004).The model was calibrated and applied to data from two temperate peatland sites in Michigan to demonstrate the capabilities of models with different complexities in simulating CH 4 effluxes.A group of sensitivity analyses were conducted to assess the need for a four-substance model with an improved ebullition algorithm.

The revised CH 4 module
The governing equation for a non-adsorbed substrate in a soil column (Fig. 1) is: where and z wt (unit: m) is the water table depth, being negative when it is above the soil surface.For substance i, the bulk concentration y i (unit: mol m −3 ) is related to its aqueous concentration y i,w (unit: mol m −3 water) and gaseous concentration y i,a (unit: mol m −3 air) through where (z,t) (unit: m 3 air m −3 soil) is air-filled porosity, α i is the Bunsen coefficient for gas i (see Appendix A for its calculation) and θ (z,t) (unit: m 3 water m −3 soil) is the volumetric soil moisture.The boundary conditions for Eq. ( 1) are y 0 (t) = y(0,t) for volatiles (4)

Chemistry involved in methane production and consumption
In wetland ecosystems, CH 4 is produced primarily through methanogenesis and consumed through methanotrophy Methanogenesis can proceed in either of the two pathways (Conrad, 1989) (Conrad, 1989).If one makes the assumption that dissolved oxygen is negligible in the aqueous phase, then the water table serves as a boundary in the soil between the oxic zone above the water table and the anoxic zone below the water table.Consequently, the equation to solve for the wetland CH 4 profile in a soil column is reduced to a system of a single substance, i.e.CH 4 , only.Such was adopted in Walter et al. (2001) and Zhuang et al. (2004), where the bubbling was modeled as a switch-on and -off process with a prescribed threshold CH 4 concentration, CH 4,max (unit: mol m −3 water).
If one considers the competition for O 2 in CH 4 oxidation and respiration processes, a third stoichiometry is involved: With such, we obtained a two-substance model considering both CH 4 and O 2 in a soil profile.Characterization of the aerobic and anaerobic zone in a soil column by the water table in the one-substance system is now revised by introducing the inhibition of O 2 on CH 4 production where P * CH 4 (unit: mol m −3 s −1 ) is the maximum CH 4 production potential when the environment is completely anoxic, and η (unit: m 3 water mol −1 ) is a parameter representing the sensitivity of methanogenesis to the concentration of dissolved oxygen y O 2 ,w in pore water.A value of 400 m 3 water mol −1 from Arah and Kirk (2000) was used for η.P * CH 4 is defined in Appendix B. Accordingly, the methanotrophy is restricted by the availability of O 2 as where Q * CH 4 (unit: mol m −3 s −1 ) is the oxidation potential when aqueous O 2 and CH 4 are not limited, and k CH 4 and k O 2 are Michaelis-Menten constants (unit: mol m −3 water) for CH 4 and O 2 .We use values of 0.44 mol m −3 water and 0.33 mol m −3 water, respectively, for k CH 4 and k O 2 (Arah and Kirk, 2000).Q * CH 4 is defined in Appendix B. The consumption of O 2 due to heterotrophic respiration and CH 4 oxidation is modeled as where V * R is the maximum rate of respiration when O 2 is not the limiting factor, k R is the Michaelis-Menten constant, using a value of 0.22 mol m −3 water (Arah and Kirk, 2000).As in Matthews et al. (2000), we assumed only the process of heterotrophic respiration competes with the process of methanotrophy for O 2 , thus V * R is twice that of P * CH 4 .We also neglected the O 2 consumption by electron acceptor reoxidation (Segers and Leffelaar, 2001;van Bodegom et al., 2001)   Carbon dioxide is produced in methanogenesis, methanotrophy and aerobic respiration: Just as for O 2 consumption, CO 2 production from electron acceptor reduction (Conrad, 1989) was also neglected here.
In the soil, consumption of CO 2 is zero, therefore, Q CO 2 = 0.
For N 2 , we assumed no production and consumption in the soil profile, therefore, P N 2 = Q N 2 = 0.

The pressure-based ebullition algorithm
We revised TEM to consider effects of hydrostacy on ebullition.Tokida et al. (2007) observed an abrupt change in the CH 4 emission rates associated with a decreasing atmospheric pressure, and the mixing ratio of CH 4 in the gas bubbles was no more than 50% (see their Fig. 2).Zona et al. (2009) found that, when the surface standing water increased, the CH 4 efflux was effectively retarded.Such behavior has not been explicitly considered and modeled in the process-based CH 4 models with the conventional algorithms of ebullition using a prescribed threshold of water dissolved CH 4 (e.g., Walter and Heimann, 2000;Arah and Kirk, 2000;Zhuang et al., 2004).Tokida et al. (2007) suggested a three-substance system, including CH 4 , CO 2 and N 2 , should be used to model the ebullition.Indeed, Bazhin (2001Bazhin ( , 2004) ) suggested that ebullition is triggered at a certain depth when the total pressure of the water-dissolved gases exceeds the hydrostatic pressure imposed at that depth by the water table and atmospheric pressure.Therefore, the simple concentrationbased threshold approach was replaced by an equation of hydrostatic equilibrium.In this study, we considered a foursubstance system, i.e.CH 4 , O 2 , CO 2 and N 2 , and ignored other possible trace gases (e.g.argon and hydrogen).We formulated the bubbling criterion (Fig. 1) as where P si is the partial pressure and H i (see Appendix A for the formula) is the Henry's law constant for gas i, p (= p/P 0 ) is the scaled atmospheric pressure, P 0 = 10 5 Pa, z 0 = 10 m, and z d = min(z − z s ,z − z wt ), and where z s is the depth of soil surface, set to 0.0; and head (unit: Pa) is the total hydrostatic pressure head imposed by atmosphere and water above depth z.The second equation in Eq. ( 15) accounts for the capillary force by considering the rate of saturation of the soil.Further we assumed bubbling only occurs below the water table, thus z d is always nonnegative.Note, in Eq. ( 14), we did not consider the effects of bubble shapes and number of bubbles, which would impact the surface tension between the bubble and water interface and consequently the bubbling criterion (Peck, 1960).We also neglected the change of water distribution caused by the ebullition (Rosenberry et al., 2003), which would cause some bubbles to be trapped and released later.With Eq. ( 14), the potential ebullition for gas i at a certain depth z is computed as where the equilibrium bulk concentration is and the equilibrium aqueous concentration is ỹi,w = head and δ(s) is the Dirac delta function.The potential ebullition computed from Eq. ( 16) can either be positive or negative, with positive implying bubble formation, and negative implying potential bubble re-dissolution.
To partly account for the fact that a fraction of the bubbles could be re-dissolved during their travel to the atmosphere (e.g., Martens and Klump, 1980), we used the algorithm in Fig. 2 to compute the ebullition.In this probabilistic algorithm, the possibility (p r ) of re-dissolution is proportional to the potential fraction of re-dissolution abs(E b (z))/[abs(E b (z)) + E(z)] (if E b (z) is negative calculated from Eq. 16).We modeled the re-dissolution as a yes/no process, if the random number u drawn from a uniform distribution U [0,1] is less than p r , dissolve the bubbles with an amount of abs(E b (z)) otherwise, bubbles continue moving upward without redissolution, and combine with possible bubbles generated at upper layers.
The algorithm was applied starting from the bottom of soil column to the level of water table.The total ebullition E is either released directly to the atmosphere or added to the soil column, depending on the location of water table.When the water table is at or above the soil surface, the gases carried in the bubbles are directly emitted to the atmosphere; otherwise, they are added to the soil layer right above the water table.
There is an alternative way to implement the above ebullition algorithm, i.e. using the volumetric criteria, such as in Kellner et al. (2006); Granberg et al. (2001), and most recently in Wania et al. (2010).Using the ideal gas law, the volume of substances in gaseous phase in equilibrium with the aqueous phase can be computed at all depths.The gas volumes are then compared with some predefined threshold to trigger the bubbles.However, such a threshold is fuzzy and varies temporally and spatially due to a group of different factors (Baird et al., 2004;Kellner et al., 2006).Our implementation relates the ebullition directly to the pressure.As such, the ebullition criteria can be determined physically using the available information on gas content and soil water elevations.Also, our algorithm does not need to make any assumption of the relative fractions of different gases in the bubbles (Kellner et al., 2006).Arguably, ebullition can even occur without the existence of CH 4 , as long as the buoyancy is greater than the weight of the bubble.There are processes that have not been accounted for in the algorithm, e.g.entrapped gas due to a wetting process from the soil surface down into the column, which could cause bubble formation (Kellner et al., 2006).Solutions should be found in future studies to address such events.It is likely that our algorithm will not always give superior results to that obtained using the volume-threshold-based method in other studies (e.g., Wania et al., 2010, and comparison is needed).However, ease of implementation will favor inclusion of this approach for other gases in future modeling.

Other transport routes and model implementations
We revised the pathways of diffusion and plant-aided transport in Zhuang et al. (2004) (see Appendix C for details).These and other processes described in previous sections gave the governing equations for CH 4 , CO 2 , O 2 and N 2 involved in the four-substance model in Appendix C. As a result, the net CH 4 efflux was computed where P OX is set to 0.5 for the one-substance model (Walter and Heimann, 2000) and 0.0, otherwise.
We used the mass balance approach to calculate the diffusive flux to avoid the ambiguity in choosing the depth for computation (Rothfuss and Conrad, 1998).The governing equation Eq. ( 1) was solved using the method of lines (Schiesser, 1991) with a first order implicit projectorcorrector method for the reaction terms.The integration was done with a time step of 2400 s.The soil column was approximated to a depth of 4 m with an exponentially stretching grid (a total of 40 nodes) that has finer grid resolution at the top and coarser grid resolution at the bottom (Oleson et al., 2004).
The revised TEM CH 4 module has three different levels of complexity: the one-substance model (S1 model hereafter) was obtained by (1) retaining the processes of methanogenesis and methanotrophy, (2) excluding processes involving other traces gases and (3) modeling ebullition with the conventional algorithm using a prescribed threshold CH 4,max equal to 1.31 mol m −3 water (at 25 • C); similarly, the twosubstance model (S2 model hence after) was obtained by considering CH 4 and O 2 simultaneously and modeling the ebullition with the concentration-based threshold approach, where O 2,max equal to 1.23 mol m −3 water (at 25 • C) and CH 4,max equal to 1.31 mol m −3 water (at 25 • C); when four gases were considered and ebullition was modeled with the new probabilistic pressure-based algorithm, a four-substance model (S4 model hereafter) was obtained.

Study sites
Two temperate peatlands located in southern Michigan on the Edwin S. George Reserve, a University of Michigan field station were used to test our revised CH 4 module.Three years of measurements from 1991 to 1993 were taken at Buck Hollow Bog and Big Cassandra Bog (42 • 27 N, 84 • 1 W).Buck Hollow Bog is an open peatland covered by a wet lawn of Sphagnum species, with a dense cover of Scheuchzeria palustris, an arrow-grass.Three flux chambers were grouped in a triangular pattern approximately 10 m apart to measure the net CH 4 flux in the Buck Hollow Bog.Measurements were taken at four sites at the Big Cassandra Bog.Sites 2 and 3 were used in this study, because these two sites are similar in terms of ecosystem conditions.The sites at the Big Cassandra Bog are dominated by Sphagnum and Polytrichum mosses and are covered by a dense stand of Chamaedaphne calyculata.Measurements of net CH 4 fluxes were made using static chambers, and gas samples were collected and analyzed within 3-4 days of collection on a Shimadzu GC-14A gas chromatograph with a flame ionization detector.Environmental variables including water table depth and soil temperature from 5 cm above the peat surface to 100 cm within the peat column were monitored and used as driving data for the model in this study.Available pore water concentration profiles (with a detection limit 0.1 µM) were also used in our assessment of the model.For a detailed description of the study sites and assessment of measurements, readers can refer to Shannon and White (1994).

Standard simulations
Standard simulations were conducted to evaluate the performances of the different models at both sites.The different model formulations were considered to include individual parameters whose values were calibrated in a modelspecific way by trial and error.Specifically, we modified the parameter values to match the simulated fluxes and pore water concentrations as closely as possible to the measurements (Table 1), so that the differences between different model simulations were mainly due to different model formulations.We used the measured water table depth and soil temperature as environmental forcing.Since no site-specific measurements of atmospheric pressure were available, we simply set total pressure to 1 atm, a standard value that has been used in other model studies (Walter and Heimann, 2000;Zhuang et al., 2004).For soil porosity, we assumed a value of 0.83 v v −1 for depths shallower than 0.5 m, linearly decreasing to 0.53 v v −1 at 0.9 m, and constant at 0.53 v v −1 to the lower boundary of 4 m.The scaled NPP data (Fig. 3) required to model CH 4 production were derived from simulations using TEM driven with monthly climate data (Mitchell et al., 2004).For the atmospheric mixing ratio of the gases involved, we assume 0.209 v v −1 for O 2 , 0.781 v v −1 for N 2 , 385 ppmv for CO 2 and 1740 ppbv for CH 4 (Forster et al., 2007).

Model sensitivity studies
To test the responses of the different models to water table dynamics, we ran our models with the time series of water table depth artificially increased or decreased over a specific time period during the emission season.As shown in Fig. 4, for the Buck Hollow site, we increased the water table by 10 cm between ordinal day 150 and 160 (1 January 1991 was set to ordinal day 1), and decreased by 10 cm between ordinal day 590 and 600; for the Big Cassandra site, we increased the water table depth by 10 cm between ordinal day 135 and 145, and decreased by 10 cm between ordinal day 689 and 699.These days were chosen such that the differences of the CH 4 effluxes between the standard simulations and sensitivity simulations were significant enough to be identified.The results were analyzed by comparing simulations with those from standard simulations.
Two sets of experiments were used to test the model response to atmospheric pressure change.First, we conducted sensitivity simulations using a time series of artificially perturbed low pressure or high pressure events (specifically 930 hPa low and 1045 hPa high) on two arbitrarily chosen days during the high-emission season in summer (Fig. 5).These two days were chosen based on the same criteria as that used in the water table sensitivity study.The effect of changing atmospheric pressure was analyzed by comparing the change in pathways of CH 4 transport with that from the standard simulation.This was used to analyze whether the response is physically consistent or not.A second test was carried out to evaluate the overall effect of atmospheric pressure variability using a time series of atmospheric pressure (Fig. 6) extracted from the European Centre for Medium-Range Weather Forecasts (ECMWF) Interim reanalysis dataset at the grid that encompassed the site for the same time period of measurement.The response was again analyzed by comparing the results to the standard simulations.

Comparisons between standard model simulations and site-level observations
All models resulted in similar CH 4 fluxes (Table 2).Specifically, for the Buck Hollow site, the S1, S2 and S4 models all captured the temporal variability of the CH 4 fluxes (Fig. 7a).Because of the probabilistic feature of the S4 model, mulitple runs were conducted.The differences from these multiple runs were indistinguishable because the emission routes at these two sites were dominated by plantaided transport, and bubble redissolution rarely occured.mean from an ensemble simulation of size four was shown for comparison.The S4 model performed best in terms of linear fitting and the root mean square error (RMSE) against the measurements (Table 2).The S1 model presented the second best results, with the simplest model structure.
Over the three-year period, the simulated mean daily fluxes were 107.00 mg CH 4 m −2 d −1 , 168.82 mg CH 4 m −2 d −1 and 114.43 mg CH 4 m −2 d −1 , respectively by the S1, S2 and S4 model.The S2 model gave almost 50% higher emission than did the S1 and S4 model.In the S1 simulation, the maximum production rate was at least an order of magnitude smaller than those of S2 and S4 models because it included no inhibition effect of O 2 on methanogenesis (Table 1).The maximum oxidation rates ( ÔCH 4 ) were similar in magnitude.
Note that the simulations (except S4) were not very sensitive to ÔCH 4 .The importance of different CH 4 transport pathways varied in different models (Fig. 7a).Diffusion played a significant role in the release of CH 4 into the atmosphere, even though plant transport remained the dominant pathway in S4 simulations.In S2 simulations, plant transport accounted for more than 90% of the efflux into the atmosphere.The modeled diffusion contribution to efflux reached a short-term maximum between day 170 and 180 in the S1 simulations.This was due to a sudden and substantial decrease in water table from above the peat surface to below the peat surface (Fig. 4).This increased the concentration gradient of the CH 4 near the peat surface and thus increased diffusion.Such short-term changes were also found in the S2 and S4 simulations, but with smaller magnitudes.All the models suggested that the efflux through ebullition was small, which agrees with measurements (Shannon and White, 1994).In the S4 simulation, the ebullition played a larger role than it did in S1 and S2 simulations.
S1 performed the best in simulating pore water concentrations; followed by the S2 model and then the S4 model (Fig. 7a).For 12 June 1993, none of the models presented a satisfying result.The discrepancy might be due to the nonlinearities of the transient simulations.For example, the pore water concentrations at a certain day were impacted by results in previous days.Uncertainty in the driving data (e.g.soil temperature) is another source of such discrepancy.
At the Big Cassandra site, the simulated effluxes were less satisfying than the Buck Hollow Bog site, though S4 model performed the best, followed by the S2 model and then S1 www.biogeosciences.net/7/3817/2010/Biogeosciences, 7, 3817-3837, 2010  model.Particularly for the second half year of 1991 and the year 1992, the results compared poorly with measurements.This underperformance may be due to an inadequate representation of the methanogenesis substrate, which was simulated in TEM but has not been specifically calibrated for wetland ecosystems.For the three-year period, the mean daily fluxes were 42.51 mg CH 4 m −2 d −1 , 50.83 mg CH 4 m −2 d −1 and 45.51 mg CH 4 m −2 d −1 , respectively, by the S1, S2 and S4 models.Still, S2 simulated the highest CH 4 emission among the three models.At this site, plant transport was the most important pathway, and ebullition was relatively unimportant.All three models presented similar results for pore water concentrations (Fig. 7a).

Sensitivity to water table change
We found that the response to a 10 cm change in surface standing water caused a change in the CH 4 efflux by as much as −50 ∼ 300 mg CH 4 m −2 d −1 (Fig. 8).Such responses depend on site characteristics and model complexity.At the Big Cassandra site, models S1 and S2 yielded a stronger response in effluxes to the water table increase than did the S4 model, while the opposite occurred for the water table decrease.At the Buck Hollow site, the responses of different models to the water table change were similar, but S1 gave a much stronger response to the water table decrease than did models S2 and S4.Nevertheless, all the models successfully predicted the retardation effect of an increase in surface standing water on CH 4 efflux (Fig. 9).This can be explained by the low diffusivity of gases in water relative to air.A higher column of surface standing water represents a longer distance of diffusion before gas can escape into the atmosphere.This phenomenon may account for CH 4 accumulation in the peat column, which in turn enhances plant mediated transport and the oxidation of CH 4 in the rhizosphere, which decreases the efflux.When water table depth decreases, the diffusion distance is reduced, and efflux to the atmosphere increases (Fig. 9).In places where emergent vascular plants are sparse, a decrease in water table depth could enhance CH 4 efflux through ebullition.This was tested by re-moving the plants and using the remaining parameters from Table 2 in the simulations.The results are shown in Fig. 10.We found that the S1 and S2 models had negligible ebullition compared to that from model S4 in response to changes in water table.The burst of ebullition predicted by model S4 was more significant at the onset of a drop in standing water level.Then it decreased as the water table continuously dropped and finally reached zero ebullition when the CH 4 accumulation was too low to support ebullition.In contrast, the S1 and S2 models predicted that diffusion was the major pathway of CH 4 efflux and greatly underestimated the CH 4 emissions through ebullition.Field data from Buck Hollow site also supported ebullition as important when the vegetation is sparser (Shannon et al., 1996).These findings suggest that while models S1 and S2 performed relatively well with careful calibrations, the positive results were fortuitous.They represent an inadequate formulation but an artful parameterization of the problem.

Sensitivity to atmospheric pressure change
For the first sensitivity test of the models to atmospheric pressure change, S1 and S2 models showed much weaker (10 1 ∼ 10 2 weaker) responses than that of the S4 model (Fig. 11).In both S1 and S2 models, a change in atmospheric pressure can only affect the atmospheric concentration of the gases.Given the small feasible range of atmospheric temperature change and atmospheric pressure change, the change in   atmospheric CH 4 concentration is small, implying a small change of the diffusion rate and thus the methane efflux in these two models.However, in the S4 model, the atmospheric pressure was further related to the ebullition fluxes.When a low atmospheric pressure occurred, the ebullition criterion became less restrictive, and bubbles were more easily formed, enhancing the ebullition.In this simulation, a decrease in atmospheric pressure could trigger an episodic increase in CH 4 efflux by as much as 120 mg CH 4 m −2 d −1 at the Buck Hollow site and as much as 80 mg CH 4 m −2 d −1 at the Big Cassandra site, comparable to the enhancement due to a decrease in surface standing water table depth.
Results from the sensitivity tests of transient atmospheric pressure were analyzed, which we found were very different with the different models (see Fig. 12).For instance, in cases of low atmospheric pressure events at the Buck Hollow site, the S4 model usually predicted higher fluxes through enhancement of ebullition.For the S1 model, response to atmospheric pressure change was negligible.The S2 model also responded significantly to the change of atmospheric pressure, but showed lower fluxes in accordance with a lower concentration of atmospheric CH 4 at the upper boundary.In cases of high pressure events, S4 yielded reduced the fluxes by suppressing ebullition, S2 yielded enhanced fluxes due to a higher atmospheric concentration of CH 4 at the upper boundary, and S1 showed little response.Similar results were found at the Big Cassandra site.The change of atmospheric pressure also changed the rate of plant aided transport and diffusion.However, in our formulation of the algorithm, ebullition is the preferred route if it is triggered (which we also believe is true in the field, e.g.Tokida et al., 2005).By analyzing the cumulative differences, we found for the threeyear period at the Buck Hollow site that the S4 model predicted around 5% more emitted CH 4 using the transient atmospheric pressure data than it did using the standard 1 atm pressure.The S2 model, in contrast, predicted around 3% www.biogeosciences.net/7/3817/2010/Biogeosciences, 7, 3817-3837, 2010 less emission when the transient atmospheric pressure was used.Similar results were found for the Big Cassandra site, with smaller differences, in accordance with the lower emission rates.

The importance of using temperature dependent parameters
In our standard and sensitivity simulations with the revised CH 4 module, we treated the phase-related parameters, such as diffusivities, Henry's law constants and Bunsen coefficients as temperature dependent.Analysis showed that, within the typical temperature range (e.g., −5 to 30 • C), the diffusivity of CO 2 in water changed ± 5%, the Henry's law constant and Bunsen coefficient changed ± 50% (results not shown).For CH 4 , its diffusivity in air changed ± 10% and in water changed ± 5%, but the Henry's law constant and Bunsen coefficient changed ± 40%.To test if fixing these parameters at a specific reference temperature could significantly affect the results, we conducted a set of simulations with the phase-related parameter values corresponding to a reference temperature of 12.5 • C. We found, for the S1 model, that the temperature dependence of these parameters did not change the efflux significantly (Fig. 13).However, compared to the standard simulation, both the S2 and S4 models predicted a higher CH 4 efflux during high emission periods.Given that soil temperatures at the two sites were often above the reference temperature during the high-emission summer season, the lower solubility computed with the reference temperature allowed less O 2 and CH 4 to be stored in the soil, given almost the same rate of CH 4 production.Further, considering the inhibition effect of O 2 on methanogenesis and the stimulus effect of O 2 on methanotrophy, the higher emissions associated with the S2 and S4 models using coefficient values at the reference temperature is then explained as a stimulus of gas transport to the atmosphere.We conclude that fixing phase-related parameters at their reference temperature values is safe for the S1 model, but that temperature-dependent parameters should be used in the S2 and S4 models.

The role of root density distribution
Previous CH 4 modeling has used a linear function for the vertical distribution of root density (e.g., Walter and Heimann, 2000;Zhuang et al., 2004; also see Eq. (D2) in Appendix D).However, root biomass is often found to be exponentially distributed (Jackson et al., 1996), and an exponential root distribution could also be used (e.g., van Huissteden et al., 2006).Here we implemented both linear and exponential root distribution functions in our revised CH 4 module to test if they make a difference in CH 4 effluxes (Tables 1-4) and Table 3. Parameters calibrated for model simulations using a linear root distribution function.
Table 4. Comparison of CH 4 efflux from the simulations using a linear root distribution function to the measurements at the two Michigan peatlands.All statistics were tested for significance and were found significant with p<0.001, except those denoted in the parentheses.pore water concentrations.In both configurations, the parameter values were obtained from calibration based on measurement data.The differences between the simulations were thus mainly due to the different model configurations.
We found the contributions from different pathways were different when two different root density distribution functions were used (Figs.7a and 14a).In the standard simulation using an exponential root density distribution function, ebullition played a minor role, whereas in the simulation employing a linear root density distribution, the ebullition was more significant, particularly in the S2 and S4 models.The ebullition was enhanced more at the Big Cassandra site than at the Buck Hollow site, suggesting that the responses of the models to two different root distribution functions are site dependent and related to the net CH 4 production characteristics of the site.In the simulation employing the linear root distribution function, model-derived pore water CH 4 concentration profiles underestimated field observations in the upper level of the peat column (Fig. 14a).The lower CH 4 concentrations ware due to a poor representation of transport in the upper portion of the soil column when a linear root density function was used.Conversely, the exponential root distribution extends smoothly down into depth of the peat column, producing more realistic pore water CH 4 concentration profiles.Therefore, though rigorous parameterization can lead to a good fit of the modeled CH 4 fluxes with respect to the measurements, the model fails to capture other aspects of the measurements when an improper formulation of the problem is used.Also, in our case, the exponential distribution is a superior representation of root density as a function of depth.

Issues for regional application of the different CH 4 models
The CH 4 models of different complexities developed in this study can be used for regional hindcast and projection of wetland CH 4 emissions provided that necessary climate forcing data are available.This is not a problem when these CH 4 models are used inside a biogeochemistry model, such as TEM, where the necessary climate forcing data to the CH 4 models can be computed explicitly when the biogeochemistry model is driven by a climate dataset including air temperature, cloud fraction, precipitation and vapor pressure (Zhuang et al., 2004).The three CH 4 models have almost the same requirement for climate forcing, except that the S4 model requires surface pressure data for a better performance.For historical simulations, the surface pressure data can easily be obtained from various climate data sources, e.g.datasets from NCEP Reanalysis and ECMWF Interim Reanalysis.For projections, GCM model outputs would be a source for the necessary climate data.In some cases, when sea level pressure data rather than surface pressure data are output from GCMs (e.g.models involved in IPCC AR4, http://www.ipcc-data.org/ar4/gcmdata.html).The surface pressure data can then be derived from a combination of sea level pressure data, information of air temperature, elevation and vapor pressure (Wallace, 2006).Although the CH 4 models developed here have different complexities, they have almost the same number of parameters that require calibration.The more complicated S2 and S4 models have even fewer parameters to calibrate.For instance, the S2 and S4 models compute the fraction of CH 4 oxidation in the rhizosphere explicitly -no parameterization is needed as in the S1 model.
The parameters of the CH 4 models should be handled carefully in regional applications.For instance, upscaling maximum CH 4 production potential ( PCH 4 ) and maximum CH 4 oxidation potential ( QCH 4 ) from the calibrated sites to a region is critical.Currently, we use the maximum monthly NPP derived from a 50-year historical TEM simulation to scale the parameter PCH 4 and the maximum monthly soil respiration to scale the parameter QCH 4 .Both NPP and soil respiration are simulated with TEM.The extrapolation is based on the fact that CH 4 productivity is usually positively related with NPP (e.g.Chanton et al., 1995), and CH 4 oxidation is positively related with respiration (e.g.Nakano et al., 2004).The scaling is based upon the vegetation cover data.The remaining model parameters derived from the calibrated site are used for our regional extrapolations.Thus, as a next step, we will test how different ways in extrapolating the sitespecific parameters to a region affect the uncertainties in the wetland CH 4 emissions quantified with the CH 4 models of different complexities.Also, an analysis of uncertainty due to equifinality will be attempted to investigate robustness of the parameterization from calibration at the limited number of sites (Tang and Zhuang, 2008).
The regional water table dynamics are another major source of uncertainty in quantifying regional wetland CH 4 emissions.Standing water depth on top of soils is also essential to a proper quantification of regional CH 4 effluxes.
In particular, when the S4 model is used in regional simulations, there are grid cells, where vegetation is sparse, emitting CH 4 mainly via ebullition.In contrast, the S1 and S2 models greatly underestimate the CH 4 emissions in such cases (e.g.Fig. 10).In these simulations, water table depths play a significant role in affecting CH 4 production, oxidation, soil pressure profile, and diffusion process.To more accurately simulate water table dynamics, we are currently testing several different algorithms (e.g.Granberg et al., 1999;Weiss et al., 2006).The methane models with different complexities will be further coupled with existing soil physics models (e.g.Zhuang et al., 2001Zhuang et al., , 2003;;Tang and Zhuang, 2010) and with the tested water table depth model to conduct regional and global analyses of wetland CH 4 emissions.

Conclusions
We revised an extant process-based biogeochemistry model, the Terrestrial Ecosystem Model to account for the effects of multiple in a soil profile on CH 4 production, oxidation, and transport.The new development allows CH 4 effluxes to be modeled with different levels of model complexity.When (O 2 , N 2 , CO 2 and CH 4 ) are considered, the inhibitory effect of O 2 on CH 4 production and the stimulatory effect of O 2 on CH 4 oxidation are well accounted for, and ebullition is modeled in a physically logical manner.When ebullition is modeled with a concentrationbased threshold approach and the inhibition effect of O 2 on CH 4 production, and the competition for O 2 between methanotrophy and heterotrophic respiration are considered, the model becomes essentially a two-substance system.If we ignore the role of O 2 , while modeling bubble ebullition with the concentration-based threshold function, the model is reduced to a one-substance system.These models were tested through a group of sensitivity analyses at two temperate peatland sites in Michigan.We showed that only the foursubstance model with the new ebullition algorithm is able to account for the effects of a sudden drop in atmospheric pressure or in water table on episodic emissions.All models simulated the retardation of CH 4 efflux after an increase in surface standing water due to inhibited diffusion and enhanced rhizospheric oxidation.We conclud that, to more accurately account for the effects of atmospheric pressure dynamics and water table dynamics on methane effluxes, the four-substance model with the probabilistic but physics-based ebullition algorithm should be used in the future to quantify global wetland CH 4 emissions.Further, to more accurately simulate the pore water gas concentrations and different pathways of CH 4 transport, an exponential root distribution function should be used and the phase-related parameters should be treated as temperature dependent.
where f (S OM (z,t)), f (T (z,t)), f (pH(z,t)), f (Eh(z,t)) are multiplier functions of methanogenesis substrate availability (modeled as a function of scaled NPP), soil temperature, pH value and redox potential, as defined in Zhuang et al. (2004).PCH 4 (unit: mol m −3 s −1 ) is a scaling parameter for model calibration.Another site specific parameter that needs calibration is the Q 10 coefficient (P Q 10 ) of f (T (z,t)).
The maximum CH 4 oxidation potential is defined as where f (T (z,t)), f (θ(z,t)), f (Eh(z,t)) are functions of soil temperature, soil moisture and redox potential (see Zhuang et al., 2004 for detailed descriptions).Parameter ÔCH 4 (unit: mol m −3 s −1 ) is calibrated for every representative site.The Q 10 coefficient for temperature effect is set to 2 throughout this study.

Appendix C
For the diffusive flux, the diffusion constant in Eq. ( 1) is defined for the bulk medium, which is conventionally computed (Stephen et al., 1998) as where subscripts a and w denote the diffusivity in air and in water (see Appendix A for ways of computation).τ is the tortuosity factor in the soil, taken as 1.5 throughout the study (Arah and Stephen, 1998).
For gas transport through the aerenchyma of wetland plants, we, following the argument in other studies (Teal and Kanwisher, 1966;Matthews et al., 2000;Segers and Leffelaar, 2001), assumed the N 2 , CO 2 and CH 4 are transported in a similar way, such that R i = R * i (y i,a − y i,atm ) = λ r L v D i,a f (t)(y i,a − y i,atm ) (C2) where λ r (unit: m air (m root) −1 ) is the specific conductivity of the root system and L v (unit: m root m −3 soil) is the root length density.A value of 3.0 × 10 −4 was used for λ r .The vertical distribution of L v in soil is assumed following the Gale-Grigal model (Jackson et al., 1996) (see the exponential model in Appendix D).The temporal variation f (t) of the root is modeled similarly to Zhuang et al. (2004) and Walter and Heimann (2000).Also, we assumed the four gases can either be transported from the atmosphere to the roots or from the roots to the atmosphere.When the onesubstance model is switched on, the oxidation of CH 4 in the rhizosphere (Beckett et al., 2001) is not considered explicitly, rather, as in Walter and Heimann (2000), we assume 50% of CH 4 is oxidized. In and for N 2 is

Appendix D
The root length density in Eq. ( C2) is defined as where β is 0.943 for Buck Hollow Bog, and 0.910 for Big Cassandra Bog.R veg is a scaling parameter needed in calibration to account for differences in conducting capabilities for different plants.Note, the integrated root distribution function f (z) from lower boundary to soil surface equals one.The alternative root distribution used in Sect.3.4, is defined as where R d is root depth, computed using the Gale-Grigal model.This is different from the formula adopted in Walter and Heimann (2000) and Zhuang et al. (2004) in that we here imposed the constraint such that the vertical integration of f (z) equals one.

Fig. 1 .
Fig. 1.Diagram showing the three major gas-transport pathways involved in the multi-substance CH 4 model.See text for details about the calculations.

Fig. 2 .
Fig. 2. The probabilistic algorithm used in the S4 model to compute ebullition.

Fig. 3 .
Fig. 3.The time series of scaled NPP (scaled with the median peak NPP at the site from a 100-year simulations) used as driving data in this study.Same data were used at both sites, with truncation into proper time periods.

Fig. 4 .
Fig. 4. The time series of water table depth used as driving data in this study at: (a) Buck Hollow site; (b) Big Cassandra site.

Fig. 5 .
Fig. 5. Synthetic atmospheric pressure used in simulations for sensitivity analysis at: (a) Buck Hollow site; (b) Big Cassandra site.

Fig. 6 .
Fig.6.Transient atmospheric pressure used in simulations for sensitivity analysis.The same time series was applied at both sites, with truncation into proper time periods.

Fig. 7a .Fig. 7b .Fig. 7c .
Fig. 7a.Methane effluxes, component-wise emissions and pore water concentration profiles from one-substance model (S1 model) in the standard simulations.(a) Panels for the Buck Hollow site.(b) Panels for the Big Cassandra site.Dashed lines indicate the level of water tables.
Fig.8.Water table sensitivities with models of different complexities.Left-hand panels are for Buck Hollow site; right-hand panels are for Big Cassandra site.The red symbol "+" indicates the period of artificially increased water table depth, and the green symbol "−" indicates the period of artificially decreased water table depth.

Fig. 9 .
Fig.9.Zoom-in plots for the water table sensitivity with models of different complexities at the Buck Hollow site.Left-hand panels were for the period where the water table depth is increased by 10 cm compared to the measured data, right-hand panels are for the period where the water table depth is decreased by 10 cm.The arrows indicate when the specified water table change starts and ends.

Fig. 10 .
Fig. 10.Methane effluxes through ebullition at the Buck Hollow site when no vegetation is present to support plant aided transport.See text for details.

Fig. 11 .
Fig. 11.Atmospheric pressure sensitivities with models of different complexities.(a), (c) and (e) panels are for Buck Hollow site; (b), (d)and (f) panels are for Big Cassandra site.The red symbol "+" indicates the day of artificially perturbed high pressure event, and the green symbol "−" indicates the day of artificially perturbed low pressure event.Note the scale of results from S4 model is 10 2 of that from S2 model, and 10 3 of that from S1 model.

Fig. 12 .Fig. 13 .
Fig. 12. Differences in the response to transient atmospheric pressure and that to constant atmospheric pressure with models of different complexities.(a) Time series of the absolute difference at the Buck Hollow site.(b) Time series of cumulative difference at the Buck Hollow site.(c) Time series of the absolute difference at the Big Cassandra site.(d) Time series of cumulative difference at the Big Cassandra site.The difference is defined by subtracting the fluxes simulated with the standard 1 atm pressure from that using the transient atmospheric pressure data.

Fig. 14a .Fig. 14b .Fig. 14c .
Fig. 14a.Methane effluxes, component-wise emissions and pore water concentration profiles from one-substance model (S1 model) in test simulations using a linear root distribution function.(a) Panels for the Buck Hollow site.(b) Panels for the Big Cassandra site.Dashed lines indicate the level of water tables.
for the moment.Since no O 2 is produced in the soil, P O 2 is set to zero.

Table 1 .
Parameters calibrated for standard model simulations.

Table 2 .
Comparison of CH 4 efflux from the standard simulations to the measurements at the two Michigan peatlands.All statistics were tested for significance and were found significant with p < 0.001.
table sensitivity with models of different complexities at the Buck Hollow site.Left-hand panels were for the period where the water table depth is increased by 10 cm compared to the measured data, right-hand panels are for the period where the water table depth is decreased by 10 cm.The arrows indicate when the specified water table change starts and ends.