Variability in methane emissions from West Siberia ’ s shallow boreal lakes on a regional scale and its environmental controls

Small lakes represent an important source of atmospheric CH4 from northern wetlands. However, spatio-temporal variations in flux magnitudes and the lack of knowledge about their main environmental controls contribute large uncertainty 20 into the global CH4 budget. In this study, we measured methane fluxes from small lakes using chambers and bubble traps. Field investigations were carried out in July-August 2014 within the West Siberian middle and south taiga zones. The average and median of measured methane chamber fluxes were 0.32 and 0.30 mgCH4 m -2 h -1 for middle taiga lakes and 8.6 and 4.1 mgCH4 m -2 h -1 for south taiga lakes, respectively. Pronounced flux variability was found during measurements on individual lakes, between individual lakes and between zones. To analyze these differences and the influences of 25 environmental controls we developed a new dynamic process-based model. It shows good performance with emission rates from the south taiga lakes and poor performance for individual lakes in the middle taiga region. The model shows that, besides well-known controls such as temperature, pH and lake depth, there are significant variations in the maximal methane production potential between these climatic zones. In addition, the model shows that variations of gas-filled pore space in lake sediments are capable to control the total methane emissions from individual lakes. The CH4 emissions exhibited 30 distinct zonal differences not only in absolute values but also in their probability density functions: the middle taiga lake fluxes were best described by a lognormal distribution while the south taiga lakes followed a power law distribution. The latter suggests applicability of self-organized criticality theory for methane emissions from the south taiga zone, which could help to explain the strong variability within individual lakes.


