Articles | Volume 17, issue 6
Development and technical paper
20 Mar 2024
Development and technical paper |  | 20 Mar 2024

Optimising CH4 simulations from the LPJ-GUESS model v4.1 using an adaptive Markov chain Monte Carlo algorithm

Jalisha T. Kallingal, Johan Lindström, Paul A. Miller, Janne Rinne, Maarit Raivonen, and Marko Scholze

The processes responsible for methane (CH4) emissions from boreal wetlands are complex; hence, their model representation is complicated by a large number of parameters and parameter uncertainties. The arctic-enabled dynamic global vegetation model LPJ-GUESS (Lund–Potsdam–Jena General Ecosystem Simulator) is one such model that allows quantification and understanding of the natural wetland CH4 fluxes at various scales, ranging from local to regional and global, but with several uncertainties. The model contains detailed descriptions of the CH4 production, oxidation, and transport controlled by several process parameters.

Complexities in the underlying environmental processes, warming-driven alternative paths of meteorological phenomena, and changes in hydrological and vegetation conditions highlight the need for a calibrated and optimised version of LPJ-GUESS. In this study, we formulated the parameter calibration as a Bayesian problem, using knowledge of reasonable parameters values as priors. We then used an adaptive Metropolis–Hastings (MH)-based Markov chain Monte Carlo (MCMC) algorithm to improve predictions of CH4 emission by LPJ-GUESS and to quantify uncertainties. Application of this method on uncertain parameters allows for a greater search of their posterior distribution, leading to a more complete characterisation of the posterior distribution with a reduced risk of the sample impoverishment that can occur when using other optimisation methods. For assimilation, the analysis used flux measurement data gathered during the period from 2005 to 2014 from the Siikaneva wetlands in Southern Finland with an estimation of measurement uncertainties. The data are used to constrain the processes behind the CH4 dynamics, and the posterior covariance structures are used to explain how the parameters and the processes are related. To further support the conclusions, the CH4 flux and the other component fluxes associated with the flux are examined.

The results demonstrate the robustness of MCMC methods to quantitatively assess the interrelationship between objective function choices, parameter identifiability, and data support. The experiment using real observations from Siikaneva resulted in a reduction in the root-mean-square error (RMSE), from 0.044 to 0.023 gC m−2 d−1, and a 93.89 % reduction in the cost function value. As a part of this work, knowledge about how CH4 data can constrain the parameters and processes is derived. Although the optimisation is performed based on a single site's flux data from Siikaneva, the algorithm is useful for larger-scale multi-site studies for a more robust calibration of LPJ-GUESS and similar models, and the results can highlight where model improvements are needed.

1 Introduction

Methane (CH4) is the second most important long-lived greenhouse gas after carbon dioxide (CO2) (Ciais et al.2013; Kirschke et al.2013). It has been reported that the global atmospheric CH4 concentration has been growing since the pre-industrial period. In 2021, it reached a value of 1908 ppb (parts per billion), nearly 2.62 times greater than its estimated value in 1750 (Dlugokencky2021). This increase in the atmospheric concentration of CH4 is responsible for around 16.5 % of the total effective radiative forcing (in W m−2) of the well-mixed greenhouse gases (IPCC AR6: Forster et al.2021).

Among the biogenic sources, wetlands contribute around 19 %–33 % of current global terrestrial CH4 emissions and are the largest and the most uncertain (Kirschke et al.2013; Saunois et al.2020; Bousquet et al.2006). Wetlands occupy around 3.8 % of the Earth's land surface and are mainly located in high-latitude regions. There is approximately 455 Pg of carbon stored in boreal and subarctic wetland peat/Histosol soils. Under long-term anaerobic soil situations, this carbon will be metabolised by the anaerobic microorganisms (called methanogens) and will eventually be emitted as CH4 back to the atmosphere (Aurela et al.2009).

In the future, climate change may cause a positive feedback on CH4 emissions from wetlands due to a warmer and wetter climate (Johansson et al.2006; Bridgham et al.2008). According to Zhang et al. (2017), at the end of the 21st century, 38 %–56 % of the CH4 production from the wetlands will be climate change induced. Increased uncertainty in CH4 emission from boreal wetlands (Christensen et al.2007) is also predicted, partly due to expected spatio-temporal changes in wetland extent (Saunois et al.2016). Considering the fragility of boreal wetlands in a changing environment (Jacob et al.2007), one way to quantify their carbon budget is to model their carbon dynamics, including CH4 emission. Realistic and optimised process-based vegetation models can be used to reach a more precise estimation of emission variability and trends. However, representation of the complex biogeochemical processes, including soil carbon turnover, vegetation dynamics, hydrology, soil thermal dynamics, and defining wetland boundaries are complex; therefore, estimating the contribution from multiple pathways for CH4 production, consumption, and release complicates wetlands CH4 modelling (Ahti et al.1968; Wania et al.2010, 2013; Susiluoto et al.2018).

The Lund–Potsdam–Jena General Ecosystem Simulator (LPJ-GUESS) (Smith et al.2014) is one of a few available process-based dynamic global vegetation models (DGVMs) that simulate local to global vegetation dynamics and soil biogeochemistry (Smith2001; Sitch et al.2003). Using information about the climate and the concentration of CO2 in the atmosphere, it predicts the structural, compositional, and functional properties of the native ecosystems of major climate zones of the Earth. Considering the complexity of LPJ-GUESS, with its large number of uncertain process parameters, the model requires a mathematically robust framework for parameter optimisation (Wramneby et al.2008). Data assimilation using Bayesian statistics can be seen as a way of combining observations with prior information (i.e. model process formulation and prior model parameter values) to derive posterior parameter and emission estimates (Susiluoto et al.2018; Ghil and Malanotte-Rizzoli1991; Dee2005; Carrassi et al.2018). The Markov chain Monte Carlo (MCMC) method (Metropolis et al.1953b) is a powerful and convenient Bayesian framework (Tarantola1987) for data assimilation, as it can combine prior information with observations to sample from the posterior distributions in complex models. This study has developed an adaptive MCMC Metropolis–Hasting (AMCMC-MH) framework (Hastings1970b; Tarantola1987) with “Rao-Blackwellised” adaptation (transformation of estimators using the Rao–Blackwell theorem) of the multivariate Gaussian random walk proposals (Andrieu and Thoms2008). The algorithm minimises the model–data misfit (i.e. a cost function) by sampling from the probability density function (PDF) of the posterior parameters. The adaptation allows the algorithm to learn the shape of the posterior, thereby improving sampling efficiency. The main objective of this paper is to evaluate the capabilities and limitations of the AMCMC-MH framework to optimise CH4 wetland emissions simulated by the LPJ-GUESS model by analysing the posterior parameter distributions, the parameter correlations, and the processes they control.

2 Data and methodology

2.1 Siikaneva wetland and measurements

The Siikaneva wetland is located at 61°49 N, 24°12 E, and 160 m a.s.l., and it is the second largest undrained wetland complex in Southern Finland (Ahti et al.1968; Rinne et al.2007). This boreal wetland complex has an area of 12 km2, including minerotrophic and ombrotrophic sites with over 6 m of peat deposition under the surface (Mathijssen et al.2016; Aurela et al.2007; Rinne et al.2007). The estimated average annual total precipitation is about 707 mm. The average temperatures for January and July are approximately −7.2 and 17.1 °C, respectively. The estimated mean annual temperature is around 4.2 °C (Korrensalo et al.2018). Total annual CH4 emissions from the Siikaneva wetland vary between 6.0 and 14 gC m−2, while net CO2 fluxes vary between −96 and 27 gC m−2 (Rinne et al.2018).

The CH4 flux data were obtained from (Alekseychik et al.2019a). Daily measurements of incoming short-wave radiation (swr), precipitation, and air temperature collected at the wetland were used as input to the model. As the meteorological data measured directly at the Siikaneva wetland have several significant gaps, which made them unsuitable as input for the model, we used precipitation and temperature data collected at a nearby station called Juupajoki – Hyytiälä (, last access: 5 February 2024), located approximately 5.5 km from Siikaneva. The swr data collected at the Hyytiälä weather station, situated around 6 km from Siikaneva, were obtained from (Aalto et al.2022). Given the short distances between these sites and Siikaneva, we assumed that the meteorological variables were representative of Siikaneva. To verify this assumption, we analysed the available data from Siikaneva and the datasets collected from the Juupajoki and Hyytiälä sites. The air temperature and precipitation from the Juupajoki and Siikaneva sites showed Pearson correlation values of 0.998 and 0.706, respectively. The swr data collected at Hyytiälä and Siikaneva showed a correlation of 0.98. However, there were some minor gaps in the swr data collected at Hyytiälä that were gap-filled using the available data collected at Siikaneva, which were obtained from (Alekseychik et al.2019b) for the corresponding periods.

Additional inputs to the model are the atmospheric CO2 concentration, as described by McGuire et al. (2001) and updated until recent years using data from (last access: 5 February 2024). Daily water table depth (wtd) and soil temperature at a 5 cm depth collected from the Siikaneva site are obtained from (Alekseychik et al.2019b) and are used to evaluate the modelled values.

2.2CH4 model description in LPJ-GUESS

Compared with version 4 of LPJ-GUESS described by Smith et al. (2014), version 4.1, which we used for this study, has more detailed representations of plant functional type (PFT) characteristics and processes in wetlands. This includes improved descriptions of peatland-specific PFTs, peatland hydrology, soil temperature estimation, and CH4 emissions. Brief descriptions of the important wetland processes in LPJ-GUESS version 4.1 are given below and in Sect. S1 in the Supplement. For a more detailed description, the reader is referred to Gustafson (2022).

Figure 1Schematic representation of the CH4 model in LPJ-GUESS with the CENTURY soil organic model implemented. Carbon for methanogens is allocated to soil layers based on the distribution of roots in each layer. The root density decreases from the top of the peat to the bottom. The assigned carbon in each layer is divided into CH4 and CO2. Oxygen (O2) either directly diffuses or is transported through plants. The availability of O2 determines the amount of CH4 in the soil, as it oxidises a fraction of CH4. Similarly, CH4 can also either directly diffuse or be transported to the atmosphere in bubbles or it can be transported by vascular plants. The equilibrium between gaseous bubbles of CH4 and dissolved CH4 in water is controlled by the maximum solubility of CH4. Any CH4 that exists in gaseous form will escape to the atmosphere via ebullition.


The active wetland peat in the LPJ-GUESS is represented by a 1.5 m deep column that is further divided into 15 layers of 0.1 m thickness each (see Fig. 1). The uppermost three layers comprise the acrotelm, within which the water table can vary. The underlying 12 layers of catotelm are permanently saturated with water (Wania et al.2009a). The decomposed organic carbon on each day (explained in Sect. S1.4) is distributed vertically in different peat soil layers in proportion to an assumed static root distribution (see Eq. 1). This carbon pool is considered to be the “potential carbon pool” for methanogenic archaea and is the basic concept behind the CH4 model in LPJ-GUESS. The total available carbon is decomposed into two components – CO2 and CH4 – depending on the availability of O2 in the soil. The dissolved CH4 concentration and the gaseous CH4 fraction are calculated based on the estimated CH4 content in each layer. A portion of the estimated CH4 is oxidised by soil O2, whereas the remainder is transported to the atmosphere by diffusion, ebullition, or plant-mediated transport. Apart from being the key factor in estimating the potential carbon pool, root biomass in each soil layer also plays a role in the transport of O2 and CH4 into and out of each layer, as this transport is mediated by the plants. From different studies of various wetland PFTs, Wania et al. (2010) observed an exponential decrease in the root biomass with depth that was proportional to the degree of anoxia, which is expressed by the following equation, also used in LPJ-GUESS:

(1) f root = C root e z / λ root ,

where froot is the fraction of root biomass at a certain depth z, λroot=0.2517 m is the decay length, and Croot=0.025 is a normalisation constant. This distribution ensures that approximately 60 % of the roots are distributed within the acrotelm, and the root fraction in the lowest soil layer is adjusted to achieve a total root distribution of 1 across all 15 soil layers.

