Upscaling Northern Peatland CO₂ Fluxes Using Satellite Remote Sensing Data

Abstract

Peatlands play an important role in the global carbon cycle as they contain a large soil carbon stock. However, current climate change could potentially shift peatlands from being carbon sinks to carbon sources. Remote sensing methods provide an opportunity to monitor carbon dioxide (CO₂) exchange in peatland ecosystems at large scales under these changing conditions. In this study, we developed empirical models of the CO₂ balance (net ecosystem exchange, NEE), gross primary production (GPP), and ecosystem respiration (ER) that could be used for upscaling CO₂ fluxes with remotely sensed data. Two to three years of eddy covariance (EC) data from five peatlands in Sweden and Finland were compared to modelled NEE, GPP and ER based on vegetation indices from 10 m resolution Sentinel-2 MSI and land surface temperature from 1 km resolution MODIS data. To ensure a precise match between the EC data and the Sentinel-2 observations, a footprint model was applied to derive footprint-weighted daily means of the vegetation indices. Average model parameters for all sites were acquired with a leave-one-out-cross-validation procedure. Both the GPP and the ER models gave high agreement with the EC-derived fluxes (R² = 0.70 and 0.56, NRMSE = 14% and 15%, respectively). The performance of the NEE model was weaker (average R² = 0.36 and NRMSE = 13%). Our findings demonstrate that using optical and thermal satellite sensor data is a feasible method for upscaling the GPP and ER of northern boreal peatlands, although further studies are needed to investigate the sources of the unexplained spatial and temporal variation of the CO₂ fluxes.

1. Introduction

Peatlands are an important ecosystem for the global carbon cycle, as they cover only 3% of the terrestrial area but contain around 40% of global soil organic carbon [1]. The waterlogged and anoxic conditions in peatlands inhibit the decomposition of plant material, creating a long-term net sink of atmospheric carbon dioxide (CO₂) and, consequently, an effective storage of organic carbon [2]. Peatlands assimilate CO₂ through photosynthesis (gross primary productivity, GPP) and release it through autotrophic respiration (from vegetation) and heterotrophic respiration (microbial decomposition). The sum of the autotrophic and heterotrophic respiration is the ecosystem respiration (ER). The difference between ER and GPP is the net ecosystem exchange (NEE), which determines whether an ecosystem is a carbon sink or source (excluding gains and losses of carbon through flowing water, biomass removal, or methane fluxes). Pristine peatlands are considered to be net carbon sinks [3], but changing environmental conditions due to global warming or land use change may increase CO₂ release and potentially shift peatlands from being carbon sinks to sources [4,5]. For that reason, it is important to estimate peatland CO₂ exchange across time and space.

A common method to measure land–atmosphere CO₂ exchange at the ecosystem level is the eddy covariance (EC) technique [6]. The EC method provides direct measurements of NEE and nighttime ER, which can be subsequently partitioned into GPP and daytime ER. As the flux partitioning relies on modelling procedures, EC-derived GPP and ER are not considered to be direct measurements but modelled estimates at the ecosystem level.

However, there is only a limited number of EC sites measuring peatland carbon exchange. Peatlands are widespread across the northern latitudes, and to enable the estimation of the CO₂ flux from all of these peatlands, it is important to develop remote sensing-based methods for upscaling their CO₂ fluxes. The estimation of peatland carbon fluxes with optical remote sensing presents a series of unique challenges due to several characteristic features of peatlands, such as the high proportion of non-vascular plants, high water content, and the effects of microtopographic variations [7]. As a result of their heterogeneity, carbon fluxes may vary widely between peatlands [8,9] and within peatlands [10,11].

Unlike CO₂ concentration, CO₂ fluxes are not directly detectable by satellites, and hence remote-sensing-derived estimates rely on the relationship between the fluxes and remotely sensed environmental variables. A widely used method to estimate GPP is the light use efficiency (LUE) model by Monteith [12], where GPP is the product of the total photosynthetically active radiation (PAR) incident on the vegetation, the fraction of photosynthetically active radiation absorbed by vegetation (fAPAR), and the conversion efficiency of the absorbed energy (ɛ). Several studies have shown that fAPAR is linearly or near-linearly correlated with spectral vegetation indices like the Normalized Difference Vegetation Index (NDVI) or the Enhanced Vegetation Index (EVI), and therefore a vegetation index often provides a proxy for fAPAR. Although NDVI is widely used in remote sensing-based ecosystem studies, EVI is less sensitive to perturbing factors like atmospheric effects (water vapor and aerosols) and soil reflectance, and more sensitive to variations in vegetation structure [13]. The two-band EVI (EVI2) has been developed to provide similar information as EVI, but it only uses the red and the near infrared (NIR) bands, whereas EVI also uses the blue band [14]. Despite the differences between the spectral indices, both NDVI and EVI have been successfully used to estimate GPP in peatland ecosystems [15,16].

To obtain the full CO₂ balance of an ecosystem, both GPP and ER need to be estimated, but as Lees et al. [7] highlight, remote sensing-derived peatland studies have mainly focused on estimating GPP with less consideration of ER. The key challenges are to estimate heterotrophic respiration with high accuracy using remote sensing data [17] and to explain the variability in ER across ecosystems [18]. Despite the challenges, ER has been successfully modelled using satellite land surface temperature (LST) [16,19], and some studies have also included vegetation productivity and soil moisture to improve the ER estimates [18,20,21,22].