Introduction
Due to its higher global warming potential, methane contributes about 20 % of the overall greenhouse effect (Cicerone and Oremland, 1988;Wuebbles and Hayhoe, 2002; Published by Copernicus Publications on behalf of the European Geosciences Union.IPCC, 2013).Lakes and wetland ponds have strong potential impacts on the methane budget (Tranvik et al., 2009) due to their anoxic sediment conditions and often high organic matter content (Zehnder, 1978).Lake methane fluxes and their temporal patterns are still poorly constrained and form a major gap in the northern C budget (Rasilo et al., 2015).Over the past decade, new evidence has demonstrated that these systems have been underestimated in their contribution to the northern carbon balance (Kortelainen et al., 2006;Juutinen et al., 2009;Holgerson and Raymond, 2016;Wik et al., 2016b).Lakes and wetland ponds can form high CH 4 emission hotspots that contribute largely to landscape-scale CH 4 budgets but create uncertainty for bottom-up regional CH 4 emission estimates due to their small size (Bubier et al., 2005).For example, a recent global methane assessment estimates methane emissions in a wide range from 37 to 112 TgCH 4 year −1 (Saunois et al., 2016).Thus, more data are needed to resolve the divergence between top-down and bottom-up estimates, including the generation of flux measurement that is more representative in time and space (Nisbet et al., 2014;Saunois et al., 2016).
Small shallow lakes have high methane emission potential for several reasons.First, methanogenesis is sensitive to temperature conditions (Zeikus and Winfrey, 1976;Dunfield et al., 1993;Schulz et al., 1997;Yvon-Durocher et al., 2014) and shallow lakes warm up quickly during the summer season.Second, the small water volume coupled with high organic carbon content promotes the formation of anoxic hypolimnion and is related to increased concentrations and fluxes of CH 4 (Juutinen et al., 2009).Third, these lakes occupy significant areas in waterlogged regions (Downing et al., 2006) where they receive large inputs of substrate for methanogenesis in the form of terrigenous dissolved organic matter (Segers, 1998).
Methane fluxes from small lakes have great spatial variability (e.g., Casper et al., 2000;Dzyuban, 2002;Kankaala et al., 2004;Bergström et al., 2007), which hampers flux upscaling and modeling of the processes required for making regional and global estimations.This high variability results from multiple environmental controls, including biological (system productivity, organic matter loading and its mineralization, methane production, and oxidation by different groups of microorganisms), physical (temperature, mixing rate, stratification, diffusion and bubble transport rate) and chemical factors (such as concentrations of methane, oxygen, inhibitors and pH) (Rudd and Hamilton, 1978;Bastviken et al., 2004;Lofton et al., 2014;DelSontro et al., 2011DelSontro et al., , 2015DelSontro et al., , 2016)).Due to these variations, the effect of different control factors is complicated and still insufficiently known.In particular we must develop robust relationships between lake CH 4 emissions and their potential controlling factors to facilitate both prediction and spatial extrapolation (Rasilo et al., 2015).
The excess water supply and flat topography with impeded drainage provides favorable conditions for wetland and lake formation in West Siberia (Terentieva et al., 2016).A few studies of methane emission from lakes of this region (Gal'chenko et al., 2001;Repo et al., 2007;Glagolev et al., 2011;Sabrekov et al., 2013) indicate that the methane flux from middle taiga lakes is 10 times lower in magnitude than from southern taiga lakes.These findings agree with data showing strong latitude gradient of methane release from permafrost lakes in the Northern Hemisphere (Holgerson and Raymond, 2016).Despite many studies addressing local spatial variability in methane emission from lakes, there is still a gap of knowledge about methane emission from lakes and, even more rare, their environmental controls on a regional scale (West et al., 2015).An understanding of the nature of this difference is important for modeling and reliably estimating regional methane emission from lakes.In this study, we have the following objectives: to estimate methane fluxes from small lakes of the middle and southern taiga to detect key environmental controls of methane emissions on both regional (between zones) and local (within each zone) scale to improve the precision of methane emission modeling.
In order to complete these tasks, the following methodology was applied.The best way to take into account the key processes is to construct a process-based model founded on established physical, chemical and biological dependences.
Regression procedures, commonly used for identification of environmental controls, cannot fully reflect complicated interactions between different controls and thus mask a lot of important details.Rather, with process-based modeling, it is possible to test which dependences and which parameters are reliable (at least, for a certain climate zone) and which are not.After taking into account well-known dependences, a comparison of predicted and measured values can help to find new potentially important controls.Since the focus of our work is on lake methane emissions, it is necessary to examine their huge temporal variability and episodic peaks attributed to bubbling, which is challenging to quantify (Walter et al., 2006;Wik et al., 2013).It can be suggested that these stochastic emissions can be modeled using self-organized criticality (SOC) theory (for details on this theory see, e.g., Bak et al., 1988;Bak, 1996;Turcotte, 1999), which can simplify scaling the flux measurements across larger areas.Bubbling is similar to systems showing SOC behavior, in which constant external force leads to rapid changes in nonlinear interactions after reaching a certain threshold (Bak et al., 1988).Systems with SOC behavior occur in many disciplines, including physics, biology and economics (Bak et al., 1988).It has been argued that earthquakes, landslides, forest fires and species extinctions are examples of SOC in nature (Bak, 1996).To the best of our knowledge, no applications of SOC to bubbling in lakes have been previously published.However, such a behavior of gas bubbles in foam (Kawasaki and Okuzono, 1996) and in different artificial systems (see, for example, Juodis et al., 2006;Petrashenko et al., 2005) is well known.Therefore, we test whether the measured flux values demonstrate SOC behavior and examine the consequences for upscaling.
2 Materials and methods

Study sites
The studied lakes are located in two different boreal zones (both with a subarctic climate according to Köppen climate classification) of the Western Siberian Lowlands (Russian Federation) (Fig. 1).The northern study area is in the middle taiga zone (referred to hereafter as MT) about 20-30 km southwest from Khanty-Mansiysk (61 • N, 69 • E).Mean (range for 1979-2014) temperature is 18.4 (14.0 to 23.1) • C in July and −18.9 (−22.9 to −14.9) • C in January; mean annual temperature and precipitation are −0.8 (−4.7 to 3.5) • C and 530 (308-762) mm, respectively (All Russian Research Institute of Hydrometeorological Information -World Data Center, Khanty-Mansiysk).The southern study area is in the southern taiga (referred to hereafter as ST) zone about 100-200 km northwest from Tomsk (57 • N, 83 • E), approximately 900 km southeast of the MT study area.Mean (range for 1979-2014) temperature is 18.7 (13.7 to 24.8) • C in July and −17.1 (−20.9 to −13.0) • C in January; mean annual temperature and precipitation are 0.9 (−3.5 to 6.2) • C and 567 (292 to 768) mm, respectively (All Russian Research Institute of Hydrometeorological Information -World Data Center, Tomsk).The climate in both regions is continental with moderate annual rainfall, long and cold winters, and warm summers.Permafrost is absent in both study regions.
The MT lakes are mostly acidic, have low ionic strength and are surrounded by wetlands (Fig. 1, Table 1).Four lakes were selected to cover the range of sediment properties in this zone.Lake Muhrino has peaty sediments with high mineral content (sandy bedrock), lakes Babochka and Lebedinoe both have mineral-free peaty sediments, and Lake Bondarevskoe has sapropel (flocculated humic material) sediments.The ST lakes are more diverse due to high groundwater mineralization.Lakes Bakchar-ryam, Plotnikovo and the three Bakchar-forest lakes (1-3 in Fig. 1) represent mesotrophic or eutrophic lakes surrounded by soils rich in clay and grasslands.Lakes Gavrilovka-1 and Gavrilovka-2 represent lakes with low nutrient concentrations influenced by groundwaters with a high pH.Lakes Bakchar-bog-1 and Bakchar-bog-2 represent acidic humic wetland lakes with low pH, low ionic strength and low nutrient concentrations.Finally, lake Ob' Floodplain represents floodplain lakes (oxbows) with extremely high nutrient concentrations.

Field measurements
Field investigations were carried out during summer 2014.We conducted 190 methane flux and 170 carbon dioxide flux measurements, with 70 and 60 in MT and 120 and 110 in the ST, respectively.The total field measurement time varied from 4 to 10 h per lake, while the average was 6 h.All measurements were carried out between 10:00 and 20:00 LT; each lake was visited one time.All measurements were conducted using a boat to prevent any influence on the lake vegetation or sediments.CH 4 emissions were measured from the lakes using closed floating chambers.CO 2 fluxes were also measured as a point of comparison.The plexiglas chambers were equipped with plastic bottles to ensure a 5 cm floating depth and had dimensions of 40 × 40 × 30 cm (length × width × height), creating www.biogeosciences.net/14/3715/2017/Biogeosciences, 14, 3715-3742, 2017 a headspace volume of 0.048 m 3 .Chambers were covered by aluminum foil to prevent changes in temperature inside due to solar heating.Four gas samples for both CH 4 and CO 2 were taken from the fan-mixed chamber headspace at 10-15 min intervals during a chamber closure period of 30-45 min with 12 or 20 mL polypropylene syringes (sfm medical devices, Germany).Prior to sampling, chamber air was used to flush the sampling tube several times.Until the chromatographic analysis, the syringes with the CH 4 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 amounts capable of affecting the measurement.The sample CH 4 concentration was measured on a calibrated gas chromatograph Crystall-5000 (Chromatec, Yoshkar-Ola, Russia) with a flame ionization detector (FID) and column (3 m) filled by HayeSep Q (80-100 mesh) at 70 • C with nitrogen as a carrier gas (flow rate 30 mL min −1 ) or on a calibrated gas chromatograph KhPM-4 (Khromatograf Co., Moscow, Russia) with an FID and column (1 m) filled by Sovpol at 40 • C with hydrogen as a carrier gas (flow rate 10 mL min −1 ) within 72 h after sampling.Uncertainty of individual concentration measurement averaged ±0.03 ppm.
The CH 4 concentrations were corrected for leakages (as described in Glagolev et al., 2011).CO 2 concentrations were measured using an infrared gas analyzer within 4 h of sampling (DX-6100; RMT Ltd., Moscow, Russia).The uncertainty of individual CO 2 concentration measurements averaged ±2.1 ppm.Each gas sample for both methane and carbon dioxide was analyzed in triplicate.Fluxes were calculated from linear regression between chamber headspace CO 2 and CH 4 concentration versus measurement time using weights inverse to the measurement of gas concentration un-certainty (Kahaner et al., 1989).Fluxes are reported following the sign convention that exchange from the landscape to the atmosphere is positive.
During the chamber measurements, near-surface water (10 cm depth) was sampled for dissolved CH 4 .The 10 mL water samples were taken with a polypropylene syringe and vigorously shaken for 3 min with a 10 mL headspace of known CH 4 concentration (taken from ambient air in which methane concentration was 1.88 ppm).The headspace CH 4 concentrations were then analyzed in a field laboratory on a gas chromatograph as described above within 24 h after lake water sampling.Samples were kept in a refrigerator before analysis.The dissolved CH 4 concentrations were calculated with Henry's law, accounting for the temperature dependence of solubility (Sander, 2015).Due to logistical problems, these measurements were only conducted for three ST lakes (Bakchar-forest-2, Gavrilovka-1 and Plotnikovo).
We tested for potential diurnal variability in stratification and vertical mixing due to the difference between day and night temperatures that can occur in shallow lakes by examining the lake temperature and oxidation-reduction potential profiles.We found that in this study's lakes, these terms do not show strong vertical gradients, and the lakes all belong to the class of continuous cold polymictic lakes (Wetzel, 2001).We are therefore confident that our single daytime measurements will not generate a bias in the measured flux estimates (Ford et al., 2002).During flux measurements there were no periods with strong thermal stratification, as temperature gradients between surface and bottom water never exceeded 2 • C. In this concern the West Siberian lakes that were studied do not correspond to the summer pattern, described by Ford et al. (2002), in which diurnal flux cycles are coupled with afternoon stratification and night mixing.Rather, they correspond to the autumn pattern, in which there are neither strong gradients nor pronounced diurnal flux variability.We also neglected the storage flux because it is important only for stratified lakes, while all study lakes were mixed and no stratification was found during field campaign (Bastviken et al., 2004).
Submerged funnel gas collectors analogous to those in other measurement campaigns (Huttunen et al., 2001;Repo et al., 2007) were used to monitor CH 4 ebullition in the lakes.The gas collectors were 20 and 50 cm diameter funnels feeding into a graduated 500 mL polypropylene cylinder (Thermo Fisher Scientific, USA), which was fitted with a PVC tube to a 20 mL polypropylene syringe for sampling.One or two gas collectors were installed randomly near the study sites on the same day as the chamber measurements for 1-2 days in MT and 3-4 days in ST lakes.The net CH 4 ebullition was determined from the released gas bubble's volume and its concentration as analyzed on the gas chromatography flame ionization detector (GC-FID) described above.The short period of sampling used in our study does not allow the obtainment of highly accurate bubble flux estimates because the ebullitive methane flux showed strong temporal and within-lake spatial variability defined by weather events (Schilder et al., 2016;Wik et al., 2016a).These data are used only for comparison in a first approximation.
At each lake site, environmental characteristics were measured at three depth levels -20 cm below water surface, the lake profile midpoint and 10 cm above sediment depth.At each level, we measured air and water temperatures (T ) with Thermochron iButton loggers (DS 1921(DS -1922, Dallas Semiconductor, USA), pH, oxidation-reduction potential (Eh) using an SG-8 (Mettler Toledo, USA) and electrical conductivity (EC) using an SG-7 (Mettler Toledo, USA).At the same three depth levels, water samples were collected using PVC tube and immediately filtered in prewashed 30 mL PP Nalgene ® flacons through single-use 0.45 µm filter units MINISART (Sartorius, acetate cellulose filter) with a diameter of 25 mm.After discarding the first 20 to 50 mL of filtrate, the filtered solutions for cation analyses were acidified (pH ∼ 2) with ultrapure double-distilled HNO 3 and stored in pre-cleaned high-density polyethylene bottles.The sample storage bottles were prepared in a clean bench room (ISO A 10 000) and blank samples were used to check the level of pollution induced by sampling and filtration.Major anion (Cl − , SO 2− 4 ) concentrations were measured using ion chromatography (high-performance liquid chromatography (HPLC), Dionex ICS-2000) with an uncertainty of 2 %.Dissolved organic carbon (DOC) and dissolved inorganic carbon (DIC) were analyzed using a carbon total analyzer (Shimadzu TOC VSCN) with uncertainty below 3 %.Special calibration of the instrument for analysis of both form of dissolved carbon in organic-rich, DIC-poor waters was performed as described elsewhere (Prokushkin et al., 2011).Major cations (Ca, Mg, Na, K), Si and trace metal (Fe, Cu, Ni, Co) concentrations were measured with an inductively coupled plasma mass spectrometry (ICP-MS) Agilent ce 7500 with In and Re as internal standards and three various external standards, measured as check samples between each run of 10 lake water samples.Further details about analysis, uncertainties and detection limits are given elsewhere (Pokrovsky et al., 2015(Pokrovsky et al., , 2016)).Sediment layer depth was determined using a peat auger (accuracy is ±0.2 m).It should be noted that bottom pH served as a pH in the whole sediment layer since no data were available for the latter.This approximation bears the risks of misinterpretation because in some cases lake water pH can strongly differ from lake sediment pH, especially for acidic lakes (Casper et al., 2003).All studied acidic lakes have secondary origin, as they have the same oligotrophic sphagnum peat in the bottom (at least in the several upper meters that are only important for lake emission) as peat in wetlands around.Bottom-water pH of acidic lakes in our study is very close to pH of surrounding wetland according to available data for ST (Kotsyurbenko et al., 2004(Kotsyurbenko et al., , 2007) ) and MT (Sabrekov et al., 2011) wetlands.This peat has very low ash content (Turunen et al., 2001) and is unlikely able to strongly regulate pH.Therefore, it could be suggested that bottom pH is a fair proxy of pH in a whole sediment layer.
It is important to notice that the obtained flux and supporting data give only a momentary snapshot of methane emission from a certain lake section.While the spatiotemporal variability in CH 4 fluxes is critical when making whole lake or annual budgets (Natchimuthu et al., 2015;Wik et al., 2016a), our target was to show only momentary variability on a regional spatial scale.As it was mentioned above, the best way to take into account complicated and nonlinear key processes is to construct a process-based model.Modern lake methane emission models tend to operate not with whole lake and seasonal budgets but on much smaller temporal (days or hours) and spatial (with subsequent averaging) scales (Stepanenko et al., 2011(Stepanenko et al., , 2016;;Tan et al., 2015) because it helps to resolve small-scale heterogeneity of such important controls as lake depth or water temperature.Considering not the whole lake but certain lake sections can also improve aquatic greenhouse gas emission estimates (Schilder et al., 2016).Therefore, it should be mentioned that in our study not the whole lake but the lake section (area about 10 m 2 ) is the studied object.

Statistical analysis
All statistical analyses were performed using the Statistica 8 software (StatSoft, USA).Ordinary least square regression (α = 0.05) is used to find the significance of the relationship between each environmental variable and the measured CH 4 flux (for average and median values across all the individual fluxes for each lake).Stepwise multiple regressions (α = 0.05) included parameters such as air temperature ( • C), lake depth (LD, m), sediment depth (m), area (ha), CO 2 flux www.biogeosciences.net/14/3715/2017/Biogeosciences, 14, 3715-3742, 2017 (mgCO 2 m −2 h −1 ) and for three water depths (surface, middle, bottom) we measured temperature ( • C), pH, Eh (mV), EC (µS cm −1 ), concentrations of DOC, DIC, SO 2− 4 , Cl − , P, Fe, Cu, Ni, Co, and K (mg L −1 ).To estimate the powerlaw distribution parameters C and α in f (x) = Cx −α (in which f (x) is a probability density function and x is a flux), a maximum-likelihood estimate was used (Newman, 2005).The minimum chi-square estimation test was used to check how different probability density functions fit the flux data.In order to test the linearity between flux rank and absolute value of this flux value in doubly logarithmic coordinates (which is typical for systems with SOC behavior; Bak et al., 1987Bak et al., , 1988;;Jensen, 1998;Turcotte, 1999), the sample of methane fluxes was sorted (where rank 1 indicates the highest magnitude flux).The Levenberg-Marquardt algorithm was used to define coefficients of empirical dependences of methane production on pH and temperature.

Model description
To analyze the zonal difference between fluxes, a processbased model reproducing the effect of the main environmental controls that are well known from literature (such as temperature, pH, lake area and depth, DOC concentration, etc.) was developed.The model is designed to couple the pro-cesses of production, consumption, and transport of methane and consumption and transport of oxygen in the water column and sediments of shallow boreal lakes.The model structure is similar to other methane emission models for wetlands (Walter and Heimann, 2000;Tian et al., 2010;Meng et al., 2012) and lakes (Stepanenko et al., 2011(Stepanenko et al., , 2016;;Tan et al., 2015).The model structure is represented in Fig. 2, and a full description is given in Appendix A. Input model parameters include the temperature profile of the lake, concentrations of DOC, total phosphorous and the pH value where these values are measured in the near bottom water.The parameters are assumed to be representative for the sediment layer, lake and sediment depth, latitude, and wind speed at the 10 m height.The model outputs are methane and oxygen concentration profiles, methane ebullition rate, and diffusive flux of methane to the atmosphere.
The model is constructed similarly to other modern models that have shown good ability to predict methane emissions from lakes and ponds (Stepanenko et al., 2011(Stepanenko et al., , 2016;;Tan et al., 2015).There are several differences between our approach and these models.First, in our model, the necessary parameters were each obtained from published literature for the appropriate climate zone (where possible) and averaged across different sources.There was no calibration of model parameters because we try to test how current sci- entific knowledge about the methane cycle in boreal lakes can simulate the chamber-measured methane fluxes.In order to avoid using different calibrated constants relating the dependence of methane production from substrates (Stepanenko et al., 2011;Tan et al., 2015), DOC was selected as a single proxy for substrate of methane emission (Tian et al., 2010).In order to avoid calibrating the strongly variable temperature dependence of methane production and to take into account its potential climatic differences, a climate-sensitive approach was used (see Appendix B).Second, unlike previous models (Stepanenko et al., 2011(Stepanenko et al., , 2016;;Tan et al., 2015) we have added the influence of pH through a nondimensional scaling factor (Appendix B).Third, gaseous molecular methane diffusion in lake sediments is included in contrast with the previous models (Stepanenko et al., 2011(Stepanenko et al., , 2016;;Tan et al., 2015).It was introduced because initial numerical experiments demonstrated that taking into account only liquid CH 4 molecular diffusion in the pore space of lake sediments leads to a concentration of dissolved methane in lake water more than an order of magnitude lower than observed in several ST lakes in this study (see Sect. 3 and Table 2) and in other temperate and boreal lakes according to the literature (see Sect. 4.2 for details).Data regarding the gas-filled pore space in shallow lake sediments are very sparse (see Appendix A).Wetland gas-filled pore space occupies from 3 to 18 % of total peat volume (Fechner-Levy and Hemond, 1996;Rosenberry et al., 2003;Strack et al., 2005).Close values (7-18 %) were obtained in laboratory experiments with muddy lake sediments (Flury et al., 2015).Therefore, we assume that the lake sediment value of this parameter has the same order of magnitude.This choice of the model framework is based on the data availability, which covers a mix of both spatial and seasonal conditions.Recent models for lake methane emissions (Stepanenko et al., 2011(Stepanenko et al., , 2016;;Tan et al., 2015) are validated mostly against seasonal time series taken at singular locations.Thus, it is not clear whether the influence of spatial variability can be explained according to modern knowledge about environmental controls of methane emission or whether there are controls that are valuable on different spatial scales but not included in models.For example, controls that are relatively stable for a single lake and on a seasonal scale (climate, lake pH and trophic state, sediment porosity) may not be relevant on greater spatial scales.Since this paper's obtained flux data cover regional and local spatial variability, we use simple empirical relationships for controls that are known to be important on these scales: temperature (on a climate-sensitive basis), pH and DOC concentration (Le Mer and Roger, 2001;Nazaries et al., 2013;Serrano-Silva et al., 2014).The microbial communities of methanogens and methanotrophs and their dynamics were not simulated (as performed, for example, in Grant andRoulet, 2002, andKettunen, 2003) despite their importance because it is currently not possible to obtain reliable estimates of microbiological parameters for lakes with different pH and trophic states.Therefore, we compromise between the model's complexity (which cannot be overly detailed due to the challenge of obtaining reliable data for validation) and data availability (i.e., that the model should describe the influence of important and measured controls on the scale of the data that are present).
The partial differential equations were solved with MAT-LAB v. 7.8.0 (Mathworks, USA).A bootstrap method (Efron and Tibshirani, 1986) was implemented to find the uncertainty bounds on the modeled fluxes, as follows.First, artificial errors were introduced for each model parameter using their given standard deviations and a normal distribution.Then, 1000 iterations of these noisy parameter values were used to generate noisy flux estimates, and the uncertainty on the predicted flux value was derived as the standard deviation of these outputs.

Results
A summary of methane flux measurements is presented in Table 3.The median methane fluxes were 0.3 and 4.1 mgCH 4 m −2 h −1 for MT and ST lakes, respectively.For MT lakes the median relative uncertainty of the individual measurements was 10 %.For ST lakes the median relative uncertainty of the individual measurements was 20 % for fluxes higher than 1 mgCH 4 m −2 h −1 and 50 % for fluxes lower than 1 mgCH 4 m −2 h −1 .The higher relative uncertainty for smaller flux observations may be caused by the fact that smaller fluxes are potentially influenced by single rare bubbles.This impact creates higher scatter in the measured gas concentrations than the relatively constant bubbling that drives higher fluxes.The median fluxes for individual lakes vary from 0.1 to 0.5 mgCH 4 m −2 h −1 for MT and from 0.8 to 7.4 mgCH 4 m −2 h −1 for ST.Methane fluxes in both zones differ significantly (Wilcoxon test, p < 0.00001).Average values of pH, EC, and P differed between the two zones with the confidence level of 0.05 in contrast to average DOC, DIC, temperature, Eh, Fe and CO 2 flux, which were not significantly distinguished.Variability amongst repeated measurements for individual ST lakes was higher than in the MT: the coefficients of variation were 1.68 and 0.71, respectively.The average dissolved methane concentration for three ST lakes (Bakchar-forest-3, Gavrilovka-2 and Plotnikovo) is 8.3 ± 6.3 mgCH 4 m −3 (see Table 2).Simple regression showed that there were several correlations between environmental variables and either average or median CH 4 flux.The average CH 4 flux for all data positively correlates with CO 2 flux (R 2 = 0.43, p = 0.012) and surface [P] (R 2 = 0.30, p = 0.042).The average CH 4 flux for ST lakes correlates negatively with pH middle (R 2 = 0.48, p = 0.022), pH bottom (R 2 = 0.51, p = 0.020) and lake depth (R 2 = 0.39, p = 0.052).Median CH 4 flux for all data positively correlates with surface [Cu] (R 2 = 0.33, p = 0.030) and CO 2 flux (R 2 = 0.30, p = 0.043).Median CH 4 flux for ST lakes only correlates negatively with bottom [Cu] (R 2 = 0.44, p = 0.037).Multiple linear regression by flux median and the average with two or more independent variables did not show any reliable dependences (p ≤ 0.05).
The ST lake methane flux is more variable (both for individual lakes and for data combined for all lakes) than the MT lake fluxes (see Table 3).For example, the median coefficient of variation for average flux values from ST lakes (0.87) is more than twice the value from MT lakes (0.36).CH 4 flux values in MT and ST lakes are also from different continuous probability density distributions (two-sample Kolmogorov-Smirnov test, p < 0.00001).The combined ST lake flux values correspond to a power-law distribution with parameters C = 0.86 ± 0.05 and α = 1.71 ± 0.07 (minimum chi-square estimation test, p = 0.46; Fig. 3a) and do not correspond to a lognormal distribution (p < 0.00001).Conversely, the fluxes sampled from MT lakes have a lognormal distribution (mean = −1.46 ± 0.19, variance = 0.72 ± 0.12, p = 0.25; Fig. 3b) and do not have a power-law distribution (p < 0.00001).A linear dependence of flux rank on the absolute flux magnitude in doubly logarithmic coordinates (with the exception of a few points at the bounds) is also observed for methane fluxes from ST lakes (Fig. 3c) and is not observed for fluxes from MT lakes (Fig. 3d).Both power-law probability distributions and linearity in doubly logarithmic coordinates are typical for systems with SOC (Bak et al., 1988).Hence, SOC can be used to describe the dynamic of methane emission from lakes.
Since multiple linear regression did not reveal statistically significant dependences with two or more independent variables from the environmental factors listed in Sect.2.2.2, the multiple effect of environmental controls is confounding.Further analysis was provided using a process-based model (see Sect. 2.2.3 and Appendix A), which reproduced the methane and oxygen production, consumption and transport in lake water and lake sediments.
The modeling results are presented in Fig. 4 and in Table 4.The predicted fluxes fit the observed values for ST lakes quite well (R 2 = 0.76); however, the model overestimates MT lake fluxes by more than 1 order of magnitude.Modeled concentrations of dissolved methane in ST lakes vary in a wide range from 0.35 to 21.52 mgCH 4 m −3 with an average value of 7.63 mgCH 4 m −3 (recall our observations, in which for three ST lakes, the dissolved methane had a mean of 8.3 ± 6.3 mgCH 4 m −3 ).The modeled fraction of oxidized methane varies from 12 to 40 % with an average value of 22 %.Linear regression analysis of the residual dif-ferences between modeled and measured fluxes did not reveal statistically significant dependences from factors listed in the Sect.2.2.2.
Methane concentrations appear to be strongly underestimated (4-6-fold) for those ST lakes in which it was measured (Table 2).Our numerical experiments showed that the CH 4 molecular diffusion in liquids within the pore space of lake www.biogeosciences.net/14/3715/2017/Biogeosciences, 14, 3715-3742, 2017 sediments by itself could not generate a surface CH 4 concentration close to the observed values.These modeled values are very sensitive to the gas-filled porosity of the sediments (Table 2).Thus, the main differences between observed and predicted methane emissions are that the model overestimated fluxes for MT lakes by more than 1 order of magnitude underestimated concentration of dissolved methane in both MT and ST lakes (4-6-fold).
Additionally, the data showed extremely high variability in fluxes from ST lakes.Without additional flux monitoring and a greater focus on the driving process mechanisms, it may be that this experimental dataset is not suitable for a model comparison or validation effort.In the discussion section we try to suggest where these discrepancies have come from and how they can be explained.

Discussion
The obtained data indicate that CH 4 fluxes are distinctly higher in the ST than MT lakes.They are also correspond well with the data reported in these West Siberian zones by other researchers.Fluxes from the MT lakes agree with Repo et al. (2007) data from near Khanty-Mansiysk sites, where medians for two individual lakes were 0.2 and 1.0 mgCH 4 m −2 h −1 .They are also similar to the previously defined flux median of 0.6 mgCH 4 m −2 h −1 (across 51 flux measurements) for MT wetland lakes and ponds (Glagolev et al., 2011;Sabrekov et al., 2013).The median flux for ST lakes corresponds to the wide interval between the first and second quartiles (4.3-23.9mgCH 4 m −2 h −1 , across 82 flux measurements) reported earlier for the wetland lakes and ponds of the same climate zone (Glagolev et al., 2011;Sabrekov et al., 2013).
Comparison with measurements from small (< 100 ha) boreal lakes beyond West Siberia is presented in Table 5.Data for MT lakes are in the range of data from other regions in all parameters.Alternatively, methane flux (8.3 versus 1.3 mgCH 4 m −2 h −1 as the average sum of diffusive and ebullitive fluxes for lakes from Table 5) and DOC concentration (30 versus 11 mg L −1 ) in ST lakes are considerably higher, but the concentration of dissolved methane is lower (8.3 versus 15.0 mgCH 4 m −3 ).A comparison between the model and observational results can produce information for which controls and processes should be measured and examined to improve predictions of boreal lake methane emissions.

Differences in methane production between ST and MT lakes
The significant differences in measured CH 4 flux between MT and ST lakes can be explained with the help of the model results.Modeled fluxes from MT lakes are overes-Biogeosciences, 14, 3715-3742, 2017 www.biogeosciences.net/14/3715/2017/timated by 1 to 1.5 orders of magnitude.However, for ST lakes the model is close to the observations in both the mean level of emissions (at which the intercept of the observed vs. predicted flux linear regression is close to zero in comparison with the average flux; Fig. 4) and the representation of controls (in which the slope is close to unity).Hence, the question is what model parameter(s) should be changed for MT lakes to also correspond well.A first choice could be the potential controlling environmental parameters of CH 4 exchange.However, pH, DOC and temperature are in the same range in both zones (although their average values vary) and if the model interprets these values correctly in one zone it is unlikely to shift their interpretation in the other zone.Thus, differences in measured emissions between zones could come from parameters that we had modeled as unchanging between zones: the maximal methane production rate (MMPR; Eq.A14 in Appendix A) and/or the maximal intensity of methane oxidation (Eq.A15).It is doubtful that the key parameters of methane oxidation are different between the zones because the methanotrophic microbial community is assumed be adaptable to any climate and pH conditions and thus the oxidation rate is not heavily influenced by these factors (Le Mer and Roger, 2001;Nazaries et al., 2013;Serrano-Silva et al., 2014).Additionally, methane oxidation is proportional to methane production in lakes in which there is no influence of plant methane and oxygen transport (Bastviken et al., 2008;Duc et al., 2010).Therefore, differences in methane flux between zones are likely caused by differences in MMPR.
We therefore focus on MMPR, which may actually be lower in MT lakes compared to ST lakes due to substrate availability.While the mean difference in DOC between zones is not significant (15 mg L −1 versus 30 mg L −1 for MT and ST lakes, respectively, p = 0.142), if the two low-DOC and low-emitting Gavrilovka lakes were excluded from the ST sample, the difference would become significant (15 mg L −1 versus 36 mg L −1 for MT and ST lakes, respectively, p = 0.028).The DOC concentration is taken into account in calculations using the Michaelis-Menten equation (see Appendix A for details) and it can also influence the MMPR.MMPR implicitly reflects the abundance of methanogenic microbia.Higher substrate availability leads to higher methanogenic biomass and consequently to higher MMPR as simulated explicitly in previous research (Grant and Roulet, 2002;Kettunen, 2003).Greater substrate availability may be caused by higher plant productivity.Methane production correlates positively with plant productivity because root exudates provide additional fresh organic substrates for methanogens (Whiting and Chanton, 1993;Aulakh et al., 2001;Mitra et al., 2005).Net primary www.biogeosciences.net/14/3715/2017/Biogeosciences, 14, 3715-3742, 2017 production (NPP) in the MT wetlands that surround MT lakes is 40 % less than in ST wetlands (Peregon et al., 2008).This mechanism is extremely important for lakes with a low level of nutrients and in which a greater fraction of organic matter is allochthonous.The lower trophic state in MT lakes may also lead to a decrease in the MMPR.It is well known that higher trophic states generate both higher methane production and emission (Bubier, 1995;Segers, 1998;Duc et al., 2010).As a result, a critical concentration for bubble formation is either not attained at all or is attained at a lower depth in lakes with a lower MMPR.Therefore, a greater fraction of methane diffuses through the water column and is oxidized, further decreasing the flux.This low-production, low-ebullition hypothesis is supported in this study by measurements with bubble traps (see Table 3): the ebullition flux in MT lakes is less or equal to the diffusion flux calculated as the difference between the flux measured by static chambers and the flux measured by bubble traps.Meanwhile, the ebullition flux in ST lakes is many times higher than the diffusive flux from both model predictions and measurements.Certainly, our field experiments covered a relatively short period and were insufficient for exhaustively estimating methane emission pathways because we lacked the time to measure more bubbles.However, data by Repo et al. (2007) obtained in similar lakes with bubble traps during a month or more in summer are in good agreement with our measurements: no bubbles caught in a lake with a peat bottom (similar to lakes Lebedinoe and Babochka in the current study in which bubbles were also not detected) and small fluxes in lakes with sandy bottoms (0.04-0.4 mgCH 4 m −2 h −1 in Repo et al. (2007) and 0.01-0.1 mgCH 4 m −2 h −1 in the current study).
One can try to estimate the impact of these two possible causes: climatic and trophic.There are no data about the MMPR in lake sediments in West Siberia but we can estimate the climatic impact driving MMPR differences using data for ST and MT wetlands.According to Kotsyurbenko et al. (2004), the methane production under optimal temperature conditions and without substrate limitation measured in ST wetlands is 110 mgCH 4 m −3 h −1 .The same parameter for MT wetlands can be estimated from Kotsyurbenko et al. (2008) as 38 mgCH 4 m −3 h −1 .Thus, taking into account that both sites have similar pH conditions (because both wetlands are acidic and pH is 4.8 and 4.4, respectively) and trophic states (both wetlands are ombrotrophic Sphagnum bogs), the climatic methane production in MT wetlands is 3 times lower than in ST wetlands.The trophic impact can be estimated using average values of methane production for two groups of Swedish lakes (Duc et al., 2010).These groups were situated within the same region (thus had no climatic effect) and approximately correspond in phosphorous concentration to our MT (in Duc et al., 2010, this group of lakes is labeled low methane formation, P = 0.011 mg L −1 ) and ST (high methane formation, P = 0.064 mg L −1 ) lakes.Since phosphorous is known as a main control of methane production in lakes (Bastviken et al., 2004;Juutinen et al., 2009), in the first approximation it represents a difference in the trophic state between groups.For the low-methaneproduction lakes at optimal temperature, methane production is approximately 4 times lower than for the high-methaneproduction lakes.
Therefore, the sum of the climatic and trophic impacts gives a 12-fold reduction for the MMPR value for MT lakes in comparison with ST lakes.If we presume that the model's MMPR value is typical for ST lakes, the MMPR for MT lakes should be 2.60 mgCH 4 m −3 h −1 .Model experiments show that MMPR fitted to the corresponding measured methane fluxes from MT lakes is in the range 1.5-3.7 mgCH 4 m −3 h −1 .Hence, it can be expected that accounting for both climatic and trophic differences can help to adequately predict MMPR in a variety of lakes and reach a correspondence between measured and modeled fluxes for the MT zone.For these calculations and extrapolations, cross-ecosystem comparisons to evaluate the potential effects of both climate and local biogeochemistry on MMPR are needed.

Effect of diffusivity in the lake sediments
Model calculations show that only on average 22 % (12-40 %) of total produced methane is oxidized (see Table 3).The latter value is lower than the experimentally measured oxidized CH 4 fraction from Bastviken et al. (2008), reaching 22-40 % for the epilimnion of stratified lakes, similar to this study's lakes, because this layer is both oxic and has high turbulence.Both modeled and measured concentrations of dissolved methane in ST lakes are also less than the literature values for different temperate and boreal lakes (Table 5), including the West Siberian MT lakes (Repo et al., 2007); at the same time total methane flux in ST was considerably higher.This situation is in disagreement with the typical pattern in which higher methane concentration correlates with greater fluxes (as the sum of diffusive and ebullition flux).Thus, both the oxidized fraction of methane and concentration of dissolved methane seem to be underestimated.
The first possible reason for these differences is that the model has underestimated methane oxidation.Indeed, a comparison of half saturation constants for methane oxidation from different studies showed that this constant for highly productive CH 4 samples was greater than for samples with low production rates (Segers, 1998).However, the model parameters of methane oxidation, estimated based on several different sources (see Appendix A), correspond to measured and literature data on concentration of dissolved methane for studied lakes.Model experiments also show that increased oxidation leads to much lower modeled concentrations of dissolved methane, moving them even further from the measured values.Methane oxidation can also vary strongly depending on water chemistry or availability of different metals, such as Ni and Cu, which are main micronutrients necessary for enzyme production in both methanogens and methanotrophs (Krüger et al., 2003;Sazinsky and Lippard, 2015).Linear regression revealed that a higher copper concentration in the bottom waters corresponds to lower fluxes in ST lakes.It could be suggested that at higher copper concentrations, methane oxidation is enhanced in bottom waters, decreasing total emissions.The known mechanism of this possible enhancement is a "copper switch": in water at a Cu concentration < 50 µg L −1 soluble methane monooxygenase dominates activity, while Cu concentration > 250 µg L −1 is necessary to fully activate much more effective (for methane oxidation) particulate methane monooxygenase (Hakemian and Rosenzweig, 2007).Measured Cu concentration in MT and ST lakes does not exceed 2.5 µg L −1 making the copper switch hypothesis irrelevant in this situation.In the observed range of bottom-water copper concentrations, no significant effect on methane oxidation is seen (Schnell and King, 1995;Scheutz and Kjeldsen, 2005;Van der Ha et al., 2013).
This pattern could also be explained by the underestimated gas-filled porosity in lake sediments, which is an important control of dissolved CH 4 concentration influencing its diffusion through sediments (see Table 2).Literature data about gas-filled porosity in shallow lake sediments are sparse, while the variability in gas-filled porosity is very high (Valsaraj et al., 1999;Brennwald et al., 2005).However, it is well known that higher silt and mineral content as well as higher bulk density sediments have lower diffusion coefficients for the same values of gas-filled porosity (Clapp and Hornberger, 1978;Moldrup et al., 2003).This mechanism may explain why ebullition was not observed in MT lakes with peat bottoms and higher diffusion fluxes (Babochka, Lebedinoe) but was observed in lakes with lower organic content in their sediments and lower diffusion fluxes (Muhrino; see Table 3).The same situation was revealed by Repo et al. (2007), who found no ebullition and high diffusive fluxes in lake MTPond with peat sediments of several meters.Alternately, the lake MTlake of non-wetland origin was characterized by higher mineral content in sediments and lower diffusive flux.As a consequence, there was significant ebullition flux in this lake.
Accumulation of free gas affects the tortuosity of the sediment and leads to an underestimated diffusion coefficient for dissolved gas (Flury et al., 2015).The gas-filled porosity influence on methane cycling in lakes can be tested with a quick numerical experiment as follows.Consider doubling the gas-filled porosity for ST lakes to 0.05.This value is still typical for natural shallow lake sediments.For example, according to Valsaraj et al. (1999) the maximal gas-filled porosity is 0.07, a value more than 2 times higher than the 0.025 used by default in our model.In this higher-porosity case, the oxidized fraction of produced methane will increase to average 49 % (over a wide range from 19 to 90 %).The concentration of dissolved methane will increase to average 27.7 mgCH 4 m −3 , becoming 4 times higher than calculated using the default value of gas-filled porosity.Linearity between the observed and predicted fluxes under this experiment still remains high: predicted = 1.02 • observed -1.47 mgCH 4 m −2 h −1 ; R 2 = 0.73.Thus, both the underestimated fraction of oxidized CH 4 and concentrations of dissolved CH 4 can reach literature and measured values through natural variability in gas-filled porosity.
It could be concluded that the natural variability in gasfilled porosity in the sediments can strongly influence the ratio between diffusive transport and ebullition and, hence, the fraction of oxidized methane and total emissions.This variation may result from the extremely nonlinear influence of relatively low values of gas-filled porosity on gas diffusivity (Sallam et al., 1984;Flury et al., 2015).This nonlinearity is related to interconnected water films causing disconnectivity in gas-filled pore space and, thus, reducing gas diffusivity (Moldrup et al., 2003).Unfortunately, data about gas-filled porosity in lake sediments are very sparse and it is difficult to provide a comprehensive analysis of this parameter's influence on methane emission from lakes.

Emission uncertainty and self-organized criticality
The power-law dynamics of methane emission from ST lakes (Fig. 3) are similar to dynamic system behavior in the SOC theory (Bak et al., 1987(Bak et al., , 1988;;Jensen, 1998;Turcotte, 1999).SOC is based upon the idea that complex behavior can develop spontaneously in certain multicomponent systems whose dynamics vary abruptly.The paper by Bak et al. (1987) contained the hypothesis that systems that (i) are driven by some external force and (ii) consist of nonlinear interactions amongst their components may generate a characteristic self-organized behavior.The self-organized state into which systems organize themselves has properties similar to equilibrium systems at their critical point; thus, they are described as having SOC behavior (Bak et al., 1987).SOC dynamics are assumed to evolve through the contribution of processes on different timescales.The processes driven externally are typically much slower than the internal relaxation processes.A prototypical example is an earthquake, driven by stress that has slowly accumulated in the Earth's crust due to tectonic activity.This slowly built stress is subsequently released very quickly (in seconds or minutes) in an earthquake (Jensen, 1998).There is an analogous situation in lake sediments, as they become saturated by methane.Methane molecules and energy input continue much longer and more continuously than the release of bubbles and relaxation to the new steady state (Scandella et al., 2011).
The separation of relevant timescales is generated by the threshold responses -which build up over time -and metastability, which awaits a triggering event.In lake sediments, the situation is generated by microorganisms that produce and emit methane molecules into the surrounding lake water.The methane concentrations increase slowly until a solubility limit is reached.In this moment a new phase in the form of a bubble is produced.Then the methane concentration inside the bubble slowly continues to increase until the www.biogeosciences.net/14/3715/2017/Biogeosciences, 14, 3715-3742, 2017 moment when pressure in the bubble is high enough to work against forces preventing its release to the atmosphere (Scandella et al., 2011).When a critical pressure is exceeded, bubbles very quickly leave sediments via the previously formed channel.The applied force -the buildup of the CH 4 concentration -has to accumulate in order to overcome the critical threshold.This buildup occurs over a much longer timescale than the short time interval it takes the bubble to be released.The release of accumulated energy is nearly instantaneous in the moment the bubble moves.If the CH 4 molecules were produced very slowly by microorganisms and diffused in water without ebullition, then no threshold for motion would exist.In this situation the dissolved methane would be continuously released and its energy dissipated at the same rate as it was produced by the system.The actual force that the generated bubble of CH 4 must overcome depends on the molecular details of how the bubble interlocks with the sediment particles.As a result, there is a multitude of states in which the bubble will remain immobile even in response to an applied force.When these forces do not pass the release threshold, these states are metastable.The forces induce strain in the sediment material, corresponding to a certain amount of stored elastic energy.Thus, despite the bubble-sediment material system existing in an apparently stable, time-independent state, the system is not actually in its lowest energy state.A small increase in applied force can lead to a number of different responses from no motion of the bubble to a large jump and removal of the bubble from the sediment matrix (Scandella et al., 2011(Scandella et al., , 2016)).
There are several practical consequences of SOC behavior of methane emission in lakes.The high values of SD in Fig. 4 show not the low accuracy of measurements but natural spatial and temporal variability in methane emissions from lakes.Short-term measurements can produce uncertainty if they are extrapolated to a long time period or season.Controls found to be important from short-term measurements may be unreliable on other spatial or temporal scales.Whole season, multiyear measurements in three lakes in northern Sweden, made by Wik et al. (2013Wik et al. ( , 2014)), confirm this hypothesis.Each season of their measurements has a unique type of seasonal dynamic with a unique pattern of peaks and falls related to temperature and atmospheric pressure dynamics (Wik et al., 2013).However, the whole season methane budget clearly linearly correlates with seasonal energy input to lakes (Wik et al., 2014).
Another practical consequence is in the upscaling of flux measurements in lakes for large regions.Once we determine that the probability distribution law is relevant across all the ST lakes, we can use it for upscaling.The mean value for a power-law distribution is (Newman, 2005): in which x mean is mean flux value, x min is the minimal flux value, x max is the maximal flux value and the other param-eters were described in the Sect.2.2.2.While x min , C and α can be easily calculated based on our flux measurement campaign, it is more complicated to give a reliable estimate of x max because we cannot be sure that in our sample set we have obtained the maximal possible flux value.In SOC theory, x max is infinite, but in real lake conditions it is a function of methane production (as a measure for the applied external force) and sediment diffusivity (as a measure of energy dissipation).The maximal measured flux in ST lakes in our previous work (Glagolev et al., 2011;Sabrekov et al., 2013) is 359 mgCH 4 m −2 h −1 .If we assume that this value is the x max , then the mean flux value according to Eq. ( 1) would be 13.2 mgCH 4 m −2 h −1 .This value is 50 % higher than the simple average and 220 % higher than the median obtained from the flux dataset for ST lakes of the current study.Therefore, use of simple average and median statistics can lead to substantial underestimation of the total methane amount emitted from lakes.Despite this stochastic behavior of emissions, our modeled flux values are in good correspondence with measured fluxes.There are several reasons for this agreement.According to the probability law distribution identified for ST, 10 or more flux measurements, as we have performed, allow detection of high-flux moments (for example, for our ST flux power-law function the probability of detecting a flux with a magnitude from 10 to 20 mgCH 4 m −2 h −1 is 0.095) and represent the total emission in the first approximation.The wide range of fluxes and relatively high number of studied lakes also help to obtain good correspondence.

Other important controls for regional model development
Comparison of observed and predicted fluxes can help to reveal other important methane emission controls on a spatial scale.There are two strong site discrepancies for our model: CH 4 emission for the Ob' Floodplain oxbow lake is strongly underestimated, while the model generates a large overestimation for the Bakchar-forest-1 lake.Both these model-data disagreements can be considered in the context of organic matter quality.In our model it is assumed that organic matter in the form of DOC has the same quality for all lakes, but in reality the quality depends on its origin.There are generally two possible sources of this organic matter -plant and algae exudations and decomposition of organic matter (dead plants, different types of peat, gyttja, sapropel, etc.); both of them can be autochthonous and allochthonous (Whiting and Chanton, 1993;Cao et al., 1996;Segers, 1998).As a rule, fresh, labile and/or nitrogen-rich organic matter leads to higher methane fluxes (Segers, 1998;Duc et al., 2010).The Ob' Floodplain lake has the highest trophic state (in terms of P concentration) in our sample (see Table 1).Therefore, it is natural to suggest that higher trophic states produce higher MMPRs (in this case approximately 50 % higher) and hence higher emissions for this lake.Phosphorous does not directly influence methane production but strongly positively correlates with chlorophyll concentration, indicating productivity of algae, and with sediment respiration, indicating higher intensity of organic matter decomposition and higher oxygen consumption by sediments, as reviewed by Pace and Prairie (2005).Higher algae productivity (West et al., 2015) and peat decomposition supply methanogenesis with fresh organic substances, while lower oxygen concentration leads to decreasing methane oxidation.Moreover, temperature dependency of CH 4 fluxes actually increases with the overall system productivity (DelSontro et al., 2016).The Bakchar-forest-1 lake is mesotrophic but has the highest DOC concentration between the studied lakes and a CO 2 flux 2 or more times lower than other lakes with a relatively high methane flux (see Table 3).DOC in lake water positively correlates with both plankton and sediment respiration (Pace and Prairie, 2005).It can be suggested that the highest DOC and the lowest CO 2 flux for the same lake together demonstrate that this DOC is formed from recalcitrant organic matter leading to lower CH 4 production.This hypothesis corresponds with the findings of Duc et al. (2010) in which the similar values of MMPR correlate not with DOC concentration but with the quality of organic matter in form of the C : N ratio.
We decided not to compare residuals for the MT lakes because of the small sample size and, as mentioned in Sect.4.2, possible differences in gas-filled porosity.The latter parameter needs special investigation since now, without further datasets, it requires near arbitrary selection.It is interesting to compare our CO 2 and CH 4 flux data with data of Repo et al. (2007) for the same region.It was obtained that the MTlake site corresponds with our dataset's CO 2 and CH 4 flux values and ebullition-to-diffusive flux ratio.At the same time, much higher diffusive CH 4 flux and no ebullition were found for lake MTPond of 0.5 ha (see Table 5).The DOC concentration for this lake is not higher than in MT lakes studied by us, but 3-fold higher CO 2 flux can indicate better quality of this substrate for methanogenesis, as mentioned earlier in this section.The latter finding can be explained by the fact that MTPond is located in a through-flow poor fen and is partly vegetated (Repo et al., 2007).This setting means that methanogenesis in this lake is supplied by both autochthonous and allochthonous organic matter.In contrast, our study's MT lakes are surrounded by a pine-shrubsphagnum community that prevents any through-flowing and are not even partly covered with any vegetation.Therefore, carbon dioxide flux can be useful to predict methane fluxes from shallow lakes, as shown in other studies (Rasilo et al., 2015).Regression dependences of CO 2 flux from different controls that are reliable on a large scale (as obtained, for example, by Kortelainen et al., 2006) can increase the precision of global models of methane emission from lakes.
Another possible important control of CH 4 emission is the presence of chemical inhibitors.It is well known that a number of alternative electron acceptors (such as dissolved NO − 3 , Fe 3+ , Mn 4+ and SO 2− 4 ) inhibit methane production (Ehrlich, 1987;Conrad, 1996;Nealson and Saffarini, 1994).However, only in one studied lake (the Ob' Floodplain) did the concentration of an inhibitor (SO 2− 4 ) in the near bottom water exceed the threshold value of 1.92 mg L −1 reported in previous research (Kuivila et al., 1989).Despite this fact, emissions from this lake are even underestimated by the model.The discrepancy can be explained by the fact that an amount of sulfate can prevent methanogenesis only in a small part of the sediment layer, as discovered previously by Kuivila et al. (1989).The mechanism of inhibition corresponds with the study by Sabrekov et al. (2016), in which two groups of wetlands were distinguished in the forest-steppe zone (next to the southern climatic zone after ST) of West Siberia: one in which pore water EC is about 600-800 µS cm −1 and CH 4 fluxes are 2-4 mgCH 4 m −2 h −1 , and one in which EC is more than 2000 µS cm −1 and fluxes are not higher than 0.1 mgCH 4 m −2 h −1 .The maximal EC in the studied lakes is 500 µS cm −1 .Therefore, we can suppose that in the humid climate of the West Siberian taiga zone, there is a small probability of methane inhibition in lakes.This mechanism can be important for regions where evaporation exceeds precipitation, leading to a higher concentration of mineral components in lake water.Certain trace elements can be beneficial for methane production (Basiliko and Yavitt, 2001), but linear regression did not reveal significant correlation between concentrations of these elements in the lake water and model residuals.

Conclusion
A study of small-size bodies of water in the non-permafrost region of West Siberia has demonstrated that lake and pond methane fluxes vary on both regional and local spatial scales.Based on the presented model's calculations, it can be suggested that it is possible to predict fluxes for individual lakes within the same climate zone with a fair agreement by taking into account such established controls as temperature, pH and substrate availability.Individual characteristics of lake origin and development, such as sediment gas-filled porosity, trophic state and organic matter quality, can also have crucial effects on methane emission.
To successfully predict CH 4 fluxes in several zones, different values of MMPR should be used.The climate and trophic state may be primary controls of MMPR variability in the interzonal scale.Searching for simple governing relationships for MMPR on all spatial scales may be the most feasible manner of improving the precision of methane emission modeling.
The constructed ab initio model is much more primitive than more complex recent models (Tan et al., 2015;Stepanenko et al., 2016), but it does not include calibrated parameters because all parameters can be adopted from the literature as average values from several literature sources for the www.biogeosciences.net/14/3715/2017/Biogeosciences, 14, 3715-3742, 2017 suitable climate zone.It can be assumed that this approach can be effective for analysis of spatial variability in methane emission, which appears to be higher than the temporal variability (Treat et al., 2007;Olefeldt et al., 2013;Sabrekov et al., 2014).Additionally, controls of spatial variability seem to have lower predictive ability (for example in terms of R 2 for regression models) than for temporal variability (Treat et al., 2007;Olefeldt et al., 2013;Wik et al., 2013Wik et al., , 2014;;Rasilo et al., 2015).
For global modeling it is important to know which lakes, and with what kind of ecological features and in what season, there exists SOC behavior.These lakes can emit significantly more methane because methane bypasses the oxidation filter through ebullition.The most interesting question in this concern is about the limits of environmental controls in time and space that define the switch between ebullitive and nonebullitive regimes.Because of the variability in MMPR and diffusivity of lake sediments, the presence of such methane emission hot spots as small shallow lakes is expected in any climate zone.However, because of their great extent in the taiga and tundra regions, small lakes in those zones are particularly relevant for the global CH 4 budget.

Appendix A: Model description
The functional forms of process controls were chosen in order to obtain reliable estimates of the governing parameters using publically available information from the appropriate climatic zone.There was no calibration of any model parameters.
Oxygen and methane dynamics in the water column from the water-atmosphere border to the lower boundary of sediments were modeled using the following equations according to Tang and Riley (2014): in which F CH 4 (z) and F O 2 (z) (mg m −2 h −1 ) are the transport terms; Prod (z), Ebul (z), and Ox (z) (mg m −3 h −1 ) are the rates of methane production by methanogens, ebullition, and consumption by methanotrophs, respectively; Resp (z) (mg m −3 h −1 ) is oxygen consumption by plankton and sediment respiration; Phot(z) (mg m −3 h −1 ) is the oxygen production via photosynthesis; z (m) is the spatial coordinate (positive downward); and t (h) is the time.The coefficient 4 reflects the stoichiometric relationship for oxidation, in which for each 1 g of methane, 4 g of oxygen are necessary, according to the equation CH 4 + 2O 2 = CO 2 + 2H 2 O. Oxygen and methane diffusion can be written as (Stepanenko et al., 2011;Tan et al., 2015) in which D CH 4 (z) is the diffusivity for methane and C CH 4 (mg m −3 ) is the methane concentration in the liquid phase.The resultant diffusion coefficient was calculated as either the sum of the molecular diffusivity coefficient within a liquid and eddy diffusivity in the water column or as the sum of molecular transport within both the liquid and gas phases in the lake sediment layer (Tang and Riley, 2014): in which D 0,liq,CH 4 (m 2 h −1 ) is the molecular diffusivity of methane in lake water at 0 • C; T (z) ( • C) is the water and/or sediment temperature from field observations; D tur (z), D mol,liq,CH 4 (z), and D mol,gas,CH 4 (z) (m 2 h −1 ) are the wind-induced turbulent diffusivity in the water column, molecular diffusivity of methane in sediment pore water, and molecular diffusivity of methane through gas-filled pore space of lake sediments, respectively; (z) (m 3 m −3 ) is the total sediment porosity; ε a (z) (m 3 m −3 ) is the gas-filled porosity; α CH 4 (nondimensional) is the Bunsen solubility coefficient for methane; and z bot (m) is the depth to lake bottom.Molecular diffusion in both the liquid and gas phases was taken into account in the sediment layer as it is performed for wetlands (Walter and Heimann, 2000;Zhu et al., 2014).The Penman equation was used for molecular diffusion in the liquid phase because the fraction of pores in sediments filled with water was high (0.6-0.9).Thus, the equation is quite precise under observed porosity (Moldrup et al., 2000): Because gas-filled porosity in sediments is very low (0.015-0.07 according to Valsaraj et al., 1999, andBrennwald et al., 2005), the Penman equation is not accurate in these conditions (Moldrup et al., 2000).As a result, we used the Millington-Quirk equation (Jin and Jury, 1996), which generates a diffusion coefficient similar to experimental measurements conducted under low gas-filled porosity (Sallam et al., 1984;Moldrup et al., 2000): in which D 0,gas,CH 4 (m 2 h −1 ) is the molecular diffusivity of methane in air at 0 • C. We neglected the solubility effect in consistency with Stepanenko et al. (2011).The depth profile of (z) was adopted from Gadzhiev and Kovalev (1982) and used for all lakes.The same value of ε a (z) = 0.025 was chosen for all lakes as an average (Valsaraj et al., 1999;Brennwald et al., 2005).Since gas-filled porosity does not vary strongly with depth for the first 1-3 m of sediments (Brennwald et al., 2005), it was assumed to be independent of depth.α CH 4 was calculated as in Tang et al. (2010): ) in which K H,CH 4 (T (z)) (mg m −3 atm −1 ) is the temperaturedependent Henry's law constant for methane, calculated as (Sander, 2015) in which K 0,H,CH 4 (mg m −3 atm −1 ) is the Henry's law constant for methane at 25   (Arah and Stephen, 1998): in which V 10 (mg m −3 h −1 ) is the maximal sediment respiration rate at 10 • C, E (J mol −1 ) is activation energy of respiration and R (J K −1 mol −1 ) is universal gas constant.Since no data about solar radiation are available, photosynthesis is not taken into account.However, numerical tests show that oxygen does not limit methane oxidation in the water column.
Photosynthesis is the only process that produces O 2 in a water column.For its calculation, parameterization from Stefan and Fang (1994) was used.This parameterization assumes the rates of biogeochemical processes to depend on temperature and photosynthetic active radiation (PAR) and to be proportional to chlorophyll a concentration.PAR was calculated for a period of chamber measurements with the help of the simple model SPLASH v. 1.0 (Davis et al., 2017) using latitude, date and cloudiness (Russian Federal Service, 2017) as input parameters.Chlorophyll a concentration was calculated from total phosphorous concentration based on empirical function from Pace and Prairie (2005).For more details, an interested reader may refer to the original paper.
Ebullition was calculated under the assumption that emitted methane bubbles immediately reach the surface (Stepanenko et al., 2011;Tan et al., 2015).
in which c e (h −1 ) and α e (nondimensional) are empirical parameters and C cr,CH 4 (z) (mg m −3 ) is the critical methane concentration of bubble formation that has been estimated according to Stepanenko et al. (2011): in which C N 2 (z) (mg m −3 ) is the nitrogen concentration in the sediments according to the depth profile from Bazhin (2001), K H,N 2 (T (z)) (mg m −3 atm −1 ) is the temperature-dependent Henry constant for nitrogen and p a (Pa) is the atmospheric pressure.Methane flux through the bubbles was calculated by integration within the sediment layer: ) in which z sed (m) is the depth of the lower bound of sediments.Gas exchange between bubbles and ambient water was neglected because its relative impact on methane transport is very small for relatively shallow lakes (Stepanenko et al., 2011;Tan et al., 2015).
As for the boundary conditions, we specify zero flux for both gases at the lower bound: At the upper bound, we specified diffusive methane flux calculated according to Riera et al. (1999), Bastviken et al. (2004) and Rasilo et al. (2015): ) in which k CH 4 (m h −1 ) is the so-called piston velocity, an empirical gas exchange coefficient, and C eq,CH 4 (mg m −3 ) is the concentration of dissolved methane, corresponding to atmospheric concentration of methane using Henry's law: C eq,CH 4 = P CH 4 ,atm • K H,CH 4 (T (0)) , (A27) in which P CH 4 ,atm (atm) is the partial pressure of methane in the atmosphere above the lakes, calculated via concentration using the ideal gas law.k CH 4 was calculated as follows (Rasilo et al., 2015): Temperature sensitivity of the Schmidt number was calculated using interpolation of experimental data from Jähne et al. (1987) using third-order polynomial function.Parameter n was assumed to be −2/3 for wind speed < 3.7 m s −1 and −1/2 for wind speed ≥ 3.7 m s −1 (Crusius and Wanninkhof, 2003).The upper boundary condition for oxygen was calculated in the same way.
Appendix B: pH and temperature effect on methane production Since data about the temperature and pH sensitivity of methane production are highly variable (Dunfield et al., 1993;Segers, 1998;Meng et al., 2012), special consideration is required for these important controls.f prod,pH is a nondimensional multiplier reflecting the decrease in methane production under real, nonoptimal pH conditions.f prod,pH was calculated using data and the functional form from Meng et al. (2012) but with other coefficient values: We preferred not to use the coefficients given in Meng et al. (2012) because they strongly underestimate production in acidic conditions.Data obtained in a number of Russian lakes (Gal'chenko et al., 2001;Kotsyurbenko et al., 2004;Sabrekov et al., 2012) have shown that both production and emission can be very high even when pH is about 4 and lower.We suppose that this underestimation is caused mainly by lack of CH 4 production data for acidic wetlands.The coefficient of determination, R 2 for this dependence, given in Meng et al. (2012), is quite low (0.44) due to scatter in data.This scatter can be explained by variations in other methane production controls.In order to avoid both of these problems, we have obtained an empirical function of CH 4 production's pH dependence using data binned into 0.5 of pH unity intervals (Fig. B1).In order to obtain a function varying from 0 to 1, the fitted function was divided by its maximal value.Therefore, we have the term f prod,pH in the form of Eq. (B1).Fitted parameters are given in Table A1.f prod,T is a nondimensional multiplier reflecting the decrease in methane production under real, nonoptimal temperature conditions.f prod,T was calculated using the empirical function suggested by O'Neill et al. (1972) (presented within Strashkraba andGnauk, 1985) for describing the effect of temperature change on photosynthesis: T opt ( • C) is optimal temperature (i.e., when the process intensity is maximal); and Q 10 (nondimensional) is a parameter showing how many times the process rate will grow for each 10 • C increase in temperature (for low temperatures).
To constrain the C 1 and C 2 parameters, we used literature data (Svensson, 1984;Moore et al., 1990;Sass et al., 1991;Dunfield et al., 1993;Parashar et al., 1993;Klinger et al., 1994;Kotsyurbenko et al., 2004) about the intensity of methanogenesis under different temperature conditions in different climatic zones.T opt was calculated using linear regression from the average number of days per year with an average day temperature higher than 10 • C (N T >10 , days): T opt = 0.055 • N T >10 + 13.08. (B7) This empirical Eq. (B7) was found using the studies cited above as well as others (Cicerone and Shetter, 1981;King et al., 1981;Whalen and Reeburgh, 1988;Frolking and Crill, 1994;Best and Jacobs, 1997).The average annual number of days with an average day temperature higher than 10 • C had better correlation with T opt (R 2 = 0.62, p = 0.002) than other climatic parameters such as average daily air temperatures, or average daily air temperature minimums and maximums of January, July or the whole year.
T max was calculated as a function of T opt using linear regression (R 2 = 0.74, p < 0.001), carried out with the help of data from a variety of sources (Van den Berg et al., 1976;Zeikus and Winfrey, 1976;Williams and Crawford, 1984;Svensson, 1984;Sass et al., 1991;Miyajima et al., 1997;Kotsyurbenko et al., 2001Kotsyurbenko et al., , 2004)): T max = 1.023 • T opt + 15.29. (B8) When C 1 = 590 • C, C 2 = 1000 • C 2 and Q 10 = 2, the empirical Eqs.(B2-B6) fit almost all the experimental data (Fig. B2).The only exception (Fig. B2c) can be explained by the fact that the emission was measured at a site with low water table depth (10 cm below moss surface), and as a consequence the influence of temperature on methane oxidation is coupled with methane production.Usually the temperature optimum of methane oxidation is lower than that for methane production (Segers, 1998).Therefore, the left shoulder of experimental data in Fig. B2c lays below the expected, modeled values.

Figure 1 .
Figure 1.Studied lakes in MT (a, shades of yellow correspond to floodplain, brown to forests, green to pine-shrub-sphagnum communities, purple to ridge-hollow complexes) and ST (b, shades of green correspond to forests and grasslands, blue to wetlands) on Landsat satellite images.

Figure 2 .
Figure2.Schematic representation of the model structure.The one-dimensional column is divided into lake water and sediments.The forcing consists of lake and sediment depths, the water T water (z) and sediment T sed (z) temperature, DOC and phosphorous concentration [P], and pH.CH 4 production occurs only in the sediments.The methane production rate R prod (z) is a function of the sediment temperature; the DOC, which is taken as a measure for substrate availability; and pH.CH 4 oxidation is calculated in the different way for sediments and for lake water.The CH 4 oxidation rate R oxid (z) follows Michaelis-Menten kinetics for both methane and oxygen and is a function of the water and sediment temperature.The water respiration rate WaterResp(z) is a function of the water temperature and the phosphorous concentration, which is taken as a measure for abundance of phytoplankton.The aquatic photosynthesis rate Phot(z) is a function of water temperature, photosynthetic active radiation and phosphorous concentration, which is used as a substitution for chlorophyll a concentration through an empirical function.The sediment respiration rate SedResp(z) is a function of the sediment temperature.Both water and sediment respiration follow Michaelis-Menten kinetics for oxygen uptake.Transport of CH 4 in lake sediments proceeds by (1) molecular diffusion through gasfilled and water-filled pores and (2) ebullition, which is the formation of gas bubbles in the sediment layer and their immediate ascent to the water surface.Transport of CH 4 in lake water proceeds by (1) wind-induced turbulent diffusion and (2) molecular diffusion in water.Transport of O 2 is the same except for ebullition.The model calculates methane fluxes to the atmosphere and CH 4 and O 2 concentration profiles in the lake water and sediments.

Figure 3 .
Figure 3. Flux data for (a) ST and (b) MT lakes on a log-log scale and probability distribution fitting: (c) power law for ST lake fluxes and (d) lognormal for MT lake fluxes.Rank 1 indicates the highest magnitude flux.Note the strong difference in the y axis scaling between the two regions.

Figure 4 .
Figure 4. Observed versus predicted values of methane flux.Whiskers denote ±1 SD.Predicted flux uncertainties are calculated from bootstrapping (see Sect. 2.2.3) and are explained by high uncertainty in model parameters adopted from literature (see Appendix A).High magnitudes of observation SDs may be explained by ebullition and SOC behavior of lakes as methane sources (see Sect. 4.3).Note that predicted flux uncertainties are higher than the magnitude of observation SDs for the MT lakes but less for the ST lakes.

Figure B1 .
Figure B1.Empirical function of relative methane production dependence on pH.The function from Meng et al. (2012) is given for comparison.Whiskers denote ±1 SD for binned values of production, according to data presented in Meng et al. (2012).

Figure B2 .
Figure B2.Empirical function of relative methane emission (a-c) and production (d-h) dependence on temperature.Black squares are experimental data and solid lines represent the fitted empirical function using Eqs.(B2)-(B8).Standard deviations are given as whiskers for investigations in which they have been presented.

Table 1 .
Basic characteristics of study lakes (mean values of three measured values -bottom, middle, surface).

Table 2 .
Surface-dissolved water CH 4 concentration (mgCH 4 m −3 ) in three ST lakes, both measured and calculated assuming a range of values of gas-filled porosity.
*Value used in model by default (see Appendix A for details).

Table 3 .
Summary of field flux observations (empty cells indicate no data).
a Number of individual flux measurements.b Measured using bubble traps; if there were two replicate traps, standard deviation is given.

Table 4 .
Summary of modeling results.
* Difference between modeled total flux (as a sum of diffusive and ebullition flux) and measured average flux.

Table 5 .
Summary for temperate and boreal lakes with an area < 100 ha.Empty cells mean no data; storage fluxes are not shown and do not have values more than 1 mgCH 4 m −2 h −1 .

Table A1 .
List of the model parameters.