Articles | Volume 16, issue 3
Research article
04 Feb 2019
Research article |  | 04 Feb 2019

Modeling anaerobic soil organic carbon decomposition in Arctic polygon tundra: insights into soil geochemical influences on carbon mineralization

Jianqiu Zheng, Peter E. Thornton, Scott L. Painter, Baohua Gu, Stan D. Wullschleger, and David E. Graham

Rapid warming of Arctic ecosystems exposes soil organic matter (SOM) to accelerated microbial decomposition, potentially leading to increased emissions of carbon dioxide (CO2) and methane (CH4) that have a positive feedback on global warming. Current estimates of the magnitude and form of carbon emissions from Earth system models include significant uncertainties, partially due to the oversimplified representation of geochemical constraints on microbial decomposition. Here, we coupled modeling principles developed in different disciplines, including a thermodynamically based microbial growth model for methanogenesis and iron reduction, a pool-based model to represent upstream carbon transformations, and a humic ion-binding model for dynamic pH simulation to build a more versatile carbon decomposition model framework that can be applied to soils under varying redox conditions. This new model framework was parameterized and validated using synthesized anaerobic incubation data from permafrost-affected soils along a gradient of fine-scale thermal and hydrological variabilities across Arctic polygonal tundra. The model accurately simulated anaerobic CO2 production and its temperature sensitivity using data on labile carbon pools and fermentation rates as model constraints. CH4 production is strongly influenced by water content, pH, methanogen biomass, and presence of competing electron acceptors, resulting in high variability in its temperature sensitivity. This work provides new insights into the interactions of SOM pools, temperature increase, soil geochemical feedbacks, and resulting CO2 and CH4 production. The proposed anaerobic carbon decomposition framework presented here builds a mechanistic link between soil geochemistry and carbon mineralization, making it applicable over a wide range of soils under different environmental settings.

1 Introduction

The northern permafrost region contains 1400–1800 Pg soil carbon (C), which is more than twice as much C as is currently contained in the atmosphere (Tarnocai et al., 2009; McGuire et al., 2012). Persistent cold and saturated soil conditions have limited C decomposition in this reservoir. However, rapid warming and permafrost thaw expose previously frozen organic carbon to accelerated microbial decomposition, potentially leading to emissions of carbon dioxide (CO2) and methane (CH4) that have a positive feedback on global warming (Zimov et al., 2006; Schuur et al., 2009, 2015). How quickly frozen soil organic matter (SOM) will be mineralized, and how much permafrost C will be released to the atmosphere following thaw, is highly uncertain. Earth system models project 27–508 Pg carbon release from the permafrost zone by 2100 under current climate forcing (Zhuang et al., 2006; Koven et al., 2015; MacDougall et al., 2012; Schaefer et al., 2014). Understanding environmental dependencies of SOM decomposition is therefore essential for reducing model uncertainties and improving predictions of future climate change.

Disagreement in model projections for the northern permafrost region could be due to differences in model structure, model initialization, or parameters used in simulations. Despite increasingly detailed process representations in many models that simulate terrestrial CO2 and CH4 fluxes, important geochemical and metabolic constraints might still be poorly represented, oversimplified, or missing in current biogeochemical models (Xu et al., 2016). The northern permafrost region is rapidly changing in response to the changing climate. Rising temperatures not only release more labile carbon from permafrost for decomposition but also create thermal and hydrological heterogeneity that further affects biogeochemical processes. Here, we examine two mechanisms that substantially affect SOM turnover in permafrost-affected soils. First, rising temperature alters the kinetics of biogeochemical reactions (Segers, 1998). This effect is more pronounced at subzero temperature (Bore et al., 2017), and the process rate increase is higher at lower temperature ranges (Davidson and Janssens, 2006). Microbial communities also change with temperature, compounding effects on process rates (Karhu et al., 2014). Models address this temperature effect using empirical functions and parameters (Tuomi et al., 2008; Xu et al., 2016), which might be highly biased depending on model assumptions and original curve-fitting techniques, generating large uncertainties. Second, heterogeneity in permafrost thaw and related hydrological responses creates geochemical gradients in soils. Models use different levels of detail to simulate effects of water saturation (Meng et al., 2012; Xu et al., 2016). Soil moisture limits gas transport, and it is often used as an implicit control on heterotrophic respiration and methanogenesis. However, the explicit processes resulting from soil oxygen depletion (e.g., soil redox status and pH dynamics) are not widely represented (Riley et al., 2011; Meng et al., 2012; Xu et al., 2015).

The extent of SOM decomposition and gas emissions depends upon soil geochemical characteristics beyond temperature and O2 availability. Among the wide range of environmental variables, pH emerges as a primary control on decomposition by regulating both microbial communities and microbial metabolic activities (Zhalnina et al., 2015; Bethke et al., 2011; Jin and Kirk, 2018). pH affects microbial metabolism by modulating the thermodynamics and kinetics of redox reactions. Redox reactions produce or consume protons, and thus their free energy yields vary with pH (Bethke et al., 2011; Jin and Bethke, 2007). The Gibbs free energy available to anaerobic microorganisms that degrade simple organic molecules generally increases (becomes less favorable) with increasing pH (Bethke et al., 2011). Notably, iron [Fe(III)] reduction is highly proton consuming and becomes less favorable at higher pH (Fig. S1 in the Supplement). Previous studies identified iron reduction as a major process in anoxic Arctic soils (Lipson et al., 2010, 2013), which increases local pH and might favor co-occurring methanogenesis (Tang et al., 2016; Wagner et al., 2017). However, the influence of iron reduction on methanogenesis rates in different soils is rarely investigated. The reactivity of iron and its pH feedback impose additional complexity on the controls of SOM decomposition and associated CH4 production.

Despite the importance of pH in controlling redox reactions and resulting C emissions, pH change is not explicitly represented in biogeochemical models. Most of the current biogeochemical models apply a single initial pH value for redox reactions without considering proton production and consumption during the processes. Traditional decomposition models use landscape position, soil moisture content, or other proxies for O2 concentration to determine the form of C release. Scalars on aerobic respiration (Riley et al., 2011; Lawrence et al., 2015) or empirical ratios of CO2 and CH4 (Koven et al., 2015) are often used to inform the extent of C decomposition and partitioning of CO2 and CH4 production. Reactions that produce or consume protons and the resulting pH changes or ion exchange reactions are not considered in these empirical models. Some process-rich models explicitly include details of methanogen populations and their interactions with substrates and other environmental factors, but these models still lack the capability to simulate pH changes during long-term carbon decomposition. Instead, constant pH is often assumed within bell-shaped pH response functions (Meng et al., 2012; Tian et al., 2010; Xu et al., 2015). Without underlying proton exchange and pH buffering mechanisms, a significant error may occur when rate calculations depend heavily upon the initial choice of a single optimal pH value for various reactions.

In this study, we developed a new anaerobic carbon decomposition model framework with explicit representation of aqueous-phase geochemistry to allow pH and thermodynamic calculations. By coupling three different models, including a thermodynamically based microbial growth model, a substrate pool-based model, and a humic ion-binding model, we built a process-rich carbon decomposition model that allows simultaneous thermodynamic and pH calculations. Results from anoxic incubations of permafrost-affected soils along a hydrological gradient were synthesized to parameterize and validate this new model framework. The main objectives of this study were to (i) examine the role of soil geochemical variables in controlling anaerobic carbon decomposition and carbon release (as both CO2 and CH4); (ii) develop a common set of parameters in the new anaerobic carbon decomposition framework to capture variabilities in CO2 and CH4 production; and (iii) evaluate model uncertainties in responses to both soil heterogeneity and model parameterization, emphasizing effects of soil saturation, pH, and temperature response.

2 Materials and methods

2.1 Anaerobic carbon decomposition model

The anaerobic carbon decomposition framework was developed with explicit representation of fermentation, methanogenesis, and iron reduction, which were identified as key mechanisms for anaerobic CO2 and CH4 production (Roy Chowdhury et al., 2015; Yang et al., 2016; Zheng et al., 2018b). The main structure of this framework included two major components: a simplified Community Land Model with prognostic carbon and nitrogen (CLM-CN) decomposition cascade (converging trophic cascade, or CTC) (Thornton and Rosenbloom, 2005) to facilitate parameterization of the upstream carbon flow entering the aqueous-phase dissolved organic carbon (DOC) pool (Fig. 1, process 1), and an aqueous phase to facilitate calculations of thermodynamics and redox-reaction-associated acid–base chemistry. An empirical approach was used to represent non-aqueous-phase soil organic carbon (SOC) decomposition. Additionally, mechanistic representations of methanogenesis and iron reduction were developed in this work based on aqueous-phase thermodynamic calculations.

