Assimilating multi-site eddy-covariance data to calibrate the wetland CH₄ emission module in a terrestrial ecosystem model

In this study, we use a data assimilation framework based on the adaptive Markov chain Monte Carlo (MCMC) algorithm to constrain process parameters in LPJ-GUESS model using CH₄ eddy-covariance flux observations from 14 different natural boreal, temperate, and arctic wetlands. The objective is to derive a single set of calibrated parameter values. The calibrated parameter values are then used in the model to validate its CH₄ flux output against independent CH₄ flux observations from five different types of natural wetlands situated in different locations, assessing their generality for simulating CH₄ fluxes from boreal, temperate, and arctic wetlands. The results show that the MCMC framework has substantially reduced the cost function (measuring the misfit between simulated and observed CH₄ fluxes) and facilitated detailed characterisation of the posterior parameter distribution. A reduction of around 50 % in RMSE was achieved, reflecting improved agreement with the observations. The results of the validation experiment indicate that for four out of the five validation sites the RMSE was successfully reduced, demonstrating the effectiveness of the framework for estimating CH₄ emissions from wetlands not included in the assimilation experiment. For wetlands above 45° N, the total mean annual CH₄ emission estimation using the optimised model resulted in 28.16 Tg C yr⁻¹ and for regions above 60 ° N it resulted in 7.46 Tg C yr⁻¹.

1 Introduction

Methane (CH₄) emissions from wetlands contribute 20 %–30 % to the total global emissions (IPCC AR6 Chap. 5: Canadell et al., 2021; Saunois et al., 2020). About one third to one half of these wetland emissions are from wetlands located at northern latitudes of North America, Europe, and Russia (Saunois et al., 2016a). According to the IPCC AR6 report, wetlands are the largest single source of uncertainty to the global CH₄ budget estimate. It is expected to have increased uncertainties in wetland CH₄ emissions in the future (Christensen et al., 2007), partly due to climate change and partly due to spatio-temporal changes in wetland extent (that in itself is partly a consequence of climate change) (Saunois et al., 2016b; Zhang et al., 2017). A key question to consider here is the extent to which these changes in emissions are occurring and how they will impact the future global greenhouse gas (GHG) budget and hence the climate. While current in situ measurement techniques such as eddy-covariance (EC) flux observations are promising for drawing assumptions at local scales, studies to date have faced difficulties in estimating wetland CH₄ emissions over large landscapes primarily due to the spatial heterogeneity of wetlands, strong temporal variability in emissions, and the limited coverage of long term flux observations (Saunois et al., 2020).

An attempt to overcome this limitation through process-based modelling of global CH₄ emissions was first initiated by Fung et al. (1991) followed by Christensen and Cox (1995) and more mechanistically by Cao et al. (1996) and Walter and Heimann (2000). These models were simple in structure and later more attention was given to improving process representation through the studies of mainly Segers and Leffelaar (2001), Gedney et al. (2004), and Zhuang et al. (2006). In the past decade, more advanced models have incorporated additional complexities such as dynamic water table fluctuations, temperature-dependent microbial activity, substrate diffusion, and oxidation in the soil column, along with broader applications across diverse wetland types and climatic conditions (e.g. Wania et al., 2010; Ringeval et al., 2010; Susiluoto et al., 2018). All of these past efforts indicate that comprehensive, process-based modelling of CH₄ emissions from wetland ecosystems is unquestionably a key way to understand the variability of wetlands and how they respond to environmental stresses and climate change (Saunois et al., 2020).

As all these models are approximations of the real world and exhibit their own uncertainties; here again the question is how to reduce the uncertainty for large scale applications. According to Kuppel et al. (2012), every terrestrial biosphere model contains uncertainties in five different ways: errors in real data used for calibration, errors in meteorological forcing, errors in process descriptions, errors in model parameter values, and inaccurate initial state of the model. The first two errors are related to measurement, while the last three are related to model formulation and are important to improve the model performance for general applications. There has been a growing effort to reduce uncertainty related to the last three sources of error factors in several ways. A popular method to reduce uncertainty in model parameters is to calibrate the model simulations against observations. Previous studies like Williams et al. (2009), Susiluoto et al. (2018) and Kuppel et al. (2012), based on different models, different data, and different parameter sets provide examples of improving model parameters and reducing uncertainties through data assimilation.

In this study, we consider uncertainties in parameter values of the CH₄ module of a global process-based ecosystem model, Lund–Potsdam–Jena General Ecosystem Simulator (LPJ-GUESS) v4.1 and aim to reduce their uncertainties. Dynamic global vegetation models (DGVMs) like LPJ-GUESS are state-of-the-art tools for studying the functioning of high latitude wetlands and estimating the dynamics in their global carbon balance (Sitch et al., 2003). In a previous study, Kallingal et al. (2024) have used EC flux observations collected at an individual site to investigate the potential of a Markov chain Monte Carlo (MCMC) type (Global Rao–Blackwellised Adaptive Metropolis, GRaB-AM) algorithm to optimise the parameters in the CH₄ module of LPJ-GUESS. The study showed that EC flux measurements of CH₄ contains useful information for optimising the CH₄ model parameters due to the high temporal resolution of the CH₄ flux measurements. However, the small spatial scale (site scale) and limited temporal extent of data collected from a single site could have over fitted parameters to the specificities of the particular site used. This points to the need for a more general approach. For example, in a study conducted by Groenendijk et al. (2011) the parameters of a photosynthesis model are optimised using EC flux observations from several Fluxnet sites. Similarly, studies like Kuppel et al. (2012) and Raoult et al. (2016) have constrained the parameters of a global ecosystem model using multi-site EC flux observations. Considering these studies and the results of Kallingal et al. (2024), we hypothesise that assimilating multi-site daily CH₄ flux observations using the GRaB-AM framework can derive a set of optimised general parameters capable of representing various types of northern wetlands.

The present study's objective is to investigate the capacity of the GRaB-AM framework developed by Kallingal et al. (2024) for calibrating selected model process parameters in a more general multi-site setup. We aim to use the calibrated model parameters to assess the extent to which this optimisation improves the model's ability to estimate CH₄ emissions from various wetlands across northern latitudes and to estimate the total mean annual budget at larger scales. We also aim to estimate the posterior process and parameter uncertainties, as well as the posterior parameter correlations. The 14 different arctic, boreal and temperate wetland sites chosen for this study were selected to encompass diverse bioclimatic and geographical characteristics of wetlands. This deliberate selection aimed to endow the optimised parameters with the ability to represent a range of wetland types, independent of their distinct climatic and geographical features. We then perform additional validation simulations to evaluate the performance of the calibrated model against five different, independent validation sites to verify our above-mentioned hypothesis. In the following, we first give a brief overview of the LPJ-GUESS model, the data used in the assimilation and the data assimilation methodology itself. The assimilation results are then presented and discussed in Sects. 3 and 4 before we end with the conclusion (Sect. 5).

2 Data and methodology

2.1 LPJ-GUESS model

LPJ-GUESS represents the structure and dynamics of terrestrial ecosystems from local to global scales (Smith, 2001; Smith et al., 2014). The model combines basic eco-physiological features with detailed vegetation dynamics and canopy structure as used in forest gap models and includes an interactive nitrogen cycle (Smith et al., 2014). In version 4.1, which we used for this study, global vegetation is grouped into 13 different co-occurring mixtures of plant functional types (PFTs) and 5 additional PFTs that can only exist on peatland stands. The model input data consists of climate parameters (mean daily air temperature, precipitation, and incoming shortwave radiation), atmospheric CO₂ concentrations, and soil properties. LPJ-GUESS simulates vegetation dynamics, ecosystem biogeochemistry, water cycling, and energy and carbon fluxes on a daily time step. The peatland module in LPJ-GUESS contains detailed representations of wetland PFT characteristics and bio-geo-chemical processes such as estimation of peat temperature, hydrology, and ecosystem exchanges, including CH₄ emissions.

2.1.1 Main process description in CH₄ module of LPJ-GUESS