The Moderate Resolution Imaging Spectroradiometer (MODIS) instrument mounted on NASA’s Terra and Aqua satellites is widely used to monitor vegetation dynamics, land cover changes, and LST at local and global scales, as it provides both optical and thermal data at a high temporal resolution (1–2 days). However, the heterogeneous character of peatland ecosystems presents a challenge for remote sensing using coarse-resolution (250 m–1 km) sensors like MODIS [23,24]. The MODIS GPP/NPP (MOD17) product at 1 km spatial resolution provides GPP estimates in different biomes [25]. However, the MOD17 product lacks a peatland category in its land cover classification and an estimate of surface moisture, which limits its applicability in peatlands.

The possibility to generate time-series data at high spatial (10–60 m) and temporal resolution (2–3 days at mid-latitudes) using the new Sentinel-2A and 2B satellites with the MultiSpectral Instrument (MSI) presents a new opportunity to map CO₂ fluxes in these landscapes. Nevertheless, even with the high frequency of data provided, the presence of clouds and seasonal snow cover cause data gaps of varying length. The TIMESAT v.4 software package provides a flexible yet robust data processing method based on box constrained double logistic functions combined with splines to gap-fill and smooth time series data [26]. The availability of high spatial resolution satellite data also raises the question of how to select a spatial sample of the data that is representative of the EC flux footprint area. Accounting for the flux footprint may be particularly important in peatlands due to their spatial heterogeneity, and can be achieved by using footprint models that can help select the satellite pixels that fall within the EC footprint area [27,28,29].

The potential of remote sensing-based models to predict peatland carbon fluxes has been demonstrated by Schubert et al. [16]. We will expand on this work by using newly available high spatial resolution (10–20 m) optical Sentinel-2 MSI data at five peatland sites across Sweden and Finland, state-of-the-art time-series processing, and footprint modeling to characterize the source area of the CO₂ fluxes. Although other carbon fluxes, such as methane (CH₄) emissions or dissolved carbon exported by water, also play an important role in peatlands, this study focuses on CO₂ fluxes, as these are typically the largest component of carbon exchange in peatlands at annual timescales [5,30].

The aim of this paper is to investigate the feasibility of developing a simple model based on remote sensing data to estimate the CO₂ balance (NEE) and its component fluxes (GPP and ER) across different peatlands types at northern latitudes, which could be used for upscaling CO₂ fluxes across northern boreal peatlands. As part of this aim, we investigate a) with what accuracy we can model the NEE flux and its components using a single parameter set for our five sites, compared to using individual parameter sets for each site, and b) with what accuracy directly modeled NEE differs from the sum of modelled GPP and ER.

2. Materials and Methods

2.1. Study Sites

The study includes five Nordic peatland sites located in Sweden and Finland: Abisko-Stordalen, Degerö, Mycklemossen, Siikaneva, and Lompolojänkkä (Figure 1,

Table 1. The study sites, their characteristics, and references for additional site descriptions. Temperature and annual precipitation values are 1981–2010 averages.

Site Name and Infrastructure	Location	Peatland Type	Vegetation Cover	Annual Precipitation and Air Temperature	Data Years	Reference
Abisko-Stordalen (SE-Sto) ICOS	68.356°N, 19.045°E	Sub-arctic ombrotrophic bog	Carex rostrata, Betula nana, Eriophorium angustifolium, Sphanum fuscum, Empetrum hermaphroditum	332 mm –0.1 °C	2017–2019	web- site 1
Lompolojänkkä (FI-Lom) ICOS	67.997°N, 24.209°E	Boreal medium rich fen	Carex rostrata, Menyanthes trifoliata, Betula nana, Salix lapponum, Sphagnum angustifolium, S. riparium, S. fallax	484 mm –1.4 °C	2017–2018	[31,32]
Degerö (SE-Deg) ICOS	64.182°N, 19.557°E	Boreal oligotrophic fen	Sphagnum balticum, S. Lindbergii, S. majus, Eriophorum vaginatum, Vaccinium oxycoccos L., Andromeda polifolia, Trichophorum caespitosum	613 mm 1.9 °C	2017–2019	[33]
Siikaneva (FI-Sii) ICOS	61.833°N, 24.193°E	Boreal oligotrophic fen	Carex chordorrhiza, C. Rostrata, Sphagnum papillosum, S. magellanicum, S. balticum, Salix phylicifolia, Betula nana	703 mm 3.5 °C	2017–2019	[5,34]
Mycklemossen (SE-Myc) SITES	58.365°N, 12.169°E	Hemi-boreal oligotrophic fen	Sphagnum rubellum L., Sphagnum fallax L., Sphagnum austinii L., Eriophorum vaginatum, Calluna vulgaris, Erica tetralix, Pinus sylvestris	803 mm 6.8 °C	2017–2018	website 2

¹https://www.icos-sweden.se/abisko-stordalen, accessed on 24 August 2020. ² https://www.fieldsites.se/en-GB/research-stations/skogaryd-32652394, accessed on 26 August 2020.

). Mycklemossen is part of the SITES (Swedish Infrastructure for Ecosystem Science) infrastructure, whereas the other sites belong to the ICOS (Integrated Carbon Observation System) infrastructure. The sites were selected to cover peatlands exposed to different climatic conditions (hemi-boreal, boreal, and subarctic), at northern latitudes from 58°N to 68°N. Four of the sites are fens, and one, Abisko-Stordalen, is classified as a bog, although it also contains areas of fen. The study period spanned between 2017 and 2019, though data availability varied among the sites. Further information on each site can be found in Table 1.

