Mechanistic modeling of microbial interactions at pore to profile scale resolve methane emission dynamics from permafrost soil

2.1 Representation of Soil Structural Features: Pores to Profile Scales

We discretized the soil profile into layers of thickness, Δz which contains bulk soil and embedded aggregates of different sizes located in bulk soil and subjected to macroscopic conditions (temperature, hydration, and oxygen) that vary with depth. The bulk soil and aggregate domains inhabiting microbial life were represented by 2-D and 3-D angular pore networks, respectively [Ebrahimi et al., 2013; Ebrahimi and Or, 2014] (see Figure 1). Similar angular pore networks have been used to model effects of hydration conditions on oxygen/carbon diffusion processes in soil aggregates of different sizes [Ebrahimi and Or, 2014, 2015]. The aggregate size distribution at each soil depth is represented by lognormal distribution as often observed in natural soils [Gardner, 1956; Ebrahimi and Or, 2016]. To keep the interpretation sufficiently general, these structural domains could also be considered as representative of generic hot spots centered around a temporally anoxic core and thus requiring no particular pore space characteristics (distinguished from bulk soil by potentially different pore properties or presence of exopolymers and particulate organic matter in their core) as often ascribed to real soil aggregates.

The soil aggregate size distribution is often a function of soil texture, organic matter content, vegetation, and freezing/thawing processes. Here we have chosen four aggregate size classes (1, 2, 3, and 4 mm); these were selected larger than 1 mm to ensure development of persistent anoxic microsites [Ebrahimi and Or, 2016]. The associated angular pore networks hosting microbial life consisted of between 10³ and 30³ bonds to simulate aggregate sizes of physical dimensions in the range from 1 to 4 mm in diameter. The bond length is derived to maintain equal porosity in all aggregate sizes (using tortuous helical bonds as described in Ebrahimi and Or [2015]). The population mean aggregate size was chosen as 2 mm to conform with common values in the literature [Schneider and Gupta, 1985; Six et al., 2000] (calculations are based on Ebrahimi and Or [2016]).The soil properties (e.g., water characteristics, porosity, and pore size distribution) used for pore network simulations are chosen to represent parameters for peat soil to enable comparisons with observations (see Table S2 in the supporting information). Peat properties widely vary due to several factors (e.g., bulk density, organic matter content, presence of plant roots, and soil compaction) [Schwarzel et al., 2002; Nimmo, 2004; Horn and Smucker, 2005]. In this study, we assumed constant bulk density and the parameters sampled from the range of typical observed values (see Weiss et al. [1998] for detailed analyses and Table S2). The mean pore sizes of the pore networks for bulk soil are estimated from the water characteristic curves of the peat soil. For the simulations, the mean pore sizes of aggregates are assumed 5 times smaller than the bulk soil. The macroscopic conditions experienced by the 2-D and 3-D pore networks at a given soil layer were kept similar (e.g., temperature, water, and oxygen at the aggregate boundary). No nutrient interactions and gas exchange among aggregates in the same layer are considered except the consumption of methane flux by methanotrophs inhabiting the aerated 2-D pore networks (representing life in aerated bulk soil). More details on the upscaling of biogeochemical gas fluxes produced in each member of a population of soil aggregates of different sizes and at different depths to fluxes at the soil surface are provided in Ebrahimi and Or [2016]. A brief description on the nutrient (carbon for methanogens and CH4 for methanotrophs) and oxygen diffusions within single soil aggregates are represented in the Text S1.

2.2 Methanogenic and Methanotrophic Microbial Activity in Thawing Permafrost

We built on past experiences with implementing individual-based models (IBM) [Kreft et al., 1998; Wang and Or, 2010; Ebrahimi et al., 2013] to simulate microbial community composition, spatial distribution, and functioning in simplified soil pore spaces. The biophysical model resembles the formulation used in Ebrahimi and Or [2015] that explicitly considers bacterial cell dispersal, growth, division, and nutrient/oxygen consumption according to physiological traits associated with methanogens and methanotrophs as explained next. In this study, we have only considered methane oxidation by aerobic methanotrophs that are known to be dominant in peatlands and wetland soils [Segers, 1998; Strous, 2010]. However, it should be mentioned that the recent evidences suggest a potentially important role of anaerobic methanotrphs (linked to microbial sulfate reduction) in various peatland ecosystems [Smemo and Yavitt, 2007, 2011]. Further experimental studies are required to investigate the significance of anaerobic methane oxidation and its relationship with other known anaerobic pathways (e.g., denitrification, sulfate, Fe, and Mn reductions) [Raghoebarsing et al., 2006; Smemo and Yavitt, 2011].

2.2.1 Microbial Growth Kinetics

A modified Monod-kinetic parameterization was used to describe the growth of methanogens and methanotrophs microbial cells (and emergent communities). The model considers growth inhibition by presence of oxygen for methanogens. Both methane and oxygen were considered as growth limiting “nutrients” for methanotrophs. The growth kinetics of individual cells in the simulated pore networks were represented as follows:

