Process-based modelling of the methane balance in periglacial landscapes (JSBACH-methane)

A detailed process-based methane module for a global land surface scheme has been developed which is general enough to be applied in permafrost regions as well as wetlands outside permafrost areas. Methane production, oxidation and transport by ebullition, diffusion and plants are represented. In this model, oxygen has been explicitly incorporated into diffusion, transport by plants and two oxidation processes, of which one uses soil oxygen, while the other uses oxygen that is available via roots. Permafrost and wetland soils show special behaviour, such as variable soil pore space due to freezing and thawing or water table depths due to changing soil water content. This has been integrated directly into the methane-related processes. A detailed application at the Samoylov polygonal tundra site, Lena River Delta, Russia, is used for evaluation purposes. The application at Samoylov also shows differences in the importance of the several transport processes and in the methane dynamics under varying soil moisture, ice and temperature conditions during different seasons and on different microsites. These microsites are the elevated moist polygonal rim and the depressed wet polygonal centre. The evaluation shows sufficiently good agreement with field observations despite the fact that the module has not been specifically calibrated to these data. This methane module is designed such that the advanced land surface scheme is able to model recent and future methane fluxes from periglacial landscapes across scales. In addition, the methane contribution to carbon cycle– climate feedback mechanisms can be quantified when running coupled to an atmospheric model.


Introduction
Knowledge of atmospheric methane concentrations is a key factor in several global-scale environmental research fields.Besides acting as a highly potent greenhouse gas and thus influencing global climate change, methane also contributes to degrading the ozone layer.Its average atmospheric lifetime is about 12.4 years, and its current atmospheric concentration in the Arctic is about 1850 ppbV (Ito and Inatomi, 2012).Concentrations have been reported as rising slowly but steadily since the onset of industrialisation, and, after a hiatus at the beginning of the 21st century, have recently been found to be rising again.These recent dynamics in the global atmospheric methane budget are still not fully explained, emphasising the fact that future trajectories of methane and its role in global climate change are also highly uncertain.The global warming potential of methane is 84 to 86 times that of carbon dioxide over an integration period of 20 years and 28 to 34 times over 100 years (Myhre et al., 2013).Accordingly, even though its absolute mixing ratios are quite low compared to carbon dioxide, it makes up for about 20 % of the radiative Published by Copernicus Publications on behalf of the European Geosciences Union.
S. Kaiser et al.: Process-based modelling of the methane balance in periglacial landscapes with JSBACH forcing from all greenhouse gases.Thus, for the radiation balance and the chemistry of the atmosphere, it is important to understand land-atmosphere exchanges of methane.
Environmental conditions are highly heterogeneous in permafrost regions, where landscapes are often characterised by small-scale mosaics of wet and dry surfaces (Sachs et al., 2008).The heterogeneous aerobic and anaerobic conditions in permafrost soils, in concert with elevated soil carbon stocks (Hugelius et al., 2014), set the conditions for large and spatially heterogeneous methane emissions in these areas (Schneider et al., 2009).Such strongly varying environmental and soil conditions as well as processes that influence the methane production and emissions are challenges in a process-based model with a bottom-up approach for methane balance estimation, simply because of the complexity of the network of processes to consider as well as their unclear interactions.
However, process-based modelling approaches are powerful tools that help to quantify recent and future methane fluxes on a large spatial scale and over long time periods in such remote areas.They can give first estimates of where field measurements are missing and help to understand the effects of climate change on permafrost methane emissions.In addition, the effect of methane emissions on climate, and hence feedback mechanisms, can be analysed using an Earth system model.For such purposes, a methane module for an Earth system model has to be process-based and working under most environmental conditions, including permafrost.
Currently existing process-based methane models have usually been developed for applications in temperate or tropical wetlands, without considering permafrost-specific biogeophysical processes, such as e.g.freezing and thawing soil processes (e.g.Zhu et al., 2014;Schuldt et al., 2013).In other cases, they are embedded within a vegetation model, which cannot easily be coupled to an atmospheric model (e.g.Schaefer et al., 2011;Wania et al., 2010;Zhuang et al., 2004).Some models have been developed only for smallscale applications (e.g.Xu et al., 2015;Mi et al., 2014;Khvorostyanov et al., 2008;Walter and Heimann, 2000) or use an empirical approach (e.g.Riley et al., 2011).Highly simplified models might be less reliable for global applications (e.g.Jansson and Karlberg, 2011;Christensen et al., 1996) because of oversimplification in simulating the complexity of the methane processes.
The aim of this study is to introduce a new methane module that is running as part of a land surface scheme of an Earth system model.Moreover, it shall be general enough for global applications, including terrestrial permafrost ecosystems.The methane module presented in this work represents the gas production, oxidation and relevant transport processes in a process-based fashion.Among other processes, this new methane module takes into account the size variation of the pore spaces in the soil column in relation to the freezing and thawing cycles, influencing directly the methane concentration in the soil.Furthermore, in this module the oxygen content is explicitly taken into account, enabling two process-based oxidation processes: bulk soil methane oxidation and rhizospheric methane oxidation.
The platform chosen to develop the methane module is the JSBACH (Jena Scheme for Biosphere Atmosphere Coupling in Hamburg) land surface scheme of the MPI-ESM (Max Planck Institute Earth System Model).The starting point was a model version that has a carbon balance (Reick et al., 2013) and a five-layer hydrology (Hagemann and Stacke, 2015) and that includes permafrost as described in Ekici et al. (2014).A parallel development by Schuldt et al. (2013) incorporated wetland carbon cycle dynamics and was also integrated into the model version presented in this work.The bases for the methane-related processes were the works by Walter and Heimann (2000) and Wania et al. (2010).Special focus was also placed on the connections with permafrost and wetland as well as the explicit consideration of oxygen.This paper describes the newly developed methane module, and, for the purpose of model evaluation, it presents an application at a typical polygonal tundra site in the north of the Sakha (Yakutia) Republic, Russia.

Site description
For the purpose of evaluation, this model has been applied at the Samoylov island site, located 120 km south of the Arctic Ocean in the Lena River Delta in Yakutia, with an elevation of 10 to 16 m above sea level.The mesorelief of Samoylov is flat, while as microrelief, there are low-centre polygons with the soil surface about 0.5 m higher at the rim than at the centre.This results in different hydrological conditions also influencing heat conduction.The average maximum active layer depth at the dryer but still moist polygonal rims and the wet polygonal centres is about 0.5 m (Boike et al., 2013).While the water table at the polygonal rims is generally well below the soil surface, the polygonal centres are often water-saturated, with water tables at or above the soil surface (Sachs et al., 2008).
The vegetation on Samoylov can be classified as wet polygonal tundra that is composed of mosses, lichens and vascular plants.According to Kutzbach et al. (2004), mosses and lichens grow about 5 cm high and cover about 95 %, while vascular plants grow about 20 to 30 cm high and cover about 30 % of the area.The most dominant vascular plant, both at the rim and at the centre, is Carex aquatilis, but with a dominance of only 8 % at the rim compared to 25 % at the centre.However, most of the species present at the rim are different from those present at the centre.According to Sachs et al. (2010), the proportions of moist and wet microsites are approximately 65 % moist and 35 % wet.The reader is referred to Sachs et al. (2010) for more details on the study site.Below, moist microsites will be referred to as rim and wet microsites as centre.

Layer structure
For a numerically stable representation of gas transport processes in soils, a much finer vertical soil structure is required than what is normally used for thermal and hydrological processes in JSBACH.Therefore, a new soil layering scheme has been implemented for the methane module.This scheme is variable and allows fine layers (of the order of a few centimetres), but still inherits the hydrological and thermal information contained in the coarse scheme.The number and height of layers can be chosen arbitrarily, also allowing nonequidistant solutions.
Internally, the module uses midpoints and lower boundaries of the layers as well as distances between midpoints.At the bottom, the layering scheme is truncated at depth to bedrock.The layers where the plant roots end, i.e. the rooting depth lies, the water table lies and the minimum daily water table over the previous year lies (permanent saturated depth) have also been determined.These layers have a specific function for methane production and various transport processes.Details will be given below in the respective sections (see also Appendix A1).
For model evaluation, fine layers with a height of 10 cm have been used.For all the layers of the new soil layering scheme, the soil temperature is interpolated linearly from the coarse JSBACH layering scheme.From these values, the previous day's mean soil temperature is also calculated.In addition to geometry and soil temperature, each layer has its own hydrological parameters, as described in the next section, and various state variables describing the different gases' concentrations.

Hydrology
For the fine layers, several hydrological values have to be determined using the relative soil moisture and ice content from the coarse JSBACH layering scheme.Fine-scale layer values are derived such that known values at common layers are kept and only those layers that span more than one input layer get values of the weighted mean of the involved coarselayer values.The relative soil water content is then defined by the sum of the relative soil moisture and ice content.
Subtracting the relative ice content from the volumetric soil porosity leads to the ice-corrected volumetric soil porosity.With this, the relative moisture content of the ice-free pores can be defined, which is calculated by division of the relative soil moisture content by the ice-corrected volumetric soil porosity.Finally, the relative air content of the ice-free pores is defined as 1 minus the relative moisture content of the ice-free pores.
The water table is calculated following Stieglitz et al. (1997).From the uppermost soil layer, the water table is located in the immediate layer above the first one with a relative soil water content of at least 90 % of field capacity.This definition was used because the current hydrology scheme in JS-BACH does not allow one to consider water content of soils higher than field capacity or standing water (Hagemann and Stacke, 2015).Instead, water content exceeding field capacity is removed by runoff and drainage.In this context, the current model implementation considers only mineral soil (field capacity: 0.435; porosity: 0.448); i.e. no peat layers exist in this version.The dimensionless but ice-uncorrected field capacity is used because the relative soil water content already includes ice.The water table depth is then defined as Here, b is the lower boundary of the soil layer of interest with height h and relative soil water content r w .fc is the field capacity.If even the uppermost layer has a relative soil water content of at least 90 % field capacity, the water Evidently, this stratification is hydrological, while the layering scheme is purely numerical.Thus, each stratum may contain several soil layers.For carbon decomposition, the mean temperatures of the previous day at the midpoints of these three strata are needed.These values are derived analogously to the temperatures in the fine layers by interpolating the mean temperatures of the previous day linearly.
With these three strata, carbon that may experience unsaturated conditions is split into an unsaturated and a saturated pool by the water table.In addition, a permanently saturated carbon pool is defined by the permanently saturated depth.This scheme is similar to what Schuldt et al. (2013) proposed.Further details about the calculation of the carbon decomposition are given in Appendix A2.