2.2. Eddy Covariance Flux Data

The eddy covariance (EC) measurement system at Lompolojänkkä included a USA-1 three-dimensional sonic anemometer (METEK GmbH, Elmshorn, Germany) and a closed-path LI-7000 CO₂/H₂O infrared gas analyzer (LI-COR, Inc., Lincoln, NE, USA). At Siikaneva and Mycklemossen a USA-1 anemometer was used with an LI-7200 (LI-COR Inc., USA) gas analyzer. At Degerö, a 1012R3 Solent (Gill Instruments, Lymington, UK) anemometer was used with an IRGA LI-6262 (LI-COR Inc., USA) gas analyzer and the Abisko-Stordalen site used a uSonic-3 Class A (METEK GmbH, Germany) anemometer with an LI-7200 gas analyzer.

From the high frequency (10 Hz) EC data, 30-min CO₂ fluxes were derived taking into account the corrections following [35]. At the ICOS Sweden sites, the open source EddyPro^® Eddy Covariance Processing Software, version 6.0 (LI-COR Biosciences, Lincoln, NE, USA), was used following the ICOS protocol for data processing [36,37]. The processing procedures at Lompolojänkkä and Siikaneva have been presented in Aurela et al. [31] and Korrensalo et al. [10], respectively. More details on the data processing at Mycklemossen are given in the Supplementary Materials. The processed EC fluxes were partitioned into GPP and ER and gap-filled using standard methods with an exponential temperature response for respiration [38] and a saturating light response for GPP [39].

Daily averages of the fluxes were calculated from the gap-filled 30-min CO₂ flux data acquired from the sites and the time series were smoothed using a spline function in TIMESAT (see Table S1 for the spline parameters). We applied the same spline parameters as used for the LST data (see Section 2.3) to reduce the high-frequency (daily) variations in the flux data that could not be captured by the remote-sensing data due to its coarser temporal resolution. The smoothing of the CO₂ flux data did not have a substantial effect on the seasonal shape of the fluxes or the annual NEE balance (see Section 3.3).

2.3. Remote Sensing Data

The models to estimate GPP, ER, and NEE were driven by vegetation indices based on Sentinel-2-retrieved data along with LST data from MODIS.

Sentinel-2 MSI images covering an area of 100 by 100 km at each site from 2017–2019 were downloaded from the European Space Agency (ESA) Copernicus Sentinels Scientific Data Hub. The Sen2Cor processor (version 2.8.0) was used to perform atmospheric corrections and obtain land surface reflectance and scene classification [40]. The vegetation index EVI2 [14], with a spatial resolution of 10 m, was calculated using red (RED) and near-infrared (NIR) reflectances:

$$\mathrm{EVI}2=2.5\frac{NIR-RED}{NIR+2.4RED+1} .$$

(1)

A vegetation index related to water content, the Normalized Difference Water Index (NDWI) [41], was calculated using 20 m resolution NIR and short-wave infrared (SWIR) bands:

$$\mathrm{NDWI}=\frac{NIR-SWIR}{NIR+SWIR} .$$

(2)

Both data sets were gap-filled and smoothed pixel-wise using a spline function in TIMESAT to produce daily values (Table S1).

In order to select the NDWI and EVI2 pixels that fell within the flux source area, the EC flux footprint was estimated at each site with the Flux Footprint Prediction (FFP) model [28]. The model was driven by 30-min micrometeorological data (collected as part of the EC measurements, see Section 2.2) to produce daily footprint climatologies. The footprint climatologies were used to calculate weighted daily means of the EVI2 and NDWI within 80% of the footprint area. Furthermore, the time-series of the mean daily EVI2 and NDWI were both scaled so that the index value was 0 when GPP was 0. This procedure enabled the combination and comparison of data from all the study sites, as the base values of the indices were usually not equal to zero and varied from site to site.

As NDWI expresses the effect of soil hydrology on the vegetation, NDWI was used to calculate a water scalar (W_s), following Xiao et al. [42]. The W_s represented the variations of moisture conditions interannually and between the sites. Other scalars that were tested are presented in Table S2. A simple version of the water scalar by Xiao et al. [43] was selected for the further analysis:

$${\mathrm{W}}_{\mathrm{s}}=1+{\mathrm{NDWI}}_{\max },$$

(3)

where NDWI_max is the 99th percentile of NDWI during the growing season. During the dormant season, W_s was set to 1. The growing season was defined for each site and year based on when the daily mean air temperature (T_air) exceeded +5 °C for seven consecutive days [9].

