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

A consistent, 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. Oxygen has been explicitly incorporated in diffusion, transport by plants and two oxidation processes, of which one uses soil oxygen, while the other uses oxygen that is available 5 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 polygonal tundra site Samoylov, Lena 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 10 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 center. The evaluation shows sufficiently good agreement with field observations despite the fact that the module has not been specifically calibrated to these data.

1 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 to rise slowly but steadily since the onset of industrialisation, and, after a hiatus at the beginning of the 21st century, have recently be found 25 to rise again. These recent dynamics in the global atmospheric methane budget are still not fully explained, emphasising the fact that also future trajectories of methane and its role in global climate change are 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 30 dioxide, it makes up for about 20 % of the radiative 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 35 often characterised by small-scale mosaics of wet and dry surfaces. 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-oriented 40 model with a bottom-up approach for methane balance estimation. However, process-based modelling approaches are powerful tools that help to quantify recent and future methane fluxes at large spatial scale and over long time periods in such remote areas. They can give first estimates 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, hence feedback mech-45 anisms, 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.

Scheme for Biosphere Atmosphere Coupling in Hamburg) of the MPI-ESM (Max Planck Institute
Earth System Model)was chosen for this work. The starting point was a model version that has a carbon balance (Reick et al., 2013), a five layer hydrology (Hagemann and Stacke, 2015) and includes permafrost as described in Ekici et al. (2014). A parallel development by Schuldt et al. (2013) 75 incorporated wetland carbon cycle dynamics and was also integrated in the model version presented in this work. The basis for the methane-related processes were the works by Walter and Heimann (2000) and Wania et al. (2010). Special focus was also put 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 80 polygonal tundra site in Yakutia ::::::: northeast :::::: Siberia.

Site description
For the purpose of evaluation, this model has been applied at the site Samoylov Island, located 120 km south of the Arctic Ocean in the Lena River Delta in Yakutia with an elevation of 10 to 16 m 85 above sea level. The mesorelief of Samoylov Island is flat, while as microrelief, there are low-center polygons with the soil surface about 0.5 m higher at the rim than at the center. 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 centers is at about 0.5 m (Boike et al., 2013). While the water table at the polygonal rims is generally well below the soil 90 surface, the polygonal centers are often water saturated with water tables at or above the soil surface . 3 The vegetation on Samoylov Island 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 95 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 center, is Carex aquatilis but with dominance of only 8 % at the rim compared to 25 % at the center. However, most of the species present at the rim are different from those present at the center. According to Sachs et al. (2010), the proportions of moist and wet microsites are approximately 65 % moist and 35 % 100 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 center.

Layer structure
For a numerically stable representation of gas transport processes in soils, a much finer vertical soil 105 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 (in the order of a few cm) but still inherits the hydrological and thermal information contained in the coarse scheme. Number and height of layers can be chosen arbitrarily, allowing also non-equidistant solutions. 110 Internally, the module uses midpoints and lower boundaries of the layers as well as distances between the midpoints. At the bottom, the layering scheme is truncated at depth to bedrock. The layers, where the plant roots end, i.e., rooting depth lies,

115
water table lies and 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 App. A1).

120
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, also the previous day's mean soil temperature is 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 coarse layer values. The relative soil water 130 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 ::: ice-free pores can be defined, which is calculated by division of the relative soil moisture content by the ice-corrected volumetric 135 soil porosity. Finally, the relative air content of the ice free pores is defined as one minus the relative moisture content of the ice free pores.

185
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 :::::: ice-free pore volume changes. Thus, concentration values also change, but 190 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 195 on the soil hydrological conditions, carbon dioxide or methane are produced from the decompos-6 ing 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 200 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 App. A3).

205
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 stratum 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 210 is set assuming all decomposed carbon above the water table and half of it below the water table gets 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.

215
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 220 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. 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).

225
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 230 oxygen and K O2 m = 2 mol m 3 , which is chosen to be the average concentration of oxygen at the 7 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 235 calculated as the sums of the respective amounts.

Ebullition
The implementation of the ebullition of methane follows largely 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 240 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 245 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 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, 250 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 proposed in Wania (2007), these two steps are iterated until stable state conditions are reached.

255
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 260 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 265 above the water table. The former is calculated after 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 270 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,

