Predicting catchment-scale methane fluxes with multi-source remote sensing

Spatial patterns of CH4 fluxes can be modeled with remotely sensed data representing land cover, soil moisture and topography. Spatially extensive CH4 flux measurements conducted with portable analyzers have not been previously upscaled with remote sensing. How well can the CH4 fluxes be predicted with plot-based vegetation measures and remote sensing? How does the predictive skill of the model change when using different combinations of predictor variables? We measured CH4 fluxes in 279 plots in a 12.4 km2 peatland-forest-mosaic landscape in Pallas area, northern Finland in July 2019. We compared 20 different CH4 flux maps produced with vegetation field data and remote sensing data including Sentinel-1, Sentinel-2 and digital terrain model (DTM). The landscape acted as a net source of CH4 (253–502 µg m−2 h−1) and the proportion of source areas varied considerably between maps (12–50%). The amount of explained variance was high in CH4 regressions (59–76%, nRMSE 8–10%). Regressions including remote sensing predictors had better performance than regressions with plot-based vegetation predictors. The most important remote sensing predictors included VH-polarized Sentinel-1 features together with topographic wetness index and other DTM features. Spatial patterns were most accurately predicted when the landscape was divided into sinks and sources with remote sensing-based classifications, and the fluxes were modeled for sinks and sources separately. CH4 fluxes can be predicted accurately with multi-source remote sensing in northern boreal peatland landscapes. High spatial resolution remote sensing-based maps constrain uncertainties related to CH4 fluxes and their spatial patterns.


Introduction
Methane (CH 4 ) is the second most important greenhouse gas (IPCC 2013), and a large part of CH 4 emissions originate from northern peatlands (Frolking et al. 2011;Turetsky et al. 2014;Abdalla et al. 2016). On the one hand, CH 4 is produced by methanogens in anoxic conditions, which are typically obtained in water-saturated conditions in which diffusion of gases is low and oxygen is depleted quickly. CH 4 is transported from soil to the atmosphere by molecular diffusion, in the form of gas bubbles (i.e. ebullition) and plant-mediated processes. On the other hand, CH 4 can be oxidized to carbon dioxide by CH 4 consuming methanotrophs in oxic conditions, typically found in the surface peat layers. A net emission into the atmosphere is observed when production and transport exceed consumption (Le Mer and Roger 2001;Lai 2009;Abdalla et al. 2016).
Important factors controlling CH 4 emissions are soil temperature, water table depth, plant community composition, and soil pH (Turetsky et al. 2014;Abdalla et al. 2016). These factors are spatially heterogenic in diverse scales in the landscape, and there are considerable uncertainties in the resulting net emission of CH 4 by the peatlands and how the emissions are distributed spatially (Turetsky et al. 2014). In some circumstances, also mineral soils can act as CH 4 sources (Lohila et al. 2016). To constrain the different uncertainties, there have been attempts to assess and upscale CH 4 flux patterns over different spatial extents ranging from a single mire (Lehmann et al. 2016) to landscapes (Dinsmore et al. 2017) and also to synthesize CH 4 controls globally (Turetsky et al. 2014).
The spatial patterns of CH 4 fluxes can be assessed with remote sensing-based maps. One approach is to monitor CH 4 fluxes with imaging spectroscopy, in particular in shortwave and thermal infrared region. However, this approach has been limited to CH 4 plumes instead of diffuse fluxes originating from peatlands (Thorpe et al. 2016). Coarse-grained estimates of CH 4 fluxes can be derived with large scale process models resembling carbon cycle (Bruhwiler et al. 2014) or more locally with the eddy covariance method (Davidson et al. 2017). However, in spatially heterogeneous landscapes, such as forest-peatland mosaics, CH 4 fluxes have fine-scale spatial heterogeneity. Therefore, a typical approach is to measure fluxes with chambers and then upscale of the observations with remote sensing data (Lehmann et al. 2016;Dinsmore et al. 2017;Morozumi et al. 2019).
Chambers can measure fluxes only from a small area (\ 0.5 m 2 ), and the number of chamber plots is limited in one study; therefore, there remain large uncertainties whether the limited number of chamber plots are representatively sampled over a specific landscape (Lehmann et al. 2016;Davidson et al. 2017). A typical approach is to upscale fluxes with categorical land cover maps (Lehmann et al. 2016;Morozumi et al. 2019), and it has been shown that more accurate upscaling can be derived with land cover maps with higher thematic resolution (Lehmann et al. 2016). However, it has been discussed that maps of continuous variables (e.g. topographic features, spectral indices and tree canopy cover) represent more accurately landscape processes such as CH 4 fluxes, which are typically not limited to crisp patches (McGarigal et al. 2009;Kedron et al. 2018). Therefore, another option is to conduct CH 4 upscaling with continuous maps of specific property such as satellitederived spectra (Dinsmore et al. 2017).
It is possible to link CH 4 fluxes to diverse remotely sensed or other data sources. In numerous mapping endeavors it has been shown that multi-source remote sensing (i.e. inclusion of multiple types of remotely sensed information in one mapping workflow) increases the accuracy (Bourgeau-Chavez et al. 2016;Karlson et al. 2019;Räsänen and Virtanen 2019). As CH 4 fluxes in northern landscapes are mainly connected to peatlands, it can be assumed that similar methods can be used for detection of peatlands and tracking spatial patterns of CH 4 . In previous studies, it has been shown that multiple remote sensing datasets are useful in wetland and peatland delineation, containing optical data depicting spectral properties of vegetation and land cover (Bourgeau-Chavez et al. 2016;Minasny et al. 2019), synthetic aperture radar (SAR) data sensitive to moisture, surface roughness and vegetation structure (Widhalm et al. 2015;Bourgeau-Chavez et al. 2016;Millard and Richardson 2018;Karlson et al. 2019;Minasny et al. 2019) and topographic data capable to model presumably wet areas based on water flow routes and flatness (Murphy et al. 2007;Karlson et al. 2019;Minasny et al. 2019;Lidberg et al. 2020). Nonetheless, the use of multi-source remote sensing data has been limited in CH 4 flux studies. Furthermore, recent technological development of portable CH 4 measurement devices has opened up possibilities for spatially dense and extensive field sampling. However, to the best of our knowledge, spatially extensive measurements have not yet been linked to remotely sensed data.
We measured CH 4 fluxes with spatially extensive chamber measurements (n = 279) in two adjacent catchments (altogether 12.4 km 2 ) in northern Finland covering all main land cover types in the area during a two-week period in early July 2019. We then upscaled the measurements to the whole study area with multisource remote sensing data including Sentinel-2 optical satellite imagery, Sentinel-1 SAR data and topographic data. Our two research questions were: How well can the CH 4 fluxes be predicted with plotbased vegetation measures and remote sensing data representing land cover, soil moisture and topography? How does the predictive skill of the model change when using different combinations of predictor variables?