The LST data was acquired from the MODIS instruments mounted on the Terra and Aqua satellites. The MODIS LST products (MOD11A1, MYD11A1) include daily daytime and nighttime LST at 1 km spatial resolution, providing a maximum of four LST measurements per day. LST data were smoothed and gap-filled in the same way as the flux data and spectral indices using TIMESAT (Table S1). Due to the coarse resolution of the MODIS data, the footprint model was not used with these data, but we extracted data from one pixel where the EC tower and the peak of the flux footprint was located at each site. Daytime LST was used in the rest of the analysis due to the higher number of high quality images and a better fit with ground LST measurements that were available from some of the sites (R² = 0.85, RMSE = 5.3 °C for Degerö 2017–2019 and R² = 0.75, RMSE = 6.1 °C for Abisko-Stordalen 2017–2019).

2.4. Empirical Regression Models for GPP, ER, and NEE

The GPP model structure was based on previous work by Schubert et al. (2010), where EVI, LST, and photosynthetic photon flux density (PPFD) were used as regression model inputs. In this study, we investigated how satellite-derived EVI2, LST, and W_s along with site-measured environmental variables (T_air, PPFD, water table depth (WTD) and annual precipitation) were related with GPP time-series data. Based on the goodness-of-fit of the regressions, EVI2, W_s, and daytime LST were chosen as input variables into the GPP model. The derived linear regression model for estimating GPP at all the sites simultaneously can be described as:

$$\mathrm{GPP}=a\times \mathrm{EVI}2\times {\mathrm{W}}_{\mathrm{s}}\times \mathrm{LST}.$$

(4)

Because of the abovementioned scaling of EVI2 and NDWI, it was possible to force the regression line through the origin, which is a physically reasonable way to estimate GPP. Equation (4) can be regarded as an empirical variation of the LUE model, where the parameter a, W_s, and LST jointly function as a proxy for PAR and the light use efficiency coefficient. Since photosynthesis does not take place in cold temperatures, only data when LST was above 0 °C were used to determine the model parameter a.

Ecosystem respiration (ER) was also modelled using an empirical model. Three existing models were tested: Lloyd and Taylor [44] (Equation (5)), Heskel et al. [45] (Equation (6)), and Gao et al. [46] (Equation (7)), who added the dependence on GPP to Lloyd and Taylor’s model:

$$\mathrm{ER}=R_{ref}{e}^{E_0\left(\frac{1}{T_{ref}-T_0}-\frac{1}{T+273.15-T_0}\right)},$$

(5)

$$\mathrm{lnER}=a+bT+c{T}^2,$$

(6)

$$\mathrm{ER}=a\times \mathrm{GPP}+R_{ref}\times e^{E_0\left(\frac{1}{T_{ref}-T_0}-\frac{1}{T+273.15-T_0}\right)},$$

(7)

where T_ref = reference temperature (set to 10 °C); T = daytime LST; T₀ = –46.02 °C (value taken from Lloyd and Taylor [44]); GPP = modelled GPP; and R_ref, E₀, a, b, and c are parameters to be estimated. To model ER, we started by testing the three models (Equations (5)–(7)) at each of the sites separately, for all available years, in order to find the best model (see results in Table S3). We also tested multiplying each model by an EVI2 scalar (formulated in the same way as W_s) or the W_s scalar (Table S3). The Gao et al. [46] model (Equation (7)) provided the best fit. Including vegetation productivity in the ER model improved the fit, because it helped better capture the magnitude and timing of the peak growing season ER. We tested modelling ER separately during the growing and dormant seasons, but it yielded no or only minor improvements in the goodness-of-fit.

Since Equation (7) is sensitive to possible errors in the modelled GPP as an input, EVI2 was used as a proxy for modelled GPP to reduce the uncertainties in the modelled ER (Equation (8)). Equation (8) was the ER model used in the rest of the analysis:

$$\mathrm{ER}=a\times \mathrm{EVI}2+R_{ref}\times e^{E_0\left(\frac{1}{T_{ref}-T_0}-\frac{1}{T+273.15-T_0}\right)}.$$

(8)

We compared two methods for modelling NEE, where NEE = ER - GPP. We followed the sign convention that NEE is negative when the ecosystem is a carbon sink (GPP > ER) and NEE is positive when the ecosystem is a carbon source (ER > GPP). NEE was estimated 1) by subtracting the modelled GPP (Equation (4)) from the modelled ER (Equations (2) and (8)) by fitting a non-linear regression model that included the GPP and ER components (Equation (9)) to the EC NEE measurements:

$$\mathrm{NEE}=b_1\times \mathrm{EVI}2+R_{ref}\times e^{E_0\left(\frac{1}{T_{ref}-T_0}-\frac{1}{T+273.15-T_0}\right)}-b_2\times \mathrm{EVI}2\times {\mathrm{W}}_{\mathrm{s}}\times \mathrm{LST}.$$

(9)

In addition to these two NEE models, we also tested a third, purely empirical, NEE model based on multiple linear regression. EVI2, NDWI, LST, PPFD, WTD, and W_s were fitted as predictor variables against EC NEE measurements. However, this approach was not selected for further analysis, because PPFD and WTD were not available at all sites and years, and it did not provide a consistent improvement compared to the non-linear NEE model (Equation (9)). See Supplementary Materials for more details.

