A process-based model of methane consumption by upland soils

This study combines a literature survey and field observation data in an ad initio attempt to construct a process-based model of methane sink in upland soils including both the biological and physical aspects of the process. Comparison is drawn between the predicted sink rates and chamber measurements in several forest and grassland sites in the southern part of West Siberia. CH4 flux, total respiration, air and soil temperature, soil moisture, pH, organic content, bulk density and solid phase density were measured during a field campaign in summer 2014. Two datasets from literature were also used for model validation. The modeled sink rates were found to be in relatively good correspondence with the values obtained in the field. Introduction of the rhizospheric methanotrophy significantly improves the match between the model and the observations. The Q10 values of methane sink observed in the field were 1.2–1.4, which is in good agreement with the experimental results from the other studies. Based on modeling results, we also conclude that soil oxygen concentration is not a limiting factor for methane sink in upland forest and grassland ecosystems.


Introduction
The field of greenhouse gas exchange has been coming into prominence since the 1960s, as the scientific community faced the need to predict climate change that is tightly linked with the evolution of the Earth's atmosphere (Solomon 2007). However, atmospheric greenhouse gas monitoring yielded information only on net planetary-scale fluxes. As a consequence, in the 1980s the realization came that reliable long-term climate projections are impossible without the knowledge of the distribution and changes in the greenhouse gases surface sources and sinks. It was mainly the necessity to estimate these changes for the need of long-term planning of human activities that has sparked high interest in the quantification of gas exchange in natural ecosystems, particularly in soils (Heimann 2011, Pachauri et al 2014. Methane (CH 4 ) is a potent greenhouse gas, and the data on net CH 4 fluxes is important for the understanding of the climate system. It strongly influences the photochemistry of the atmosphere (Ramanathan et al 1987, Cao et al 1995. In the recent decades, the attention to methane budgets has been growing, as it was found that the radiative forcing of the atmospheric methane is second only to CO 2 (Myhre et al 2013).
The increase in atmospheric CH 4 is caused by an excess of sources over sinks, amounting on average to 5-48 TgCH 4 year −1 (1 Tg=10 12 g) in recent decades (Khalil and Rasmussen 1985, Cicerone and Oremland 1988, Dlugokencky et al 2003, Dlugokencky et al 2009. In terms of their effect on the atmospheric CH 4 budget, soils can act in two different ways: (1) flooded soils can act as sources of methane which is produced under strictly anoxic conditions (Arah and Stephen 1998); (2) upland soils can act as sinks of tropospheric CH 4 and as a biofilter for microbially produced CH 4 in anoxic soil or sediment compartments, which reduces the CH 4 emission into the atmosphere (Bender and Conrad 1994). Therefore, assessment of CH 4 fluxes from/into soils should play an important role in the prediction of trends in atmospheric CH 4 concentration.
One possible way to estimate and predict CH 4 fluxes at the regional level is to develop a processbased model to simulate CH 4 flux rates in various ecosystems. It should implement an assessment of CH 4 fluxes by linking climatic, edaphic and biological controls and provide a mechanistic basis for spatial analysis and for future change projections at the regional and global scales (Cao et al 1995).
Wetlands Methane is largely removed from the atmosphere by tropospheric oxidation by hydroxyl radical (OH), which accounts for 85%-90% of the estimated annual mean sink of 570 TgCH 4 year −1 (Bousquet et al 2006). The remaining sink is thought to be split in roughly equal parts between the stratospheric removal by OH and O 1 (D) and biological consumption in surface soils (Smith et al 2000, Curry 2007. Terrestrial environments are the only known net biological sinks for atmospheric methane. Soils are considered an important component of global methane dynamics (Adamsen and King 1993), consuming 9 to 45 Tg CH 4 year −1 (Ehhalt et al 2001, Curry 2007, an amount comparable to 1 to 7% of the total global emission. However, statistical upscaling from the distribution of actual measurements leads to a much wider range of uncertainty, 7-120 TgCH 4 year −1 (Smith et al 2000, Curry 2007. Currently, the level of understanding of the soil sink of atmospheric methane is inferior to that of atmospheric sink, and few attempts to model it have been made (Striegl 1993, Potter et al 1996, Curry 2007. Net methane emission is a sum of production and oxidation. The latter might be rather high, e.g. methane oxidation in the surface was shown to reduce methane emissions from saturated wetland soils by 10%-90% (King et al 1990, Epp and Chanton 1993, Schipper and Reddy 1996. Therefore, any adequate CH 4 emission model must contain a module describing methane consumption. Besides, on the regional scale, the areas that emit methane of geological origin, e.g. the 'mud volcanoes', (Etiope 2009, Etiope et al 2009, are typically overlaid by efficiently oxidizing ecosystems as forests and grasslands (Belova et al 2013, Oshkin et al 2014). Thus, a methane oxidation module is a requirement for any regional-scale model as well.
Methane consumption in upland soils varies strongly on local, global, seasonal and interannual scales (see table 1). To explain this variability, such wellknown CH 4 sink controls as ground surface temperature, water table depth, above-ground biomass (Arah andStephen 1998, Cao et al 1998), soil properties (e.g., bulk density, porosity) (Striegl 1993, Ridgwell et al 1999, Curry 2007) are nowadays accounted for in field measurements and laboratory experiments. The emerging relationships are attractive to global climate modelers, but they mask a lot of potentially important detail. A process-based model should indicate what underlies the correlations obtained by measurement, and under what circumstances they are susceptible to alteration (Arah and Stephen 1998). Last but not least, the ability to build a model founded on basic principles is the best test of our understanding of the process.
Over the recent years, estimation of methane emission from the Russian territory has been our general goal (Glagolev et al 2011, Glagolev et al 2012, Sabrekov et al 2013, Sabrekov et al 2014. As mentioned earlier, such an estimate cannot be deduced without having a methane sink model for upland soils. Existing methane sink models are largely empirical, particularly in regard to their treatment of biological oxidation, with rare exceptions (Grant 1998, Saggar et al 2007, Zhuang et al 2013. However, even those models do not account for certain specific features of methane sink in soils, such as methane consumption by microorganisms living on plant roots. Therefore, the model in development had to satisfy the two main requirements. First, the model must be a process-based so that it can well reproduce the process of methane sink by the land biomes based on the known biochemical and physical processes. Second, it must contain only those parameters that can be obtained for all types of soils and biomes over the Russian territory. Due to these requirements, we had to only use the average parameter values found for the respective biome types and soils in literature. This study presents an ad initio attempt to construct a process-based model of methane sink in upland soils including both the biological and physical aspects of this problem without any calibration of model parameters. Since we do not consider seasonally or permanently waterlogged soils, methane production is assumed to be negligible.

Test sites
The field experiments were carried out during the 2014 summer period at 3 sites in the south taiga zone of Western Siberia, in one forest and two grassland sites. The forest site (FS) is a coniferous spruce-pine-fir forest, the grasslands sites include a mesophilic grassland (G1) and a mesophilic grassland with sparse birch cover (G2). Detailed plant community descriptions are presented in appendix A. Basic soil properties of the test sites are shown in table 2.
The water table depth was located at a 3-4 m depth. Soil water content did not exceed 70% of the maximum water holding capacity. Methane concentration measurements at different depths revealed that, in such conditions, methane production by soils is negligible (Whalen et al 1992, Adamsen and King 1993, Priemé and Christensen 1997. Therefore, methane production is assumed to be negligible in our model.

Study methods
2.2.1. CH 4 and CO 2 flux measurements The measurements of CH 4 flux were performed using the static chamber method (Hutchinson and Mosier 1981) following the methodology described earlier by Glagolev et al (2011). The chamber consisted of two parts: a permanent square stainless steel collar (40 cm×40 cm, embedded 15 cm into the soil surface), and a removable plexiglass box (30 or 40 cm height). To minimize the changes of chamber temperature during measurement, the plexiglass box was covered with reflecting aluminum fabric. The air inside the chamber was circulated by a battery-operated internal fan; a water channel on the chamber rim acted as a lock against leaks into or out of the chamber. Four gas samples were taken at 20 min intervals into 12 ml nylon syringes (SFM, Germany). The total chamber closure time was 60 min After sampling, the syringes were immediately sealed with rubber stoppers and delivered to the laboratory. The CH 4 concentrations were corrected for leakages as described in (Glagolev et al 2011). Heather grassland 52.00°N, 5.78°E −0.028±0.020 b July −0.030±0.015 b August a SD is usually given for temporal replicates. b SD is given for spatial replicates. Until the chromatographic analysis, the syringes with the samples had been kept in salt solution to prevent methane leakage. Boiled water was used for this purpose, because it does not contain methane in the amounts capable of affecting the measurement. Methane concentrations were measured with a modified gas chromatograph CPM-4 ('Chromatograph', Moscow, Russia) having a flame ionization detector of a chromatograph LHM-80 ('Chromatograph', Moscow, Russia), 1 m stainless steel column (2.5 mm o.d.) filled with Sovpol (80-100 mesh) at 40°C with hydrogen as a carrier gas (flow rate 5 ml min −1 ). The loop volume was 0.5 ml, the volume of injected sample was 3-4 ml. Each sample of gas from a syringe was analyzed three times; the mean of the three concentrations was used for the flux calculation. The gas chromatograph was calibrated with standard gases (1.99±0.01, 5.00±0.01 and 9.84±0.01 ppmv methane in a synthetic air) prepared at the National Institute for Environmental Studies, Japan. The R 2 values for the linear correlation between signal (area of peak) and concentration in the standard were 0.998 and higher. The error of concentration measurement (standard deviation as percent of the mean of 6 to 10 daily repetitions of the standard) was typically 0.5% for the 1.99 ppmv CH 4 standard. Units were converted from ppmv to mgCH 4 m −3 using the ideal gas law.
Carbon dioxide fluxes were measured in the same way, except that in that case, four samples were taken at 3 min intervals over a period of 9 min. The aboveground biomass was not removed within the chamber collars. CO 2 concentration in the samples was measured not later than a few hours after sampling using an infrared gas analyzer DX-6100 (RMT Ltd, Russia). This device was calibrated with a standard gas (357±5 and 708±10 ppmv CO 2 in a synthetic air), prepared at VNIIEM Corporation, Russia. The relative error of concentration measurement (standard deviation as percent of the mean of 3 to 5 daily repetitions of standard) was typically 3% for both standards.
Fluxes were calculated from the linear regression (Kahaner et al 1989) for CO 2 emission and exponential regression for CH 4 uptake (see appendix B), with weights inverse to concentration measurement uncertainty for the chamber headspace concentration versus measurement time. The minimal detectable fluxes (corresponding to a chamber headspace concentration change on magnitude of concentration measurement error over chamber closure time (Wang and Wang 2003)) were 0.003 gC-CO 2 m −2 h −1 for CO 2 emission and −0.004 mgCH 4 m −2 h −1 for methane uptake. Throughout the manuscript, the convention that fluxes to the atmosphere are positive is adopted. Since the NDFE-corrected CH 4 and CO 2 fluxes (see (Livingston et al 2010) for details) were just 2% and 6% higher in terms of absolute values, respectively, than the fluxes calculated with the above method, the uncorrected values were used instead.

Environmental characteristics
Air and soil temperatures were measured during flux measurement by the temperature loggers TERMO-CHRON iButton DS 1921-1922. The frequency of measurements was once per minute, the accuracy of individual measurement being 0.125°C. The obtained temperature data were averaged across all replicates for each chamber site. After the flux measurement at each site, soil samples were taken from three depths (5, 10 and 15 cm), three replicates each, in order to measure their water content. Soil water content was measured gravimetrically by oven-drying at 105°C. Soil samples for physical analyses 0.5 kg each were randomly picked from each soil horizon. Soil physical properties were determined as described in (Shein 2015). Bulk density was measured by drying a known volume of field soil sample until its weight stabilized, whereafter solid phase density was determined by displacement of water by the known mass of soil. The soil clay content was determined by the pipette method following pretreatment to remove soluble salts, organic matter and iron oxides. Soil organic content was determined by CHNS-analyzer PE-2400 (PerkinElmer, USA). Soil pH measured in a 1:2 ratio (soil:distilled water) using a glass electrode.
A WRB soil classification scheme was used to classify the soils (IUSS 2014). Botanical descriptions of the vegetation communities within each chamber site were conducted.

Data analysis
For the comparison of modeled against measured fluxes, measurements on eight different chamber sites (six on FS site, one on G1 and one on G2) were provided. The total number of CH 4 flux measurements is 40, CO 2 flux-38. For each chamber site (in the same point in space) 2 to 12 temporal replicates of methane flux and 2 to 10 temporal replicates of total ecosystem respiration were taken in a row (i.e., within several hours) were obtained. For further calculations and comparisons, we use the weighted median of methane flux and median of total ecosystem respiration across all replicates for each chamber site. Weights were assigned in inverse proportion to squares of individual flux uncertainty. Weighted median was calculated as described by Cormen (2009). The solution of partial differential equations, numerical integration for calculation of root biomass were performed with MATLAB v. 7.8.0 (MathWorks, USA). Uncertainty of model predictions and sensitivity of model to uncertainty in certain parameters were calculated using bootstrapping as described in appendix C.

Model validation against data from other studies
For the model to be fully validated, datasets originating from other sites were also used. Unfortunately, the model runs require rather large datasets, so only two other datasets were deemed fit for comparison with the model values and the observed sink rates.
The first of the comparison datasets was obtained in the mixed hardwood forest site in College Woods (43.13°N, 71.95°W, New Hampshire, USA) during the growing seasons of 1989-1993. However, while only the 1989-1990 data have been published (Crill 1991), the 1991-1993 data were obtained using the same methods at the same sites. In order to achieve higher precision in the live root biomass and in comparison of modeled versus observed methane sink magnitudes, soil respiration, methane consumption by soil methanotrophs (bothper gram of dry soil) and their temperature sensitivities were adopted from the incubation experiment data obtained at the same site (Crill 1991).
The seasonal soil moisture profile was derived from the averaged multiannual data (2001)(2002)(2003)(2004)(2005)(2006)(2007) from the Hemlock tower site in Harvard Forest, which is similar to the College Wood site in terms of climatic characteristics and vegetation cover (Harvard Forest Data Archive 2009). The data on soil clay fraction was adopted from the field sampling results for the same region (Finzi 2009). The rest of the required data were adopted from the original source (Crill 1991). In order to reduce the importance of individual outliers and measurement bias, monthly averages for April-October were used (30-37 individual CH 4 and CO 2 flux measurements per month). Since there was virtually no vegetation inside the chamber collars during the measurement (i.e. the entire root biomass was contributed by the nearby trees), above-ground biomass was not considered in the calculation of root biomass.
The second dataset was obtained in the native tallgrass prairie and experimental agricultural field sites in The Konza Prairie Research Natural Area (39.08°N, 95.58°E, Kansas, USA) over the growing season of 1990 (Tate and Striegl 1993). For more precise calculation of the live root biomass, the soil respiration magnitudes (per gram of dry soil) and its temperature sensitivity were taken from the The Konza Prairie Research Natural Area soil incubations (Fierer et al 2006). The soil clay content (Boutton et al 1998, Nippert et al 2012), root-to-shoot ratio (Ojima et al 1994, Fay et al 2003, Nippert et al 2012 and soil bulk density and solid phase density (Grahammer et al 1991, Shaver et al 2002 of the Konza Prairie Research Natural Area soil were used (as averages). Where unavailable for the same area, the soil data from the same region were used.

Model description
A number experimental studies show that the following factors are of importance for methane consumption: • soil moisture-(Adamsen and King 1993, Castro et al 1995, Gulledge and Schimel 1998; • diffusivity in the soil pore space-(Dörr et al 1993, Potter et al 1996, Ball et al 1997; • methanotrophy substrate concentrations (oxygen and methane) and the capacity of methanotrophs of various soils to consume them-(Bender and Conrad 1992, 1994, Knief et al 2003, Knief and Dunfield 2005. The model is designed to couple the processes of consumption and transport of gaseous oxygen and methane in pore space of one-dimensional column of upland soils. Influence of the factors listed above was taken into account. Since the soil at the site is not waterlogged (soil moisture is less than 70% of maximum water holding capacity), CH 4 production is assumed to be negligible (Crill 1991, Whalen et al 1992, Adamsen and King 1993, Priemé and Christensen 1997. The model assumes that methane is consumed by two groups of methanotrophs: those living on plant roots and those inhabiting the soil, but not associated with the rhizosphere (termed 'rhizospheric' and 'soil' methanotrophs from now on, correspondingly). The CH 4 consumption rate by both rhizospheric and soil methanotrophs follows Michaelis-Menten kinetics for both methane and oxygen, and is also a function of soil temperature and moisture.
The model describes respiration of both the plant roots and the microorganisms inhabiting the soil. Soil respiration rate is the function of soil temperature and soil carbon content. Root respiration rate is the function of soil temperature and root biomass. Both soil and root respiration follows Michaelis-Menten kinetics for oxygen. Transport of both CH 4 and O 2 in soil is by molecular diffusion through the air-filled soil pore space. The model calculates methane fluxes to the atmosphere. Detailed description of the model is given in appendix C (table C1). Input parameters include air temperature, CO 2 flux measured by dark chambers (total respiration, TR), soil profiles of temperature, moisture, bulk density, solid phase density, carbon content and clay content.
The model is formulated similarly to the other modern models or model blocks predicting methane consumption in wetlands (Walter et al 1996, Arah and Stephen 1998, Grant 1998) and upland soils (Zhuang et al 2013). However, there are notable differences between those prototype models. First, all the necessary parameters were obtained from literature for appropriate climate zone (if it was possible) and averaged across all sources. The model parameters were not calibrated, as the aim was to examine how modern knowledge of methane consumption in upland soils can reproduce the values of methane fluxes observed in chamber measurements. Second, in the above-mentioned models, rhizospheric methane consumption was not considered. Current model introduces that process, taking into account root biomass and root density distribution in the soil profile (see appendix C for details). In order to estimate root biomass, the balance approach was used. Since total respiration is a sum of soil, below-ground and above-ground plant biomass respiration, root biomass can be estimated from their difference. Soil, root and shoot respiration rates per unit soil/plant mass and rootto-shoot ratio required for these calculations were borrowed from literature (see appendix C for details).

Results and discussion
The values of the weighted medians (WM) of methane flux, medians of TR and the mean magnitudes of various ecological parameters are presented in table 3 (see appendix E (table E1) for the full overview of results).
The results of model runs are presented in figure 1. The simulated methane uptake was generally in good agreement with the chamber flux data, although underestimation did occur in two of the chamber sites (2 and 8 from table 3). Statistical analysis of the differences in methane oxidation between the sites is presented in appendix F (table F1). It doesn't show a lot of significant differences. However, according to literature data reported in table 1 it appears to be a typical pattern. Methane fluxes obtained in our investigation are in the same range and have the same level of uncertainty. It means that for analysis of predictive ability of process-based models more precise data on methane consumption are necessary. Of course, the agreement between modeled and experimental values could be worse for individual sites if model parameters adopted only from literature are used. However, at the regional and global scale this approach could be functional because of the averaging of model parameters from several sources. Figure 1 shows standard deviations of methane fluxes simulated by the model (vertical error bars). The uncertainty follows from the precision of the model parameters, and on average is 0.04 mgCH 4 m −2 h −1 . Since all model parameters were adopted from literature data, it is possible to identify parameters crucial for making methane uptake modeling more accurate. It was done by bootstrapping as described in appendix C. Since several model parameters are dependent and/or were measured simultaneously, they were grouped in categories. This analysis is presented in table 4. It shows that, on average, uncertainties introduced by soil and rhizospheric methanotrophy components are close and are higher than uncertainty from the other groups of parameters. While there are no studies about variability of rhizospheric methanotrophy in upland soils, several studies indicate that parameters of free soil methanotrophy are very variable between different soil and ecosystems types (Czepiel et al 1995, Knief et al 2003, Nazaries et al 2013b. For improving the soil methane consumption models it is necessary to reveal controls of this variability.

Model sensitivity
A model simulation was conducted to determine the effect of soil temperature on methane fluxes. Numerical tests have shown that a 10°C increase in soil temperature leads to a 18%-40% increase in methane sink (in other words, the flux becomes more negative). This result is in good agreement with the data from experimental studies on temperature sensitivity of CH 4 sink in soils and confirm the idea that during growing season methane consumption is mostly limited not by temperature but by diffusivity of CH 4 in soil pore space (Born et al 1990, King 1997, Bowden et al 1998, Gulledge and Schimel 2000. It is also noteworthy that the temperature increase produced a greater effect when the soil moisture was lower, just as expected (Priemé and Christensen 1997).
Modeled oxygen concentrations were high enough in the soil even at the highest possible soil moisture contents (i.e. the lowest diffusivity), implying that oxygen did not limit methanotrophy. It is confirmed by the numerical experiments: a twofold increase of both root and soil respiration reduced the CH 4 flux by less than 0.1%.

The role of rhizospheric methanotrophy
The model describes two methane sink components: consumption by the soil methanotrophs and consumption by the rhizospheric methanotrophs. While the former is well known (e.g., Holmes et al 1999, Kolb et al 2003, Kolb et al 2005, the latter is much more obscure, having been mostly studied in plant communities of wet ecosystems (Gerard andChanton 1993, King 1994). We conducted a quantitative examination of the importance of rhizospheric methane consumption components (table 5, model experiments 1-5). It appeared that the assumption of zero activity of the rhizospheric methanotrophs significantly reduces the match between the model and the observations (table 5, model experiment 1 and 2). Thus without rhizospheric methanotrophy it is impossible to explain local spatial variability of methane flux into the soil.
As we use the parameters of soil methane oxidation taken from the literature data, they would not, most likely, match those typical for the investigated soils. However, even if those values are in fact higher or lower, model experiments show that it would not result in a much better R 2 for the observed versus predicted methane fluxes with excluded rhizospheric methanotrophy (table 5, model experiment 3). Neither do the values of rhizospheric methane oxidation parameters affect the correspondence between measured and modeled fluxes. Because they were obtained for wet ecosystems, they could overestimate rhizospheric methane consumption in upland soils. Model calculations reveal that R 2 does not become much worse if rate of rhizospheric methane oxidation becomes 2 or 4 times smaller (table 5, model experiment 4). Of course, the relative contributions of the rhizospheric and soil methanotrophs to methane consumption cannot be reliably partitioned based solely on the present data. But numerical experiments show that improvement effect due to taking rhizospheric methanotrophy into account does not depend on values of microbiological parameters used in the model. , 4 B r ) 0.0277 Sensitivity of methanotrophy to soil moisture (a 0 , a 1 , a 2 ) 0.0105 Sensitivity of methanotrophy to soil temperature All other parameters 0.0083 It is also interesting to test how methane consumption would change if rhizospheric methanotrophy is not influenced by soil moisture. Microorganisms on the root surface usually live in optimal moisture conditions thanks to the presence of root exudates and specific microclimate and do not depend on moisture of surrounding soil (Philippot et al 2009). Numerical experiment shows that this assumption leads to better R 2 for the observed Versus predicted methane fluxes (table 5, model experiment 5).
The data presented above reflect the within-ecosystem local spatial variability of methane sink in natural ecosystems. Modeling results suggest that root biomass might be the factor explaining a substantial fraction of this variability. However, the spatial variability of the sink might also be explained by some other factors. As long as the CH 4 consumption rate by soil free methanotrophs used in the model was the same for all chamber sites, it could not explain this variability in the model output.
The data from Tate and Striegl (1993) allow the performance of simultaneous model validation for both spatial and temporal methane flux variability. In that study, CH 4 fluxes were measured in tallgrass prairie ecosystem sites during the growing season of 1990. Calculations were done only for May-July period because since August the respiration was suppressed by drought so that a reliable root biomass estimate could not be obtained. The effect of rhizospheric methanotrophy inclusion is shown in figure 2 and table 5 (model experiments 6-8). Without rhizospheric methanotrophy, the model could not satisfactorily explain CH 4 uptake variability (table 5, model experiment 7). Conversely, the inclusion of rhizospheric methanotrophy leads to the significance of the predicted versus observed flux linear dependence (table 5, model experiments 6, 8).
The role of rhizospheric methanotrophy can be clearly demonstrated by the example of agricultural ecosystems. Having been planted, the cultivated plants intensively grow their root systems. Comparison with the model may be drawn based on the sorghum field data of the first month after planting from (Tate and Striegl 1993). Calculations performed both with and without rhizospheric methanotrophy are presented in figure 3 and table 5 (model experiments 9-11). Again, the inclusion of rhizospheric methanotrophy grants a statistical significance to the linear regression of predicted versus observed fluxes. In contrast, the same correlation would be negative and not significant if rhizospheric methanotrophy is omitted (table 5, model experiment 10). Use of multiannual (1989)(1990)(1991)(1992)(1993) seasonal dynamic of methane sink in the College Woods (Crill 1991 and  unpublished data by P Crill for the same site) allows to avoid influence of spatial uncertainty on model validation. Since measurements were conducted seasonally in the same locations spatial variability could not affect the CH 4 consumption rates. The root biomass has pronounced seasonal dynamic (Vogt et al 1995, Pregitzer et al 2000 and can influence on methane consumption on a seasonal scale. Since the profile data on methane consumption rate by soil free methanotrophs are available for this site (Crill 1991), it is also possible to provide a quantitative estimation of the role of rhizospheric methanotrophy.
The results of model calculations are presented in figure 4 and table 5 (model experiments 12-14). Free soil methanotrophy cannot explain the methane flux from the atmosphere to the soil on its own. Predicted fluxes without rhizospheric methanotrophy are out of 99% confidence interval for average observed flux for all months. Only the inclusion of rhizospheric methanotrophy lets one achieve quantitative correspondence between predicted and observed fluxes.

Perspectives of the model development
Several improvements may be introduced in the model. First, several more factors affecting the methanotrophy activity may be considered: In the current model version, inhibition could not be simulated because data on soil inhibitor concentrations were not available. The soil pH effect on methanotrophy is ambiguous. According to the results presented elsewhere, methanotrophs can live within a wide pH range (Serrano-Silva et al 2014). To our knowledge, no data showing independent pH effect on upland soil methanotrophy were obtained. One may surmise that, by virtue of great species diversity and adaptation capacity, methanotrophs may form consortia that efficiently oxidize soil methane at any pH level possible in upland soils.
Another possibility of improvement lies in the incorporation of microbial community in the model. No current models account for the role of various microbial communities in soil CH 4 processes. In part, this is because, for all the extreme diversity of the soil microbial communities, they also show high spatial and temporal variability. Few microbes are cultivable under current laboratory regimes and their physiological capabilities therefore remain unknown. The above issues make the parameterization of microbial data challenging (Nazaries et al 2013a). Nowadays, all the information on microbial communities is implicitly contained in the other model parameters. Besides, a number of methanotrophy-related effects in soils can only be explained using microbiological data. For example, the decreasing of maximal soil free methanotrophs CH 4 oxidation rate with depth (Crill 1991, Czepiel et al 1995) can be associated with the reduction of the methanotrophs biomass (Bender and Conrad 1992). This can be related to the downward reduction of methane concentrations in the soil profile, i.e. lower substrate amount leads to lower methanotroph biomass.
The differences in the maximal rate of CH 4 oxidation by soil free methanotrophs between different soil types may be related to the variability in the methano-  providing the lack of experimental data, the best approach would be to search for significant correlations between the methanotroph communuity structure and the climatic, soil and vegetation community features. Here, one should eliminate the effects of the other factors (temperature, soil moisture and root biomass). If significant dependencies are found, they can be used to forecast the response of the CH 4 consumption on global climatic changes.

Conclusions
The factors controlling CH 4 consumption show effects on different temporal and spatial scales and confound each other. The necessity to take into account the combined effects of different interacting controls on soil methane consumption motivated the formulation of the new process-based model in the current study.
Both the data obtained by the authors and those adopted from literature equivocally indicate that the inclusion of rhizospheric methanotrophy significantly improves correlation between observed and predicted methane fluxes. Numerical experiments show that this improvement does not depend on the values of microbiological parameters used in the model. Comparison with the Crill (1991) data showed that, without root methanotrophy, one cannot achieve quantitative correspondence between the observed and predicted fluxes. Important limitations in model validation originated from the root biomass not being estimated from field data and the maximal rate of CH 4 oxidation by soil free and rhizospheric methanotrophs not being measured directly in the studied ecosystems. Therefore, one cannot state that exactly the rhizospheric methanotrophy explains the observed variability in the CH 4 uptake. Nevertheless, it should be noted that the process explaining this variability is related to the other, ecologically different, group of methanotrophs.
The key difference is that activity of this methanotrophs group correlates with the plant root biomass, which is not the case with the soil free methanotrophs.
The most obvious explanation here is the activity of rhizospheric methanotrophs.
The revealed relationship may have important consequences in the future research into upland soil methane sink. First, collars for methane flux chamber measurements are often placed in the areas completely devoid of or sparsely covered with vegetation, to simplify their installation. That might lead to underestimation of consumption. Second, global climatic changes might show effects on the methane sink in soils also via vegetation community changes. Third, the rhizospheric methanotrophs live in different ecological conditions compared with the free soil methanotrophs. For that reason, their response to the change in soil properties would be different, too. Fourth, the methanotrophic consortium associated with the plant roots would differ from that of the soil free methanotrophs in terms of species composition. Consequently, that consortium might have different ecological characteristics. The latter is important for the perspective inclusion of microbial data in models for improved prediction of CH 4 flux from terrestrial ecosystems.

Acknowledgments
We thank all participants of the 2014 summer field expedition. This study SSM acknowledges support from the MEXT Arctic GRENE (ID 5) program. We also want to thank P Crill and R Varner for giving access to College Woods site data. AFS thanks M Semenov, A Zakharenko and T Yorke for their contributions.

Appendix A. Test site descriptions
The forest site (FS) is a the coniferous spruce-pine-fir forest with Picea abies, Pinus sibirica and Abies sibirica in woody layer, Prunus padus, Caragana arborescens and Rubus idaeus in a shrub layer and Equisetum sylvaticum, Figure A1. Map of test sites on a Landsat satellite image. Aegopodium podagraria, Oxalis acetosella, Agrimonia pilosa, Urtica dioica, Glechoma hederacea, Vicia sylvatica, Aconitum septentrionale, Maianthemum bifolium and other species in an understory. Site G1 is mesophilic grassland with Cirsium vulgare, Galium spurium, Phleum pratense, Stellaria media, Ranunculus acris, Achillea millefolium as dominant species. Site G2 is a mesophilic grassland with a sparse birch (Betula pendula) cover, Caragana arborescens and Rubus idaeus in the shrub layer and Chamerion angustifolium, Filipendula ulmaria, Cirsium arvense, Rumex confertus, Dactylis glomerata, Phleum pratense and others comprising the grass layer. All sites are shown on figure A1.

Appendix B. Exponential regression for methane consumption
When measuring the rate of soil-to-atmosphere gas emission, it is usually assumed that gas concentration in the chamber increases linearly (Boeckx et al 1996, Alm et al 1999, Altor and Mitsch 2008, Sabrekov et al 2014. However, the situation changes dramatically when gas sink in soil is concerned. At high oxygen concentration, which is true within the chamber during the measurement, the rate of methane oxidation (R oxid , mg m −3 h −1 ) in the chamber headspace can be described by Michaelis-Menten-type equation: where k=V max /K m . Thus, first-order kinetics is obtained, which is characterized by exponential decrease of concentration with time (t): is the initial concentration in the chamber headspace. King and Adamsen (1992) and King (1994) provided experimental evidence of methane oxidation being described by first-order kinetics in pure cultures of methanotroph Methylomoas rubra and in methanotroph associations on roots and rhizomes of aquatic vegetation. This was true even when methane oxidation equalled 56 times the mean atmospheric level. Boeckx et al (1996) showed the same for soil, albeit at only 6 times the atmospheric level. Chan and Parkin (2000) state that the first-order kinetics is still observed at methane concentrations exceeding 7000 times the mean atmosperic level; however, it is not readily seen in the experimental results they provide. It was shown for Methylomoas rubra that the kinetics approach zero-order type already at the 5600-fold excess. Moreover, King (1994) showed that in associations of methanotrophs on roots and rhizomes of aquatic vegetation, zero-order kinetics is obserrvable at as low as 560-fold excess. Firstorder kinetics in methane oxidation is also postulated in Jensen andOlsen (1998), Curry (2007) and many other studies. Therefore, as the exponential law of concentartion decrease in a chamber should be used, regardless of the time interval length, non-linear equation parameters С 0 and k from (B3) should be determined. However, determination of the (B3) parameters does not directly lead to the surface flux density magnitude. The net methane flux in soil ( F, mg m −2 h −1 ) may be found as is the chamber volume, S (m 2 ) the chamber base area (or collar area), C A (mg m −3 ) atmospheric methane concentration. C A is substituted into the linear approximation of R oxid , as microbial consumption of methane in real field measurements occurs at this value of methane concentration. In order to find the parameters in (B3), a logarithmic transformation can be performed to obtain After a substitution у=ln(С) the parameters of (B5) can again be found using the linear least squares method (Ryan 1997). In (B5) the sum of squared residuals for the transformed values ln(С) and the calculated values ln(С о )−k·t is minimized. Note that it is not the sum of squared residuals for the measured values С and the corresponding calculated ones. Thus the resulting parameter values do not satisfy the least squares principle (Vernin and Chanon 1986) and may serve as only the first approximation.
Improved parameter estimates can be obtained by introducing the weights (w i ). Substituting it appears reasonable to introduce the weights as is chamber headspace concentrations measured in moments t i . Using (B6-B7), weights for (B5) can be found as However, even as the weights are introduced, there is no guarantee that the weighted linear regression for (B5) would yield correct parameter values of the original nonlinear equation (B3). In order to avoid this Coefficients of the dependency of maximal methane oxidation rate on soil temperature (T s ): exp(b 2 T s 2 +b 1 T s +b 0 )/b max −3.695±0.463 Glagolev (2004Glagolev ( , 2006; b max was chosen so that the function ranges between 0 and 1   uncertainty we use a nonlinear regression (function nlinfit in MATLAB) to find parameter values in (B3).

Appendix C. Detailed model description
The model is designed to couple the processes of consumption and transport of oxygen and methane in the gaseous phase in the pore space of aerated soils. The parameters of the below relationships are presented in table C1. They were adopted from previous studies describing identical or similar ecosystems or soil types (see the references in table C1). If alternative values were reported in several studies, their average value±standard deviation was used. It is assumed that the soil pore space does not become totally saturated with water, and the water films are not interconnected (probably, with the exception of the snowmelt period which is not included in this study). Thus, we do not consider gas transport in the liquid phase, as it is negligibly small compared with gas-phase transport in the conditions described above (Moldrup et al 2000). Methane was allowed to be consumed by two groups of methanotrophs: those living on plant roots and those inhabiting the soil, but not associated with the rhizosphere (termed 'rhizospheric' and 'soil' methanotrophs from now on, correspondingly).
We formulate the general equation describing the transport and sink of methane in gaseous phase as: (g m −3 h −1 ) the rates of methane consumption by root-associated and soil free methanotrophs respectively, z (m) the spatial coordinate (positive downward) and t (h) the time. The equation for oxygen (C2) is fully analogous: diffusion in gaseous phase is only considered, and uptake is allowed by both the plant roots and the microorganisms inhabiting the soil.
(g m −3 h −1 ) the rates of root and soil respiration, respectively. For the top boundary (soil surface, 0 m), Dirichlet-type boundary condition was assumed for both gases: Diffusion is the only gas transport mechanism that is important in aerated soils (Striegl 1993), which is reflected in (3) and (4): (m 2 h −1 ) are the diffusivities for methane and oxygen respectively. Those are traditionally obtained by multiplying the gaseous tracer diffusivity in air (with the air temperature dependency) by the reducing coefficient describing the soil properties. Many ways to calculate this coefficient exist, depending on different soil properties (Moldrup et al 2000). However, more detailed coefficient calculation schemes would require an amount of field measurement data that are usually unavailable. Therefore, calculation of diffusivity according to Moldrup et al (2003) was used as in (C5) and (C6): where d 1 and d 2 (unitless) are empirical coefficients presented in table C1. Air filled porosity and total porosity were calculated according to Shein (2015) from the water density d w (kg m −3 ), experimental data on soil density ( )  (8) and (9): The effect of oxygen and methane concentrations on the rate of methane consumption was described with a Michaelis-Menten function. Methane consumption rate by the rhizospheric methanotrophs was calculated as in (C10): where V r,CH 4 (g CH 4 kg −1 root dry matter h −1 ) is the maximal rate of CH 4 oxidation by rhizospheric methanotrophs, ( ( )) f T z s (unitless) the temperature dependency function of methane oxidation varying from 0 to 1, ( ( )) f m z s (unitless) the soil moisture dependency function of methane oxidation varying from 0 to 1, RD(z) (m −1 ) the root density distribution in the soil profile calculated according to Jackson et al (1996), B r (kg root dry matter m −2 ) is the total biomass of the living roots, K r,CH 4 and K ox,O 2 (g m −3 ) the Michaelis constants (the methane and oxygen concentrations at which the methane oxidation rate by rhizospheric methanotrophs is at half-maximum). The total mass of the living roots was calculated from the root-to-shoot ratio p (kg root dry matter kg −1 shoot dry matter) and the aboveground biomass B s (kg shoot dry matter m −2 ): where TR (gC-CO 2 m −2 h −1 ) is the measured total ecosystem respiration (for FS-total respiration of soil and moss-grass layer), c (gC-CO 2 g −1 O 2 ) the recalculation coefficient from gO 2 to gC-CO 2 (ratio of molar masses), i the soil horizon index, n the number of horizons, z i and z i+1 the upper and the lower boundary of a horizon, correspondingly, V , l,resp V r,resp and V s,resp (gO 2 kg −1 dry matter or dry soil h −1 ) the maximal rates of aboveground biomass, root biomass and soil respiration respectively at 20°C, Q , (unitless) the temperature dependency coefficients estimating change of aboveground biomass, root biomass and soil respiration respectively as a consequence of increasing the temperature by 10°C, T a (°C) the measured air temperature. The numerator in (C11) contains the difference between the total ecosystem respiration and the calculated respiration of soil microorganisms and roots at the measured soil temperature. The denominator represents the calculated respiration rate by unit of the aboveground biomass at the measured air temperature. The final expression for B r is: The soil moisture dependency for both components of methane consumption was taken into account as a dimensionless coefficient ranging between 0 and 1: where a 0 , a 1 , a 2 , a max (unitless) are the empirical coefficients presented in table C1 and m whc (unitless) the relation between the soil moisture m s (gH 2 O g −1 soil) and the maximal soil water holding capacity WHC (gH 2 O g −1 soil): Maximal soil water holding capacity was calculated according to Shein (2015) from the total porosity, the soil bulk density and the water density: The temperature dependence for both methane consumption groups was also derived as a dimensionless coefficient ranging between 0 and 1: ,resp 10, 20 10 ,resp 10, 20 10 . Methane oxidation by the free soil methanotrophs was taken into account in a similar way as the oxidation by the rhizospheric methanotrophs: where V s,CH 4 (gCH 4 kg −1 dry soil h −1 ) is the maximal methane consumption rate by soil methanotrophs, K s,CH 4 (gCH 4 m −3 ) the Michaelis constant (the methane concentration at which the oxidation rate by soil free methanotrophs is at half-maximum). V s,CH 4 was adopted individually (see table C1) for the litter V lit,CH 4 (gCH 4 kg −1 dry soil h −1 ) and the mineral soil layer V min,CH 4 (gCH 4 kg −1 dry soil h −1 ): where z lit (m) is the lower boundary of the litter layer and z max (m) the upper boundary where significant methane oxidation first occurs according to literature data. Due to lack or ambiguity of data on vertical profile of maximal methane sink in soil, it was described with a step-function: V min,CH 4 is constant above z max and zero below z max . The rate of root respiration was also calculated using the above expression for live root biomass and the data on root depth distribution: where K resp,O 2 (gO 2 m −3 ) is the Michaelis constant for respiration. The soil respiration rate, or, strictly speaking, the heterotrophic respiration of soil microbial community was calculated in a similar fashion: Maximal soil respiration rate at 20°C was adopted individually for the litter and mineral soil layers: lit,resp lit min,resp org lit where, V lit,resp (gO 2 kg −1 dry soil h −1 ) is the maximal litter respiration rate at 20°C, V min,resp (gO 2 kg −1 soil C h −1 ) the maximal soil respiration rate at 20°C and C org (z) (kg soil C kg −1 dry soil) the soil organic content at depth z.
The calculations were performed for the soil layer between 0 m and 1 m depths. Methane flux (mgCH 4 m −2 h −1 ) was calculated using gradient approach from the methane concentration profile in the soil (Walter et al 1996, Zhuang et al 2004: The physical and chemical properties of the soil profiles involved in the analysis are presented in appendix D (table D1). Measured values of total ecosystem respiration (for FS-total respiration of soil and moss-grass layer), soil temperature and soil moisture are shown in table 3. Bootstrap method (Efron and Tibshirani 1986) was used to estimate the uncertainty of modeled fluxes. Artificial errors were introduced into each model parameter using their standard deviations obtained from literature data (see table C1). Then, the 'noisy' fluxes were calculated using these 'noisy' parameter values. Uncertainty was estimated based on 1000 such iterations for each individual predicted flux value. Model sensitivity to uncertainty in groups of parameters was calculated in the same way: artificial errors were introduced into each model parameter of a certain group while other parameters remain constant.
Appendix D