Figure 1Conceptual diagram showing key processes in the anaerobic carbon decomposition framework, beginning with plant material and coarse wood debris (CWD). The numbers indicate different processes. 1. SOM degradation from soil organic carbon pools with increasing turnover time produces dissolved organic carbon (DOC) and CO2. 2. Fermentation of DOC into organic acids, H2, and CO2. 3. Methanogenesis from organic acids or H2. 4. Fe(III) reduction from organic acids or H2. 5. Fe(OH)3 dissolution.


The simplified CTC cascade included four SOM pools to represent bulk SOC with different levels of complexity. Changes of these SOM pools followed modified first-order decay (see the Supplement for details). We modified the original respiration fraction (Thornton and Rosenbloom, 2005; Koven et al., 2013) into direct and indirect fractions. Thus, for each SOM pool, the direct respiration fraction represented CO2 lost as originally defined, while the indirect respiration fraction was labile C produced from bulk C entering the aqueous-phase carbon pool (DOC pool, Fig. 1).

The large biomolecules in the DOC pool went through multiple hydrolysis and fermentation steps to produce low-molecular-weight organic acids that would further respire into CO2 and CH4 (Boye et al., 2017; Zheng and Graham, 2018; Yang et al., 2016; Roy Chowdhury et al., 2015). Under anoxic conditions, hydrolysis of polysaccharides was considered the rate-limiting step for downstream methanogenesis (Glissmann and Conrad, 2002). Polysaccharide hydrolysis has a favorable free energy, due to increased entropy, but cannot be readily coupled to biological energy transduction outside of the cell. We previously measured a rapid decrease in reducing sugar concentrations in pore water during tundra soil incubations, which indicated that hydrolysis limits decomposition (Yang et al., 2016). At low temperatures (below 15 C), the microbial degradation of cellulose was considerably diminished, while other polymers, such as starch or proteins, were degraded much faster at low temperature, resulting in the accumulation of organic acids, primarily acetic, propionic, and butyric acids (Kotsyurbenko, 2005; Yang et al., 2016). These low-molecular-weight organic acids further fueled microbial mineralization reactions that led to production of CH4 and CO2. Given that most anaerobic lignocellulose degraders also fermented sugars following hydrolysis (Blumer-Schuette et al., 2014), we assumed the turnover of DOC into low-molecular-weight organic acids was a single lumped fermentation process (Fig. 1, process 2), in which labile DOC (C6H12O6) was fermented into acetate, H2, and CO2 (Appendix A, Reaction AR1, and Table S1 in the Supplement). This assumption gave a fixed stoichiometry ratio: one-third of the fermented carbon was oxidized to CO2.

Redox reactions including methanogenesis and iron reduction were represented using a thermodynamically based approach (Istok et al., 2010), with unique microbial growth kinetics incorporated into energy-yielding redox reactions. In this thermodynamically based approach, the growth equations of methanogens and iron reducers were derived from paired electron donor (acetate or H2) and electron acceptor half reactions and a biomass synthesis equation (Istok et al., 2010). Using a constant molecular formula as biomass (C5H7O2N) and ammonium (NH4+) as the nitrogen source for biosynthesis, we derived the growth equations for methanogenesis and iron reduction (Appendix A, Reactions AR2–AR5). Rate calculations followed the generalized Monod rate law (Jin and Bethke, 2007) with an additional thermodynamic factor representing the thermodynamic driving force. The thermodynamic factor f(G) is calculated using Eq. (1) (Jin and Bethke, 2003):

(1) f G = 1 - exp - - Δ G - m Δ G p χ R T ,

where ΔG (kJ (mol reaction)−1) is the free energy change of the redox reaction. ΔG depends on the standard Gibbs free energy change (ΔG0) and the concentrations of chemical species involved in the reaction. ΔGp is the phosphorylation potential, i.e., the energy required to synthesize adenosine triphosphate (ATP) from adenosine diphosphate (ADP) and dihydrogen phosphate in the cell's cytoplasm. ΔGp is about 45 kJ (mol ATP)−1 (Jin and Kirk, 2018). m is the number of ATP molecules synthesized per redox reaction. χ is the average stoichiometric number (Jin and Kirk, 2018), R is the gas constant (kJ mol−1 K−1), and T is the absolute temperate (K). The factor f(G) ranges from 0 to 1, where the reaction is thermodynamically favorable when f(G)>0.

Both methanogenesis and iron reduction contribute to pH change. Reactions such as ferrihydrite reduction substantially increase alkalinity (Appendix A, Reactions AR4, AR6). Furthermore, the solubility of CO2 and the composition of dissolved CO2 and bicarbonate vary significantly over typical soil pH values, affecting all C mineralization processes. In the organic-rich soils modeled here, SOM rather than minerals provides most buffering capacity. Therefore, we used the humic ion-binding model to describe pH buffering during carbon decomposition. A simplified parameterization of proton binding is available in the Windermere Humic Aqueous Model (WHAM; Tipping, 1994, 1998), which has been extensively calibrated to represent the acid–base chemistry of “average” humic and fulvic acids, and benchmarked with heterogeneous natural organic matter (Atalay et al., 2009). We adopted the WHAM parameterization to represent proton-binding characteristics (pH buffering) provided by SOM (Tang et al., 2016). Using representative binding constants provided by WHAM, the pH buffering capacity can be directly adjusted by altering the number of proton-binding sites, which is assumed to be linearly correlated with the total amount of SOM (see the Supplement for details).

2.2 Model implementation and initialization

The above model structure was implemented using the open-source geochemical program PHREEQC 3.0 (Charlton and Parkhurst, 2011) with a new database describing SOC decomposition cascade, redox reaction kinetics, and pH buffering (redox.dat, available at, last access: January 2019). This model assumed thermodynamic equilibrium of aqueous chemical speciation, mineral dissolution/precipitation, and ion sorption/desorption based on the updated PHREEQC thermodynamic database (phreeqc.dat; Charlton and Parkhurst, 2011). The database was modified to include WHAM pH buffering and reaction kinetics for SOM pools decay and reaction kinetics for fermentation, methanogenesis, and iron reduction. The kinetic rate constants and microbial biomass growth and decay rates were adopted from former thermodynamically based studies (Istok et al., 2010) and previously tested with low-center polygon Arctic soils (Tang et al., 2016).

The model initialization was based on both the incubation conditions and soil geochemical characterizations (Fig. S2). The initial partitioning of SOM pools was assumed to be at fixed ratios due to the limitation of short-term incubation data. Under the experimental conditions, we assumed SOM1 and SOM2 pools with relatively shorter turnover rates (τ=14 days and 70 days, respectively) were most relevant in the model. On the other hand, SOM3 and SOM4 pools were relatively inert (τ>2 years for both pools). We started with the relative fractions of SOM pools at approximately 10 %, 40 %, 10 %, and 40 % of SOC in organic and mineral soils for SOM1–4. We further assessed the bias of this assumption with sensitivity analysis. Sizes of SOM1 and SOM2 pools were reduced by 90 % for permafrost to better account for the overall low levels of carbon degradation.

All other variables required were initialized using measurements based upon 10 to 15 g of wet soil incubated in 60 to 70 mL sealed bottles. Total SOC, total water (TOTW), total organic acid carbon (TOAC), pH, and the initial concentration of Fe(II) were specified in the model based on measurements (Table S2). The DOC pool in the model was initialized using the measured water-extractable organic carbon (WEOC) expressed as a fraction of SOC (fdoc). On average, WEOC accounts for approximately 2 % of SOC based on our synthesized data (see Sect. 3.1). This value is consistent with previous long-term incubations, which suggested less than 5 % of SOC was fast-decomposing carbon in permafrost-affected soils at a standardized temperature of 4–5 C (Knoblauch et al., 2013; Schädel et al., 2014). The starting biomass of methanogens and iron reducers was assumed to be within the range of 10−3 to 10−5 g C g SOC−1 for organic soils, 10−4 to 10−6 g C g SOC−1 for mineral soils, and 10−7 to 10−9 g C g SOC−1 for permafrost (Table S2). This stratified microbial biomass distribution was used to represent the vertical gradient in the relative abundance of microbial communities (Yang et al., 2017).

The lumped fermentation process was the rate-limiting step in the model and was fitted individually with data from each soil microcosm. Based upon reaction stoichiometry, the fermentative conversion of each mole of labile C led to two-thirds of mole of organic acids and one-third of a mole of CO2. Organic acids were mineralized via methanogenesis or iron reduction to convert approximately 49 % to 88 % of C in organic acids into CO2. This estimation was based on reaction stoichiometry of Reactions (AR2) and (AR4), and a fraction of the C was incorporated into microbial biomass. Therefore, the percentage of respired C would be less than 100 % even if all organic acids were respired as CO2. If we assume all fermentation products were mineralized into CH4 and CO2, we could estimate the fermentation rate (Rfer) from measured CO2 production. Thus, Rfer was estimated using the initial CO2 production rate in the incubation data,and further optimized using the least squares method by fitting with observed CO2 production values (Table S3).