Study area
We studied landscape-level CH 4 fluxes in two adjacent catchments (4.5 km 2 and 7.9 km 2 respectively) in Pallas area, located in Muonio municipality, in northern Finland (67°57 0 -68°01 0 N; 24°10 0 -24°15 0 E, Fig. 1). The area belongs to northern boreal vegetation zone and subarctic climate zone with the annual average temperature being -1.3°C, precipitation 547 mm, and the number of growing degree days above 5°C 698 (years 1981-2010, data from nearby Kittilä Pokka station). The area is characterized by gently undulating terrain with the elevation above sea level ranging between 267 and 564 m. According to our land cover classification, coniferous and mixed coniferous-deciduous forests (68%) and forested and open peatlands (17%) dominate the whole area, while its western part includes mountain tundra and treeline vegetation (14%) (Fig. 1, Table 1; methods described in Supplementary material). Dominating tree species include Norway spruce (Picea abies), Scots pine (Pinus sylvestris) and Downy birch (Betula pubescens). The western part of the study area belongs to Pallas-Yllästunturi National Park and includes old-growth forests, mostly with mixed tree species composition (Fig. 1). The eastern part contains mostly coniferous dominated forests, as in the latter half of the twentieth century, deciduous trees were selectively logged, and Scots pines were planted in some locations. Furthermore, some of the peatlands in the eastern part have been drained for forestry purposes in the same period. Atmospheric, ecological and hydrological research within and nearby the study area has been extensive during the past decades (Lohila et al. 2015).