A detailed description of the wetland and CH₄ emissions module is given in Wania et al. (2010) and in Kallingal et al. (2024). Here, we only briefly summarise the most important aspects of the module. The wetland peat in LPJ-GUESS is 1.5 m deep and is divided into an acrotelm with a thickness of 0.3 m with varying water table depth (wtd) and a permanently saturated catotelm. Peat hydrology and peat temperature in this layered structure depend on the composition of each layer and prevailing meteorological conditions. The five types of PFTs implemented in the wetlands are Sphagnum mosses, C₃ graminoids, evergreen and deciduous shrubs, and a generic herbaceous cushion lichen moss. Shade mortality, inundation stress, and daily desiccation stress are limiting factors for the existence and productivity of these PFTs. The basic concept of the CH₄ module in LPJ-GUESS is a soil carbon pool distributed in proportion to the root distribution. This “potential carbon pool” serves as the substrate for methanogens to produce CH₄. A portion of the soil carbon get transformed to CH₄ and/or CO₂ depending on the hydrological conditions. A fraction of the produced CH₄ is in dissolved form and the remainder is in gaseous form. A part of this CH₄ is oxidised by the oxygen in the soil and the other part is eventually transported to the atmosphere through either diffusion, plant-mediated transport, or ebullition. As shown in Eq. (1) the sum of the emissions through these three pathways constitutes the total CH₄ flux from the soil to the atmosphere (Wania et al., 2010; Kallingal et al., 2024):

$$\begin{array}{}\text{(1)}& F_{{\mathrm{CH}}_{\mathrm{4}}}={{\mathrm{CH}}_{\mathrm{4}}}_{\text{diff}}+{{\mathrm{CH}}_{\mathrm{4}}}_{\text{plant}}+{{\mathrm{CH}}_{\mathrm{4}}}_{\text{ebul}},\end{array}$$

where $F_{{\mathrm{CH}}_{\mathrm{4}}}$ is the total CH₄ flux, ${{\mathrm{CH}}_{\mathrm{4}}}_{\text{diff}}$ is the CH₄ flux component from diffusion, ${{\mathrm{CH}}_{\mathrm{4}}}_{\text{plant}}$ is the CH₄ flux component from plant-mediated transport, and ${{\mathrm{CH}}_{\mathrm{4}}}_{\text{ebul}}$ is the CH₄ flux component from ebullition.

2.1.2 Parameter selection

The parameters selected for optimisation in this study are shown in Table 1. For this study we have considered 10 out of 11 parameters calibrated by Kallingal et al. (2024) in their single-site optimisation. The parameter w_tiller, which is the plant tiller weight, is removed because it showed high correlation with r_tiller in Kallingal et al. (2024).

Table 1 Parameters selected for the multi-site assimilation. Model prior values, prior standard deviation (std), units, and description of the parameters are given.

2.2 Flux sites and climate data

As mentioned in Sect. 1, for this study we selected 14 natural wetland sites for the assimilation and five additional wetland sites for validation (see Fig. 1, Table 2 and Table 4). The selection criteria for the sites were (1) that they are located above 40° N, (2) that they include at least 3 years of consecutive CH₄ (at least the summer) measurements and meteorological measurements available for the sites used for assimilation (same criteria but at least 2 years of measurements available for the sites used for validation), and (3) that they represent arctic, boreal, or temperate ecosystems. We did not consider lakes, uplands, etc. as they are beyond the scope of the present study. The 19 sites are representative of a range of wetland types including fens, mires, bogs, marsh, and a wet tundra (for validation).

Figure 1 Location of measurement sites selected for the study. The 14 sites used for assimilation are indicated as circles and the 5 sites used for validation are indicated as stars. The numbers in the left and right side legends correspond to the numbers assigned to the sites in Tables 2 and 4, respectively.

In total we have used 18 437 data points of daily CH₄ EC flux observations for assimilation, spanning approximately 93 measurement years. We did not consider any sites where the climate data had gaps of more than 14 d. In case of gaps smaller than 14 d the data are backfilled (by copying the values from the preceding cells) to ensure dataset continuity, as data with gaps cannot be used in the model as input. For validation, we utilised 5111 data points spanning a total of 14 measurement years. For all sites, only the available CH₄ flux observations are compared with the model both for assimilation and validation; i.e. we do not use any gap-filled CH₄ fluxes. It should be noted that for the assimilation sites Att, Bon, Sco, Uoa, and Zak, most of the data from winter months were missing, while for the validation sites Lgt and Wpt, most of the data from summer months were missing. See Appendix A for the details regarding data collection.

Łakomiec et al. (2021)Todd and Humphreys (2018)Ueyama et al. (2020)Euskirchen and Edgar (2020)Granberg et al. (2001)Koebsch and Jurasinski (2020)Lohila et al. (2020)Aurela et al. (2015)Desai and Thom (2020)Bohrer and Morin (2020)Sonnentag and Helbig (2020)Rinne et al. (2018)Iwata et al. (2020)Scheller et al. (2021)Sachs et al. (2020)

Table 2 Site information and data references of the 14 natural wetland sites used for assimilation. MAT refers to the mean annual temperature and MAPr to the mean annual precipitation. For Bib, Deg, Lom, Los, Ole, Sco, Sii, and Uao, MAT and MAPr are extracted from the corresponding grid cells of the site locations of the WorldClim 2.0 gridded product (Fick and Hijmans, 2017) and for the rest of the sites the information is taken from their references. The table also includes type, coordinates, and climate zones of the wetlands, as well as the time period of data availability.

2.3 Data assimilation system

To find an optimal posterior parameter set we used an adaptive Rao–Blackwellised Markov chain Monte Carlo Metropolis–Hastings (MCMC-MH) algorithm (Andrieu and Thoms, 2008) to iteratively reduce the model–observation misfit in terms of a so-called cost function (see Eq. 2) that compares the modelled observable with the field measurements. As mentioned in Sect. 1, details of this search algorithm and its application to a single-site optimisation have been described in Kallingal et al. (2024). Efficient sampling from the target distribution requires a proposal distribution that correctly represents the dependence structure of the target; to avoid manual tuning of the proposal they used an adaptive MCMC to tune the proposal distribution, where the Rao–Blackwellisation improves the adaptation step. The tuning improves the MCMC convergence speed and avoids cases of incomplete convergence (Andrieu and Thoms, 2008), especially for a complex non-linear model like LPJ-GUESS.

We assumed errors in observation and parameters in the form of Gaussian distributions yielding the cost function J(x),