In the first stage of the analysis, the free parameters for the GPP, ER, and NEE models were estimated separately for each site using all the available years of flux data (Table S4), and in the second stage, the parameters were estimated for all the sites together (Tables S5–S7). To obtain one parameter set for all the sites, leave-one-out cross-validation (LOOCV) was used, i.e., the models were fitted using data from all the sites except for one site-year, which was used to validate the fitted models. After all the LOOCV runs, the mean model parameter values were applied to all the sites.

3. Results

3.1. Relationships between GPP, ER and Remote Sensing Variables

Daily mean EVI2 from Sentinel-2 showed a clear relationship with flux tower derived GPP, R² = 0.61 for all sites and years together (Figure 2a). Stronger relationships were found when studying each site and year separately, and differences in the relationship between GPP and EVI2 during the spring and autumn were also clearer at some sites (Figure S1a). These results indicate that it would be possible to use different regression models for the spring green-up and the autumn senescence, as the relationships seemed to be logarithmic and exponential, respectively, instead of linear. However, since it was not possible to define general functional relationships that could be used across all sites, we used linear relationships between EVI2 and GPP. In addition to EVI2, we used LST and W_s as proxies for the environmental variables controlling GPP, reducing the tendency to underestimate high GPP and increasing the strength of the relationship (R² = 0.72, Figure 2b).

There was also a strong exponential relationship between ER and LST, R² = 0.63 at all sites, and years together (Figure 3). Similarly to the relationship between EVI2 and GPP, some sites showed clear seasonal variation in the ER–LST relationship between the spring and autumn (Figure S1b).

3.2. GPP and ER Models

For GPP, the LOOCV parameterized linear regression model (Equation (4)), with the average coefficient a = 0.70 (Table S5), resulted in high agreement between modelled and EC-derived GPP at all the sites (R² = 0.54–0.78, RMSE = 0.42–0.98 µmol m⁻² s⁻¹,

Table 2. Goodness-of-fit statistics for the predicted GPP, ER, and NEE, using the LOOCV parameterized models, for all available site years. NRMSE is RMSE normalized using the range (maximum–minimum) of the EC-derived flux.

Site	Flux	R2	RMSE (µmol m−2 s−1)	NRMSE (%)
SE-Sto	GPP	0.76	0.42	10
	ER	0.23	0.41	19
	NEE (Equation (9))	0.59	0.26	10
	NEE (ER–GPP)	0.16	0.37	15
FI-Lom	GPP	0.78	0.98	12
	ER	0.68	0.62	14
	NEE (Equation (9))	0.57	0.79	12
	NEE (ER–GPP)	0.59	0.77	12
SE-Deg	GPP	0.68	0.48	13
	ER	0.56	0.41	16
	NEE (Equation (9))	0.34	0.31	11
	NEE (ER–GPP)	0	0.50	18
FI-Sii	GPP	0.73	0.59	15
	ER	0.85	0.31	10
	NEE (Equation (9))	0.33	0.39	15
	NEE (ER–GPP)	0	0.54	20
SE-Myc	GPP	0.54	0.93	20
	ER	0.51	0.98	18
	NEE (Equation (9))	0	0.41	15
	NEE (ER–GPP)	0	0.51	19
	GPP	0.70	0.68	14
Average	ER	0.56	0.54	15
	NEE (Equation (9))	0.34	0.43	13
	NEE (ER–GPP)	0	0.54	17

). Based on the normalized root mean square error (NRMSE), the best fit was found at Abisko-Stordalen, where the error was 10% of the range of the EC-derived GPP, whereas the worst fit was found at Mycklemossen (NRMSE = 20%).

For ER, the agreement between the LOOCV parameterized model fit (Equation (8), see Table S6 for parameters) and the EC-derived ER differed between the sites more than the GPP model fit: R² varied between 0.23–0.85, and RMSE was between 0.31–0.98 µmol m⁻² s⁻¹ (Table 2). The NRMSE of the ER model fit was similar to that of the GPP model, varying from 10% to 20%. The ER model performed best at Siikaneva, whereas the weakest relationships were found at Abisko-Stordalen and Mycklemossen. Using site-specific parameters (Table S4), the model NRMSE was reduced to between 7–12% of the flux range for both the GPP and ER estimates (

Table 3. Goodness-of-fit statistics for the site-specific GPP, ER and NEE models, for all available site years. NRMSE is RMSE normalized using the range (maximum–minimum) of the EC-derived flux.

Site	Flux	R2	RMSE (µmol m−2 s−1)	NRMSE (%)
SE-Sto	GPP	0.85	0.33	8
	ER	0.86	0.17	8
	NEE (Equation (9))	0.70	0.22	9
	NEE (ER–GPP)	0.55	0.27	11
FI-Lom	GPP	0.89	0.69	8
	ER	0.93	0.30	7
	NEE (Equation (9))	0.75	0.60	9
	NEE (ER–GPP)	0.64	0.73	11
SE-Deg	GPP	0.69	0.47	12
	ER	0.80	0.27	11
	NEE (Equation (9))	0.42	0.29	11
	NEE (ER–GPP)	0.06	0.38	14
FI-Sii	GPP	0.88	0.39	10
	ER	0.91	0.24	7
	NEE (Equation (9))	0.56	0.32	12
	NEE (ER–GPP)	0.24	0.42	16
SE-Myc	GPP	0.82	0.58	12
	ER	0.81	0.61	11
	NEE (Equation (9))	0.03	0.39	15
	NEE (ER–GPP)	0	0.62	23
	GPP	0.83	0.49	10
Average	ER	0.86	0.32	9
	NEE (Equation (9))	0.49	0.37	11
	NEE (ER–GPP)	0.01	0.48	15

).