Field inventory and vegetation data
We gathered the vegetation cover data and conducted the CH 4 flux measurements between July 3 and 13, 2019 in twelve 0.5-2.0 km long transects located within the catchments so that all major land cover types were included (Fig. 1). Along each transect, we sampled circular plots with 10 m diameter in 100 m intervals. Within most of the circular plots, we sampled three subplots with 0.5 m diameter located at a 2 m distance to the north, east, and west from the plot centroid. In total, we sampled 131 larger plots and 300 subplots. We collected the location of each plot with Trimble R8 GPS device with ± 0.05 m accuracy.
For each plot, we calculated a number of plot-based vegetation measures representing vegetation composition and structure. In the larger plots, we measured the height, crown length and diameter at breast height of each tree located within the plot. Within the subplots, we estimated %-cover of each vascular plant and moss species and measured the average height of each vascular plant species. For each larger plot, we estimated the ground and field layer biomass and leafarea index for four vascular plant functional types (evergreen shrubs, deciduous shrubs, forbs, and graminoids) and biomass for mosses as well as total biomass and leaf-area index (Supplementary material). Additionally, we estimated tree biomass and leaf-area index for spruce, pine, coniferous trees in total, deciduous trees in total, and all trees in total (Supplementary material). For each subplot, we defined membership values for fuzzy plant community clusters (Supplementary material).
At each of the 279 subplots, we conducted chamber-based measurements of CH 4 fluxes once during the field inventory period. The chamber closure time was always 5 min and measurements took place between 9:15 AM and 5:30 PM. We used portable LI-COR trace gas analyzer (LI-7810, LI-COR Environmental, Lincoln, NE, USA) and a closed opaque chamber (height 30.5 cm, diameter 31.5 cm). During the flux measurements, we measured the air temperature inside the chamber headspace. We measured soil Fig. 1 Land cover map of the study area located in northern Finland. The classification accuracy was 76% and methods description is given in Supplementary material Non-vegetated Non-vegetated areas, such as rock, sand and gravel roads 0.7 0 temperature and moisture next to the chamber at the surface (WET-2 sensor, Delta-T Devices Ltd, Cambridge, UK). Additionally, we recorded 5 cm soil temperature (Pt-100 RTD Thermometer HH376, Omega, Taiwan). CH 4 flux (F) was calculated according to ideal gas law: where dC t ð Þ dt t¼0 is the time derivative (ppm s -1 ) of exponential regression (e.g. Korkiakoski et al. 2017) at the beginning of the closure, M is the molecular mass of CH 4 (16.042 g mol -1 ), P is the air pressure (Pa) measured from the nearby official weather station (Kenttärova, \ 2 km from the flux measurements), R is the universal gas constant (8.314 J mol -1 K -1 ), T is the chamber headspace temperature (K), V and A are the chamber volume (m 3 ) and base area (m 2 ), respectively. Here, a micrometeorological sign convention is used where a positive flux indicates a flux from the ecosystem to the atmosphere (CH 4 emission/source) and a negative flux indicates a flux from the atmosphere to the ecosystem (CH 4 uptake/sink). To control for diurnal variation and changing weather conditions, we standardized the CH 4 fluxes to 10°C soil temperature. We measured fluxes in 84 (30%) of the subplots 9-10 times between June and August and in the rest only once between July 3 and 13, 2019. For the subplots with repeated measurements, we estimated F CH4;10c by using Arrhenius model: where F CH4;meas is the measured CH 4 flux, T is the soil temperature (K) at 5 cm, F CH4;10c is the CH 4 flux at 10°C and Q 10 is the coefficient describing flux sensitivity to temperature. We kept Q 10 and F CH4;10c as free parameters during fitting. For the subplots with no repeated measurements, we solved F CH4;10c directly from Eq. (2) by using median Q 10 estimated from repeated measurements (2.149 for F CH4;meas \ 0 and 2.386 for F CH4;meas [ 0) and measured T and F CH4;meas .

Remote sensing data
We acquired diverse remote sensing datasets to map land cover, moisture and topography patterns within the study area (Table 2). For predicting CH 4 flux patterns, we were primarily interested in analyzing the applicability of optical Sentinel-2 satellite images and Sentinel-1 SAR data (European Space Agency, ESA) as well as a lidar-based digital terrain model (DTM, National Land Survey of Finland, NLS). In addition, we included canopy and intensity metrics from aerial lidar data, spectral and textural features from an aerial orthophoto (NLS), and spectral indices from PlanetScope satellite images (Planet Labs Inc., San Francisco, CA, USA) in the analysis. All optical remote sensing data (Sentinel-2, orthophoto, PlanetScope) were cloud-free. For each field inventory plot, we calculated features from each dataset. Sentinel-2 satellites carry optical instruments detecting reflected solar radiation from the Earth's surface at 13 different wavelength bands (Drusch et al. 2012). Optical data can be used to differentiate vegetation and land cover types, due to diverging reflectance properties of landscape elements. We used the bottom-of-atmosphere reflectance Sentinel-2 image (Level-2A; atmospherically corrected with ESA's Sen2Cor algorithm) taken during our fieldwork (July 8). We included 11 bands (excluding bands 8a and 10, i.e. narrow near-infrared and cirrus respectively) in the analysis and also calculated 29 different widely used spectral indices (Table S1) with RStoolbox (Leutner et al. 2019) in R. These included indices that capture vegetation greenness and plant community composition such as normalized difference vegetation index (Rouse et al. 1974;Pettorelli et al. 2005;McPartland et al. 2019), and soil moisture such as water indices (Gao 1996;McFeeters 1996;Arroyo-Mora et al. 2018).
Sentinel-1 satellites carry C-band SAR instruments emitting microwave pulses at a central frequency of 5.405 GHz and measuring their backscatter characteristics (Torres et al. 2012). The backscatter is sensitive to land surface roughness, vegetation structure and in particular moisture and water (Widhalm et al. 2015;Bourgeau-Chavez et al. 2016;Millard and Richardson 2018;Karlson et al. 2019;Minasny et al. 2019). Backscatter is also dependent on incidence angle, and it thus has been recommended to use data with multiple incidence angles (O'Grady et al. 2014;Millard and Richardson 2018). We acquired ground range Sentinel-1 data from three consecutive days (July 7-July 9), each taken in the morning but with different incidence angles (mean values 34.2°, 39.4°a nd 44.1°respectively). We corrected the datasets geometrically with the help of the DTM and nearest neighborhood resampling and resampled the data to 10 m spatial resolution. We used both VH (vertical European Space Agency 10 m 30: VH and VV polarization GRD images from three consecutive days, generalized using mean filter with five different window sizes (3 9 3, 5 9 5, 7 9 7, 9 9 9, 11 9 11) Sentinel-2 indices are specified in Table S1 B blue, DTW depth to water, G green, GLCM grey-level co-occurrence matrix, GRD ground range detected, GRVI green-red vegetation index, NDVI normalized difference vegetation index, NDWI normalized difference water index, NIR near infra-red, R red, SD standard deviation, TPI topographical position index, TWI topographic wetness index, VH vertical transmit, horizontal receive, VV vertical transmit, vertical receive, l mean transmit, horizontal receive) and VV (vertical transmit, vertical receive) polarized bands, and reduced speckle using mean filter with five different window sizes (3 9 3, 5 9 5, 7 9 7, 9 9 9, 11 9 11) ( Table 2). Topography is a major factor controlling soil moisture and hydrological patterns within landscapes (Murphy et al. 2007;Lidberg et al. 2020). From the topographic DTM, we calculated following layers with SAGA-GIS (Conrad et al. 2015): elevation, slope, topographic position indices with different neighborhood distances (Guisan et al. 1999), topographic wetness index (TWI, Böhner and Selige 2006) and two depth to water (DTW) indices with stream networks derived from DTM with minimum 4 ha catchment area (DTW-DTM) and NLS topographic database (DTW-NLS) respectively (Murphy et al. 2007). Of these predictors, in particular TWI and DTW have been shown to be useful in predicting spatial patterns of wetness (Murphy et al. 2007;Lidberg et al. 2020) and position index has been found to be function well in locating wetlands in depressions (Riley et al. 2017).
As forest structure differs between forests on mineral soil and peatlands ), we also calculated seven canopy metrics from the lidar point cloud for each larger plot using LAStools (Isenburg 2018) (Table 2). We also included two lidar intensity layers (all returns and ground returns only) whose functionality in wet area mapping has been documented (Stevens and Wolfe 2012).
For each orthophoto band, we calculated eight different grey-level co-occurrence matrix texture layers, since the quantification of pixel value neighborhood patterns (i.e. texture) has shown to increase distinctiveness between different land covers (Hall-Beyer 2017). The layers were energy, entropy, correlation, inverse difference moment, inertia, cluster shade, cluster prominence and Haralick correlation (Haralick et al. 1973) calculated with Orfeo ToolBox (Grizonnet et al. 2017).
Wet and dry areas have disparate seasonal cycles in their spectral properties, and it has been suggested that multi-temporal imagery should be used (Halabisky et al. 2018). Therefore, we used six PlanetScope images, each from a different phenological stage. The images were delivered as orthocorrected surface reflectance products. To correct the minor geometrical errors and match the data with the orthophoto, we geometrically corrected the images with 12 reference points by using a combination of polynomial transformation and a triangulated irregular network interpolation as well as nearest neighbor resampling. For each image, we calculated two vegetation indices capable of capturing vegetation phenological dynamics, i.e. normalized difference vegetation index (Rouse et al. 1974;McPartland et al. 2019), and green-red vegetation index (Gitelson et al. 2002;Anderson et al. 2016), and one water index, i.e. normalized difference water index (McFeeters 1996) (Table 2).