$$\begin{array}{}\text{(2)}& \begin{array}{rl}J\left(x\right)& ={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\sum _{i=\mathrm{1}}^n(Y_i-M_i(x\left){)}^t{R}_i^{-\mathrm{1}}\right(Y_i-M_i\left(x\right))\\ & +{\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}(x-x_{\mathrm{p}}{)}^t{\mathbf{B}}^{-\mathrm{1}}(x-x_{\mathrm{p}}),\end{array}\end{array}$$

where Y_i, M_i(x), and R_i are the CH₄ flux observations, model simulations, and the covariance matrix of the observation errors, respectively at the ith site, x_p are the expected prior parameters, and B is the prior parameter error covariance matrix. Thus, the first term represents the model–data misfit weighted by the observation error covariances and the second term represents the prior information on the parameters weighted by the parameter error covariances.

Samples are generated by drawing x_prop from a proposal distribution and then either accepting the proposed state (x_i=x_prop) or keeping the current state ($x_i=x_{i-\mathrm{1}}$) based on the posterior probabilities. The probability of accepting the proposed state (α) is generally computed as

$$\begin{array}{}\text{(3)}& \mathit{\alpha }=min\left(\mathrm{1},{\displaystyle \frac{P\left(x_{\text{prop}}\right)}{P\left(x_{i-\mathrm{1}}\right)}}\right).\end{array}$$

Here, P(x_prop) is the posterior probability of the proposed state and P(x_i−1) is the posterior probability of the current state, both computed using the cost function, Eq. (2). The acceptance probability ensures a balanced exploration of the parameter space, accepting states that improve the fit while allowing occasional exploration of less favourable regions (see Andrieu and Thoms, 2008, for technical details and Kallingal et al., 2024, for implementation).

Table 3 Data availability and threshold estimated for the base error values. The number of available flux observations from each site is also provided.

For optimisation we used the GRaB-AM with a chain length of 100 000 iterations, where each iteration involves one complete model run for all 14 sites. As mentioned above, daily averages of flux observations collected from the above-mentioned 14 sites are assimilated simultaneously. For each site, only the actual CH₄ flux observations; i.e. non-gap-filled data are used to calculate the cost function. Having multiple sites in this framework, one crucial challenge was scaling the cost function to maintain a balanced representation of each site's contribution to the overall model–data misfit. This process is particularly relevant when sites exhibit variations in the magnitude of their individual cost functions or when the number of observations at each site differs significantly. Here, the scaling factors are carefully chosen to ensure an approximately equal representation of all sites in the cost function, regardless of their individual characteristics, and to ensure that each site has an equal influence on the optimisation outcome.

Considering the difficulty of calculating error correlations in the flux observations, we only considered errors in individual observations; i.e. we did not consider off-diagonal elements in specifying the observational error covariance matrix R in the cost function (Eq. 2). Estimating the exact observation error for each site is again challenging. Assigning a constant percentage error for all measured values could introduce a bias, as it would result in high error values for measurements with high magnitudes and very low error values for observations with small magnitudes. To overcome this challenge, we followed the procedure introduced by Knorr and Kattge (2005) for the case of assimilating CO₂ eddy-covariance observations and assign a threshold value set at 5 % of the variance of the distributions of observations, calculated separately for each site. Values below this threshold are identified and a uniform error is assigned to them (see Table 3). An error of 30 % is estimated for the observations greater than the threshold values. For the matrix B in the Eq. (2), for each parameter a standard deviation of 40 % of their possible range is assumed based on the expertise of LPJ-GUESS modellers.

2.4 Posterior estimation

After completing a full run of the MCMC chain, the posterior parameter error covariance matrix (B_p) is estimated from the prior error covariance matrices of the observations, R, and of the parameters, B, and the linearisation of the model at the minimum of the cost function, J(M_∞), as described in Tarantola (1987):

$$\begin{array}{}\text{(4)}& {\mathbf{B}}_{\mathrm{p}}=[M_{\mathrm{\infty }}^t{R}^{-\mathrm{1}}{M}_{\mathrm{\infty }}+B^{-\mathrm{1}}{]}^{-\mathrm{1}}.\end{array}$$

B_p is then used to estimate the level of optimisation of each parameter and the sensitivity of the cost function to them. The posterior parameter uncertainties have been estimated from the square root of the diagonal elements of B_p. Large absolute values of posterior error correlations indicate that the observations do not provide independent information to distinguish the effects of a given parameter pair (Tarantola, 1987).

From the 100 000 samples yielded by the GRaB-AM framework 75 % of this chain was discarded as burn-in. The remaining part of the chain, which we consider converged to its stationary distribution, was used for calculating the maximum a posteriori estimation (MAP) and posterior mean estimations. The posterior distributions of parameters are classified as “well-constrained”, “poorly constrained”, and “edge-hitting” parameters. The well-constrained parameters are characterised by a clearly defined unimodal distribution with a low standard deviation. Conversely, poorly constrained parameters exhibit a relatively flat multimodal distribution with a large standard deviation. For a more precise estimation, we classified posterior parameter distributions as poorly constrained if the standard deviation exceeded 20 % of the total range. Edge-hitting parameters cluster near one of the edges of their prior range, as described by Kallingal et al. (2024) and Braswell et al. (2005).

2.5 Experimental setup

The model was spun up for 500 years using Climate Research Unit (CRU) meteorological forcing data (Harris et al., 2014), which includes daily measurements of air temperature, precipitation, and incoming shortwave radiation, to bring the model's state variables, i.e. the various carbon pools, to initial equilibrium. After spinning up, the model was run for the 14 assimilation sites using local, daily meteorological forcing data collected directly at the sites. We have bias-corrected the CRU data for the grid cells in which the sites are located to enforce agreement with monthly mean values of the site-specific meteorological data. For this we have used at least 2 years of meteorological data that are recorded prior to the time period of the site flux observations that we used for the assimilation. Atmospheric CO₂ concentration, as described in McGuire et al. (2001) and updated until recent years using data from the NOAA Global Monitoring Laboratory (NOAA-GML, https://gml.noaa.gov/ccgg/trends, last access: 29 April 2025) is used as the CO₂ concentration input to the model.

Most of the CH₄ flux observations at the sites were available with a half-hourly resolution, but for the assimilation we used daily mean values corresponding to the LPJ-GUESS temporal resolution. Also, using daily values reduces the complexity of error correlations of half-hourly data and is better suited for a broad range timescale assimilation (Lasslop et al., 2008). Days with less than 50 % of half-hourly CH₄ data availability were removed from the assimilation. To test whether the resulting posterior parameters can enhance model performance for other individual wetlands within the study area, we designed a validation experiment that uses additional flux observations not included in the assimilation experiment. Independent flux observations from five different wetlands located in various parts of the study area were used for the experiment (see Sect. 2.2 and Table 4).

Merbold et al. (2020)Jacotot et al. (2020)White et al. (2023)Schmid and Klatt (2020)Chen and Chu (2020)

Table 4 Site information and data references of five natural wetland sites used for validation. MAT refers to the mean annual temperature and MAPr to the mean annual precipitation collected from their references. The table also includes the time period of data collected, the availability of data and the type and climate zones of the wetlands.

Additionally, we conducted a mean annual budget estimation for all wetlands above 45° N and for all wetlands above 60° N using the optimised parameters. The model was run using the peat fraction map PEATMAP (Xu et al., 2018) and daily gridded climate data from the Climatic Research Unit-Japanese Reanalysis (CRU-JRA) (Harris et al., 2020) as model inputs. CRU-JRA data can be accessed at https://rda.ucar.edu/datasets/ds628.0/ (last access: 29 April 2025). The results were then compared with the output of the JSBACH-HIMMELI model, as described in Zhang et al. (2023) (hereafter referred to as JSBACH-H), for regions above 45° N (for easier comparison with other literature) as well as with the estimates from the Global CH₄ Budget of the Global Carbon Project (GCP) (Saunois et al., 2025) for regions above 60° N.

3 Results

3.1 Posterior parameter distributions

The PDFs of the parameters from the MCMC chains after the burn-in are displayed in Fig. 2. All parameters, except for ϕ_tiller, f_air, and por_acro, are well-constrained with only one peak in the PDF. However, the parameters ϕ_tiller, f_air, and por_acro are rather poorly constrained, exhibiting some clustering around multiple peaks in the PDF and having large standard deviations. The R_moist, f_air, and λ_root showed high skewness and kurtosis. However, the smaller kurtosis values of f_oxid, por_cato, and $R_{{\mathrm{moist}}_{\mathrm{an}}}$, along with their low skewness, indicate that they closely resemble Gaussian distributions. The remaining four parameters showed low skewness, suggesting agreement with a Gaussian distribution, but with very low kurtosis indicating somewhat flatter distributions than a Gaussian distribution. The figure shows that, except for λ_root, none of the parameters exhibited edge-hitting behaviour, indicating that the hypothetical boundaries assigned for each parameter align well with the model structure. The parameter λ_root finds its solution nearly at the lower edge. Most parameters displayed posterior solutions far from their prior values, i.e. the prior values where outside the posterior mean estimate ±1σ, except for ϕ_tiller. For ϕ_tiller, both the MAP and posterior mean appeared close to the prior value, i.e. within 1σ of the posterior mean estimate.

Figure 2 Probability distribution functions (PDFs) of the posterior obtained after the GRaB-AM experiment. The blue curves are smoothed Gaussian kernel estimates on the posterior histograms, while the black curves represent the prior distributions. The dotted vertical lines in green, orange, and black correspond to the posterior mean, MAP, and prior means, respectively. The shaded grey area in the distributions represents the 1σ error estimate of the PDFs. Skewness and kurtosis values for each posterior distribution are provided in insets.

3.2 Posterior parameter estimates and cross-correlation

The posterior parameter values estimated from the MAP and the posterior mean estimate, along with their standard deviations, are presented in Table C1 and Fig. 2. The MAP and posterior mean estimates of all parameters remained close to each other, except for ϕ_tiller, r_tiller, and f_air. The parameters λ_root and ϕ_tiller exhibited high standard deviations in their posterior distributions, indicating greater uncertainty. In contrast, the parameters ${\mathrm{CH}}_{\mathrm{4}}/{\mathrm{CO}}_{\mathrm{2}}$, r_tiller, por_cato, and r_most showed low standard deviations, suggesting more precise posterior estimates.

A cross-correlation plot (see Fig. 3) shows the correlation between all 10 parameter pairs after optimisation to examine potential optimisation issues due to parameter correlation. High positive or negative correlations suggest that these parameters may convey similar information and one of them might be redundant in further studies. The results show that not many parameters have strong positive or negative correlations, except for the correlation between ${\mathrm{CH}}_{\mathrm{4}}/{\mathrm{CO}}_{\mathrm{2}}$ and $R_{{\mathrm{moist}}_{\text{an}}}$, which has a high negative correlation of −0.82. Many slight positive correlations are observed, with a few pairs like por_acro and por_cato, por_cato and f_air having comparatively higher values. A detailed discussion of the parameter correlation and their possible impacts on the posterior estimates is given in Sect. 4.1.

Figure 3 Posterior correlations between parameters derived from the GRaB-AM optimisation. In the upper triangle of the figure, negative correlations are depicted in blue and positive correlations are shown in red. The numerical labels on the upper triangle correspond to values of Pearson's correlation coefficient. The diagonal panels exhibit 1-D histograms for each model parameter. The lower triangle displays two-dimensional marginal distributions for each parameter. The grey dots on the marginal distributions represent the parameter values obtained from the posterior GRaB-AM chain. The ranges of the distributions are labelled on the left and bottom of the figure.

3.3 Performance of the optimisation

The site-wise data–model misfit, presented in terms of RMSE both before and after optimisation, along with the average RMSE for all sites combined is presented in Fig. 4. Most sites demonstrate a substantial reduction in RMSE, with many achieving over 50 % improvement. The unweighted cost function values for each of the 14 sites individually and their collective sum, along with the corresponding RMSE and reduced χ² values, are presented in Table C2 in the appendix. The estimated prior cost function value was 1 763 294.9. After optimisation the value changed to 79 296.4 (around 5 % of the prior) for the posterior mean estimate. The total χ² value is 8.6, accounting for measurement and parameter uncertainty. Notably, none of the sites exhibit overconfidence, with individual χ² values larger than 1, except for the site Hue, which shows a χ² value slightly below but close to 1. The median RMSE reduction is 52.2 % and the χ² value is 6.8. Taking these values as thresholds, the sites Abi, Att, Sii, and Lom show significant reductions in RMSE but with high χ² values. In contrast, sites such as Bib, Deg, Hue, Sco, Zak, and Zar exhibit low RMSE reductions but low χ² values. The sites Bon, Los, Ole, and Uoa display both high RMSE reductions and low χ² values. None of the sites are characterised by low RMSE reduction and a high χ² value.

Figure 4 Prior and posterior Root Mean Square Error (RMSE) estimates are provided for each of the 14 sites individually, along with the combined average values. In the figure, purple, cyan, and grey bars represent the RMSE corresponding to the prior estimate, maximum a posteriori (MAP) estimate and posterior mean estimate, respectively.

3.4 Impact of the optimisation

The observed CH₄ flux data available from all 14 sites collectively provide a total of 888.8 g C m⁻². The prior sum estimate for the corresponding data was 1957 g C m⁻². This value reduced to 850.9 g C m⁻² after optimisation. A reduction of approximately 56.52 % from the prior to the posterior CH₄ flux was observed. Compared to the observation, the posterior resulted only a slight underestimation of −38.1 g C m⁻². Time series of the annual sums of CH₄ emissions at 4 of the 14 sites used in this study are shown in Fig. 6a. The time series of the remaining sites are shown in Appendix C3, Fig. C3 for convenience. All four sites shown in Fig. 6a exhibited a better fit to the annual budget after the optimisation. Particularly noteworthy is the site Bib, which had a prior estimation in 2016 significantly deviating from the observed values. This discrepancy was corrected in the posterior estimate and the annual posterior CH₄ emissions at all four sites align well with the flux observations after the optimisation. Figure C3 shows that the sites Att, Los, Uoa, and Zar have shown improvement in all the years assimilated. Deg improved in more than half of the years assimilated. Hue improved in 2013 but did not improve in the other years. Bon improved in 2 out of 3 years but was slightly overestimated in 2017. For Zak, although the prior and posterior merged in some years, it showed improvement in almost all years. Lom showed improvement in 2010, 2012, 2013, and 2014. Sco showed improvement in 2016 and 2017, which accounts for half of the total assimilated years. Figure 6b displays the mean annual sums of CH₄ estimated at all 14 sites. The figure illustrates that the highest contribution of CH₄ came from the sites Hue and Zar, which have the highest mean annual temperature (MAT) as compared to other sites. The lowest contributions are from Zak and Abi, which have below-zero MAT. Abi, Att, Los, Ole, Sco, Sii, Uoa, Zak, and Zar showed improvement in the mean annual budgets after the optimisation. The remaining sites did not show an improvement.

Table 5 presents the total uncertainty for each site and the total uncertainty estimated for all sites together (see Sect.  B in the appendix for technical details regarding the uncertainty estimation). The total uncertainty of the posterior CH₄ flux estimates for all the sites together was 0.19 g C m⁻² d⁻¹, whereas for the prior fluxes it was 0.36 g C m⁻² d⁻¹. This results in a reduction of the total uncertainty of around 50 % after the optimisation. Comparing the prior and posterior RMSE (Fig. 4) and the uncertainty reduction, it can be concluded that the more constrained sites, such as Abi, Att, Bon, and Uoa, exhibited high uncertainty reduction. Notably, Abi and Att, which had the highest prior RMSE, showed a reduction of uncertainty of around 95 % and 82 %, respectively. A low reduction in uncertainty was mainly observed in sites that demonstrated a low reduction in RMSE. In contradiction to this, even though the RMSE reduction observed in the case of Zak is very small, this site showed an uncertainty reduction of around 33 %.

Table 5 Prior and posterior total uncertainty estimates (σ (g C m⁻² d⁻¹)). Total parameter and model uncertainty estimates separately for each site and collectively for all sites combined are shown. Uncertainty is abbreviated to “unc.”.

Figure 5 illustrates the time series of observed, prior, and posterior fluxes for all 14 sites considered in the assimilation. The posterior model continued to overestimate emissions at the Abi, Att, Los, Ole, Sii, and Uoa sites. For Bib (except for 2015), Deg, Hue, and Zak emissions were consistently underestimated. The remaining four sites demonstrated reasonably good agreement of the posterior estimates with the flux observations. The model adeptly captures the seasonal cycles of CH₄ emissions from all wetlands, both for flux estimates using the prior and the posterior parameter values. Generally, no significant phase shift is observed in either the prior or the posterior estimates. However, for Bib, both prior and posterior estimates exhibit a slight phase shift to early summer, while the sites Hue, Los and Zar show a similar shift to late summer.

Figure 5 The CH₄ simulation from the LPJ-GUESS model from 14 different wetland sites (purple dots) after optimising with the GRaB-AM algorithm. The black dots are the real CH₄ flux observations from corresponding wetlands. The teal dots are the prior simulation with the prior model parameters used to start the MCMC chain; 3 d running averages are calculated from the original time series and from most of the figure a few outliers on the vertical axis have been removed for better visualisation.

Figure 6 Panel (a) displays the annual sum estimation of CH₄ from 4 out of 14 sites used in the study. The sites are represented in different colours with distinct markings to distinguish between observation, prior, and posterior. Panel (b) presents the mean annual sums of CH₄ estimated for all 14 sites used in this study as bar charts. The averages of the flux observations are calculated with only the available daily averages used for the assimilation.

Summer and winter anomalies

Figure 7 illustrates the mean annual summer and winter emissions for all sites and their corresponding standard deviations. Figure C2 in the appendix shows their “summer” and “winter” anomalies. It should be noted that the winter mean and anomaly estimation for the sites Att, Bon, Sco, Uoa, and Zak was conducted with only a very limited number of available data points, as most of them were missing. Conversely, it should be taken into account that proper winter measurements were not carried out at these sites, given the almost negligible emission estimates during the winter months due to their extremely cold temperatures. For all these sites, the mean annual temperature (MAT) was estimated to be below zero (refer to Table 2 and the corresponding site references).

Figure 7 Mean annual summer (a) and winter (b) emissions for all sites, along with their corresponding standard deviations. The dots represents mean seasonal values and error lines indicates their 1 standard deviation.

After the optimisation, both mean annual summer and winter emission estimates for sites Abi, Att, Los, Sco, Sii, and Zar exhibited improved agreement with the flux observations. The sites Uoa and Bib showed improvement in summer but not in winter, whereas the sites Bon, Hue, Lom, and Ole improved in winter but not in summer. The sites Zak and Deg did not show any significant improvement in either season. For the site Deg, the prior estimates were closer to the flux observations than the posterior and for the site Zak no significant changes were observed. Summer emissions in Zak were underestimated by both the prior and posterior models. This may be attributed to the relatively low (below zero) mean annual temperature (MAT) of −8.6 at this location, variation of summer months with latitude, and the seasonality of decomposition. In contrast, for the sites Abi and Att with negative MAT, the model tends to overestimate summer emissions but the observed average summer emissions were comparatively lower. However, despite a MAT of −2.9 °C at the site Sco and −2.6 °C at the site Uoa, both these sites demonstrated relatively high summer emissions, which were better captured by the model using the posterior parameter values. Sites with a higher MAT, such as Hue, Bib, and Zar, exhibited the highest summer emission values. Although, the site Ole, which has the highest MAT of 12.1 °C, displayed comparatively lower summer emissions. This difference could be influenced by the substantial MAPr of 1120 mm at the Ole site.

For most of the sites, the observed winter annual mean was very close to zero, except for Bib, Deg, Hue, Ole, and Zar. For Hue, Ole, and Zar, the posterior estimate showed better performance in capturing the seasonal trends in flux observations than the prior. The site Hue exhibited high winter emissions. When the air temperature input for this site was estimated, it showed a mean value of 3.8 °C in winter months. Overall, although the majority of sites showed improved estimation of winter and summer emissions, some of the sites remained unchanged or did not perform better than before. The estimation of the standard deviation for summer and winter months showed a reduction for all sites after optimisation.

3.5 Validation of optimised parameters and large scale simulation

Table 6 and Fig. 8 show the results of the validation experiment explained in Sect. 2.5. A collective reduction of 58.8 % in RMSE was observed across all the validation sites. The RMSE decreased from 0.4 in the prior simulation to 0.16 in the posterior simulation. The total sum of flux observations for all five sites was 319.4 g C m⁻². Corresponding posterior simulation resulted in a total sum of 421.3 g C m⁻², compared to 521.9 g C m⁻² in the prior simulation. This represents a 19.3 % decrease in the posterior estimate compared to the prior. Similarly, the prior σ was 0.61, which was reduced to 0.17, representing a 72.13 % reduction.

Table 6 RMSE reduction and the changes in the total emission estimate (in %) for the validation sites, along with their prior and posterior uncertainty (σ) estimates.

The posterior estimates for four of the five sites showed improved RMSE (Fig. 8). For the sites Sfn, Myk, and Che, reductions of 81.1 %, 59.8 %, and 31.2 % were observed, respectively. The RMSE improvement for the site Lgt was negligible, with a value of 0.6 %. The site Wpt, a temperate marsh (note that marshes were not included in the sites used for assimilation), exhibited a very slight increase in RMSE. The total prior model–data mismatch of CH₄ estimated at this site during the time period was 72.5 g C m⁻², which increased to 78.98 g C m⁻² after optimisation. Including Wpt, the posterior estimates of all these sites appeared improved in terms of σ reduction. Wet tundra was not used for assimilation; however, the site Che, a wet tundra used for validation, demonstrated a remarkable 31.2 % reduction in RMSE, along with a significant reduction in total posterior σ, from 0.29 to 0.19. For all the sites except Che, the prior was overestimating the posterior flux. For Che, there was an increase of 21.5 % in the posterior estimate.

Figure 8 Scatterplot comparing prior and posterior RMSE values for each site. Each point represents one of the study sites. The 1:1 dashed line indicates no change in RMSE. Points below the line indicate improved model performance after assimilation.

Figures 9 and 10 show the results of large scale simulations of CH₄ emissions from wetlands located above 45° N. Figure 9 includes prior and posterior simulations and their differences, as well as the fluxes from JSBACH-H simulations for the same domain. The maps in the top row indicate that both the prior and posterior simulations identify hotspots of emission in North America, Canada, Russia, and Europe, though the total posterior estimate is clearly smaller compared to the prior. However, examining the differences of prior and posterior emission illustrated in the bottom left map reveals that the smaller posterior emissions are not systematic over the whole simulation domain. Essentially the magnitude of the emissions at the hotspot locations is reduced in the posterior, whereas surrounding areas of low emissions in the prior are slightly enhanced in the posterior. The bottom right map displays the CH₄ simulation from JSBACH-H, with a resolution 1.875° × 1.875°, for comparison.

It is evident that both LPJ-GUESS and JSBACH-H exhibit similar patterns in high and low emissions. The total mean emission over 2010–2020 was 32.32 Tg C yr⁻¹ for the LPJ-GUESS prior simulation and 28.16 Tg C yr⁻¹ for the posterior simulation. For JSBACH-H, the mean emission over 2010–2016 was 15.66 Tg C yr⁻¹. The posterior simulation of LPJ-GUESS is closer to the result from JSBACH-H but still much larger. When comparing CH₄ emissions from wetlands above 60° N, the agreement between the posterior LPJ-GUESS estimate (7.46 Tg C yr⁻¹) and the JSBACH-H estimate (7.62 Tg C yr⁻¹, mean over 2010–2016) is much closer. However, the LPJ-GUESS prior estimate (8.69 Tg C yr⁻¹) aligns better with the Global CH₄ Budget estimate from the GCP (9.0 Tg C yr⁻¹, mean over 2010–2019).

Figure 9 Total mean annual CH₄ estimates above 45° N. (a) Prior CH₄ emissions from LPJ-GUESS averaged from 2010 to 2020 with a spatial resolution of 0.5° × 0.5°. (b) Posterior CH₄ emissions from LPJ-GUESS averaged from 2010 to 2020 with a spatial resolution of 0.5° × 0.5°. (c) Difference between the prior and posterior emissions from LPJ-GUESS. (d) CH₄ emissions from JSBACH-H averaged from 2010 to 2016 with a spatial resolution of 1.875° × 1.875°.

4 Discussion

In this study we have conducted detailed examination of the posterior results to assess parameter distributions, uncertainties, error correlations, changes in flux components, and model–observation fit as well as the changes in summer and winter anomalies. The well-constrained seven parameters discussed in Sect. 3.1 with single-peaked PDFs indicate strong convergence and effective estimation, suggesting the model is sensitive to the data for these parameters. However, the poorly constrained parameters with multiple peaks highlight areas where the model might need refinement or where additional types of observations could improve parameter estimation. The ranges of posterior parameter distributions indicate a successful search within the permitted parameter ranges (Fig. 2). The non-Gaussian behaviour of the parameters ϕ_tiller, f_air, and por_acro described in Sect. 3.1 can still be considered acceptable, as the posterior distribution is expected to be perfectly Gaussian only under the assumption of a linear model. Beyond the inherent nonlinearity of LPJ-GUESS, this behaviour also suggest that the assimilated fluxes provide more complex information about the wetland processes than initially assumed or that there are complex interactions between the model parameters. This complexity may warrant further investigation into the model structure. It also indicates that the model could benefit from accounting for errors that place less weight on extreme observations, allowing for a smoother response curve from the model.

4.1 Parameter correlations and their implications

R_moist and $R_{{\mathrm{moist}}_{\text{an}}}$: these parameters are related to the moisture response in the acrotelm and catotelm, respectively. A smaller posterior value of R_moist and $R_{{\mathrm{moist}}_{\mathrm{an}}}$ would result in a slower soil carbon turnover time in the posterior model, leading to slightly less carbon available for CH₄ production. Although this could decrease the total decomposed carbon in the soil, the strong negative correlation between $R_{{\mathrm{moist}}_{\mathrm{an}}}$ and ${\mathrm{CH}}_{\mathrm{4}}/{\mathrm{CO}}_{\mathrm{2}}$, and the weak negative correlation between R_moist and ${\mathrm{CH}}_{\mathrm{4}}/{\mathrm{CO}}_{\mathrm{2}}$, indicate that the decrease in these parameters has influenced the increase in the CH₄ fraction from this reduced amount of decomposed soil fraction. This indicates that other factors such as water table depth, availability of oxygen, soil temperature, etc., might have influenced the CH₄ production.

${\mathrm{CH}}_{\mathrm{4}}/{\mathrm{CO}}_{\mathrm{2}}$: this parameter represents the ratio of CH₄ to CO₂. The parameter was slightly increased to a value of 0.14 from 0.085 following optimisation. This indicates a comparatively higher CH₄ emission fraction from the total decomposed carbon. All the parameters except for R_moist and $R_{{\mathrm{moist}}_{\mathrm{an}}}$ discussed above showed weak positive or neutral correlations with ${\mathrm{CH}}_{\mathrm{4}}/{\mathrm{CO}}_{\mathrm{2}}$. The weak positive correlations observed between most parameters and the ${\mathrm{CH}}_{\mathrm{4}}/{\mathrm{CO}}_{\mathrm{2}}$ ratio suggest that, while these factors may not be the primary drivers of CH₄ production, they still contribute to the increases in the ratio.

f_oxid: the fraction of oxidised CH₄, utilising available oxygen in the soil, is represented by the parameter f_oxid. The posterior parameter value (0.76) is increased compared to the prior (0.5). This indicates that a substantial fraction of the produced CH₄ will get oxidised, while the remaining CH₄ (24 %) will get transported to the atmosphere. This parameter has shown very low correlation with the other parameters and zero correlation with r_tiller and por_cato.

ϕ_tiller and r_tiller: the posterior parameter values estimated for ϕ_tiller and r_tiller are higher than the prior values, 0.77 and 0.0081, respectively. With aerenchyma tissues having more porous space and a larger radius, the plant-mediated transport of CH₄ to the atmosphere is facilitated. However, through the same spacious aerenchyma tissues, plants also have the potential to transport more O₂ to the soil. This potential increase in the transport of O₂ to the soil could be a reason for increase in f_oxid, considering the slight positive correlation observed between f_oxid and ϕ_tiller.

f_air, por_acro, and por_cato: these three parameters are related to soil composition. The posterior values of f_air increased compared to the model prior indicate a higher fraction of air in the soil. The decrease in por_acro and por_cato indicates decreased porosity in the acrotelm (which can contain both water or air) and catotelm (which can contain only water) respectively. f_air and por_acro are positively correlated, indicating more air in a more compact acrotelm environment with less water. A higher amount of air in the acrotelm can have a positive effect on diffusion. In the model, the diffusivity of CH₄ in air is estimated to be four orders of magnitude larger than in water. On the other hand, lower porosity in the catotelm can reduce ebullition due to the limited porous space in the soil, which holds less water. This decreases the capacity to retain excess CH₄ and release it when the solubility threshold is reached. (see Kallingal et al., 2024; Wania et al., 2010, for details).

λ_root: λ_root played a crucial role in this optimisation. After optimisation, this parameter got a significantly lower value of 10.25 cm compared to the prior. It seems that the optimisation, when generally trying to reduce the emission from the model, has a tendency to reduce the decay length of root biomass in the soil. The posterior parameter value closely aligns with the values reported in Kallingal et al. (2024) (10.47 cm). The optimisation results indicate a much shallower soil profile for the majority of root decay activities and CH₄. Given that most of the peat decomposition activities are assumed to occur in acrotelm, the reduction in the magnitude of λ_root could substantially facilitate diffusion. λ_root has not shown any strong correlation with any of the other parameters, but its weak negative correlation with r_moist, r_moistan, f_oxid, and f_air suggests that soils with shallow root depth exhibit increased moisture response and air fraction, with higher oxygen availability.

4.2 Posterior model estimates

Overall, the optimisation successfully reduced the model–data misfits, with some variation in the degree of improvement among different sites. This improvement is evident from the substantial magnitude change in the total cost function, RMSE and uncertainty from the prior to the posterior, and the estimated χ² values (see Tables C1 and C2, Fig. 4, and Table 5). It should be noted that this successful improvement occurred even when assimilating data from multiple sites with diverse climatic conditions.

The approximately 95 % improvement in the fit, as measured by the cost function, demonstrates the effectiveness of the GRaB-AM algorithm in sampling from high-probability regions. In general, such a large reduction in the magnitude of the cost function value can be interpreted as the model becoming overly tuned to the specific dataset, which risks a loss of generalisability. However, the high magnitude of the posterior cost, 79 296.4, along with the high chi-square value of 8.6, indicates a slight underfitting. A portion of this underfitting in posterior estimate can be attributed to the inability of GRaB-AM to capture the collective variability in such a large and diverse EC flux observations (see the Fig. 5). GRaB-AM assumes normally distributed fluxes, but in reality, observed fluxes often deviate from normality, exhibiting outliers, skewness, and heavy tails. However, the simplification to normality is a direct result of applying the commonly used quadratic loss function in data assimilation, and aids tractability of the MCMC-implementation. Introducing more complex error distributions would require estimation of the parameters of those distributions, significantly increasing both the number of parameters to estimate and the computational expense. Such an expansion, while possible and interesting future work, is outside the scope of this study.

Although LPJ-GUESS is a well-developed model, there are still several processes affecting the calculation of CH₄ emissions that are still not fully included in the model. For example, the assumption of zero wind speed above wetlands, the lack of detailed representation of ebullition, and the inadequate representation of wintertime emissions are some of the process in need of improvement (Kallingal et al., 2024). Another main reason for this underfitting could be the variability in the model input data. It has been observed that the model is highly sensitive to precipitation and temperature inputs. Temperature alone can explain a large proportion of the variation in the CH₄ simulations (Aalto et al., 2025), which lies beyond the constraints of GRaB-AM. And even though we used the meteorological measurements at the sites, there may still be biases in the representativeness of these measurements compared to the flux footprint.

Examining the sites individually, the four sites with a significant reduction in RMSE and a high χ² value imply enhanced accuracy in the model predictions but indicate similar underfitting as discussed before. More attention should be given to these sites, as there may be potential for further improvement in capturing the variability of the assimilated data. On the other hand, the seven sites showing the smallest RMSE reduction and low χ² values indicate less model improvement, but the model captures the variability in the assimilated data. The three sites with a high reduction in RMSE and low χ² values indicate good performance in model prediction and in capturing the variability. No sites have appeared with low RMSE reduction and high χ² values, which indicates that the predictions after optimisation have improved in at least one way or the other. No systematic trends are observed in these matrices, which indicates that the reason a few sites show considerable improvement in both RMSE and χ² is more due to the external factors such as model structure and assimilated or input data rather than the assimilation framework. Even though no correlations are observed between the types of wetlands and their locations, we note that the sites that showed a considerable reduction in RMSE, such as Abi, Att, Bon, Deg, Lom, Los, Sii, Uoa, and Zar, are those that are boreal or arctic in nature, missing only Sco and Zak. Except for Los, all sites located in the temperate region, namely Hue, Ole, Bib, and Zar, showed comparatively lower reductions in RMSE (see Table 2 and Fig. 4).

Compared to the total assimilated flux observations, the posterior estimate showed a slight underestimation of −38.1 g C m⁻² (see Sect. 3.4). The assimilated sites did not exhibit any consistent patterns of over or underestimation (Fig. 5). Testing the GRaB-AM algorithm at the Siikaneva site, Kallingal et al. (2024) observed systematic underestimation in LPJ-GUESS posterior over many years. Employing multiple sites with varying climatic variability has proven beneficial in resolving this issue, as sites like Bib, Deg, Sco, etc., exhibit both overestimation and underestimation in consecutive years. On the other hand, Kallingal et al. (2024) observed a considerable reduction in RMSE for Siikaneva in their single-site experiment, whereas in this study, being one among the 14 sites, the optimisation for Siikaneva was comparatively limited. This is rather expected in a multi-site setup, as the model must balance variability across multiple flux observations, which can degrade the fit at any individual site in favour of overall performance and generalisability.

Emission of CH₄ from wetlands can vary significantly across latitudes, largely influenced by the growing season and climatic conditions. When assimilating data from a single site into the LPJ-GUESS model, Kallingal et al. (2024) encountered difficulties in accurately capturing winter-time emissions. We faced a similar issue, where the model emitted zero CH₄ when temperatures dropped below freezing point. However, observational data showed that wetlands continue to emit CH₄ even under sub-zero conditions (Ito et al., 2023; Treat et al., 2018; Aalto et al., 2025). This discrepancy highlights a notable limitation of the LPJ-GUESS model: its high sensitivity to temperature inputs for CH₄ emissions and microbial CH₄ production and consumption being strongly inhibited in cold and frozen soils. The study by Ito et al. (2023), in which they compared cold-season CH₄ fluxes simulated by 16 models, shows that many similar models exhibit the issue with underestimating wintertime emissions. The same was observed in the study by Aalto et al. (2025), in which six models were compared. In both studies, LPJ-GUESS was a participant. One of the main observations by Ito et al. (2023) is that LPJ-GUESS is one of the models that discretely suppressed emissions under sub-zero temperatures and set a clear temperature threshold for CH₄ production. This could limit the GRaB-AM framework from adjusting the model's wintertime emissions and capturing the overall variability. Consequently, the algorithm compensates for the winter model–data mismatch by adjusting summer values. One potential model modification could involve incorporating mechanisms to simulate microbial activity under frozen conditions, snowpack insulation, and detailed representation of soil temperature dynamics, allowing for CH₄ emissions even when surface temperatures drop below zero. It has been observed that models representing detailed microbial activities and soil temperature profiles, such as the WETMETH model (Nzotungicimpaye et al., 2021), could simulate wintertime emissions on a comparable scale. A future study should involve using the GRaB-AM framework on such a somewhat simpler model to compare and verify its performance.

4.3 Validation of the posterior model and budget estimation of northern wetlands

Through selecting 14 different wetlands as described in Sect. 2.2, we aimed to equip the optimised parameters with the capability to accurately represent different types of wetlands, irrespective of their specific climatic and geographical features. Our validation analysis suggests that the optimised parameters achieved the goal to a large extent (see Sect. 3.5). In general, the optimised parameters perform better in representing different types of wetlands, especially bogs. A mismatch of 102 g C m⁻² is observed between the total flux observation and the total posterior model estimate, which is a large value compared to the result from the optimisation. Given this observation and the less constrained nature of sites like Wpt and Lgt, which are temperate, future studies might consider the necessity of different parameter sets for different wetland types. This is evident from the study of Treat et al. (2018), which shows that the behaviour of wetlands can be completely different based on their type and location, especially during the non-growing season. Here it is worth noting, as mentioned in Sect. 4.2, that the majority of temperate sites used for optimisation also exhibited lower fit improvement in terms of cost function.

There is a significant shortage of reliable, long-term CH₄ flux observations from wetlands, which is a major challenge for studies like ours. Many of the wetland sites used in this analysis only provided measurements during summer months, or contained substantial data gaps. This limitation makes it difficult for any assimilation algorithm to fully capture the natural variability in the data and could potentially lead to biased misinterpretations. Having used all the main sites with long-term measurements for assimilation, it was challenging to identify additional measurement sites for validation.

Figure 10 Time series of the total annual CH₄ estimates above 45° N. Prior and posterior CH₄ estimates from LPJ-GUESS, from 2010 to 2020 are shown in blue and green, respectively. CH₄ estimate from JSBACH-H, from 2010 to 2016, is shown in purple.

Both the assimilation and validation experiments showed a reduction in the posterior compared to the prior. The large scale simulation above 45° N also showed a similar tendency, with around a 12.87 % reduction in the posterior compared to the prior (see Figs. 9 and 10 and Sect. 3.5). We compared our results against CH₄ wetland emissions from JSBACH-H. This comparison primarily indicates that the change in the LPJ-GUESS simulated emissions after optimisation is in the direction towards better agreement with JSBACH-H fluxes, even though the comparison with JSBACH-H still suggests that LPJ-GUESS overestimates emissions. However, other studies, such as the one by Nzotungicimpaye et al. (2021), reported 24.825 Tg C yr⁻¹ from wetlands north of 45° N, which aligns much better with the posterior LPJ-GUESS estimate. Also, the study by Aalto et al. (2025) shows that the 16 wetland models they used simulated annual CH₄ emissions from wetlands of 45–90° N during 2000–2020 as 20.1 Tg C yr⁻¹. We also compared the resulting posterior to the observations quantified for the same domain. The gridded data product of northern wetland CH₄ emissions, based on upscaling EC flux observations estimated by Peltola et al. (2019), resulted in an average value of 25.98 Tg C yr⁻¹ (23.77–28.2 Tg C yr⁻¹). The results indicate that the LPJ-GUESS posterior is close to their quantified observation. However, the large uncertainty in various estimates of wetland CH₄ emissions makes it difficult to compare the results with other studies. For example, as mentioned in Sect. 3.5, when compared with the GCP estimate of wetland CH₄ emissions above 60° N, the prior emissions from LPJ-GUESS show better agreement. Nevertheless, the GCP estimate also reports a wide range, from 4.5 to 13.5 Tg C yr⁻¹, which would also cover the LPJ-GUESS posterior.

4.4 Possibilities and limitations of GRaB-AM

The GRaB-AM algorithm incorporates the adaptation mechanism with Rao–Blackwellisation, which recursively updates the covariance of the proposal distribution to capture the dependence among different parameters. Such a recursive update of the proposal distribution can improve the efficiency of MCMC by allowing the search algorithm to take larger steps in the parameter space, while still accounting for parameter inter-dependencies. This is particularly important for high-dimensional and correlated parameter spaces. In the optimisation process laid out here, with flux observations from multiple sites and with a relatively high-dimensional parameter space of a highly non-linear model, the GRaB-AM algorithm is particularly beneficial because it enhances the exploration of a wider parameter space while adapting the proposal distribution over iterations. The algorithm has the potential to incorporate other types of natural data distributions and account for temporal correlations in the data. Based on the results of this study and the results form Kallingal et al. (2024), it is evident that the algorithm is useful for optimising complex models but can also be easily applied to simpler, more linear models such as WETMETH (see Sect. 4.2 for the discussion). Additionally, there is potential for the algorithm to be used in optimising other model variables, such as CO₂ fluxes, vegetation, or soil dynamics, where long, reliable time series data are available.

While GRaB-AM is a powerful technique for highly non-linear models, there are several considerations and potential challenges when applying it to multi-site assimilation. (a) GRaB-AM can be computationally intensive, especially in high-dimensional parameter spaces or when dealing with a large number of sites. In this study, it took around 480 computational hours to complete 100 000 iterations on an AMD Ryzen Threadripper processor. (b) Adaptation mechanisms in GRaB-AM need to be carefully tuned to balance parameter exploration. Improper tuning may lead to suboptimal exploration of the parameter space. Additionally, the algorithm requires tuning hyperparameters that control the speed of adaptation, and its performance is sensitive to these choices, making it challenging to find appropriate hyperparameter values to begin with. (c) GRaB-AM assumes that observational data are Gaussian distributed, which is not always true when considering natural variability in multi-site flux data. The data distribution may be Lorentzian, Bernoulli, heavy-tailed, or take other forms. (d) Temporal correlations in the data are not addressed in the current form of GRaB-AM. (e) We encountered difficulties in efficiently minimising the total misfit, i.e. simultaneously reducing the misfit across multiple sites. This is due to the largely different cost function values for individual sites, caused by the varying conditions at each location. As a result, scaling was applied independently to each site to homogenise the contribution (weight) of each site in the optimisation process. However, this weighting of each site in a multi-site assimilation is a general challenge in multi-site data assimilation experiments and does not pertain to the GRaB-AM technique.

5 Conclusions

This study aimed to optimise the simulation of CH₄ emissions from natural wetlands in the LPJ-GUESS DGVM using eddy-covariance flux measurement data obtained from 14 diverse natural wetlands, characterised by variations in temporal, spatial, and/or climatic features. Ten selected model process parameters with the greatest influence on wetland CH₄ flux simulation were optimised using the Global Rao–Blackwellised Adaptive MCMC (GRaB-AM) algorithm within a Bayesian framework as a follow-up study of Kallingal et al. (2024). Following the optimisation, the study used flux observations from five different wetlands, which again differ in their temporal, spatial, and bioclimatic features, to validate the results of the optimisation. The optimisation results showed a substantial enhancement in the model's capacity to align with observed CH₄ fluxes, with a total reduction of approximately 50 % in RMSE and an approximately 53 % reduction in total uncertainty. The discrepancy between the modelled and observed values decreased from 1068.5 to 38.1 g C m⁻². For wetlands above 45 ° N, the total mean annual emission from the posterior LPJ-GUESS estimate is 7.46 Tg C yr⁻¹. Validation results demonstrate that four out of five sites reduced RMSE, contributing to an overall reduction of approximately 58.8 %. Given the remaining mismatches between observations and simulations and the presence of less constrained sites, future investigations will focus on individual sites, and grouping them based on their bio-geo-climatic characteristics, to examine if they need to be parameterised with different sets of parameters. Additionally, further studies are planned to quantify CH₄ emissions from boreal and temporal wetlands on large spatial scales, using the optimised parameters, and to validate them against independent atmospheric observations, i.e. atmospheric CH₄ observations provided by the European ICOS observation network. Another intended outcome of this study is to make use of the error correlation derived from the study as prior input to the atmospheric CH₄ inversion model, such as Lund University Modular Inversion Algorithm (LUMIA) (Monteil and Scholze, 2021).

Appendix A: Data source description

Among the sites used for assimilation, Bib (JP-BBY), Bon (US-BZB), Deg (SE-Deg), Hue (DE-Hte), Los (US-Los), Ole (US-ORv), Sco (CA-SCB), Uoa (US-Uaf), and Zar (DE-Zrk) were collected from Fluxnet datasets (the IDs in the bracket corresponds to Fluxnet) (https://fluxnet.org/data/fluxnet-ch4-community-product/, last access: 7 April 2025, Delwiche et al., 2021; Pastorello et al., 2020). For the site Abi, CH₄ data was collected from ICOS Sweden (2023) (https://doi.org/10.18160/Q6H6-B94B and the climate data was obtained from SMHI (https://www.smhi.se/data/meteorologi/, last access: 29 April 2025). Data for Att was collected from Ameriflux (https://doi.org/10.17190/AMF/1902822; Todd and Humphreys, 2025). For the site Lom (FI-Lom), climate observations of 2006 were taken from the Fluxnet site mentioned above. Observations for the remaining years were obtained from the station Principal Investigator (PI). Precipitation and temperature data for Sii were taken from FMI (https://en.ilmatieteenlaitos.fi/download-observations, last access: 7 April 2025), and CH₄ data and short-wave radiation data for Sii were collected from AVAA-SMEAR (https://smear.avaa.csc.fi/download, last access: 29 April 2025) (see Kallingal et al., 2024, for details). For the site Zak, the data was taken from GEM CH₄ (https://doi.org/10.17897/430P-DS31; Greenland Ecosystem Monitoring, 2025a) and GEM climate (https://doi.org/10.17897/9RY6-ZB75; Greenland Ecosystem Monitoring, 2025b). The data for Che, Lgt, Sfn, and Wpt, used for validation, were collected from the Fluxnet datasets mentioned above. Climate data for Myk was obtained from SITES (https://www.fieldsites.se/sv-SE, last access: 29 April 2025) and the CH₄ data were obtained from station PIs.

Appendix B: Statistical metrics

The following statistical metrics are used, among others, to analyse the results. The total reduced chi-square (χ²) is used to assess the goodness of the fit between the flux data and the model outputs. We calculated it by dividing twice the value of the cost function by the number of observations as

$$\begin{array}{}\text{(B1)}& {\mathit{\chi }}^{\mathrm{2}}={\displaystyle \frac{\mathrm{2}J\left(x\right)}{N}},\end{array}$$

where J is the cost function, x corresponds to the parameters, and N is the number of observations.

The root mean square error (RMSE) is estimated to quantify the average error between predicted and observed flux data as

$$\begin{array}{}\text{(B2)}& \text{RMSE}=\sqrt{{\displaystyle \frac{\mathrm{1}}{N}}\sum _{i=\mathrm{1}}^N(Y_i-M_i{)}^{\mathrm{2}}},\end{array}$$

where N is the number of observations, Y_i is the observed value, and M_i is the predicted value.

An uncertainty estimation was conducted assuming independence between parameter uncertainty and model uncertainty, using the following equation:

$$\begin{array}{}\text{(B3)}& {\mathit{\sigma }}_{\mathrm{total}}=\sqrt{{\mathit{\sigma }}_{\text{model}}^{\mathrm{2}}+{\mathit{\sigma }}_{\text{param}}^{\mathrm{2}}},\end{array}$$

where σ_model is the model structural uncertainty estimated from the standard deviation of the prior and posterior residuals (Desroziers et al., 2005). σ_param represents the contribution of parameter uncertainty to the overall uncertainty in observation space, estimated from the 95 % credible interval of the parameters and the standard deviation of total sums of the model prediction by taking into account both the parameter uncertainty from the MCMC sampling and the variability in the model predictions. The calculation is performed as follows:

$$\begin{array}{}\text{(B4)}& {\mathit{\sigma }}_{\text{param}}={\displaystyle \frac{{\mathit{\sigma }}_{\text{predic}}({\text{CI}}_{\text{upp}}-{\text{CI}}_{\text{low}})}{\mathrm{1.96}}},\end{array}$$

where σ_predic is the standard deviation of total sums of the model prediction over MCMC runs, CI_upp and CI_low are the upper and lower bounds of the credible intervals of the parameters, and the factor 1.96 is the conversion factor to convert the 95 % credible interval to a standard deviation assuming a Gaussian distribution.

In addition, we used skewness and kurtosis estimates of the posterior parameter probability density functions (PDFs) to describe their structure. We also calculated summer and winter anomalies by subtracting the mean values for April through September (summer) and October through March (winter) from the corresponding time series to estimate the variability in the seasonal cycle of the observed and modelled CH₄ fluxes.

Appendix C: Result of optimisation

C1 Changes in component contribution

Changes in the component-wise estimation of ebullition, diffusion, and plant-mediated transport before and after optimisation is illustrated in Fig. C1. The inner circles represent the priors, while the outer circles represent the posterior model estimates. The optimised parameters are constrained differently for different components across sites. There was no general trend of any one of the three transport mechanisms being dominant after optimisation. Regarding the prior estimate, all sites but Hue and Att had significant contributions from all three emission components. Hue and Att had a very minor contribution from plant-mediated transport. After optimisation, zero contributions from plant-mediated transport were estimated for both these sites. Interestingly, for the site Los, the majority of the prior was contributed by plant-mediated transport and ebullition. However, in the posterior, ebullition contributed very little and was taken over by diffusion. Furthermore, many sites showed the dominance of only two components after optimisation. It was consistently observed that the third component was suppressed, regardless of the nature or climatic conditions of the site.

Table C1 Posterior parameter estimates from the GRaB-AM optimisation. The “start prior values” are the initial random values sampled from the prior distributions. Maximum a posteriori (MAP) estimates, posterior means, and posterior standard deviations (std) are provided. The cost function values corresponding to the prior and posterior estimates are also given.

Table C2 A comprehensive overview of un-weighted (uw) prior and posterior cost values, RMSE reduction in percentages, and the χ² values estimated for all sites individually and together.

Figure C1 Component-wise percentage contribution of CH₄ to the total modelled emission for all 14 sites is presented separately. The inner circle represents the prior estimate, and the outer circle represents the posterior estimate.

Another interesting observation pertains to the site Zak, an arctic fen with a very low mean annual temperature (MAT), where nearly equal contributions from all three components were observed in the prior. The posterior, however, showed that nearly all emissions were from diffusion, with very little contributions from the two other components. The RMSE estimate for this site indicates a very low reduction compared to other sites, suggesting that the optimisation did not perform well in constraining this site.

The dominance of only one or two components after optimisation indicates possibilities of biases. Verifying and resolving this issue might be achievable through component-wise assimilation into the model using data from all three components and local hydrology observations. However, this will be challenging due to the unavailability of data, especially of the ebullition. Measuring ebullition fluxes poses significant challenges, primarily attributed to the pronounced spatio-temporal variability. Ecosystems exhibit rapid, momentary surges in fluxes, reaching exceptionally high levels within seconds, interspersed with prolonged periods of negligible ebullition (Canadell et al., 2021)

C2 Summer and winter anomaly estimation

Summer and winter anomalies of flux observations, prior, and posterior estimated separately for all 14 sites used can be seen in Fig. C2. The figure also provides details about the years in which the model either underestimated or overestimated the emissions. It is clear that neither the simulation nor the observation follows any common seasonal patterns or trends. This indicates that CH₄ emissions from wetlands are generally highly dependent on the variabilities in the underlying climatic variables, and the same holds for the model. A detailed analysis of the correlation and sensitivity between the model's CH₄ emission and input climatic variables can be seen in Kallingal et al. (2024).

High deviations were observed in the summer anomaly compared to the winter anomaly at all the sites. In general, in the majority of cases, the model was capable enough to capture trends shown by the observational anomaly, though there were differences in magnitude. For example, for the sites Bon and Zak, the model was successful in capturing all the summer and winter trends of the observation. Notably, the high positive anomaly of Abi in 2016 and of Ole in 2015, etc., and the high negative anomaly of Att in 2015, of Deg in 2017, of Low in 2017, etc., were also captured by the model.

C3 The annual sum estimation of CH₄

Figure C3 shows the annual sum estimation of CH₄ from 10 out of 14 sites used in the study. The remaining four sites are illustrated in Sect. 3.4 of the main paper. The figure illustrates that, after optimisation, most sites exhibited improved annual CH₄ estimation throughout the year. However, for the site Hue, the model consistently failed to capture the observation pattern in most years, except for 2013. For the site Att, particularly for the year 2016, the model displayed shortcomings. On the other hand, Att in 2016 completely aligned with the observed value for both prior and posterior estimations.

Figure C2 Summer and winter anomalies were estimated from the averages of the summer months (April to September) and winter months (October to March). The black, green, and purple dashed lines represent flux observations, prior, and posterior values, respectively. Dots and plus signs denote summer and winter data points of the season, respectively.

Figure C3 The annual sum estimation of CH₄ from 10 out of 14 sites used in the study. The sites are represented in different colours with distinct markings to distinguish between observation, prior, and posterior.

Appendix D: Time series estimation of validation sites

Figure D1 The prior and posterior CH₄ simulation from the LPJ-GUESS compared with flux observations from five different wetland sites used for validation. The black dots represent the CH₄ flux observations. The teal dots depict the prior simulation using the prior model parameters and the purple dots represent the posterior simulation; 3 d running averages are calculated from the original time series. In most of the figures, a few outliers on the vertical axis have been removed for better visualisation.