2.2.1CH4 production

Due to its wide range, the CH4/CO2 ratio from decomposition is challenging to predict. For example, Segers (1998) observed high variation in the CH4-to-CO2-production molar ratio of between 0.001 and 1.7 under anaerobic conditions. Hence, it is defined in the model as an adjustable parameter weighted by the degree of anoxia (α): α=1-(Fair+fair), where Fair is the fraction of air in the soil layers and fair is the fraction of air in peat (Wania et al.2009a) (see Sect. S1 for details).

The production of CH4 on each day in each layer is determined as follows:

(2) CH 4 prod = α ( z ) × f root ( z ) × CH 4 / CO 2 × R h ,

where α(z) is the degree of anoxia at depth z, froot(z) is the fraction of roots in the peat at depth z, CH4/CO2=0.085 (prior value in the model) is the tuning parameter for the CH4-to-CO2-production ratio, and Rh is the daily heterotrophic respiration. Note that the model is set to CH4prod=0 when Fwater<0.1, assuring no CH4 production in frozen and/or dry soils (i.e. the model assumes there is no water when the water is frozen and, hence, Fwater is 0).

2.2.2CH4 oxidation

As mentioned above, the CH4 fraction that is oxidised (represented by the parameter foxid=0.5 as the prior value in the model) depends on the availability of O2 in the soil. A part of the O2 transported to the soil will be consumed by the plant roots and non-methanotrophic microorganisms. The remaining part is then used to oxidise CH4. The oxidised CH4 is added to the CO2 pool, and the remainder stays in the CH4 pool.

2.2.3 Total CH4 flux

Diffusion, ebullition, and plant-mediated transport are the three pathways via which CH4 is transported to the atmosphere. The total CH4 flux from high-latitude wetland patches in the model is represented as follows:

(3) F CH 4 = CH 4 diff + CH 4 plant + CH 4 ebul ,

where CH4diff is the CH4 flux component from diffusion, CH4plant is the CH4 flux component from plant-mediated transport, and CH4ebul is the CH4 flux component from ebullition. As the daily CH4 production in each layer is dependent on Rh (Eq. 2), FCH4 is subtracted from Rh before saving the daily heterotrophic respiration. Any CO2 generated, whether from heterotrophic respiration or CH4 oxidation, is released into the atmosphere.


The fractions of CH4, CO2, and O2 that are transported to the atmosphere and from the atmosphere through diffusion are calculated by solving the gas diffusion equation within the peat layers using a Crank–Nicolson numerical scheme with a time step of 15 min. The molecular diffusivities of these gases in soil depend on temperature, soil porosity, and the water and air content of the soil. Diffusivity in water is derived by fitting a quadratic curve to observed diffusivities at different temperatures, as described in Broecker and Peng (1974); diffusivity in the air and its temperature dependency are derived from the values taken from Lerman (1979); and diffusivity in soil and its temperature dependency are estimated from the Millington–Quirk model, as described in Millington and Quirk (1961). A detailed description can be found in Wania et al. (2010).

At the water–air surface, the gas diffusivities change by a minimum of 4 orders of magnitude; hence, at the water–air boundary, the flux is calculated by the following equation:

(4) J = - ψ C surf - C eq ,

where Csurf is the surface water gas concentration; Ceq is the concentration of gas in equilibrium with the atmospheric partial pressure, estimated using Henry’s law; and ψ is the gas exchange coefficient, also called the piston velocity, which is usually difficult to estimate for different gases. In this case, the piston velocities of CH4, CO2, and O2 are calculated by relating them to the known piston velocity of SF6 via the following equation:

(5) ψ * = ψ 600 S c * 600 n ,

where ψ600=2.07+0.215×U101.7 is the piston velocity of SF6 normalised to a Schmidt number of 600 (subjected to the wind speed U10 at 10 m from the ground, which is considered to be 0 in the model), Sc* represents the Schmidt number of the gas under consideration, and n=-1/2 (see Wania et al.2010, for details).

As mentioned above the diffusion through the soil is affected by soil porosity, i.e. by the value of Fair(t,z). When Fair≤0.05 in a given soil layer, the diffusivity values of water are used; otherwise (Fair>0.05), the diffusivity values of air, which are 4 orders of magnitude larger than those of water, are used. At each daily time step, before diffusion is calculated, the gas flux J at the boundary is used to update the dissolved gas content. The surface concentration Csurf of CH4 will mostly be greater than Ceq; hence, J will be negative, denoting a flux to the atmosphere, although it is possible for CH4 to diffuse into the soil in small amounts if the concentrations at the surface are suitable. The resulting daily flux of CH4 is determined as the total CH4diff.


Ebullition depends on the solubility of CH4 at a given temperature and pressure and occurs when the water table reaches the surface during periods of high CH4 emission. Following Wania et al. (2010), in LPJ-GUESS, the best-fitted curve is represented as follows:

(6) S B = 0.05708 - 0.001545 T + 0.00002069 T 2 ,

where SB is the Bunsen solubility coefficient (i.e. the volume of gas dissolved per volume of liquid at atmospheric pressure and a given temperature).

The CH4 in each layer is converted to a maximum allowable dissolved mass, and this limit is used to separate the CH4 in the form of dissolved and gaseous components. If there is any CH4 that exceeds the maximum solubility of a layer, it will be released into the atmosphere. The CH4ebul is calculated by adding the ebullition fluxes from all layers.

Plant-mediated transport

Plant-mediated transport of CH4 occurs via the aerenchyma (the gas-filled tissues) of vascular plants, either through concentration gradient or active pumping from the soil to the atmosphere. Only the passive mechanism (through a concentration gradient) is considered in the model, as it is the most dominant one (Cronk and Fennessy2016). Abundance, biomass, phenology, and the rooting depth of aerenchymatous plants are considered to calculate this. Only the flood-tolerant C3 graminoids are considered for plant-mediated gas transport in the model (see Table S2 in the Supplement); hence, plant-mediated transport of O2 and CH4 can only occur when C3 graminoids are present in a simulated patch.

The transport depends on the cross-sectional area of plant tillers (all of the secondary shoots produced by grasses like Poaceae or Gramineae) in each soil layer, assuming that a significantly high percentage of CH4 is oxidised in the highly oxic zone near the roots, where methanotrophs flourish, before they enter into the plant tissue.

The mass of the tiller is calculated as follows:

(7) m tiller = b graminoid × P ( leaf ) ,

where bgraminoid is the leaf biomass of graminoids, and P represents the daily phenology, which is the fraction of potential leaf cover that has reached full development. To calculate the number of tillers (ntiller), the total weight of tillers (mtiller) is divided by the average weight of an individual tiller (wtiller). The cross-sectional area of tillers (Atiller) can then be obtained as follows:

(8) A tiller = n tiller × ϕ tiller × π r tiller 2 ,

where rtiller is the tiller radius and ϕtiller is the tiller porosity. Based on the optimisation of McGuire et al. (2012), Tang et al. (2015), and Zhang et al. (2013), the value of rtiller is estimated as 0.0035 m, and based on Wania et al. (2010), the values of ϕtiller and Wtiller are estimated as 70 % and 0.22 gC per tiller, respectively. Each soil layer is allocated a fraction of the total cross-sectional area of tillers based on the root fraction estimated in that layer. The CH4plant is estimated by adding the plant-mediated CH4 fluxes from all layers.

2.3 Parameters selected for optimisation

Parameter values related to the processes of CH4 emission in LPJ-GUESS are mostly adopted from the parameter values described in Wania et al. (2010). As Wania et al. (2010) had difficulties finding the optimal parameter values for many of the parameters, they performed some preliminary analyses for seven uncertain parameters for which there were little or no data available. They performed a simple initial sensitivity test by taking four sets of values for each of the seven parameters, followed by a parameter fitting exercise with three sets of values for each of the seven parameters. They ran the model with all 2187 different combinations for 7 sites for 1 year. As a result, they got an RMSE range of between 226.4 and 18.3 (mg CH4 m−2 d−1) for the different sites, which clearly indicates loosely fitted parameters with a high degree of uncertainty.

In this study, parameters for the optimisation are selected based on their sensitivity to the model output (CH4) and expert opinion. We used a simple method to calculate the percentage difference in output (single simulation) when varying only one input parameter at a time from its permitted minimum value to its maximum (Hoffman and Miller1983; Bauer and Hamby1991). The “sensitivity index” (SI) is calculated using the following equation:

(9) SI = D max - D min D max ,

where Dmin and Dmax represent the model output values corresponding to the respective minima and maxima of the corresponding parameter range. The values are taken based on expert opinion.

We considered 5 of the 7 parameters that Wania et al. (2010) tested in their sensitivity analysis (2 parameters related to the root exudate decomposition are not used in LPJ-GUESS) as well as 11 other parameters used in LPJ-GUESS. Based on their high SI values, we chose 11 of them for the optimisation (Fig. 2, Table 1).

Among the five eliminated parameters, agfrac is the fraction of annual net primary production (ANPP) used to calculate number of tillers, DTmin is the minimum temperature (°C) for heterotrophic decomposition, pororg is the porosity of organic material, Ccon is the carbon content of biomass, and U10 is the possible constant value of wind speed at 10 m height. Among the selected parameters related to the CH4 production, Rmoist and Rmoistan are the response of the soil moisture content to the soil organic matter decomposition in the acrotelm and catotelm, respectively (see Eq. 8 in Supplement); CH4/CO2 is the CH4-to-CO2 ratio under the anaerobic conditions (Eq. 2); fair is the fraction of air in peat (Sect. 2.2.1 and Eq. 2); poracro and porcato are the porosity in the acrotelm and catotelm, respectively (see Sect. S1); and λroot is the decay length of root biomass in peat (Eq. 1). Among the selected parameters related to the CH4 transportation, foxid is the fraction of oxidised CH4 (Sect. 2.2.2) related to the CH4 oxidation, wtiller is the average weight of an individual tiller, rtiller is the tiller radius, and ϕtiller is the tiller porosity (Eq. 8).

Figure 2Tested parameters for the optimisation and their SI values. Red and blue colours indicate an increase or decrease in the total CH4 flux, respectively, when the value of the parameter increases.


Table 1Selected model parameters for the assimilation. The prior values, units, prior standard deviation (SD), and parameter description are given.

Download Print Version | Download XLSX

2.4 Parameter optimisation framework

We assumed Gaussian PDFs to depict both the prior distributions of the parameters and the deviation between model and observations. The resulting model can be formulated as follows:

(10) Y | x N ( M ( x ) , R ) , x N ( x p , B ) .

Here, Y represents the observations; M(x) is the LPJ-GUESS output given parameters x; xp represents the prior values of the parameters; and R and B are error covariance matrices describing the uncertainty in observations and priors, respectively.

The prior uncertainties (B) are based on expert opinion and were kept relatively large to reduce the prior's influence on the posterior parameter estimates. We have assumed prior variance for each parameters as 40 % of their expected range (see Table 1). The parameters are also assumed to be a prior uncorrelated, due to the lack of good and consistent expert opinions regarding dependence.

2.4.1 Cost function

Using the Bayesian framework, the posterior for the parameters becomes

(11) P ( x | Y ) = P ( Y | x ) p ( x ) p ( Y ) P ( Y | x ) p ( x ) ,

which, on a log scale, results in the following quadratic cost function (Tarantola1987):

(12) log P ( x | Y ) = - J ( x ) + const. J ( x ) = 1 2 ( Y - M ( x ) ) t R - 1 ( Y - M ( x ) ) + 1 2 ( x - x p ) t B - 1 ( x - x p ) ,

where “const.” represents normalising constants not depending on the unknown parameters. The two terms in J(x) represent model–data misfit and the prior information on the parameters. A number of experiments aim to achieve the smallest cost function values to locate the optimal parameter set within the parameter space.

