Monthly gridded data product of northern wetland methane emissions based on upscaling eddy covariance observations

Abstract

Natural wetlands constitute the largest and most uncertain source of methane (${\mathrm{CH}}_{\mathrm{4}}$) to the atmosphere and a large fraction of them are found in the northern latitudes. These emissions are typically estimated using process (“bottom-up”) or inversion (“top-down”) models. However, estimates from these two types of models are not independent of each other since the top-down estimates usually rely on the a priori estimation of these emissions obtained with process models. Hence, independent spatially explicit validation data are needed. Here we utilize a random forest (RF) machine-learning technique to upscale ${\mathrm{CH}}_{\mathrm{4}}$ eddy covariance flux measurements from 25 sites to estimate ${\mathrm{CH}}_{\mathrm{4}}$ wetland emissions from the northern latitudes (north of 45${}^{\circ }$ N). Eddy covariance data from 2005 to 2016 are used for model development. The model is then used to predict emissions during 2013 and 2014. The predictive performance of the RF model is evaluated using a leave-one-site-out cross-validation scheme. The performance (Nash–Sutcliffe model efficiency $=\mathrm{0.47}$) is comparable to previous studies upscaling net ecosystem exchange of carbon dioxide and studies comparing process model output against site-level ${\mathrm{CH}}_{\mathrm{4}}$ emission data. The global distribution of wetlands is one major source of uncertainty for upscaling ${\mathrm{CH}}_{\mathrm{4}}$. Thus, three wetland distribution maps are utilized in the upscaling. Depending on the wetland distribution map, the annual emissions for the northern wetlands yield 32 (22.3–41.2, 95 % confidence interval calculated from a RF model ensemble), 31 (21.4–39.9) or 38 (25.9–49.5) Tg(${\mathrm{CH}}_{\mathrm{4}}$) yr${}^{-\mathrm{1}}$. To further evaluate the uncertainties of the upscaled ${\mathrm{CH}}_{\mathrm{4}}$ flux data products we also compared them against output from two process models (LPX-Bern and WetCHARTs), and methodological issues related to ${\mathrm{CH}}_{\mathrm{4}}$ flux upscaling are discussed. The monthly upscaled ${\mathrm{CH}}_{\mathrm{4}}$ flux data products are available at10.5281/zenodo.2560163 (Peltola et al., 2019).

1 Introduction

Methane (${\mathrm{CH}}_{\mathrm{4}}$) is the second most important anthropogenic greenhouse gas (GHG) in terms of radiative forcing after carbon dioxide (${\mathrm{CO}}_{\mathrm{2}}$): 34 times (GWP${}_{\mathrm{100}}$, including climate-carbon feedbacks) as strong as ${\mathrm{CO}}_{\mathrm{2}}$ (Ciais et al., 2013). ${\mathrm{CH}}_{\mathrm{4}}$ has contributed $\sim \mathrm{20}$% to the cumulative GHG-related global warming (Etminan et al., 2016). Deriving constraints on ${\mathrm{CH}}_{\mathrm{4}}$ sources and sinks is thus of utmost importance. The net atmospheric ${\mathrm{CH}}_{\mathrm{4}}$ budget is well constrained by precise ${\mathrm{CH}}_{\mathrm{4}}$ mole fraction measurements around the globe, yet the contribution of individual sources and sinks to this aggregated budget remains poorly understood. This is primarily due to lack of data to constrain the modelling results (Saunois et al., 2016). In order to make more accurate predictions of the atmospheric ${\mathrm{CH}}_{\mathrm{4}}$ budget in a changing climate, the response of the various sources and sinks to different drivers needs to be better identified and quantified.

Natural wetlands are the largest and quantitatively most uncertain source of ${\mathrm{CH}}_{\mathrm{4}}$ to the atmosphere (Saunois et al., 2016). An ensemble of land surface models estimated global ${\mathrm{CH}}_{\mathrm{4}}$ emissions from wetlands for the period 2003–2012 to be 185 Tg(${\mathrm{CH}}_{\mathrm{4}}$) yr${}^{-\mathrm{1}}$ (range 153–227 Tg(${\mathrm{CH}}_{\mathrm{4}}$) yr${}^{-\mathrm{1}}$) and for the same period inversion models estimated it to be 167 Tg(${\mathrm{CH}}_{\mathrm{4}}$) yr${}^{-\mathrm{1}}$ (range 127–202 Tg(${\mathrm{CH}}_{\mathrm{4}}$) yr${}^{-\mathrm{1}}$) (Saunois et al., 2016). This discrepancy between bottom-up (process model) and top-down (inversion model) estimates, as well as the range of variability, exemplifies the large uncertainty of the current estimate for natural wetland ${\mathrm{CH}}_{\mathrm{4}}$ emissions. Sources of this uncertainty can be roughly divided into two categories: (1) uncertainty related to the global areal extent of wetlands (e.g. Petrescu et al., 2010; Bloom et al., 2017a; Zhang et al., 2016) and (2) uncertainties related to the key ${\mathrm{CH}}_{\mathrm{4}}$ emission drivers and responses to these drivers (e.g. Bloom et al., 2017a; Saunois et al., 2017). Evaluation of the emission estimates is thus urgently needed, and results from these efforts will lead to refined process models. Process model improvements will also directly affect the uncertainty of inversion results since they provide important a priori information to the inversion models (Bergamaschi et al., 2013).

Boreal and arctic wetlands comprise up to 50 % of the total global wetland area (e.g. Lehner and Döll, 2004) and the wetlands in these northern latitudes substantially contribute to total terrestrial wetland ${\mathrm{CH}}_{\mathrm{4}}$ emissions (ca. 27 %, based on the sum of regional budgets for boreal North America, Europe and Russia in Saunois et al., 2016). In wetlands, ${\mathrm{CH}}_{\mathrm{4}}$ is produced by methanogenic Archaea under anaerobic conditions, and hence the production takes place predominantly under water-saturated conditions (e.g. Whalen, 2005). The microbial activity and the resulting ${\mathrm{CH}}_{\mathrm{4}}$ production is thus controlled by the quality and quantity of the available substrates, competing electron acceptors, and temperature (Le Mer and Roger, 2001). Once produced, the ${\mathrm{CH}}_{\mathrm{4}}$ can be emitted to the atmosphere via three pathways: ebullition, molecular diffusion through soil matrix and water column, or plant transport. If plants capable of transporting ${\mathrm{CH}}_{\mathrm{4}}$ are present, plant transport is generally the dominating emission pathway (Knoblauch et al., 2015; Kwon et al., 2017; Waddington et al., 1996; Whiting and Chanton, 1992). A large fraction of ${\mathrm{CH}}_{\mathrm{4}}$ transported via molecular diffusion is oxidized into ${\mathrm{CO}}_{\mathrm{2}}$ by methanotrophic bacteria in the aerobic layers of wetland soils and hence never reaches the atmosphere (Sundh et al., 1995), whereas ${\mathrm{CH}}_{\mathrm{4}}$ transported via ebullition and plant transport can largely bypass oxidation (Le Mer and Roger, 2001; McEwing et al., 2015). Furthermore, processes related to permafrost (e.g. active layer, thermokarst) and snow cover dynamics (e.g. snow melt, insulation) have an impact on ${\mathrm{CH}}_{\mathrm{4}}$ flux seasonality and variability (Friborg et al., 1997; Helbig et al., 2017; Mastepanov et al., 2008; Zona et al., 2016; Zhao et al., 2016). Hence, wetland ${\mathrm{CH}}_{\mathrm{4}}$ emissions to the atmosphere largely depend on interplay between various controls, including water table position, temperature, vegetation composition, methane consumption, availability of substrates and competing electron acceptors.