275
with k H CH4 being the Henry solubility, c CH4 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 of 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 280 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 285 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 290 concentrations. For soil as a porous medium, moreover with changing pore volumes because of different contents of ice, the ice-corrected 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 295 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 hour, 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 300 diffusion coefficients.

9
As 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.

305
If there are at least 5 cm of snow, no flux is assumed, and therefore Neumann type is applied also at the top. However, if there are less than 5 cm of snow, ambient air conditions are assumed to hold at the boundary, and therefore Dirichlet type with gas concentration in free air is applied, Here, v p is the volumetric soil porosity, c denotes the gas concentration, t is the time, x is the depth,

310
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 App. 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.

315
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 with their respective pore volumes as usual and summing the resulting amounts over the layers gives the total fluxes of methane and oxygen.

320
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 simulated :::::::: 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 325 efficient barrier against gas exchange.
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.,340 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 · ⇡, 345 multiplied by plant density, q p = t ph tp . 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 LAI, using t ph = max(t m ) · LAI max(LAI) , with t m being the number of tillers per square metre. Please see also App. A5.

350
The root tissue is assumed to be distributed equally between all root-containing layers, thus A rl r = A ges r · h P 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 O2 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.

360
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 CarexA.
dom V ascularP. . 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 365 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 center, 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 370 The variables definitions are the same as for oxygen and c CH4 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, 375 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, 380 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 center, 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 385 roots in the layer that are able to transport gases, and c O2 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.

410
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 415 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 simulated ::::: ruled ::: out ::: that ::::::::: modelled methane processes may be biased ::::::: modified : through the chosen order under certain conditions.

495
For additional results concerning modelled physical conditions, such as soil moisture and ice content as well as soil temperatures, the reader is referred to App. B1 to B3.

Modelled methane flux in summer and winter
The modelled methane fluxes at rim and center 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.

500
The spread of the flux is greater for the center than for the rim in both summer and winter. While the majority of flux values in summer is positive at the center, 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 center, there are some rare events when uptake takes place. In the mixed approach, which means 65 % rim and 35 % center, the overall mean emission is about 0.0813 mgC m 2 h 1 505 for the summer period. The overall higher emissions at the center are due to higher moisture and thus more favourable conditions for methane production in concert with lower methane oxidation rates.
This is due to the drier conditions at the rim that allow a fast diffusion through air, while ebullition is only possible with a minimum of soil moisture ::::: under :::::::: conditions ::: of :::: high ::::: water :::::: content. Because in the model, higher soil moisture is calculated from the middle to the end of the thawing season, most of the emissions by ebullition and plant transport :::: occur : at the center occur then (Fig. 3B).

530
In the mixed approach, only the diffusion of the rim alters :::::::::: substantially : the pattern of the emissions substantially (Fig. 3C). In total, the polygon center 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 535 center. Emissions at the rim are highest during spring, while they are highest at the center during the mid and late season (Fig. 3D).
When comparing the total fluxes of the center 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 540 seven times higher total methane emissions at the center than at the rim. While diffusion at the rim :: At ::: the ::: rim :::::::: diffusion is more than 13 times as high as plant transportat the rim, the diffusion : , ::::: while at the center : it : is just slightly higher than the plant transportthere :::: plant :::::::: transport. Ebullition is about 4.5 times as high as plant transport both at the rim and at the center. These differences are again due to the differences in soil moisture content, which allow more production under higher soil moisture and 545 thus also lead :::: leads 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 , but :: of ::: gas, :::: and under wetter conditions , plant transport may increase relative to diffusion. With reduced soil air, the remaining velocity of the diffusion is almost at the same order of magnitude like ::: than : the overall 550 velocity of plant transport, in contrast to the velocity of diffusion mainly through air.
Still it seems, that the plant transport in the model is too low compared to the total flux. While the diffusion flux to the atmosphere only happens at the soil surface, the surface area of the gas transporting roots is the relevant boundary for plant transport. The value of this is not well-known, 555 so the module might need further adjusting of parameters connected to plant root surface area to improve the share of plant transport. Furthermore, ebullition needs substantial amounts of soil moisture, and this is more common at the center than at the rim. Consequently, substantially more ebullition is found at the center 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,

Production versus oxidation
Methane oxidation follows the pattern of methane production as long as enough oxygen is available 600 (Fig. 5A). Production, and hence also oxidation, is higher during times of more moist conditions for both, the rim and the center, and also higher for the center than for the rim (Fig. 5B). At the center, 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 605 rim and at the center 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 center and their dynamics determine these results.