Temperature and pH response functions were used to further constrain model simulations (Fig. S2). A temperature effect was parameterized using the CLM-CN temperature response function (Appendix B, Eq. B1). Additional temperature response functions were evaluated by sensitivity analysis (see Sect. 2.4). The effect of pH on biological reaction rates is modulated by bell-shaped pH response functions (Tang et al., 2016; Xu et al., 2016). Here, we used the Dynamic Land Ecosystem Model (DLEM) pH response function (Appendix B, Eq. B5), since it generated the least variation in parameter perturbation tests (Tang et al., 2016).

2.3 Incubation data synthesis for model validation

Incubation data from Utqiaġvik (Barrow) Alaska soil cores that represent the microtopographic heterogeneity of polygonal tundra were synthesized to validate the new anaerobic carbon decomposition model. The selected datasets represent fine-scale variabilities in thermal and hydrological regimes across the gradient of soil microtopographic positions (Herndon et al., 2015). The synthesized data contain complete sets of soil geochemical descriptions for organic, mineral, transition zone (if identified), and permafrost layers from each microtopographic feature (see the Supplement for details). Levels of total soil organic carbon, WEOC, and TOAC were available before and after soil incubation. Besides CH4 and CO2 production during low-temperature soil decomposition, data on Fe(II) concentrations and pH changes were also available for model initialization and validation.

2.4 Model parameter uncertainty

This model was designed as a generic framework to simulate anaerobic carbon decomposition across a range of soil physiochemical conditions. Two types of sensitivity analysis were conducted to evaluate model performance. First, possible bias and variations associated with model initialization variables (soil geochemical attributes) were assessed using perturbation simulations. Variations of ±25 % and ±50 % (+100 % and 200 % for some variables) were applied to these variables, and the resulting changes in cumulative CO2 and CH4 production were evaluated by comparison with reference simulations. This evaluation helps to identify critical measurements needed for initializing the model. Second, parameters specifically benchmarked in this study and parameters adopted from empirical relationships were also evaluated with perturbation simulations. This test helps to apportion the model prediction uncertainties into different sources, including model input, parameters, or model structure.

3 Results

3.1 Meta-analysis to validate model assumptions

Incubation data used in this study were generated from soils representing different microtopographic features with a wide range of moisture and SOC contents and reported elsewhere (Roy Chowdhury et al., 2015; Zheng et al., 2018b). Correlation analysis revealed a close relationship between soil moisture and organic carbon pools (measured as SOC, WEOC, and TOAC) among examined soil microtopographic features and across soil depth (p<0.01, Table S4). All these soil properties significantly correlated with cumulative CO2 and CH4 production (p<0.05), suggesting the important role of initial soil geochemical properties in controlling carbon degradation.

Although various levels of carbon mineralization were measured as CO2 and CH4 production during incubations, changes in WEOC and TOAC were consistent among treatments with distinct patterns. WEOC represents 0.3 % to 2.6 % of total SOC among all test soils, and this ratio remained constant before and after anoxic incubations (Fig. S3). On the other hand, TOAC showed much more dynamic changes among different soils and different incubation temperatures. TOAC generally increased in soils from the organic layer, transition zone, and permafrost. In contrast, TOAC drastically decreased by up to 90 % in mineral soils. These results indicate that WEOC was in a steady state among examined soils, while TOAC varied substantially due to microbial mineralization processes, supporting the model assumption of lumped fermentation (the conversion of WEOC to TOAC) as the rate-limiting step.

Both CO2 and CH4 production rates responded strongly to rising incubation temperature (p=0.02 and p=0.04, respectively; Fig. S4, Tables S5 and S6). The estimated Q10 values of CO2 production showed a relatively narrow range, while methanogenesis had much larger variation in estimated Q10 values ranging from 1.6 to 48.1. Using Q10 values to simulate the temperature dependence of processes might work for CO2 production but could generate significant errors in predicting CH4 production.

Figure 2Comparison between modeled and observed production of CO2 (a) and CH4 (b). Averaged measurements of triplicate microcosms at each time point from each incubation temperature were calculated as observed values.


3.2 Modeled CO2 and CH4 production using observed parameters

The model performed well in simulating CO2 and CH4 dynamics across a range of moisture and SOC gradients and among different soil types (Figs. S5 and S6). Variations in gas production among different conditions, including microtopographic features, soil layers, and different incubation temperatures, were well captured (Fig. S7). The comparisons between modeled and observed CO2 and CH4 production are shown in Fig. 2. The model slightly underestimates CO2 production towards the end of the incubations but still maintains a good agreement between modeled and observed CO2 production (R2=0.89). The underestimation of CO2 production is likely due to substrate limitations caused by the initial distribution of different carbon pools. Model-predicted CH4 production also showed good agreement with observations (R2=0.79). More variation between modeled and observed CH4 production suggests a systematic pattern in the model parameterization of methanogenesis: the model underestimates CH4 production at 4 and 8 C, and overestimates CH4 production at −2C.

Figure 3Perturbations of initial soil geochemical conditions differentially affected model predictions (including CH4, CO2, Fe(II), TOAC, WEOC, and pH) during anaerobic carbon decomposition. For example, when the initial pH decreased by 8 % and 17 %, CH4 production decreased by 40 % and 80 %, respectively. Normalized changes in model output were calculated as the ratio of changes caused by perturbation simulations (differences between perturbation and reference runs) to reference simulation output after 60 days of anaerobic decomposition at 8 C. To test model sensitivity in response to initial pH, the reference run started with pH 6, and up to 1 pH unit changes were applied in perturbation simulations to represent a realistic pH range for soils. Reference simulations were based on soils with 30 % SOC (water content of 2 g g−1 dwt, and pH = 5).


To assess the model sensitivity to initial model inputs, we compared model predictions in response to varying initial model inputs via perturbation simulations. First, we examined the influence of the partitioning of different carbon pools. Significant changes in model predictions of CO2 and CH4 were observed in response to perturbations of initial input of SOC and WEOC but not TOAC (Fig. 3). SOC determines the size of different carbon pools in the model, and it further influences the predictions of WEOC, TOAC, CO2, and CH4. For example, predicted CO2 and CH4 production increased by about 200 % when +200 % changes were applied to initial SOC input. This trend is consistent with correlation analysis of incubation results, described above (Table S4). Perturbations in initial WEOC strongly altered the predictions of TOAC and CO2, consistent with the model assumption of the conversion of WEOC to TOAC (fermentation process) as the rate-limiting step. The model also predicted increases in CH4 and Fe(II) accumulation in response to lower WEOC. Lower WEOC significantly reduced organic acid accumulation and thus increased system pH and accelerated rates of both methanogenesis and iron reduction. The starting level of TOAC showed minimal influence on model predictions of CO2 and CH4, suggesting other factors rather than substrate availability were limiting carbon mineralization. The initial sizes of SOM1 and SOM2 pools showed very slight changes in model predictions of WEOC and CO2, and minimal influence on CH4 prediction, further supporting the assumption that downstream fermentation is the rate-limiting step in the model. Additional soil geochemical factors, including soil moisture, Fe(II), and pH, also significantly influence model output. In particular, initial soil pH showed a dramatic effect on predicted CO2 and CH4 production. With initial soil pH increasing from 5 (reference simulation) to 6, the model predicted 160 % and 308 % increases in CO2 and CH4 production, respectively. Perturbations in initial soil pH had the strongest effect on the prediction of CH4 by assigning different values in fpH that were directly proportional to the methanogenesis rates. The above results of perturbation simulations demonstrated the high sensitivity of this model in response to varying soil geochemical properties.

Figure 4Simulated changes in model predictions (including CH4, CO2, Fe(II), TOAC, WEOC, and pH) during anaerobic carbon decomposition in response to perturbations of (a) fermentation rate and (b) fermentation stoichiometry (acetate : CO2 is 1:1 for reference simulation). Normalized changes in model output were calculated as the ratio of perturbation simulation output to reference simulation output after 60 days of anaerobic decomposition at 8 C. Reference simulations were based on soils with 30 % SOC (water content of 2 g g−1 dwt, and pH = 6).


3.3 Model sensitivity to parameterization uncertainties

To further validate the model, we performed additional sensitivity analysis to justify model assumptions and estimate the uncertainties generated from model parameterizations. One major assumption of this modeling framework is to lump multiple fermentation processes into one reaction stoichiometry, controlled by one reaction rate constant. It is critical to evaluate how this simplified structure influences model performance and contributes to model output uncertainties. The model parameter sensitivity analysis indicated the TOAC pool was most sensitive to changes in the fermentation rate (Rfer) and reaction stoichiometry (Fig. 4). Downstream reactions were less affected by the uncertainties of the two tested parameters. These results supported our assumption of lumped fermentation with fixed stoichiometry, indicating the robustness of the model structure presented here.