During the past 2 decades, eddy covariance (EC) measurements of wetland ${\mathrm{CH}}_{\mathrm{4}}$ emissions have become more common, due to rapid development in sensor technology (e.g. Detto et al., 2011; Peltola et al., 2013, 2014). The latest generation of low-power and low-maintenance instruments are rugged enough for long-term field deployment (Nemitz et al., 2018; McDermitt et al., 2010); thus, the number of sites where ${\mathrm{CH}}_{\mathrm{4}}$ flux measurements have been made is increasing. Due to this progress, EC ${\mathrm{CH}}_{\mathrm{4}}$ flux synthesis studies have been emerging (Petrescu et al., 2015; Knox et al., 2019). Similar progress was made with ${\mathrm{CO}}_{\mathrm{2}}$ and energy flux measurements in the 1990s and now these measurements form the backbone of the global EC measurement network FLUXNET (https://fluxnet.fluxdata.org/, last access: 6 August 2019), whose data have provided invaluable insights into terrestrial carbon and water cycles. Some of the most important results have been obtained by upscaling FLUXNET observations using machine-learning algorithms to evaluate terrestrial carbon balance components and evapotranspiration (Beer et al., 2010; Bodesheim et al., 2018; Jung et al., 2010, 2011, 2017; Mahecha et al., 2010). These results are now widely used by the modelling community to evaluate process model performance (e.g. Wu et al., 2017) and to validate satellite-derived carbon cycle data products (e.g. Sun et al., 2017; Y. Zhang et al., 2017).

In this study, we synthesized EC ${\mathrm{CH}}_{\mathrm{4}}$ flux data from 25 EC ${\mathrm{CH}}_{\mathrm{4}}$ flux sites and developed an observation-based monthly gridded data product of northern wetland ${\mathrm{CH}}_{\mathrm{4}}$ emissions. We focus on northern wetlands (north of 45${}^{\circ }$ N) due to their significance in the global ${\mathrm{CH}}_{\mathrm{4}}$ budget and relatively good data coverage and process understanding, at least compared to tropical systems (Knox et al., 2019). High-latitude regions are projected to warm during the next century at a faster rate than any other region, which will likely significantly impact the carbon cycling of wetland ecosystems (Tarnocai, 2009; Z. Zhang et al., 2017) and permafrost areas of the Arctic boreal region (Schuur et al., 2015). To date, ${\mathrm{CH}}_{\mathrm{4}}$ emission estimates for northern wetlands are typically based on process models (Bohn et al., 2015; Bloom et al., 2017a; Chen et al., 2015; Melton et al., 2013; Stocker et al., 2013; Wania et al., 2010; Watts et al., 2014; Zhang et al., 2016) or inversion modelling (Bohn et al., 2015; Bruhwiler et al., 2014; Spahni et al., 2011; Thompson et al., 2017; Thonat et al., 2017; Warwick et al., 2016), yet scaling of existing chamber measurements to the northern wetland area has also been published (Zhu et al., 2013). However, ${\mathrm{CH}}_{\mathrm{4}}$ emission estimates obtained with the former two approaches are not independent since the attribution of ${\mathrm{CH}}_{\mathrm{4}}$ emissions derived using inversion models to different emission sources (e.g. wetlands) depends largely on a priori estimates of these emissions (i.e. process models for wetland emissions), highlighting the tight coupling between these two approaches (Bergamaschi et al., 2013; Spahni et al., 2011). Hence, the main objective of this study is to produce an independent data-driven estimate of northern wetland ${\mathrm{CH}}_{\mathrm{4}}$ emissions. This product could be used as an additional constraint for the wetland emissions and hence aid in process model refinement and development. Additionally, the drivers causing ${\mathrm{CH}}_{\mathrm{4}}$ flux variability at the ecosystem scale are also evaluated and methodological issues are discussed, which will support future ${\mathrm{CH}}_{\mathrm{4}}$ wetland flux upscaling studies.

2 Materials and methods

Data from flux measurement sites (Fig. 1) were acquired and used together with forcing data to estimate ${\mathrm{CH}}_{\mathrm{4}}$ emissions from northern wetlands with a monthly time resolution using a random forest (RF) modelling approach. Both in situ measurements and remote sensing are utilized in this study. In this section, the RF approach is briefly introduced (Sect. 2.1) and data selection, quality filtering, gap filling and aggregation to monthly values are described (Sect. 2.3). We identified 40.7 site years available for analysis, measured between 2005 and 2016. To perform upscaling to all wetlands north of 45${}^{\circ }$ N, gridded data products of the flux drivers and wetland distribution maps were needed. These products are presented in Sect. 2.4 and 2.5, respectively. Finally, the upscaled wetland ${\mathrm{CH}}_{\mathrm{4}}$ emissions are compared against process model outputs, with the models briefly described in Sect. 2.6.

2.1 Random forest algorithm

Random forest (RF) is a machine-learning algorithm that can be used for classification or regression analyses (Breiman, 2001). In this study the RF models consist of a large ensemble of regression trees. Each individual regression tree is built by training it with a random subset of training data and the trees are trained independently of each other. The RF model output is then the average of all the predictions made by individual regression trees in the forest. Hence, the RF algorithm applies the bootstrap aggregation (bagging) algorithm and takes full advantage of the fact that ensemble averaging decreases the noise of the prediction. In addition to random selection of training data, the predictor variables used in split nodes are also selected from a random sample of all predictors, which minimizes the possible correlation between trees in the forest (Breiman, 2001) and decreases the possibility of overfitting. The predictor variables can be either categorical or continuous. The variables are then used in the split nodes to divide the data into two (e.g. categorical variable true or false or a continuous variable, such as temperature above or below 5 ${}^{\circ }$C).

Performance of RF algorithms to predict ${\mathrm{CO}}_{\mathrm{2}}$ and energy fluxes across FLUXNET sites have been compared against other machine-learning algorithms, such as artificial neural networks and multivariate regression splines, by Tramontana et al. (2016), who showed that differences between methods were negligible. We anticipate a similarly negligible effect of machine-learning algorithm choice for ${\mathrm{CH}}_{\mathrm{4}}$ fluxes. For a thorough description of the RF algorithm for flux upscaling purposes, the reader is referred to Bodesheim et al. (2018) (and references therein).

In this study, the RF models were developed using the MATLAB 9.4.0 (R2018a) TreeBagger function with default values similar to those of Bodesheim et al. (2018). These settings included a minimum of five samples in a leaf node and used mean squared error (MSE) as a metric for deciding the split criterion in split nodes. Each trained forest consisted of 300 randomized regression trees.

2.1.1

Description of input variables for RF model development for upscaling. Data were aggregated to monthly values (see text) unless otherwise noted below.

[Figure omitted. See PDF]

We further tested whether information about wetland type or sedge cover would improve the model performance even though these categorical variables were not available in gridded format and hence were not usable for upscaling. Including the sedge flag increased the NSE to 0.53, although the increase in Pearson correlation was not statistically significant ($p>\mathrm{0.05}$, comparison of correlation coefficients using Fisher's $r$ to $z$ transformation). Also, wetland type did not have a statistically significant influence on the model performance ($p>\mathrm{0.05}$ and NSE $=\mathrm{0.49}$ if type included). Using too many categorical variables in a RF model may be problematic because each site may end up with a unique combination of categorical variables.

The most important predictor for the model was temperature, similar to numerous studies showing that wetland ${\mathrm{CH}}_{\mathrm{4}}$ emissions are strongly correlated to soil temperature (Christensen et al., 2003; Helbig et al., 2017; Jackowicz-Korczyński et al., 2010; Rinne et al., 2018; Yvon-Durocher et al., 2014; Knox et al., 2019). Selection of LSTn as the primary driver instead of the other temperature variables was likely an outcome of the available data and the algorithm used to select the drivers. With slightly different dataset (more sites) other temperature variables (e.g. $T_{\mathrm{air}}$) might have been more important drivers for the ${\mathrm{CH}}_{\mathrm{4}}$ flux variability. Estimating apparent $Q_{\mathrm{10}}$ from the RF model LSTn dependence yielded a value of $\mathrm{1.90}\pm \mathrm{0.03}$ and for validation data it was slightly higher ($\mathrm{1.97}\pm \mathrm{0.06}$) (Fig. 3). These values are comparable to the ones reported in Turetsky et al. (2014) for ${\mathrm{CH}}_{\mathrm{4}}$ chamber measurements at bog and fen sites. The temperature dependence of ${\mathrm{CH}}_{\mathrm{4}}$ production is modelled in many process models with the parameter $Q_{\mathrm{10}}$ value close to 2 (Xu et al., 2016b), which agrees with the ${\mathrm{CH}}_{\mathrm{4}}$ emission temperature dependence shown here. However, one should note that ${\mathrm{CH}}_{\mathrm{4}}$ oxidation also depends on temperature and the derived apparent $Q_{\mathrm{10}}$ value describes the temperature dependence of surface ${\mathrm{CH}}_{\mathrm{4}}$ emission, which is always a combination of ${\mathrm{CH}}_{\mathrm{4}}$ production and oxidation.

Figure 3

Dependence of monthly mean ${\mathrm{CH}}_{\mathrm{4}}$ emissions on monthly mean land surface temperature at night (LSTn) derived from MODIS data. Eddy covariance measurements are shown with filled markers (unique colour for each site) and random forest model predictions for each site are given with black dots.

[Figure omitted. See PDF]

3.2 Model agreement with validation data

The overall systematic bias (BE) between the RF predictions and validation data was negligible (Fig. 4), whereas the spread of the data (RE) was more pronounced (Fig. 4). Following Moffat et al. (2010), RE was analysed further by binning the data based on ${\mathrm{CH}}_{\mathrm{4}}$ flux magnitude and calculating RE for each bin. RE was clearly correlated with flux magnitude ($\text{RE}=(\mathrm{0.52}\pm \mathrm{0.06}){\mathrm{FCH}}_{\mathrm{4}}+(\mathrm{3.3}\pm \mathrm{2.0})$ nmol m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$, where ${\mathrm{FCH}}_{\mathrm{4}}$ denotes ${\mathrm{CH}}_{\mathrm{4}}$ flux), indicating that the relative random error of the RF model prediction was nearly constant and approximately 50 % for high fluxes. The systematic error BE did not show a clear dependence on flux magnitude. The RF model performance was worse on a site mean level than with monthly data. When comparing site means, NSE and R${}^{\mathrm{2}}$ were both 0.25 and RE and BE were 27.0 and 1.5 nmol m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$, respectively. Possible drivers causing the remaining ${\mathrm{CH}}_{\mathrm{4}}$ flux variability not captured by the RF model (i.e. the scatter in Fig. 4) are discussed in Sect. 4.2.1.

Figure 5

Time series of modelled ${\mathrm{CH}}_{\mathrm{4}}$ emissions (red lines) together with validation data (circles) at four example sites: (a) Siikaneva oligotrophic fen in Finland, (b) Lost Creek shrub fen in Wisconsin, USA, (c) Atqasuk wet tundra in Alaska, USA, and (d) Chersky wet tundra in northeastern Siberia, Russia. Dashed vertical lines denote a new year. Note the changes in $y$ axis scales. Site-specific model performance metrics are also included.

[Figure omitted. See PDF]

The mean annual cycle of ${\mathrm{CH}}_{\mathrm{4}}$ emission predicted by the RF model agrees well with the mean annual cycle calculated from the validation data (not shown). During the non-growing season the RF model slightly overestimates the fluxes (15 % overestimation) but such differences were negligible during the rest of the year ($<\mathrm{1}$ %). However, for individual sites ${\mathrm{CH}}_{\mathrm{4}}$ emission seasonality agrees less. For instance, at US-Los the modelled ${\mathrm{CH}}_{\mathrm{4}}$ emissions start to increase 1 month earlier in the spring (Fig. 5b). The non-growing season fluxes are overestimated at four example sites (FI-Sii, US-Los, US-Atq and RU-Ch2; Fig. 5). The mean flux magnitude is modelled well at FI-Sii (Fig. 5a), whereas at US-Los (Fig. 5b) and US-Atq (Fig. 5c) the RF model overestimates and at RU-Ch2 (Fig. 5d) underestimates the ${\mathrm{CH}}_{\mathrm{4}}$ emissions. The flux bias had a relatively large impact on site-specific NSE. For example, for US-Atq NSE was $-\mathrm{1.85}$, meaning that the observation mean would be a better predictor for this site than the RF model (see the NSE definition in Sect. 2.2). The RF model is not able to replicate the between-year differences in ${\mathrm{CH}}_{\mathrm{4}}$ emissions at the example sites. Capturing interannual variability has also been difficult in previous upscaling studies of ${\mathrm{CO}}_{\mathrm{2}}$ and energy fluxes (e.g. Tramontana et al., 2016).

In general, the RF model performance was better for permafrost-free sites than for sites with permafrost ($r=\mathrm{0.66}$ and $r=\mathrm{0.51}$, respectively; $p<\mathrm{0.05}$), which is likely related to the fact that at sites with permafrost the MODIS LSTn is not as directly related to the soil temperature than at sites without permafrost. Hence, LSTn is not as good proxy for the temperature that is controlling both ${\mathrm{CH}}_{\mathrm{4}}$ production and consumption and this results in a worse performance than at sites without permafrost.

3.3

Upscaled ${\mathrm{CH}}_{\mathrm{4}}$ fluxes

The RF model developed in this study was used together with the gridded input datasets (Sect. 2.4) and wetland distribution maps (Sect. 2.5) to estimate ${\mathrm{CH}}_{\mathrm{4}}$ emissions from northern wetlands in 2013 and 2014. The mean ${\mathrm{CH}}_{\mathrm{4}}$ emissions of the 2 years from the RF model are plotted in Fig. 6 together with ${\mathrm{CH}}_{\mathrm{4}}$ wetland emission maps from the process model LPX-Bern and model ensemble WetCHARTs. Differences between the process model estimations and upscaled fluxes are shown in Fig. 7. In general, the spatial patterns are similar among emission maps, which is not surprising given that the spatial variability is largely controlled by the underlying wetland distributions. One noteworthy difference is that WetCHARTs, RF-PEATMAP (i.e. RF modelling with PEATMAP) and RF-GLWD show higher emissions from western Canada than LPX-Bern or the upscaled fluxes using the wetland map from that process model (RF-DYPTOP). The other difference is that RF-GLWD show negligible emissions from Fennoscandia (Fig. 6c). These differences are related to differences in the underlying wetland maps. While the wetland maps differ, there is no consensus on which is more accurate, so comparisons indicate the uncertainty in upscaling emanating from uncertainties in wetland distribution.

Figure 6

Mean annual ${\mathrm{CH}}_{\mathrm{4}}$ wetland emissions during years 2013–2014 estimated by upscaling EC data using the RF model and three wetland maps (a, b, c) and process models (d, e). Grid cells with low ${\mathrm{CH}}_{\mathrm{4}}$ wetland emissions (below 0.1 g(${\mathrm{CH}}_{\mathrm{4}}$) m${}^{-\mathrm{2}}$ yr${}^{-\mathrm{1}}$) are shown in grey. The flux rates refer to total unit area in a grid cell.

[Figure omitted. See PDF]

Figure 7

Difference in mean annual ${\mathrm{CH}}_{\mathrm{4}}$ wetland emissions during the years 2013–2014 estimated by upscaling EC data using the RF model with different wetland maps and process models. All the ${\mathrm{CH}}_{\mathrm{4}}$ emission maps were aggregated to 1${}^{\circ }$ resolution before comparison. The flux rates refer to total unit area in a grid cell.

[Figure omitted. See PDF]

Three statistical metrics (NSE, $R^{\mathrm{2}}$ and RE) were calculated between RF-DYPTOP and LPX-Bern for each grid cell (Fig. 8). The figure illustrates how well the temporal variability of ${\mathrm{CH}}_{\mathrm{4}}$ emissions estimated by RF-DYPTOP and LPX-Bern agree in each grid cell. NSE values are low in areas where the systematic difference between RF-DYPTOP and LPX-Bern was high (compare Figs. 8a and 7a) since the bias strongly penalizes NSE. The $R^{\mathrm{2}}$ values are high throughout the study domain, likely due to the fact that the seasonal cycle of ${\mathrm{CH}}_{\mathrm{4}}$ emissions dominated the temporal variability in most of the grid cells and the seasonal cycles were in phase between RF-DYPTOP and LPX-Bern. RE values calculated between RF-DYPTOP and LPX-Bern were high in areas where the emissions estimated by RF-DYPTOP were also high (compare Figs. 8c and 6a). This is likely due to the fact that, even though the seasonal cycles were in phase, their amplitudes were different which increased the variability between LPX-Bern and RF-DYPTOP (i.e. increase in RE).

Figure 8

NSE, $R^{\mathrm{2}}$ and RE calculated between RF-DYPTOP and LPX-Bern. Grid cells with low ${\mathrm{CH}}_{\mathrm{4}}$ wetland emissions (below 0.1 g(${\mathrm{CH}}_{\mathrm{4}}$) m${}^{-\mathrm{2}}$ yr${}^{-\mathrm{1}}$) are shown in grey. RE values refer to total unit area in a grid cell.

[Figure omitted. See PDF]

The uncertainties of the upscaled fluxes were estimated from the spread of predictions made with the ensemble of 200 RF models (Fig. 9). The uncertainty mostly scales with the flux magnitude (compare Fig. 6a–c with Fig. 9a–c), meaning that grid cells with high fluxes tend to also have high uncertainties. However, the relative flux uncertainty does have some geographical variation (Fig. 9d–f). The highest relative uncertainties are typically at the highest and lowest latitudes of the study domain. In these locations the dependencies between the predictors and the ${\mathrm{CH}}_{\mathrm{4}}$ flux are not as well defined as in the locations with lower uncertainties leading to larger spread in the ensemble of RF model prediction. For instance, at low latitudes LSTn may go beyond the range of LSTn values in the training data (see the range in Fig. 3), and hence the RF model predictions are not well constrained in these situations. On the other hand, lower relative uncertainties are typically obtained for locations close to the measurement sites incorporated in this study (compare Figs. 1 and 9), since the dependencies between the predictors and the ${\mathrm{CH}}_{\mathrm{4}}$ flux are better defined.

Figure 9

Absolute (a–c) and relative (d–f) uncertainties of the upscaled ${\mathrm{CH}}_{\mathrm{4}}$ fluxes using different wetland maps. Uncertainty is estimated as 1$\mathit{\sigma }$ variability of the predictions by 200 RF models developed by bootstrapping the training data (Sect. 2.1.2). Grid cells with low ${\mathrm{CH}}_{\mathrm{4}}$ wetland emissions (below 0.1 g(${\mathrm{CH}}_{\mathrm{4}}$) m${}^{-\mathrm{2}}$ yr${}^{-\mathrm{1}}$) are shown in grey. The absolute uncertainties refer to total unit area in a grid cell.

[Figure omitted. See PDF]

The seasonalities of the upscaled fluxes and ${\mathrm{CH}}_{\mathrm{4}}$ fluxes from process models are similar with the highest ${\mathrm{CH}}_{\mathrm{4}}$ emissions in July–August and the lowest in February. This seasonal pattern is consistent throughout the study domain (Fig. 10). Warwick et al. (2016) and Thonat et al. (2017) showed that the northern wetland ${\mathrm{CH}}_{\mathrm{4}}$ emissions should peak in August–September in order to correctly explain the seasonality of atmospheric ${\mathrm{CH}}_{\mathrm{4}}$ mixing ratios and isotopes measured across the Arctic. Hence, the wetland ${\mathrm{CH}}_{\mathrm{4}}$ emissions presented here are peaking approximately 1 month too early to perfectly match with their findings. ${\mathrm{CH}}_{\mathrm{4}}$ flux magnitude agrees well between WetCHARTs and the upscaled flux during spring and midsummer (April–July), whereas LPX-Bern estimates lower fluxes (0 % and 26 % difference, respectively). During late summer and autumn (August–October) both process models estimate slightly lower fluxes than the upscaled estimate (17 % and 19 % difference, respectively). The upscaled fluxes also show somewhat higher emissions during the non-growing season (November–March) than the two process models (27 % and 35 % difference; see Table 2), and the upscaled estimates of non-growing season emissions are relatively close to a recent model estimate (Treat et al., 2018). This result promotes the recent notion that process models might be underestimating non-growing season fluxes at high latitudes (e.g. Treat et al., 2018; Xu et al., 2016a; Zona et al., 2016).

Figure 10

Monthly time series of zonal mean ${\mathrm{CH}}_{\mathrm{4}}$ fluxes. The upscaled fluxes with different wetland maps are shown in (a, b, c) and wetland ${\mathrm{CH}}_{\mathrm{4}}$ emissions estimated with the two process models are given in (d, e).

[Figure omitted. See PDF]

Table 2

Annual ${\mathrm{CH}}_{\mathrm{4}}$ wetland emissions in different subdomains (Hudson Bay lowlands and western Siberian lowlands; see Fig. 1) and time periods. The values are given in Tg(${\mathrm{CH}}_{\mathrm{4}}$) yr${}^{-\mathrm{1}}$. Note that estimates from some reference studies are not for the same period as the one studied here (2013–2014). For WetCHARTs the mean of the model ensemble together with the range (in parentheses) is given, whereas for the upscaling results the 95 % confidence intervals for the estimated emissions are given.

	Reference	Hudson Bay	Western	Non-growing	Annual
		lowlands	Siberian	season fluxes	emissions
			lowlands	from northern	north of
				wetlands	45${}^{\circ }$ N
				(November–March)
Inversion	Bohn et al. (2015), WETCHIMP-WSL		$\mathrm{6.06}\pm \mathrm{1.22}$
models	Bruhwiler et al. (2014)${}^{\mathrm{a}}$				23
	Kim et al. (2011)		$\mathrm{2.9}\pm \mathrm{1.7}$ and $\mathrm{3.0}\pm \mathrm{1.4}$
	Miller et al. (2014)	$\mathrm{2.4}\pm \mathrm{0.3}$
	Spahni et al. (2011)				$\mathrm{28.2}\pm \mathrm{2.2}$
	Thompson et al. (2017)	2.7–3.4	$\mathrm{6.9}\pm \mathrm{3.6}$
Process	Bohn et al. (2015), WETCHIMP-WSL		$\mathrm{5.34}\pm \mathrm{0.54}$
models	Chen et al. (2015)${}^{\mathrm{b}}$	$\mathrm{3.11}\pm \mathrm{0.45}$			$\mathrm{35.0}\pm \mathrm{6.7}$
	Melton et al. (2013), WETCHIMP${}^{\mathrm{c}}$	$\mathrm{5.4}\pm \mathrm{3.2}$
	Pickett-Heaps et al. (2011)${}^{\mathrm{d}}$	$\mathrm{2.3}\pm \mathrm{0.3}$
	Treat et al. (2018)${}^{\mathrm{e}}$			$\mathrm{6.1}\pm \mathrm{1.5}$	$\mathrm{37}\pm \mathrm{7}$
	Watts et al. (2014)				53
	Zhang et al. (2016)${}^{\mathrm{f}}$	$\mathrm{5.5}\pm \mathrm{1.1}$	$\mathrm{4.6}\pm \mathrm{0.6}$		$\mathrm{30.3}\pm \mathrm{5.4}$
	This study, LPX-Bern	2.5	4.4	4.5	24.7
	This study, WetCHARTs	2.8 (0.5–8.7)	4.2 (1.6–9.4)	5.1 (0.6–17.0)	29.7 (8.7–74.0)
Flux	Glagolev et al. (2011)		$\mathrm{3.9}\pm \mathrm{1.3}$
measurement	Zhu et al. (2013)				44.0–53.7
upscaling	This study, RF-PEATMAP	4.8 (3.3–6.3)	6.6 (4.9–8.4)	6.7 (4.9–8.5)	31.7 (22.3–41.2)
	This study, RF-DYPTOP	4.6 (3.1–6.0)	7.0 (5.2–8.8)	6.2 (4.6–7.8)	30.6 (21.4–39.9)
	This study, RF-GLWD	4.9 (3.4–6.5)	6.8 (5.0–8.5)	8.0 (5.8–10.2)	37.6 (25.9–49.5)

${}^{\mathrm{a}}$ Approximately north of 47${}^{\circ }$ N. ${}^{\mathrm{b}}$ Approximately north of 45${}^{\circ }$ N. ${}^{\mathrm{c}}$ Mean annual ${\mathrm{CH}}_{\mathrm{4}}$ emissions from eight models $\pm$1$\mathit{\sigma }$ of interannual variation in the model estimates for the period 1993–2004. ${}^{\mathrm{d}}$ Process model tuned to match atmospheric observations. ${}^{\mathrm{e}}$ North of 40${}^{\circ }$ N. ${}^{\mathrm{f}}$ Mean $\pm$1$\mathit{\sigma }$ over the LPJ-wsl model results using different wetland extents for the period 1980–2000.

Treat et al. (2018) adjusted WetCHARTs model output so that it matches with their estimates of non-growing season ${\mathrm{CH}}_{\mathrm{4}}$ emissions and then estimated annual wetland ${\mathrm{CH}}_{\mathrm{4}}$ emissions north of 40${}^{\circ }$ N to be $\mathrm{37}\pm \mathrm{7}$ Tg(${\mathrm{CH}}_{\mathrm{4}}$) yr${}^{-\mathrm{1}}$ using this adjusted model output. The estimates derived here for the annual emissions using the three wetland maps are similar (see Table 2), especially when considering our slightly smaller study domain (above 45${}^{\circ }$ N). The two process models included in this study estimated slightly lower mean annual emissions than the upscaled fluxes (11 % and 26 % difference between the mean upscaled estimate and WetCHARTs and LPX-Bern, respectively; see also Table 2). However, given the uncertainties in upscaling as well as in process models, this can be regarded as relatively good agreement. Different process models may be driven with different climate forcing data and they may have discrepancies in the underlying wetland distributions, in addition to the different parameterizations and descriptions of the processes behind the ${\mathrm{CH}}_{\mathrm{4}}$ emissions. These sources of uncertainty should be recognized when models are compared against each other or against upscaling products.

In order to further evaluate the agreement between the upscaled fluxes and process models we focused on two specific regions: Hudson Bay lowlands (HBL) and western Siberian lowlands (WSL) (see locations in Fig. 1). The upscaled fluxes indicate higher annual emissions for both subdomains compared to the two process models or previously published estimate (Table 2). For WSL the upscaled estimates are within the range of variability observed between process models and inversion modelling in WETCHIMP-WSL (Bohn et al., 2015) and close to Thompson et al. (2017). The upscaled estimates by Glagolev et al. (2011) might underestimate ${\mathrm{CH}}_{\mathrm{4}}$ emissions from the WSL area (Bohn et al., 2015). Furthermore, the process models in Bohn et al. (2015) are likely underestimating the non-growing season ${\mathrm{CH}}_{\mathrm{4}}$ emissions which might partly explain the discrepancy to the upscaled estimates in this study. Hence, the upscaled ${\mathrm{CH}}_{\mathrm{4}}$ emission estimates for the WSL area, while large, are still in a reasonable range.

For HBL, the discrepancy between upscaled emission estimates and the estimates based on process models or previous studies is larger (Table 2). The upscaling results agree with Zhang et al. (2016) and Melton et al. (2013) but show emissions that are twice as large for the HBL than the other estimates (Table 2). This cannot be explained by wetland mapping since the difference also holds when the DYPTOP wetland map is used in upscaling. There are only few long-term EC flux studies conducted in the HBL area and the only one found (Hanis et al., 2013) showed on average 6.9 g(${\mathrm{CH}}_{\mathrm{4}}$) m${}^{-\mathrm{2}}$ annual emissions at a subarctic fen located in the HBL. If the upscaled ${\mathrm{CH}}_{\mathrm{4}}$ emissions are downscaled back to ecosystem level in the HBL area with wetland maps, we get on average 11.0 g(${\mathrm{CH}}_{\mathrm{4}}$) m${}^{-\mathrm{2}}$ annual ${\mathrm{CH}}_{\mathrm{4}}$ emission for the HBL area based on the RF model output, which is 1.6 times larger than the estimate by Hanis et al. (2013). While Hanis et al. (2013) studied only one wetland during different years than those used here (years 2008–2011 in Hanis et al., 2013, here 2013–2014), it is still noteworthy that the relative difference between Hanis et al. (2013) and this study is similar to the discrepancy between this study and the inversion estimates (Pickett-Heaps et al., 2011; Thompson et al., 2017) at the whole HBL scale. Pickett-Heaps et al. (2011) and Thompson et al. (2017) show near zero ${\mathrm{CH}}_{\mathrm{4}}$ emissions during October–April and onset of ${\mathrm{CH}}_{\mathrm{4}}$ emissions in mid-May or even June, largely dependent on when the ground was free of snow and unfrozen. This is somewhat surprising given the fact that only 32 % of wetlands in the area are underlain by permafrost (based on amalgam of PEATMAP and the permafrost map), and hence the soils are likely not completely frozen and some non-growing season ${\mathrm{CH}}_{\mathrm{4}}$ emissions are likely to occur in such conditions (e.g. Treat et al., 2018). The upscaled non-growing season ${\mathrm{CH}}_{\mathrm{4}}$ emissions show on average 1.1 Tg(${\mathrm{CH}}_{\mathrm{4}}$) yr${}^{-\mathrm{1}}$ emissions for the HBL area. This partly, but not completely, explains the discrepancy between the ${\mathrm{CH}}_{\mathrm{4}}$ emission estimates for the HBL area. All these results suggest that the upscaled product likely overestimates ${\mathrm{CH}}_{\mathrm{4}}$ emissions from the HBL area.

4 Discussion

4.1 Comparing the RF model predictive performance to previous studies

The RF model performance was worse when compared against independent validation data than what has been achieved in previous upscaling studies for GPP and energy fluxes ($R^{\mathrm{2}}>\mathrm{0.7}$) and ecosystem respiration ($R_{\mathrm{eco}}$; $R^{\mathrm{2}}>\mathrm{0.6}$) (e.g. Jung et al., 2010; Tramontana et al., 2016). However, the RF model performance for monthly ${\mathrm{CH}}_{\mathrm{4}}$ emissions was comparable to net ecosystem exchange of ${\mathrm{CO}}_{\mathrm{2}}$ (NEE) ($R^{\mathrm{2}}<\mathrm{0.5}$) (e.g. Jung et al., 2010; Tramontana et al., 2016). Likely reasons for this finding include, for instance, that for other fluxes there is simply more data available from several sites spanning the globe. For example, the La Thuile synthesis dataset used by Jung et al. (2010) and Tramontana et al. (2016) consists of 965 site years of data from over 252 EC stations. Here we have data from 25 sites with ${\mathrm{CH}}_{\mathrm{4}}$ fluxes. Furthermore, the drivers (or proxies for the drivers) of, for example, GPP and energy fluxes are more easily available from remote-sensing (e.g. MODIS) and weather forecasting re-analysis datasets (e.g. WFDEI). In contrast, ${\mathrm{CH}}_{\mathrm{4}}$ emissions are more related to below-ground processes, thus drivers for these processes are more difficult to measure remotely. Also, there are temporal lags between changes in drivers (e.g. LSTn) and ${\mathrm{CH}}_{\mathrm{4}}$ fluxes in response to these changes. Consequently, training a machine-learning model such as RF on such data is difficult since the RF model assumes a instantaneous relationship between the change and response. However, one should also note that GPP or $R_{\mathrm{eco}}$ are never directly measured with the EC technique, they are always at least partly derived products (Lasslop et al., 2009; Reichstein et al., 2005). Hence, direct functional relationships between GPP and $R_{\mathrm{eco}}$ and their environmental drivers are inherently included in these flux estimates, whereas NEE and ${\mathrm{CH}}_{\mathrm{4}}$ emissions are directly measured without additional modelling. Also, both NEE and ${\mathrm{CH}}_{\mathrm{4}}$ emissions are differences between component fluxes (NEE: GPP and $R_{\mathrm{eco}}$; ${\mathrm{CH}}_{\mathrm{4}}$ flux: production and oxidation). Therefore, GPP and $R_{\mathrm{eco}}$ upscaling algorithms show better correspondence with validation data than for NEE or ${\mathrm{CH}}_{\mathrm{4}}$ emissions and the results for NEE would be the correct point of reference for the RF model performance presented here.

While the RF model performance in this study was inferior to previous upscaling studies for other fluxes when evaluated using different statistical metrics, it was still comparable to what has been shown before for several process models for ${\mathrm{CH}}_{\mathrm{4}}$ emissions (McNorton et al., 2016; Wania et al., 2010; Zürcher et al., 2013; Zhu et al., 2014; Xu et al., 2016a). For instance, McNorton et al. (2016) validated the land surface model JULES against ${\mathrm{CH}}_{\mathrm{4}}$ flux data from 13 sites and found $R^{\mathrm{2}}=\mathrm{0.10}$ between the validation data and the model. Wania et al. (2010) found on average RMSE $=\mathrm{29}$ nmol m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$ and RMSE $=\mathrm{42}$ nmol m${}^{-\mathrm{2}}$ s${}^{-\mathrm{1}}$ with and without tuning their model LPJ-WHyME against ${\mathrm{CH}}_{\mathrm{4}}$ flux data from seven sites, respectively. Zürcher et al. (2013) found the time-integrated ${\mathrm{CH}}_{\mathrm{4}}$ flux to be well represented by LPX-Bern model across different sites. A tight correlation ($R^{\mathrm{2}}=\mathrm{0.92}$) is found between simulated and measured cumulative site emissions after calibrating the model against the measurements. While Xu et al. (2016a) did not explicitly show any statistical metrics, their model (CLM4.5) comparison against site level ${\mathrm{CH}}_{\mathrm{4}}$ flux data seemed to be somewhat better than in Wania et al. (2010) or McNorton et al. (2016). Xu et al. (2016a) emphasize the importance of non-growing season emissions and the fact that their model was clearly underestimating these emissions. Zhu et al. (2014) calibrated their model (TRIPLEX-GHG) for each measurement site by changing, e.g. the $Q_{\mathrm{10}}$ for ${\mathrm{CH}}_{\mathrm{4}}$ production and ${\mathrm{CH}}_{\mathrm{4}}$ to ${\mathrm{CO}}_{\mathrm{2}}$ release ratio to be site-specific and found on average $R^{\mathrm{2}}=\mathrm{0.64}$ when comparing the calibrated model against measurements at 17 ${\mathrm{CH}}_{\mathrm{4}}$ flux measurement sites. However, their findings are not directly comparable to the RF model agreement with validation data shown here due to their model calibration against data before comparison. Nevertheless, their results show that, even after calibration, the process models are not fully able to capture the ${\mathrm{CH}}_{\mathrm{4}}$ flux variability in measurements. Miller et al. (2014) argued that the structure of some of the process models is so complex that the required forcing variables may not be reliable at larger spatial scales. All of these five models (JULES, LPJ-WHyME, LPX-Bern, CLM4.5 and TRIPLEX-GHG) are contributing to the global ${\mathrm{CH}}_{\mathrm{4}}$ budget estimation within the Global Methane Project (Saunois et al., 2016), highlighting that these results summarize the agreement between state-of-the-art process models and field measurements.

4.2 Methods to improve RF model predictive performance

4.2.1 Missing predictors

In this study a statistical model was developed using the RF algorithm, and the model was able to yield $R^{\mathrm{2}}=\mathrm{0.47}$ against monthly ${\mathrm{CH}}_{\mathrm{4}}$ flux validation data. Our upscaling using the RF model focused on 2013–2014, as these were the years with the largest overlap of collected data. However, all data from all the years (2005–2016) were used to develop and validate the model. The incomplete match between the RF model and validation data is likely caused by the fact that not all the possible drivers causing inter- and intra-site variability in ${\mathrm{CH}}_{\mathrm{4}}$ emissions were included in the analysis, and hence all the variability could not be explained by the model.

Christensen et al. (2003) were able to explain practically all the variability ($R^{\mathrm{2}}=\mathrm{0.92}$) in annual ${\mathrm{CH}}_{\mathrm{4}}$ emissions in their multi-site chamber study with only two predictors: temperature and the availability of substrates for ${\mathrm{CH}}_{\mathrm{4}}$ production. Also, Yvon-Durocher et al. (2014) speculate that the amount of substrates for microbial ${\mathrm{CH}}_{\mathrm{4}}$ production explains across-site variability of ${\mathrm{CH}}_{\mathrm{4}}$ fluxes in their data. However, gridded data on spatially explicit substrate information are currently nonexistent. Hence, proxies for the substrates available for methanogenesis are needed. The current paradigm on wetland ${\mathrm{CH}}_{\mathrm{4}}$ emissions is that most of the emitted ${\mathrm{CH}}_{\mathrm{4}}$ is produced from recently fixed carbon being used as precursors for the ${\mathrm{CH}}_{\mathrm{4}}$-producing Archaea (e.g. Chanton et al., 1995; Whiting and Chanton, 1993). Most process models are based on the premise that a certain fraction of ecosystem net primary productivity (NPP) is available and used for ${\mathrm{CH}}_{\mathrm{4}}$ production or alternatively a fraction of heterotrophic respiration is allocated to ${\mathrm{CH}}_{\mathrm{4}}$ emissions (e.g. Xu et al., 2016b). Thus, NPP (or GPP) could potentially be included as a predictor for the RF model and used as a proxy for the amount of substrates available for ${\mathrm{CH}}_{\mathrm{4}}$ production. However, the RF model performance in this study was not enhanced if variables closely related to NPP (EVI and the product of EVI and LSTd) were included as predictors. Also, Knox et al. (2019) did not find GPP as an important predictor of ${\mathrm{CH}}_{\mathrm{4}}$ emission variability in their multi-site synthesis study.

Using NPP (or proxies for it) for the RF model development might be an oversimplification, since it has been shown that the deep-rooted sedges and their NPP are especially important for ${\mathrm{CH}}_{\mathrm{4}}$ production (Joabsson and Christensen, 2002; Ström et al., 2003, 2012; Waddington et al., 1996). Hence, information about plant functional types (PFTs) would be needed to better explain the ${\mathrm{CH}}_{\mathrm{4}}$ flux variability (Davidson et al., 2017; Gray et al., 2013). Furthermore, the fraction of the fixed carbon allocated to the roots and released as root exudates (hence, available for ${\mathrm{CH}}_{\mathrm{4}}$ production) varies between species and root age (Proctor and He, 2017; Ström et al., 2003), further complicating the connection between NPP and ${\mathrm{CH}}_{\mathrm{4}}$ emissions. The sedges also act as conduits for ${\mathrm{CH}}_{\mathrm{4}}$, allowing the ${\mathrm{CH}}_{\mathrm{4}}$ produced below water level to rapidly escape to the atmosphere and bypass the oxic zone in which the ${\mathrm{CH}}_{\mathrm{4}}$ might have otherwise been oxidized (Waddington et al., 1996; Whiting and Chanton, 1992). Besides sedges, Spaghnum mosses are also important because methanotrophic bacteria that live in symbiosis with these mosses significantly decrease the ${\mathrm{CH}}_{\mathrm{4}}$ emissions to the atmosphere when they are present (Larmola et al., 2010; Liebner et al., 2011; Parmentier et al., 2011b; Raghoebarsing et al., 2005; Sundh et al., 1995). In a modelling study, Li et al. (2016) showed that it was essential to consider the vegetation differences between sites when modelling ${\mathrm{CH}}_{\mathrm{4}}$ emissions from two northern peatlands. Hence, ideally one should have gridded information on wetland species composition and associated NPP across the high latitudes to significantly improve the upscaling results. Unfortunately, such information is not yet available and therefore modelled estimates could be used (e.g. LPX-Bern, which includes several peatland-specific PFTs allowed to freely evolve during the model run) (Spahni et al., 2013). However, in such cases the upscaled ${\mathrm{CH}}_{\mathrm{4}}$ emission estimates would no longer be independent of the model and therefore would be less suitable for model validation. We also note that many process models have only one PFT per wetland.

Different variables related to water input to the ecosystem (i.e. $P$, $P_{\mathrm{ann}}$) or surface moisture (SRWI) did not enhance the RF model predictive performance, not only reflecting that water table depth (WTD) is not solely controlled by input of water via precipitation but also that evapotranspiration and lateral flows affect wetland WTD, data that were missing from our study. These findings are consistent with previous studies (e.g. Christensen et al., 2003; Rinne et al., 2018; Pugh et al., 2018 and Knox et al., 2019), who showed only a modest ${\mathrm{CH}}_{\mathrm{4}}$ flux dependence on WTD in wetlands and peatlands. In contrast, several chamber-based studies have shown a positive relationship between WTD and ${\mathrm{CH}}_{\mathrm{4}}$ fluxes (Granberg et al., 1997; Olefeldt et al., 2012; Treat et al., 2018; Turetsky et al., 2014). In general, chamber-based studies often show spatial dependency of ${\mathrm{CH}}_{\mathrm{4}}$ flux on WTD, whereas studies done at an ecosystem scale with EC generally do not show temporal WTD dependency, albeit there are exceptions (e.g. Zona et al., 2009). This might indicate that WTD controls metre-scale spatial heterogeneity of ${\mathrm{CH}}_{\mathrm{4}}$ flux between microtopographical features (e.g. Granberg et al., 1997) but not temporal variability on the ecosystem scale, provided that WTD stays relative close to the surface. Also, the chamber studies tend to observe spatial variation, which can be indirectly influenced by WTD via its influence on plant communities, whereas EC studies observe typically temporal variation in sub-annual timescales. However, the effect of WTD might be masked by a confounding effect caused by plant phenology, since vegetation biomass often peaks at the same time as the WTD is at its lowest. While the variables related to WTD did not increase the RF model performance, WTD might still play a role in controlling ecosystem-scale ${\mathrm{CH}}_{\mathrm{4}}$ variability when it is exceptionally high or low. For instance, the year 2006 was exceptionally dry at the Siikaneva fen, and hence ${\mathrm{CH}}_{\mathrm{4}}$ emissions during that year were lower than on average (see Fig. 5a). However, in order to accurately capture such dependencies with machine-learning techniques (such as RF), they should be frequent enough so that the model can learn these dependencies.

RF model performance was better at permafrost-free than at sites with permafrost, which might indicate that the LSTn might not be an appropriate proxy for the temperature controlling the ${\mathrm{CH}}_{\mathrm{4}}$ production and oxidation rates at sites with permafrost. Also, no information on the development of the seasonally unfrozen, hydrologically and biogeochemically active layer was included in the RF model. Furthermore, Zona et al. (2016) showed strong hysteresis between soil temperatures and ${\mathrm{CH}}_{\mathrm{4}}$ emissions at their permafrost sites in Alaska, whereas Rinne et al. (2018) show a synchronous exponential dependence between soil temperature and ${\mathrm{CH}}_{\mathrm{4}}$ emissions at a boreal fen without permafrost. The hysteresis observed in Zona et al. (2016) could be explained by the fact that part of the produced ${\mathrm{CH}}_{\mathrm{4}}$ at these permafrost sites is stored below ground for several months before it is being emitted to the atmosphere, causing a temporal lag between soil temperature and observed surface flux. In any case, more knowledge of soil processes (soil thawing and freezing, ${\mathrm{CH}}_{\mathrm{4}}$ production and storage) is needed before the ${\mathrm{CH}}_{\mathrm{4}}$ emissions from these permafrost ecosystems can be extrapolated to other areas with greater confidence.

It should be emphasized that the drivers causing across-site variability in ecosystem-scale ${\mathrm{CH}}_{\mathrm{4}}$ emissions are, in general, unknown since studies comparing EC ${\mathrm{CH}}_{\mathrm{4}}$ fluxes from multiple wetland sites have only recently been published (Baldocchi, 2014; Knox et al., 2019; Petrescu et al., 2015). Most previous ${\mathrm{CH}}_{\mathrm{4}}$ synthesis studies were based on plot-scale measurements (Bartlett and Harriss, 1993; Olefeldt et al., 2012; Treat et al., 2018; Turetsky et al., 2014). However, the ${\mathrm{CH}}_{\mathrm{4}}$ flux responses to environmental drivers and their relative importance might be different at an ecosystem scale since ${\mathrm{CH}}_{\mathrm{4}}$ fluxes typically show significant spatial variability at sub-metre scale (e.g. Sachs et al., 2010). Furthermore, the temporal coverage of plot-scale measurements with chambers is usually relatively poor, whereas EC measurements provide continuous data on ecosystem scale. This study and Knox et al. (2019) show that temperature is important when predicting ${\mathrm{CH}}_{\mathrm{4}}$ flux variability in a multi-site ${\mathrm{CH}}_{\mathrm{4}}$ flux dataset, but a significant fraction of ${\mathrm{CH}}_{\mathrm{4}}$ flux variability is still left unexplained. It remains a challenge for future EC ${\mathrm{CH}}_{\mathrm{4}}$ flux synthesis studies to discover the drivers explaining the rest of the variability.

4.2.2

Quality and representativeness of ${\mathrm{CH}}_{\mathrm{4}}$ flux data

The RF model performance may improve if instrumentation, measurement setup and the data processing are harmonized across sites, since these discrepancies between flux sites might have caused spurious differences in ${\mathrm{CH}}_{\mathrm{4}}$ fluxes. These differences would have created additional variability in the synthesis dataset, which would in turn (1) influence the training of RF model and (2) decrease, for example, NSE values obtained against validation data, since there would be artificial variability in the validation data, which is not related to the predictors. In this study, the site PIs processed the data themselves using different processing codes, albeit the gap filling was done centrally in a standardized way.

While these issues mentioned above could impact the upscaling results shown here, prior studies have shown that the usage of different instruments or processing codes does not significantly impact ${\mathrm{CH}}_{\mathrm{4}}$ flux estimates. For instance, Mammarella et al. (2016) showed that the usage of different processing codes (EddyPro and EddyUH) resulted, in general, in a 1 % difference in long-term ${\mathrm{CH}}_{\mathrm{4}}$ emissions. On the other hand, ${\mathrm{CH}}_{\mathrm{4}}$ instrument cross comparisons have shown small differences (typically less than 7 %) between the long-term ${\mathrm{CH}}_{\mathrm{4}}$ emission estimates derived using different instruments (Goodrich et al., 2016; Peltola et al., 2013, 2014). While these studies show consistent ${\mathrm{CH}}_{\mathrm{4}}$ emissions, they also stress that the data should be carefully processed to achieve such good agreement across processing codes and instruments. In addition, many issues related to, for example, friction velocity filtering and gap filling of ${\mathrm{CH}}_{\mathrm{4}}$ fluxes are still unresolved, and the role of short-term emission bursts, which are common in methane flux time series, needs to be further investigated (e.g. Schaller et al., 2017). However, recently Nemitz et al. (2018) advanced these issues by proposing a methodological protocol for EC measurements of ${\mathrm{CH}}_{\mathrm{4}}$ fluxes used to standardize ${\mathrm{CH}}_{\mathrm{4}}$ flux measurements within the ICOS measurement network (Franz et al., 2018).

A total of 25 flux measurement sites were included in this study and they were distributed across the Arctic boreal region (see Fig. 1). The measurements were largely concentrated in Fennoscandia and Alaska, whereas data from, for example, the HBL and WSL areas, were missing. Long-term EC ${\mathrm{CH}}_{\mathrm{4}}$ flux measurements are largely missing from these vast wetland areas, casting uncertainty on wetland ${\mathrm{CH}}_{\mathrm{4}}$ emissions from these areas. The location of a flux site is typically restricted by practical limitations related to, for example, ease of access and availability of grid power. Hence, open-path instruments with low power requirements potentially open up new areas for flux measurements (McDermitt et al., 2010), yet they need continuous maintenance, which is not necessarily easy in remote locations. However, one could argue that the geographical location of flux sites is not vital for upscaling, more important is that the available data represents well the full range of ${\mathrm{CH}}_{\mathrm{4}}$ fluxes across the northern latitudes and more importantly the ${\mathrm{CH}}_{\mathrm{4}}$ flux responses to the environmental drivers. Also, sites should ideally cover all different wetlands with varying plant species composition, whereas geographical representation is not necessarily as important. ${\mathrm{CH}}_{\mathrm{4}}$ flux site representativeness could be potentially assessed in the same vein as in previous studies for other measurement networks (Hargrove et al., 2003; Hoffman et al., 2013; Papale et al., 2015; Sulkava et al., 2011). However, before such analysis can be done, the main drivers causing across-site variability in ecosystem-scale ${\mathrm{CH}}_{\mathrm{4}}$ fluxes should be better identified.

Most of the ${\mathrm{CH}}_{\mathrm{4}}$ flux data here and in the literature have been recorded during the growing season when the ${\mathrm{CH}}_{\mathrm{4}}$ fluxes are at a maximum, whereas year-round continuous ${\mathrm{CH}}_{\mathrm{4}}$ flux measurements are not as common. This is likely due to the harsh conditions in the Arctic during winter that make continuous high-quality flux measurements very demanding (e.g. Goodrich et al., 2016; Kittler et al., 2017a) but also in part since the large-scale importance of non-growing season emissions has just recently been recognized (Kittler et al., 2017b; Treat et al., 2018; Xu et al., 2016a; Zona et al., 2016). For upscaling year-round ${\mathrm{CH}}_{\mathrm{4}}$ emissions, continuous measurements are vital to accurately constrain also the non-growing season emissions and their drivers.

5 Data availability

The presented upscaled ${\mathrm{CH}}_{\mathrm{4}}$ flux maps (RF-DYPTOP, RF-PEATMAP and RF-GLWD), their uncertainties and the underlying ${\mathrm{CH}}_{\mathrm{4}}$ flux densities are accessible via an open-data repository Zenodo (Peltola et al., 2019). The datasets are saved in netCDF files and they are accompanied by a readme file. The dataset can be downloaded from 10.5281/zenodo.2560163.

6 Conclusions

Methane (${\mathrm{CH}}_{\mathrm{4}}$) emission data comprising over 40 site years from 25 eddy covariance flux measurement sites across the Arctic boreal region were assembled and upscaled to estimate ${\mathrm{CH}}_{\mathrm{4}}$ emissions from northern ($>\mathrm{45}$${}^{\circ }$ N) wetlands. The upscaling was done using the random forest (RF) algorithm. The performance of the RF model was evaluated against independent validation data utilizing the leave-one-site-out scheme, which yielded value of 0.47 for both the Nash–Sutcliffe model efficiency and $R^{\mathrm{2}}$. These results are similar to previous upscaling studies for the net ecosystem exchange of carbon dioxide (NEE) but worse for the individual components of NEE or energy fluxes (e.g. Jung et al., 2010; Tramontana et al., 2016). The performance is also comparable to studies where process models are compared against site ${\mathrm{CH}}_{\mathrm{4}}$ flux measurements (McNorton et al., 2016; Wania et al., 2010; Zürcher et al., 2013; Zhu et al., 2014; Xu et al., 2016a). Hence, despite the relatively high fraction of unexplained variability in the ${\mathrm{CH}}_{\mathrm{4}}$ flux data, the upscaling results are useful for comparing against models and could be used to evaluate model results. The three gridded ${\mathrm{CH}}_{\mathrm{4}}$ wetland flux estimates and their uncertainties are openly available for further usage (Peltola et al., 2019).

The upscaling to the regions $>\mathrm{45}$${}^{\circ }$ N resulted in mean annual ${\mathrm{CH}}_{\mathrm{4}}$ emissions comparable to prior studies on wetland ${\mathrm{CH}}_{\mathrm{4}}$ emissions from these areas (Bruhwiler et al., 2014; Chen et al., 2015; Spahni et al., 2011; Treat et al., 2018; Watts et al., 2014; Zhang et al., 2016; Zhu et al., 2013) and hence, in general, support the prior modelling results for the northern wetland ${\mathrm{CH}}_{\mathrm{4}}$ emissions. When compared to two validation areas, the upscaling likely overestimated ${\mathrm{CH}}_{\mathrm{4}}$ emissions from the Hudson Bay lowlands, whereas emission estimates for the western Siberian lowlands were in a reasonable range. Future ${\mathrm{CH}}_{\mathrm{4}}$ flux upscaling studies would benefit from long-term continuous ${\mathrm{CH}}_{\mathrm{4}}$ flux measurements, centralized data processing and better incorporation of ${\mathrm{CH}}_{\mathrm{4}}$ flux drivers (e.g. wetland vegetation composition and carbon cycle) from remote-sensing data needed for scaling the fluxes from the site level to the whole Arctic boreal region.

Appendix A Table A1

Description of eddy covariance sites included in this study.

Name	Site ID	PI	Latitude,	Amount of	Reference	Permafrost	Sedges as	Biome	Wetland	Time
			longitude	monthly ${\mathrm{CH}}_{\mathrm{4}}$		present	dominant	based on	type	resolution
				flux data		(true/false)	vegetation	Olson et		of data
				available			type	al. (2011)
Schechenfilz	DE-SfN	Janina Klatt,	47.8064,	2	Hommeltenberg	false	false	temperate	bog	30 min
Nord		Hans Peter Schmid	11.3275		et al. (2014)
Chokurdakh	RU-Cok	Albertus J. Dolman	70.8291,	5	Parmentier et	true	true	tundra	wet	30 min
			147.4943		al. (2011a)				tundra
Vorkuta	RU-Vor	Thomas Friborg	67.0547,	5	Marushchak et	true	false	tundra	wet	30 min
			62.9405		al. (2016)				tundra
Stordalen	SE-St1	Thomas Friborg	68.3542,	6	Jammet et	false	true	tundra	fen	30 min
			19.0503		al. (2017)
Stordalen	SE-Sto	Janne Rinne	68.3560,	55		true and	false	tundra	bog	30 min
(ICOS)${}^{*}$			19.0452			false
Siikaneva 1	FI-Sii	Timo Vesala,	61.8327,	104	Rinne et	false	true	boreal	fen	30 min
		Ivan Mammarella	24.1928		al. (2018)
Siikaneva 2	FI-Si2	Timo Vesala,	61.8375,	26	Korrensalo et	false	false	boreal	bog	30 min
		Ivan Mammarella	24.1699		al. (2018)
Lompolojänkkä	FI-Lom	Annalea Lohila	67.9972,	59	Aurela et	false	true	boreal	fen	30 min
			24.2092		al. (2009)
James Bay	CA-JBL	Daniel F. Nadeau	53.6744,	3	Nadeau et	false	false	boreal	bog	daily
lowlands			$-\mathrm{78.1706}$		al. (2013)
Lost Creek	US-Los	Ankur R. Desai	46.0827,	30	Pugh et	false	false	temperate	fen	30 min
			$-\mathrm{89.9792}$		al. (2018)
Atqasuk	US-Atq	Donatella Zona	70.4696,	11	Zona et	true	true	tundra	wet	30 min
			$-\mathrm{157.4089}$		al. (2016)				tundra
Barrow	US-Beo	Donatella Zona	71.2810,	16	Zona et	true	true	tundra	wet	30 min
Environmental			$-\mathrm{156.6123}$		al. (2016)				tundra
Observatory
Biocomplexity	US-Bes	Donatella Zona	71.2809,	16	Zona et	true	true	tundra	wet	30 min
Experiment			$-\mathrm{156.5965}$		al. (2016)				tundra
South tower
Ivotuk	US-Ivo	Donatella Zona	68.4865,	15	Zona et	true	true	tundra	wet	30 min
			$-\mathrm{155.7502}$		al. (2016)				tundra
Western	CA-WP1	Lawrence B.	54.9538,	5	Long et	false	false	temperate	fen	30 min
peatland 1		Flanagan	$-\mathrm{112.4670}$		al. (2010)
Mer Bleue	CA-Mer	Elyn Humphreys	45.4094,	16	Brown et	false	false	temperate	bog	daily
			$-\mathrm{75.5186}$		al. (2014)
Chersky	RU-Ch2	Mathias Göckede	68.6169,	21	Kittler et	true	true	boreal	wet	daily
reference			161.3509		al. (2017)				tundra
Rzecin	PL-wet	Janusz Olejnik	52.7622,	4	Kowalska et	false	true	temperate	fen	30 min
			16.3094		al. (2013)
Degerö	SE-Deg	Mats B. Nilsson,	64.1820,	22	Nilsson et	false	true	boreal	fen	30 min
Stormyr		Matthias Peichl	19.5567		al. (2008)
Seney	US-Sen	Thomas Pypker	46.3167,	5	Pypker et	false	true	temperate	fen	daily
			$-\mathrm{86.0500}$		al. (2013)
Scotty Creek	CA-SCC	Oliver Sonnentag	61.3000,	14	Helbig et	false	false	boreal	bog	30 min
			$-\mathrm{121.300}$		al. (2016)
Samoylov	RU-Sam	Torsten Sachs	72.3667,	11	Sachs et	true	true	tundra	wet	30 min
			126.5000		al. (2008)				tundra
Imnavait	US-ICh	Eugenie S.	68.6060,	7		true	true	tundra	wet	30 min
Creek		Euskirchen	$-\mathrm{149.3110}$						tundra
Bonanza	US-BCF	Eugenie S.	64.7040,	16	Euskirchen et	false	true	boreal	fen	30 min
Creek, fen		Euskirchen	$-\mathrm{148.3130}$		al. (2014)
Bonanza	US-BCB	Eugenie S.	64.7000,	14	Euskirchen et	false	false	boreal	bog	30 min
Creek, bog		Euskirchen	$-\mathrm{148.3200}$		al. (2014)

${}^{*}$ Data from this site is divided into two since data from two wind directions differ from each other (with and without permafrost).

Author contributions

OP, TA and TV designed the study and YG contributed further ideas for the study. OP did the data processing and analysis. OR prepared the PEATMAP map for the study. PA, MA, BC, ARD, AJD, ESE, TF, MG, MH, EH, GJ, JK, NK, LK, AL, IM, DFN, MBN, WCO, MP, TP, WQ, JR, TS, MS, HPS, OS, CW and DZ provided ${\mathrm{CH}}_{\mathrm{4}}$ fluxes and other in situ data for the study. FJ and SL did the LPX-Bern model runs. OP wrote the first version of the manuscript and all authors provided input.

Competing interests

The authors declare that they have no conflict of interest.

Acknowledgements

Lawrence B. Flanagan is acknowledged for providing data from CA-WP1 site. Lawrence B. Flanagan acknowledges support from the Natural Sciences and Engineering Council of Canada and Canadian Foundation for Climate and Atmospheric Sciences. Olli Peltola is supported by the postdoctoral researcher project (decision 315424) funded by the Academy of Finland. Olle Räty is supported by the Academy of Finland IIDA-MARI project (decision 313828). Financial support from the Academy of Finland Centre of Excellence (grant nos. 272041 and 307331), Academy Professor projects (grant nos. 312571 and 282842), ICOS-Finland (grant no. 281255) and CARB-ARC project (grant no. 285630) is acknowledged. Sara H. Knox and Robert B. Jackson acknowledge support from the Gordon and Betty Moore Foundation through grant GBMF5439 “Advancing Understanding of the Global Methane Cycle”. Ankur R. Desai acknowledges support of the DOE Ameriflux Network Management Project. Albertus J. Dolman acknowledges support from the Netherlands Earth System Science Centre, NESSC). Torsten Sachs was supported by the Helmholtz Association of German Research Centres (grant no. VH-NG-821). Ivan Mammarella and Timo Vesala thank the EU for supporting the RINGO project funded by the Horizon 2020 Research and Innovation Programme (grant no. 730944). The EU-H2020 CRESCENDO project (grant no. 641816) is also acknowledged. Fortunat Joos and Sebastian Lienert are thankful for support from the Swiss National Science Foundation (grant no. 200020_172476). Mats B. Nilsson and Matthias Peichl acknowledge support from the National Research Council (VR 2018-03966) SITES and ICOS-Sweden.

Financial support

This research has been supported by the Academy of Finland (grant nos. 315424, 313828, 312571, 282842, 281255, 285630, 272041 and 307331), the Gordon and Betty Moore Foundation (grant no. GBMF5439), the Helmholtz Association (grant no. VH-NG-821), and Horizon 2020 (RINGO (grant no. 730944) and CRESCENDO (grant no. 641816)).

Review statement

This paper was edited by David Carlson and reviewed by two anonymous referees.