Comparison to chamber measurements
Although the number of available field data is small and from a different year than the meteorologi-

Comparison to eddy measurements
Eddy covariance data had the best available data coverage of field measurements (light grey areas 635 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 mgC m 2 h 1 with mean 0.609 mgC m 2 h 1 . The range of 640 modelled summer methane emissions in this time frame is -0.023 to 30.4 mgC m 2 h 1 with mean 0.0813 mgC 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 shows higher extremes.

655
Still, Fig. 7 also shows some patterns that are present in both model results and observations, e.g.
For additional results concerning modelled oxygen uptake, such as mixed daily sum, seasonally 670 split and cumulative sums as well as transport process split, see App. B4.

755
Moreover, the modelled soil climate and hydrology are not the same as those to be found in the field.
The model application for remote permafrost areas may also be limited by the availability of longterm and complete observations of meteorological data to be used as model forcing. Forcing data 770 and methane fluxes are required for the same time period, which optimally lasts over one 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-oriented ::::::::: data-driven : upscaling of eddy covariance or chamber based observations (e.g. Christensen et al., 1995;Marushchak et al., 2015).
Despite our mean values are located towards the lower end, our minimum, mean and maximum values fit well within the given range, that shows a widespread of possible observations. Bartlett et al.

Conclusions
The aim of this study was to develop a consistent, :::: 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 land surface scheme JSBACH of the global Earth system model MPI-ESM has been enhanced for this purpose. The 950 new methane module : of :::::::::::::::: JSABCH-methane represents methane production, oxidation and transport.

960
When combined with a module for oversaturated :::::::::::: water-saturated : soil conditions like TOPMODEL (e.g. Kleinen et al., 2012), such 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-oriented modelling can 965 contribute to understand 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 highly integrated also 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 970 with site level field measurements, sufficiently good agreements could be shown, despite the module has not been adjusted to site specific processes or features. Coupling such 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. The model code used in this work is available upon request for academic and non-commercial use. 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 similar to Schuldt et al. (2013) though appropriate tem- Though the rate of organic matter decomposition at the evaluation site is not known, the presentday amount of carbon in the soil is known (Sect. 2.2.3). Considering short time scales only, the 1000 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 1005 of the soil profile by Zubrzycki et al. (2013). The values of the six measurement depths were averaged over the sixteen different center respectively six rim cores. These resulting averages have been After 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 respective 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 1035 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 pres-1040 sure. They are calculated after Massman (1998) from values at the soil surface with over 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 soil surface are calculated from values at 0 C and 1 atm (Massman, 1998).

1045
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 after Matthaus as quoted by Kukulka et al. (1987), Here, r m is the relative soil moisture content, v p the ice-corrected 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 respective water pressure in atm and k H the Henry constant, all of 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 1055 oxygen are defined as D a CH4 (0,1) = 1.952 · 10 5 m 2 s 1 , 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.

1060
To establish the boundary conditions for the system properly, for both the upper and lower boundary 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 in the same distance from the upper respective lower boundary as the first respective last inner point are needed.
These points have as properties their location and diffusion coefficient only, which are the same as 1065 those of the first respective last layer. The layer heights are used as weights for the weighted harmonic means of the diffusion coefficients at the borders between the layers. Just if boundary points are involved, half of the layer heights are used as weights.

A5 Plant transport -setup details
The thickness of the exodermis is set to 0.06 mm (Kutzbach et al., 2004). The number of tillers per 1070 square metre for rim and center 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. 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. 8A), soil moisture increases in the upper soil part in spring but decreases with ongoing thawing season. In contrast, at the center (Fig. 8B), 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 1080 greater amount of ice in the soil, which thaws slowly. On the other hand, the greater input of water to the center than to the rim as soon as there is runoff created at the rim is a continuous additional supply of soil moisture to the center later in the thawing season. With this, the rim is more moist than the center in the beginning of the thawing season but drier in the middle and at the end of it (Fig.   8C). Just in the deeper layers, rim has a little bit more liquid water during the whole thawing season.

1085
In winter, however, the amount of liquid water is negligible both at the rim and at the center. Thus, differences may only be seen in the timing of changes due to thawing respective freezing, which both happen earlier at the rim than at the center. Consequently, they result in earlier wetting of the rim's soil during spring as well as earlier drying of it during freezing.