CH 4 regressions and upscaling with maps
We evaluated how much variance of measured CH 4 fluxes can be explained with different explanatory predictor sets and what kind of differences there are in the upscaled landscape-level fluxes and their spatial patterns when using 20 different CH 4 maps. We constructed the first map with the land cover classification so that we assigned the mean measured CH 4 flux estimate for each land cover type (classification accuracy 76%, Table 1, Fig. 1, Supplementary material). The following 19 maps were generated with regression-based techniques, in which predictors were remote sensing features and plot-based vegetation measures. We produced 15 maps in which we included all CH 4 measurement data (n = 279) in one regression (Table 3). Initially, we tested 14 different predictor set options, and we further tested if a smaller number of predictors would notably increase the regression performance ([ 3%-point increase in explained variance) in each predictor set option. Consequently, we included a regression with three DTM-based predictors (elevation, TWI, DTW-NLS) as it had higher performance than regression with all DTM predictors. The subsequent four maps were produced so that we first divided the CH 4 measurement data before regressions into sinks (n = 194) and sources (n = 85). We tested two options for separating sinks and sources in the landscape: separation of sinks and sources based on land cover types ( Table 1) and separation of sinks and sources with a binary classification in which training data consisted of CH 4 flux measurement plots. In the sink and source regressions, we tested two different predictor sets: all predictors and remote sensing predictors (i.e. options 1 and 3 in Table 3).
To produce CH 4 maps other than land cover-based maps, we conducted random forest regressions (Breiman 2001) and binary classifications with 500 trees. Random forest is insensitive to overfitting and can thus be used in remote sensing analyses with several predictors (Belgiu and Dragut 2016). We tested onethird of all predictors at each tree node (square root of all predictors in binary classifications). We evaluated model performance with random forest out-of-bag assessment in which two-thirds of the data was used for training and the rest for validation for each tree. This method has been reported to result in a slightly conservative estimate of model performance (Clark et al. 2010). As the result of the random forest may slightly change between different random forest runs, we conducted 100 random forest runs for each predictor set option. For each regression, we calculated mean, minimum and maximum explained variance (random forest pseudo R 2 = 1 -(mean squared error)/variance(response)), root mean square error (RMSE) and normalized RMSE (nRMSE; i.e. RMSE divided by the range of observed values). For binary classifications, we calculated overall classification accuracy.
To remove irrelevant predictors from the final regressions and binary classifications, we reduced the number of predictors with Boruta (Kursa and Rudnicki 2010), which is a random forest wrapper predictor selection algorithm. We used mean decrease accuracy variable importance measure. We chose those variables to final random forest models, which were not rejected during 999 random forest runs in Boruta (i.e. their importance was not statistically significantly lower than that of randomly shuffled shadow variables). We also calculated mean predictor importance scores over the Boruta runs, and Pearson correlation coefficients between the predictors and CH 4 flux. Analyses were conducted in R with packages ran-domForest (Liaw and Wiener 2002) and Boruta (Kursa and Rudnicki 2010).
For those remote-sensing based maps that included predictors only from one dataset and for maps combining Sentinel-1 and Sentinel-2 datasets, we produced pixel-based maps with the spatial resolution being as it was in the input data (Table 2). For the other regression-based maps, we utilized a geographic object-based image analysis approach with Baatz and Schäpe (2000) segmentation in TerraView (Câmara et al. 2008), and calculated similar predictors to each segment as for each field inventory plot (Table 2). For regressions including plot-based vegetation predictors, we utilized maps of those predictors (map production described in Supplementary material). We segmented the aerial orthophoto with the following parameter information: minimum segment size 100 (25 m 2 ), similarity threshold 0.15, weight for color 0.75, and weight for compactness 0.5. The final average segment size was 50 m 2 . Finally, to quantify differences in the spatial patterns between maps, we resampled all maps into 10 m resolution and calculated a standard error of the mean map and pairwise Pearson correlation coefficients.

Results
The measured CH 4 fluxes at the catchment during the campaign in the first half of July 2019 varied between sinks and sources (-662-14,211 lg CH 4 m -2 h -1 ). In general, peatlands were sources (on average 2573 lg CH 4 m -2 h -1 ) while forests and tundra on mineral soil were sinks (on average -154 and -133 lg CH 4 m -2 h -1 respectively) (Fig. 2, Table S11). However, some of the plots on drier peatlands (bogs and paludified forests) were weak sinks, and some of the plots on mineral soil were weak sources (Fig. 2).
Combined regressions for sinks and sources explained 59-76% of the variance in CH 4 fluxes, while sink regressions explained 21-24% and source regressions 58-59% (Table 4). In the binary sink and source classifications, mean classification accuracy was 96.5% (minimum to maximum 95.7-97.1%) and 96.6% (96.1-97.1%) in classification with all predictors and classification with remote sensing predictors, respectively.
In the regressions with all data, most of the regressions with remote sensing predictors only had higher performance than regressions with plot-based vegetation predictors. Nevertheless, the regression combining remote sensing and plot-based vegetation predictors had the highest performance (Table 4). Of the different remote sensing regressions, regression omitting DTM had the highest performance, and  regression omitting both Sentinel datasets and DTM had the worst performance. Relatively good performance was obtained in regressions with a single remote sensing data source, with Sentinel-1 having the highest performance followed by DTM and Sentinel-2. However, explanatory capacities were slightly higher when multiple data sources were included in a regression (Table 4). The importance of Sentinel-1 and DTM in predicting CH 4 flux patterns was supported by predictor importance scores in which predictors calculated from these data were dominating among the top-ranked predictors in the regression including all remote sensing predictors (Fig. S1). For Sentinel-1, meanfiltered VH predictors from different incidence angles and window sizes were the most important ones. In the case of DTM, important predictors included those predicting wetness (TWI, DTW) as well as elevation and slope. For Sentinel-2, most important predictors included water vapor, blue and coastal aerosol bands and some spectral indices. In contrast, the role of other remote sensing-based predictors was small in the regression with all remote sensing predictors. In the correlation analysis, the presence of wetter fen plant community cluster, TWI, Sentinel-1 VH features, graminoid biomass and Sentinel-2 blue band had the highest absolute correlation with CH 4 (Table S12).
Although the proportion of area with net CH 4 emission was in minority within the landscape, the land area in total acted as a net source of CH 4 (Table 4). There was a considerable variation in the proportion of sources between different maps (12-50%), while the variation was in relative terms smaller in landscape-level mean flux estimations (253-502 lg CH 4 m -2 h -1 ).
When visually interpreting the spatial patterns of CH 4 , there were notable differences between maps produced with different predictor sets (Figs. 3, 4 and Supplementary material). Most of the regressionbased maps erroneously predicted CH 4 sources to mountain tundra in the western part of the study area. This pattern was the most evident with the Sentinel-1based map (Fig. 3d), but not visible in the DTM-based maps (Fig. 3f), in land cover classification-based maps (Figs. 1, 2 and 3g) nor in binary classification-based maps (Fig. 3h). In the DTM-based maps, CH 4 sources were predicted to forests on mineral soil nearby peatlands in the middle part of the study area. Moreover, many other regression-based maps predicted CH 4 sources to forests on mineral soil that were CH 4 sinks according to our fieldwork data (Fig. 2).
According to the standard error map (Fig. 3i), differences between maps were the largest in the peatland areas, mountain tundra areas and roadsides. When looking at pairwise map comparisons, it could be seen that differences between maps were the largest between maps with predictors calculated from single remote sensing data source and similarities increased when maps with predictors calculated from multiple data sources were compared (Fig. 4).
In comparison to the land cover classification-based map, most of the regression-based maps overestimated landscape-level fluxes, and all regression-based maps considerably overestimated the proportion of sources. In contrast, separate sink and source regression-based maps slightly underestimated landscape-level fluxes Minimum and maximum values are given in parentheses. In addition, predicted landscape-level mean CH 4 flux (landscape minimum and maximum in parentheses) and the proportion of area with net CH 4 emissions c Fig. 3 Maps of predicted CH 4 fluxes based on a both remote sensing and vegetation predictors, b vegetation predictors, c all remote sensing predictors, d Sentinel-1, e Sentinel-2, and f digital terrain model, g sinks and sources derived from land cover classification and regressions with all predictors, h sinks and sources derived from binary classification and regressions with remote sensing predictors, i standard error map over all CH 4 maps. Positive values indicate sources and negative values sinks and binary classifications proportion of sources (Table 4). Moreover, sink and source regression-based maps had the highest correlations with land coverbased map (Fig. 4).

Discussion and conclusions
According to our results, CH 4 fluxes can be predicted with a range of different remote sensing datasets, and remote sensing works better than plot-based vegetation measures in explaining CH 4 patterns ( Table 4). The latter is a somewhat surprising result as vegetation has been shown to be a major controlling factor of CH 4 fluxes (Turetsky et al. 2014;Abdalla et al. 2016). Furthermore, we assessed vegetation patterns in the plots where CH 4 fluxes were measured; thus, vegetation predictors were more directly linked to plot characteristics than remote sensing predictors. In upscaling, vegetation predictors and concurrent maps of vegetation have double uncertainty when compared to upscaling solely with remote sensing data as the vegetation maps need first to be produced with remote sensing data. Therefore, our results suggest that maps of vegetation may be an unnecessary step in an upscaling process. However, they are more useful when the amount of sampled vegetation plots is larger than the number of chamber plots, which typically is the case when portable CH 4 measurement devices are not used.
In regressions with all data, explanatory capacities were relatively high (Table 4) which is probably linked to high spatial heterogeneity in terms of topography and land cover within the study area and also due to considerable range in CH 4 fluxes. In previous studies, relatively high explanatory capacities have been obtained when modeling CH 4 flux patterns with remote sensing data in study areas containing both peatlands and forests on mineral soil (Dinsmore et al. 2017). This explanation is further supported by the observation that the explanatory capacities were lower in the separate sink and source regressions in which within-regression landscape heterogeneity and range of fluxes were smaller than in the regressions with all data. Explanatory capacities were higher in source than in sink regressions which is logical as the range of fluxes were higher in source regressions and as peatlands-where the CH 4 sources are mainly located-are spatially heterogenic in terms of vegetation, topography and CH 4 fluxes. High classification accuracies in vegetation mapping have been previously obtained within peatlands (Lehmann et al. 2016;Räsänen and Virtanen 2019). To test how well upscaling can be performed within a single peatland, more studies employing spatially dense and extensive flux measurements should be conducted. Also, datasets should have ultra-high spatial resolution, preferably \ 0.5 m grain size, due to fine-scale heterogeneities .
The landscape-level mean fluxes were relatively similar in the different maps (i.e. similar landscape composition), and the landscape was a net source of CH 4 according to all maps. However, the spatial pattern differed strikingly between the maps (i.e. the difference in configuration) ( Table 4, Figs. 3 and 4). The differences in spatial patterns were notably different in particular when looking at the threshold between sinks and sources. When compared to land cover and binary classification-based maps, regression-based maps overestimated the proportion of sources within the landscape and the mean CH 4 flux. On the one hand, the threshold between sinks and sources is probably most accurately predicted in the land cover or binary sink-source classification. This is because both land cover classification discriminating peatlands (mostly sources) from mineral soils (mostly sinks) and binary classification had high classification accuracies. On the other hand, it is more difficult to judge which one of the landscape-level mean flux estimates is the most accurate. Although land cover classification-based map models well the threshold between sinks and sources, regression-based maps can assess the variation of CH 4 fluxes within specific land cover types. Therefore, the separate sink and source regression-based maps could provide the most accurate estimates of spatial patterns and landscape-level mean fluxes although the regression performance was lower in the separate sink and source regressions than in the regressions with all data. As the classification accuracy in the binary classification was high, we recommend using it in regression-based CH 4 map production instead of a land cover map as it requires smaller and simpler training data. Nevertheless, probably the most feasible option is to produce and compare multiple different maps to evaluate the composition and configuration of different processes and biogeochemical fluxes within various landscapes.
The most important remote sensing predictors included Sentinel-1 and topography. This is in line with earlier studies in which the locations of peatlands have been mapped with SAR and topographic data with relatively high accuracy (Murphy et al. 2007;Widhalm et al. 2015;Karlson et al. 2019;Lidberg et al. 2020). The highest accuracies were obtained when remote sensing predictors were included from multiple sources, but contrary to our expectations, boosts in explanatory capacities with multi-source data were relatively small (Table 4). Even with a single data source, in particular Sentinel-1, high accuracy was obtained, but accuracy dropped considerably when Sentinel-1, Sentinel-2 and topography were omitted from the models. Nonetheless, in pairwise map comparisons, the spatial patterns were considerably different between maps with predictors from a single data source. Similarities increased when maps with more predictors were compared (Fig. 4), which suggests that the prediction of spatial patterns is more robust with multi-source remote sensing. Therefore, in line with previous literature (Bourgeau-Chavez et al. 2016;Karlson et al. 2019;Räsänen and Virtanen 2019), it is advisable to use multiple data sources to get both high explanatory capacities and realistic predictions of spatial patterns of CH 4 fluxes.
When looking at the importance of different topography predictors, the regression model with only three topographic predictors (TWI, elevation, DTW-NLS) had a higher explanatory capacity than a model with more topographic predictors emphasizing the importance of these three predictors in CH 4 upscaling (Table 4). Of these predictors, TWI and DTW should in principle model soil wetness (Murphy et al. 2007;Lidberg et al. 2020) which is a key control of CH 4 fluxes (Turetsky et al. 2014;Abdalla et al. 2016).
When looking at Sentinel-1 predictors, VH predictors were more important than VV predictors which is in line with previous research (Baghdadi et al. 2001;Karlson et al. 2019). Our results also showed that multiple different VH predictors (i.e. VH meanfiltered with different window sizes and acquired with different incidence angles) were among the most important predictors. This suggests that it is preferable to test, and, if necessary, include multiple diverse SAR datasets and preprocessing options in peatland mapping. Additionally, as SAR is sensitive to temporal moisture conditions (Bourgeau-Chavez et al. 2013), data taken on different time-points, even consecutive days, can give divergent results. Further boosts in regression performance could have been obtained if we had included SAR datasets from early and late growing season and also from winter (Widhalm et al. 2015;Karlson et al. 2019). This could have alleviated problems related to mountain tundra, which was erroneously modeled to be a significant emitter of CH 4 in the Sentinel-1-based map. Nevertheless, it has been found that SAR backscatter is more sensitive to vegetation structure than soil moisture in peatland environments (Millard and Richardson 2018). Moreover, in a study that incorporated multi-temporal SAR data in peatland detection, open peatlands were mixed up with mountain tundra, and the misclassifications were reduced with the help of topographic predictors (Karlson et al. 2019).
Interestingly, of the Sentinel-2 predictors, the most important ones were water vapor, blue and coastal aerosol bands instead of predictors, which are usually important in vegetation mapping, such as greenness indices (Pettorelli et al. 2005), or water indices sensitive to soil moisture (Arroyo-Mora et al. 2018). Our results suggest that in CH 4 flux studies, four spectral bands is not sufficient, and higher spectral resolution should be included. This finding is further backed by the low importance of PlanetScope data with four spectral bands in the regressions. When assessing CH 4 flux patterns within peatlands with finescale spatial heterogeneity, preferred datasets could be hyperspectral or multispectral sensors with at least ten spectral bands aboard unmanned aerial vehicles. However, the importance of different spectral features should be further researched; in particular because the most important spectral bands according to our study are primarily designated for atmospheric corrections (in particular water vapor and coastal aerosol bands) or are sensitive to atmospheric scattering even after the ESA Sen2Cor atmospheric correction (Drusch et al. 2012;Li et al. 2018).
There were some limitations in our study which should be addressed in future studies. Firstly, our measurement campaign covered only a snapshot situation during peak summer, and we did not quantify the seasonal or diurnal trends nor long-term balances of CH 4 . The broad spatial flux patterns in fluxes would probably be more or less similar also when evaluating a longer time span, i.e. mineral soils and peatlands would act as sinks and sources, respectively, but there would be changes in the magnitude of the fluxes within the catchment as well as within the vegetation and land cover types. This is backed by previous studies, which have reported that seasonal variation in CH 4 flux is largely controlled by variation in soil and air temperature and moisture (Chi et al. 2020), whereas diurnal variation is correlated in addition with changes solar radiation and latent heat flux (GaŽ ovic et al. 2010;Helbig et al. 2017;Long et al. 2010). These variables mainly affect the temporal trends in CH 4 fluxes, and spatial variation can mostly be explained by vegetation and land cover types (Forbrich et al. 2011;Davidson et al. 2016;Tuovinen et al. 2019;Chi et al. 2020); however, also soil temperatures may vary considerably between land cover types, in particular during summer (Helbig et al. 2017). Secondly, we omitted the fluxes originating from trees in our analysis. It has been shown that trees on mineral soil emit CH 4 ; for instance, emissions of an average Scots pine tree has been measured to be 0.8% of uptake by forest floor (Machacova et al. 2016). Therefore, future studies should concentrate on developing remote sensing-based spatiotemporal models for CH 4 and carbon balances that account for fluxes originating from open peatlands, forest floors and trees.
Recently developed portable carbon analyzers and modeling approaches presented in this article enable to measure and model landscape-level variation of CH 4 dynamics in such detail that spatially explicit maps can be produced. The maps produced with high resolution remote sensing constrain uncertainties related to CH 4 fluxes and their spatial patterns. The approach presented is applicable at various spatial extents, but field observation data should be collected so that it is representative for the whole region of interest, and at smaller extents, finer resolution remote sensing data acquired for example from unmanned aerial vehicles should be used.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.
Funding Open Access funding provided by University of Helsinki including Helsinki University Central Hospital.
Data availability The datasets generated during and/or analyzed during the current study are available from the corresponding author on reasonable request.