Figure 5Comparison of simulated and observed temperature responses for the production of CO2 (a) and CH4 (b). Results were all normalized to CO2 or CH4 production rates at 8 C for direct comparison. Observations at −2 and 4 C were plotted in black dots and the median values were marked in red. The shaded area represents output uncertainties generated from rate estimations within 60±5 days. Reference simulations were based on soils with 30 % SOC (water content of 2 g g−1 dwt, and pH = 5).


The selection of temperature response functions represents one of the major sources of model uncertainties. A sensitivity analysis was performed by comparing four different temperature response functions (Appendix B). In our simulations, the quadratic temperature response function proposed by Ratkowsky et al. predicted much higher CO2 and CH4 production rates at higher temperature, and the lowest rates of both CO2 and CH4 at temperatures below 0 C, giving the highest temperature response among tested response functions (Fig. 5). In contrast, the Arrhenius equation predicted much lower temperature response for both CO2 and CH4. Empirical functions used in the CLM-CN and CENTURY models gave similar temperature responses for both CO2 and CH4. Variations in low-temperature CO2 production are well constrained by established temperature response functions, while CH4 production at −2C showed a much wider range of temperature response, and the median value is best simulated using the Ratkowsky function. This sensitivity analysis is consistent with model output of CO2 and CH4 production, where CO2 is well constrained by the model, but CH4 is significantly overestimated at −2C using the CLM-CN temperature response function. A unified temperature response function for all reactions under different biotic or abiotic constraints substantially contributes to the disagreement between model output and observations.

Figure 6Simulated changes in model predictions (including CH4, CO2, Fe(II), TOAC, WEOC, pH, and fpH) during anaerobic carbon decomposition in response to perturbations of (a) pH buffering capacity and (b) pH response function. Normalized changes in model output were calculated as the ratio of perturbation simulation output to reference simulation output after 60 days of anaerobic decomposition at 8 C. Reference simulations were based on soils with 30 % SOC (water content of 2 g g−1 dwt, and pH = 5).


Redox reactions contribute to proton production or consumption, and the resulting pH alters the value of the pH response function (fpH) that directly controls reaction kinetic functions, creating a feedback loop. pH buffering capacity (BC) provided by SOM with proton-binding sites and fpH represent two major sources of uncertainties in this feedback loop. Thus, we performed perturbation simulations to characterize the sensitivity of model output to variations in BC and fpH (Fig. 6). Higher BC stabilized system pH during prolonged incubations, while lower BC permitted a pH increase by up to 0.71 pH unit compared to the reference simulation. This 14 % pH increase led to a 123 % increase in fpH, accelerating both methanogenesis and Fe(III) reduction rates substantially. Perturbations on the pH response function were directly reflected in the slopes of pH response curves (Fig. S8). We found up to 372 % change in the value of fpH during a 60-day simulation, as a steeper increase in fpH accelerated both methanogenesis and iron reduction (Reactions AR2–AR5), which contributed to pH rise that further accelerated fpH increase. Correspondingly, both CH4 and Fe(II) increased by more than 100 % after the simulation. While BC is an important factor controlling both redox reactions and pH fluctuations, a unified fpH for all reactions may impose significant variations in model output.

Figure 7Temperature response of CH4 and Fe(II) production rates at varying soil pH buffering capacities (BCs). Varying BCs with respect to the reference simulation (BC = 1) creates strong feedback to rates of methanogenesis and iron reduction. Reference simulations were based on soils with 30 % SOC (water content of 2 g g−1 dwt, and pH = 5).


Figure 8Cluster analysis of soil geochemical properties related to CO2 and CH4 production using Ward's linkage method. (a) Cluster analysis of measured soil geochemical characteristics and observed CO2 and CH4 production (n=42); (b) cluster analysis of modeled results (n=42). Model-simulated CO2, CH4, and Fe(II) production, and final pH values are labeled as M_CO2, M_CH4, M_Fe, and M_pH, respectively. Biomasses of methanogens and iron reducers were tracked in the model and labeled M_Meb and M_Feb, respectively.


BC is an intrinsic soil property simulated with a simplified linear relationship to soil SOM. However, it generates a strong nonlinear response in the simulations of methanogenesis and Fe(III) reduction (Fig. 7a). Simulations with varying soil BC revealed dynamic pH change at lower BC (Figs. 8 and 9, with BC = 1 as the reference simulation) and stabilized pH at higher BC. At constant temperature, rates of both methanogenesis and Fe(III) reduction increased significantly at lower BC due to pH control. At lower BC, when pH change is not well buffered, higher pH accelerated CH4 and Fe(II) production rates (Fig. 7), giving much higher apparent temperature responses, while at higher BC with stabilized pH in the system, apparent temperature responses of these redox processes were significantly lower than the reference simulation (BC = 1). Variations in pH buffering capacity generated large variations in apparent temperature responses of methanogenesis and Fe(III) reduction due to this pH feedback loop.

4 Discussion

4.1 Synthesized soil geochemistry and model validation

Soil geochemical characteristics represent important abiotic controls on anaerobic carbon decomposition and subsequent CO2 and CH4 production. SOC content, soil pH, water table position, C:N ratio, and landscape position were all suggested to contribute to the variability in anaerobic CO2 and CH4 production (Lee et al., 2012; Schädel et al., 2014; Treat et al., 2015). We synthesized incubation data for gelisol soils from different pedons and soil moisture regimes representing heterogeneity across the Barrow Environmental Observatory (BEO). This coordinated dataset allowed us to focus on individual factors and their roles in relation to anaerobic CO2 and CH4 production.

Carbon released as CO2 and CH4 during anoxic incubations decreased with depth. Permafrost was associated with low levels of CO2 production and very low CH4 production, consistent with a previous synthesis (Treat et al., 2015). Nevertheless, permafrost TOAC, WEOC, and SOC concentrations were all comparable to organic soils, suggesting high substrate availability but low microbial activity. This trend is consistent with previous studies (Walz et al., 2017; Treat et al., 2015), where highest microbial abundance and diversity were observed in surface soil, and permafrost contained low microbial abundance (Treat et al., 2014; Waldrop et al., 2010). Among surface soils, higher moisture in low-centered polygon soils significantly promoted CO2 and CH4 production and the accumulation of fermentation products (measured as TOAC), emphasizing the importance of soil SOC content and moisture as strong environmental drivers for carbon decomposition. Given the bias in correlation analysis created by the skewed distribution of CO2 and CH4 production in our dataset, additional cluster analysis was performed based on data similarity rather than correlations. High similarity of soil attributes (depth, moisture, pH, C:N ratio, SOC, TOAC) with CH4 production (Fig. 8a) was found, suggesting methanogenesis is potentially controlled by a set of soil geochemical characteristics in the local microenvironment.

These synthesized observations support the major assumptions of our model development: (1) the coupled hydrolysis and fermentation processes converting macromolecular SOM into low-molecular-weight organic acids are the rate-limiting step; and (2) different rates of CO2 and CH4 production from different soil layers can be attributed to variations in microbial activity manifested as differences in initial microbial biomass or growth rates. Additional observations of substantial Fe(III) reduction and associated pH increases during anaerobic decomposition (Fig. S9) confirmed the need to simulate pH variations associated with redox reactions and corresponding microbial responses. This anaerobic carbon decomposition framework adequately modulated the involved biotic and abiotic interactions by splitting the carbon flow to different redox reactions and simulating pH buffering capacity to mediate associated changes in acidity or alkalinity.

The model presented here identified fermentation, acetoclastic methanogenesis, and acetotrophic iron reduction as key mechanisms for anaerobic CO2 and CH4 production (Vaughn et al., 2016; Lipson et al., 2010). Although denitrification, ammonification, and sulfate reduction are all thermodynamically more favorable, low nitrate and sulfate concentrations in BEO soils limit flux through these pathways (Newman et al., 2015). We performed another cluster analysis on the model output (Fig. 8b), where we not only simulated fermentation, methanogenesis and iron reduction rates, and associated pH changes but also tracked the biomass of methanogens (M_Meb) and iron reducers (M_Feb). A dendrogram depicting data similarity showed four distinct clusters consisting of WEOC, CO2 (CO2 prediction), ferrous (Fe(II) prediction), and CH4 (CH4 prediction) that closely associated with soil geochemical properties and incubation temperature. This result is similar to the cluster analysis of synthesized data, demonstrating that the proposed model structure captured major relationships between carbon mineralization and soil geochemical attributes. Predicted CH4 production is strongly influenced by incubation temperature, soil pH, soil moisture, and depth which determine the size of the methanogen population. This model prediction is consistent with previous studies on the vertical distribution of methanogen population (Waldrop et al., 2010). Environmental factors, such as labile organic matter, water table depth, and soil redox status, soil alkalinity, and salinity (Wachinger et al., 2000; Rivkina et al., 2007; Høj et al., 2006; Yang et al., 2017) are all likely to contribute to the variabilities in the distribution and abundance of methanogens and subsequent methane production.

