Calibrating a wetland methane emission model with hierarchical modeling and adaptive MCMC

Methane (CH4) emission estimation for natural wetlands is complex and the estimates contain large uncertainties. The models used for the task are typically heavily parametrized and the parameter values are not well known. In this study we perform a Bayesian model calibration for a new wetland CH4 model to improve quality of the predictions and to understand the limitations of such models. The detailed process model that we analyze contains descriptions for CH4 production from anaerobic respiration, CH4 5 oxidation, and gas transportation by diffusion, ebullition, and the aerenchyma cells of vascular plants. The processes are controlled by several tunable parameters. We use a hierarchical statistical model to describe the parameters and obtain the posterior distributions of the parameters and uncertainties in the processes with adaptive MCMC techniques. For the estimation, the analysis utilizes measurement data from the Siikaneva flux measurement site in Southern Finland. The model parameters are calibrated using six different modeled peat column depths, and the hierarchical modeling allows 10 us to assess the effect of the parameters on an annual basis. The results of the calibration and their cross validation suggest that the early spring net primary production and soil temperatures could be used to predict the annual methane emissions. The modeled peat column depth has an effect on how much the plant transport pathway dominates the gas transport, and the optimization moved most of the gas transport from the diffusive pathway to plant transport. This is in line with other research, highlighting the usefulness of algorithmic calibration of biogeochemical models. 15 Modeling only 70 cm of the peat column gives the best flux estimates at the flux measurement site, while the estimates are worse for a column deeper than one meter or shallower than 50 cm. The posterior parameter distributions depend on the modeled peat depth. At the process level, the flux measurement data is able to constrain CH4 production and gas transport processes, but for CH4 oxidation, which is an important constituent of the total CH4 emission, the determining parameter is not identifiable. 20 1 Geosci. Model Dev. Discuss., doi:10.5194/gmd-2017-66, 2017 Manuscript under review for journal Geosci. Model Dev. Discussion started: 21 April 2017 c © Author(s) 2017. CC-BY 3.0 License.


Introduction
Methane is the second most important gas in the atmosphere in terms of its capacity to warm the climate, currently with the radiative forcing power of 0.97 Wm −2 .This is a sizable part of the total effect of well-mixed greenhouse gases, which is approximately 3.0 Wm −2 .According to IPCC (2013), the amount of CH 4 in the atmosphere has risen to its highest level in at least the last 800000 years due to human activity, and based on ice core measurements, also its growth rate is presently very likely at its highest level in the last 22000 years.
The sources of CH 4 are both anthropogenic and natural.In years 2003-2012, 60% of the global emissions were anthropogenic (range 50-65 %) and about one third came from natural wetlands.The most important source of uncertainty in the global methane budget is attributable to emissions from wetlands and other inland waters.Combining top-down and bottom-up estimates, natural wetland emissions range from 127 to 227 Tg CH 4 yr −1 (Saunois et al., 2016).Anthropogenic sources include rice paddies, landfills, enteric fermentation and manure, incomplete combustion of hydrocarbons, and natural gas leaks (Ciais et al., 2013).
The methane from wetlands is produced by prokaryotic archaea under anaerobic conditions.The main sink for atmospheric CH 4 is its oxidation in troposphere by OH ions and the average lifetime of a CH 4 molecule in the atmosphere is 9.1 ± 0.9 years (Prather et al., 2012;IPCC, 2013).
The wetlands in the boreal zone are a significant contributor to the total CH 4 emissions from wetlands (Kirschke et al., 2013), and for this reason the CH 4 emissions from them have been intensively studied, also with models, during the past years (Wania et al., 2010;Kaiser et al., 2016;Petrescu et al., 2015).
The need for improved wetland methane emission modeling is amplified by the fact that although annual mean precipitation is projected to increase in the boreal zone (Ruosteenoja et al., 2016), changes in the frequency and duration of severe drought may follow an alternate path (Lehtonen et al., 2014), manifesting the need to study wetland responses to extreme events.
Changes to hydrological conditions such as draining or recurring low water table depth can alter the balance of greenhouse gas emissions (Frolking et al., 2011;Petrescu et al., 2015).Modeling and calibrating for such exceptional events can be difficult, as was found for instance by Mäkelä et al. (2016).
The HelsinkI Model of MEthane buiLd-up and emIssion for peatlands (HIMMELI) is a relatively full-featured wetland/peatland CH 4 emission model and it is described in detail in Raivonen et al. (2017).The model contains process descriptions for CH 4 production from anaerobic respiration, O 2 consumption and CO 2 production from oxic respiration, and gas transport processes via diffusion, ebullition, and plant transport.Modeling the concentrations of CH 4 , O 2 , and CO 2 in the peat column is explicitly included.The peat column depth can be set at any desired value, and the water table movement determines the part of the peat column that is favorable for CH 4 production.The version of HIMMELI in this work has additional processes, described in Sec.3.1, and the modified model is referred to as sqHIMMELI (square root HIMMELI), as it contains a description of CH 4 production from root exudates.
Computer models describing environmental processes accumulate error at many levels (Sanso et al., 2007).The sources include time-and space discretization, incomplete physics and biochemistry descriptions, insufficient information about the does not accurately and appropriately describe the observations.On the other hand, with different parameters for all the years, the parameters are easily overfitted, meaning that while the resulting model fits the data well, it does not accurately predict future fluxes (Gelman et al., 2013).Hierarchical modeling provides a solution for these problems.
In the present study, the sqHIMMELI model is calibrated using adaptive Markov chain Monte Carlo (MCMC) techniques to evaluate a hierarchical statistical model for the model parameters.The calibration is done for the boreal Siikaneva site.This study complements the work in Raivonen et al. (2017) in describing the effects of various parameters on the processes and fluxes, and analyzing what kinds of configurations best describe the studied boreal wetland.
Classical optimization is often misleading due to the multi-modality of the objective function, as for example Müller et al. (2015) reported in a study where they used a surrogate model to calibrate the parameters of the CH 4 model component of the Community Land Model.This multi-modality can be accommodated for by using MCMC techniques.Utilizing MCMC methods for optimizing environmental models and studying their uncertainties is not new (Laine, 2008;Ricciuto et al., 2008;Hararuk et al., 2014), but to our knowledge they have not been used for wetland CH 4 model parameter estimation before.
Moreover, the research that the authors are aware of does not investigate the interannual variability of parameters, as is done in this study.
The main objective of this work is to analyze the capabilities and limitations of a modern featureful wetland CH 4 model by looking into the shape of the posterior parameter distributions, parameter correlations, and the roles, identifiabilities, interdependencies, and interconnections of the parameters and the processes they control.The simulations and the analyses are performed with six different modeled peat depths, which allows for assessing how the modeled peat column depth affects the model behavior, and how deep a peat column is optimal based on the flux measurement data used.
2 Siikaneva wetland flux measurement site and model input data Methane flux measurements were needed for estimating the model parameters, and for that purpose observational data from the Siikaneva peatland flux measurement site in southern Finland (61 • 50'N, 24 • 12'E) were used.The site is a boreal oligotrophic fen with a peat depth of up to four meters.The data collection started in 2005, and in this work eddy covariance (EC) CH 4 flux measurements from years 2005 to 2014 were used.In the current application of the EC method, the flux was calculated from the wind speed and direction, and CH 4 concentration information, both of whose sampling frequency was 10 Hz.During the measurement period several different instruments were used for methane concentration measurements: Campbell TGA-100 (2005)(2006)(2007) and 04/2010-08/2010), Los Gatos RMT-200 (01/2008RMT-200 (01/ -02/2014)), Picarro G1301-f (04/2010-10/2011) and Los Gatos FGGA (2014).The wind velocity vector was analyzed by a USA-1 acoustic anemometer by METEK (Rinne et al., 2007).All the EC-data were post-processed in a consistent manner using an in-house software EddyUH (Mammarella et al., 2016).
The required inputs for sqHIMMELI are daily soil temperatures, water table depths (WTD), net primary production (NPP), and leaf area indexes (LAI).The soil temperature profile for the grid used was generated by interpolating between the measurement depths (-5 cm, -10 cm, -20 cm, -35 cm and -50 cm) and assuming that at -3 meters and below the temperature is a constant +7°C.This was the mean temperature of all the years at -50 cm depth.The WTD data used was available as measurement data, and where data was missing, it was gap-filled by repeating the previous measured value.Net primary production cannot be measured in a direct way, and hence modeled values for it were used.Also for LAI, a simple model was used for obtaining the input.For more details, see Appendix E. A summary of the data used is given in Table 1.