Production
Initial values of methane and oxygen concentrations have been derived using reported gas concentrations in free air for oxygen and methane.For oxygen, the global mean value for 2012 is used (8.56 mol m −3 , http://cdiac.ornl.gov/tracegases.html).The value for methane is defined as the March 2012 value (77.06 µmol m −3 , http://agage.eas.gatech.edu/data.htm).
The initial gas concentrations in the soil profile are determined assuming equilibrium condition between free ambient air as well as the air and moisture in the soil pore space.Thus, Henry's law with the dimensionless Henry constant is applied.The dimensionless Henry constant is defined as the ratio of the concentration of gas in moisture to its concentration in air (Sander, 1999).The chosen temperature dependence values, which are d(ln = 1700 K, as well as the Henry constants at standard temperature, which are k 25 H CH 4 = 0.0013 mol dm −3 atm −1 and k 25 H O 2 = 0.0013 mol dm −3 atm −1 , are all from Dean (1992).
The calculated initial values for methane and oxygen concentrations in the soil profile can be transformed into gas amounts and vice versa.During methane transport process calculation, concentration values are widely used.In between time steps, however, the volume of ice is recalculated and therefore the relative ice-free pore volume changes.Thus, concentration values also change, but only the gas amounts stay constant.Therefore, at the beginning of each methane module execution, the total gas amounts that have been saved at the end of the previous time step are divided by the current relative ice-free pore volume to recalculate the current concentration values.
The final products of the decomposition of soil carbon are carbon dioxide and methane.Depending on the soil hydrological conditions, carbon dioxide or methane are produced from the decomposing carbon pools that belong to the three strata described above.These decomposition results are distributed over fine-scale layers of the whole soil column.Because no direct vertical information about the amount of decomposing carbon is available, equal decomposition velocity in all layers of one stratum is assumed.Thus, once the decomposed amount of carbon per stratum is known, the decomposed amount of carbon per layer per stratum depends on the amount of available carbon in that layer only.And the carbon content in the soil layers for Samoylov has been prescribed from measurements by Zubrzycki et al. (2013), Harden et al. (2012) and Schirrmeister et al. (2011), taking local horizontal variations of polygonal ground (Sachs et al., 2010) into account (see Appendix A3).
The initial amount of carbon in the pools is obtained from the sum of carbon in each layer of the strata.In this case, the first and second strata share one carbon pool which is split after calculation of the mean water table over the previous day.The amount of carbon per layer is divided by the amount of carbon per stratum.These weights are used for distributing the amounts of decomposed carbon from strata to layers.In addition, the share of initially produced carbon dioxide and methane is set assuming all decomposed carbon above the water table and half of it below the water table get carbon dioxide: Here, sl means all layers in the stratum, and C s is the decomposed carbon in the stratum.f C is the soil carbon content of the layer with height h, and v p is the ice-corrected volumetric soil porosity.Mass conservation is done if the stratum is too small to get a layer assigned, so that the associated carbon is not neglected.The gas fluxes for methane and carbon dioxide are calculated via the sums of the respective amounts, and the produced gases are added to their respective pools in the layers.

Bulk soil methane oxidation
Only part of the oxygen in the soil is assumed to be available for methane oxidation.In layers above the mean water table over the previous day, available oxygen is reduced by the amount that corresponds to the amount of carbon dioxide which is produced by heterotrophic respiration but not more than 40 % of the total oxygen content.An additional 10 % of oxygen is assumed to be unavailable and also reduced.In layers below the water table, the amount of oxygen is reduced by 50 %.This approach is similar to Wania et al. (2010).
For methane oxidation itself, a Michaelis-Menten kinetics model is applied.The Q 10 temperature coefficient is similar to the one used by Walter and Heimann (2000), but with a reference temperature of 10 • C rather than the annual mean soil temperature.Reaction velocities of both, methane and oxygen, are taken into account by using an additional equivalent term with the concentration of oxygen and K O 2 m = 2 mol m −3 , which is chosen to be the average concentration of oxygen at the water table.Furthermore, methane and oxygen follow a prescribed stoichiometry: (3) c denotes the concentration of oxygen or methane in the layer.T is the soil temperature in the layer, and dt is the time step.The total gas fluxes for methane, oxygen and carbon dioxide are again calculated as the sums of the respective amounts.

Ebullition
The implementation of the ebullition of methane largely follows the scheme from Wania (2007).Ebullition is the transport of gas via bubbles that form in liquid water within the soil and transport methane rapidly from their place of origin to the water table.The amount of methane to be released through ebullition is determined by that amount of the present methane that can be solute in the present liquid water.This amount depends on the overall amount of methane present in the layer, but also on the storage capacity of the present liquid water.
In a first step, the concentration of methane in soil air is assumed to be in equilibrium with the concentration in soil water.Thus, by application of Henry's law, the present methane can be partitioned into the potentially ebullited methane concentration in soil air and the potentially solute methane concentration in soil water.The dimensionless Henry solubilities at current soil temperature conditions are used for this.As an initial approximation, all methane is assumed to be in soil air and potentially ebullited.Thus, first, the potentially solute methane in soil water can be determined, but it will also be overestimated because of this approximation.Therefore, second, an updated potentially ebullited concentration of methane in soil air is determined by subtracting the potentially solute methane from the total methane.Unlike what was proposed in Wania (2007), these two steps are iterated until stable-state conditions are reached.
In a second step, to calculate the maximal amount of methane that can be soluble in the present soil water, the Bunsen solubility coefficient from Yamamoto et al. (1976) is applied.By considering the available pore volume, this gives the volume of methane that can maximally be dissolved.The ideal gas law results in the maximally soluble amount of methane.For that, the soil water pressure in layers below the water table needs to be derived.This is determined from soil air pressure and the pressure of the water column, using the basic equation of hydrostatics.For this, the specific gas constant of moist air and the soil air pressure in layers above the water table are required.For the air pressure calculation, the barometric formula is used.Hereby, the first layer uses the air pressure at the soil surface and deeper layers use the above layer's soil air pressure.The specific gas constant of moist air finally needs the saturation vapour pressure and relative soil air moisture, both in layers above the water table.The former is calculated following Sonntag and Heinze (1982), and the latter is set to 1 if the relative water content is at least at the wilting point and to 0.9 elsewhere.Now, the maximally soluble concentration of methane is derived by dividing the maximally soluble amount of methane by the available pore volume.Thus, the concentration of methane that is solute and in equilibrium with methane in the air is the lesser of the following two concentrations: the potentially solute methane that was calculated in the first step, and the maximally soluble methane that was calculated in the second step.Finally, the actually ebullited methane is the difference between all methane and solute methane, with k H CH 4 being the Henry solubility, c CH 4 gas the methane concentration that can potentially be ebullited, β the Bunsen solubility coefficient, p w the soil water pressure and T the soil temperature.All these variables relate to the layer, and R is the gas constant.
The ebullited methane is removed from the layers and, if the water table is below the surface, added to the first layer above the water table.In this case, the ebullition flux to atmosphere is zero, and the methane is still subject to other transport or oxidation processes in the soil.Otherwise, if the water table is at the surface and if snow is not hindering, it is added to the flux to atmosphere.Snow is assumed not to hinder if snow depth is less than 5 cm.If, finally, the water table is at the surface but snow is hindering, ebullited methane is put into the first layer and the ebullition flux to atmosphere is zero like in the first case.

Diffusion
For the diffusion of methane and oxygen, Fick's second law with variable diffusion coefficients is applied.The possibility of a non-equidistant layering scheme is specifically taken into account.Diffusion is a molecular motion due to a concentration gradient, with a net flux from high to low concentrations.For soil as a porous medium, moreover with changing pore volumes because of different contents of ice, the icecorrected soil porosity of the layers also has to be accounted for in the equation system directly as a factor (Schikora, 2012).The discretisation of the computational system is done with the Crank-Nicholson scheme with weighted harmonic means for the diffusion coefficients.While ice is treated as non-permeable for gases, the diffusion is allowed to continue if the soil is frozen but not at field capacity; i.e. there is no simple cut at 0 • C.During every model time step of 1 h, two half-hourly diffusion steps are calculated to prevent instabilities like oscillations or unrealistic behaviour like negative concentrations.The diffusion-specific time step can be decreased further if necessary and if an adjustment of the layering scheme is not desired.The possibility of these effects results from the tight connection between layering scheme, time step and diffusion coefficients.
As an initial condition, free ambient air, soil air and moisture phase are assumed to be in equilibrium.The boundary condition at the bottom of the soil column is always of Neumann type; i.e. no flux is assumed.At the top of the soil column, boundary conditions are assumed to depend on snow depth.If there are at least 5 cm of snow, no flux is assumed, and therefore the Neumann type also is applied at the top.However, if there are less than 5 cm of snow, ambient air conwww.geosci-model-dev.net/10/333/2017/Geosci.Model Dev., 10, 333-358, 2017 ditions are assumed to hold at the boundary, and therefore a Dirichlet type with a gas concentration in free air is applied: (5) Here, v p is the volumetric soil porosity, c denotes the gas concentration, t is the time, x is the depth, D denotes the diffusion coefficient, D is the boundary with Dirichlet type boundary conditions and N is the boundary with Neumann type boundary conditions.For details on how the diffusion coefficients are determined, see Appendix A4.The solution of the diffusion equation system is obtained by the tridag_ser and tridag_par routines from Press et al. (1996) in Numerical Recipes.By subtracting the gas concentrations after diffusion from those before for methane and vice versa for oxygen, concentration changes are derived with positive values for lost methane and gained oxygen.Multiplying the concentration changes by their respective pore volumes as usual and summing the resulting amounts over the layers gives the total fluxes of methane and oxygen.

Plant transport
Gas transport via plants is first calculated for oxygen entering the soil.Then, another oxidation mechanism with this newly gained oxygen takes place (see Sect. 2.2.8).After that, the transport of methane via plants is modelled.The transport via plants happens through the plant tissue that contains big air-filled channels, the aerenchyma, to foster aeration of the plant's roots.However, because plants need the oxygen that reaches their roots for themselves, their root exodermis acts as an efficient barrier against gas exchange.
In this model configuration, gas transport by plants is assumed to happen only via the phenology type grass with a C3 photosynthetic pathway.The contribution to methane emissions due to the degradation of labile root exudates is not taken into account here.The potential role of this process is reviewed in the discussion section.Furthermore, the gas transport via plants will occur only if snow is not hindering, i.e. if there are less than 5 cm of snow.This is justified by the consideration of snow crinkling the culms such that transport is not possible anymore.A diffusion process from aerenchyma through the root tissue to soil is assumed as a key process, and it is described by Fick's first law.Gas transport is fast inside the air-filled aerenchyma; hence, atmospheric air conditions can be assumed there.
The diffusion flux via the plants is determined from the oxygen concentration gradient between ambient air and the root zone soil layers.The diffusion coefficients of methane and oxygen in the exodermis are unknown but can be assumed to be slightly lower than in water (e.g.Kutzbach et al., 2004;Končalová, 1990).Therefore, their values are set to be 80 % of their respective values in soil water at the given soil temperatures and pressures, D r = 0.8 • D w .
The oxygen flux entering the soil is furthermore constrained by the surface area of root tissue, A ges r = A r • q p , which is determined from the surface area of a single plant's roots, A r = l r • d r • π, multiplied by plant density, q p = t ph t p .Here, l r is the root length, d r the root diameter, both in metres, t ph the number of tillers per square metre depending on phenology and t p the number of tillers per plant.Finally, the number of tillers per square metre is influenced by plant phenology, which is determined from the leaf area index (LAI), using t ph = max(t m ) • LAI max(LAI) , with t m being the number of tillers per square metre.Please also see Appendix A5.
The root tissue is assumed to be distributed equally between all root-containing layers; thus, A rl r = A ges r • h rl h , with h denoting the layer height and rl all layers with roots.The travel distance, dx, is set to the thickness of the exodermis in metres because this is the limiting factor.The plant transport per layer is thus modelled as Here, c O 2 air is the concentration of oxygen in free air and dt the time step length.For every soil layer, the resulting amount of oxygen is converted into concentration and added to the oxygen pool.As usual, the flux of oxygen into the soil is calculated by the total soil column balance.
After plant transport of oxygen, additional methane can be oxidised by the amount of oxygen that leaves the roots (Sect.2.2.8).The remaining methane is then available for plant transport, which is modelled exactly as for oxygen, with one exception: it is necessary to account for the fraction of roots able to transport gases, f r = dom Carex A.
dom Vascular P. .This can be thought of as a measure of distance between the methane and the transporting roots.With increasing amounts of roots being able to transport gases, the distance for methane to travel to them is getting smaller and transport is generally enhanced.To account for that, f r is set for rim and centre, respectively, as the fraction of the dominance measure for Carex aquatilis divided by the dominance of vascular plants (Kutzbach et al., 2004).The plant transport of methane is thus modelled as The variables' definitions are the same as for oxygen and c CH 4 air is the concentration of methane in free air.A similar effect will be taken into account for oxygen when it is allowed to oxidise only methane near the transporting roots.To determine the flux out of the soil, the differences of methane concentrations in the soil subtracted by the concentration in ambient air are used.For every layer, the amount of methane is converted into concentration and removed from the methane pool.Again, the total methane flux out of the soil is calculated by summing up individual layer balances.

Rhizospheric methane oxidation
The oxygen gained by the transport via plants is assumed to foster methane oxidation next to their roots.Thus, if oxygen is leaving these roots, the same oxidation routine as described above in Sect.2.2.4 is applied to calculate how much additional methane is oxidised by this oxygen.Obviously, only gas concentrations in layers with roots will be influenced.Because the amount of vegetation with roots that are able to supply oxygen varies between rim and centre, the dominance measure (f r from Sect.2.2.7) is applied again as a factor to account for the distance to these roots: The variables' definitions are the same as for the bulk soil methane oxidation, f r is the fraction of roots in the layer that are able to transport gases, and c O 2 plant is the concentration of oxygen transported by plants.Carbon and oxygen pools are adjusted accordingly.The total exchange with the atmosphere is determined by summing the total amount of gas that is calculated by multiplying the concentrations by their pore space.

Simulation set-up
As a global land surface scheme, the JSBACH model is set up for spatially explicit model runs at larger scales.Accordingly, many assumptions behind the model structures are only valid at large spatial scales.One prominent example here is the hydrology scheme, which works exclusively vertically, and therefore cannot represent lateral water flow from rim to centre, which is a process of major importance in polygonal tundra sites.Other examples include assumptions regarding e.g. the modules for radiation scheme and energy balance (no south-versus north-facing slopes, etc.).Since our ultimate target is to provide a new methane module that can be integrated into global-scale JSBACH simulations, accordingly the structure of our methane module also needs to target spatially explicit experiments.Thus, the site-level runs presented here are landscape-scale spatial runs with a grid cell size of 0.5 • using input data representing a very small domain.
To still facilitate site-level simulations that capture the general hydrologic characteristics of a polygonal tundra site, we split the model experiments into two separate runs, one for rim and one for centre.A redistribution of excess water from the rim area into polygon centres was added in order to mimic lateral flow.In more detail, the performed experiment consisted of two simulation runs with different settings for rim and centre.The polygon rim is assumed to be a normal upland soil, and a standard JSBACH simulation run was performed.For the polygon centre, runoff and drainage of the rim have been collected and added to centre precipitation.Additionally, for the centre run, runoff and drainage have been switched off until the soil water content reached field capacity.
The sequence of methane processes executed in the module is identical to the above-described order within Sect.2.2.1 to 2.2.8, and has been sorted according to the velocity of the specific processes, from fast to slow.The impact of changing this sequence on total and component methane flux rates was tested in a separate sensitivity study (not shown).These tests indicated only a minor influence of the sequence to the partitioning of the fluxes between the transport processes compared to the influence of hydrology or the definition of the processes themselves.Still, it cannot be excluded that modelled methane processes may be modified through the chosen order under certain conditions.
The carbon pools for rim and centre were initialised using data from Zubrzycki et al. (2013) and information from Harden et al. (2012), Schirrmeister et al. (2011) and Sachs et al. (2010).The used values for rim and centre for Samoylov are 627.61 and 731.94 mol m −2 for the upper carbon pool (i.e. the zone that is made up of the unsaturated and temporarily saturated soil layers) and 16 355 and 25 424 mol m −2 for the lower carbon pool (i.e. the permanently saturated zone).Because of the lack of information on how the modelled soil carbon from these two pools is distributed vertically, a depth distribution is applied to the decomposed carbon instead.For all layers within one stratum, equal decomposition velocity is assumed.The relative amounts of measured carbon are applied as a distribution aid for the decomposed carbon.The layers used were 10 cm in height.The only further settings varying between rim and centre are two vegetation parameters required for the process of plant transport, i.e. the number of tillers per square metre and the dominance of Carex aquatilis.Beyond the definitions cited above, the model has not been calibrated to site-specific processes or properties.
To initialise hydrological conditions, a spin-up of 100 years was done using 1 single year of climate data with average conditions from the period of observations.Starting in year 41 of this spin-up, the methane processes were activated.This set-up was chosen to stabilise the hydrological conditions before the methane processes were included.After finalising the spin-up, the time period of interest has been calculated with actual climate data.

Parameter sensitivity study
We reviewed the list of parameters that are required to run the new methane module of JSBACH and categorised them by relevance and available information to support the chosen settings.Based on this survey, we identified a shortlist of 10 parameters, which are listed in Table 2. To allow for a uniform processing of all parameters in this list, we assumed an uncertainty range of ±10 % for each of these settings.Changing each parameter by these percentages and performing for www.geosci-model-dev.net/10/333/2017/Geosci.Model Dev., 10, 333-358, 2017 each of those an individual model run yielded a range of resulting methane emissions according to the influence of each parameter.Model sensitivity towards the setting of the chosen parameters was evaluated through changes in the cumulative methane emissions over the modelled time period that followed the variation of the parameter.

Forcing and evaluation data
The climate forcing data used in the simulations are the same as in Ekici et al. (2014), spanning from 14 July 2003 to 11 October 2005.The climate input consists of air temperature, precipitation, atmospheric relative humidity, short-and long-wave downward radiation and wind speed, all at hourly resolution.
For model evaluation, data from chamber measurements have been used.These data were collected over 39 days from July to September 2006 by Sachs et al. (2010), resulting in 55 single measurements for the rim and 48 for the centre.In addition, eddy-covariance-based fluxes from Wille et al. (2008) have been used, integrating rim and centre.From this, 3340 data points were available for the simulation time period.

Modelled water table and permanent saturated depth
The modelled depth of permanent saturation for both, rim and centre, is always at the same level of 31.9 cm.In contrast, the modelled water table changes during the seasons for rim and centre differently (Fig. 1).In general, it is higher at the centre than at the rim, though there are few cases in early spring when the rim has a higher water table than the centre.This results from the different soil water contents at the rim and at the centre, which were forced by adding runoff and drainage from the rim to the centre as precipitation and prohibiting runoff and drainage at the centre until the soil water content reached field capacity.Still, in the early part of the thawing season, the water tables at the rim and at the centre are similar.While in general, at the rim, the water table is highest during the early thawing season, at the centre, there is a tendency to high values towards the end of the thawing season.But if the rim shows a high water table, there will generally also be a high water table at the centre.Overall, the water table in the model is changing relatively quickly, due to the quick changes in modelled soil water conditions.However, JSBACH does not allow one to model soil water content higher than field capacity or standing water at the surface.Thus, the maximal soil water content in the model is field capacity.It is obvious that there is a mismatch with the real situation in the field, where the centre is often water-saturated, with water tables at or above the soil surface.While measurements of the water table at the rim give values between 35 and 39 cm (Kutzbach et al., 2004), the mean summer value in the model is 30.88 cm.For the centre, measurements give values between −12 and 17 cm (Sachs et al., 2010), while the mean summer value in the model is 24.52 cm.Hence, the model tends to have a slightly higher water table at the rim, but the calculated water table is too low at the centre.Still, this water table has been calculated using the unsaturated soil water content.For the interpretation of the methane module results, it therefore has to be taken into consideration that JSBACH is currently not capable of filling the entire pore space up to saturation with water; i.e. a realistic representation of saturated water content like in the field is not possible.
For additional results concerning modelled physical conditions, such as soil moisture and ice content as well as soil temperatures, the reader is referred to Appendix B1 to B3.Table 1 gives the maximal values.
Table 1.Maximal values of the cumulative sums of modelled methane flux over the modelled time period for rim, centre and a mixed approach of 65 % rim plus 35 % centre for the different transport processes and combined in g C m −2 , rounded to three non-zero digits.

Rim
Centre

Modelled methane flux in summer and winter
The modelled methane fluxes at the rim and centre are different for the different seasons (Fig. 2).While most of the modelled flux is positive (i.e.emission to the atmosphere), there are also uptake events.The spread of the flux is greater for the centre than for the rim in both summer and winter.While the majority of flux values in summer is positive at the centre, it is more balanced at the rim.In winter, the methane flux is almost always zero, following the assumption that snow may hinder the exchange.However, at the centre, there are some rare events when uptake takes place.In the mixed approach, which means 65 % rim and 35 % centre, the overall mean emission is about 0.0813 mg C m −2 h −1 for the summer period.The overall higher emissions at the centre are due to higher moisture and thus more favourable conditions for methane production in concert with lower methane oxidation rates.

Role of different transport processes
During most of the year, the diffusive methane flux is rather small at the rim (Fig. 3a) and sometimes slightly negative at the centre (Fig. 3b).The largest methane emissions, both at the rim and at the centre, occur during spring.In this season, the methane that is produced in the topsoil from late autumn on and accumulated during winter is released in the form of so-called spring bursts upon snow thaw.Methane transport via plants is smaller than via diffusion, but more pronounced at the centre than at the rim (Fig. 3a  and b), because plant transport was defined as being slower than diffusion in water and should thus lead to lower emissions under less wet conditions.Despite the exodermis being a very thin layer, it is an efficient barrier against gas exchange, maintaining gases such as oxygen that are necessary for metabolic processes inside the roots.Thus, the diffusion rate through roots is slower than through water, and in turn, diffusion in water is much slower than diffusion in air.Moreover, the soils in the centre were not water-saturated in the model, promoting diffusive methane release through coarse pores.However, the wetter the soil, the more plant transport relative to diffusion should occur, because the more water, the more diffusion is slowed down.While ebullition is the most important process at the centre (Fig. 3b), it is diffusion at the rim (Fig. 3a).This is due to the drier conditions at the rim that allow a fast diffusion through air, while ebullition is only possible under conditions of high soil moisture.Because in the model higher soil moisture is calculated from the middle to the end of the thawing season, most of the emis-sions by ebullition and plant transport at the centre occur in this period (Fig. 3b).
In the mixed approach, only the diffusion of the rim alters the pattern of the emissions substantially (Fig. 3c).In total, the polygon centre accounts for a 6.8 times as large fraction of emissions as the rim due to the higher methane production under wetter conditions (Fig. 3d).This means a total share of 78.6 % of the methane emissions in the mixed approach is coming from the centre.Emissions at the rim are highest during spring, while they are highest at the centre during the mid and late season (Fig. 3d).
When comparing the total fluxes of the centre to the ones of the rim, diffusion is almost doubled, plant transport is 19 times as high, and ebullition is 18 times as high (Table 1).This results in almost 7 times higher total methane emissions at the centre than at the rim.At the rim, diffusion is more than 13 times as high as plant transport, while at the centre, it is just slightly higher than plant transport.Ebullition is about 4.5 times as high as plant transport both at the rim and at the centre.These differences are again due to the differences in soil moisture content, which allow more production under higher soil moisture and thus also lead to more methane emissions.On the other hand, plant transport is in principle a slower transport process than diffusion in water, but diffusion in water is much slower than diffusion in air.Thus, under drier conditions, diffusion in air will transport the main portion of gas, but under wetter conditions, plant transport may increase relative to diffusion.With reduced soil air, the remaining velocity of the diffusion is almost of the same order of magnitude as the overall velocity of plant transport, in contrast to the velocity of diffusion mainly through air.
Splitting the total methane flux into several transport processes not only allows one to evaluate the relative contribution of each process linked to rim or centre characteristics, but it is also possible to analyse differences in temporal patterns (Fig. 4a).As noted above, at the rim, the fluxes are much lower than at the centre (Fig. 4b), because less methane is produced under drier conditions, or methane becomes oxidised in the soil column.Ebullition makes up a large portion of the total budget at both microsites at isolated time steps, reflecting the nature of this process, while its total amount for the rim is rather small over longer time frames.At the rim, diffusion represents the second largest methane release but also substantial uptake during the season (Fig. 4a).The smallest flux portion at the rim is delivered by plant transport, which also shows some uptake.In contrast, at the centre, plant transport plays a much more pronounced role, and diffusion fluxes are more negative.All these effects occur in the different hydrological regimes at the rim and at the centre.
Furthermore, ebullition can only take place in soils with high soil moisture content, and this is more common at the centre than at the rim.Consequently, substantially more ebullition is found at the centre than at the rim.In the mixed approach, diffusion accounts for about 2.5 times of the emissions of plant transport, while ebullition accounts for 4.5 times of it.Overall, 0.588 g of carbon are emitted by each square metre during the modelled time period from 14 July 2003 to 11 October 2005.

Production versus oxidation
Methane oxidation follows the pattern of methane production as long as enough oxygen is available (Fig. 5a).Production, and hence also oxidation, is higher during times of more moist conditions for both, the rim and the centre, and is also higher for the centre than for the rim (Fig. 5b).At the centre, a substantial amount of methane is oxidised in the rhizosphere with oxygen that enters the soil via plant transport.This happens when a high amount of methane is produced, which is rather rare at the rim due to lower soil moisture (Fig. 5a).During spring, bursts of oxidation occur both at the rim and at the centre because methane produced during the winter and stored below the snow gets in contact with oxygen.The different moisture and temperature regimes at the rim and the centre and their dynamics determine these results.

Parameter sensitivity study
Results of the parameter sensitivity study are summarised in Table 2 and indicate that just one of the chosen parameters, fracCH4Anox, has a major influence on the cumulative methane emissions when varied within a 10 % range.
FracCh4Anox represents the fraction of methane produced under anoxic conditions compared to the total decomposition flux.For two more parameters, fracO2forOx+fracO2forPh and KmO2, the net effect was still larger than 1 %.FracO2forOx+fracO2forPh influences the available amount of oxygen for the methane oxidation, whereas KmO2 influences the oxidation as the Michaelis-Menten constant for oxygen.For all remaining parameters, only negligible effects on the cumulative methane emissions were found.

Comparison to chamber measurements
Although the number of available field data is small and from a different year than the meteorological forcing data, the field measurements and model results are of the same order of magnitude (Fig. 6).Observations and model results show higher centre values compared to the rim, but the model seems to underestimate occasional uptake events.For the rim, the model gives methane fluxes to the atmosphere of between −0.0237 and 39.3 mg C m −2 h −1 , with a mean of 0.0267 mg C m −2 h −1 , while the available field measurement values range from −0.111 to 0.881 mg C m −2 h −1 , with a mean of 0.154 mg C m −2 h −1 .For the centre, the model gives values between −0.0189 and 86.8 mg C m −2 h −1 , with a mean of 0.231 mg C m −2 h −1 , while the available field measurement values range from −0.0584 to 1.22 mg C m −2 h −1 , with a mean of 0.327 mg C m −2 h −1 .Besides higher mean values, the extremes are thus lower for the field measurements.This is due to the observation period excluding spring time, when the model calculates the highest emissions in the form of spring bursts.
One should also take into account that JSBACH is a global model; therefore, it requires input parameters from global fields.Furthermore, other modules of JSBACH, like the hydrology or the carbon decomposition, are adjusted for global applications.Therefore, JSBACH integrates processes over much larger grid cell areas than what chamber measurements may represent.Hydrological conditions and other processes are highly variable in polygonal tundra environments and are of crucial importance for methane processes.Still, they may not be represented with the required detail by the model, so that the modelled conditions are the same as those at the measurement site.Thus, it is obvious that with coarser and different hydrological conditions, the modelled methane fluxes per square metre for a 0.5 • grid cell cannot be identical to the point measurements of chambers.In particular, the low soil moisture in the hydrological conditions of the model may explain the lower mean modelled methane fluxes compared to what is reported by chamber data.

Comparison to eddy measurements
Eddy covariance data had the best available data coverage of field measurements (light grey areas in Fig. 7).Overall model results are of the same order of magnitude as observations, but there are also seasonal shifts between model results and measurements.This is due to a mismatch between the real soil conditions at the measurement site and the modelled soil climate and hydrology that cannot be expected to be the same as those in the field.The range of available measurements in the modelled period is 0.0233 to 4.59 mg C m −2 h −1 , with a mean of 0.609 mg C m −2 h −1 .The range of modelled summer methane emissions in this time frame is −0.023 to 30.4 mg C m −2 h −1 , with a mean of 0.0813 mg C m −2 h −1 .If less than 5 cm of snow are on the ground, this is defined as summer time.Besides lower mean values, the model also shows higher extremes.
For this comparison, the same constraints hold like for the comparison to chamber data.The modelled fluxes differ from field measurements because of differences in thermal or hydrological conditions.Critical are periods where observations show substantial methane emissions while at the same time model results show only minor emissions, e.g. in autumn 2003 or spring 2004.During these periods, modelled soil temperature values below zero and snow cover result in modelled methane fluxes of virtually zero, while in reality soils might be warmer and gas diffusion through snow might be possible (see Sect. 4).
Still, Fig. 7 also shows some patterns that are present in both model results and observations, e.g.periods with increasing fluxes that are followed by a sudden decline in the fluxes in a cyclic manner during a single season.These patterns are linked to the changing soil moisture content.Unfortunately only the first season is covered well by field measurements, while the second misses the later part, and the third covers just a part within.The model shows the largest methane emissions during spring upon snow thaw for both rim and centre in the form of burst.There is still little evidence in field measurements of the occurrence and magnitude of spring bursts, and to our knowledge no published data on this effect exist for Samoylov.In Sect.4, we briefly review the evidence of spring bursts in other northern wetland areas to evaluate the representativeness of these events in the model results.
For additional results concerning modelled oxygen uptake, such as the mixed daily sum, seasonally split and cumulative sums, as well as transport process split, see Appendix B4.

Discussion
This paper aims to present the structure of a newly developed methane module for the JSBACH land surface scheme and evaluate its general performance against field observations.The new module itself is completely integrated into the larger framework of the JSBACH model; therefore, sensitivity tests can only be conducted using the full model and a clean separation between the existing structure and new components is not always possible.The interpretation and discussion of all findings should therefore consider that the functioning of the new methane module is to a large extent dependent on, and in many aspects limited by, the performance of the JSBACH model as a whole.
The presented methane module determines production, oxidation and transport of methane to the atmosphere.All of these key processes are heavily dependent on soil moisture status as well as the quality and quantity of carbon in different soil pools.Both of these aspects, i.e. soil hydrology and carbon decomposition, are handled by existing JS-BACH modules which were not modified in the context of the presented study.With an exclusive focus on simulating processes at site-level scale, it may even be possible to upgrade these modules and add some features that would be relevant for the methane processes.However, since our scope was to provide a methane extension for JSBACH that can be applied globally, certain limitations regarding the representation of site-level observations need to be taken into account.This situation is even aggravated due to the use of parameter settings from global fields, i.e. with a coarse spatial resolution that aggregates conditions over larger areas and thus naturally cannot provide the exact details for the field site where the reference fluxes were measured.Such systematic deviations in modelling framework and parameter configurations will generate systematic differences between model output and site-level measurements.Accordingly, modelled hydrological conditions and amounts of decomposed carbon need to be considered when comparing modelled methane fluxes to the site-level observations and interpreting the spatiotemporal differences.
As mentioned above, the JSBACH hydrology module has been designed for global applications and is not capable of capturing conditions in complex landscapes such as polygonal tundra.Therefore, for the Samoylov site, which we used for this site-level analysis, the modelled soil climate and hydrology systematically deviate from those found in the field (Beer, 2016).We still chose to work at this site, because a highly valuable interdisciplinary dataset could be provided to evaluate different facets of the model output.To adapt the model to represent the complex hydrology, a mixed approach of combining two different model runs was applied.This approximation implies a very simplified representation of the real hydrological conditions and cannot fully offset all site-level differences between model simulations and observational datasets.Accordingly, systematic biases need to be considered when interpreting the findings.However, through this approach, we could demonstrate the paramount importance of realistic hydrologic boundary conditions for simulations of the methane balance.In many aspects, details in the behaviour of the methane processes are tightly linked to the spatiotemporal variation of hydrological conditions; therefore, biases in hydrology are directly projected onto the methane processes.
Still, the authors believe that the comparison of methane simulations against selected site-level measurements is an important first step to evaluating the overall performance of the new module.It is obvious that the limitations of the observational database employed herein, i.e. using just one single observation site and focusing on the growing season alone, cannot allow for a comprehensive assessment of the newly implemented algorithms.Accordingly, the limited amount of available field measurements from chamber and eddy-covariance-based fluxes requires a careful interpretation when compared to model results, particularly regarding the evaluation of JSBACH as a process-based global biosphere model.For the Arctic domain, methane emissions during shoulder and winter seasons have been shown to add considerably to the full annual budget, an aspect that we cannot evaluate based on the given database.Moreover, the question of temporal and spatial representativeness is complicated by the discontinuous nature of the methane fluxes (e.g.Tokida et al., 2007a;Jackowicz-Korczyński et al., 2010;Tagesson et al., 2012).To overcome these limitations, in follow-up studies the authors plan to conduct model evaluations based on longer-term flux measurements, covering full annual cycles for multiple Arctic sites.
Even though eddy-covariance-based fluxes are regarded as the most reliable reference data source for longer-term sitelevel model evaluation, the influence of microsite variability in the area surrounding the tower clearly poses a challenge here.Particularly with respect to methane fluxes, pronounced variability in the distribution of soil organic matter and water content may lead to a mosaic of different source strengths.For the Samoylov domain, which is characterised by polygonal structures, the apparent differences between wet (centre) and moist (rim) areas were mimicked through the execution of two model runs with different settings.Still, the footprint composition of the eddy covariance tower might not match the mixed approach of 65 % rim and 35 % centre used for modelling (Sachs et al., 2010).Even though this mixture generally captures the composition of the larger area surrounding the tower, particularly when footprints are smaller during daytime, the reduced field of view of the sensors might focus on areas that are wetter or drier than the average.The concept of combining two separate model runs has to be regarded as an approximation to cope with the hydrological constraints of a global model on the one hand and the complex landscape on the other.
The model application for remote permafrost areas may also be limited by the availability of long-term and complete observations of meteorological data to be used as model forcing.Forcing data and methane fluxes are required for the same time period, which optimally lasts over 1 or more years.When going towards regional to global applications, this new model might be additionally compared to regional or global atmospheric inversion results (e.g.Bousquet et al., 2011;Berchet et al., 2015) or data-driven upscaling of eddycovariance-or chamber-based observations (e.g.Christensen et al., 1995;Marushchak et al., 2016).
Within the methane module presented in this work, the discretisation as well as the pore volume are variable.This requires that the time step of calculation and the diffusion coefficients must fit to the thicknesses of the soil layers.If not set up properly, instabilities like oscillations or unrealistic behaviour like negative concentrations may occur.However, because the new methane module has been designed to be flexible in this respect, adjustments can easily be made in case numerical problems arise.
A parameter sensitivity study (Sect. 3.5) shows that the uncertainty of the resulting methane emissions scales linearly only for one parameter with the uncertainty of that parameter.This parameter represents the amount of methane produced under anoxic conditions compared to the total decomposition flux of carbon dioxide and methane combined  (Segers, 1998), this parameter was chosen to be 0.5 in Eq. ( 2).In other models, this parameter is used as an effective parameter and has been tuned to match ultimate methane and carbon dioxide emissions from soil to the atmosphere in the absence of an explicit representation of oxygen and hence methanotrophy (Wania et al., 2010).
Regarding the assumptions concerning fluxes during winter time or plant transport, according to recent findings (Zona et al., 2016;Marushchak et al., 2016), the settings chosen within the context of this work might be oversimplifying the actual processes in the field.The implemented mechanism that prevents gas exchange with the atmosphere once the snow cover reaches a depth of 5 cm is a very crude approximation of the snow cover influence.It resulted from biases in the modelled hydrological conditions in winter, where freezing of relatively dry soils led to oxic soil conditions that facilitated methane transport into the soil.The next iteration of the model development will include a more sophisticated, process-based representation of methane diffusion through snow.This upgrade, however, needs to be coupled to a major restructuring of several model components and thus cannot be reconciled with the model version presented within the context of this study.
The implementation of the plant transport follows a mechanistic approach, but its definition is limited by the availability of observational evidence on e.g.diffusion velocities.Therefore, the parameter settings used in this study are subject to high uncertainty.The value for the diffusion coefficient in the exodermis was chosen to be 80 % of the diffusion coefficient in water (C.Knoblauch, personal communication, 2014).The subsequent gas transport within the aerenchyma is assumed to be as quick as diffusion in air.With this setup, the effective barrier of the root exodermis will limit the plant transport efficiency and therefore act as a dominant control for this emission pathway.The thickness of this barrier has a large influence on plant transport as well; i.e. a thinner root exodermis would lead to increased plant transport.While this parameter is relatively easy to define, the cumulative surface area of all gas transporting roots in the soil column is difficult to constrain.Considering our basic assumption that plant transport is slower than diffusion in water, the general patterns of flux processes and soil moisture for rim and centre conditions appear plausible.Regarding the quantitative flux rates, however, the fraction of the total flux emitted through plant transport in the model tends to be too low.With larger root surface leading to increased plant transport, we therefore could use this setting as a tuning parameter to improve this issue.However, the oxygen available to consume methane also plays another modulating role, particularly for plant transport.Accordingly, new observational evidence would certainly improve the associated uncertainties; therefore, this issue is subject to ongoing investigations.With the new methane module, designed to be flexible re-garding these kinds of settings, parameter adjustments with respect to newer findings can easily be implemented.
The contribution of labile root exudates to methane production and emission has been largely neglected in existing model implementations and is also not considered in this model configuration.This is also an understudied process in field experiments and can only be estimated indirectly.The rate of root exudates is linked to the nutrient availability in soils, with more root exudates present in plants located in nutrient-poor wetland soils (Koelbener et al., 2010).The wetland soils in Arctic tundra are known to be nitrogen-limited (Melle et al., 2015;Gurevitch et al., 2006).The plant growth in the polygonal lowland tundra of Indigirka, Russia, is colimited by nitrogen and phosphorus, and only about 5 % of the total nitrogen soil content is active in the biological fraction (Beerman et al., 2015).The presence of vascular plants in Arctic wetlands supports the production of highly labile low molecular weight carbon compounds which can promote methane emissions through their methanogenic decomposition (Ström et al., 2012).Indirect evidence of the role of root exudates in methane production in polygonal ponds and water-saturated soils in Samoylov is presented by Knoblauch et al. (2015).The authors found almost 4-fold higher potential methane production rates in vegetated sites compared to the non-vegetated ones, both with the same C and N soil concentrations.Thus, the contribution to methane emissions from wetland soils in Arctic tundra due to the decomposition of root exudates should be taken into account in models.This will allow the understanding of the role of root exudates under present climate conditions.On the other hand, the potential nutrient mobilisation in soils due to permafrost degradation under climate change (Kuhry et al., 2010) may reduce the role of root exudates in methane emissions.However, the current JSBACH configuration lacks a full soil nutrient cycle, and the assimilation of nutrients by plant roots, as well as the contribution of root exudates to the total methane emissions, cannot be modelled at this point.
In Samoylov, the minimum of modelled daily sums of methane emissions during summer is smaller and the maximum much higher for rim and centre compared to measurements published by Kutzbach et al. (2004).However, these observations do not include spring bursts with very short but also very high emissions or even dry phases with small uptake.On the other hand, the mean of those measurements is 3 times as high for the rim and 3.5 times as high for the centre compared to the modelled daily sums in summer (Table 3).But such high modelled emissions are rather rare, and the general level of modelled values is lower than in observations (Fig. 7).
When comparing our model results at Samoylov to published results from other high-latitude regions, reasonable agreement is found.Our modelling results are about 40 to 60 % lower than measurements for BOREAS, Canada, and Abisko, Sweden, (Wania et al., 2010).The Lena River Delta region is much colder and drier compared to these sites, suggesting that lower flux rates are indeed reasonable.Furthermore, the Samoylov site is characterised by mineral soils containing substantially lower organic carbon as a substrate for methane production than the organic soils at the BOREAS site and the mire in Abisko.Compared to measurements done by Desyatkin et al. (2009) on a thermokarst terrain at the Lena River near Yakutsk, our mean results are well within the measurement range when comparing our rim to the drier sites, our centre to the wetter sites, and our mixed approach to the entire ecosystem (Table 4).However, climate and environmental conditions in this study were very different from those observed in Samoylov; thus, this comparison can only be regarded as a rough guideline.Nakano et al. (2000) measured methane fluxes at Tiksi near the mouth of the Lena River.While our mean value at the rim is 4.5 times as high as the mean measurements in Tiksi, the mean at the centre is 5.5 times as high as our mean value (Table 3).The modelled minimum is lower for the centre but comparable for the rim.
The large methane spring burst simulated by the model at both the rim and centre may represent the release of methane that has been accumulated during winter in the topsoil below the snow layer.To our knowledge, there is no observational reference of spring bursts measured in Samoylov.However, evidence of these events have been presented for other wetland areas using chambers and eddy covariance measurements, e.g. in northern Sweden (Jammet et al., 2015;Friborg et al., 1997), in Finland (Hargreaves et al., 2001), in northern Japan (Tokida et al., 2007b) and in Northeast China (Song et al., 2012).These studies suggest the presence of spring thaw emissions of methane that occur sporadically over short periods in the form of bursts.The magnitude of the spring bursts can exceed the mean summer fluxes by a factor of 2 to 3.Although spring emissions can account for a large share of the total annual fluxes, their occurrence, duration and magnitude are still uncertain.To adequately characterise the spring bursts in Samoylov, it is necessary to perform dedicated field measurements during the spring thaw period.These results will then help to evaluate the representativeness of the modelled spring bursts.In future model iterations, the spring bursts will also be evaluated for larger spatial scales.
In Zona et al. (2009), several measurements of methane emissions in the Arctic tundra are presented.Despite our mean values being located towards the lower end, our minimum, mean and maximum values fit well within the given range.Bartlett et al. (1992) measured methane fluxes near Bethel in the Yukon-Kuskokwim Delta, Alaska.The provided values for upland tundra compare well to our mean and minimum values.However, the model maximum fluxes are higher than the measurement values for upland tundra, but still well within the range of measured values for wet meadow, which has higher moisture contents than upland tundra.In fact, the highest values are calculated if soil moisture is highest, so despite being more on the lower end of this waterlogged landscape type's emissions, they also fit well therein.In summary, the variability of results of this pan-Arctic survey indicates that methane budgets within all these places are influenced by different conditions in terms of weather, hydrology and carbon pools.Accordingly, the good agreement of our modelled values with these references confirms that our results are within a plausible range at the greater picture, but a detailed evaluation cannot be performed without in-depth analysis of the site-level conditions.
Regarding the general structure of the JSBACH model, other parts of the land surface scheme require advancements before its application with the methane module at a global scale and over long time periods can be suggested.For example, soil organic matter should be represented as vertically resolved (Braakhekke et al., 2011(Braakhekke et al., , 2014;;Koven et al., 2015;Beer, 2016), with different soil carbon pools and a moisturedependent decomposition.Furthermore, the site hydrology should include soil moisture contents above field capacity and standing water above the surface (Stacke and Hagemann, 2012).We are also aware, however, that it is not the best approach to calculate an empirical water table depth following Stieglitz et al. (1997) under unsaturated soil water conditions.Together with the water table depth, the soil moisture content itself is of great importance to the presented methane module.Still, with this model version, the importance of different processes, their interplay and the influence of climatic or hydrologic drivers can be studied at site level, which is a major step forward.Furthermore, this process-based implementation can be applied at other sites or with another hydrology, and still, the methane-related processes will only depend on the soil conditions.In order to improve the hydrological www.geosci-model-dev.net/10/333/2017/Geosci.Model Dev., 10, 333-358, 2017 scheme of the current model version, it would be desirable to use other approaches like TOPMODEL (e.g.Kleinen et al., 2012) that would allow one to represent the fraction of the inundated area in a model grid cell based on the topography profile.This would provide a modelled wetland extent and a representation of the water table depth in saturated soils, especially for large-scale applications.This step has been considered and will be included in future model iterations.
Despite being a complex process model, the interplay of the processes is consistent.Thus, the influence of climate and hydrology on methane fluxes can be studied in detail.Knowing the dominating processes and environmental conditions provides useful information about the complex behaviour of the methane dynamics in permafrost soils.To summarise, a lot of information can be gained from using this model that may all help understand the complex network of drivers, influencing factors and constraints that govern the methane balance in periglacial landscapes.

Conclusions
The aim of this study was to develop a more detailed and consistent process-based methane module for a land surface scheme which is also reliable in permafrost ecosystems.
Based on previous work by Wania et al. (2010) and Walter and Heimann (2000), the JSBACH land surface scheme of the MPI-ESM global Earth system model has been enhanced for this purpose.The new methane module of JSBACHmethane represents methane production, oxidation and transport.Methane transport has been represented via ebullition, diffusion and plant transport.Oxygen can be transported via diffusion through soil pores and plant tissue (aerenchyma).Two methane oxidation pathways are explicitly described: one takes the amount of soil oxygen into account and the other explicitly uses oxygen that is available via roots (rhizospheric oxidation).All methane-related processes respond to different environmental conditions in their specific ways.They increase or decrease according to their requirements with changing soil moisture, temperature or ice content.The differences between the processes, seasonal differences as well as differences between the rim and centre microsites have been shown.When combined with a module for water-saturated soil conditions like TOPMODEL (e.g.Kleinen et al., 2012), such a methane-advanced land surface scheme can be used to estimate the global methane land fluxes, including for periglacial landscapes.These regions are rich in soil carbon (Hugelius et al., 2014) and show good conditions for methane production (Schneider et al., 2009).However, they are often remote and rather hard to investigate.Thus, process-based modelling can contribute to understanding the role of methane emissions as long as widespread and long-term measurements remain scarce.In addition, the role of methane for future permafrost carbon feedbacks to climate change can be studied.For these reasons, the module in this study is also highly integrated with permafrost and wetland processes, e.g.changing pore space in the soil because of freezing and thawing or changing water table depths due to changing soil water content.In a first comparison with site-level field measurements, sufficiently good agreements could be shown, despite the module not having been adjusted to site-specific processes or features.Coupling such a land surface scheme to atmosphere and ocean schemes in an Earth system model will provide the basis for studying methane-related feedback mechanisms to climate change.

Code availability
The model code used in this work is available upon request for academic and non-commercial use.
Appendix A: Additional methods A1 Layer structure -specific layer determination Specific layers are determined by comparing the midpoints of the layers to rooting depth, water table or minimum daily water table over the previous year, respectively.If one of these lies between two layer midpoints, the layer with the upper midpoint is chosen to be the specific layer for that.If the depth under consideration and the midpoint of a layer are the same, the corresponding layer is chosen.

A2 Hydrology -decomposition of carbon
The decomposition of carbon is determined similarly to Schuldt et al. (2013), though appropriate temperatures are used for each of the three strata.Furthermore, the decomposition times for the three carbon pools have been adjusted to ensure that the two pools under partially oxic conditions are relatively stable, neither accumulating nor decomposing great portions within a few years, and the last pool slowly accumulating.In numbers, the former two pools change only about 1 mol m −2 each within the calculation period from 14 July 2003 to 11 October 2005.The decomposition timescales used are 80, 400 and 30 000 years for the unsaturated, currently saturated and permanently saturated stratum's carbon pool.
Though the rate of organic matter decomposition at the evaluation site is not known, the present-day amount of carbon in the soil is known (Sect.2.2.3).Considering short timescales only, the above-described approach should give reasonable amounts of decomposed carbon in the three strata.This way, the input to our methane routine, the amount of decomposed carbon per time step in each stratum, is provided daily.

A3 Production -soil carbon per layer
The amount of soil carbon per layer has been prescribed based on measurements for the first metre of the soil profile by Zubrzycki et al. (2013).The values of the six measurement depths were averaged over the 16 different centre and 6 rim cores.These resulting averages have been interpolated to 1 cm values for rim and centre accordingly.The means of the corresponding 1 cm values are then used for the modelling layers within the first metre of the soil profile.
As Zubrzycki et al. (2013) only give values for the first metre, additional information for the rest of the soil profile is needed.Schirrmeister et al. (2011) give an estimate for Lena Delta soil carbon content of 553.33 kg m −2 with a soil depth of 18.25 m, which is converted in a volumetric estimate of 30.32 kg m −3 .Harden et al. (2012) give quantitative information about the depth distribution of soil carbon up to 3 m.Horizontal variations are accounted for by a partitioning into 65 % rim and 35 % centre (e.g.Sachs et al., 2010).
Using this information, values are assigned to the remaining layers so that the overall mean over all layers, rim and centre mixed in the proposed partitioning, gives the volumetric estimate gained from Schirrmeister et al. (2011).Hereby, the information from Harden et al. (2012) about the variability over depth, which is a slight decrease until 1.7 m and a slight increase thereafter, is taken into account.
As the uppermost values for this, at a depth of 1.05 m, the means of the deepest measured values are taken as 21.24 kg m −3 for the rim and 35.00 kg m −3 for the centre.As values at the turning point, at depths of 1.65 to 1.75 m, the ceiled mean values of the first metre are used, which are 20 kg m −3 for the rim and 34 kg m −3 for the centre.In between, the values are interpolated linearly, and then, towards the depth, extrapolated linearly, to meet the criterion of overall fitting to the value of Schirrmeister et al. (2011) as mentioned above.

A4 Diffusion -diffusion coefficients
Following Collin and Rasmuson (1988), the diffusion coefficients of methane and oxygen in the soil layers are calculated by adding the diffusion coefficients in soil moisture times the dimensionless Henry solubility to the diffusion coefficients in soil air.Both are weighted by the relative pore moisture or air content, and the ice-corrected soil porosity of the modelling layers is also considered.The exponents for this are estimated with Newton's method.For fast convergence, an appropriate starting value has been chosen that was found to be 0.62.The dimensionless Henry solubilities for methane and oxygen at the current soil temperatures are applied, and the diffusion coefficients in soil air and moisture are derived.
The diffusion coefficients in soil air can be seen as such in free air at soil temperature and pressure.They are calculated following Massman (1998) from values at the soil surface with depth-variable soil temperature and pressure.The latter one arises from soil air and water pressure.The values of diffusion coefficients in free air at the soil surface are calculated from values at 0 • C and 1 atm (Massman, 1998).
The diffusion coefficients in soil moisture can be seen as such in free water at soil temperature and pressure.They are calculated differently for the two gas species.For methane, Jähne et al. (1987) is used, whereas for oxygen, Boudreau (1996) is used with the calculation of the dynamic viscosity of water following Matthaus as quoted by Kukulka et al. (1987), Here, r m is the relative soil moisture content, v p the icecorrected volumetric soil porosity, a and w the exponents from Collin and Rasmuson (1988) for air and water, T the soil temperature, p s the soil air or water pressure in atm and k H the Henry constant.All these variables relate to the layer.D a (0,1) is the diffusion coefficient in free air at T 0 = 273.15K and standard pressure p 1 = 1 atm, and D w is the diffusion coefficient in water under the conditions of the layer.The latter two for methane and oxygen are defined as with A and E a from Jähne et al. (1987), and R being the gas constant.T is once more the temperature and µ the dynamic viscosity of water, both of the layer.
To establish the boundary conditions for the system properly, for both the upper and lower boundaries of the soil column, one additional computational point has to be added to the computational system.Also for the boundary conditions, but just for computational reasons, two virtual points at the same distance from the upper or lower boundary as the first or last inner point are needed outside the computational domain.These points have as properties their location and diffusion coefficient only, which are the same as those of the first or last layer.The layer heights are used as weights for the weighted harmonic means of the diffusion coefficients at the borders between the layers.If just boundary points are involved, half of the layer heights are used as weights.

A5 Plant transport -set-up details
The thickness of the exodermis is set to 0.06 mm (Kutzbach et al., 2004).The number of tillers per square metre for rim and centre are given by Kutzbach et al. (2004).The number of tillers per plant is set to one.While the mean accumulated root length of one plant is derived from Shaver and Billings (1975) to be 0.739 m, the root diameter is derived from Kutzbach et al. (2004) to be 1.9 mm.

B1 Modelled relative soil moisture content
The modelled soil moisture content changes seasonally very much.However, because soil water content is restricted to field capacity, there is also a limit for soil moisture content at field capacity.At the rim (Fig. B1a), soil moisture increases in the upper soil part in spring but decreases with the ongoing thawing season.In contrast, at the centre (Fig. B1b), soil moisture increases only slowly in spring, but this increase is ongoing until almost the end of the thawing season.This is due to the greater amount of ice in the soil, which thaws slowly.On the other hand, the greater input of water to the centre than to the rim as soon as there is runoff created at the rim is a continuous additional supply of soil moisture to the centre later in the thawing season.With this, the rim is more moist than the centre in the beginning of the thawing season but drier in the middle and at the end of it (Fig. B1c).Just in the deeper layers, the rim has a little bit more liquid water during the whole thawing season.In winter, however, the amount of liquid water is negligible both at the rim and at the centre.Thus, differences may only be seen in the timing of changes due to thawing or freezing, which both happen earlier at the rim than at the centre.Consequently, they result in earlier wetting of the rim's soil during spring as well as earlier drying of it during freezing.

B2 Modelled relative soil ice content
The modelled soil ice content, in contrast, is almost always higher at the centre than at the rim.Only during freezing in autumn is there a short period when there is more ice in the uppermost soil part at the rim than at the centre.During the thawing season, there is generally very little ice in the upper part of the rim's soil (Fig. B2a), while at the centre, small amounts of ice may also occur in this period (Fig. B2b).Both rim and centre show substantial amounts of ice below 30 cm, even during the summer.Furthermore, during spring, while the uppermost part of the soil at the centre is already thawed, an accumulation of new ice takes place right below, which thaws shortly after.In general, the upper soil part gets its ice thawed and frozen more slowly and later at the centre than at the rim because there is more ice at the centre.Below 30 cm, the difference in ice content between rim and centre increases in summer (Fig. B2c).However, this levels off during freezing until it reestablishes at a lower level in winter.In winter, the soil part with the least amount of ice is not on top, but between 10 and 30 cm both at rim and centre.

B3 Modelled soil temperature
The modelled soil temperatures show deeper thawing and higher temperatures during the thawing season at the rim compared to the centre (Fig. B3a).In addition, rim temperatures reach lower values in winter.Moreover, the thawing season starts earlier and ends later for the rim than for the centre (Fig. B3b).These effects are due to the generally drier soil at the rim compared to the centre.Water dampens the amplitude of the temperature change, and, in addition, the phase change takes up energy.While the warming to 0 • C occurs quickly, the phase change takes time and the soil can only warm further after the phase change is completed.During freezing, the reverse occurs.The cooling then is faster and to lower temperatures at the rim compared to the centre.In general, deeper layers react more slowly and are dampened compared to layers close to the surface.At the rim as well as at the centre, there are short periods with tempera- tures below 0 • C even during summer.The highest temperature differences occur during early spring when there is more ice in the ground at the centre than at the rim.Thus, the rim can reach the zero curtain easier (Fig. B3c).

B4 Modelled oxygen uptake B4.1 Mixed daily sum
The overall pattern of oxygen uptake shows big portions during the early and late thawing season, with a reduced uptake during the mid season (Fig. B4).This is the most moist part of the season, and water effectively reduces oxygen diffusion into the soil.There is also some daily variation in the amount of uptake during the thawing season that is connected to the soil moisture content.The wetter the soil, the less oxygen can enter.Because there is high uptake at the beginning and the end of the thawing season, the overall transport of oxygen is more similar for the rim and the centre, in contrast to methane, where the centre dominates.In winter, no uptake takes place because snow hinders the exchange.

B4.2 Seasonal split
The modelled oxygen uptake at the rim and at the centre is different for the different seasons (Fig. B5).In summer, the uptake is purely positive and greater for the rim than for the centre.Also, the spread of uptake is greater for the rim than for the centre.This is again due to the drier conditions that www.geosci-model-dev.net/10/333/2017/Geosci.Model Dev., 10, 333-358, 2017 allow more diffusion through air, which is quicker and can thus lead to higher uptake compared to diffusion in water or via plants under the wetter conditions at the centre.In winter, the uptake is zero, following the assumption that snow hinders the exchange.In the mixed approach, the overall mean uptake is about 2.21 g O 2 m −2 h −1 .

B4.3 Cumulative sums
At the rim, diffusion delivers a much greater portion of oxygen than plant transport (Fig. B6a).At the centre, both processes provide almost the same amount of oxygen (Fig. B6b).
There are no such pronounced bursts during spring as for methane.While plant transport is smaller than diffusion for both, rim and centre, the difference is much bigger at the rim.At the centre, there is more plant transport but less diffusion than at the rim.Diffusion at the rim and plant transport at the centre are increasing towards the end of the thawing season.In contrast, diffusion at the centre and plant transport at the rim show decreasing contributions towards the end of the thawing season.
In the mixed approach, rim and centre add to a relatively uniform increase in oxygen flux by diffusion over the whole thawing season.For plant transport, the mid season increase is highest, with smaller contributions at the beginning and the end of the thawing season (Fig. B6c).This results from the different timings of high soil moisture content at the rim and at the centre that compensate each other for diffusion.Furthermore, the wetter the soil, the more plant transport relative to diffusion should occur, because the more water, the more diffusion is slowed down.If, moreover, these conditions occur towards the end of the growing season, which is the case at the centre, the effect is bigger than if this happens in spring, which is the case at the rim.Still, diffusion ac-  B1 gives the maximal values.
Table B1.Maximal values of the cumulative sums of modelled oxygen uptake over the modelled time period for rim, centre and a mixed approach of 65 % rim plus 35 % centre for the different transport processes and combined in kg O 2 m −2 , rounded to three non-zero digits.counts for a larger proportion of uptake than plant transport because plant transport was defined as being slower than diffusion in water, while diffusion in air is rather quick.It might still be that the plant transport is too low compared to the total uptake because the root surface might have been chosen too small, like the results for the methane emissions suggest.In total, the rim accounts for more oxygen uptake than the centre (Fig. B6d), but the difference is not as high as for the methane emissions.While the late season is slightly more important at the rim, it is the early season for the centre.
When comparing rim and centre total uptake, diffusion gets reduced to about a third at the centre compared to the rim, and plant transport gets almost 4 times as high (Table B1).This results in a reduction to less than two-thirds of the overall uptake at the centre compared to the rim.While at the rim, diffusion is almost 12 times as high as plant transport, they are almost at the same level at the centre.These differences are again due to the differences in soil moisture content.In the mixed approach, diffusion accounts for about 4.5 times the uptake of plant transport.Overall, 16 kg of oxygen are taken up by each square metre in the course of the modelled time period.

B4.4 Transport process split
Splitting the overall oxygen uptake into the transport processes shows differences in the amount of their contribution per process, depending on location, but also differences in the pattern (Fig. B7a).The uptake is split into different portions between the processes that are more equal for the centre (Fig. B7b), but differ a lot for the rim.There, diffusion is responsible for the majority of the uptake.At the centre, this is only true in the early season and at the freezing.In the mid season, plant transport is much higher than diffusion.While the diffusion part is lower at the centre than at the rim, the opposite is the case for plant transport.In spring, big amounts of oxygen are taken up both at the rim and at the centre.In the late season, some small emissions via diffusion also occur at the centre.In general, uptake through diffusion is greater when soil is drier, which is the case for the rim in the late season and for the centre in the early season.While plant transport is more steady at the rim, there are pronounced peaks at the centre when the soil is wettest.In spring, when the soil is wettest at the rim, plants are not yet that far developed that plant transport could increase to similarly high values as at the centre during the respective times with high soil moisture content.Geosci.Model Dev., 10,2017 www.geosci-model-dev.net/10/333/2017/

Figure 1 .
Figure 1.Modelled water table at rim and centre.Solid lines indicate 1 January; dashed lines indicate 1 April, 1 July and 1 October of the respective year.Only the summer periods are shown, which means less than 5 cm of snow are on the ground.

]Figure 2 .
Figure 2. Modelled methane flux out of soil at the rim, centre and a mixed approach of 65 % rim plus 35 % centre, split into summer and winter.Summer means less than 5 cm of snow are on the ground; winter is the remainder.Because of the wide spread of values, from −0.0747 mg C m −2 h −1 to as high as 86.8 mg C m −2 h −1 , a portion of 4.66 % values was cut to provide a reasonable picture.

Figure 3 .
Figure 3. Modelled methane flux out of soil at (a) the rim, (b) the centre, (c) a mixed approach of 65 % rim plus 35 % centre, split into the different transport processes, and at (d) the rim, the centre and a mixed approach of 65 % rim plus 35 % centre combined, as a cumulative sum.Solid lines indicate 1 January; dashed lines indicate 1 April, 1 July and 1 October of the respective year.Please note the different scales.Table1gives the maximal values.

Figure 4 .
Figure 4. Modelled methane flux out of the soil at the (a) rim and (b) centre, split into the different transport processes.Solid lines indicate 1 January; dashed lines indicate 1 April, 1 July and 1 October of the respective year.Please note the different scales.Because of the wide spread of high values, to as high as 39.3 (a) and 86.6 (b) mg C m −2 h −1 , a portion of 0.108 % (a) and 0.0609 % (b) values was cut to provide reasonable pictures.The minima of the values are −0.0234(a) and −0.158(b) mg C m −2 h −1 .

Figure 5 .
Figure 5. Modelled methane amounts that get produced and oxidised at the (a) rim and (b) centre, split into the different processes.Solid lines indicate 1 January; dashed lines indicate 1 April, 1 July and 1 October of the respective year.Please note the different scales.The maxima of the values are 0.670 (a) and 1.02 (b) mgC m −2 h −1 .

]Figure 6 .
Figure 6.Modelled methane flux out of soil at the rim and centre compared to chamber measurements.Modelled values are only from the summer periods 2003 to 2005, which means less than 5 cm of snow are on the ground.Field measurements took place on 39 days from July to September 2006.Because of the wide spread of high modelled values, to as high as 86.8 mg C m −2 h −1 , a portion of 0.347 % modelled values was cut to provide a reasonable picture.The minimum of the modelled values is −0.0237 mg C m −2 h −1 .