4.2 Temperature and pH response of anaerobic carbon decomposition

Rising temperature promotes anaerobic carbon decomposition, resulting in increased rates of anaerobic CO2 and CH4 production (Treat et al., 2014; Lupascu et al., 2012). It is widely recognized that methanogenesis is more sensitive to temperature than respiration (Yvon-Durocher et al., 2012, 2014), and it is usually associated with large variations. Segers estimated the Q10 value of methanogenesis ranged from 1.5 to 28 among 1043 incubation experiments using wetland soils (Segers, 1998). Our data synthesis revealed higher temperature sensitivity than other reported values. High estimated temperature sensitivity across the freezing point of water has previously been documented (Waldrop et al., 2010) and further attributed to limited water availability for microbial activities at subzero temperature (Tilston et al., 2010). Ratkowsky et al. proposed a quadratic relationship for the temperature dependence of microbial growth rates that modeled low-temperature growth better than the Arrhenius law (Ratkowsky et al., 1982). Our simulations suggest better prediction of methanogenesis with this temperature response function, possibly due to a more suitable representation of growth limitation of methanogens at subzero temperature. Methanogenesis rates are also influenced by the availability of alternative electron acceptors and carbon source. Processes contributing to the accumulation or consumption of carbon substrates and competing electron acceptors may respond differently to temperature change, which could further complicate the temperature sensitivity of methanogenesis. Current modeling approaches heavily depend upon empirical temperature response functions, which may be associated with large uncertainties due to variations in the selection of data and curve-fitting methods. Extrapolation of carbon decomposition rates, particularly methanogenesis rates, into a future warmer climate remains uncertain. More accurate simulations will require additional information on geochemical properties that contribute to the variations of methanogens distribution and methanogenesis activity.

pH values impose fundamental physiological restrictions on microbial activities. Soil pH ranges from acidic to circumneutral (pH 4–7.5) in northern Alaska and varies substantially through the soil profile and along the microtopographic gradient. Accumulation of organic acids in anoxic soils leads to pH decline (Jones et al., 2003), while consumption of organic acids by methanogenesis and iron reduction increases the alkalinity of the system via the production of HCO3- and OH (Drake et al., 2015; Roy Chowdhury et al., 2015; Howell et al., 1998). The interplay of these processes leads to strong nonlinear pH feedbacks in the system, and previous studies have observed up to 1–2 pH unit changes during short-term anoxic incubations (Xu et al., 2015; Drake et al., 2015; Roy Chowdhury et al., 2015). These relationships between pH and organic carbon decomposition can vary in sign and magnitude. Our model simulations with mechanistic pH evolution indicate that constant pH assumed in previous models may cause significant errors in simulating long-term anaerobic CO2 and CH4 production. The intrinsic soil pH buffering capacity plays a large role in stabilizing soil pH and may be heterogeneous depending upon solution acidity or alkalinity, cation exchange capacity and residual acidity or mineral dissolution. These properties derive from SOM characteristics, moisture, mineral content, and additional geochemical properties, leading to complex correlations between soil pH and SOC decomposition rate that require future investigation.

4.3 Fast-decomposing carbon pool

Substrate availability is a primary determinant of potential CO2 and CH4 production (Lee et al., 2012; Schuur et al., 2015; Tarnocai et al., 2009). Total SOC is composed of heterogeneous C pools characterized by different turnover times. Carbon release during short-term incubation originates from the C pool with relatively rapid turnover. The size and turnover time of this quickly metabolized carbon pool are usually estimated by two-pool or three-pool conceptual models with a maximum likelihood solution using time series of CO2 data (Schädel et al., 2013). A previous study on Siberian permafrost soils using a two-pool model estimated a turnover time of 0.26 years for the fastest-responding pool (Knoblauch et al., 2013). A three-pool model was applied using more extensive incubation datasets collected from 23 high-latitude ecosystems, yielding an estimate of a 0.35-year mean turnover time for the fastest-responding carbon pool (Schädel et al., 2014).

In our synthesis study, we directly quantified WEOC and assumed it represented the fast-decomposing labile carbon pool. The size of the labile carbon pool is constant during anaerobic decomposition, while total CO2 and CH4 release represent up to 194 % of the labile carbon pool, indicating continuous replenishment of the labile carbon pool from non-labile carbon pools within the hierarchy. The replenishment of the labile carbon pool can be attributed mostly to decomposition of SOM1 and SOM2 pools with faster turnover (Koven et al., 2013). Overall, we estimated the fast-decomposed carbon pool is approximately 2 %–4 % of total SOC, similar to previous estimates. The turnover time calculated from the fermentation rate was comparable to estimates of the turnover time of the fastest-responding carbon pool in previous studies (Fig. 9), suggesting these quantifications and parameterization in the anaerobic carbon decomposition framework apply broadly.

Figure 9Model-estimated turnover rates of the fastest-decomposing carbon pool. Organic, mineral, and permafrost labels represent estimations from our model simulations (rates estimated at 4 C). Schadel data represent turnover rates estimated via a three-pool model from pooled anaerobic incubations with normalized incubation temperature of 5 C (tags 1, 2, and 3 represent pool estimation from different soil types: 1. organic, 2. mineral <1 m, 3. mineral >1 m). Knoblauch data are rate estimates (at 4 C) made via a two-pool model (Schädel et al., 2014; Knoblauch et al., 2013). Open symbols represent the average values, and the vertical lines represent the estimated range.


4.4 Key features of the anaerobic model framework and future considerations

Here, we present an anaerobic carbon decomposition framework by combining three well-known modeling approaches developed in different disciplines. A pool-based model to represent upstream carbon transformations and replenishment of a DOC pool, a thermodynamically based model to calculate rate kinetics and biomass growth for methanogenesis and Fe(III) reduction, and a humic ion-binding model for aqueous-phase speciation and pH calculation are implemented into the open-source geochemical model PHREEQC (Charlton and Parkhurst, 2011). The model framework presented here has several unique features. First, this model is built upon a thermodynamically based approach, which allows consistent parameterization of individual reactions along the redox ladder. Such a model structure is particularly useful in circumstances when function-specific microbial growth is difficult to quantify and parameterize. Second, calculations of free energy changes of redox couples are used to modulate redox reaction hierarchy. Considering the difficulty in obtaining growth-associated parameters for every functional group, a thermodynamically based approach significantly decreases the number of parameters that are difficult to measure. In addition, proton production and consumption during redox reactions are incorporated into a dynamic pH calculation, allowing various simulations on aqueous solubility and reactivity of different elements. The anaerobic carbon decomposition framework presented here holds a significant advantage over traditional models in simulating carbon decomposition process within a wide range of environmental settings.

In permafrost-affected regions, studies consistently identify iron reduction, denitrification, and sulfate reduction (Lipson et al., 2010, 2013; Ernakovich et al., 2017; Hansen et al., 2007) as alternative anaerobic pathways, which are recognized as energetically more favorable processes than methanogenesis. Fe reduction makes a significant contribution to total respiration (Roy Chowdhury et al., 2015; Herndon et al., 2015), and adding Fe reduction simulations to a baseline model (without Fe reduction or dynamic pH calculations) caused faster decreases in TOAC and WEOC pools and increased CO2 production as expected (Fig. S10). More indirect feedbacks were revealed when dynamic pH calculation was enabled. With dynamic pH simulation during anaerobic decomposition, the model revealed strong pH dynamics that are counterbalanced by Fe reduction. By including both Fe reduction and dynamic pH calculations, the model accurately reproduced the initial pH drop and subsequent pH rise during incubations, which were commonly observed in permafrost-affected soils (Roy Chowdhury et al., 2015; Herndon et al., 2015). The new model framework presented here provides a basis for a deeper understanding of carbon decomposition under oxygen-limited conditions where the importance of accounting for alternative election acceptors and pH feedbacks becomes more pronounced. Future fine-scale experiments on carbon decomposition using alternative electron acceptors would be beneficial for more comprehensive parameterization of this model framework. Additional observations on temperature and pH sensitivity of specific redox reactions would also be quite useful in reducing large uncertainties generated by the current representation of temperature and pH responses. Application of such a modeling framework at the field scale requires close coupling with hydrology models to facilitate estimations of aqueous-phase concentrations. Additional assumptions on vertical mixing and gas diffusion in the soil column should also be considered.

5 Conclusions