The time-series of the modelled and EC-derived fluxes (Figure 4 for GPP, Figure 5 for ER) show that the LOOCV parameterized models are able to capture the seasonal variation of the fluxes. However, some noticeable under- and overestimations of GPP and ER occur during the peak of the growing season. The LOOCV parameterized ER model also tended to underestimate the wintertime ER fluxes, especially at Lompolojänkkä and Mycklemossen (Figure 5c,i). Using the models with site-specific parameters helped better capture the peak of the fluxes, with less improvement in the general seasonal shape. The time-series of the LOOCV parameterized modelled GPP for Mycklemossen and Lompolojänkkä show a delayed start of the growing season compared to the EC-derived GPP (Figure 4c,i). Similar effects can be observed at Abisko-Stordalen in 2019 and Degerö in 2017 (Figure 4a,e). As noted above, the LOOCV parameterized GPP and ER models performed poorly at Mycklemossen, and Figure 4i and Figure 5i show that both fluxes were significantly underestimated at this site. Several environmental variables related to GPP and ER were tested in order to find a scaling factor that could help to improve the ability of the models to capture the growing season peak of the fluxes (Table S2), but no strong correlations were found. The NDWI-based water scalar (W_s) was selected for GPP modelling, as it improved the agreement between modelled and EC-derived GPP at Lompolojänkkä and Mycklemossen. The influence of this variable was small at the other sites, and it did not improve the accuracy of modelled ER (Table S3).

3.3. NEE Models

Out of the two NEE models, the non-linear regression model (Equation (9)) generally performed better than calculating the difference of modelled GPP (Equation (4)) and modelled ER (Equation (8)) when compared against EC-measured NEE (see Table 2 for LOOCV parameterized model prediction estimates and Table 3 for site-specific estimates). The site-specific model parameters (Table S4) gave higher agreement than the average parameters from the LOOCV (Table S7). Overall, the performance of the NEE models was weaker than the GPP and ER models, although at Abisko-Stordalen and Lompolojänkkä the non-linear regression NEE model with site-specific parameters gave similar agreement to that of GPP and ER, R² > 0.70, and NRMSE around 9% of the range of measured NEE (Table 3).

The time-series of modelled and measured NEE show that modelled NEE corresponds fairly well to measured NEE, and the seasonal shape of NEE was captured at most of the sites and years (Figure 6). However, the models have difficulties capturing the growing season peak, as seen for instance at Abisko-Stordalen in 2017 (Figure 6a), and the periods of positive NEE (when ER is higher than GPP) at the start and end of the growing season. The latter can be seen clearly at Lompolojänkkä in both years (Figure 6c). The poor performance of the GPP and ER models at Mycklemossen are also shown in the NEE time-series, where the modelled NEE corresponds weakly to the measured NEE (Figure 6i).

The comparison of the annual cumulative modelled NEE and measured NEE (Figure 7) emphasizes the importance of accurate flux modelling, as small over- or underestimations, not only during the peak growing season, but during the whole year, can result in a significant difference in the annual CO₂ budget. For example, the NEE models at Abisko-Stordalen in 2017 underestimated the NEE growing season peak equally (Figure 6a), but result in very different annual accumulation of CO₂. The model with the site-specific parameters suggests that the site is a CO₂ sink, whilst the model with the average parameters suggests that the CO₂ budget is close to zero (Figure 7). Table S8 shows that the modelled annual NEE differed largely from the measured annual NEE at most of the sites and years. The model was able to estimate the annual NEE accurately (maximum 3% difference to the measured NEE) at Abisko-Stordalen 2017 (site-specific estimate), Abisko-Stordalen 2018 (predicted with joint parameters), and Siikaneva 2017 (site-specific estimate). The linear regression model with site-specific parameters performed slightly better than the model with jointly estimated parameters at eight out of 13 site-years.

3.4. Upscaling GPP to the Peatland Scale

In order to demonstrate the applicability of the simple regression model with high resolution satellite data for upscaling, we applied the GPP model (Section 3.2, Equation (4)) with the average parameters from the LOOCV (Table S5) to Sentinel-2 data for the Lompolojänkkä site (Figure 8). We selected Sentinel-2 pixels for upscaling based on the wetland land cover classification in Finland (SYKE, 2018), so that only pixels in the class “open peatbog” were used. Figure 8 shows the estimated average GPP during the growing season in 2017 and 2018. GPP varies from 1 up to 5 µmol m^-2 s^-1 within the flux footprint area, and similar spatial variation can be seen also beyond the flux footprints in both years.

4. Discussion

