Modelling Daily Gross Primary Productivity with Sentinel-2 Data in the Nordic Region–Comparison with Data from MODIS

: The high-resolution Sentinel-2 data potentially enable the estimation of gross primary productivity (GPP) at ﬁner spatial resolution by better capturing the spatial variation in a heterogeneous landscapes. This study investigates the potential of 10 m resolution reﬂectance from the Sentinel-2 Multispectral Instrument to improve the accuracy of GPP estimation across Nordic vegetation types, compared with the 250 m and 500 m resolution reﬂectance from the Moderate Resolution Imaging Spectroradiometer (MODIS). We applied linear regression models with inputs of two-band enhanced vegetation index (EVI2) derived from Sentinel-2 and MODIS reﬂectance, respectively, together with various environmental drivers to estimate daily GPP at eight Nordic eddy covariance (EC) ﬂux tower sites. Compared with the GPP from EC measurements, the accuracies of modelled GPP were generally high (R 2 = 0.84 for Sentinel-2; R 2 = 0.83 for MODIS), and the differences between Sentinel-2 and MODIS were minimal. This demonstrates the general consistency in GPP estimates based on the two satellite sensor systems at the Nordic regional scale. On the other hand, the model accuracy did not improve by using the higher spatial-resolution Sentinel-2 data. More analyses of different model formulations, more tests of remotely sensed indices and biophysical parameters, and analyses across a wider range of geographical locations and times will be required to achieve improved GPP estimations from Sentinel-2 satellite data.