For methanogens:

(1)

and for methanotrophs:

(2)

where [C], [CH₄] and [O₂] are carbon, methane, and oxygen concentrations (mg/L), respectively; V_bond [L] is the aqueous volume of each bond of pore network and s[mg] = ρv_cell is the cell dry mass of the individual bacteria calculated from multiplying the cell density, ρ[mg/L] by cell volume, v_cell[L]. The terms ν_[C] and are dissolved organic carbon and methane uptake rates, respectively. (maximum specific growth rate/growth yield of methanotrophs) and (maximum specific growth rate/growth yield of methanogens) are the maximum specific nutrient uptake rates of methanotrophs and methanogens, respectively. For the special conditions of permafrost, we represent the maximum specific growth rate, μ_max as a function of ambient temperature, T [K] according to a modified Ratkowsky model [Ratkowsky et al., 1982; Zwietering et al., 1991],

(3)

where b₃ and c₃ are fitting parameters obtained experimentally. We note that equations 1 and 2 are applied for each individual cell; hence, the growth dynamics of the respective populations would reflect spatial and temporal limitations imposed by nutrients (carbon/methane) or oxygen that are expected to limit the size and position of these interacting communities. The physiological parameters of the kinetic microbial model are summarized in Table S1.

The simulations were initialized by introducing 1000 methanotrophs and 1000 methanogen cells uniformly inoculated within bonds of each pore networks. After inoculation of individual cells, they disperse within the entire network, produce enzymes, consume nutrients (carbon for methanogens and CH₄ and O₂ for methanotrophs), grow, and divide to new daughter cells. The net growths of individual cells (μ_net) are calculated based on the actual uptake rate of nutrient ( ) and with respect to the maintenance of individual cells (m_i). In this study, a proportional relationship between the maintenance rate and the maximum growth rate (m = αμ_max) is considered [Kim and Or, 2016]; the net growth in individual cell mass (μ_net) can be then described as follows:

(4)

The minimum volume of the cell (V_d , min) for division to new daughter cells could be estimated from the descriptive Donachie model ( ) represented by Kreft et al. [1998]. is median volume of the individual cell. The individual cell is considered to be dead if the cell volume becomes smaller than the threshold minimum volume (V_min) which is assumed to be one fifth of V_d , min [Kreft et al., 1998].

2.2.2 Microbial Motility and Dispersal

For a range of hydration conditions, flagellated motility is the main dispersal mechanism for soil microbial cells [Dechesne et al., 2010; Wang and Or, 2010; Tecon and Or, 2016]. We employed a previously developed motility model [Wang and Or, 2010; Ebrahimi and Or, 2014] that considers the balance between self-propulsion force for flagellated cells (F_M) in bulk water with interfacial drag forces (cell-wall hydrodynamic interactions, F_λ, and capillary pinning, F_Ca) as a function of water film thickness of bond, δ. The resulting local cell velocity within a single bond in the pore network is assumed to be proportional to the cell velocity in the bulk water, V₀ (considered here as 1 μm s⁻¹ reference) and expressed as

(5)

The water film thickness, δ, for each bond in the network is defined as a function of matric potential (water pressure, ψ) and bond angularity, α [Ebrahimi and Or, 2014].

An important factor for self-organization is the assumption that microbial motility is guided by chemotaxis toward higher nutrient concentrations and/or to move away from unfavorable chemical components [Brosilow et al., 1996; Olson et al., 2004; Ford and Harvey, 2007]. In this study, we applied chemotactic motion to guide microbial communities to self-organize in bulk soil and aggregates in which methanotrophs are attracted for regions with higher oxygen and methane concentrations while movement of methanogens are biased to avoid high oxygen concentrations. Localization of methanotrophs and methanogens along the radius of aggregates (hot spots) is shown in Figure S2 in the supporting information.

2.3 Variations of Active Layer Depth With Temperature: Freezing/Thawing Processes

2.3.1 Soil Thawing Model

We expressed the dynamics of the thawed active layer thickness in permafrost soil as a function of soil surface temperature and its water content using the model of Kurylyk et al. [2014]. The model considers different thermal properties and two distinct transient heat conduction equations as follows:

For thawed zone:

(6)

For frozen zone

(7)

where T and are the temperature distribution in the thawed and frozen zones (°C). The α and α_f are the bulk thermal diffusivity of the thawed and frozen zones (m² s⁻¹), respectively, and Z is the position of the thawing front. The initial condition is uniform (constant) temperature over the soil profile (T(z, t = 0) = T_i), as shown in Figure 2. The bottom boundary condition is assigned as constant temperature condition (T(t) = T_f) at the interface of frozen and thawed zones (Figure 2).