The results of this study showed that purely remote sensing (RS) derived data can be used in regression models to estimate seasonal CO₂ uptake and release from northern peatlands. The product of EVI2, LST, and an NDWI-based water scalar (W_s) explained a large proportion of the variability in GPP at all the sites. In addition, ER was strongly related to daytime LST, and based on these findings, NEE could be estimated from LST in combination with EVI2 and W_s. However, the accuracy of the estimated fluxes was inconsistent across the sites and years. It should be noted that the EC-derived fluxes vary at a high frequency that satellite remote sensing data cannot reproduce. Smoothing the time-series of EC-derived NEE, ER and GPP made them more comparable to the remote sensing data. Nevertheless, the RS models rely on LST to capture the daily variations in the fluxes, because LST is the only remote sensing dataset that varies significantly at a daily scale, whereas EVI2 and NDWI vary more slowly over the growing season. This means that short-term variability can only be captured at a coarse spatial scale.

Some sites showed clear variations in the ER-LST and the GPP-EVI2 relationships between the spring green-up phase and the senescence phase, suggesting that each phase should be modelled separately. However, during initial testing, we found that separating the CO₂ fluxes between different seasonal stages was not straightforward. A simple split of a flux time series between the green up phase and the senescence phase can cause problems with modelling the peak season accurately, and therefore we decided to model the whole year at once. Further analyses may reveal solutions to this challenge.

The linear regression model for GPP that we developed is an empirical variation of the formal LUE model by Monteith [12]. LUE models usually include PAR. In this study, PAR was not included in the GPP model, as it showed very weak relationships with daily GPP at all the sites. Instead, the regression parameter, the water scalar W_s, and LST jointly provided a proxy for PAR and the LUE term, but these factors might not be sufficient to explain all spatial and temporal variations of the photosynthetic efficiency within a peatland ecosystem. Alternative ways to estimate the LUE coefficient using remote sensing methods, such as the photochemical reflectance index (PRI) [47], may solve this issue [48]. Unfortunately, the wavelength bands to calculate PRI are not available from Sentinel-2. Another remote sensing-derived variable that can be used to estimate GPP is solar induced fluorescence (SIF) [49], which will be operational at the ecosystem scale in future from the forthcoming Fluorescence Explorer (FLEX) satellite mission by ESA.

LST explained ER well at all sites. Including productivity in the ER model (using EVI2 or modelled GPP) improved the model fit compared to the original Lloyd and Taylor [44] model. Gao et al. [46] argue that the Lloyd and Taylor equation represents heterotrophic respiration, and including GPP is a way of including information related to autotrophic respiration in the estimate. It has been shown that ER is tightly linked to vegetation productivity [50,51,52], and a productivity or greenness index has also been used in remote sensing models of ER before [18]. Our results demonstrate that using EVI2 instead of modelled GPP (as proposed by Gao et al. [46]) gives similar model accuracy whilst also making the model simpler and more robust. However, the LST and EVI2 based model might be less suitable to estimate ER in winter, since low soil surface temperatures and vegetation productivity predict low fluxes, while higher temperatures in deeper layers of peat enable decomposition to continue [53]. Wintertime ER can be an important part of the annual CO₂ balance in northern peatlands that is possibly underestimated by surface temperature driven models [54]. Separate ER model fits for the whole dormant season did not improve ER estimates, but future work could investigate this using separate model fits specifically during periods with snow cover (similar to Jägermeyr et al. [18]) or using LST as input to model peat temperature during the winter.

Of the two types of NEE models tested, the non-linear regression model performed generally better than the model that was calculated by subtracting GPP from ER. This might be due to error propagation [16]. It is challenging to model NEE directly, as it is the difference between two large fluxes driven by separate processes with different seasonal dynamics. Although the largest under- or overestimations were found during the peak growing season, our results indicate that small differences during the whole year greatly influence the annual NEE as well.

The modelling results showed that site-specific parameterization usually gave higher agreement with the EC fluxes than the average parameters from the LOOCV runs. However, for upscaling purposes, an average parameter set is more applicable than parameters that are strictly associated with one site. To reduce the importance of the regression parameters, other variables, that are able to explain the differences between the years within a site and also between the sites, could be included in the models. We tested several variables to scale GPP and ER in order to capture differences in the growing season fluxes between the sites, but none of them were strongly related to peak GPP or ER.

The NDWI-based water scalar was selected for GPP modelling, as in general, WTD is one of the most important variables regulating the function of peatland ecosystems [55]. Measures of soil moisture availability or precipitation have also been used previously to model ER [21], but our initial tests found that including the NDWI or W_s in the ER model did not improve the fit. NDWI is a spectral index expressing the influence of soil moisture on the vegetation rather than an actual soil moisture index, and it does not fully capture the seasonal variations in water table depth within a site (see Figure S2) or differences between sites (Table S2). Further work is needed to find remote sensing variables better able to capture variations in water table depth or soil moisture at peatland sites; for this purpose, hyperspectral and synthetic aperture radar data have shown promising results [56,57]. Alternatively, the CO₂ flux modelling could be combined with hydrological modelling to further investigate how, for example, water table depth changes or snow melt timing affect the fluxes. For instance, Huang et al. [58] assimilated regional soil moisture network data, remote sensing data, and high-resolution land surface parameters to develop an empirical soil moisture model.