Introduction
Atmospheric carbon dioxide (CO 2 ) from anthropogenic emissions is one of the key drivers of climate change [1].The techniques used to quantify land-atmosphere carbon exchanges are essential for the understanding of the global carbon cycle [2].Especially Northern land ecosystems play an important role in the global carbon cycle due to the high amount of carbon stored in boreal ecosystems [3] and the increasing seasonal CO 2 exchange [4].Micrometeorological eddy covariance (EC) measurements, conducted at fixed towers, can provide in situ estimates of CO 2 exchange at the ecosystem level, from which Gross Primary Productivity (GPP) can be estimated [5].However, these measurements lack the spatial coverage needed for the characterization of CO 2 fluxes at larger regional extents (10 3 -10 6 km 2 ), and to represent different age classes and species diversity [6].
An efficient way to estimate GPP across larger areas is to use observations from Earth observation satellites.Satellite-derived approaches for carbon modelling are usually based Remote Sens. 2021, 13, 469 2 of 18 on the light use efficiency (LUE) model [7], empirical models with site-or ecosystem-specific environmental variables [8], or mechanistic models with remotely sensed inputs [9][10][11].LUE models express GPP as a product of the photosynthetic light use efficiency (ε), and photosynthetically active radiation absorbed by the vegetation (APAR).Although APAR is nonlinearly correlated to GPP at the daily time step [12,13], the relationship can be considered linear over monthly or annual periods [14][15][16].APAR can be derived by multiplying PAR (photosynthetically active radiation) by fAPAR (fractional PAR absorbed by the vegetation).Several studies have shown that satellite-derived spectral vegetation indices have linear or near-linear relationships with fAPAR within different vegetation types and climatic conditions and are therefore used as proxies of fAPAR [17][18][19][20].The LUE factor, on the other hand, is more difficult to estimate.LUE models usually rely on biome-specific maximum LUE factors that are reduced due to environmental constraints such as air temperature and water content [21][22][23].However, meteorological data and a biome-specific maximum LUE factor might be too coarse to accurately describe spatial and temporal variations of the LUE factor [11,22,24]. Due to this issue, empirical models based on remotely sensed input have been developed for estimating terrestrial GPP [25][26][27].In these models the regression parameter represents the LUE factor since it scales fAPAR to the level of carbon assimilation.The simplest empirical models use only a vegetation greenness index such as the Normalized Difference Vegetation Index (NDVI) proposed by Tucker [28] or the Enhanced Vegetation Index (EVI) proposed by Huete et al. [20] to estimate GPP.In comparison to the widely used NDVI, EVI is more adjustable for reducing soil and atmospheric effects and more sensitive to structural variations in vegetation [20].The further developed two-band Enhanced Vegetation Index (EVI2) provides similar information on vegetation properties to EVI but only uses the red and near-infrared bands, whereas EVI also requires the blue band [29].However, Sims et al. [24] showed that a model based entirely on EVI only poorly explains the seasonal variation of GPP at sites with evergreen vegetation that experience summer drought.Especially in boreal forests the air temperature (T a ) and the availability of light (PAR) regulate the seasonal photosynthetic activity [30].In addition, summer drought is typically the result of high T a and increased vapor pressure deficit (VPD) [31].Consequently, including PAR, T a and VPD variables in the model along with a remotely sensed vegetation index is assumed to improve the GPP estimation.
When calibrating remote sensing models against EC CO 2 data, it should be noted that the EC data represent large exchange areas, roughly covering a distance of 10 to 100 times the height of the measurements over displacement height [32,33].The area influencing the measured flux is called the source area [34], the contribution of each part to the measured flux (source weight) is called the flux footprint, and it depends on wind direction, wind speed and turbulence, leading to the variation of footprint from one half-anhour flux averaging period to another [35].The rough approximation above disregards this variation in flux footprint, thus introducing scale and weather-dependent approximation errors [36,37].For accurate spatial representation of flux variations, a footprint model should be used [38,39].This type of model can be used to assess which pixels belong to the footprint and hence should be selected for the calibration.
To model GPP on local and global scales, data from the Moderate Resolution Imaging Spectroradiometer (MODIS) sensor on board NASA's Terra and Aqua satellites have been widely used.Due to its high temporal resolution, MODIS provides data for monitoring vegetation dynamics and changes in land cover over time.However, MODIS acquires data at three relatively low spatial resolutions: 250 m, 500 m, and 1 km.The widely used MODIS GPP/NPP (MOD17, collection 6) product at 500 m spatial resolution has been questioned by many, mainly due to a lack of sufficient flux measurements for validation and calibration of the algorithm [40][41][42].
The spatial resolution of satellite products has been identified as a critical factor that limits spatially detailed carbon flux estimates, especially in heterogeneous landscapes [43].The new Sentinel-2A and 2B satellites with the MultiSpectral Instrument (MSI) have great potential for studying CO 2 fluxes in terrestrial ecosystems due to improved spatial resolution (10 m, 20 m, and 60 m) and high temporal resolution (2-5 days' revisit frequency depending on latitude).A few studies have explored the performance of the Sentinel-2 products for CO 2 flux modelling, e.g., Lin et al. [44] applied 30 m Sentinel-2 data from the Harmonized Landsat and Sentinel-2 product to estimate the GPP (R 2 = 0.77 and RMSE = 0.81 g C m −2 day −1 comparing to EC GPP) of three evergreen broadleaf forest sites and two grassland sites in southeast Australia.In the study by Wolanin et al. [45], Sentinel-2 and Landsat 8 data were used to estimate GPP (R 2 = 0.92 and RMSE = 1.38 g C m −2 day −1 ) on five crop fields (four in the USA and one in Germany).Given that Sentinel-2 has been shown to be capable of accurate representation of phenology [46,47] and has the capability to provide detailed estimates of biophysical and biochemical vegetation parameters [48,49] it can be hypothesized that Sentinel-2 MSI data has the potential to enhance the performance of GPP models also in northern landscapes.However, the impact of the lower temporal frequency of Sentinel-2 compared to MODIS in these high-cloudiness areas is a factor that could have negative impact.
In this study, we analyze eight study sites within the Nordic region (Denmark, Finland and Sweden, Figure 1) to investigate the accuracy of GPP estimates from Sentinel-2 data in comparison with MODIS data.We base the comparison on empirical regression models as they are commonly used for estimating GPP, and are functionally similar to the LUE type models.Following this, we examine the relationships between EC-derived GPP and EVI2-derived GPP using EVI2 and meteorological variables.In investigating the potential of Sentinel-2 vegetation index data for modeling carbon uptake, the study aims to (1) examine the effects of using Sentinel-2 MSI 10-m data compared to MODIS 250 m and 500 m data on modeling GPP from EC measurements, and (2) determine the importance of the input data for modeling accuracy in comparison to the choice of environmental driving variables.
The spatial resolution of satellite products has been identified as a critical factor that limits spatially detailed carbon flux estimates, especially in heterogeneous landscapes [43].The new Sentinel-2A and 2B satellites with the MultiSpectral Instrument (MSI) have great potential for studying CO2 fluxes in terrestrial ecosystems due to improved spatial resolution (10 m, 20 m, and 60 m) and high temporal resolution (2-5 days' revisit frequency depending on latitude).A few studies have explored the performance of the Sentinel-2 products for CO2 flux modelling, e.g., Lin et al. [44] applied 30 m Sentinel-2 data from the Harmonized Landsat and Sentinel-2 product to estimate the GPP (R 2 = 0.77 and RMSE = 0.81 g C m −2 day −1 comparing to EC GPP) of three evergreen broadleaf forest sites and two grassland sites in southeast Australia.In the study by Wolanin et al. [45], Sentinel-2 and Landsat 8 data were used to estimate GPP (R 2 = 0.92 and RMSE = 1.38 g C m −2 day −1 ) on five crop fields (four in the USA and one in Germany).Given that Sentinel-2 has been shown to be capable of accurate representation of phenology [46,47] and has the capability to provide detailed estimates of biophysical and biochemical vegetation parameters [48,49] it can be hypothesized that Sentinel-2 MSI data has the potential to enhance the performance of GPP models also in northern landscapes.However, the impact of the lower temporal frequency of Sentinel-2 compared to MODIS in these high-cloudiness areas is a factor that could have negative impact.
In this study, we analyze eight study sites within the Nordic region (Denmark, Finland and Sweden, Figure 1) to investigate the accuracy of GPP estimates from Sentinel-2 data in comparison with MODIS data.We base the comparison on empirical regression models as they are commonly used for estimating GPP, and are functionally similar to the LUE type models.Following this, we examine the relationships between EC-derived GPP and EVI2-derived GPP using EVI2 and meteorological variables.In investigating the potential of Sentinel-2 vegetation index data for modeling carbon uptake, the study aims to (1) examine the effects of using Sentinel-2 MSI 10-m data compared to MODIS 250 m and 500 m data on modeling GPP from EC measurements, and (2) determine the importance of the input data for modeling accuracy in comparison to the choice of environmental driving variables.

Study Sites and Environmental Variables
The study area contains eight EC flux sites located in Denmark, Finland and Sweden, at latitudes from 55 • N to 68 • N (Figure 1).The flux sites belong to the ICOS infrastructure (https://www.icos-ri.eu)and follow the ICOS measurement protocols [50].We analyzed data from 2016 to 2017, except at the SE-Svb site where we used only the data from 2016.There is an air-temperature gradient from south to north across the eight sites, with climate zones ranging from temperate in the south, through hemiboreal and boreal, and up to the northernmost subarctic landscapes.The sites cover diverse and major Nordic land cover types: evergreen forest, deciduous forest, wetlands, and agriculture (Table 1).At the ICOS sites, auxiliary meteorological parameters, such as photosynthetically active radiation (PAR), air temperature (T a ) and water vapor pressure deficit (VPD) are measured using standardized methodology [51,52].PAR is measured as instantaneous photosynthetic photon flux density (PPFD, µmol m −2 s −1 ) at the sites and was summed to daily total PPFD (µmol m −2 d −1 ).The measurements of T a ( • C), relative humidity, wind components and gas concentrations were located above the canopies at the forest sites and at around 2 m above the ground at the agriculture and wetland sites.VPD (hPa) is estimated from T a and relative humidity.The VPD and T a data used in the study were daily averages.
The net ecosystem exchange of CO 2 (NEE) at ICOS sites is measured by the EC method, with standardized methods and data processing [53,54].GPP (g C m −2 day −1 ) was estimated by partitioning the NEE, with the help of auxiliary in situ environmental variables, such as global radiation, air or soil temperature, and VPD, by using the REd-dyProc tool [55,56], implemented as an online tool (https://www.bgc-jena.mpg.de/bgi/index.php/Services/REddyProcWeb).

Eddy Covariance Influence Area
To estimate the approximation of the climatological footprint of the EC measurements, we used the Matlab implementation of the Flux Footprint Prediction (FFP) model of Kljun et al. [39], which provides position, size, and contribution details of surface source areas.Most of the input variables were collected as a part of the EC systems at the sites.The displacement height was defined as 0.67 times the canopy height.The boundary layer height was not measured directly at the sites but estimated using the friction velocity, Obukhov length and latitude of the site for near-neutral and stable conditions following Kljun et al. [39].For convective conditions the boundary layer was set to be 1500 m.
Through input half-hourly measured data (meeting the condition that the friction velocity was higher than 0.1 m s −1 ) to the FFP model, we obtained annual flux footprint climatology contours (Figure 2) and contribution weights in 10 m resolution pixel (matched to Sentinel-2 pixel).Following the suggestion in Kljun et al. [39], we used the 80% flux footprint contribution contour line as the maximum footprint area since it is sufficient for estimating the main source area in most cases.
Remote Sens. 2020, 12, x FOR PEER REVIEW 5 of 18 Obukhov length and latitude of the site for near-neutral and stable conditions following Kljun et al. [39].For convective conditions the boundary layer was set to be 1500 m.Through input half-hourly measured data (meeting the condition that the friction velocity was higher than 0.1 m s −1 ) to the FFP model, we obtained annual flux footprint climatology contours (Figure 2) and contribution weights in 10 m resolution pixel (matched to Sentinel-2 pixel).Following the suggestion in Kljun et al. [39], we used the 80% flux footprint contribution contour line as the maximum footprint area since it is sufficient for estimating the main source area in most cases.

Daily EVI2 Trajectories of Satellite Pixels
The eight ICOS sites are covered by seven Sentinel-2 MSI tiles.A total of 1312 available Sentinel-2 images acquired in 2016 and 2017 for the seven tiles were downloaded from the ESA Copernicus Sentinels Scientific Data Hub [57].The Sen2Cor processor (version 2.5.5) was applied for atmospheric correction to obtain Level-2A surface reflectance and scene classification layer (SCL) [58].EVI2 of each Sentinel-2 pixel (10m resolution) was calculated by red reflectance ( ) and near-infrared reflectance ( ) using the following equation: Through the SCL, pixels recorded as vegetation and soil were marked as high quality, whereas those recorded as cloud, snow, and water, were marked as low quality.Two different spatial resolution MODIS EVI2 (250 m and 500 m) datasets for 2016 and 2017 were generated from MODIS/Terra MOD09Q1 (250 m) and MOD09A1 (500 m) [59] to be comparable to Sentinel-2 EVI2.Each MODIS pixel's quality flag was created from the binary MODIS quality (QA) data.

Daily EVI2 Trajectories of Satellite Pixels
The eight ICOS sites are covered by seven Sentinel-2 MSI tiles.A total of 1312 available Sentinel-2 images acquired in 2016 and 2017 for the seven tiles were downloaded from the ESA Copernicus Sentinels Scientific Data Hub [57].The Sen2Cor processor (version 2.5.5) was applied for atmospheric correction to obtain Level-2A surface reflectance and scene classification layer (SCL) [58].EVI2 of each Sentinel-2 pixel (10m resolution) was calculated by red reflectance (ρRED) and near-infrared reflectance (ρN IR) using the following equation: Through the SCL, pixels recorded as vegetation and soil were marked as high quality, whereas those recorded as cloud, snow, and water, were marked as low quality.Two different spatial resolution MODIS EVI2 (250 m and 500 m) datasets for 2016 and 2017 were generated from MODIS/Terra MOD09Q1 (250 m) and MOD09A1 (500 m) [59] to be comparable to Sentinel-2 EVI2.Each MODIS pixel's quality flag was created from the binary MODIS quality (QA) data.
For generating daily EVI2 time-series from Sentinel-2 data and MODIS data, we applied a time-series fitting algorithm, box constrained separable least squares fitting to double logistic model functions [60], so that the datasets could be compared.This is a timeseries fitting method that reduces noise, gap-fills the data, and generates output a regular time step (Figure 3).The inputs of the method were EVI2 time-series and quality flags.Only highest quality data were used in the function fitting.To avoid errors due to unequal time sampling of MODIS [61], the actual acquisition dates of the MODIS observations were used as a time vector of the MODIS EVI2 time-series.The default procedure and settings as described by Jönsson et al. [60] were applied to produce daily EVI2 time-series.The fitting method is similar to the commonly used double logistic model functions [62] but with the following three differences: (1) the box constrained method applies a constant base level, which was estimated from a 5% low threshold of the clear observation histogram for the full time-series; (2) through a shape prior, i.e., a common season estimated for the full time-series through which the box constrained method improves the reconstruction of the seasons also in cases when the amount of observation data is insufficient; and (3) where the observation data was sufficiently large over a growing season, the estimation of parameters of the double logistic function for the season was constrained within a set of ranges ('box') in the box constrained method.
Remote Sens. 2020, 12, x FOR PEER REVIEW 6 For generating daily EVI2 time-series from Sentinel-2 data and MODIS data, w plied a time-series fitting algorithm, box constrained separable least squares fitting to ble logistic model functions [60], so that the datasets could be compared.This is a series fitting method that reduces noise, gap-fills the data, and generates output a re time step (Figure 3).The inputs of the method were EVI2 time-series and quality Only highest quality data were used in the function fitting.To avoid errors due to un time sampling of MODIS [61], the actual acquisition dates of the MODIS observa were used as a time vector of the MODIS EVI2 time-series.The default procedur settings as described by Jönsson et al. [60] were applied to produce daily EVI2 time-s The fitting method is similar to the commonly used double logistic model function but with the following three differences: (1) the box constrained method applies a con base level, which was estimated from a 5% low threshold of the clear observation gram for the full time-series; (2) through a shape prior, i.e., a common season estim for the full time-series through which the box constrained method improves the r struction of the seasons also in cases when the amount of observation data is insuffi and (3) where the observation data was sufficiently large over a growing season, the mation of parameters of the double logistic function for the season was constrained w a set of ranges ('box') in the box constrained method.

Daily EVI2 of Flux Footprint
The pixels' daily EVI2 trajectories, which were generated from the previous se were further used to estimate daily EVI2 over the flux footprint.To compute MODI m EVI2 averages over the flux footprint area, denoted as EVI2250m, a maximum o pixels around the tower were used, and 500 m EVI2 averages (EVI2500m) were comp from the 500 m pixel(s) where the towers are located (Figure 4).Sentinel-2 weighted ages over the flux footprint areas, denoted as EVI210m, were thus computed from 450 pixels depending on the estimated footprint size.Since 10 m Sentinel-2 pixels are smaller than flux footprints, we calculated the footprint EVI2 more precisely by

Daily EVI2 of Flux Footprint
The pixels' daily EVI2 trajectories, which were generated from the previous section, were further used to estimate daily EVI2 over the flux footprint.To compute MODIS 250 m EVI2 averages over the flux footprint area, denoted as EVI2 250m , a maximum of four pixels around the tower were used, and 500 m EVI2 averages (EVI2 500m ) were computed from the 500 m pixel(s) where the towers are located (Figure 4).Sentinel-2 weighted averages over the flux footprint areas, denoted as EVI2 10m , were thus computed from 450-5800 pixels depending on the estimated footprint size.Since 10 m Sentinel-2 pixels are much smaller than flux footprints, we calculated the footprint EVI2 more precisely by where EV I2 10m is the footprint EVI2; n is the number of Sentinel-2 pixels covering the flux footprint; EV I2 i and w i (0-1) are EVI2 of the pixel i and the corresponding contribution weight to the total flux.

Empirical Regression Models for GPP Estimation
For the eight stations, regression models were used to estimate daily GPP (g C m −2 day −1 ) following previously tested model formulations for our study region by Schubert et al. [27], one model for each ICOS site: where the variable 2 is daily averaged EVI2 reconstructed by the box constrained method from Sentinel-2 MSI data and MODIS data. is the EC tower measured daily sum of PPFD (µmol m −2 day −1 ), the daily average air temperature Ta (°C), the daily average VPD (hPa) or the product of them.Available light and temperatures above the freezing point are requirements for photosynthesis.Therefore, only dates with both a daily sum PPFD larger than 0 µmol m −2 day −1 and a daily average Ta above 0 °C were used in the further model calculation.
in the model was smoothed by a moving average with a 7day moving window for reducing the noise of the data.Although there would be six possible model formulas based on the combinations between satellite-derived EVI2 and the different environmental variables at each site, a single model formula was required to enable a direct comparison between Sentinel-2 MSI EVI2 and MODIS EVI2 at each site.The selection of this single model formula was achieved by assessing the coefficient of determination values (R 2 ) and root mean square error (RMSE) between the EC GPP and model estimated GPP.
To evaluate the EVI2-driven GPP models, we computed the root mean square error (RMSE), mean absolute error (MAE) and R 2 between EVI2-driven daily GPP ( ) and EC-driven daily GPP ( ) for each GPP site:

Empirical Regression Models for GPP Estimation
For the eight stations, regression models were used to estimate daily GPP (g C m −2 day −1 ) following previously tested model formulations for our study region by Schubert et al. [27], one model for each ICOS site: where the variable EV I2 is daily averaged EVI2 reconstructed by the box constrained method from Sentinel-2 MSI data and MODIS data.E is the EC tower measured daily sum of PPFD (µmol m −2 day −1 ), the daily average air temperature T a ( • C), the daily average VPD (hPa) or the product of them.Available light and temperatures above the freezing point are requirements for photosynthesis.Therefore, only dates with both a daily sum PPFD larger than 0 µmol m −2 day −1 and a daily average T a above 0 • C were used in the further model calculation.E in the model was smoothed by a moving average with a 7-day moving window for reducing the noise of the data.Although there would be six possible model formulas based on the combinations between satellite-derived EVI2 and the different environmental variables at each site, a single model formula was required to enable a direct comparison between Sentinel-2 MSI EVI2 and MODIS EVI2 at each site.
The selection of this single model formula was achieved by assessing the coefficient of determination values (R 2 ) and root mean square error (RMSE) between the EC GPP and model estimated GPP.
To evaluate the EVI2-driven GPP models, we computed the root mean square error (RMSE), mean absolute error (MAE) and R 2 between EVI2-driven daily GPP (GPP EV I2 ) and EC-driven daily GPP (GPP EC ) for each GPP site: where n is the number of valid days for modelling GPP at the site.In order to further clarify whether the RMSE, MAE, and R 2 at different spatial resolutions (10 m, 250 m and 500 m) are significantly different, we used the Wilcoxon signed-rank test to perform pairwise tests with a null hypothesis of zero medians for the difference of RMSE/MAE/R 2 between paired sample sites.

Relationships between GPP, EVI2 and Environmental Variables
Daily EVI2 from Sentinel-2 MSI and MODIS all showed linear relationships with flux tower-derived GPP (R 2 = 0.45-0.93,Table 2).There was no tendency for stronger relationships with higher resolution data.MODIS EVI2 showed stronger relationships than Sentinel-2 at four sites; however, at the evergreen forest site SE-Htm, the relationships with EC GPP were considerably weaker than for Sentinel-2 (R 2 = 0.45 and 0.56).Comparing the two MODIS resolutions, the higher spatial resolution EVI2 250m provided stronger relationships with flux tower GPP at the agricultural site SE-Lnn.At this site and the deciduous forest site DK-Sor, daily EVI2 at all spatial resolutions from both Sentinel-2 MSI and MODIS showed very strong relationships with flux tower GPP (R 2 > 0.85).The environmental variables showed lower correlations with GPP than EVI2 for five of the eight sites.T a had higher R 2 than PPFD at most sites, expect at the agricultural site SE-Lnn and the evergreen forest site SE-Htm.For PPFD, R 2 values were lower at high latitudes than at low latitudes, and values at the forest sites were higher than values for the agricultural and wetland sites.VPD showed the weakest relationships among all three environmental variables at all sites.In DK-Sor and SE-Nor, however, the relationships between VPD and GPP were relatively strong, with R 2 values of 0.57 and 0.53, respectively.Spring 2017 was relatively dry in the study areas and for SE-Nor summer 2017 remained dry which increased the sensitivity to VPD.Otherwise, shallow soil depth at these sites might affect the water storage capacity.Despite this, the influence of VPD seemed to be low for most of the sites so only EVI2, T a and PPFD were chosen for the final regression models, because of their stronger relationships with GPP.

Multiple Linear Regression Models of GPP
Based on the general relationships shown in Table 2, and with the aim of fitting models that had the same formulation for the three satellite data sets at each site, T a was chosen as input environmental variables into the final linear GPP regression model: We tested three models at all sites, and the RMSEs, MAEs, and R 2 of satellite-derived GPP and flux data-derived GPP were calculated (Table A1).We found that the differences of RMSE, MAE, and R 2 due to the differences in the model were greater than those due to the difference in the satellite data.Referring to the RMSEs from the three models at each site (Table A1), we used the model that generated the highest R 2 at each site for the final GPP estimation results (Table 3).Therefore, we applied the model driven by T a (Equation ( 5)) at all sites.The slopes and intercepts of the regression models are similar within the same site (Table 3).The largest differences of slopes and intercepts among the different satellite EVI2 data were found at the three southern sites, SE-Lnn, SE-Htm, and DK-Sor.RMSEs, MAEs, and R 2 had similar patterns across sites, satellite data sets and model formulations, which means that the low RMSE values usually coincided with low MAE and high R 2 values.There was no single satellite dataset that always generated the lowest RMSE, lowest MAE, or highest R 2 .RMSE, as well as MAE and R 2 , differed mostly across sites, and the smallest differences were found between the EVI2 variables from the different satellite datasets (Table A1).All the p-values shown in Table 4 from Wilcoxon signed-rank tests are over 0.10, which means that the test failed to reject the null hypothesis of zero medians for the difference of RMSE/MAE/R 2 between paired sample sites.Separately averaging RMSE, MAE, and R 2 across sites, it was seen that there was almost no difference in the results generated from Sentinel-2 MSI data and MODIS data (Table 5).
The modelled GPP time-series were very close between Sentinel-2 and 250 m MODIS data at all sites (Figure 5).The shape and magnitude of the daily GPP time-series between satellite driven GPP and flux data-derived GPP were similar (Figure 5).Peak season GPP decreased from the south to the north.The subarctic and boreal mire sites (SE-Sto and SE-Deg, respectively) had much smaller peak GPP (less than 5 g C m −2 d −1 ) compared to other sites.The GPP peak at the coniferous forest sites varied between 10-15 g C m −2 d −1 , whereas the deciduous and agricultural sites had the peak around 15-20 g Cm −2 d −1 .The agriculture site SE-Lnn had high peak GPP but a relatively short growing season with a narrow shape, mainly controlled by the activities of cultivation.The GPP time-series pattern at the deciduous and evergreen forest sites, DK-Sor and SE-Htm, was much wider than at SE-Lnn, which indicates a much longer growing season (Figure 5).The modelled GPP time-series were very close between Sentinel-2 and 250 m MODIS data at all sites (Figure 5).The shape and magnitude of the daily GPP time-series between satellite driven GPP and flux data-derived GPP were similar (Figure 5).Peak season GPP decreased from the south to the north.The subarctic and boreal mire sites (SE-Sto and SE-Deg, respectively) had much smaller peak GPP (less than 5 g C m −2 d −1 ) compared to other sites.The GPP peak at the coniferous forest sites varied between 10-15 g C m −2 d −1 , whereas the deciduous and agricultural sites had the peak around 15-20 g Cm −2 d −1 .The agriculture site SE-Lnn had high peak GPP but a relatively short growing season with a narrow shape, mainly controlled by the activities of cultivation.The GPP time-series pattern at the deciduous and evergreen forest sites, DK-Sor and SE-Htm, was much wider than at SE-Lnn, which indicates a much longer growing season (Figure 5).

Discussion
The results showed that for all sites and all satellite products, the relationship between EVI2 and the flux tower modelled GPP was generally strong (Table 2).This means that seasonal GPP variations can be estimated at different resolutions, from 500 down to 10 m spatial resolution, implying that remotely sensed data has the potential to provide highly detailed spatial estimates of carbon uptake dynamics across northern boreal landscapes.This tool can contribute to future work on carbon-friendly forest management and climate mitigation.
The model selection using climatological flux footprint led to a linear regression model, driven by T a , in combination with satellite reconstructed daily EVI2 (Equation ( 6)).However, the selected model is not always optimal for each site (Table A1).The agricultural sites and deciduous forest sites had a larger seasonal variation of greening in comparison with the evergreen coniferous forest sites, which helps to explain the strong relationship between EVI2 and GPP in SE-Lnn and DK-Sor.Air temperature was the most important environmental variable for estimating GPP at the northern sites (SE-Sto, SE-Svb, SE-Deg and FI-Hyy) whereas at the southern sites (SE-Nor, SE-Lnn, SE-Htm and DK-Sor) the differences between models were small.This is explained by the fact that in the northern part of our study area the growing season begins when the temperature exceeds the temperature limit, which occurs well after the amount of light has reached sufficient levels for photosynthesis [63,64].This finding indicates that, although EVI2 is the most important driving variable overall, choosing the correct model formulation is important for GPP estimation across a study area traversing such a large latitudinal range.
The method of estimating GPP used in this study is similar to the method proposed by Schubert et al. [27] regarding the methodology, model design and study area, but we used finer spatial and temporal resolution satellite data, dynamic flux footprint, and groundbased measured environmental variables.In comparison with the study by Schubert et al. [27], we used four common sites: SE-Sto, FI-Hyy, SE-Nor and DK-Sor.The estimation of GPP has small RMSE and higher R 2 at site SE-Sto, but it has larger RMSE and lower R 2 at site FI-Hyy.The RMSEs are larger but R 2 are lower at sites SE-Nor and DK-Sor.One of the reasons for this inconsistency is that we used a different time period in the study.

Similar Performance of Sentinel-2 and MODIS on GPP Modelling
Despite its finer spatial resolution, Sentinel-2 data did not enhance the performance of the empirical GPP models.This result is mainly explained by the homogeneous surrounding vegetation of the EC areas in combination with the quite large area of flux footprint; however, other factors that may have contributed to this result are, e.g., (1) the function fitting model used eliminated minor fluctuations of EVI2 time-series, which might have masked some differences among the input datasets; (2) the use of linear regression models reduced the magnitude difference between Sentinel-2 EVI2 and MODIS EVI2 data.These points are further discussed below.
A factor that strongly reduced the differences in GPP relationships between Sentinel-2 and MODIS was the homogeneity of the EC sites.We applied the FFP model to more precisely estimate the EC measurements' footprint [35,39].The coverage areas of MODIS pixels and Sentinel-2 pixels that are used to represent the flux footprint are quite different at some sites, e.g., SE-Lnn.However, the footprint EVI2 in Sentinel-2 and MODIS scales did not show large differences due to the homogeneity of the surrounding areas of EC measurement.
As previously mentioned, the satellite daily EVI2 time-series were generated based on double logistic functions, which are robust at handling data gaps and reducing data noise but lead to a loss of daily fluctuations of the EVI2.Therefore, the day-by-day fluctuations in the output GPP were all mainly generated by T a or PPFD.Since the same environmental datasets were applied to the GPP models for all satellite data, the differences in the GPP outputs between the satellite datasets at each site only reflect the differences in fitted EVI2.Thus, the seasonal shape and magnitude of EVI2 make up the principal differences between the GPP of the different satellite datasets (Figure A1 in Appendix A).Our previous experiences show that Sentinel-2 MSI EVI2 time-series perform better than MODIS in capturing ground multispectral measured EVI2 time-series [65].However, in the current study, the effects of the differences in seasonal shape and magnitude of EVI2 time-series on GPP modelling at each site were reduced by the intercepts (a) and slopes (b) in the linear regression models.The intercepts reduced the bias of estimations, while the slopes reduced the differences in the amplitudes of the estimations.
In a recently published study, Wang et al. [11] applied a very high resolution (VHR) hyperspectral sensor on a UAV platform, to study the effects of spatial resolution on the model performance in a willow field at a 30 minutes' flux basis.They found that the performance of footprint weighted simulations to represent the measured flux started to decline as resolutions exceeded 10 m, while performance was high and relatively constant from 0.03 m to 10 m.They compared their results with the typical scales of the canopy structure and concluded that the spatial resolution should be in the same order of magnitude as the typical crown sizes and distances between the plant rows.While this high requirement might be typical for a homogeneous study site as the willow plantation, the rule of thumb might be general [11].Yet the spatial resolutions of satellite-based remote sensing techniques are currently much coarser.

Challenges and Outlook
In this study, we did not separate the data into training data and evaluation data.This was because of the limited number of flux sites and the short coverage time of Sentinel-2 in Northern Europe, and because our main concern was the ability of Sentinel-2 and MODIS to drive the GPP models but not the absolute accuracy of the prediction model.With more data in the future, it will be possible to fine-tune an accurate model by dividing the data into training data and evaluation data.The high p-values in the test indicated that the conclusion of no significant differences between the data sets may change with a larger sample size.
Although the empirical regression models were generally accurate, further GPP model development may need to be carried out when developing upscaling methodology from Sentinel-2 MSI.The fact that model selection had more impact on the accuracies than the choice of input data highlights the importance of further studying the input environmental drivers and the model formulations for estimating GPP, e.g., Tagesson et al. [66].To further improve the match between flux tower EC data and remote sensing, applying daily rather than annual footprint climatologies is a possibility that should be pursued in future studies [38,39].This is a rather complex analysis and given the small differences observed between the different data sets it is not likely that this would radically change the conclusions in the current study.For upscaling across large areas, a practically applicable methodology needs to be developed (such as the abovementioned study by Wang et al. [11]).It may also be fruitful to further investigate the drivers behind the spatial variation of GPP, e.g., canopy openness or leaf area index, and find suitable remote sensing variables that directly or indirectly represent these drivers.Another factor to consider is the contribution of the understory vegetation to the ecosystem carbon exchange, which is not usually taken into account in satellite-derived GPP models.This affects remotely sensed data [67] and could account for errors in GPP estimates at the beginning and the end of the growing season [68].Along with EC towers, flux chambers may be useful for mapping spatial heterogeneity in GPP in low vegetation areas [69].
With regard to satellite data processing, we suggest analyzing more efficient vegetation indices, e.g., the Plant Phenology Index, which has shown to correlate more strongly with GPP than EVI or NDVI [70,71].Though efficient, this index requires angular effects to be corrected for, which precluded us from using the index with the available Sentinel-2 data.

Conclusions
In this paper, we examined the effect of the spatial resolution of daily two-band enhanced vegetation index (EVI2, 10m from Sentinel-2 MSI, 250m and 500m from MODIS) on the accuracy of gross primary productivity (GPP) estimations with different regression models at eight Nordic eddy covariance flux sites.The results gave high and consistent coefficients of explanation (R 2 = 0.84 for Sentinel-2; R 2 = 0.83 for MODIS) verifying the potential of remotely sensed data to provide accurate estimates of the seasonal carbon uptake dynamics across northern boreal landscapes at both the coarse and fine spatial scale.Sentinel-2 opens the possibility of very detailed GPP estimation across the northern landscape, which is of benefit for climate mitigation and adaptation work.The analysis showed that the quantitative differences in GPP estimation due to the different input data were small between Sentinel-2 MSI GPP and MODIS GPP; smaller than the effect of different model formulations.This can be explained by the combined effects of the higher temporal frequency of MODIS compensating for its lower spatial resolution, the time-series fitting method used for reconstructing daily EVI2, the homogeneity of the surrounding areas of the EC measurements, and the use of linear regression models.To further strengthen the development of GPP models based on Sentinel-2 data will require more analyses of different model formulations (including environmental drivers), more tests of remotely sensed indices and biophysical parameters, and analyses across a wider range of geographical locations and time.

Figure 1 .
Figure 1.Study area with location of flux sites and corresponding Sentinel-2 tiles.Figure 1. Study area with location of flux sites and corresponding Sentinel-2 tiles.

Figure 1 .
Figure 1.Study area with location of flux sites and corresponding Sentinel-2 tiles.Figure 1. Study area with location of flux sites and corresponding Sentinel-2 tiles.

Figure 2 .
Figure 2.An example of the contribution contours of eddy covariance flux measurement.The black dot represents the location of the flux tower.The red curves represent the flux contribution of the area.The background image (green) shows Sentinel-2 EVI2 at the flux tower site SE-Lnn on day of year 182 in 2017.The 80% flux footprint contribution contour was used to define the footprint area of the site.

Figure 2 .
Figure 2.An example of the contribution contours of eddy covariance flux measurement.The black dot represents the location of the flux tower.The red curves represent the flux contribution of the area.The background image (green) shows Sentinel-2 EVI2 at the flux tower site SE-Lnn on day of year 182 in 2017.The 80% flux footprint contribution contour was used to define the footprint area of the site.

Figure 3 .
Figure 3. Fitted daily EVI2 time-series on a Sentinel-2 pixel (10 m × 10 m) at site SE-Lnn (Sentinel-2 tile 33VUE; row 889 and column 3190).The black dots represent the raw EVI2 time-series of the pixel.The red and black circles are point marked as vegetation and soil from SCL.The red curve is the fitted daily EVI2 time-series by box constrained separabl least squares fitting to double logistic model functions.
EVI2; is the number of Sentinel-2 pixels coverin flux footprint;2 and (0-1) are EVI2 of the pixel and the corresponding c bution weight to the total flux.

Figure 3 .
Figure 3. Fitted daily EVI2 time-series on a Sentinel-2 pixel (10 m × 10 m) at site SE-Lnn (Sentinel-2 tile 33VUE; row 8891 and column 3190).The black dots represent the raw EVI2 time-series of the pixel.The red and black circles are points marked as vegetation and soil from SCL.The red curve is the fitted daily EVI2 time-series by box constrained separable least squares fitting to double logistic model functions.

Figure 4 .
Figure 4.An example of the corresponding Sentinel-2 and MODIS pixels and the flux footprint site SE-Lnn.Sentinel-2 10 m resolution pixels (red) cover the region that contributed 80% of the total flux in 2017.The flux footprint region is covered by two MODIS 250m pixels (hatched) and one MODIS 500m pixels (black).The background image (green) shows Sentinel-2 EVI2 at the flux tower site SE-Lnn on day of year 182 in 2017.

Figure 4 .
Figure 4.An example of the corresponding Sentinel-2 and MODIS pixels and the flux footprint site SE-Lnn.Sentinel-2 10 m resolution pixels (red) cover the region that contributed 80% of the total flux in 2017.The flux footprint region is covered by two MODIS 250m pixels (hatched) and one MODIS 500m pixels (black).The background image (green) shows Sentinel-2 EVI2 at the flux tower site SE-Lnn on day of year 182 in 2017.

Figure 5 .
Figure 5. Time-series of daily satellite estimated GPP and EC GPP for 2017-2018.The red curve and dashed black line are GPP-derived from 10 m Sentinel-2 EVI2 and from 250 m MODIS EVI2 We have omitted the 500 m GPP curve in the figure for better visualization.

Figure 5 .
Figure 5. Time-series of daily satellite estimated GPP and EC GPP for 2017-2018.The red curve and dashed black line are GPP-derived from 10 m Sentinel-2 EVI2 and from 250 m MODIS EVI2 We have omitted the 500 m GPP curve in the figure for better visualization.

Table A1 .
RMSEs, MAEs, and R 2 between EC driven GPP (g C m −2 day −1 ) and EVI2 driven GPP from three linear models in 10 m, 250 m, and 500 m resolutions at eight Nordic sites.The three linear models are model I: GPP = a•T a •EV I2 + b, model II: GPP = a•PPFD•EV I2 + b, and model III: GPP = a•T a •PPFD•EV I2 + b.

Table 1 .
The eight ICOS flux tower sites used in this study.For more site information, refer to the national ICOS information pages (https://www.icos-ri.eu/).

Table 2 .
Coefficient of determination, R 2 , of environmental variables/satellite EVI2 versus EC derived GPP in daily time step.No VPD data were available at FI-Hyy site.

Table 3 .
Linear regressions for GPP (g C m −2 day −1 ) at each site driven by the products of T a and EVI2 in 10 m from Sentinel-2, and 250 m, 500 m from MODIS.

Table 4 .
Paired p-values of RMSE, MAE and R 2 over 8 sites from Wilcoxon signed-rank test.All p-values of 0.10 over fail to reject the null hypothesis.

Table 5 .
Average errors and coefficients of determination across sites.

Table 5 .
Average errors and coefficients of determination across sites.
Note: the numbers in bold type are the smallest RMSEs, the smallest MAEs, and the largest R 2 in the three models.