1090
The modelled soil ice content, in contrast, is almost always higher at the center than at the rim.
Only during freezing in autumn, there is a short period when there is more ice in the uppermost soil part at the rim than at the center. During the thawing season there generally is very little ice in the upper part of the rim's soil (Fig. 9A), while at the center, small amounts of ice may also occur in this period (Fig. 9B). Both, rim and center, show substantial amounts of ice below 30 cm even 1095 during the summer. Furthermore, during spring, while the uppermost part of the soil at the center 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 center than at the rim because there is more ice at the center. Below 30 cm, the difference in ice content between rim and center is increasing in summer (Fig. 9C). However, this levels off during freezing, until it 1100 reestablishes in winter at a lower level. 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 center.

B3 Modelled soil temperature
The modelled soil temperatures show deeper thawing and higher temperatures during the thawing season at the rim compared to the center (Fig. 10A). In addition, rim temperatures reach lower 1105 values in winter. Moreover, the thawing season starts earlier and ends later for the rim than for the center (Fig. 10B). These effects are due to the generally drier soil at the rim compared to the center.
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 1110 then is faster and to lower temperatures at the rim compared to the center. In general, deeper layers react more slowly and dampened compared to layers close to the surface. At the rim as well as at the center, there are short periods with temperatures below 0 C even during summer. The highest temperature differences occur during early spring when there is more ice in the ground at the center than at the rim. Thus, the rim can reach the zero curtain easier (Fig. 10C).

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. 11). 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 1120 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 center, in contrast to methane, where the center is dominating. 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 center is different for the different seasons (Fig.   12). In summer, the uptake is purely positive and greater for the rim than for the center. Also, the spread of uptake is greater for the rim than for the center. This is again due to the drier conditions that allow more diffusion through air, which is quicker and can thus lead to higher uptake compared 1130 to diffusion in water or via plants under the wetter conditions at the center. 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 gO 2 m 2 h 1 .

B4.3 Cumulative sums
At the rim, diffusion delivers a much greater portion of oxygen than plant transport (Fig. 13A). At 1135 the center, both processes provide almost the same amount of oxygen (Fig. 13B). There are no such pronounced bursts during spring as for methane. While plant transport is smaller than diffusion for both, rim and center, the difference is much bigger at the rim. At the center, there is more plant transport but less diffusion than at the rim. Diffusion at the rim and plant transport at the center are increasing towards the end of the thawing season. In contrast, diffusion at the center and plant 1140 transport at the rim show decreasing contributions towards the end of the thawing season.
In the mixed approach, rim and center add to a relatively uniform increase of 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. 13C). This results from the different timing of high soil moisture content at the rim and at the center that compensate each other for diffusion. Furthermore, the wetter the soil, the more plant transport relatively to diffusion should occur, because the more water the more is diffusion slowed down. If, moreover, these conditions occur towards the end of the growing season, which is the case at the center, the effect is bigger than if this happens in spring, which is the case at the rim. Still, diffusion accounts for 1150 a larger proportion of uptake than plant transport because plant transport was defined to be 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 center (Fig. 13D), but the difference is not as high as for the methane emissions. While the 1155 late season is slightly more important at the rim, it is the early season for the center.
When comparing rim and center total uptake, diffusion gets reduced to about a third at the center compared to the rim, and plant transport gets almost 4 times as high (Table 5). This results in a reduction to less than two-thirds of the overall uptake at the center compared to the rim. While at 1160 the rim, diffusion is almost 12 times as high as plant transport, they are almost at the same level at the center. These differences are again due to the differences in soil moisture content. In the mixed approach, diffusion accounts for about 4.5 times of 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.

1165
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. 14A).
The uptake is split into different portions between the processes, that are more equal for the center ( Fig. 14B) but differ a lot for the rim. There, diffusion is responsible for the majority of the uptake. At the center, this is only true in the early season and at the freezing. In the mid season, plant transport 1170 is much higher than diffusion. While the diffusion part is lower at the center 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 center. In the late season, also some small emissions via diffusion occur at the center.
In general, uptake through diffusion is greater when soil is drier, which is the case for the rim in the late and for the center in the early season. While plant transport is more steady at the rim, there are 1175 pronounced peaks at the center 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 center during the respective times with high soil moisture content.