Microbial processes are the driving forces for biogeochemical cycling of soil carbon and are subjected to environmental constraints beyond temperature and organic substrate availability. The present study incorporated microbial redox reactions and mechanistic pH evolution to simulate anaerobic carbon decomposition in Arctic soils with depth and across soil moisture gradients. Our data synthesis and modeling results quantify direct effects of temperature on anaerobic carbon decomposition, as well as indirect effects of soil geochemistry that cause strong redox reaction–pH feedback. We identified substantial pH feedbacks on the predicted CO2 and CH4 production. The anaerobic carbon decomposition framework presented in this study provided the essential model structure to incorporate redox reactions of alternative electron acceptors for accurate simulation of CO2 and CH4 production. Soil geochemistry imposes critical constraints on SOM decomposition and further regulates permafrost carbon feedback in response to changing climate.

Code and data availability

PHREEQC (version 3) is publicly available at (last access: January 2019).

The model is archived at (Zheng et al., 2018c), with a detailed description of model implementation, input files, and various sensitivity analyses described in this paper.

Datasets used in this work can be found at (Herndon et al., 2017), (Zheng and Graham, 2018), and (Zheng et al., 2017), and a synthesis of the incubation data is available at (Zheng et al., 2018a).

Appendix A: Anaerobic carbon decomposition model

This section lists reactions used in the anaerobic carbon decomposition model. Under anaerobic conditions, dissolved organic carbon is converted to low-molecular-weight organic acids via fermentation. One simplified fermentation reaction is used to represent this lumped fermentation process, where one-third of the fermented organic carbon is converted to CO2 (Tang et al., 2016; Xu et al., 2015):

(AR1) C 6 H 12 O 6 + 4 H 2 O 2 CH 3 COO - + 2 HCO 3 - + 4 H + + 4 H 2 .

This fermentation reaction generates protons and decreases pH in the system. Fermentation products acetate and H2 are further consumed via methanogenesis and iron reduction. The growth equations of methanogenesis and iron reduction were derived for each group using a thermodynamically based approach, in which biomass synthesis is included in paired electron donor and electron acceptor half reactions. A general molecular formula (C5H7O2N) is used for microbial biomass, and the growth equations are written as (Istok et al., 2010)





In addition, Fe(III) concentration was calculated based on the dissolution of representative amorphous ferric hydroxides (Reaction AR6), with a solubility constant Ks0=103.96. This process consumes many protons and contributes to pH increases.

(AR6) Fe ( OH ) 3 ( s ) + 3 H + Fe 3 + + 3 H 2 O

A complete set of rate constants used in this model can be found in Table S1.

Appendix B: Temperature and pH response functions

We used the CLM-CN temperature response function (Eq. B1) in our simulations (Thornton and Rosenbloom, 2005). Additional tested temperature response functions included Eq. (B2), which is used by the CENTURY model (Grosso et al., 2005), the Arrhenius equation (Eq. B3) used in ecosys (Grant, 1998), and the quadratic Eq. (B4) (Ratkowsky et al., 1982). Tref is set at 25 C, Ea is the activation energy (J mol−1), and R is the universal gas constant (J K−1 mol−1). Tm used in Ratkowsky's model represents a conceptual temperature of no metabolic significance and is set at −8C in this study.

(B1) ln f T = 308.56 × 1 71.02 - 1 T - 227.13

(B2) f T = 0.56 + 0.465 arctan [ 0.097 T - 15.7 ]

(B3) f T = e - E a R 1 T - 1 T ref

(B4) f T = T - T m T ref - T m 2

The discontinuous bell-shaped pH response function from the DLEM model was used here (Eq. B5; Tian et al., 2010):

(B5) f pH = 1.02 1.02 + 10 6 exp ( - 2.5 pH ) ( 0 < pH < 7 )

(B6) f pH = 1.02 1.02 + 10 6 exp ( - 2.5 14 -pH ) ( 7 < pH < 14 ) .

The supplement related to this article is available online at:

Author contributions

DG, SDW, and BG conceived and organized the research study; PET, SLP, DEG, and JZ built the conceptual model framework; JZ preformed all model simulations; JZ and DEG drafted the manuscript. All authors contributed revisions to the manuscript and have given approval to the final version of the manuscript.

Competing interests

The authors declare that they have no conflict of interest.


We appreciate comments and suggestions on earlier versions of this paper offered by Ethan Coon and manuscript reviewers. The Next-Generation Ecosystem Experiments in the Arctic (NGEE Arctic) project is supported by the Biological and Environmental Research program in the US Department of Energy (DOE) Office of Science. Oak Ridge National Laboratory is managed by UT-Battelle, LLC, for the DOE under contract no. DE-AC05-00OR22725.

Edited by: Anja Rammig
Reviewed by: Ali Ebrahimi and one anonymous referee


Atalay, Y. B., Carbonaro, R. F., and Di Toro, D. M.: Distribution of proton dissociation constants for model humic and fulvic acid molecules, Environ. Sci. Technol., 43, 3626–3631,, 2009. 

Bethke, C. M., Sanford, R. A., Kirk, M. F., Jin, Q., and Flynn, T. M.: The thermodynamic ladder in geomicrobiology, Am. J. Sci., 311, 183–210,, 2011. 

Blumer-Schuette, S. E., Brown, S. D., Sander, K. B., Bayer, E. A., Kataeva, I., Zurawski, J. V., Conway, J. M., Adams, M. W. W., and Kelly, R. M.: Thermophilic lignocellulose deconstruction, FEMS Microbiol. Rev., 38, 393–448,, 2014. 

Bore, E. K., Apostel, C., Halicki, S., Kuzyakov, Y., and Dippold, M. A.: Microbial metabolism in soil at subzero temperatures: adaptation mechanisms revealed by position-specific 13C labeling, Front. Microbiol., 8, 946,, 2017. 

Boye, K., Noël, V., Tfaily, M. M., Bone, S. E., Williams, K. H., Bargar, John R., and Fendorf, S.: Thermodynamically controlled preservation of organic carbon in floodplains, Nature Geosci., 10, 415,, 2017. 

Charlton, S. R. and Parkhurst, D. L.: Modules based on the geochemical model PHREEQC for use in scripting and programming languages, Comput. Geosci., 37, 1653–1663,, 2011. 

Davidson, E. A. and Janssens, I. A.: Temperature sensitivity of soil carbon decomposition and feedbacks to climate change, Nature, 440, 165–173,, 2006. 

Drake, T. W., Wickland, K. P., Spencer, R. G. M., McKnight, D. M., and Striegl, R. G.: Ancient low–molecular-weight organic acids in permafrost fuel rapid carbon dioxide production upon thaw, P. Natl. Acad. Sci. USA, 112, 13946–13951,, 2015. 

Ernakovich, J. G., Lynch, L. M., Brewer, P. E., Calderon, F. J., and Wallenstein, M. D.: Redox and temperature-sensitive changes in microbial communities and soil chemistry dictate greenhouse gas loss from thawed permafrost, Biogeochemistry, 134, 183–200,, 2017. 

Glissmann, K. and Conrad, R.: Saccharolytic activity and its role as a limiting step in methane formation during the anaerobic degradation of rice straw in rice paddy soil, Biol. Fertility Soils, 35, 62–67,, 2002. 

Grant, R. F.: Simulation of methanogenesis in the mathematical model ecosys, Soil Biol. Biochem., 30, 883–896,, 1998. 

Grosso, S. J. D., Parton, W. J., Mosier, A. R., Holland, E. A., Pendall, E., Schimel, D. S., and Ojima, D. S.: Modeling Soil CO2 Emissions from Ecosystems, Biogeochemistry, 73, 71–91,, 2005. 

Hansen, A. A., Herbert, R. A., Mikkelsen, K., Jensen, L. L., Kristoffersen, T., Tiedje, J. M., Lomstein, B. A., and Finster, K. W.: Viability, diversity and composition of the bacterial community in a high Arctic permafrost soil from Spitsbergen, Northern Norway, Environ. Microbiol., 9, 2870–2884,, 2007. 

Herndon, E., Yang, Z., and Gu, B.: Soil Organic Carbon Degradation during Incubation, Barrow, Alaska, 2012, Next Generation Ecosystem Experiments Arctic Data Collection,, 2017. 

Herndon, E. M., Yang, Z., Bargar, J., Janot, N., Regier, T. Z., Graham, D. E., Wullschleger, S. D., Gu, B., and Liang, L.: Geochemical drivers of organic matter decomposition in arctic tundra soils, Biogeochemistry, 126, 397–414,, 2015. 

Høj, L., Rusten, M., Haugen, L. E., Olsen, R. A., and Torsvik, V. L.: Effects of water regime on archaeal community composition in Arctic soils, Environ. Microbiol., 8, 984–996,, 2006. 