Figure 7 .
Figure 7. Modelled methane flux out of soil in a mixed approach of 65 % rim plus 35 % centre compared to eddy covariance measurements.Light grey background indicates measurement data coverage.X axes indicate the first day of the respective month of the year.Dashed lines indicate 1 July and 1 October of the respective year.Please note the cutouts in between the different years.Because of the wide spread of high modelled values, to as high as 30.4 mg C m −2 h −1 , a portion of 0.0507 % modelled values was cut to provide a reasonable picture.The minimum of the modelled values is −0.0235 mg C m −2 h −1 .

Figure B1 .
Figure B1.Modelled relative soil moisture content of the uppermost metre at the (a) rim and (b) centre as well as (c) the difference centre minus rim in several depths.Solid lines indicate 1 January; dashed lines indicate 1 April, 1 July and 1 October of the respective year.The scale maximum for (a) and (b) is field capacity, ceiled to two digits.

Figure B2 .
Figure B2.Modelled relative soil ice content of the uppermost metre at the (a) rim and (b) centre as well as (c) the difference centre minus rim at several depths.Solid lines indicate 1 January; dashed lines indicate 1 April, 1 July and 1 October of the respective year.The scale maximum for (a) and (b) is the field capacity, ceiled to two digits.The scale for (c) is the same as for the difference of the modelled relative soil moisture content.