2.5 Adaptive Metropolis–Hastings (MH) algorithm

To search for the optimal posterior parameters, we used an MCMC-MH algorithm (Metropolis et al.1953a; Hastings1970a). A standard MH algorithm generates samples from a target distribution by, in each iteration, drawing from a proposal distribution and then either accepting the new state or copying the old state. For a Gaussian random walk, proposals (x^) are generated by adding a mean-zero normal random variable to the current value (xt):

(13) x ^ = x t + ϵ , ϵ N ( 0 , λ Σ ) ,

where λΣ is a scaling and covariance matrix describing the spread of the added random variable. The new value is accepted with a probability that compares the likelihood (or cost function) of the proposed and old sample:

(14) α = min 1 , P ( x ^ | Y ) P ( x t | Y ) = min 1 , exp ( - J ( x ^ ) + J ( x t ) ) .

If the value is accepted, set xt+1=x^, otherwise keep the previous value, xt+1=xt. The resulting sequence of states will represent dependent samples from the target distribution.

The adaptive MH algorithm used here contains three key concepts that are explained in the following sections: transformed proposals, providing a natural way of including limits on the parameters; the adaptive random walk, allowing the algorithm to estimate λ and Σ from the target distribution; and tempering of the target distribution, to reduce the effects of local maxima and allow better exploration of the target.

2.5.1 Transformed proposals

The standard proposal in Eq. (13) does not include any restrictions on the parameters. To handle parameter limits, we transformed the parameters, resulting in an adjusted random walk proposal:

(15) z t = g - 1 ( x t ) , z ^ = z t + ϵ , ϵ N ( 0 , λ Σ ) , x ^ = g z ^ .

A list of possible limits and corresponding functions are given in Table 2. Note that different functions can be applied to each parameter in x.

Having a transformed proposal requires an adjustment of the acceptance probability (Hastings1970a) in Eq. (14) to