The sqHIMMELI model
The HIMMELI (HelsinkI Model of MEthane buiLd-up and emIssion for peatlands) model (Raivonen et al., 2017) is a detailed model for estimating CH 4 emissions from wetlands.It was developed at the University of Helsinki in collaboration with the Finnish Meteorological Institute and the Max Planck Institute for Meteorology in Hamburg.The model is designed to be used as a submodel for different larger modeling environments, such as regional and global biosphere models.It contains processes describing the production of CH 4 and CO 2 , the loss of CH 4 and O 2 , and transport of CH 4 , O 2 , and CO 2 between the soil and the atmosphere.The CH 4 transport can take place by diffusion in peat (in water and in the air), by ebullition (transport by bubble formation), and by diffusion in the porous aerenchyma tissues in vascular plants.The model is driven with peat temperature, WTD and LAI of the aerenchymatous plants.The process descriptions are mainly adopted from previous wetland CH 4 models such as Arah and Stephen (1998), Wania et al. (2010) and Tang et al. (2010).The version of the model used here differs slightly from that presented in (Raivonen et al., 2017), and is therefore called with the different name of sqHIMMELI to avoid confusion.
The model simulates the processes in a discretized peat column.The number and thickness of the peat layers can be varied, but in this work a variable number of 10 cm layers is used, similarly to e.g.Kaiser et al. (2016).Effectively, the total depth of the peat column changes, not the thickness of the layers.The water table divides the column into water-filled and air-filled parts, and CH 4 is produced only in the inundated anoxic layers.In the present configuration, the NPP-related CH 4 production is allocated into the layers according to the vertical distribution of the root mass, described in Sect.3.2.The internal time resolution of the model is dynamically adjusted depending on the model state, and the output interval is set to one day.HIMMELI itself, as presented in Raivonen et al. (2017), does not simulate carbon uptake (photosynthesis) or peat carbon pools but instead it takes as input the rate of anoxic respiration.The differences between HIMMELI and sqHIMMELI are described below in Sec.3.1 and 3.2 and in Appendix A3.
For each modeled process in sqHIMMELI, there are parameters regulating the process, affecting the concentrations of CH 4 , O 2 and CO 2 in the peat column, and the wetland methane emissions.The equations describing the physics relevant to the optimized parameters are listed in section 3.4.Other relevant model equations are listed in Appendix A.

Root exudates and peat decomposition
Methanogens prefer recently assimilated fresh carbon as their energy source, for instance the root exudates of vascular plants (Joabsson and Christensen, 2001).A connection between ecosystem productivity and CH 4 emission has been observed in several wetland studies (Bellisario et al., 1999;Whiting and Chanton, 1993).However, anoxic decomposition of litter and Geosci.Model Dev. Discuss., doi:10.5194/gmd-2017-66, 2017 Manuscript under review for journal Geosci.Model Dev. Discussion started: 21 April 2017 c Author(s) 2017.CC-BY 3.0 License.older peat also produces CH 4 (Hornibrook et al., 1997).Many models form CH 4 substrates by extracting directly a fraction of the net primary production (van Huissteden et al., 2009;Wania et al., 2010), and some rely on heterotrophic peat respiration only (Riley et al., 2011).In sqHIMMELI both primary production and anaerobic peat decomposition were included.
The modified sqHIMMELI model contains an exudate pool description, from which it produces methane via Eq. 3 and A5.
The exudate pool itself is described by Eq. 4, detailing how the modeled NPP turns into root exudates.Effectively, a fraction of NPP determined by the parameter ζ exu produces root exudates, which are then distributed as anaerobic respiration according to the root distribution into the peat column at the rate determined by the model parameter τ exu .The part ending up under the water table produces CH 4 and CO 2 , depending on the oxygen content of the water, and above the water table the exudates are respired into CO 2 .
The second source of anaerobic respiration, the anaerobic peat decomposition, is modeled in sqHIMMELI with a simple Q 10 -model adopted from Schuldt et al. (2013).The peat under the water table is prescribed a turnover time, based on which anaerobic respiration and CH 4 are produced according to Eq. 5 and A6.

Root distributions
The sqHIMMELI model differs from HIMMELI in the details regarding the root distribution model.Compared to measurement data of root distributions of aerenchymatous sedges from Saarinen (1996), the original root distribution π(z), adopted from Wania et al. (2010) and described by does not describe the distribution of roots well.Here z is depth, and λ root is a parameter describing the steepness of the decaying exponential curve.This formula is replaced with With the Gaussian shape, the new root density decreases faster with depth.Without this change, the optimization process calibrates the model to have very high root masses below 50 cm underground.The other difference between the models is that in the original model there are vanishingly few roots below the depth of one meter, but according to Saarinen (1996), sedge roots can reach to as low as 2.3 m under the surface.The term C 1 in Eq. 2 was added to remedy this.
Before starting the optimization, the parameters C 0 , C 1 , and z 0 were fitted to data from Saarinen (1996), resulting in values of C 0 = 215, C 1 = 6, and z 0 = 0.105.The different root distributions are shown in Fig. 1.

Peat depth
Methane is produced from anaerobic peat decomposition at all peat depths in the sqHIMMELI model, and its transport and oxidation affect the modeled CH 4 emission.The homogeneous model description of the peat column is highly idealized, as in reality the peat column varies from place to place with respect to CH 4 production rate, production depth, and gas transport.
Increasing peat depth in the model is a liability, since the deeper the column, the more expensive the model is to run (see Sect.F2).The model calibration is run for the peat depths of 40, 50, 70, 100, 150, and 200 cm in order to find the optimal peat column depth for the model.