Figure B3 .
Figure B3.Modelled soil temperature of the uppermost metre at the (a) rim and (b) centre as well as (c) the difference rim minus centre at several depths.Solid lines indicate 1 January; dashed lines indicate 1 April, 1 July and 1 October of the respective year.

Figure B4 .]Figure B5 .
Figure B4.Modelled oxygen flux into soil in a mixed approach of 65 % rim plus 35 % centre as the daily sum.Solid lines indicate 1 January; dashed lines indicate 1 April, 1 July and 1 October of the respective year.The range of the modelled values is −0.00184 to 87.6 g O 2 m −2 d −1 .

Figure B6 .
Figure B6.Modelled oxygen flux into soil at (a) the rim, (b) the centre, (c) a mixed approach of 65 % rim plus 35 % centre, split into the different transport processes, and at (d) the rim, the centre and a mixed approach of 65 % rim plus 35 % centre combined, as a cumulative sum.Solid lines indicate 1 January; dashed lines indicate 1 April, 1 July and 1 October of the respective year.Please note the different scales.TableB1gives the maximal values.

Figure B7 .
Figure B7.Modelled oxygen flux into soil at the (a) rim and (b) centre, split into the different transport processes.Solid lines indicate 1 January; dashed lines indicate 1 April, 1 July and 1 October of the respective year.Because of the wide spread of high values, to as high as 16.3 (a) and 14.4 (b) g O 2 m −2 h −1 , a portion of 0.0254 % (a) and 0.0178 % (b) values was cut to provide reasonable pictures.The minima of the values are −0.00185(a) and −0.136(b) g O 2 m −2 h −1 .

Table 2 .
Change in the cumulative methane emissions over the modelled time period in %, when the parameter was modified by ±10 % compared to its default setting.

Table 3 .
Modelled daily methane flux for the summer periods 2003 to 2005 for the rim, centre and a mixed approach of 65 % rim plus 35 % centre in mg CH 4 m −2 d −1 , rounded to three non-zero digits.Summer means less than 5 cm of snow are on the ground.Please note the different unit here.

Table 4 .
Modelled hourly methane flux for the summer periods 2003 to 2005 for the rim, centre and a mixed approach of 65 % rim plus 35 % centre in mg C m −2 h −1 , rounded to three non-zero digits.Summer means less than 5 cm of snow are on the ground.