Our modelling approach investigated to what extent a single set of parameters is sufficient to upscale the CO₂ fluxes from a variety of northern peatlands. Future work could examine which improvements may be gained by modelling the fluxes of different peatland types (i.e., bogs and fens) separately. For example, our ER model with average parameters performed worst at Abisko-Stordalen, the only ombrotrophic bog in our study. Similarly, Kross et al. [15] found differences in the relationships between remote sensing vegetation indices and peatland CO₂ fluxes depending on the presence of trees and peatland type. Bogs and fens have also been shown to respond differently to changes in water table and air temperature [59,60] due to their contrasting vegetation assemblages.

Improving the GPP and ER estimations separately would also improve the NEE model estimation. An alternative approach to estimate NEE could be a multivariate linear regression model including all remote sensing and environmental variables that are available at a site. This method lacks the more mechanistic scheme that our NEE model (Equation (9)) has, but the advantage is the high number of free parameters that makes the model adjustable. This method could be used to estimate NEE at an individual site that differs from other sites in its characteristics, e.g., at Mycklemossen, which is more light-dominated than the other sites in our study (Figure S3).

During the summer of 2018, northwestern Europe experienced a severe drought, which increased air temperatures, lowered WTD, and thus reduced the CO₂ uptake at several northern peatlands [61]. Rinne et al. [61] found that the 2018 drought did not affect all peatlands equally, as local hydrological features can make a peatland less sensitive to climatic variations, as was the case at Lompolojänkkä. Similarly, compared to 2017 and/or 2019, the EC-derived NEE we present was lower in 2018 at Abisko-Stordalen, Degerö and Siikaneva, but not at Lompolojänkkä. The decrease in NEE at these sites mainly originated from reduced GPP, although increased ER in 2018 was observed at Siikaneva. Lower GPP and ER were found at Degerö in 2019, which was also a dry year at that site. However, the dry conditions in 2019 did not influence the spectral properties of the vegetation, as the EVI2 in 2019 is at the same level as in 2017, which causes the mismatch between modelled and measured GPP. It is also possible the 2018 drought has a legacy effect in 2019, i.e., the vegetation did not fully recover from the drought in the previous year, and therefore both GPP and ER were lowin 2019 at Degerö. Mycklemossen experienced droughts both in 2017 and 2018 [61], and was also the southern-most (i.e., warmest and most light-dominated) site in the study, which may partly explain why our models with the average parameters performed poorest here. Overall, with only two or three years of data, it is difficult to draw strong conclusions about the temporal variations of the fluxes, especially if a severe drought has been experienced at a site.

Despite the limitations discussed above, our findings suggest that using satellite data in regression models is a simple yet powerful tool for modelling and upscaling GPP and ER. Due to its high spatial resolution, Sentinel-2 data has great potential to upscale GPP from the flux footprint to larger areas. The main input of the ER model, MODIS LST, has a coarse spatial resolution (1 km), which makes it less feasible to study heterogeneous ecosystems like peatlands. However, in this study, the MODIS LST was chosen over the 30 m resolution thermal Landsat data because of the higher temporal resolution and because it is a readily available LST product. To improve the spatial resolution of MODIS data, an appropriate downscaling procedure from MODIS reflectance to the Landsat 30 m resolution could be used [62]. Preserving the high temporal resolution while increasing the spatial resolution would enable the mapping of the spatial variation of LST within a peatland and within the flux tower footprint and may lead to improved ER and GPP models.

5. Conclusions

In this study, we demonstrated that regression models driven by remote sensing-derived data can be used to estimate CO₂ fluxes in northern peatlands (Abisko-Stordalen, 68°N; Lompolojänkkä, 68°N; Degerö, 64°N; Siikaneva, 62°N; Mycklemossen, 58°N) at relatively high accuracy. A strong relationship was found between GPP and the product of Sentinel-2 EVI2, Sentinel-2-derived water scalar (W_s), and daytime LST from MODIS. ER was strongly related to MODIS LST, and including EVI2 as a proxy for GPP improved the agreement with eddy covariance-derived ER. NEE was successfully explained by MODIS LST in combination with Sentinel-2 EVI2 and W_s at Abisko-Stordalen and Lompolojänkkä. At other sites, however, the relationship between modelled NEE and EC-measured NEE was weak. The regression models performed better when the model parameters were estimated separately for each site in comparison to the average parameters acquired with a Leave-One-Out-Cross-Validation (LOOCV) procedure, but the LOOCV parameterization is more appropriate for upscaling the CO₂ fluxes to a regional scale.

We conclude that there is great potential to upscale GPP and ER from a site to the regional level using solely remote sensing-derived data. With high-resolution Sentinel-2 data, we were able to take advantage of a flux footprint model that maps the flux source area and thus minimize the uncertainty from mismatches in the scale of the EC and remote sensing data. It also enabled us to map the spatial variability of GPP within a peatland. In addition, there is an opportunity to develop high spatial and temporal resolution models of ER by downscaling coarse MODIS data to finer resolution or using the forthcoming Landsat Collection 2 LST product. Overall, modelling peatland NEE at regional scale still poses a challenge, and in order to fully utilize the potential of remote sensing data to upscale peatland CO₂ fluxes, further studies would be needed to examine the yet unexplained spatial and temporal variation of these fluxes.