Parameter descriptions for sqHIMMELI
The parameters for the optimization were chosen to constrain the processes most important for the CH 4 emission.Of the optimized parameters, all but ζ exu and Q 10 are constant for all years.However, ζ exu and Q 10 change year to year to reflect the changes in the relative CH 4 input to the system from peat decomposition and NPP-based production.This will allow to analyze the year to year changes in relative importances of the production pathways.The setup is natural, as for example Bergman et al. (2000) report the Q 10 -values changing from measurement date to another, even within a single year.
The parameters and their physical meanings are CH 4 production-related parameters 1. τ exu : Controls the decay rate of exudates, ν, from the root exudate pool P exu , 2. ζ exu : Fraction of NPP carbon that goes to the root exudate pool.
where ψ t is the rate of NPP at time t, P exu is size of the root exudate pool, and ν was given by Eq. 3.
3. τ cato : Controls the base rate of peat decomposition into CH 4 in Eq. 5.
4. Q 10 : Controls the temperature dependence of the rate of peat decomposition into CH 4 in anaerobic conditions via factor k cato , given by the equation Oxidation and respiration parameters 5. V R0 : Respiration parameter controlling the rate of heterotrophic respiration, which consumes O 2 and produces CO 2 .
This affects the rate of temperature dependent heterotrophic respiration, V R (z), given by Here ∆E R is a parameter affecting the temperature dependence of the heterotrophic respiration, R is the universal gas constant, and T (z) is temperature at depth z.
Here ∆E oxid is a parameter affecting CH 4 oxidation that is not part of the optimization.
Gas transport-related parameters 7. λ root : Controls how the root mass is distributed.See Eq. 2.
9. τ : Root tortuousity parameter affecting the root conductance K R .A tortuousity of 1 means that the roots are not decreasing the conductance via their curvedness.The equation for the conductance is where π(z) is the root mass density as a function of depth, over which the sum of the density is one, and m is the total root mass per square meter, set to be proportional to LAI.
10. f D,a : Fraction of the diffusion rate in air-filled peat divided by the diffusion rate in free air.The parameter affects the diffusion flux in the model: the higher this parameter is, the more there is diffusion as it takes a shorter time for the CH 4 to exit the peat reducing the possibility of oxidation and increasing the concentration gradient driving diffusion.The equation is where D air is the diffusion rate in air-filled peat, D 273 air is the diffusion base rate at 273K, and T is the temperature.This parameter is also present in Eq. 8.
11. f D,w : Same as above, but in water.The equation describing the peat-water diffusion rate is where the terms are analogous to the ones in the previous equation.

Model calibration and MCMC
The model calibration consisted of two steps: optimization and MCMC.Both of these steps were run separately for each different number of peat soil layers assessed (4, 5, 7, 10, 15, and 20 layers, each layer corresponding to 10 cm of peat).In the following, an experiment refers to one of these MCMC runs.

Calibration algorithms
For optimization, an initial parameter vector was first drawn from the prior, and the parameters were then optimized against the costfunction described in Eq. 11.The algorithm used was the simplex-based BOBYQA, described in Powell (2009).In our tests, it was significantly faster to converge than NEWUOA (Powell, 2004), L-BFGS or Nelder-Mead (Nelder and Mead, 1965).For each experiment, the model was optimized by running 350 model simulations with the minimization algorithm, which was enough for finding a local minimum to start the MCMC sampling from.
At the points obtained in the optimization, the model was linearized and from the Jacobian a suitable initial proposal covariance matrix for MCMC was estimated.After this the MCMC sampling was performed to estimate the posterior distribution.
In order to be able to assess the annual parameter and CH 4 transport pathway changes, a hierarchical description for two of the parameters was used.These parameters were Q 10 controlling the temperature dependence of the peat decomposition rate, and ζ exu , regulating the production of root exudates from NPP.
The hyperparameters are the means and variances defining the Gaussian priors of the hierarchical parameters Q 10 and ζ exu .
They were updated using fixed Gaussian hyperpriors with Gibbs sampling.The sampling distribution depends on the current values of the hyperparameters.The role of the hyperprior is to constrain the distribution from which the hyperparameters are sampled.
Technically, a Metropolis-within-Gibbs-method Gelman et al. (2013) for sampling the hierarchical parameters, non-hierarchical parameters, and the hyperparameters was used, presented briefly in Appendix D. The model parameters (i.e.everything except the hyperparameters) were sampled with the Adaptive Metropolis (AM) MCMC algorithm (Haario et al., 2001), which uses a Gaussian proposal distribution, whose covariance matrix is adapted as the chain evolves, and over time the acceptance rate gets closer to an optimal value, which is 0.23 for Gaussian targets in large dimensions (Roberts et al., 1997).If the algorithm proposes values outside the hard parameter limits listed in table 2, the model will not be evaluated and the value is rejected.
Our empirical data for the hierarchical model were the nine years from 2006 to 2014, meaning that for each of these years there were corresponding ζ exu and Q 10 parameters in the optimization.The model needed to be spun up for each annual flux estimation in order to have a realistic column of gas concentrations available.For this reason, the previous year was always also simulated, and for the likelihood only the residuals from the latter year were included in the calculations.Therefore year 2005 did not contribute directly to the values of the objective function.The different years were run in parallel to save execution time.

Objective function
As in many practical MCMC applications, a major part of the parameter estimation problem is the proper definition of the objective function.It is defined here based on a priori information about the measurement uncertainties, based on information from the model residuals, and based on annual flux estimates.Additionally prior information about the parameter values is also utilized.

Model residuals and error model
The second component of the objective function contains the daily CH 4 flux measurements.It is assumed that the daily flux estimate uncertainty is dependent on a fraction α of the flux measurement (Richardson et al., 2006) and some constant error, γ (e.g.measurement device precision).The model error is expected to follow a similar form, and hence α and γ contain the contributions from both the model and measurement errors.
The combined error is described by a Laplace distribution.The flux observations are reported to follow a distribution of this type, rather than a Gaussian distribution (Richardson et al., 2006).
When determining the parameters γ and α, the resulting residuals end up being autocorrelated.Therefore they are treated as such with the AR(1)-model, described e.g. in Chatfield (1989).Applying it, a set of Laplace-distributed residuals r * is obtained.The error model is explained in more detail in Appendix B.

Annual CH 4 fluxes
For future climate projections, the annual total emission of the snow-free period is one of the important quantities, as cumulative emissions are what determine the radiative forcing.Therefore a Gaussian-distributed observation, G obs , is included for the total annual flux, with error variance σ 2 G , where σ G is set to 10% of the annual total flux for each year.The modeled total annual flux is denoted by G M .This term keeps the annual emission estimates reasonable in the early stages of the sampling.As the annual flux estimate errors are small when sampling parameters close to the posterior mode, this term has only a minor effect at the late stages of the MCMC.

Prior information
The parameters affecting the CH 4 production of the wetland model are not known well, but despite this, not setting any prior distributions on parameters can lead to nonphysical parameter values in the posterior distribution.
The parameter priors are set to zero outside prescribed bounds.Within these bounds, most parameters are assigned Gaussian priors, and for the others the priors are set to be flat.The prior values are based on both literature and expert knowledge and the information regarding the parameter values is summarized in Table 2.

The objective function
The objective function for the parameter optimization, J(θ θ θ), is the negative logarithm of the value of the posterior probability density function at θ θ θ.It combines our statistical knowledge of flux observations, annual flux estimates, and parameter priors presented in Sec.4.2.1 -4.2.3, and is given by: Here |r * t | are the AR(1)-transformed Laplace-distributed residuals, G M , G obs and σ 2 G are the components of the annual flux term, and the last term is the prior contribution, where θ i is the proposed parameter value, µ i is the prior mean, and σ 2 i is its variance.For technical details, see Appendix B.

Results and discussion
The experiments yielded an MCMC chain for each modeled peat depth and the final number of model simulations varied from 78000 to 391000.To look at statistics, 50% of values from the start of each MCMC chain were discarded as warm-up.The posterior covariance structures of the chains were found mostly to be similar to each other.The posterior distribution from the experiment with 70 cm of peat is shown in Fig. 2, and the correlation features for all peat depths are shown in the upper triangle of that figure.For the different processes, Fig. 11 shows an example of the posteriors and the process correlations.
For each MCMC chain, three different estimates for the parameters and fluxes were looked at: the maximum a posteriori ("MAP") estimate, posterior mean estimate ("PM"), and a "non-hierarchical" posterior mean estimate ("NHPM"), where the mean values of the parameters ζ exu and Q 10 over the different years were used.