(16) α = min 1 , p ( x ^ | Y ) p ( x t | Y ) i q x t ( i ) x ^ ( i ) q x ^ ( i ) x t ( i ) = = min 1 , exp ( - J ( x ^ ) + J ( x t ) + i log q x t ( i ) x ^ ( i ) q x ^ ( i ) x t ( i ) .

Here, xt(i) denotes the ith parameter in the xt vector, and the q terms are given in Table 2. Note that each transformation results in one adjustment for that parameter and that all adjustments have to be multiplied together.

Table 2Summary of transformations and corresponding adjustments to the acceptance probability for the three cases of variables with a lower limit, variables with an upper limit, and variables constrained to an interval.

Download Print Version | Download XLSX

2.5.2 Adaptive random walk

It has been shown that the optimal MH algorithm behaviour is obtained when about 20 %–30 % of samples are accepted (Gelman et al.1996; Roberts et al.1997). For the proposal in Eq. 13, this is achieved when Σ corresponds to the posterior covariance matrix of the target distribution, P(x|Y), and the scaling is λ=2.382/d, where d is the number of parameters in x. The key idea of the adaptive algorithms suggested in Andrieu and Thoms (2008) is to recursively estimate both Σ and λ from previous samples. An important note for the transformed proposal in Eq. (15) is that Σ and λ relate to the transformed variable zt and not xt.

A Rao-Blackwellised update of the covariance matrix will consider both the proposal, z^, and the previous value, zt, weighted according to the acceptance probability, α, computed in Eq. (16). The expectation and covariance matrix are updated recursively as follows:


where γt is an adaptation factor and we have used γt=t-0.51.

The global adaptive scaling then updates the scaling factor λ as follows:

(18) log λ t + 1 = ( 1 - γ t + 1 ) log λ t + γ t + 1 α - α target .

The update is on a log scale to ensure that λt stays positive; the second part of the equation compares the current acceptance probability with a target probability, αtarget=0.234, and adjusts λt to obtain this desired overall acceptance probability.

To limit the effect of initial values, we first ran 5000 steps of the chain without adaptation and with Σ as 10−3 times an identity matrix (i.e. very small initial steps). Thereafter, a covariance matrix Σ was estimated based on the initial samples, and the chain was run for a further 15 000 steps adapting only λ and not Σ. Finally, for the last 80 000 steps, both λ and Σ were updated as described in Andrieu and Thoms (2008). We call the resulting framework the global Rao-Blackwellised adaptive metropolis (GRaB-AM) algorithm.

2.5.3 Tempering the target distribution

For large numbers of data, the cost function J(x) can have very deep local minima, causing the MH algorithm to get stuck, even with the two already outlined adjustments. To reduce the scale of the cost function, we temper the target distribution (Jennison1993).

(19) P ̃ ( x | Y ) = P ( x | Y ) 1 / T = exp ( - J ( x ) / T ) ,

where T is a suitably large value.

Having run an MH algorithm for N samples, the first Nb samples are discarded as “burn-in”, and expectations or variances can be computed as averages of the remaining samples.

However, with a tempered target distribution, we have samples from P̃(x|Y) and need to use importance sampling to adjust for the difference in distributions (Jennison1993), resulting in weighted averages


where the weights are given by


and Jmin=miniJ(xi). The subtraction by Jmin is for numerical stability to avoid cases of 0/0 when all J(xi) are very large. For the case of no tempering, i.e. T=1, the weights simplify to wi=1/(N-Nb), resulting in unweighted averages.

2.6 Experimental design

We have designed a set of twin experiments and a real-data experiment. For both the twin and real-data experiments, we generated MCMC chains with a length of 100 000 samples, as mentioned in Sect. 2.5.2. MCMC approaches are computationally intensive and time-consuming. In this study, each model simulation took approximately 9 s to complete (using an AMD Ryzen Threadripper processor). As a result, for the 100 000 iterations, it consumed nearly 250 computational hours. However, it should be noted that this study involves the model running for a single site, and the computational speed is highly dependent on the performance of the processors being used.

Figure 3An example of posterior PDFs from the twin experiment. Prior and posterior distributions are illustrated with solid black and solid green lines, respectively. True parameter values (Ztrue), the prior mean, the posterior mean, and the MAP are shown in red, black, lime, and orange colours, respectively.


Twin experiment

Twin experiments are designed to assess the performance of the developed GRaB-AM and its ability to recover the parameter values. The daily CH4 output simulated by LPJ-GUESS using randomly chosen true parameter values (Ztrue) within their permitted range of variation is used as the synthetic observation. As the synthesised observation conforms completely to the model, any potential errors in the model have not influenced the parameter optimisation process, ensuring unbiased posteriors. The observations have been assigned uncertainties to resemble the real-data experiment, as described below. It is expected that the assimilated parameters will converge to the true values as the MCMC chain progresses in time. To freely recover the Ztrue values, the prior parameter values (xp) in the cost function (Eq. 12) are set as Ztrue. Two scenarios are considered for the twin experiment to test the identifiability of the parameters under different conditions: scenario 1, which has a shorter temporal scale from 2005 to 2014 (10 years), and scenario 2, which has a longer temporal scale from 1901 to 2015 (115 years). Scenario 1 is more realistic and is chosen to mimic the real data at Siikaneva, whereas scenario 2 constitutes an ideal, hypothetical case with observations over the entire simulation period. In each set of the scenarios, the optimisation started from a different initial point in the parameter space that was randomly selected from their prescribed ranges.

Real-data experiment

To estimate the posterior parameter values, an experiment using the real observation from Siikaneva is designed. The observed daily averages are compared with the model simulation in the cost function only when more than 90 % of the hourly observation are available each day. When there are gaps in the daily observation, we eliminate them and their corresponding modelled values from the cost function calculation. In principle, the error covariance matrix R should include both observation uncertainties and their correlations. As the latter is difficult to estimate, we neglected them, and the observation uncertainties are estimated as 30 % for the daily observations greater than 0.01 gC m−2 d−1 and a floor value of 0.3 for the observations less than 0.01 gC m−2 d−1.

2.7 Parameter value estimation

For all of the experiments conducted in this study, the first 75 % of the GRaB-AM chains are discarded as the burn-in. The PDFs generated after the burn-in are used to estimate each parameter's maximum a posteriori probability (MAP), posterior mean, and standard deviation (SD). Following the idea used in Braswell et al. (2005), the parameter distributions are grouped into three categories: “well-constrained”, “poorly constrained”, and “edge-hitting” parameters. The well-constrained parameters are those that exhibit a well-defined unimodal distribution with a low SD. The poorly constrained parameters are those that exhibit a relatively flat multi-modal distribution with a large SD. To be more precise with the estimation, for posterior parameter distributions that appeared multi-model, if the SD of the distribution was greater than 20 % of the total range, we classified them as poorly constrained. The edge-hitting parameters are those that cluster near one of the edges of their prior range.

2.8 Posterior resampling experiment

To examine the effect of parameter optimisation on flux components, we designed a resampling experiment from the posterior parameter distributions. From the experiment conducted using site observations, 1000 sets of parameters were randomly selected and used to run the model to simulate the CH4 flux components. The outputs from each simulation of the experiment were used to analyse the process correlations and process–parameter relationships.

3 Results

3.1 Twin experiment using GRaB-AM

We ran a set of four different twin experiments for each of the scenarios mentioned in Sect. 2.6. Each of them shows a reasonably good convergence, regardless of their chosen initial values. In scenario 1, all parameters except for CH4/CO2 and λroot showed good convergence (see Fig. S1 in the Supplement; posterior parameter correlations of the first experiment shown in Fig. S1 are provided in Fig. S2). The results of scenario 2 are not shown, as they followed the same pattern. The resulting PDFs of experiment 1 in scenario 1 after the burn-in are represented in Fig. 3, which shows the mean and MAP values as well as the SD of the parameters. The parameter values and related statistics are given in Table 3.

The parameter retrieval capacity of the GRaB-AM algorithm is estimated as the “retrieval score” by dividing the posterior mean estimates of the parameters from all of the chains in each scenario by the Ztrue parameter values. The idea behind the retrieval score is that, in an ideal case of complete recovery, the posterior parameter estimate and the Ztrue value are the same; hence, the retrieval score would be 1. Figure 5 shows the retrieval scores obtained for each parameter and their 1σ value. In scenario 1, the ϕtiller, poracro, porcato, wtiller, and λroot are well retrieved with a low SD. Scenario 2 performed better in parameter retrieval compared with scenario 1. In scenario 2, the majority of the parameters except the CH4/CO2, rtiller, wtiller, and λroot showed good retrieval scores but with comparatively high SDs (see Fig. 5). The overall mean retrieval score estimation is based on the ratio of the estimated and true values, given a value of 0.95 with an SD of 0.19 for scenario 1 and a value of 1 with an SD 0.21 for scenario 2 (see Fig. 5).

Figure 4Time series estimates of the twin experiments in scenario 1. The simulations, using four sets of posterior parameter values obtained from the twin experiments, are plotted against the twin observation used.


Table 3The posterior mean, MAP, and SD of retrieved parameters for selected twin experiments in scenario (Sc) 1 and 2 and the parameter classes estimated from analysing the distributions of all four chains. The parameter classes include well-constrained (WC) and poorly constrained (PC) parameters.

Download Print Version | Download XLSX

In general, the twin experiments have resulted in well-constrained and poorly constrained parameter classes. Examples of the different classes of the distributions for experiment 1 of scenario 1 are shown in Fig. 3. Based on the posterior distributions estimated from all four GRaB-AM chains, the parameters Rmoist, CH4/CO2, foxid, rtiller, fair, poracro, porcato, and λroot are well constrained in scenario 1, whereas the parameters Rmoist, CH4/CO2, foxid, rtiller, fair, poracro, porcato, wtiller, and λroot are well constrained in scenario 2 (see Table 3).

Table 4Cost function reduction observed in scenario (Sc) 1 and 2 of the twin experiments. The prior and posterior cost function values obtained from four sets of experiments for each scenario are provided. The misfit between observed and expected (zero) cost function values are represented as the reduced χ2 value.

Download Print Version | Download XLSX

The reduced posterior cost function values and their χ2 values are given in Table 4. Here, the reduced χ2 values are calculated by dividing twice the cost function by the number of observations used in the assimilation. Overall, the χ2 values indicate a statistically robust reduction in the cost function in all of the experiments, although a systematic behaviour of being below 1 is observed.

Figure 5Twin experiment results in terms of the mean retrieval score based on the ratio of the estimated and Ztrue values of the parameters. The horizontal red lines indicate a complete retrieval, and the error bar shows the SD obtained from different chains in different scenarios. Sc1 and Sc2 indicate the two respective scenarios.


3.2 Real-data experiments and optimised parameters

For the experiment with the real data, the observations collected at the Siikaneva wetland are assimilated using the GRaB-AM algorithm. The MCMC trace plots are exemplified in Fig. 6.

Figure 6An example of the GRaB-AM chains for the experiment with real observations showing all 100 000 values in the chain. The first 75 % of values were discarded as burn-in and are greyed out in the figures. The remaining 25 % (from the red vertical lines) of values were used for the analyses.


3.2.1 Optimised parameter values and distributions

The posterior parameter PDFs are shown in Fig. 7. In contrast to the twin experiments, the parameters fell into three categories: well-constrained, poorly constrained, and edge-hitting parameters. The classifications are given in Table 5. The PDFs for the Rmoist, CH4/CO2, ϕtiller, fair, poracro, wtiller, and λroot parameters are classified as well-constrained distributions; the PDFs for rtiller, porcato, and Rmoistan are classified as poorly constrained distributions; and that for foxid is classified as an edge-hitting distribution. In both the well-constrained and poorly constrained parameters, high kurtosis is observed. The values of foxid, which is the edge-hitting parameter, lay near the higher bound of the edges of the prior range, and most of the retrieved values were clustered near this edge. This parameter also exhibited large positive kurtosis and negative skewness. The posterior parameter values and their 1σ SDs along with the prior values are shown in Table 5. The MAP and posterior mean estimates agree on the value for the CH4/CO2, fair, and poracro parameters. For foxid, ϕtiller, and rtiller, both the MAP and posterior mean estimates stayed out of one-third of the 1σ range of the posterior distribution, which we consider a large difference. For the remaining parameters, the MAP and posterior mean estimates stayed within one-third of the 1σ of their posterior distribution; hence, we consider this a small difference. For the Rmoist and CH4/CO2 parameters, the posterior values appeared close to but below the prior values. The posterior values of Rmoistan appeared very close to but above the prior values. For the ϕtiller parameter, the MAP estimate appeared very close to but above the prior value, whereas the posterior mean estimate appeared very close to but below the prior value. For these four parameters, the posterior mean stayed within one-third of the 1σ range of the assumed prior uncertainty. The posterior values of the foxid, rtiller, and fair parameters appeared slightly above the prior values but outside of one-third of the 1σ range of the prior uncertainty. The prior and posterior values of the poracro parameter remained the same. In contrast, the porcato, wtiller, and λroot parameters appeared very distant from and below the prior values, out of one-third of the 1σ range of prior uncertainty, but stayed within the prior range (see Sect. 4.2 for details).

3.2.2 Posterior parameter correlation

The 2-D distributions of the posterior parameters and their correlations are illustrated in Fig. 9. Overall, the majority of the parameters showed weak positive or negative correlations with a few exceptions with extreme correlations (the values and corresponding colour code in the triangle above depict this). For example, Rmoistan showed a high negative correlation with Rmoist, and poracro showed a high positive correlation with fair. The 2-D marginal distributions (scatterplots), illustrated in the lower triangle, showed a general tendency toward high clustering within the 1σ range for all of the parameters; in general, the 1-D histograms (on the diagonal; also shown in Fig. 7) appeared as well-constrained unimodal distributions. For further details, the reader is referred to Sect. 4.2.1.

3.2.3 Cost function reduction

The prior and posterior cost function values corresponding to the prior and posterior parameter values are listed in Table 5. The prior cost function value calculated with the default model parameters showed a high-cost value of 48424.4, with a model overestimation of around 4 times the observed flux. After the optimisation, the cost function value was reduced to 2959.8 (reduced χ2=3.82) with the MAP estimate of parameters and to 3002.6 (reduced χ2=3.88) with the posterior mean estimate of parameters.

Figure 7Posterior PDFs of parameters from the GRaB-AM real-data experiment. The green curves shown are the smoothed Gaussian kernel estimates of the posterior distribution on the posterior histograms, whereas the black curves are the prior distributions. The dotted vertical black, lime, and orange lines represents the prior mean, posterior mean, and MAP, respectively. The shaded green area of the distributions represents the 1σ error estimate of the PDFs.


Figure 8Schematic summery of the resampling experiment. Panel (a) shows the process–parameter correlation and regression slope. Three different flux components of CH4 as well as the total flux are labelled on the vertical axis, whereas the parameters are labelled on the horizontal axis. The different colours of the circles represent the regression slopes (β) scaled between −1 and 1 (in 11 steps). The blue colour indicates a steeper negative slope (and hence a strong decrease) in CH4 fluxes with the increasing parameter value, whereas the red colour indicates steeper positive slopes (and hence a strong increase) in CH4 fluxes with the increasing parameter value. The coefficient of determination (R2) scaled between 0.05 and 1 (in 11 steps) is represented by the size of the circles, with larger circles indicating higher R2 values. Panel (b) shows the process–process correlations. Numeric labels on the upper triangle correspond to Pearson's correlation coefficient values. The diagonal of the matrix shows the 1-D histogram for each flux component and the total flux. The 2-D marginal distributions of the sum of the processes and total flux are represented in the lower triangle with contours to indicate the 1σ, 2σ, and 3σ confidence levels. The black dots in the plots indicate the sums of flux components. Ranges of the distributions are labelled on the left and bottom of the panel.


3.2.4 Flux components of the CH4 simulation and parameter values

To understand how and the magnitude to which each optimised parameter influences the flux components and the total flux, the result of the resampling experiment (see Sect. 2.8) is examined using correlation and regression analyses. The Pearson correlation coefficients and regression slopes are calculated for all 1000 parameter sets and their corresponding sums of the flux components and total flux. Figure 8a shows a schematic summary of the correlation coefficients and regression slopes between the 11 parameters, the flux components, and the total flux. For both the total flux and diffusion, all parameters, except foxid and ϕtiller, showed a similar regression pattern, although with slight differences in magnitude. This similarity is not surprising, as diffusion is the most dominant among the process components. The total flux showed the highest correlations with CH4/CO2 and λroot and the lowest correlation with foxid. A detailed discussion of the process–parameter relations can be found in Sect. 4.2.1.

Table 5Parameter values estimated after the GRaB-AM real-data optimisation. The prior values, MAP, posterior mean, posterior SD, and parameter classes are shown. The parameter classes include well-constrained (WC), poorly constrained (PC), and edge-hitting (EH) parameters. The cost function values corresponding to the parameter values obtained with the prior, MAP, and posterior mean estimates are also shown.

Download Print Version | Download XLSX

Figure 8b shows the correlations between the sums of flux components resulting from the resampling experiment. The 2-D distributions in the lower triangle show a strong positive relation between diffusion and total flux. Almost all of the parameter residuals are observed within 3σ deviation without many outliers. Except for the correlation between diffusion and total flux, the analysis showed no other strong positive or negative correlations between the components, as can be seen in the correlation plot illustrated in the top triangle of the aforementioned figure panel.

Figure 9A posteriori correlations between the parameters from the GRaB-AM real-data optimisation. The blue and red colours in the upper triangle represent strong negative and positive correlations, respectively. The numerical labels on the upper triangle are the Pearson's correlation coefficient values. The panels on the diagonal show the 1-D histogram for each model parameter with a dashed red vertical line to indicate the mode. The vertical blue lines are the 0.16, 0.5, and 0.84 quantiles of the distributions. The lower triangle represents the 2-D marginal distributions of each parameter with contours to indicate the 1σ, 2σ, and 3σ confidence levels, whereas the blue dots on the marginal distributions are the posterior parameter values. The ranges of the distributions are labelled on the left and bottom of the figure.


Yearly variations in the fractional contributions of flux components simulated using prior and posterior parameter estimates are examined to understand the impact of the optimisation on the composition of the inter-annual emissions. The time series of the annual sums of flux components as a function of their total flux (in %) are shown in Fig. 10a. The result shows that, among the flux components, diffusion contributes the most to the total CH4 flux both in prior and posterior estimates, with a slightly higher contribution in the posterior estimate, followed by plant-mediated transport (see Sect. 4.2.3 for detailed discussion).

Figure 10Flux component fractions and the model–data difference (in %). Panel (a) shows the proportions of annual flux components plotted as a function of the total yearly flux. The different flux components are represented with solid lines of different colours, and the symbols on the lines correspond to each flux component. Panel (b) shows the annual model–observation mismatch (in %) with respect to the yearly total CH4 observation.


The time series model–observation mismatch of prior and posterior estimates for the annual total fluxes can be seen in Fig. 10b; the values are displayed as a percentage of the observed CH4 flux. The prior estimate showed a mismatch of around 600 % for the first 2 years. Furthermore, a considerably high mismatch is observed in the years 2011, 2012, and 2014. The MAP estimate remained near 0, while the posterior mean estimate exhibited slightly negative values, indicating an underestimation of the flux. Interestingly, the MAP followed the same pattern as the prior estimation by showing an increase whenever the prior increased and a decrease whenever the prior decreased; however, the posterior mean estimate did not show this relation.

The percentage fractions of the annual errors in the flux components are shown in Fig. 11. The effect of optimisation on the individual contributions of each component can be seen from the annual means (solid dots) of their fractional contribution to the total flux. Among the prior estimates of flux components, the prior plant-mediated transport showed the largest error (22.5 %), whereas the ebullition showed the smallest error (9.1 %). In the MAP estimate, ebullition showed the highest error, with a value of 12.3 %, followed by diffusion and ebullition, with similar error values of 6.9 % and 6.8 %, respectively. For the estimate using posterior mean values, diffusion and plant-mediated transport showed around the same errors, 7.5 % and 7.4 %, and ebullition showed the least error (2.6 %). On the right-hand side of the aforementioned figure, the fourth column displays the mean and errors for the inter-annual variation in the total fluxes obtained by prior parameter values and posterior estimates. The prior total estimate showed an error of 4.2 %, while the mean and MAP showed an error of 0.66 % and 0.72 %, respectively.

Figure 11The first three columns of the figure show the fractions of the annual fluxes from process components of the total fluxes. The vertical solid lines represent the 1σ error bars of each component, whereas the dots represent the mean of the annual fluxes. The fourth column (which corresponds to the right y axis) shows the annual mean and annual errors for the inter-annual variation in the total fluxes.


Figure 12Total CH4 simulation from the LPJ-GUESS model (red dots) after optimising with the GRaB-AM algorithm. The black dots are the real CH4 observations from Siikaneva with prior observation error (grey shading). The light red shading around the posterior model simulation is the 95 % confidence interval (CI) of the simulations. The blue dots are the prior simulation with the prior default model parameters. A few outliers above 0.3 gC m−2 on the vertical axis have been removed from the figure for better visualisation. While most of the observations fall within the confidence intervals, it is important to note that the effects of parameter variations in the posterior are part of these confidence intervals.


3.3 Fit to the observation

Figure 10b illustrates the percentage model–data misfit, and Fig. 12 shows the time series of the assimilated observations as well as the model prior and posterior estimates with their uncertainties. The total RMSE estimated between the prior and observations was 0.044 gC m−2 d−1, which was reduced to a value of 0.023 gC m−2 d−1 for the posterior case. This result indicates that most of the mismatch between the prior model estimates and observations was contributed by the large overestimation in the initial years. This overestimation disappeared in the posterior, showing a better agreement with the observation. There are years for which the observations show large peaks during the summer (such as 2010, 2012, and 2013), and the posterior estimates succeeded in capturing these peaks to a large extent, although not completely (see Sect. 4.5 for details).

4 Discussion

4.1 Twin experiment

A common problem with the adaptive MH algorithm is its inability to widely explore the target distribution if the set-up is not well tuned. This can then result in a poor approximation of the target distribution and, hence, poor adaptation. The resulting trace plots shown in Figs. 6 and S1 depict a set of well-explored parameters in their permitted space ranges during the progression of the random walk, indicating a well-tuned assimilation set-up. The use of adaptive Rao-Blackwellised learning of the posterior distribution appeared beneficial during the transients of the chains whenever the acceptance probability was dropped to low values in low-probability regions of the parameter space.

Even though twin experiments have retrieved all true parameters in different experiments, individual experiments reveal instances where one or two parameters may not fully converge to their true values due to factors such as parameter correlation or equifinality (Figs. 5 and S1). Given the complexity and non-linearity of the model, it is unsurprising that not all parameters converged completely. It is also unsurprising that different chains estimated slightly different posterior solutions for the parameters. However, most poorly retrieved parameters still have their true values within the 1σ range of the Gaussian PDFs of the optimised values. Even when the parameters are slightly off from the Ztrue values, Fig. 4 shows the capability of the twin experiments with respect to capturing the structure of the observations, including the observed spikes.

The occasional lack of convergence mentioned above likely introduces a slight bias in the posterior estimates, i.e. compensating effects due to the equifinality in subsequent real-data analyses, leading to suboptimal parameter estimates. This equifinality is essentially based on the high non-linearity of the model and a missing observational data constraint that would help distinguish between certain parameter combinations. Nevertheless, in contrast to other similar studies, like Santaren et al. (2014), which have encountered challenges in achieving complete recovery of many true parameter values, our GRaB-AM algorithm showed better performance with respect to recovering true parameter values. The analysis of the cost function reduction (Table 4), the ability to constrain the parameters (Fig. 3), the ability to capture the structure of the model (Fig. 4), and the parameter retrieval ability of the twin experiments (Fig. 5) showed that the developed GRaB-AM algorithm is still capable of optimising the process parameters related to CH4 emissions in LPJ-GUESS, given the caveats mentioned above on the equifinality.

Similarly, the systematic low χ2 values observed for the twin experiment do not necessarily affect the framework's ability to be set up for the real-data experiment (Sect. 3.1). As the twin experiments here are under the assumption of an “idealised model”, meaning that the model perfectly reproduces the observations without any errors or uncertainty, and “error-free data”, where the data perfectly represent the environmental conditions without any random errors, it is expected to have χ2 values systematically below 1. Furthermore, the χ2 value is highly sensitive to the number of observations and parameters. Having 3650 observations in scenario 1, 41 610 observations in scenario 2, and only 11 parameters can lead to low χ2 values. However, the comparatively smaller χ2 values for sets 1 and 3 in scenario 1 and set 3 in scenario 2 indicate a tendency toward overfitting the results and being overconfident in the estimated posterior values and uncertainties.

4.2 Parameter estimation using real observations

As described in Sect. 3.2.1, the experiment using real-data resulted in three poorly constrained parameters and one edge-hitting parameter. Poorly constrained or edge-hitting parameters, however, are not uncommon in MH parameter searches and are somewhat expected with a complex and highly non-linear model such as LPJ-GUESS. The correlation of parameters with other parameters can affect the result, i.e. the number of parameters that can be optimised within this data assimilation framework is limited. Although the twin experiments showed good parameter retrieval capabilities, the assimilation of real-world observations into a complex ecosystem model like LPJ-GUESS is expected to have parameter retrieval and equifinality problems. This is another reasons for selecting a small subset of the parameters associated with wetland CH4 flux simulations for this study. Considerable changes have occurred in the prior parameter values after optimisation. Here, it should be considered that, in this study, the assimilation aims to reduce the magnitude of the prior CH4 flux simulation to minimise its misfit with the observed data, which is nearly half of the prior model estimate (see Table 6).

4.2.1 Posterior correlation estimates

The following discusses the possible impact of the posterior parameters (shown in Table 5) on CH4 flux simulation, the interactions between the optimised parameters and the component fluxes (shown in Fig. 8), and the parameter–parameter correlations (shown in Fig. 9). (We distinguish between strong, >0.5, and weak, <0.2, parameter correlations but focus on the strong ones here.)

The very slight reduction, i.e. within one-third of the 1σ error, observed in the posterior mean estimate of Rmoist (Table 5) indicates a slight decrease in the moisture response under aerobic conditions, which would likely result in a slower soil carbon turnover time with a slight decrease in CH4 emission. The weak R2 value and the weak positive β value of Rmoist with all of the flux components indicate that a decrease in this parameter value decreases the emission and explains some of the variances in the flux components (Fig. 8).

In contrast, Rmoistan obtained a higher posterior value compared with the prior (Table 5) with a slightly asymmetric multi-modal distribution (see Fig. 7), indicating an increase in the moisture response in the anaerobic catotelm. Together with this, the strong negative correlation (−0.8) observed between Rmoist and Rmoistan (Fig. 9) indicates reduced decomposition in the acrotelm and increased decomposition in the catotelm. Rmoistan had a positive effect on diffusion and a negative effect on plant-mediated transport (Fig. 8). An increase in Rmoistan could enhance CH4 production in the catotelm. As the catotelm has a low plant root abundance, this increase would lead to more diffusion and a reduction in total plant-mediated transport. The increase in Rmoistan contributed very little to ebullition. This is most likely due to the negligible contribution of ebullition to the overall flux, with no contribution most of the time.

The posterior CH4/CO2 parameter, which is the CH4-to-CO2 ratio in an anaerobic environment, was found to be lower compared with the prior (Table 5). This indicates a high fraction of CO2 production from the peat compared with CH4 production. The very high R2 value of CH4/CO2 for diffusion and plant-mediated transport (which represent the two diffusive pathways) indicates that a significant portion of the variance in these pathways can be explained by this parameter (Fig. 8). Similarly, the high, positive β value for diffusion and plant-mediated transport indicates a substantial linear increase in emissions via these pathways if the parameter is increased. The increase in ebullition is marginally less than the other fluxes, most likely because ebullition is limited by the availability of the gaseous fraction of CH4. The dissolved CH4 will first emit via diffusive fluxes; hence, there could be very little CH4 left in the gaseous phase for ebullition. The CH4/CO2 ratio is negatively correlated with the Rmoistan and λroot parameters (Fig. 9). This indicates a lower CH4 fraction produced by decomposition in deep soil.

The prior parameter value for fair was 0, which means that there is no “permanent” gas fraction in peat (Table 5). The posterior value of fair was slightly positive (0.032), indicating a small air fraction in the peat. fair showed a very high positive correlation with poracro, which can simply be explained by the fact that more porous soil allows for more air in the soil (Fig. 9). An increase in the fair value would increase all of the flux components, with a notably larger effect on diffusion (Fig. 8). As stated in Sect. 2.2.3, the diffusivity of CH4 in air is 4 orders of magnitude greater than that in water, indicating that a higher fraction of air in the soil results in the rapid and easy transport of CH4 to the atmosphere. The larger increase in diffusion can be directly attributed to fair, as this parameter directly controls diffusion.

The fraction of available oxygen utilised for CH4 oxidation is determined by the foxid parameter, which has a higher value after optimisation (Table 5). The high values of foxid and fair, indicating a high available air fraction and, hence, high O2 concentration in the soil, result in the conversion of most of the available carbon into CO2. This could explain the above-mentioned reduction in CH4/CO2 as a balancing effect (Eq. 2). foxid showed a negative β value for diffusion and ebullition and a slight positive β value with a comparatively high R2 value for plant-mediated transport (Fig. 8). A decrease in ebullition can be explained by the increased availability of oxygen for CH4 oxidation, resulting in less CH4 being emitted via ebullition. A significant decrease occurs in diffusion because the diffusive flux cannot bypass the top layer, into which oxygen diffuses. Directly explaining the increase in plant-mediated transport is difficult due to the complex process formulation in the model; however, it can be assumed that the aerenchyma could transport a part of the oxygen deep down to the soil layers where it plays less of a role in oxidation but contributes more to the total gas pressure, which can escalate the passive plant-mediated transport to the atmosphere.

The optimisation of plant-related parameters depends on the specific plant species present in the wetland. A slight reduction in the posterior mean estimate of ϕtiller suggests that the tillers may be slightly more compact, with reduced porosity for CH4 transport. A considerable reduction, more than one-third of the prior uncertainty, is observed in wtiller, indicating lower leaf biomass (Table 5). This reduction in potential leaf cover would lead to less carbon added to the potential carbon pool for methanogens, resulting in lower CH4 emissions. A decrease in tiller weight would add less organic carbon to the soil, leading to a less-compact peat accumulation in the bottom soil layers with more porosity for water. This could explain the negative correlation between wtiller and porcato (Fig. 9). In contrast to the values of the two above-mentioned CH4 transport-related parameters, rtiller, which represents the tiller radius of plants, showed a value more than twice the prior value (Table 5). This increase indicates more cross-sectional area of tiller for a given biomass, resulting in an increase in plant-mediated CH4 transport (see Eq. 8). These three parameters related to plant-mediated transport showed strong positive correlation with each other. They also exhibited positive β values in relation to plant-mediated transport (Fig. 8). These parameters can have two effects on emissions: on the one hand, having aerenchyma cells with more porous space, radius, and biomass enhances CH4 transport to the atmosphere; on the other hand, via the same spacious aerenchyma cells, it is also possible for plants to transport more O2 to the soil. This enhanced O2 transport to the soil could explain the slight reduction in diffusion and ebullition observed in the cases of ϕtiller and wtiller.

The posterior value for the porosity in the catotelm (porcato) was significantly lower than the prior (Table 5), suggesting a more compact catotelm with less water (as it is assumed to be saturated). This change could have a dual effect on CH4. Variations in water content can slightly affect soil temperature, potentially leading to an increase in the flux if the temperature rises or to a decrease in the flux due to the compact peat. As described in Sect. 3.2.1, poracro remained unchanged, indicating no changes in acrotelm porosity (Table 5). The positive kurtosis observed in the PDF of this parameter indicates a well-constrained single solution, while the negative skewness indicates a more probabilistic region below the posterior estimate. Similar to fair, poracro also exhibited positive β values for all flux components, although with a relatively low R2 value (Fig. 8). This positive relationship may be attributed to the increased presence of air in the acrotelm soil, which could facilitate CH4 emissions. In contrast, an increase in porcato could lead to a slight reduction in ebullition. This could be because more water can potentially occupy the pores of permanently saturated catotelm which will indirectly affect ebullition through phase change and by affecting soil temperature.

The posterior value for λroot is estimated to be much smaller than the prior (higher than the value reported in Wania et al.2010, and in Susiluoto et al.2018), i.e. more than one-third of 1σ of the prior estimate (Table 5). This small posterior value for λroot indicates a low decay length of root biomass in the soil, meaning that more of the decomposition and CH4 production occurs in the acrotelm and that less occurs in the catotelm. The emission of CH4 produced mainly by peat decomposition in the acrotelm would be facilitated by a low posterior value for λroot, with around 60 % in the first layer of the acrotelm followed by 22 % and 8 % in the respective second and third layers of the acrotelm. λroot played a key role in this optimisation. Figure 8 shows that λroot has a strong negative β value for diffusion and a weak positive β value for plant-mediated transport, both with relatively strong R2 values (Fig. 8). As most of the peat decomposition activities are assumed to occur in the acrotelm, the reduction in the magnitude of λroot could facilitate diffusion, especially as it is the largest component. On the other hand, plant-mediated transport may be reduced, as it limits the distribution of roots in soil depths.

4.2.2 Posterior flux components

In Fig. 13, the time series of process components is shown for the posterior mean estimate. In general, the optimisation of the model parameters leads to an approximately 50 % decrease in the production of CH4 compared with the prior, with a significant reduction in the plant-mediated and ebullition components, leaving diffusion as the dominant component. Diffusion is reduced by around 30 %, and plant-mediated transport is reduced by approximately 86 %.

The low contribution of plant transport is mainly due to the reduced value of the root-depth-controlling parameter λroot, which decreased from 25.17 to 10.58. This smaller proportion of plant-mediated transport is somewhat surprising for a fen wetland site like Siikaneva, which features a significant aerenchymatous leaf area throughout the growing season. The result is contradictory to the results obtained from optimising the HIMMELI model (Susiluoto et al.2018), in which the largest fraction of CH4 is contributed by the plant-mediated transport. However, field experiments conducted to estimate plant-mediated transport by Korrensalo et al. (2022) have observed a smaller proportion of the ecosystem-scale CH4 flux attributable to plant CH4 transport at the Siikaneva fen site. This observation aligns well with the results that we obtained.

The largest reduction, however, was for ebullition (by around 92 %). This result is not surprising, as Wania et al. (2010), who provide the basic foundation of the CH4 model in LPJ-GUESS, also reported almost virtually no ebullition to the surface at several sites. Figure 13 shows that no ebullition is estimated during the years 2008, 2010, and 2012. Here, it should be considered that the representation of ebullition in LPJ-GUESS is somewhat simplified using a curve-fitting equation to calculate the solubility and using the ideal gas law to convert the volume of CH4 per volume of water into the corresponding number of moles. Due to this lack of detail and its fast timescale occurrence (mostly depends on the physical parameters such as temperature and pressure) and with no relevant parameters in the control vector, the optimisation could not alter the ebullition component directly. However, ebullition is indirectly controlled by parameters related to CH4 production and transport when there is a high amount of saturated CH4 available in the soil water. Thus, the optimisation can indirectly affect the ebullition. The total observed CH4 flux from Siikaneva during the period of 2005 to 2014 was 56.0 gC m−2, while the prior model estimate was 98.5 gC m−2 (Table 6). After the optimisation, with the posterior mean estimate of parameter values, the model estimated a flux of 53.5 gC m−2 with an estimated posterior uncertainty of ±4.82 gC m−2. This shows a reduced model–data error after optimisation with a difference of only 2.5 gC m−2.

Figure 13Time series for diffusion, ebullition, and plant transport using parameter values from the posterior mean estimate. A few outliers above 0.08 gC m−2 on the vertical axis have been removed from the figure for better visualisation.


Table 6Total emissions from flux components estimated from the MAP, posterior mean, and prior parameter values for all 10 years of the optimisation, given in grams of carbon per square metre.

Download Print Version | Download XLSX

4.2.3 Posterior process–process correlation

After the optimisation, the air fraction in the peat was increased, which is likely the cause of the enhanced diffusion. Diffusion is estimated in the model based on the soil porosity, water, temperature, and air fractions in the soil. Correlating the diffusion with the ebullition showed a negative result, i.e. illustrating the dominance of diffusion over ebullition (see Fig. 8b). A larger air fraction in the soil can also lead to an increase in plant-mediated emission, as the passive diffusion of air through the plant tissues depends on the amount of air in the soil/peat water (see Sect. 2.2.3). This can be seen in Fig. 8b as a comparatively high correlation between diffusion and plant-mediated transport. The increased tiller radius (rtiller) in plants increases the Atiller value (Eq. 8), thereby favouring faster diffusion through the aerenchyma cells. Ebullition is positively correlated with plant-mediated transport, indicating the occurrence of both of these components when there is a high concentration of CH4 in the soil. This occurs when the water table is located close to the surface and when there is a higher density of graminoids. An increase in plant-mediated transport of gases to the soil increases the net pressure imposed by the gases in soil/peat water, which could also lead to increased ebullition.

4.3 Model error and fit to the observation

Except for ebullition, all of the prior process components exhibited larger variances in the annual errors compared with the posterior estimates (Fig. 11). The plant-mediated transport is the component with the largest error in the prior estimate. The posterior error estimates for this component showed nearly equal values, with a slightly higher value for the posterior mean estimate. A similar pattern can also be seen for diffusion. In contrast, the MAP error estimate for ebullition showed a higher value compared with the posterior mean error but interestingly also to the prior. The posterior mean error estimate for ebullition showed the lowest value.

The annual sums of flux components mentioned above are illustrated in Fig. 10a. It is clear from this figure that the prior process components had large inter-annual variance, especially for the first 3 years and the last year. A considerable reduction in variance is observed for both the MAP and posterior mean estimates. The reduction in the variance observed in posterior estimates is not proportional to the prior; however, the posterior estimates showed a comparatively high variance in the first and last years. In Fig. 10b, as described in Sect. 3.2.4, the posterior mean estimate shows a comparatively high variance (with respect to the MAP estimate) for the annual errors with a negative bias throughout the time period. In contrast, the MAP estimate showed a positive bias throughout the time period. Compared with the posterior mean estimate, the MAP estimate has considerably larger parameter values for ϕtiller and rtiller; this could possibly be interpreted as slightly more CH4 emission through the increased tillers of plants and, hence, the reason for the positive bias in the MAP estimate. Figure 11 also indicates a high percentage of annual plant-mediated emissions for the MAP estimate. The negative bias in the posterior mean estimate could be due to the additional wintertime emissions from the real-world wetlands, which are not captured in the model. In the model, the emissions start around early summer, once the soil is not frozen anymore. In addition, the large daily variability in the observations of the summertime fluxes is also not represented in the model. Overall, the posterior estimates of the annual fluxes are in good agreement with the observations, leading to a small model–data mismatch for both MAP and the posterior mean estimates.

4.4 Model inputs and uncertainty

As mentioned in Sect. 3.3, a somewhat pronounced systematic underestimation of emissions was observed in the years 2010, 2012, and 2013. None of the twin experiments exhibited these systematic errors, which indicates that the issue could be attributed to a structural model error (see Fig. 4). While the CH4 module within LPJ-GUESS is relatively comprehensive when compared with many other similar models, the model's process description and parameterisation remain incomplete. For instance, in the real world, wind plays a crucial role in CH4 emissions and its atmospheric concentration; however, wind speed is set to 0 for modelling convenience in LPJ-GUESS, thereby presenting a significant limitation. Similarly, the lack of the representation of air pressure, the simplified representation of ebullition (see Sect. 4.2.2), and the simplified representation of CH4 production (see Eq. 2) are also a major limitations. Another reason for this mismatch could be the variations in the input climate data. The correlation plots of the input environmental variables of LPJ-GUESS and the CH4 residuals (Fig. S3) indicate that, in these years, CH4 emissions showed a comparatively high correlation with swr and air temperature, although precipitation did not show any significant relation. The results of the sensitivity study indicate that both the prior and posterior model estimates are significantly sensitive to the input variables. However, the posterior model estimate exhibited a considerable reduction in sensitivity, especially for swr and precipitation (see Sect. S5 and Fig. S5).

As mentioned in Wania et al. (2010), the flux components are determined by complex processes that depend on changes in many environmental factors. The model is unable to represent peak emissions caused by these microenvironmental changes. For instance, ebullition (one of the more complex CH4 emissions processes in LPJ-GUESS, as explained in previous sections) depends on the volumetric content of wind and various gases as well as on hydrostatic and atmospheric pressure. However, the model does not use them as forcing variables. Ebullition is also affected by the concentration of CH4 and the density of nucleation sites, which are difficult to represent in the model.

It should be noted that the incomplete state vector used for optimisation might have affected the optimised model result. The process representations in LPJ-GUESS are complex and interconnected, with a multitude of parameters directly or indirectly linked to CH4 fluxes. Representing the indirect parameters can be intricate, as they may depend on other fluxes or model components. For instance, the soil module in LPJ-GUESS is intricately connected to the CENTURY model, featuring 10 soil compartments. Introducing a parameter related to soil temperature or the wtd into the framework would necessitate consideration of the intricacies of the CENTURY model. Furthermore, it might require the inclusion of additional flux species, such as net primary production (NPP), soil temperature profiles, or wtd. This would significantly increase the complexity of the problem, exceeding the scope of this paper. Given these caveats, the small negative biases obtained for the posterior mean estimates when compared against the observations (see Fig. 10b) are reasonable considering the quality and uncertainty in the input data used (see Sect. 2.1) and the complexity and structural issues of LPJ-GUESS.

4.5 Optimised simulation from LPJ-GUESS

A detailed time series distribution of prior and posterior model simulations plotted against the observations is shown in Fig. 12. The posterior model predictions were adjusted by the optimisation to fit the observations with considerable adjustment to the summer peaks. For example, the large peaks in the modelled emissions in 2005 and 2006, which largely contributed to the prior cost function, disappeared in the posterior emissions. In the following years, 2007 and 2008, the prior model simulations underestimated the observation, which also got corrected in the posterior. Furthermore, the posterior emissions largely capture the comparatively high peaks in the observations for the years 2010 and 2012, although the model still underestimated the observations. In 2013, the observations were high and the optimisation failed to capture this peak; rather, it tried to compensate for the underestimation by releasing a sudden high spike at the end of the summer that year. In winter months, the model simulated very low or almost zero fluxes (as discussed before), whereas, the observations showed a small emission (around 8.3 % of the assimilated total), often with some small spikes (possibly from the ebullition). This inability of the model to capture the wintertime emission has contributed to the posterior model uncertainty and model–data misfit.

A significant mismatch in soil temperature and wtd has been observed between the model and observations, especially during the wintertime (see Sect. S4 and Fig. S4). The model tends to overestimate wintertime temperatures and underestimate the wintertime wtd, indicating completely frozen soil with a very low wtd. This discrepancy could be a reason for the complete suspension of model CH4 emissions during the winter. Observations show many days with a wtd above ground level, both in summer and winter. As the wtd is a key factor that can affect all flux components, this error could contribute to the misfits observed after the optimisation, although it cannot explain the systematic under estimation observed in the posterior. It should be noted that there are no considerable differences observed between the prior and posterior soil temperature or wtd.

As discussed in Sect. 4.2.3, the contribution of ebullition to the posterior estimate is comparatively negligible. Compared with the posterior, there were many emissions spikes observed in the prior estimate, especially during the beginning and the end of the summer months. Apart from these spikes, the prior CH4 estimates during the summer were a bit low in most of the years. The posterior estimate has considerably reduced these high spikes and adjusted the summer peaks to match the observations better. On the other hand, while compromising for the summer peaks in the observations, the model with the optimised parameters often failed to capture the abrupt high fluxes in the daily observation and simulated them at slightly wrong times. The spike shown at the end of 2013 is an example of this mistiming. This is likely caused by errors in the meteorological input data and missing wind and pressure representations.

It can be seen from the Fig. 12 that the majority of the observations lie within the 95 % confidence interval of the posterior estimate. Often, the observation uncertainty overlaps the confidence interval, except for the summer peak times of 2010, 2012, and 2013 in which the observations showed strong peaks compared with the average values. The few outliers in the observations are not captured by the model; these could likely be measurement artefacts and/or due to environmental forcing not considered here, such as wind speed or air pressure.

4.6 Merits and shortcomings of the GRaB-AM framework

The advantages of the developed GRaB-AM framework include the ability of the MCMC method to escape local maxima or minima, making it more robust in this respect than gradient-descent-based methods. Furthermore, the MCMC method is derivative-free, thereby avoiding issues with computing gradients for very rough (non-continuous) functions. The adaptive part of the MCMC “learns” about parameter correlations and utilises these dependencies in the proposal, allowing it to better explore the parameter space. Potential issues for the MCMC are that very uneven cost functions can lead to the chain getting stuck (essentially, the local minima is too deep for the algorithm to escape). Here, we have alleviated this issue by tapering the cost function. We also acknowledge that the cost function could be improved, as the squared cost function essentially assumes Gaussian errors with equal variance for all values and ignores any temporal correlations. The observations are measured as time series, and concentration/flux observations are known to have spatio-temporal error correlations. Adapting the cost function to account for these factors would be of interest and is likely to have a larger impact on the posterior parameter uncertainties than on the estimates of the parameters.

In principle, the twin experiment should be able to recover all of the parameters completely. However, the twin experiment in this study did not completely recover the true values of some parameters, especially the CH4/CO2 parameter. This could be attributed to one or more of the following issues: (1) due to the complex and extremely non-linear nature of LPJ-GUESS, there is a high possibility of equifinality; (2) the high dimensionality of the problem and limited variability in the data could be a reason for the poor convergence of some parameters; and (3) the model is not fully constrained by the limited dimensions in both the data and the parameter space. Achieving complete convergence may require incorporating (a) additional complementary observations capable of constraining additional model processes and (b) additional parameters from different modules of the model that represent different model processes.

The failure observed in convergence is expected to yield suboptimal parameter estimates with overly optimistic posterior uncertainty estimates due to the equifinality problems in the real-data experiments. It may also result in slightly corrupted relationships between variables, leading to changes in the covariance structure of the posterior parameters estimated from the real data. Considering the high degree of non-linearity in the model, a single type of twin flux from a single site might not be sufficient to identify all of the parameters. Therefore, additional sites spanning a large climatic variability and/or different species of twin fluxes from the same site may be necessary. For instance, regarding the CH4/CO2 parameter, one reason for the poor recovery could be that the study only optimises the CH4 component, leaving the CO2 fluxes significantly off. Incorporating both CH4 and CO2 fluxes to address both sides of the ratio would lead to a more accurate convergence. However, a challenge in this approach would involve formulating cost functions for different observations in a manner that properly weights them to equally represent the information without diminishing their significance. These considerations are beyond the scope of this paper.

5 Conclusions

This study represents an initial effort to optimise the model process parameters controlling the simulation of wetland CH4 fluxes within the LPJ-GUESS model using a Rao-Blackwellised adaptive MCMC technique based on Bayesian statistics. The assimilation framework has been shown to be able to retrieve true parameter values within the uncertainty limits discussed above by performing a set of twin experiments. Furthermore, we used eddy-covariance flux measurement data from a boreal wetland to calibrate the LPJ-GUESS model parameters for a site-specific simulation. The results demonstrate that the fit to the observation of the CH4 simulation of a complex terrestrial DGVM like LPJ-GUESS can be systematically enhanced with a Bayesian parameter calibration. The results also show that the modelled processes and the estimated parameters were well constrained by the observations, leading to a substantial reduction in the posterior uncertainty in the simulated CH4 emissions. However, we note here that, due to the equifinality problems identified in the identical twin experiments, the posterior parameter values contain compensating effects caused by the non-linearity of the model and the lack of an additional observational data constraint to uniquely identify each individual parameter. The results of the resampling experiment and the parameter and process correlations indicate that there are no redundant processes in the model description.

The robustness of the assimilation framework developed in this study calls for further application of the framework using observations from multiple sites in a simultaneous assimilation. Further validation of the framework's performance is necessary to confirm its applicability to other sites with diverse PFTs and climatic conditions. The relatively strong roughness in the shape of the cost function observed in this study is expected to be reduced in a multi-site assimilation experiment, as has been observed by Kuppel et al. (2012); this would provide additional observational data constraints to allow the retrieval of the global minimum of the cost function more easily and reduce the equifinality problem. These further applications are beyond the scope of this paper and will be investigated in future studies.

Code and data availability

The AMCMC code and CH4 flux data used for this article are available at (Kallingal et al.2022). The LPJ-GUESS model code can be obtained at (Nord et al.2021).


The supplement related to this article is available online at:

Author contributions

JTK designed the GRaB-AM assimilation framework with help from MS and JL. PAM developed the CH4 model in LPJ-GUESS based on the work of Wania et al. (2010). MS and PAM provided knowledge and advice about peatland CH4 processes as well as model parameters. JR and MR helped with collecting the CH4 observations. MS and JL provided supervision and technical advice with respect to programming the algorithm and analysing results. JTK prepared the manuscript with contributions from all co-authors.

Competing interests

At least one of the (co-)authors is a member of the editorial board of Geoscientific Model Development. The peer-review process was guided by an independent editor, and the authors also have no other competing interests to declare.


The views and opinions expressed are those of the author(s) only and do not necessarily reflect those of the European Union. Neither the European Union nor the granting authority can be held responsible for them.

Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims made in the text, published maps, institutional affiliations, or any other geographical representation in this paper. While Copernicus Publications makes every effort to include appropriate place names, the final responsibility lies with the authors.


The authors are grateful to Guillaume Monteil, Adrian Gustafson, Stefan Olin, and Johan Nord for discussions on the computational set-up and the manuscript. The Institute for Atmospheric and Earth System Research (SMEAR-INAR) at the University of Helsinki is acknowledged for the CH4 flux data collected at Siikaneva. We would also like to acknowledge the Finnish Meteorological Institute (FMI), Finland; the Research Council of Finland (grant nos. 336573, 341753, and 325169); and ICOS-Finland for open-data availability.

Financial support

This research has been partly supported by the Strategic Research Area “Biodiversity and Ecosystem services in a Changing Climate” (BECC), Lund University (grant no. DnrV 2018/467), and is a contribution to the Strategic Research Area “ModElling the Regional and Global Earth system” (MERGE). BECC and MERGE are funded by the Swedish government. This research has also been partly supported by the “Greenhouse Gas Fluxes and Earth System Feedbacks” (GreenFeedBack) project from the European Union’s Horizon Europe Framework Programme for Research and Innovation (project no. 101056921), funded by the European Union.

Review statement

This paper was edited by Julia Hargreaves and reviewed by two anonymous referees.


Aalto, J., Aalto, P., Keronen, P., Kolari, P., Rantala, P., Taipale, R., Kajos, M., Patokoski, J., Rinne, J., Ruuskanen, T., Leskinen, M., Laakso, H., Levula, J., Pohja, T., Siivola, E., Kulmala, M., and Ylivinkka, I.: SMEAR II Hyytiälä forest meteorology, greenhouse gases, air quality and soil, University of Helsinki, Institute for Atmospheric and Earth System Research,, 2022. a

Ahti, T., Hämet-Ahti, L., and Jalas, J.: Vegetation zones and their sections in northwestern Europe, in: Annales Botanici Fennici, pp. 169–211, JSTOR, 1968. a, b

Alekseychik, P., Peltola, O., Li, X., Aurela, M., Hatakka, J., Pihlatie, M., Rinne, J., Haapanala, S., Laakso, H., Taipale, R., Matilainen, T., Salminen, T., and Levula, J.: SMEAR II Siikaneva 1 wetland eddy covariance, University of Helsinki, Institute for Atmospheric and Earth System Research,, 2019a. a

Alekseychik, P., Kolari, P., Rinne, J., Haapanala, S., Laakso, H., Taipale, R., Matilainen, T., Salminen, T., Levula, J., and Tuittila, E.-S.: SMEAR II Siikaneva 1 wetland meteorology and soil, Universiy of Helsinki, Institute for Atmospheric and Earth System Research,, 2019b. a, b

Andrieu, C. and Thoms, J.: A tutorial on adaptive MCMC, Stat. Comput., 18, 343–373, 2008. a, b, c

Aurela, M., Riutta, T., Laurila, T., Tuovinen, J.-P., Vesala, T., Tuittila, E.-S., Rinne, J., Haapanala, S., and Laine, J.: CO2 exchange of a sedge fen in southern Finland-The impact of a drought period, Tellus B, 59, 826–837, 2007. a

Aurela, M., Lohila, A., Tuovinen, J.-P., Hatakka, J., Riutta, T., and Laurila, T.: Carbon dioxide exchange on a northern boreal fen, Boreal Environ. Res., 14, 699, 2009. a

Bauer, L. and Hamby, D.: Relative sensitivities of existing and novel model parameters in atmospheric tritium dose estimates, Radiat. Prot. Dosim., 37, 253–260, 1991. a

Bousquet, P., Ciais, P., Miller, J., Dlugokencky, E. J., Hauglustaine, D., Prigent, C., Van der Werf, G., Peylin, P., Brunke, E.-G., Carouge, C., and Langenfelds, R.: Contribution of anthropogenic and natural sources to atmospheric methane variability, Nature, 443, 439–443, 2006. a

Braswell, B. H., Sacks, W. J., Linder, E., and Schimel, D. S.: Estimating diurnal to annual ecosystem parameters by synthesis of a carbon flux model with eddy covariance net ecosystem exchange observations, Glob. Change Biol., 11, 335–355, 2005. a

Bridgham, S. D., Pastor, J., Dewey, B., Weltzin, J. F., and Updegraff, K.: Rapid carbon response of peatlands to climate change, Ecology, 89, 3041–3048, 2008. a

Broecker, W. S. and Peng, T.-H.: Gas exchange rates between air and sea, Tellus, 26, 21–35, 1974. a

Carrassi, A., Bocquet, M., Bertino, L., and Evensen, G.: Data assimilation in the geosciences: An overview of methods, issues, and perspectives, Wiley Interdisciplinary Reviews: Climate Change, 9, e535,, 2018. a

Christensen, J. H., Hewitson, B., Busuioc, A., Chen, A., Gao, X., Held, I., Jones, R., Kolli, R. K., Kwon, W.-T., Laprise, R., and Magaña Rueda, V.: Regional climate projections, in: Climate Change 2007: the physical science basis Chapter 11, Cambridge University Press, 847–940, 2007. a

Ciais, P., Sabine, C., Bala, G., Bopp, L., Brovkin, V., Canadell, J., Chhabra, A., DeFries, R., Galloway, J., Heimann, M., and Jones, C.: Carbon and other biogeochemical cycles. Climate change 2013: the physical science basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change, Comput. Geom., 18, 95–123, 2013. a

Cronk, J. K. and Fennessy, M. S.: Wetland plants: biology and ecology, CRC press,, 2016. a

Dee, D. P.: Bias and data assimilation, Quarterly Journal of the Royal Meteorological Society: A journal of the atmospheric sciences, Applied Meteorology and Physical Oceanography, 131, 3323–3343, 2005. a

Dlugokencky, E.: NOAA/GML (gml. noaa. gov/ccgg/trends_ch4/), NOAA/GML,, 2021. a

Forster, P., Storelvmo, T., Armour, K., Collins, W., Dufresne, J.-L., Frame, D., Lunt, D., Mauritsen, T., Palmer, M., Watanabe, M., Wild, M., and Zhang, H.: The Earth’s Energy Budget, Climate Feedbacks, and Climate Sensitivity, Cambridge University Press, Cambridge, United Kingdom and New York, NY, USA, 923–1054,, 2021. a

Gelman, A., Roberts, G., and Gilks, W.: Efficient Metropolis jumping rules, in: Bayesian Statistics 5, edited by: Bernardo, J. M., Berger, J. O., Dawid, A. P., and Smith, A. F. M., 599–607, Oxford University Press,, 1996. a

Ghil, M. and Malanotte-Rizzoli, P.: Data assimilation in meteorology and oceanography, in: Advances in geophysics, 33, 141–266, Elsevier,, 1991. a

Gustafson, A.: On the role of terrestrial ecosystems in a changing Arctic, Media-Tryck, Lund University, Sweden, 2022. a

Hastings, W.: Monte Carlo sampling methods using Markov chains and their applications, Biometrika, 57, 97–109, 1970a. a, b

Hastings, W. K.: Monte Carlo sampling methods using Markov chains and their applications, Biometrika, 57, 97–97, 1970b. a

Hoffman, F. O. and Miller, C. W.: Uncertainties in environmental radiological assessment models and their implications, Tech. rep., Oak Ridge National Lab.,, 1983. a

Jacob, D., Bärring, L., Christensen, O. B., Christensen, J. H., de Castro, M., Deque, M., Giorgi, F., Hagemann, S., Hirschi, M., Jones, R., and Kjellström, E.: An inter-comparison of regional climate models for Europe: model performance in present-day climate, Clim. Change, 81, 31–52, 2007. a

Jennison, C.: Discussion on the meeting on the Gibbs sampler and other Markov chain Monte Carlo methods, J. Roy. Stat. Soc. B, 55, 54–56, 1993. a, b

Johansson, T., Malmer, N., Crill, P. M., Friborg, T., Åkerman, J. H., Mastepanov, M., and Christensen, T. R.: Decadal vegetation changes in a northern peatland, greenhouse gas fluxes and net radiative forcing, Glob. Change biol., 12, 2352–2369, 2006. a

Kallingal, J. T., Scholze, M., A. Miller, P., Lindstrom, J., Rinne, J., and Raivonen, M.: Optimising the CH4 simulations from the LPJ-GUESS model using Adaptive MCMC, Zenodo [code, data set],, 2022. a

Kirschke, S., Bousquet, P., Ciais, P., Saunois, M., Canadell, J. G., Dlugokencky, E. J., Bergamaschi, P., Bergmann, D., Blake, D. R., Bruhwiler, L., and Cameron-Smith, P.: Three decades of global methane sources and sinks, Nat. Geosci., 6, 813–823, 2013. a, b

Korrensalo, A., Männistö, E., Alekseychik, P., Mammarella, I., Rinne, J., Vesala, T., and Tuittila, E.-S.: Small spatial variability in methane emission measured from a wet patterned boreal bog, Biogeosciences, 15, 1749–1761,, 2018. a

Korrensalo, A., Mammarella, I., Alekseychik, P., Vesala, T., and Tuittila, E.: Plant mediated methane efflux from a boreal peatland complex, Plant Soil, 471, 375–392, 2022. a

Kuppel, S., Peylin, P., Chevallier, F., Bacour, C., Maignan, F., and Richardson, A. D.: Constraining a global ecosystem model with multi-site eddy-covariance data, Biogeosciences, 9, 3757–3776,, 2012. a

Lerman, A.: Geochemical processes. Water and sediment environments, John Wiley and Sons, Inc., 1979. a

Mathijssen, P. J., Väliranta, M., Korrensalo, A., Alekseychik, P., Vesala, T., Rinne, J., and Tuittila, E.-S.: Reconstruction of Holocene carbon dynamics in a large boreal peatland complex, southern Finland, Quaternary Sci. Rev., 142, 1–15, 2016. a

McGuire, A., Sitch, S., Clein, J. S., Dargaville, R., Esser, G., Foley, J., Heimann, M., Joos, F., Kaplan, J., Kicklighter, D., and Meier, R.: Carbon balance of the terrestrial biosphere in the twentieth century: Analyses of CO2, climate and land use effects with four process-based ecosystem models, Global Biogeochem. Cy., 15, 183–206, 2001. a

McGuire, A. D., Christensen, T. R., Hayes, D., Heroult, A., Euskirchen, E., Kimball, J. S., Koven, C., Lafleur, P., Miller, P. A., Oechel, W., Peylin, P., Williams, M., and Yi, Y.: An assessment of the carbon balance of Arctic tundra: comparisons among observations, process models, and atmospheric inversions, Biogeosciences, 9, 3185–3204,, 2012. a

Metropolis, N., Rosenbluth, A., Rosenbluth, M., Teller, A., and Teller, E.: Equations of state calculations by fast computing machines, J. Chem. Phys., 21, 1087–1092, 1953a. a

Metropolis, N., Rosenbluth, A. W., Rosenbluth, M. N., Teller, A. H., and Teller, E.: Equation of state calculations by fast computing machines, J. Chem. Phys., 21, 1087–1092, 1953b. a

Millington, R. and Quirk, J.: Permeability of porous solids, T. Faraday Soc., 57, 1200–1207, 1961. a

Nord, J., Anthoni, P., Gregor, K., Gustafson, A., Hantson, S., Lindeskog, M., Meyer, B., Miller, P., Nieradzik, L., Olin, S., Papastefanou, P., Smith, B., Tang, J., Wårlind, D., and and past LPJ-GUESS contributors: LPJ-GUESS Release v4.1.1 model code (4.1.1), Zenodo [code],, 2021. a

Rinne, J., Riutta, T., Pihlatie, M., Aurela, M., Haapanala, S., Tuovinen, J.-P., Tuittila, E.-S., and Vesala, T.: Annual cycle of methane emission from a boreal fen measured by the eddy covariance technique, Tellus B, 59, 449–457, 2007. a, b

Rinne, J., Tuittila, E.-S., Peltola, O., Li, X., Raivonen, M., Alekseychik, P., Haapanala, S., Pihlatie, M., Aurela, M., Mammarella, I., and Vesala, T.: Temporal variation of ecosystem scale methane emission from a boreal fen in relation to temperature, water table position, and carbon dioxide fluxes, Global Biogeochem. Cy., 32, 1087–1106, 2018. a

Roberts, G., Gelman, A., and Gilks, W.: Weak convergence and optimal scaling of random walk Metropolis algorithms, Ann. Probab., 7, 110–120, 1997. a

Santaren, D., Peylin, P., Bacour, C., Ciais, P., and Longdoz, B.: Ecosystem model optimization using in situ flux observations: benefit of Monte Carlo versus variational schemes and analyses of the year-to-year model performances, Biogeosciences, 11, 7137–7158,, 2014. a

Saunois, M., Jackson, R., Bousquet, P., Poulter, B., and Canadell, J.: The growing role of methane in anthropogenic climate change, Environ. Res. Lett., 11, 120207,, 2016. a

Saunois, M., Stavert, A. R., Poulter, B., Bousquet, P., Canadell, J. G., Jackson, R. B., Raymond, P. A., Dlugokencky, E. J., Houweling, S., Patra, P. K., Ciais, P., Arora, V. K., Bastviken, D., Bergamaschi, P., Blake, D. R., Brailsford, G., Bruhwiler, L., Carlson, K. M., Carrol, M., Castaldi, S., Chandra, N., Crevoisier, C., Crill, P. M., Covey, K., Curry, C. L., Etiope, G., Frankenberg, C., Gedney, N., Hegglin, M. I., Höglund-Isaksson, L., Hugelius, G., Ishizawa, M., Ito, A., Janssens-Maenhout, G., Jensen, K. M., Joos, F., Kleinen, T., Krummel, P. B., Langenfelds, R. L., Laruelle, G. G., Liu, L., Machida, T., Maksyutov, S., McDonald, K. C., McNorton, J., Miller, P. A., Melton, J. R., Morino, I., Müller, J., Murguia-Flores, F., Naik, V., Niwa, Y., Noce, S., O'Doherty, S., Parker, R. J., Peng, C., Peng, S., Peters, G. P., Prigent, C., Prinn, R., Ramonet, M., Regnier, P., Riley, W. J., Rosentreter, J. A., Segers, A., Simpson, I. J., Shi, H., Smith, S. J., Steele, L. P., Thornton, B. F., Tian, H., Tohjima, Y., Tubiello, F. N., Tsuruta, A., Viovy, N., Voulgarakis, A., Weber, T. S., van Weele, M., van der Werf, G. R., Weiss, R. F., Worthy, D., Wunch, D., Yin, Y., Yoshida, Y., Zhang, W., Zhang, Z., Zhao, Y., Zheng, B., Zhu, Q., Zhu, Q., and Zhuang, Q.: The Global Methane Budget 2000–2017, Earth Syst. Sci. Data, 12, 1561–1623,, 2020. a

Segers, R.: Methane production and methane consumption: a review of processes underlying wetland methane fluxes, Biogeochemistry, 41, 23–51, 1998. a

Sitch, S., Smith, B., Prentice, I. C., Arneth, A., Bondeau, A., Cramer, W., Kaplan, J. O., Levis, S., Lucht, W., Sykes, M. T., and Thonicke, K.: Evaluation of ecosystem dynamics, plant geography and terrestrial carbon cycling in the LPJ dynamic global vegetation model, Glob. Change Biol., 9, 161–185, 2003. a

Smith, B.: LPJ-GUESS-an ecosystem modelling framework, Department of Physical Geography and Ecosystems Analysis, INES, Sölvegatan, 12, 22362, 2001. a

Smith, B., Wårlind, D., Arneth, A., Hickler, T., Leadley, P., Siltberg, J., and Zaehle, S.: Implications of incorporating N cycling and N limitations on primary production in an individual-based dynamic vegetation model, Biogeosciences, 11, 2027–2054,, 2014. a, b

Susiluoto, J., Raivonen, M., Backman, L., Laine, M., Makela, J., Peltola, O., Vesala, T., and Aalto, T.: Calibrating the sqHIMMELI v1.0 wetland methane emission model with hierarchical modeling and adaptive MCMC, Geosci. Model Dev., 11, 1199–1228,, 2018. a, b, c, d

Tang, J., Miller, P. A., Persson, A., Olefeldt, D., Pilesjö, P., Heliasz, M., Jackowicz-Korczynski, M., Yang, Z., Smith, B., Callaghan, T. V., and Christensen, T. R.: Carbon budget estimation of a subarctic catchment using a dynamic ecosystem model at high spatial resolution, Biogeosciences, 12, 2791–2808,, 2015. a

Tarantola, A.: Inversion of travel times and seismic waveforms, Seismic tomography, 135–157,, 1987. a, b, c

Wania, R., Ross, I., and Prentice, I. C.: Integrating peatlands and permafrost into a dynamic global vegetation model: 1. Evaluation and sensitivity of physical land surface processes, Global Biogeochem. Cy., 23,, 2009a. a, b

Wania, R., Ross, I., and Prentice, I. C.: Implementation and evaluation of a new methane model within a dynamic global vegetation model: LPJ-WHyMe v1.3.1, Geosci. Model Dev., 3, 565–584,, 2010. a, b, c, d, e, f, g, h, i, j, k, l, m

Wania, R., Melton, J. R., Hodson, E. L., Poulter, B., Ringeval, B., Spahni, R., Bohn, T., Avis, C. A., Chen, G., Eliseev, A. V., Hopcroft, P. O., Riley, W. J., Subin, Z. M., Tian, H., van Bodegom, P. M., Kleinen, T., Yu, Z. C., Singarayer, J. S., Zürcher, S., Lettenmaier, D. P., Beerling, D. J., Denisov, S. N., Prigent, C., Papa, F., and Kaplan, J. O.: Present state of global wetland extent and wetland methane modelling: methodology of a model inter-comparison project (WETCHIMP), Geosci. Model Dev., 6, 617–641,, 2013.  a

Wramneby, A., Smith, B., Zaehle, S., and Sykes, M. T.: Parameter uncertainties in the modelling of vegetation dynamics–effects on tree community structure and ecosystem functioning in European forest biomes, Ecol. Model., 216, 277–290, 2008. a

Zhang, W., Miller, P. A., Smith, B., Wania, R., Koenigk, T., and Döscher, R.: Tundra shrubification and tree-line advance amplify arctic climate warming: results from an individual-based dynamic vegetation model, Environ. Res. Lett., 8, 034023,, 2013. a

Zhang, Z., Zimmermann, N. E., Stenke, A., Li, X., Hodson, E. L., Zhu, G., Huang, C., and Poulter, B.: Emerging role of wetland methane emissions in driving 21st century climate change, P. Natl. Acad. Sci. USA, 114, 9647–9652, 2017. a

Short summary
By unlocking the mysteries of CH4 emissions from wetlands, our work improved the accuracy of the LPJ-GUESS vegetation model using Bayesian statistics. Via assimilation of long-term real data from a wetland, we significantly enhanced CH4 emission predictions. This advancement helps us better understand wetland contributions to atmospheric CH4, which are crucial for addressing climate change. Our method offers a promising tool for refining global climate models and guiding conservation efforts