Howell, J., Donahoe, R., Roden, E., and Ferris, F.: Effects of microbial iron oxide reduction on pH and alkalinity in anaerobic bicarbonate-buffered media: implications for metal mobility, Goldschmidt Conference, Tolouse, 657, 1998. 

Istok, J. D., Park, M., Michalsen, M., Spain, A. M., Krumholz, L. R., Liu, C., McKinley, J., Long, P., Roden, E., Peacock, A. D., and Baldwin, B.: A thermodynamically-based model for predicting microbial growth and community composition coupled to system geochemistry: Application to uranium bioreduction, J. Contam. Hydrol., 112, 1–14,, 2010. 

Jin, Q. and Bethke, C. M.: A New rate law describing microbial respiration, Appl. Environ. Microbiol., 69, 2340–2348,, 2003. 

Jin, Q. and Bethke, C. M.: The thermodynamics and kinetics of microbial metabolism, Am. J. Sci., 307, 643–677,, 2007. 

Jin, Q. and Kirk, M. F.: pH as a primary control in environmental microbiology: 1. Thermodynamic perspective, Front. Environ. Sci., 6,, 2018. 

Jones, D. L., Dennis, P. G., Owen, A. G., and van Hees, P. A. W.: Organic acid behavior in soils – misconceptions and knowledge gaps, Plant Soil, 248, 31–41,, 2003. 

Karhu, K., Auffret, M. D., Dungait, J. A. J., Hopkins, D. W., Prosser, J. I., Singh, B. K., Subke, J.-A., Wookey, P. A., Agren, G. I., Sebastia, M.-T., Gouriveau, F., Bergkvist, G., Meir, P., Nottingham, A. T., Salinas, N., and Hartley, I. P.: Temperature sensitivity of soil respiration rates enhanced by microbial community response, Nature, 513, 81–84,, 2014. 

Knoblauch, C., Beer, C., Sosnin, A., Wagner, D., and Pfeiffer, E.-M.: Predicting long-term carbon mineralization and trace gas production from thawing permafrost of Northeast Siberia, Global Change Biol., 19, 1160–1172,, 2013. 

Kotsyurbenko, O. R.: Trophic interactions in the methanogenic microbial community of low-temperature terrestrial ecosystems, FEMS Microbiol. Ecol., 53, 3–13,, 2005. 

Koven, C. D., Riley, W. J., Subin, Z. M., Tang, J. Y., Torn, M. S., Collins, W. D., Bonan, G. B., Lawrence, D. M., and Swenson, S. C.: The effect of vertically resolved soil biogeochemistry and alternate soil C and N models on C dynamics of CLM4, Biogeosciences, 10, 7109–7131,, 2013. 

Koven, C. D., Lawrence, D. M., and Riley, W. J.: Permafrost carbon-climate feedback is sensitive to deep soil carbon decomposability but not deep soil nitrogen dynamics, P. Natl. Acad. Sci. USA, 112, 3752–3757,, 2015. 

Lawrence, D. M., Koven, C. D., C., S. S., Riley, W. J., and Slater, A. G.: Permafrost thaw and resulting soil moisture changes regulate projected high-latitude CO2 and CH4 emissions, Environ. Res. Lett., 10, 094011,, 2015. 

Lee, H., Schuur, E. A. G., Inglett, K. S., Lavoie, M., and Chanton, J. P.: The rate of permafrost carbon release under aerobic and anaerobic conditions and its potential effects on climate, Global Change Biol., 18, 515–527,, 2012. 

Lipson, D. A., Jha, M., Raab, T. K., and Oechel, W. C.: Reduction of iron (III) and humic substances plays a major role in anaerobic respiration in an Arctic peat soil, J. Geophys. Res.-Biogeosci., 115, G00I06,, 2010. 

Lipson, D. A., Raab, T. K., Goria, D., and Zlamal, J.: The contribution of Fe(III) and humic acid reduction to ecosystem respiration in drained thaw lake basins of the Arctic Coastal Plain, Global Biogeochem. Cy., 27, 399–409,, 2013. 

Lupascu, M., Wadham, J. L., Hornibrook, E. R. C., and Pancost, R. D.: Temperature sensitivity of methane production in the permafrost active layer at Stordalen, Sweden: A comparison with non-permafrost northern wetlands, Arct. Antarct. Alp. Res., 44, 469–482,, 2012. 

MacDougall, A. H., Avis, C. A., and Weaver, A. J.: Significant contribution to climate warming from the permafrost carbon feedback, Nature Geosci., 5, 719–721,, 2012. 

McGuire, A. D., Christensen, T. R., Hayes, D., Heroult, A., Euskirchen, E., Kimball, J. S., Koven, C., Lafleur, P., Miller, P. A., Oechel, W., Peylin, P., Williams, M., and Yi, Y.: An assessment of the carbon balance of Arctic tundra: comparisons among observations, process models, and atmospheric inversions, Biogeosciences, 9, 3185–3204,, 2012. 

Meng, L., Hess, P. G. M., Mahowald, N. M., Yavitt, J. B., Riley, W. J., Subin, Z. M., Lawrence, D. M., Swenson, S. C., Jauhiainen, J., and Fuka, D. R.: Sensitivity of wetland methane emissions to model assumptions: application and model testing against site observations, Biogeosciences, 9, 2793–2819,, 2012. 

Newman, B. D., Throckmorton, H. M., Graham, D. E., Gu, B., Hubbard, S. S., Liang, L., Wu, Y., Heikoop, J. M., Herndon, E. M., Phelps, T. J., Wilson, C. J., and Wullschleger, S. D.: Microtopographic and depth controls on active layer chemistry in Arctic polygonal ground, Geophys. Res. Lett., 42, 1808–1817,, 2015. 

Ratkowsky, D. A., Olley, J., McMeekin, T. A., and Ball, A.: Relationship between temperature and growth rate of bacterial cultures, J. Bacteriol., 149, 1–5, 1982. 

Riley, W. J., Subin, Z. M., Lawrence, D. M., Swenson, S. C., Torn, M. S., Meng, L., Mahowald, N. M., and Hess, P.: Barriers to predicting changes in global terrestrial methane fluxes: analyses using CLM4Me, a methane biogeochemistry model integrated in CESM, Biogeosciences, 8, 1925–1953,, 2011. 

Rivkina, E., Shcherbakova, V., Laurinavichius, K., Petrovskaya, L., Krivushin, K., Kraev, G., Pecheritsina, S., and Gilichinsky, D.: Biogeochemistry of methane and methanogenic archaea in permafrost, FEMS Microbiol. Ecol., 61, 1–15,, 2007. 

Roy Chowdhury, T., Herndon, E. M., Phelps, T. J., Elias, D. A., Gu, B., Liang, L., Wullschleger, S. D., and Graham, D. E.: Stoichiometry and temperature sensitivity of methanogenesis and CO2 production from saturated polygonal tundra in Barrow, Alaska, Global Change Biol., 21, 722–737,, 2015. 

Schädel, C., Luo, Y., David Evans, R., Fei, S., and Schaeffer, S. M.: Separating soil CO2 efflux into C-pool-specific decay rates via inverse analysis of soil incubation data, Oecologia, 171, 721–732,, 2013. 

Schädel, C., Schuur, E. A. G., Bracho, R., Elberling, B., Knoblauch, C., Lee, H., Luo, Y., Shaver, G. R., and Turetsky, M. R.: Circumpolar assessment of permafrost C quality and its vulnerability over time using long-term incubation data, Global Change Biol., 20, 641–652,, 2014. 

Schaefer, K., Lantuit, H., Romanovsky, V. E., Schuur, E. A. G., and Witt, R.: The impact of the permafrost carbon feedback on global climate, Environ. Res. Lett., 9, 085003,, 2014. 

Schuur, E. A. G., Vogel, J. G., Crummer, K. G., Lee, H., Sickman, J. O., and Osterkamp, T. E.: The effect of permafrost thaw on old carbon release and net carbon exchange from tundra, Nature, 459, 556–559,, 2009. 

Schuur, E. A. G., McGuire, A. D., Schädel, C., Grosse, G., Harden, J. W., Hayes, D. J., Hugelius, G., Koven, C. D., Kuhry, P., Lawrence, D. M., Natali, S. M., Olefeldt, D., Romanovsky, V. E., Schaefer, K., Turetsky, M. R., Treat, C. C., and Vonk, J. E.: Climate change and the permafrost carbon feedback, Nature, 520, 171–179,, 2015. 

Segers, R.: Methane production and methane consumption: a review of processes underlying wetland methane fluxes, Biogeochemistry, 41, 23–51,, 1998. 

Tang, G., Zheng, J., Xu, X., Yang, Z., Graham, D. E., Gu, B., Painter, S. L., and Thornton, P. E.: Biogeochemical modeling of CO2 and CH4 production in anoxic Arctic soil microcosms, Biogeosciences, 13, 5021–5041,, 2016. 