Parameter values and modeled peat depth
The parameter values of the MAP and PM MCMC optimizations are shown in Table 3.The catotelm carbon pool turnover time, τ cato grows with the peat depth, and the reduction factor for the diffusion coefficient in air-filled peat, f D,a , also grows slightly, increasing the diffusive permeability of the dry part of the column hence increasing conductance.The root conductance gets larger with the increasing peat depth, by the influence of the parameter ρ, which grows slightly, and the decrease in root tortuousity given by parameter τ .
The MAP estimates of the different experiments disagree on the value of the water-diffusion rate coefficient f D,w , and the posteriors shown in Fig. 3 (k) are wide, especially for the experiments with less peat, meaning that especially in those cases this parameter is highly uncertain.However, with increasing peat depth, the mass of the posterior distribution mass moves closer to 1 compensating the decreased conductance caused by the longer distance to the surface.The air diffusion rate coefficient f D,a shows similar behavior, but with lower values as it is constrained by the prior.For the 40-cm peat column optimization, the MAP estimate for f D,w is far from the others.This can be explained by that in that particular case there is not much water-filled peat in the model leaving the parameter with less effect on the results.
The root distribution parameter, λ root , is optimized larger than expected, and is closer to the prior value only in the optimization with 15 and 20 layers.This is also true for the MAP estimates implying that the model optimizes best when the CH 4 produced from the photosynthesis-induced exudate production goes relatively far below the surface: with a value of 0.3, 49% of the roots are deeper than 25cm, 15% of the roots are deeper than 50cm, and just 2.5% are deeper than 75cm, see Fig. 1.In relation to these numbers, the water table depth is most of the time above the depth of -20 cm.Additionally, a larger λ root will facilitate the emission of the CH 4 produced by peat decomposition in the catotelm.The small values for the two experiments with the thickest peat column make the model behave differently from how it functions with 40-100 cm of peat.
The results in Table 3 reveal that the parameter regulating the exudate pool turnover time, τ exu , slightly decreases with the peat column depth implying a shorter period between photosynthesis and methane emission with more peat.However, ζ exu controlling the amount of methane produced from exudates gets smaller until 100 cm, and the values for 150 and 200 cm of peat are markedly larger.This implies, that CH 4 production from exudates is closely linked to the depth of the root mass and The non-hierarchical parameter V O0 controlling the amount of CH 4 oxidation taking place does not show a trend with respect to the modeled peat column depth in the PM estimate, but there is a clear trend in the MAP estimates, shallow peat depths favoring larger parameters and inducing more CH 4 oxidation.The effect of the parameter on the total CH 4 oxidation is substantial, which is evident from part (a) of Fig. 4 With all peat depths, the chains traverse in regions of both high and low V O0 as shown in Fig. 3 (f).Another parameter indirectly affecting CH 4 oxidation, the heterotrophic respiration parameter V R0 , drifts in all experiments close to its minimum value reducing the amount of heterotrophic respiration taking place.
Table 3 shows that the hierarchical parameter Q 10 , controlling the temperature dependence of the CH 4 production from peat decomposition, increases with peat column depth: the more peat there is, the stronger the peat decomposition process responds to soil temperature changes.Contrasting with this, the parameter τ cato controlling the peat decomposition rate in the catotelm increases with peat depth compensating for changes in total peat volume and keeping the production volumes reasonable.Part (d) of Fig. 4 shows how increasing the peat depth and keeping parameters constant drastically increases the total methane production.
The annual variability of Q 10 is similar across all peat depths.For instance years 2013 and 2014 (and sometimes 2011) are years of high Q 10 , whereas in 2006 the parameter gets its lowest or second-lowest value for all depths in the PM simulations.
For the other hierarchical parameter, ζ exu , these patterns do not exist.