The energy balance at the interface between the frozen and thawed zones matches the conductive heat flux from the thawed zone to the amount of energy required for thawing as the enthalpy of fusion and the conductive heat flux to the frozen zone [Kurylyk et al., 2014] (see Figure 2 for the illustrations); therefore, the energy balance can be formulated as

(8)

where S_wf is the liquid water saturation in the thawed zone (originally frozen), ρ_w is the liquid water density, ε is the soil porosity, and L_f is the enthalpy of fusion for water. The λ and λ_f are the bulk thermal conductivity of the thawed and frozen zones, respectively. The Neumann solution for the interface (equation 8) is typically expressed as follows

(9)

where m is the coefficient of proportionality and can be calculated by equaling Y₁ and Y₂ that are extracted from the solution of equation 3 as given

(10)

The temperature profile in the thawed zone can be then calculated from solution of equation 6 with respect to the top and bottom boundary conditions. The top and bottom boundary conditions are defined as the soil surface temperature, T_s, and the temperature of frozen permafrost at the interface, T_f. For the purpose of simulations, T_f is assumed to be zero where the observation data are not available. The top boundary condition (T_s) is taken from field observations (mean daily values) neglecting diurnal temperature variations. To evaluate the performance of the model, in Figure S1 the thaw penetrations in the frozen permafrost for four different surface temperatures are plotted. The original model has been successfully validated against experimental and numerical simulations [see Kurylyk et al., 2014]. The model is solved one dimensional over the soil depth, and lateral fluxes are not considered. For simplicity, homogenous soil properties based on bulk soil is considered and the contribution of soil aggregates on freezing/thawing process is not applied. In the simulations water content profile is assumed to be under quasi-hydrostatic conditions parameterized by the van Genuchten model [van Genuchten, 1980] (see Figure S1). The scale-appropriate estimation of water table dynamics requires the implementation of an extensive hydrological model at larger scales [Pomeroy et al., 2007; Swenson et al., 2012] that should include the role of rapid changes in the landscape of thawing permafrost (e.g., ice-wedge degradation) [Liljedahl et al., 2016], and it is well beyond the scope of this study. For estimation of seasonal variations in methane emissions (sections 16 and 17), the water table level is assigned from field observations using daily average values (Figure S3).

2.3.2 Soil Refreezing Model

The refreezing process during the arctic fall season was represented by a simplified solution of the energy and mass equations for unsaturated soil. It was then integrated with the proposed biophysical model to quantify impact of refreezing on methane emission rates. The water state in the soil profile was represented by 1-D mass conservation model given as

(11)

where M_w and M_i (kg) are the mass of water and ice, respectively. S_w is a sink term (water loss due to evaporation or root water uptake), and J_w (m s⁻¹) is the water flux obtained by the Richards equation. Here for simplicity, we assumed that S_w and J_w are equal to zero (no drainage is considered). The energy conservation in a soil volume is expressed as

(12)