Tarnocai, C., Canadell, J. G., Schuur, E. A. G., Kuhry, P., Mazhitova, G., and Zimov, S.: Soil organic carbon pools in the northern circumpolar permafrost region, Global Biogeochem. Cy., 23, GB2023,, 2009. 

Thornton, P. E. and Rosenbloom, N. A.: Ecosystem model spin-up: Estimating steady state conditions in a coupled terrestrial carbon and nitrogen cycle model, Ecol. Model., 189, 25–48,, 2005. 

Tian, H., Xu, X., Liu, M., Ren, W., Zhang, C., Chen, G., and Lu, C.: Spatial and temporal patterns of CH4 and N2O fluxes in terrestrial ecosystems of North America during 1979–2008: application of a global biogeochemistry model, Biogeosciences, 7, 2673–2694,, 2010. 

Tilston, E. L., Sparrman, T., and Öquist, M. G.: Unfrozen water content moderates temperature dependence of sub-zero microbial respiration, Soil Biol. Biochem., 42, 1396–1407,, 2010. 

Tipping, E.: WHAM – A chemical equilibrium model and computer code for waters, sediments, and soils incorporating a discrete site/electrostatic model of ion-binding by humic substances, Comput. Geosci., 20, 973–1023,, 1994. 

Tipping, E.: Humic Ion-Binding Model VI: An Improved Description of the Interactions of Protons and Metal Ions with Humic Substances, Aquatic Geochem., 4, 3–47,, 1998. 

Treat, C. C., Wollheim, W. M., Varner, R. K., Grandy, A. S., Talbot, J., and Frolking, S.: Temperature and peat type control CO2 and CH4 production in Alaskan permafrost peats, Global Change Biol., 20, 2674–2686,, 2014. 

Treat, C. C., Natali, S. M., Ernakovich, J., Iversen, C. M., Lupascu, M., McGuire, A. D., Norby, R. J., Roy Chowdhury, T., Richter, A., Šantrůčková, H., Schädel, C., Schuur, E. A. G., Sloan, V. L., Turetsky, M. R., and Waldrop, M. P.: A pan-Arctic synthesis of CH4 and CO2 production from anoxic soil incubations, Global Change Biol., 21, 2787–2803,, 2015. 

Tuomi, M., Vanhala, P., Karhu, K., Fritze, H., and Liski, J.: Heterotrophic soil respiration – Comparison of different models describing its temperature dependence, Ecol. Model., 211, 182–190,, 2008. 

Vaughn, L. J. S., Conrad, M. E., Bill, M., and Torn, M. S.: Isotopic insights into methane production, oxidation, and emissions in Arctic polygon tundra, Global Change Biol., 22, 3487–3502,, 2016. 

Wachinger, G., Fiedler, S., Zepp, K., Gattinger, A., Sommer, M., and Roth, K.: Variability of soil methane production on the micro-scale: spatial association with hot spots of organic material and Archaeal populations, Soil Biol. Biochem., 32, 1121–1130,, 2000. 

Wagner, R., Zona, D., Oechel, W., and Lipson, D.: Microbial community structure and soil pH correspond to methane production in Arctic Alaska soils, Environ. Microbiol., 19, 3398–3410,, 2017. 

Waldrop, M. P., Wickland, K. P., White III, R., Berhe, A. A., Harden, J. W., and Romanovsky, V. E.: Molecular investigations into a globally important carbon pool: permafrost-protected carbon in Alaskan soils, Global Change Biol., 16, 2543–2554,, 2010. 

Walz, J., Knoblauch, C., Böhme, L., and Pfeiffer, E.-M.: Regulation of soil organic matter decomposition in permafrost-affected Siberian tundra soils – Impact of oxygen availability, freezing and thawing, temperature, and labile organic matter, Soil Biol. Biochem., 110, 34–43,, 2017. 

Xu, X., Elias, D. A., Graham, D. E., Phelps, T. J., Carroll, S. L., Wullschleger, S. D., and Thornton, P. E.: A microbial functional group-based module for simulating methane production and consumption: Application to an incubated permafrost soil, J. Geophys. Res.-Biogeosci., 120, 1315–1333,, 2015. 

Xu, X., Yuan, F., Hanson, P. J., Wullschleger, S. D., Thornton, P. E., Riley, W. J., Song, X., Graham, D. E., Song, C., and Tian, H.: Reviews and syntheses: Four decades of modeling methane cycling in terrestrial ecosystems, Biogeosciences, 13, 3735–3755,, 2016. 

Yang, S., Liebner, S., Winkel, M., Alawi, M., Horn, F., Dörfer, C., Ollivier, J., He, J.-s., Jin, H., Kühn, P., Schloter, M., Scholten, T., and Wagner, D.: In-depth analysis of core methanogenic communities from high elevation permafrost-affected wetlands, Soil Biol. Biochem., 111, 66–77,, 2017. 

Yang, Z., Wullschleger, S. D., Liang, L., Graham, D. E., and Gu, B.: Effects of warming on the degradation and production of low-molecular-weight labile organic carbon in an Arctic tundra soil, Soil Biol. Biochem., 95, 202–211,, 2016. 

Yvon-Durocher, G., Caffrey, J. M., Cescatti, A., Dossena, M., Giorgio, P. D., Gasol, J. M., Montoya, J. M., Pumpanen, J., Staehr, P. A., Trimmer, M., Woodward, G., and Allen, A. P.: Reconciling the temperature dependence of respiration across timescales and ecosystem types, Nature, 487, 472–476,, 2012. 

Yvon-Durocher, G., Allen, A. P., Bastviken, D., Conrad, R., Gudasz, C., St-Pierre, A., Thanh-Duc, N., and del Giorgio, P. A.: Methane fluxes show consistent temperature dependence across microbial to ecosystem scales, Nature, 507, 488–491,, 2014. 

Zhalnina, K., Dias, R., de Quadros, P. D., Davis-Richardson, A., Camargo, F. A., Clark, I. M., McGrath, S. P., Hirsch, P. R., and Triplett, E. W.: Soil pH determines microbial diversity and composition in the park grass experiment, Microb. Ecol., 69, 395–406,, 2015. 

Zheng, J., RoyChowdhury, T., and Graham, D. E.: CO2 and CH4 Production and CH4 Oxidation in Low Temperature Soil Incubations from Flat- and High-Centered Polygons, Barrow, Alaska, 2012, Next Generation Ecosystem Experiments Arctic Data Collection,, 2017. 

Zheng, J. and Graham, D. E.: Soil Organic Carbon Degradation in Low Temperature Soil Incubations from Flat- and High-Centered Polygons, Barrow, Alaska, 2012–2013, Next Generation Ecosystem Experiments Arctic Data Collection,, 2018. 

Zheng, J., RoyChowdhury, T., Herndon, E., Yang, Z., Gu, B., Wullschleger, S., and Graham, D. E.: Synthesis of Soil Geochemical Characteristics and Organic Carbon Degradation from Arctic Polygon Tundra, Barrow, Alaska, Next Generation Ecosystem Experiments Arctic Data Collection,, 2018a.  

Zheng, J., RoyChowdhury, T., Yang, Z., Gu, B., Wullschleger, S. D., and Graham, D. E.: Impacts of temperature and soil characteristics on methane production and oxidation in Arctic tundra, Biogeosciences, 15, 6621–6635,, 2018b. 

Zheng, J., Thornton, P., Painter, S., Gu, B., Wullschleger, S., and Graham, D. E.: Modeling Anaerobic Soil Organic Carbon Decomposition in Arctic Polygon Tundra: Insights into Soil Geochemical Influences on Carbon Mineralization: Modeling Archive, Next Generation Ecosystem Experiments Arctic Data Collection,, 2018c. 

Zhuang, Q., Melillo, J. M., Sarofim, M. C., Kicklighter, D. W., McGuire, A. D., Felzer, B. S., Sokolov, A., Prinn, R. G., Steudler, P. A., and Hu, S.: CO2 and CH4 exchanges between land ecosystems and the atmosphere in northern high latitudes over the 21st century, Geophys. Res. Lett., 33, L17403,, 2006. 

Zimov, S. A., Schuur, E. A. G., and Chapin, F. S.: Permafrost and the global carbon budget, Science, 312, 1612–1613,, 2006. 

Short summary
Arctic warming exposes soil carbon to increased degradation, increasing CO2 and CH4 emissions. Models underrepresent anaerobic decomposition that predominates wet soils. We simulated microbial growth, pH regulation, and anaerobic carbon decomposition in a new model, parameterized and validated with prior soil incubation data. The model accurately simulated CO2 production and strong influences of water content, pH, methanogen biomass, and competing electron acceptors on CH4 production.
Final-revised paper