Costfunction values and annual discrepancies
The minimum costfunction values and annual biases provide information about how well the different configurations of the model performed in the model calibration task.Table 3 lists the costfunction values for the MAP estimates, and the annual errors for the MAP, PM, and NHPM estimates are shown for each MCMC experiment in Fig. 5.Among the MAP estimates, the costfunction value is lowest with 7 layers of peat (259) and then gradually higher for 10 (262), 15 (267), 20 (269), 5 ( 274 for peat column depths of at least 50 cm.Compared to the default parameters, the estimates are much improved.This is to be expected as the results shown are not for an independent validation dataset.Rather, the motivation is to see, how the model fit looks like for optimized parameters and how the features differ from the unoptimized ones. A cross validation of the annual errors is shown for the experiments with 70 cm and 100 cm peat columns in Fig. 6 (g) and 5 (b).While the annual estimates are worse than the estimates with the optimized parameters, compared to the NHPM estimate and the default parameter values the errors are smaller.Additionally the overall biases are as good as with the optimized parameters.The cross validation is described in Sec.5.7.
Almost all years of hierarchically optimized experiments show at least a small negative annual bias when compared to the available observations.This can be due to the high day to day variability of the summertime fluxes, which dominate yearround total fluxes, and the fact that the model can not, without data about the fine structure and heterogeneity of the wetland, match the high variability fluxes.The proportional model-data residual error component α y t (Appendix B) allows the model to underestimate the high peaks more than the low flux values.The error model favors the baseline of the lower values during periods when observed variance is very high, for instance in the peak emission season of 2010.This is also true for periods of increased ebullition, and such fluxes are very difficult to fit into.These periods contribute to both the costfunction values and the underestimation of the total flux.Any temporal shifts of peaks of seasons are penalized heavily, and the optimized parameter values rather produce less peaks than right size peaks at a slightly wrong time.

Parameter values and processes in sqHIMMELI
The sqHIMMELI model produces the CH 4 from anaerobic respiration that originates from peat decay and the decay of root exudates.These production components, along with the different output pathways, CH 4 oxidation and model residuals, are plotted as functions of water table depth in Fig. 8  In the following, all ebullition refers to any ebullition in the peat column regardless to whether the bubbles reach the peat column surface.Ebullition refers to the part of "all ebullition" which reaches the surface.Most of the time the water table is under the peat surface, and at those times "ebullition" is zero, although "all ebullition" can be substantial.In that case the ebullition flux does not go directly into the atmosphere, but into the first air-filled peat layer above the WTD, and continues from there via other pathways.The reason for this separation comes from implementation details of HIMMELI.In all experiments, ebullition reaching the surface is minor fraction of the total CH 4 emission.For the PM estimate with the 100 cm column depth, the flux components and oxidation are shown as time series in Fig. 9.
Having only four layers of peat leads to peat decay being inhibited when the water table is low as the volume of the modeled catotelm decreases (Fig. 8 (f)).The effect is seen also in the plant transport component in Fig. 8 (b).Plant transport becomes proportionally more important with increasing depth of the peat column with MAP, PM, and NHPM estimates (Fig. 5 (ac)), even though the differences get quite small and the system seems to mostly stabilize already at 7 layers.For the default parameters, however, the trend is opposite (Fig. 5 (d)) as increasing the peat depth dramatically increases CH 4 production and as the default parameter set favors ebullition and diffusion over plant transport.
Comparing results from simulations with optimized parameters to results using the default parameter values (prior mean values, shown in Table 2) shows that the optimization drastically increases the role of the plant transport pathway at the expense of the diffusion pathway.Diffusion and all ebullition fluxes are closely tied to each other, as can be seen in Fig. 5 (a-d), in that in all cases their values are close to each other.This is also visible in the flux component time series in Fig. 9.

Methane production and oxidation
Figures 4 and 11 show, that there is considerable annual variation in the production of CH 4 from both of the production processes.Year 2007 has a high amount of production from peat decomposition, whereas year 2006 shows a lot less, even though the ζ exu -controlled proportion does not change much.This is not a general trend, though, and instead in years of high emissions the amount of CH 4 from both of the production sources is increased.The shape of the NPP input, shown in Fig. 7, does not change remarkably from year to year, but the emissions change considerably, as the model state and input affect the production non-linearly.For example in times of low WTD in the peak emission season, the root exudates do not contribute to CH 4 production, as much of the roots are located in the dry part of the peat column and the exudates are deposited there.
Another explanation for changes in CH 4 production comes through the production-determining parameters, whose variation is in Sect.5.7 found to be related to the springtime temperature and NPP.
The NPP-based CH 4 production controlled by the parameter ζ exu is constrained by its hyperprior as can be seen in Fig. 3 (b) and the values accepted in the chain are from the higher side of the Gaussian hyperprior, whose mean and standard deviation are both 0.2.The posterior means in table 3 are between 0.27 and 0.38, with standard deviations of 0.17-0.18.For the MAP values the deviations are even higher, leading to wider fluctuations in the characteristics of the modeled wetland.As mentioned in Table 2, the prior for ζ exu was set quite low, and actually even these values obtained here are on the low side of the spectrum reported by Walker et al. (2003), who gives a range of 0.2-0.84 in terms of our ζ exu -parameter.Our result hence agrees with that a relatively large portion of the photosynthesized sugar is respired into methane.
The methane produced by the action of ζ exu is distributed according to the root distribution, whose form is determined by λ root .The posterior means reveal, that that the contribution of the prior component of λ root to the costfunction is large.Its values might well be larger with a wider prior and more permissive prior, but in regard to how root distributions are in reality (Fig. 3.2), larger values for the parameter would make its interpretation difficult.This parameter affects both how exudates are allocated in the column and how deep the fast plant transportation reaches.Clearly there is a need to reach further down, implying that the model performs more optimally when it transports CH 4 faster to the atmosphere.to 70 cm of peat the amount of oxidation is higher (Fig. 4 (b)).In our analysis no easy explanation was found for this feature and it is suspected that multiple processes involving the parameters governing the root distributions and the availability of oxygen are behind the phenomenon.The MAP simulations (Fig. 4 (a)), however, show that the CH 4 oxidation reduces gradually with the modeled peat depth, and that the simulations with the deepest modeled peat columns have parameter estimates with the lowest V O0 -parameters.The correlation between oxidation and the value of V O0 is high, at 0.8 for the year 2008 with 70 cm peat.With more modeled peat, the methane transportation time from the lower parts of the catotelm increases, and to compensate for this V O0 in the MAP covaryingly reduces CH 4 oxidation so that the amount coming to the surface varies only a little (Table 3).This is also supported by that methane oxidation and production from root exudates covary negatively, as shown in Fig. 11.That figure also shows, that the exudate and peat decomposition based methane production terms are strongly negatively correlated, and that either of the terms can dominate production of CH 4 within the 50% confidence interval, even though in 2008 the 70 cm experiment shows overall dominance of the peat decomposition process.
The production and oxidation related parameters τ cato , and V O0 correlate (Fig. 2), and ζ exu and Q 10 are affected via their correlations with τ cato .These parameters covary producing a total emission that minimizes the likelihood but this yields a posterior, where some parameters like V O0 have wide marginal distributions (Fig. 3 (f)), as in the presence of several covarying parameters any of those covarying ones can to some degree compensate for the movement of the others.The prior bounds of V O0 were tight and for example Segers (1998) reports that potential CH 4 oxidation can vary across three orders of magnitude.
Hence, higher proportions of CH 4 oxidation could have been seen with a more permissive prior.This would have then resulted in wider posteriors also for the covarying parameters.
Curiously, parameter V R0 controlling heterotrophic respiration correlates negatively in all experiments weakly with CH 4 production via ζ exu , and positively with parameter τ cato .Removing oxygen from the column reduces CH 4 oxidation and in order to maintain the overall level of CH 4 emission, production is reduced.The wide posterior of V R0 in Fig. 3 implies that the day to day variation in the emissions and the combined effects of other parameters dominate.In the MAP estimates, however, the parameter is close to the lower bound.

Plant transport
The amount of plant transport in the calibrated models, shown in Fig. 5, is close to 90% which is slightly over the upper end of the range of 68-85% reported in Wania et al. (2010)  The parameter posteriors of λ root in the MCMC with 15 and 20 peat layers are apart from the others in Fig. 3, implying that an optimal rooting depth is an ambiguous notion.The root distribution depths also correlate differently with other parameters in different depths -there is a negative correlation between ζ exu and root depth (Fig. 2), which gets stronger with increasing peat column depth suggesting that more exudates are needed for shallow roots, which is reasonable since the exudates above the water table are respired aerobically.The other parameters affecting plant transport, ρ and τ both are included in Eq. 8 and one could expect them to be tightly coupled.In the posterior, however, they are only slightly correlated, with coefficients in Fig. 2 from 0.16 to 0.31.The strict priors may play a role, as root tortuousity cannot go below the value of 1.

Diffusion
The masses of the diffusion coefficient parameters f D,a and f D,w in the posterior distributions (Fig. 3 (j) and (k)) are above the priors.This is true especially for f D,w , which optimizes to close to the upper limit of one, specifically for the experiments with 100 -200 cm of peat.Kaiser et al. (2016) note that these parameters are not well known, and use for both of them the value of 0.8, in which light the prior for f D,a looks narrow.The PM estimates for f D,a in Table 3 are between 0.50 and 0.65 for the depths of 40 -100 cm.
The parameter f D,a is correlated negatively with the root-ending area parameter ρ.This is because the air diffusion parameter also affects the speed of CH 4 transport in plant stems via Eq. 8, and by negatively correlating the two parameters the model can compensate for one of them by moving the other.Additionally, a smaller root conductance implies that more of the CH 4 needs to come out via the diffusive flux, which is also seen in the negative correlation of λ root and f D,w , especially in experiments with peat depths of 70 and 100 cm.
Diffusion is correlated strongly with peat decay-based CH 4 production and negatively with exudate-based production (Fig. 11), and these correlations extend to the hierarchical parameters defining the CH 4 production (not shown).This is related to the strong connection between diffusion and ebullition, and that decaying peat produces CH 4 lower in the peat column than decaying exudates, production from which is more likely transported by plants (Fig. 11).
In general, the calibrations tend to end up facilitating the total CH 4 transport as the depth increases, by the action of the parameters ρ and τ affecting plant transport, and f D,a , and f D,w affecting diffusion, implying that there is a regime of optimal conductance.Pertaining to this, Fig. 6 shows how the more modeled peat there is, the less important the diffusive component becomes.Going deeper down, plant transport becomes more competitive compared to diffusion in the MAP estimates.

Ebullition
Ebullition is very strongly tied to diffusion in the flux estimates with parameters from the posterior, as is shown for the 70 cm experiment in Fig. 11.The flux component timeseries in Fig. 9 shows that ebullition to the surface is a small fraction (circa 0-3% with optimized parameters), of the total flux, and Fig. 6 shows, that the more there is peat, the less important the ebullition flux is, including the part emitted as part of the diffusive flux.Similarly, Wania et al. (2010) report almost virtually no ebullition to the surface.This result is highly dependent on the type of the wetland as for instance Kaiser et al. (2016) report high ebullition fluxes for a polygonal tundra in the Siberian permafrost region, where the ice-free soil layer reaches only about 30 cm depth during summer.Variation between different sites is very large and depends on whether the water reaches the surface at times of high CH 4 emission.
Contrasting with this, in the simulations with the default parameters and 100 cm of peat, a major part of the diffusive flux is transported by ebullition (Fig. 6) and diffusion is the dominating flux component, even though ebullition to the surface accounts for only 5% of the total flux.Since ebullition is a fast timescale process, it was not directly constrained in the optimization with parameters, as preliminary tests revealed that daily data resolution would not be sufficient for this.While finer time resolution data would have been available, using it would not have been feasible as there is not enough knowledge about the fine structure of the wetland and micrometeorological conditions affecting the footprint area of the flux tower.It is reasonable to believe that the deviations from the daily averaged fluxes at a finer time resolution would only look like noise in the residuals not improving our parameter posterior.Despite this, ebullition is controlled indirectly by letting CH 4 production and transport parameters control when the water column has enough CH 4 available for ebullition.This happens when the sum of the partial pressures of dissolved gases is larger than the sum of atmospheric and hydrostatic pressures as shown in Eq.A13.The high ebullition-related proportion of the diffusive flux strengthens the argument that the likelihood formulation results in model optimizing towards parameter values that support rapid CH 4 transport.
The results show that with deep roots and high root conductances the wet part of the peat column rarely creates the conditions for ebullition to happen.Hence with less peat the amount of "all ebullition" increases, (Fig. 5 (b), 6 (a), and 8 (d)), as the produced CH 4 needs to be stored in a smaller volume increasing its concentration.This way the modeled peat depth has a major effect on how the model transports gases.

Parameter and process identifiability
The priors of the hierarchical CH 4 production-related parameters Q 10 and ζ exu in Fig. 3 (b) and (d) are constrained by the data, as are the hierarchical parameters themselves, shown in Fig. 10.The priors of these distributions are wider than their posteriors, which is also the case for the other production-related parameters τ exu and τ cato .The different posteriors of the catotelm peat turnover time τ cato are disjoint from each other as increasing simulated peat column depth is compensated for by reducing the peat decomposition rate per volume.Both process descriptions for obtaining the anaerobic respiration are clearly needed for a good model fit, because the parameter posteriors do not have remarkable mass in the regions minimizing either of these processes (hierarchical parameters at the lower bounds or turnover rate parameters τ exu and τ cato at the upper bound).The covariances in Fig. 11 and Fig. 2 show that the two production processes covary, and hence they are partly interchangeable.
Reasonable identifiability of the Q 10 -parameters is not obvious, as for example Müller et al. (2015) optimizing a corresponding parameter end up with the parameter at the lower bound of their prescribed range.
The posterior distributions of V R0 show, that sqHIMMELI performs better when the heterotrophic respiration is close to being minimized, but still away from the lower bound.With 100 cm of peat, the parameter has a clear mode further from the lower bound, suggesting that the flux measurement data used also constrains this process, and that the prior does not indisputably rule out the best values.However, the oxidation parameter V O0 is not identifiable, and as the strong correlation with the peat decay related parameter τ cato shows, its function is partly to calibrate the total CH 4 output of the the model and to spread the posteriors of the covarying parameters.With this model and data, methane oxidation rates at the Siikaneva site cannot be estimated without further constraints for e.g. the CH 4 production.
All the parameters affecting root conductivity are constrained by the data to maximize the conductance.The root tortuousity parameter τ has narrow posteriors close to the lower bound of one, the root depth parameter λ root is above its prior, and the root-ending area parameter ρ optimizes to very high values compared to the prior distribution (Fig. 3 (g-i)).The diffusionrelated parameters f D,a and f D,w are optimized to high values and identifiable with the exception that with the shallowest peat depths the water diffusion rate coefficient has little role and a wide posterior spanning all values from 0 to 1. Transport pathway shares are stable between the MAP and PM optimizations in Fig. 6 (a), and their annual variation is small, implying that the existence of the different pathways helps to optimize the model fit.

Low WTD in 2006, 2010, and 2011
The sqHIMMELI model is not able to estimate the CH 4 flux correctly in times of low water table in the 40 and 50 cm peat depth configurations.This is not unique to this particular model -also other studies have indicated the challenges in parametrizations of emission models in response to the water table depth (e.g.Zhu et al. (2014)).Figures 8 (b) and (f) reveal that in July 2006 the plant transport component and CH 4 production from peat decay go effectively to zero during the period of the lowest WTD.This effect is not visible in the simulations with more peat, and already the PM estimate with 7 layers gives a very nice fit for July 2006 (not shown) as the deeper peat column provides for more freedom for flux adjustment.The underestimation of the emission is visible comparing the plant transport pathways and total model residuals with respect to the water table depth at different depths in Fig. 8 (b), (f), and (h).
Other extended periods of low water table occur during the years 2010 and 2011, which explains why those years tend to be accentuatedly underestimated with respect to the observed flux with shallow simulated peat columns, as is shown in Fig. 5 (f), even though curiously the 50 cm MAP estimate also performs well.The lowest water table depth of the simulation period is in 2006, when on the 26 th of August the water table drops to 38 cm below the surface.In 2011 the WTD goes below 20 cm for a total of 77 days in a row, and in 2010 it recedes to -17 cm for a period of 71 days, the average of the period being -23 cm.Years 2010 and 2011 have the strongest tendency to underestimate the total annual CH 4 emission compared with the observations.A sufficiently deep peat column to accommodate for CH 4 emissions during the low WTD periods is needed for making accurate predictions .

Optimal modeled peat depth for sqHIMMELI
Even though most of the parameters and processes are identifiable, all of the parameter posteriors vary with peat depth, the most striking example of which is τ cato (Fig. 3 (c)).For this reason, the validity and meaning of the parameter values must be understood in each particular model setting.
The objective function incorporating prior knowledge can be used to evaluate what peat column depth best represents the data and still retains the physical interpretation of the parameters, information about which is in the prior parameter distributions.
In the MAP estimations, the costfunction values (Table 3) and the annual flux estimate errors (Fig. 6) are smaller starting with the depth of 70 cm and especially the 40 cm optimizations are systematically worse in this respect than the others, due to worse handling of periods of low WTD (Fig. 8).These problems do not exist with the 70 or 100 cm simulations, and are less pronounced already with 50 cm of peat.
With 150 or 200 cm of peat the correlations of the parameters shown in Fig. 2 show markedly different patterns from the correlations with shallower modeled columns.For these thick peat columns the costfunction values are higher, the correlations are not easy to explain, annual negative biases are not better, and model integration is more costly in terms of CPU time.For these reasons there is no reason to believe that modeling deeper peat columns than 100 cm in sqHIMMELI would be superior.
Rather, the optimal thickness lies between 50 and 100 cm.

Predicting emissions with sqHIMMELI
Modeled CH 4 flux estimates may have large errors as was shown in Fig. 6 (b) with the default parameter set.The negative biases of less than ten percent in the calibration phase that were found with the PM estimates are reasonable since the quality of the modeled input data from e.g. a land surface scheme will also contribute to the uncertainty in the model predictions.
Compared to the estimate with the optimized annual variations of the Q 10 parameters, the posterior mean estimate without the hierarchical parameters (NHPM) does not produce very good good flux estimates over the assessment period (Fig. 5 (g)).
With all peat depths, the total CH 4 emission of the first years is overestimated by up to 30 percents and for the last years there is a similar negative bias.The hierarchical posterior mean (PM) on the other hand does produce very steady estimates of the CH 4 flux, compared with observations, and for these estimates the model dynamics are similar between the estimates for the peat depths of up to 100 cm.
In order to be able to utilize the PM estimates for the future prediction of CH 4 emissions, the values of the hierarchical parameters need to be estimated for the simulation years.A simple regression analysis of the hierarchical variables with respect to relevant input data was performed in order to find out if such estimation is possible.As the explaining variables, means, minimums, and maximums of NPP, water table depth, and soil temperature at different depths and over different periods of time were looked at.These time periods were June, July, August, and various different amounts of days from the start of the year.
The analysis revealed that the mean soil temperature of the first 10 weeks of the year at the depth of 30-40 cm, denoted here by T 70 30−40 , is the best single-variable predictor of the Q 10 -value for that year, and for ζ exu , it is the sum of NPP from the first For the MCMC experiment with 40 cm of peat the p value of the regression is better when looking at the 20-30 cm average soil temperature (p = 0.075, r 2 = 0.38), than with the 30-40 cm temperature (p = 0.13, r 2 = 0.30).It is understandable that estimating Q 10 gets more difficult with a shallow peat column depth, because the parameter has less effect as in the summer a majority of the whole peat column is dry.For the simulations with 150 or 200 cm of peat the regression for Q 10 does not give meaningful results and the p values are large, implying that a deep active peat column might pose an additional degree of difficulty for performing the CH 4 emissions.
The best experiments in terms of the predictability of the hierarchical variables are those with 70 and 100 cm of peat, which also were the best performing peat column depths among the model calibration results (Table 3).For those depths, the hierarchical parameters Q 10 and ζ exu for each year can be estimated with ζ 7l exu = −56000N P P 130 + 0.468 (13) where the upper indexes 7l and 10l refer to the number of 10 cm peat layers, temperatures are in • C, and the units of NPP are mol m −2 s −1 .
The lower p and higher r 2 values for the 70 and 100 cm models suggest that also in terms of predictive skill these configurations are superior as the hierarchically varying variables can be more robustly estimated.A leave one out-cross validation (LOO-CV, see e.g.Gelman et al. ( 2013)) of the predicted fluxes was therefore performed on the 70 cm and 100 cm models by optimizing them with respect to the costfunction in Eq. 11 leaving one year at a time out, calculating the estimates for the hierarchical parameters based on the results obtained, and predicting the CH 4 emissions for the year that was left out.The algorithm (BOBYQA) and the number of iterations completed (350) were the same as before the MCMC and for the hierarchical parameters the priors were defined by the values defining the hyperprior.The results of the cross validation are shown in Fig. 5 (g) and 6 (b).Compared to the NHPM estimate and the default parameter values the annual errors were reduced and the mean annual errors were -4.95% and -3.01%with the standard deviations of 13.6% and 13.4 %, for the 70 cm and 100 cm peat column depths respectively.In comparison to the NHPM estimate the annual errors were reduced for six out of the nine years.
These results are promising, and as the analysis performed was extremely simple, there is room for further development.
where T X (t, z) describes transport of gas X containing the diffusion, ebullition, and plant transport components, and R stands for production or consumption.The different terms in the equations are described below.

A2 Anaerobic respiration producing CH 4
The equations presented in this section are specific to the version of HIMMELI used in this study.The version in Raivonen et al. (2017) takes the rate of anaerobic decomposition of carbon as input and does not treat the different sources of that carbon separately.
The carbon for methane production in this model version comes from two sources: root exudates and anaerobic peat composition.The methane production from anaerobic respiration of that carbon is given by the terms R exu CH4 and R peat CH4 described by: where in Eq.A5 ν is the decay rate of root exudates from Eq. 3, η is an oxygen inhibition parameter, C O2 (z) is the oxygen concentration at depth z, and π(z) is the normalized proportion of the total anaerobic root mass, also at depth z, given in an unnormalized form in Eq. 2. The parameter f exu CH4 is a constant determining what fraction of root exudates in anaerobic conditions will turn into CH 4 .Equation A5 is only used below the water table.In Eq.A6, g Q10 m is the proportion of the anaerobic peat decomposition process producing CH 4 , ρ cato is the peat density in the catotelm, f Ccato is the fraction of carbon in catotelm peat, and M C is mass of carbon.The parameter k cato = Q (T −273.15) 10 10 /τ cato is described in Eq. 5, and is zero above water table.
The equations for CO 2 are similar: and the meanings of the symbols are analogous to the ones in equations for CH 4 .

A3 Peat respiration and methane oxidation
Peat respiration (aerobic respiration) is described with an equation of the Michaelis-Menten form where α is a dimensionless Henry solubility constant for oxygen above the water table, and one below it, see Tang et al. (2010).
The factor C x O2 refers to C w O2 below the water table, and to C a O2 above it.Here w and a refer to whether the concentration is in the gaseous or in the liquid phase.Parameter K R is the Michaelis-Menten constant of the process, and V R (z) is given by Eq.

6.
Methane oxidation is controlled by dual-substrate Michaelis-Menten kinetics, and here α factors similar to the one in Eq.A9 have been absorbed into the concentration terms -otherwise the terms are analogous to those in Eq.A9, except for that the term V O (z) is described by Eq. 7.

A4 CH 4 transport
The transport term T X (t, z) in Eq.A1 consist of the following terms: The first of these is the diffusion, where the diffusion coefficients D are given by Eq. 9 and 10, and "medium" refers to either air or water.Due to coding mistake, the f D,a and f D,w coefficients in the aforementioned equations were set to 0.1 for gases other than CH 4 in this work.
The second equation is for plant transport, with ρ and τ described in context of Eq. 8, π(z) is the normalized root distribution mentioned above, and C atm X refers to the atmospheric partial pressure of gas X.LAI stands for the leaf area index, given as input, and SLA is the specific leaf area.The note above regarding the f D,a values is also valid for plant transport, as it is a factor determining D X air .The third equation is the ebullition component of the gas transport, where pp i refers to the partial pressure of different gases indexed with i, R is the universal gas constant, k is an ebullition rate constant, and σ is the peat porosity.The parameters P atm and P hyd (z) refer to the atmospheric pressure and hydrostatic pressure at depth z, respectively.
Table 5 shows the parameters that are used in the equations above but not optimized in this work, along with their values.The unknown hyperparameters µ i and σ 2 i have probabilistic models where µ 0 and τ 2 0 define the mean and variance of the hyperprior of µ i , n 0 ∈ N defines the number of degrees of freedom of the Inv-χ 2 distribution, and σ 2 0 is the expected value of the scaled Inv-χ 2 distribution.
In Gibbs sampling the full conditional posterior distributions of the hyperparameters and the parameters θ i are sampled in turns.Due to the conjugacy of the normal distribution and the scaled Inv-χ 2 distribution, closed form expressions exists for sampling from p(µ i |σ 2 , µ 0 , σ 2 0 , θ i ) and p(σ 2 i |σ 2 0 , n 0 , θ i ), where µ is the current mean of the parameters θ i and σ 2 is their variance.The Gibbs sampling therefore consists of three steps: , and (D5) 3. draw the parameters θ i (and the non-hierarchical parameters) with MCMC, since closed-form expression for p(θ θ θ|φ,y y y), where φ denotes all the different hyperparameters, is not available.
In this work, the value of the parameter τ 2 0 was set to the value of σ 2 0 , n i is the number of years, and the value of n 0 was set to 9. The means and variances obtained this way describe the interannual variability of the parameters, and not including them as parameters in the MCMC sampling reduces the dimension of space that the MCMC sampler needs to explore, speeding up convergence of the posterior distribution.by averaging the values reported by (Wilson et al., 2007) for the vascular species abundant at Siikaneva.For the growing season peak LAI we used the maximum LAI observed at the eddy covariance footprint area, viz.approximately 0.4 m 2 m −2 (Riutta et al., 2007b).We also included a constant wintertime LAI since a significant green sedge biomass may overwinter, approximately 15% of the maximum (Saarinen, 1998;Bernard and Hankinson, 1979).The overwintering LAI at Siikaneva would thus be 0.05 m 2 m − 2 .The same LAI was used for all the years and this LAI also was given as the input for the CH 4 transport model.
The daily averages of P n were calculated by subtracting R a from P g .The models were run with measured meteorological data.We determined the photosynthetically active seasons based on snowmelt dates in spring or arrival of snowcover in autumn from the reflected PAR data, or based on air temperature (permenently greater than 5 • C assumed to be the growing season).
After the calculation, we compared the resulting P n of vascular vegetation of year 2005 to eddy covariance CO 2 fluxes from Siikaneva.We used the GPP derived from the measured NEE by (Aurela et al., 2007).This was the only available year of processed CO 2 flux data.The GPP was on average 4.5-fold compared with our P n , with a R 2 of 0.9.GPP also includes the photosynthesis of Sphagnum mosses as well as CO 2 released in autotrophic respiration.Sphagnum accounted for 20-40% of the GPP in the study by (Riutta et al., 2007a) and autotrophic respiration has been observed to be roughly 50% of GPP (Gifford, 1994).Consequently, the NPP of vascular vegetation can be estimated by multiplying the GPP with 0.7 × 0.5.This estimate was still 1.56-fold compared with the P n for the year 2005.Since the P n also was lower than generally reported for peatlands, we chose to trust the eddy covariance measurement and scaled the P n of all the years upwards by multiplying with 1.56.For further details, please consult Raivonen et al. (2017).
Appendix F: Supplementary information F1 Details regarding the AM algorithm usage In order to infer about the posterior distribution, the MCMC chain needs to be long enough, and converged to produce the right statistics.The MCMC chains driven by the AM algorithm mixed well, example of which as can be seen in Fig. 12 showing the chain from the experiment with 100 cm of peat.The proposal distribution of the AM algorithm was adapted when the iteration number was a square of an integer, and for the adaptation 20% from the start of the chain was discarded.In the early stages of each experiment, the initial approximation for the proposal covariance, calculated from the Jacobians of the model, was allowed to dominate until after accepting enough proposed points there was sufficient data to start the proposal covariance adaptation procedure.

F2 Computational requirements
Even though the model runs fast, in around five to thirty seconds for the ten-year period on a multicore laptop, due to the large number of simulations, the MCMC experiments needed to be performed on a CRAY XC-20 supercomputer using a single node for a single MCMC chain and running all the experiment in the RAM of the computer minimizing hard drive utilization.Q10 1.7 16.0 -5.9 2.0 (Juottonen, 2008;Gedney et al., 2004;Bergman et al., 2000) ζexu 0.01 0.99 -0.2 0.2 Lower bound of (Walker et al., 2003) Table 2.
), Geosci.Model Dev.Discuss., doi:10.5194/gmd-2017-66,2017 Manuscript under review for journal Geosci.Model Dev. Discussion started: 21 April 2017 c Author(s) 2017.CC-BY 3.0 License.and 4 layers (295), the last of these being significantly worse than the first ones.Figure 5 suggests that for the MAP and PM estimates, the annual total CH 4 flux estimates are less steady with only 40 or 50 cm of peat.In figure Fig. 5 (g) the NHPM estimate shows a very large variance of the annual errors, with early years having a sizable positive bias, and later years having a similar negative bias.Incidentally the averages over the whole period for NHPM are small as Fig. 6 (b) indicates.The model estimates of the annual fluxes are good and the variance of its errors is small for both MAP and PM experiments, especially for the optima of various MCMC experiments, and for the unoptimized default parameter values.The process correlations and covariances are shown for the year 2008 from the experiment with 70 cm of peat in Fig. 11.
Geosci.Model Dev.Discuss., doi:10.5194/gmd-2017-66,2017 Manuscript under review for journal Geosci.Model Dev. Discussion started: 21 April 2017 c Author(s) 2017.CC-BY 3.0 License.The methane production from decomposition of peat in anaerobic conditions is aided by the rather strongly correlated parameters Q 10 and the catotelm carbon decay half-life τ cato .Unlike with ζ exu , the prior means of Q 10 are mostly inside the 1-σ bounds of the hyperprior, and the temperature dependence of the anaerobic respiration from peat decomposition is close to what was a priori expected.All years of the PM simulations have very little oxidation taking place with over 70 cm of peat and in the experiments with up in a study simulating CH 4 emissions for seven boreal peatlands.This is opposite to what was obtained with sqHIMMELI with the default (prior mean) parameter values and 100 cm modeled peat, where the simulation routed 71% of the flux via diffusion.Geosci.Model Dev.Discuss., doi:10.5194/gmd-2017-66,2017 Manuscript under review for journal Geosci.Model Dev. Discussion started: 21 April 2017 c Author(s) 2017.CC-BY 3.0 License.The high optimized share of plant transport is due to the deep roots and high root conductance from the high values of the root depth controlling parameter λ root and the root ending cross section parameter ρ, and the low values of root tortuousity parameter τ .These parameters are close to the limits of what the priors allow, and are the reason for that plant transport dominates the gas transport.Wania et al. (2010) used the parametrization from Eq. 1 with λ root = 0.2517, and the root distributions from the PM estimates are shown alongside that distribution in Fig. 1.Compared with measurements from Saarinen (1996), the amount of roots at 20-60 cm is exaggerated by all of the optimized parameter values.The model provides a better fit to the data when the root conductance is maximized.
Geosci.Model Dev.Discuss., doi:10.5194/gmd-2017-66,2017 Manuscript under review for journal Geosci.Model Dev. Discussion started: 21 April 2017 c Author(s) 2017.CC-BY 3.0 License.130 days of the year, denoted by N P P 130 .This is hardly surprising, since the peat decomposition process regulated by the parameter Q 10 is driven by soil temperature, and the anaerobic respiration from exudates controlled by the parameter ζ exu is driven by the NPP input.The p values summarizing the reliabilities of the regressions and the r 2 values, which are the coefficients of determination of the fit, are presented in table 4. The r 2 values explain what fraction of the variance of the dependent (predicted) variable is explained by the independent (explaining) variables.