and .where U is the sum of internal energy of the solid particles, ice, and liquid phases; q_c, J_adv., and S_loss are the conduction, convection heat fluxes, and a sink term due to energy losses to soil depth, respectively; and λ_T is the macroscopic thermal conductivity of the soil a function of its water and ice contents and temperature [Hansson et al., 2004; Dall'Amico et al., 2011]. For an unfrozen profile, the two primary unknowns are the water content, θ, and temperature profiles that are obtained directly from mass and energy conversion equations (equations 11 and 12). The mass balance equation (equation 11) is solved for the initial (before starting of the freezing process) conditions of water content profile driven from van Genuchten model by assuming quasi-hydrostatic conditions. The temperature profile is updated by considering the soil surface temperature, T_s as a top boundary condition, and the temperature at the interface of frozen and nonfrozen zones as a bottom boundary condition (assumed T_f = 0°C in the simulations). Similar to the thawing process, the energy and mass balance equations are solved for 1-D variations over the soil depth with consideration of bulk soil hydraulic and thermal properties.

Changes in the soil ice content as a function of temperature during freezing require additional closure relation. The generalized Clapeyron equation is used to represent the variation of water pressure (matric potential, ψ) as a function of temperature during phase changes [Dall'Amico et al., 2011],

(13)

Since the surface tension at air-water interfaces depresses the water melting temperature, T^*, relative to that of bulk water at atmospheric pressure (T_m = 273.15 K), we calculated the water pressure of the saturated and unsaturated soils separately by integration of equation 13:

For the saturated domain in the soil profile,

(14)

where p_a is the relative atmospheric pressure which can be set to zero. By mathematical simplifications, the water pressure as function of temperature is expressed as

(15)

The melting temperature of unsaturated soil is obtained by integration of equation 13 from T_m to T^* and matric potential, ψ, from atmospheric pressure p_a to ψ₀:

(16)

ψ₀ corresponds to the matric potential of the total water content (summation of liquid and ice contents). Using certain mathematical approximations [Dall'Amico et al., 2011], T^* at unsaturated conditions can be expressed as

(17)

T* defines the threshold temperature that below this temperature, the soil is under freezing conditions. Again, by integrating equation 13, pore water pressure is calculated as

(18)

Using van Genuchten model, the water pressure can be converted to equivalent liquid water content and thus deduce the ice content.

2.4 Scaling Up Methane Fluxes From Individual Aggregates to Emissions From the Soil Surface

The proposed biophysical model for quantifying methane production and consumption in each soil layer (consisting of bulk and aggregate soil fractions) is upscaled to simulate methane production, consumption, and net transport over the soil profile. The model is coupled with the permafrost thermal (freezing/thawing) and hydrological processes (primarily variations of water table level and liquid saturation along the soil depth) to provide the depth-dependent macroscopic boundary conditions for oxygen, water, and temperature over the soil profile. A daily time step for the thermal boundary conditions and water table level was considered, while microbial activity and nutrient diffusion processes (at the aggregate and profile scales) are simulated at constant time steps of 30 min. Each layer within soil profile (thickness of Δz) contains a volume fraction of soil aggregated (f_agg.) and a bulk soil fraction (f_bulk = 1 − f_agg.). The thickness Δz was selected to accommodate the different aggregate sizes provided from an aggregate probability distribution function. The total gas production or consumption rates at each layer are obtained by integrating the production/consumption rates (S_agg.) over aggregates of different sizes and considering the fraction of bulk soil (S_bulk) to give the gas production/consumption rate from a specific layer.

The primary transport mechanisms of methane involve gas diffusion through the unsaturated soil, gas bubble ebullition from supersaturated soil layers and aggregates, and plant-mediated passive transport. The dominant transport mechanism vary with soil structure (e.g., texture, aggregate, and layering), environmental factors (e.g., temperature and water table position), and vegetation type and density [Blodau, 2002; Chanton, 2005; Klapstein et al., 2014]. In this study, the presence of vegetation as transport mechanism has not been considered, and estimates are based on diffusion and ebullition as main transport processes.

The roles of vegetation and the plant root distributions would affect soil structure (e.g., biopores, aggregate stability, and size distribution) and rates of organic matter accumulation [Six et al., 1998; Or and Ghezzehei, 2002; Davidson et al., 2016; Ebrahimi and Or, 2016]. Additionally, vegetation type serves as an indicator of microbial community composition and a potential pathway for CH4 transport in permafrost soil [King et al., 1998; Christensen et al., 2004; Olivas et al., 2010; McEwing et al., 2015; Davidson et al., 2016]. However, mechanistic understanding of the processes that shape microbial communities under different vegetation types remain sketchy, and we thus opted not to include plant root processes in the modeling and upscaling presented in this study. More work is needed to establish mechanistic understanding of the links between microbial processes at the soil pore scale and characteristic properties of plant roots and their rhizospheres.

2.4.1 Methane Diffusive Fluxes in Soil

The diffusion of CH₄ through the soil is an important pathway that supports aerobic methanotrophic activity near the permafrost surface oxidizing methane to CO₂. To estimate the net methane diffusive flux, we solve a diffusion-reaction equation over each soil layer (d is the layer index):

(19)

where is the flux accumulation term within the layer.

C_gas and C_liq. are methane concentrations in gas and liquid phases, respectively. It is assumed that the gas and liquid phases are in equilibrium, and thus, Henry's law determines the partitioning between these phases, C_gas = HC_liq.. The total methane diffusive flux from aggregates and bulk soil for a given layer S_{d , diff.} is given as

(20)

where f(R) follows lognormal distributions found in experimental studies [Gardner, 1956; Ebrahimi and Or, 2016] and R is the aggregate size. It is assumed that gaseous fluxes are one-dimensional. The calculated flux for each layer (S_{d , diff.}) is then used as sink or source term for gas diffusion-reaction equation (equation 19) along the soil profile to quantify the emitted diffusive gas fluxes from the soil surface.

From Fick's law, fluxes are written as

(21)

and

(22)

where A_liq. and A_gas are cross-sectional area of soil profile containing liquid and gas phases, respectively; this can be calculated from liquid (φ_liq.) and gas (φ_gas) filled pore space fractions and total cross sectional areal of soil profile (A_tot.):

(23)

φ_liq. and φ_gas for each soil depth are calculated from the water characteristic curve. Finally, equation 19 can be then explicitly solved by finite difference method to get the transient gas diffusion along the soil profile.

2.4.2 Methane Bubble Ebullition

Gas bubble ebullition is considered a significant transport pathway for emitted methane fluxes [Schütz et al., 1989; Sharkey et al., 1991]. The low value of methane solubility in water combined with limited diffusion rates in soil solution (or saturated pores) results in frequent bubble nucleation within soil pores and aggregates. The methane concentration in each bond in the pore network is evaluated at each time step, and when the concentration exceeds a threshold value, a methane bubble nucleates and grows within bulk or aggregated soil pores. The threshold methane concentration is assumed to be methane solubility in water that varies exponentially with temperature (T):

(24)

The mechanics of bubble release from soil pores (buoyancy forces overcoming capillary resistive forces) are too complex to consider explicitly within the modeled pore networks [Chen and Slater, 2016; Painter et al., 2016]. Recent studies highlight the important role of pore structure (of peat) in regulating bubble storage and release dynamics with characteristically intermittent ebullition events in time and heterogeneous in space [Glaser, 2004; Green and Baird, 2013; Ramirez et al., 2016]. Evidence suggests that changes in overlaying water pressure (e.g., water table or tides) affect rates of gas bubble ebullition from the soil or sediments [Chanton et al., 1989]. Clearly, additional mechanistic models are needed to better describe effects of pore structure and pathways on bubble entrapment, storage, and release. Additionally, external pressure (water table and barometric fluctuations) and ambient temperatures are known to affect the extent and frequency of bubble expansion, release thresholds, and subsequent motion [Algar and Boudreau, 2010; Ramirez et al., 2015; Chen and Slater, 2016]. We note that the impact of soil freezing on bubble ebullition remains poorly understood and affects emission predictions during fall season and will be elaborated in the results section.

As first approximation for representing an intermittent ebullition process, we have used field observations that suggest that ebullition is associated with a certain range of bubble sizes represented macroscopically by a threshold (methane) volume fraction in soil pores [Tokida et al., 2005; Ramirez et al., 2016]. The release threshold for an ebullition event was taken from a uniform distribution of threshold values in the of 0.08 to 0.2 [Rothfuss and Conrad, 1998; Tokida et al., 2005] for both aggregates and bulk soil, creating a stochastic link between methane gas storage in the soil and release dynamics. For methane storage in a model domain that exceeds the release threshold (stochastically assigned from the range above) the entire storage methane gas is instantaneously released, travels to the soil surface, and is emitted to the atmosphere. Hence, the total methane flux at the soil surface includes a relatively steady contribution of a diffusive flux (calculated from equation 11) and superimposed abrupt bursts of methane bubble released. We consider the travel of methane bubbles bypassing the aerobic zone too fast to be intercepted and used by the aerobic methanotrophs. The contribution of entrapped methane bubbles to the methane diffusion process is neglected.

3.1 The Role of Water Table Position on Methane Production and Consumption

We begin with a test of features of the biophysical model considering methane emissions from peatlands to illustrate effects of hydration conditions (represented by water table position) in soil profile and stoichiometric constraints (C:N ratio) on methanogen and methanotroph community size and activity and resulting methane emissions from the soil surface. The simulations were motivated by laboratory experiments of peatland soil columns conducted by Moore and Dalva [1993]. A range of static water tables from 0 to 60 cm are studied, and soil column length is considered 80 cm, similar to the experiments. Water content profile is assumed to be under quasi hydrostatic conditions (the hydraulic properties of the peat soil are given in Table S2), and simulations are performed under constant temperature at 22.6°C. The model assumes that 30% of the soil matrix is aggregated. Nitrogen content (determined by the ratio of C:N) is considered to be the limiting factor for enzymatic activity and resulting microbial growth. The peat soil is assumed to be carbon rich in the simulations.

Figure 3 illustrates the key effects of water table position on the rates of methane emissions. The simulations are done for three different C:N ratios (10, 20, and 80) which covers broad range of observed C:N values in different biomes [Aitkenhead and McDowell, 2000] including peatland soils that are known to have high values of C:N ratios [Klemedtsson et al., 2005]. The simulation results show that lowering the water table level and increasing the oxygenated fraction of the soil profile reduce the net methane emissions from soil surface. This is attributed to changes in relative size and activity of methanogenic and methanotrophic microbial communities in response to lowering of water table position. As indicated in Figure 3, the methane emission rate is significantly influenced by C:N ratio. This highlights the role of nitrogen availability that adds a stoichiometric constraint for enzyme activity and thus limits microbial growth.

The roles of chemical controls (e.g., C:N and nutrient sources) and the stoichiometric constraints imposed on both microbial and enzyme activities are rarely represented in current models, primarily due to mathematical difficulties of considering such biogeochemical process dynamics and associated diffusion-reactions at the relevant scales for microbial activity [Zhuang et al., 2004; Bridgham et al., 2013]. The reasonable quantitative agreement obtained between our simulation results and laboratory experiments inspires confidence in implementing the model to study the more challenging questions of the role of soil structure (e.g., aggregate) on microbial (methanogen and methanotroph) community dynamics and the resulting methane production rates.

3.2 Soil Structure Effects on Microbial Community Dynamics and Methane Emissions

We evaluated microbial community dynamics under similar conditions (soil type and boundary conditions) as in previous section with a fixed C:N ratio of 10 and soil aggregate fraction of 30% in the peatland soil (depicted in Figure 4). The simulations are conducted for a range of water table positions from 0 to 0.8 m. The purpose of this model test was to illustrate shifts in the composition and spatial distribution of methanogens and methanotrophs and the sources and sinks of methane production and consumption with changes in water table position (considering persistent anoxic conditions within aggregates even in a fully aerated soil [see Ebrahimi and Or, 2016]).

Figure 4 depicts methane production and consumption rates over the soil profiles as a function of water table level for bulk and aggregated fractions of soil. The results are for the stationary phase of microbial activity. The evolution of the microbial methanogens and methanotrophs from numerical inoculation time is shown in the subset plots of Figure 4 for the water table at depth of 0.4 m. The simulation results suggest that methane production and consumption are higher in aggregated soil for a range of water table positions. Methane production in the aerated soil (low water table position) is completely disappeared in bulk soil while it can be still observed in the aggregates. In the bulk soil, methanogenic activity is only limited to the saturated zone of the soil profile where oxygen content is limited, and therefore, the methane production significantly drops in the aerated bulk soil. Similarly, the methanotrophic activity of bulk soil is only flourished in a narrow zone of soil profile near the interface between water-saturated and -unsaturated soil depth where the oxygen is available and water content is high enough for the purpose of nutrient (dissolved methane) diffusion. Subset of Figures 4a and 4b show the spatial distribution of methanogens and methanotrophs over the soil depth (bulk versus aggregate), respectively.

The 3-D pore structure associated with soil aggregates (or elevated microbial activity around plant residues embedded in soil matrix—a hot spot) gives rise to formation of oxic and anoxic microsites for a range of hydration conditions. Such microsites support coexistence of aerobic (e.g., respiration and methane oxidation) and anaerobic (e.g., denitrification and methane production) microbial activity even in aerated soils [Tiedje et al., 1984; Ebrahimi and Or, 2015]. As Figure 4 shows, methanotrophs and methanogens spatially organize themselves along the aggregate radius where methanotophs are located near the aggregate surface to efficiently consume oxygen (diffusing from the surface) and methane (produced by methanogens in the anoxic core of the aggregates). Due to this interesting feature of soil aggregates, the methanogenic activity is more resistant to changes of the water table and shows relatively homogenous distribution over the soil depth (see subset of Figure 4a). A similar effect can be seen in the oxic layer of soil profile for uniform distribution of methanotrophs within soil aggregates, while methanotrophs in bulk soil is localized to a narrow zone of soil depth (see subset of Figure 4b). These simulation results offer mechanistic explanations to experimentally observed variations in the methanotophic and methanogenic community compositions over the soil depth and for bulk soil and aggregates [Macalady et al., 2002; Shrestha et al., 2008; Lee et al., 2015]. Additionally, 3-D pore networks (aggregates) could be used to represent rhizospheric processes and improve understanding of microbial activity in such environments (e.g., rice paddies) known to host more methanogens and produce more methane than in bulk soil [Großkopf et al., 1998; Ramakrishnan et al., 2000].

3.3 Seasonal Variations of Methane Emission

The successful evaluation of the new features of the biophysical model related to microbial community size and spatial dynamics, hydration effects on methane production, and stoichiometric constrains for activity in carbon-rich soil permits study of the core questions of the study, namely, methane emissions from active layers of permafrost soil. It well established that soil methane emissions are strongly influenced by the seasonal dynamics of air temperature, active layer thickness, and water table position. Yet little is known regarding the production, oxidation, and storage mechanisms that directly shape methane fluxes during the spring melt and fall freezeup [Hargreaves et al., 2001; Mastepanov et al., 2008; Hanis et al., 2013].

The modeling framework was tested to predict seasonal variations in methane emissions based on two reported case studies. Both case studies were located at high-arctic regions in Zackenberg valley [Meltofte and Rasch, 2008], the Northeast Greenland National Park (74°30′N, 21°00′W). We have used field measurements for 2007 and 2010 growing seasons and subsequent freeze-in periods [Mastepanov et al., 2008, 2013]. The study site is a patterned fen often characterized by alternating regions of high and low topographies (perpendicular to slope or groundwater flow) with the hollows (flarks) often wet and inundated (see Tagesson et al. [2012] and Meltofte and Rasch [2008] for more information). The average monthly air temperature of this site is less than −20°C in winter and between +3 to +7°C for summer period, and the area has experienced a significant warming of 2.25°C, between 1991 and 2005 [Hansen et al., 2008]. The seasonal variations of active layer depth and water table levels are shown in Figure S3. The variations of soil surface temperature (taken from field observations [Mastepanov et al., 2013]) were used as top boundary conditions for the heat conduction equation (equation 6) and yield seasonal variations in soil temperature over depth, as represented in Figure S4. The thaw depth and the water table level were additional input parameters for our simulations, and these were taken from observations of Mastepanov et al. [2013], as shown in Figure S3. We have used peat soil thermal and hydraulic properties to perform the simulations (Table S2) and considered 30% of each soil layer to be composed of aggregates (or anoxic hot spots).

We note that for some measurements reported in the literature, the observed peak in methane emissions upon thawing in the early spring is often attributed to the release of stored methane from the previous year and methane production under low winter temperatures [Heyer et al., 2002; Song et al., 2012]. In this study, we assumed that simulations begin in the first day of thawing with depleted soil methane in the profile. Additional studies are needed to capture early spring emissions by updating initial conditions for methane concentration profile from previous year or by running multiyear simulations.

Simulation results for methane emission rates for the two case studies as affected by seasonal variations in environmental conditions (e.g., temperature, water table, and thaw depth) are depicted in Figure 5. The simulation results show reasonable agreement with field observations. Two different seasons were considered, where case study 1 represents high water table level (associated with a “wet” year), whereas case study 2 represents lower water table level (“dry” year) during the growing season. The simulation results are consistent with literature observations that report positive correlation between the increase in soil thawing depth and methane production rates [Hodgkins et al., 2014; Corbett et al., 2015]. We examined the sensitivity of the simulations to the assumed mean aggregate sizes in the profile. The results suggest (Figure 5) that increasing mean aggregate size (could be interpreted as mean anoxic hot spot size) from 2 to 3 mm resulted in a significant increase in methane emissions (especially during the fall refreezing in the dry year). The simulation results suggest great sensitivity of simulated methane emissions to the size and extent of anoxic volumes in the thawed soil (presently, a poorly defined variable).

Figure 6 shows that for case study 2 (dry year), methane production shifted deeper into the soil profile due to limited water content and establishment of oxic conditions at the topsoil. However, the naturally lower temperatures at deeper soil layers tend to suppress overall microbial activity and resulted in lower methane production rates (Figure 6b). Variations in the water table significantly modify the methane concentration profile and influence the evolution of methane storage rate in the form of gas bubbles (as depicted in Figure 6d). Drier soil promotes gas diffusion process which in turn reduces the methane storage rate and thus methane ebullition rate. For the simulated case studies, results indicate that the methanotrophic activity does not play a significant role on the rate of methane emission, since their activity is limited to a narrow time window where the water table is near surface and temperature is high enough. Their activity is suppressed when the water table drops to lower soil layers with low temperature (see Figure 6c for methane consumption rates due to methanotrophic activity).

3.4 Quantifying the Late Season Emissions of Methane From Permafrost Soil

The limited field measurements of methane fluxes during the arctic fall season show peaks of methane emissions that in some cases contribute significantly to the annual methane emissions [Mastepanov et al., 2008; Wille et al., 2008; Zona et al., 2016]. Late fall and early spring season methane emission peaks are attributed to the release of stored methane that has been produced during the growing season [Friborg et al., 1997; Heyer et al., 2002; Tokida et al., 2007a; Song et al., 2012] and to prolonged methanogenic activity in the refrozen soil supported by premeltwater films [Rivkina et al., 2000, 2007; Steven et al., 2007; Mackelprang et al., 2011]. The limited field data preclude definitive inferences of the primary mechanisms involved in the methane production and release during the fall season [Tokida et al., 2007b; Mastepanov et al., 2008; Zona et al., 2016]. In this study, we aim to introduce the effects of thawing and freezing of soil on methane production and transport processes that are missing in the current process-based models and to provide some insight to characterize the late season emissions.

To quantify late season peak methane emissions, we assume that the freezing of the liquid phase expels dissolved methane into gas form (methane hydrate formation is unlikely under atmospheric pressure [Maekawa et al., 1995]). Next, the accumulation of methane nucleates into gas bubbles within the corresponding soil layer (Figure 6d); hence, the freezing of the liquid phase promotes formation of methane bubbles. An inventory of gas bubbles within each aggregate and in the bulk soil of a corresponding soil layer is maintained. Methane gas bubbles are released to the atmosphere at a rate determined by bubble volume fraction threshold (in aggregates or bulk soil). The threshold values was based on the observations and was chosen from a uniform distribution in the range of 8 to 20% of the soil volume [Tokida et al., 2005] (see section 15). Note that in the absence of direct observations, we have used the same volume fraction thresholds for bulk and aggregate fractions in a soil layer. Figure 5 depicts the late fall emissions for the wet and dry scenarios. The results are consistent with the field observations despite the stochastic emissions due to methane ebullition events which add uncertainty on the observation data. The results show that the late emission rate is largely dependent on the storage of dissolved methane over the growing season that due to soil freezing process and reduction in the liquid phase content converts to gas bubbles and transfers to the atmosphere. Interestingly, because the solubility of methane increases with lower temperatures (see equation 21, based on Henry's law), the storage of dissolved methane increases before the freezing period (due to temperature drop) and thus promote methane late peak emissions.

Our findings provide mechanistic understating to this phenomenon that has been attributed to various processes such as (i) a moving freezing front that induces high pressures and “squeezed” gas bubbles out through vascular plant tissue [Mastepanov et al., 2013; Pirk et al., 2015] to (ii) double freezing fronts in which the unfrozen soil layer is sandwiched between two freezing fronts during which methanogenic activity persists and culminate in the release of a large methane peak [Zona et al., 2016]. Our simulation results suggest that a large portion of the methane is dissolved in the liquid phase before the onset of the freezing season (Figure 6d); hence, the notion of squeezed gas bubbles through plant tissue (while plausible) would require special considerations of capturing methane bubbles by plants during the freezing period. We argue not only that the freezing of unsaturated soil would not produce squeezing pressure in the unfrozen layer (gas pathways to the atmosphere remain open), but we expect methanogenic activity to be suppressed with lower temperatures (as Figure 6b suggests) and thus contribute minimally to late season peak emissions (in contrast to postulates by Zona et al. [2016]). Interestingly, late season emissions were more pronounced during the “wet year” due to generally low gas diffusion rates through the wet soil profile that may have resulted in higher accumulation of dissolved methane in the liquid phase (Figure 6a) relative to “dry year.” This occurs on top promotion of methanogenic activity due to persistence of anaerobic conditions during the wet year (for similar reasons that limit oxygen diffusion).

Field observations suggest that considerable amounts of methane could be stored in the soil during refreezing in the fall [Walter et al., 2006; Rinne et al., 2007; Merbold et al., 2009] without exhibiting a late season emission peak observed elsewhere and simulated in this study [Mastepanov et al., 2013; Zona et al., 2016]. These and other studies suggest that the two scenarios of large release of methane pulse in the fall or entrapment of large amounts of methane upon soil refreezing are possible. The conditions for these two different outcomes with respect to methane emissions seem to be determined by the position of the water table during the onset of soil refreezing and by seasonal temperature dynamics that determine the evolution of methanogenic activity and methane storage capacity. Evidence suggests that if refreezing process starts while the soil is flooded (high water table near the surface), a double freezing front from the soil surface propagates toward the permanent permafrost bottom layer leads to methane entrapment during the fall. The methane would then only releases with the following spring thaw [Rinne et al., 2007; Merbold et al., 2009]. The field data that motivated simulations in this study paint a different scenario where the soil profile becomes unsaturated (Figure S3) and the advancing soil freezing front permits escape of methane gas bubbles (released from various pockets and hot spots) contributing to the large peak emission during the fall season.

To evaluate the model performance at the regional scale, we tested the model with data obtained from Eddy covariance measurements [Zona et al., 2016] on the north slope of Alaska, Barrow, as shown in Figure S5. We have used averaged thaw depth, water table, and soil surface temperatures as inputs for the simulations based on field observations by Zona et al., 2016. The hydrological and thermal properties of peat soil are used for the simulations as listed in Table S2. The simulation results are in fair agreement with field observations despite differences in profile scale variations in methane gas fluxes not captured in averaged regional-scale observations. The model results are consistent with the field observations in showing the importance of late season emissions that could account for about 50% of the annual emissions. The role of late season emission in arctic regions remains unclear requiring better and more systematic understanding for improving current process-based models [Tang et al., 2016; Xu et al., 2016].

Before attempting to explain discrepancies between model simulations and measurements, we should describe the inherent limitations of the model and the uncertain boundary conditions imposed relative to field conditions with local variations in hydration, hydraulic and thermal properties, surface energy balance, vegetation, and other local variations. Nevertheless, it is instructive to separate the discrepancies to the growing season and the fall refreezing. During the growing season, the model captured the field observations reasonably (especially before the water table drops significantly). Following changes in the water table position, the model generally underestimates methane emissions rates. This may be attributed to underestimation of the extent of anoxic conditions, either due to under representation of large aggregates and other anoxic hot spots, or due to local variations in water table position relative to the areal mean value used in the simulations. For example, increasing the modeled mean aggregate size from 2 to 3 mm while keeping everything else constant resulted in considerable increase in methane emission rates during the growing season and fall refreezing, as seen in Figure 5. In other words, a better account of the extent of anoxic conditions (water table depth and aggregate/hot spot sizes) is critical for properly capturing the magnitude and dynamics of methane emissions (particularly during the growing season).

The discrepancies between model predictions and observation during the fall are likely due to model oversimplification of unsaturated soil freezing processes and interactions with diffusion and other methane transport pathways. Most models (including this one) represent the stochastic process of methane bubble ebullition very crudely; this limitation becomes accentuated during soil-freezing when various processes such as changes in liquid water content, temperature, hydrostatic pressure, and bubble transport pathways interact to affect bubble formation and ebullition. The hope is that better informed mechanistic models could improve measurement interpretation and guide future efforts for data gathering with different scale of measurements (i.e., chamber method versus eddy covariance).