Figure 1 .
Figure1.The different root distribution descriptions.The original description is shown as the decaying exponential, and the graph with discrete steps shows measurement data fromSaarinen (1996).The new root distribution curve with optimized parameters are shown along with the curves resulting from the MCMC optimization.The original distribution gives more root mass to depths of 50-80cm, than the MCMC-optimized curves of the new root distribution.All curves are normalized to the same total root mass.

Figure 2 .
Figure 2. Posterior distributions of the parameters from the MCMC.The two-dimensional marginal distributions of the posterior distribution from the experiment with 70 cm peat depth are shown in the triangle on the lower left (labels on the left and at the bottom), and the correlations between parameters are shown in the upper triangle on the right (labels on the left and on the top).The images in the lower left triangle show the 90% (black) 50% (red), and 10% (blue) contours, and roughly a one-thousandth of the points sampled.The first 50% of each chain was discarded as a warm-up period.On the upper right, each plot shows correlation coefficients between parameters for all experiments, color-coded to show negative correlations in blue and positive in red.Left-to-right and top to bottom in each color-coded square, the depths are 40, 50, 70, 100, 150, and 200 cm of peat.The values from the posterior of the 7-layer experiment, referring to 70 cm peat and shown also in the lower left part of the figure, are marked in boldface.

Table 3 .
Parameter values obtained in the optimization of the models with the different peat depths.Values shown are for the maximum a posteriori (MAP) / posterior mean (PM) estimates.The horizontal line in the middle separates the hierarchically optimized parameters (including their priors) from the others.

Figure 3 .Figure 5 .
Figure 3. Posterior marginal and prior distributions for all MCMC experiments and all parameters: (a-d) are the production-related, (e-f) the respiration and oxidation related, and (g-k) the gas transport related parameters.The curves shown for the MCMC experiments are smoothed slightly using Gaussian kernel estimates for readability.To make these figures, 50% from the start of each chain was discarded as warm-up.The dotted vertical lines show the default parameter values.For the parameters ζexu (b) and Q10 (d), the prior distribution drawn is the hyperprior.