Spatial, Phenological, and Inter-Annual Variations of Gross Primary Productivity in the Arctic from 2001 to 2019

: Quantifying the spatial, seasonal (phenological), and inter-annual variations of gross primary productivity (GPP) in the Arctic is critical for comprehending the terrestrial carbon cycle and its feedback to climate warming in this region. Here, we evaluated the accuracy of the MOD17A2H GPP product using the FLUXNET 2015 dataset in the Arctic, then explored the spatial patterns, seasonal variations, and interannual trends of GPP, and investigated the dependence of the spatiotemporal variations in GPP on land cover types, latitude, and elevation from 2001 to 2019. The results showed that MOD17A2H was consistent with in situ measurements (R = 0.8, RMSE = 1.26 g C m − 2 d − 1 ). The functional phenology was also captured by the MOD17A2H product (R = 0.62, RMSE = 9 days) in the Arctic. The spatial variation of the seasonal magnitude of GPP and its interannual trends is partly related to land cover types, peaking in forests and lowest in grasslands. The interannual trend of GPP decreased as the latitude and elevation increased, except for the latitude between 62 ◦ ~66 ◦ N and elevation below 700 m. Our study not only revealed the variation of GPP in the Arctic but also helped to understand the carbon cycle over this region.


Introduction
Climate change is causing permafrost melting [1], shrub cover expansion, growing season lengthening, and consequently, carbon flux changes in the Arctic [2]. Furthermore, the carbon cycle is also influenced by changes in vegetation phenology [3]. GPP, which is considered the biggest carbon flux of terrestrial ecosystems [4], not only plays a vital role in offsetting the concentration of greenhouse gases and mitigating global warming to a certain extent [5] but also builds a bridge between terrestrial and air carbon. In the context of the Arctic, the rate of climate warming is almost twice the global average, a phenomenon known as Arctic amplification [2,[6][7][8][9][10]. Therefore, quantifying the spatial, seasonal (phenological), and inter-annual variations of Arctic GPP is critical for comprehending the carbon cycle and its feedback to climate warming.
Quantifying global or local GPP has received a great deal of attention in recent studies. Utilizing satellite-based near-infrared reflectance (NIRv) as the proxy of GPP and the revised light-use-efficiency model (i.e., EC-LUE model), Wang et al. [11] and Zheng et al. [12] explored the global spatial patterns of GPP with a spatial resolution of 0.05 degrees. However, the annual average estimates of GPP were not consistent during the same period. Wang et al. [11] reported a range of 128.3 ± 4.0 Pg C year −1 while Zheng et al. [12] reported a range of 106.2 ± 2.9 Pg C yr −1 . Some studies have detected the GPP in the Arctic, but most paid attention to specific ecosystems (e.g., streams and moss communities) [13,14] and few efforts [12,15] have been devoted to investigating the specific situation of the

FLUXNET Data
FLUXNET 2015 is the latest version of the FLUXNET dataset. Compared with previous datasets, FLUXNET v.2015 improves the protocols of data quality and the pipeline of data processing [28]. The Net Ecosystem Exchange (NEE) in the FLUXNET 2015 dataset was gap-filled with the marginal distribution sampling (MDS) method [29]. It was then partitioned into Ecosystem Respiration (RECO) and GPP using the daytime fluxes method (_DT) [30] and the nighttime fluxes method (_NT) [29]. The quality flags in FLUXNET 2015 are given values ranging from 0 to 1, indicating the percentage of high quality gap-filled and measured data [12]; 1 represents the highest quality and 0 represents the poorest quality [31].
As shown in Table 1, there are 12 sites in the FLUXNET datasets located in the study area, including 7 wetland sites, 2 forest sites, a grassland site, a shrublands site, and a permanent snow and ice site. This study used GPP (GPP_NT_VUT_REF) estimated from the night-time method with a daily temporal scale. The selection of GPP followed two criteria: (1) the quality flags were larger than 0.5; and (2) the difference between the GPP derived using the night-time method (GPP_NT_VUT_REF) and the day-time method (GPP_DT_VUT_REF) was lower than 50%. After filtering the data based on the two criteria, daily GPP were temporally aggregated to generate the 8-day averaged GPP, matching the temporal resolution of MOD17A2H.

FLUXNET Data
FLUXNET 2015 is the latest version of the FLUXNET dataset. Compared with previous datasets, FLUXNET v.2015 improves the protocols of data quality and the pipeline of data processing [28]. The Net Ecosystem Exchange (NEE) in the FLUXNET 2015 dataset was gap-filled with the marginal distribution sampling (MDS) method [29]. It was then partitioned into Ecosystem Respiration (RECO) and GPP using the daytime fluxes method (_DT) [30] and the nighttime fluxes method (_NT) [29]. The quality flags in FLUXNET 2015 are given values ranging from 0 to 1, indicating the percentage of high quality gap-filled and measured data [12]; 1 represents the highest quality and 0 represents the poorest quality [31].
As shown in Table 1, there are 12 sites in the FLUXNET datasets located in the study area, including 7 wetland sites, 2 forest sites, a grassland site, a shrublands site, and a permanent snow and ice site. This study used GPP (GPP_NT_VUT_REF) estimated Remote Sens. 2021, 13, 2875 4 of 21 from the night-time method with a daily temporal scale. The selection of GPP followed two criteria: (1) the quality flags were larger than 0.5; and (2) the difference between the GPP derived using the night-time method (GPP_NT_VUT_REF) and the day-time method (GPP_DT_VUT_REF) was lower than 50%. After filtering the data based on the two criteria, daily GPP were temporally aggregated to generate the 8-day averaged GPP, matching the temporal resolution of MOD17A2H.

Satellite Data
MOD17A2H (Collection 6) is a standard satellite product with a spatial resolution of 500 m and a temporal resolution of 8 days. It is calculated based on the light use efficiency (LUE) approach by Monteith [31]: where ε, fPAR, and PAR denote the radiation use efficiency coefficient (RUE), the fraction of incident PAR absorbed by the surface, and photosynthetically active radiation, respectively [32]. According to Running et al. [32], the GPP values of MOD17A2H refer to the sum of the GPP during an 8-day period. In this study, we averaged the total GPP to generate the 8-day averaged GPP.

Land Cover
The MCD12Q1 product provides global land cover type data at a spatial resolution of 500 m at an annual time step from 2001 to 2019. It is based on the supervised classification of MODIS reflectance data with six different classification schemes, including the IGBP (Annual International Geosphere-Biosphere Programmer), which was widely utilized due to its high accuracy and widespread acceptance [33]. Thus, the IGBP classification method was utilized in this study. The land cover data from 2001 to 2019 were chosen to produce a spatially continuous dataset via mosaic. The filling data of MCD12Q1 was removed to reduce their effect on the results. Evergreen needleleaf forests, evergreen broadleaf forests, deciduous needleleaf forests, deciduous broadleaf forests, and mixed forests were grouped into forests. Closed shrublands and open shrublands were grouped together as shrublands. Woody savannas and savannas were combined into savannas.

DEM (Digital Elevation Model)
Multi-Error-Removed Improved-Terrain (MERIT) DEM is an improvement of SRTM3 (Shuttle Radar Topography Mission v.3) DEM, with a spatial resolution of 3 arc-second (~90 m). It removes multiple error components from the SRTM3 DEM, including stripe noise caused by the sensor error, speckle noise of surface reflectance, absolute bias derived from the limited control points of the ground, and tree height bias where the canopies were incorrectly classified as the land surface [34,35]. MERIT was chosen because its accuracy is higher than that of SRTM and NASADEM (NASA Digital Elevation Model) [36] and because of the data availability in the Arctic. In order to match the resolution of MOD17A2H, the DEM dataset was resampled to 500 m using the bilinear method.

Accuracy Assessment
Although the validation based on in situ leaves issues of scale unresolved, which might introduce uncertainties to the verification, it is still an important method in regions lacking long-term validation data [37][38][39]. Here, the direct comparison method was utilized because it is simple and easy to implement. To avoid the influence of data noise, geometric mismatch, and spatial heterogeneity on the validation results, the average of the 3 × 3 pixels of MOD17A2H centered around tower coordinates in situ was used to match with in situ as suggested by Ueyama et al. [40].
Statistical indices, including the coefficient of correlation (R), root mean square error (RMSE), and Bias were used to indicate the accuracy of MOD17A2H [41]. R measures the consistency between MOD17A2H and in situ data, while RMSE measures the average absolute error of the MOD17A2H over a single in situ. Bias describes the average deviation between MOD17A2H and in situ measurements: where P l and O l are the MOD17A2H and in situ-based GPP on the l th time period, respectively. P and O are the averaged value of the MOD17A2H and in situ-based GPP time series, respectively. h is the total number of time periods. The assessment was twofold. First, the performances of the MOD17A2H product were assessed separately over each site. Second, by combining the data points of all sites within each specific land cover type, the accuracies of MOD17A2H over different land cover types were assessed and compared.

Comparison of Phenological Patterns between In Situ and MOD17A2H
The phenological patterns in ecosystem GPP are important in the terrestrial carbon cycle and have significant ecological implications. To understand if the MOD17A2H satellite GPP product can capture the functional phenology, which has been defined as the interaction and close association between plant functional traits and phenology [42], of in situ GPP, we first filled the data gap in the original 8-day GPP time series using a linear interpolation method. The gap-filled GPP time series was further smoothed using the Savitzky-Golay (SavGol) filter with a window size of 9 time steps and a second-order polynomial, which not only eliminated noise but also preserved the basic phenological attributes [43]. Finally, we extracted the timing of maximum GPP (day of year, DOY) during the photosynthetically active period for each site from both the in situ and MOD17A2H data. Agreement between the peak timings extracted from the in situ GPP and those extracted from the MOD17A2H was used as an indicator of the performance of the MOD17A2H GPP product in representing the functional phenology patterns of arctic ecosystems.

The Spatial Distribution Characteristics Identification and Trend Detection
The spatial distribution of GPP was explored. First, the pixel-wise multiyear averaged monthly GPP was calculated to check the GPP spatial distribution and the phenology patterns in different months. Second, the annual-maximum and annual-averaged GPP were identified to detect the distribution of GPP in the Arctic.
A pixel-based simple linear regression, in which time is the independent variable and GPP is the dependent variable, was applied to detect the trend of GPP. In addition, the significance of the interannual trend was evaluated utilizing the Mann-Kendall (MK) test [44], and the significant trends (p < 0.025) of the Arctic were retained.
In this study, both MOD17A2H data and the linear regression function were provided by Google Earth Engine (GEE), which is a cloud-based computing platform for planetary-scale data analysis, mapping, and modeling, providing free access to numerous global datasets and advanced computational capabilities [45]. GEE was employed for the following reasons: it provides easy access to the MOD17A2H datasets and other related datasets such as land cover types and elevation; it enables rapid exploration of long time series datasets without downloading them; and it provides a library of functions such as linear regression function, which are applied for data analysis and result display.

Validation MOD17A2H Based on In Situ
The results of in situ and MOD17A2H GPP were compared in different land cover types. Figure 2 shows the time series of MOD17A2H and the in situ-based GPP over wetlands; their scatter plots are presented in Figure 3. As shown in Figure 2, missing data occurs frequently in the time series, especially in winter and early spring. This is attributed to the weak photosynthetic activity of vegetation and the lower data coverage during this period. Both MOD17A2H and in situ-based GPP show reasonable seasonal and annual variability over wetlands ( Figure 2). However, the agreement between them is significantly different from site to site.    The comparison results over wetlands can be divided into two groups according to the performance of MOD17A2H: (1) FI_Lom, US_Atq, GL_ZaF, and SJ_Adv; and (2) GL_NuF, RU_Che and US_Ivo. MOD17A2H was found to underestimate in group (1) and overestimate in group (2). Nevertheless, the degree of underestimation varies significantly from site to site and is low for FI_lom and US_Atq, both in magnitude and temporal variation trend (Figure 2a,c), with RMSEs of 1.12 and 0.69 g C m −2 d −1 and R of 0.9 and 0.86, respectively (Figure 3a,c). For GL_ZaF and SJ_Adv, the agreement is worse, as the RMSEs are 2.01 and 1.79 g C m −2 d −1 , respectively (Figure 3e,g), and the bias is very large, with values of −1.54 and −1.21, respectively. The large discrepancy between MOD17A2H and in situ GPP may be caused by their spatial scale mismatch. Furthermore, it should be noted that the data points (less than 50) over these two sites are very limited. For the second group, MOD17A2H is generally consistent with in situ measurements, with the RMSEs of 0.74, 1.22, and 1.0 g C m −2 d −1 , and R of 0.86, 0.82, and 0.86 over GL_NuF, RU_Che, and US_Ivo, respectively.
As indicated by Figure 3, the MOD17A2H performs well over wetlands when excluding the sites with limited data points (i.e., GL_ZaF and SJ_Adv). One interesting observation was that the latitudes of GL_ZaF and SJ_Adv are higher than 74 • N while those of the other sites are lower than 74 • N. In addition, the GL_NuF site where the best agreement between MOD17A2H and in situ occurs has the lowest latitude of all the wetland sites. Therefore, it can be inferred that the accuracy of MOD17A2H may be related to latitude, as it is higher over low latitudes but lower over high latitudes. When it comes to the overall performance of MOD17A2H over wetlands (Figure 3h), the overall RMSE and R are 1.17 g C m −2 d −1 and 0.76, respectively. MOD17A2H slightly underestimates GPP, with a bias of −0.19. Based on the results above, MOD17A2H can be considered capable of revealing the spatial distribution characteristics and the temporal trend of GPP over wetlands in the Arctic.  (Table 1), indicating that the accuracy of MOD17A2H is also influenced by other factors. As shown in Figure 5c, although MOD17A2H slightly underestimates GPP over forests (Bias = −0.36), overall, MOD17A2H was consistent with site measurements over forests, with RMSE and R of 1.51 g C m −2 d −1 and 0.79, respectively.
.     Figures 6 and 7 present the comparison results over grasslands, shrublands, and permanent snow and ice. MOD17A2H generally underestimates GPP over these land cover types but the degree of underestimation varies with sites, as it is relatively weak for RU_Cok (RMSE = 1.02 g C m −2 d −1 , Bias = −0.44) but strong for GL_ZaH and SJ_Blv (with the RMSEs of 0.61 and 0.19 g C m −2 d −1 and the Bias of −0.5 and −0.17, respectively). It is important to remember that the latitudes of GL_ZaH and SJ_Blv are higher than at RU_Cok, which further demonstrates that the accuracy of MOD17A2H is related to the latitude of the sites.
When combining the data points of all the sites (Figure 8), MOD17A2H slightly underestimates GPP, with the bias of −0.32. The overall accuracy of MOD17A2H is reasonable over the Arctic, with RMSE of 1.26 g C m −2 d −1 and R of 0.8, respectively. These indicators demonstrate that MOD17A2H is able to capture the spatiotemporal variation characteristics of GPP in the Arctic.
From the validation results based on in situ measurements, it is shown that the MOD17A2H generally underestimates GPP over these land cover types. Nevertheless, depending on the location of the sites, it may underestimate or overestimate GPP within each land cover type. This demonstrates that the accuracy of MOD17A2H is also influenced by other factors in addition to land cover types. For instance, the latitude seems to be associated with the accuracy of MOD17A2H given that the accuracy of MOD17A2H tends to be higher over low latitudes but lower over high latitudes (>74 • N), which might be due to the actual maximum radiation conversion efficiency (ε max ) of vegetations being quite different from the given ε max in high latitude. In fact, the misclassification of a pixel is also responsible for the inconsistency of the accuracies derived from different sites within the same land cover type. Because the classification scheme adopted by IGBP is too general, it cannot reveal the detailed categories carrying out photosynthesis. Moreover, MOD17A2H is calculated based on the concept of light use efficiency (LUE), which assumes a fixed maximum radiation conversion efficiency (ε max ) of each land cover type [19]. This treatment also introduced errors when misclassification occurred. Another cause of the inconsistency is the different degrees of surface heterogeneities within the satellite pixel, which cause the sites to be more or less representative of the satellite pixel. Last but not least, since missing data of in situ may occur unequally during each 8-day period, the temporal representativeness varies across these sites. Despite these uncertainties, the validation results still suggest that the performance of MOD17A2H is better over shrublands and wetlands than it is over forests.
Remote Sens. 2021, 13, x FOR PEER REVIEW 10 of 22 Figures 6 and 7 present the comparison results over grasslands, shrublands, and permanent snow and ice. MOD17A2H generally underestimates GPP over these land cover types but the degree of underestimation varies with sites, as it is relatively weak for RU_Cok (RMSE = 1.02 g C m −2 d −1 , Bias = −0.44) but strong for GL_ZaH and SJ_Blv (with the RMSEs of 0.61 and 0.19 g C m −2 d −1 and the Bias of −0.5 and −0.17, respectively). It is important to remember that the latitudes of GL_ZaH and SJ_Blv are higher than at RU_Cok, which further demonstrates that the accuracy of MOD17A2H is related to the latitude of the sites.   Figures 6 and 7 present the comparison results over grasslands, shrublands, and permanent snow and ice. MOD17A2H generally underestimates GPP over these land cover types but the degree of underestimation varies with sites, as it is relatively weak for RU_Cok (RMSE = 1.02 g C m −2 d −1 , Bias = −0.44) but strong for GL_ZaH and SJ_Blv (with the RMSEs of 0.61 and 0.19 g C m −2 d −1 and the Bias of −0.5 and −0.17, respectively). It is important to remember that the latitudes of GL_ZaH and SJ_Blv are higher than at RU_Cok, which further demonstrates that the accuracy of MOD17A2H is related to the latitude of the sites.   When combining the data points of all the sites (Figure 8), MOD17A2H slightly underestimates GPP, with the bias of −0.32. The overall accuracy of MOD17A2H is reasonable over the Arctic, with RMSE of 1.26 g C m −2 d −1 and R of 0.8, respectively. These indicators demonstrate that MOD17A2H is able to capture the spatiotemporal variation characteristics of GPP in the Arctic. From the validation results based on in situ measurements, it is shown that the MOD17A2H generally underestimates GPP over these land cover types. Nevertheless, depending on the location of the sites, it may underestimate or overestimate GPP within each land cover type. This demonstrates that the accuracy of MOD17A2H is also influenced by other factors in addition to land cover types. For instance, the latitude seems to be associated with the accuracy of MOD17A2H given that the accuracy of MOD17A2H tends to be higher over low latitudes but lower over high latitudes (>74° N), which might be due to the actual maximum radiation conversion efficiency (εmax) of vegetations being quite different from the given εmax in high latitude. In fact, the misclassification of a pixel is also responsible for the inconsistency of the accuracies derived from different sites within the same land cover type. Because the classification scheme adopted by IGBP is too general, it cannot reveal the detailed categories carrying out photosynthesis. Moreover, MOD17A2H is calculated based on the concept of light use efficiency (LUE), which assumes a fixed maximum radiation conversion efficiency (εmax) of each land cover type [19]. This treatment also introduced errors when misclassification occurred. Another cause of the inconsistency is the different degrees of surface heterogeneities within the satellite pixel, which cause the sites to be more or less representative of the satellite pixel. Last but not least, since missing data of in situ may occur unequally during each 8-day period, the temporal representativeness varies across these sites. Despite these uncertainties, the validation results still suggest that the performance of MOD17A2H is better over shrublands and wetlands than it is over forests.

Evaluation of the Phenological Characteristics of MOD17A2H
According to Figure 9, MOD17A2H does express the phenological characteristics of GPP in the Arctic, with cross-sites R and RMSE of 0.62 and 8.9 days, respectively. In particular, as shown in Figure 9c, the peak GPP timing DOY of grasslands was the most consistent with that extracted from in situ data (R = 0.86, RMSE = 4 days). By contrast,

Evaluation of the Phenological Characteristics of MOD17A2H
According to Figure 9, MOD17A2H does express the phenological characteristics of GPP in the Arctic, with cross-sites R and RMSE of 0.62 and 8.9 days, respectively. In particular, as shown in Figure 9c, the peak GPP timing DOY of grasslands was the most consistent with that extracted from in situ data (R = 0.86, RMSE = 4 days). By contrast, the discrepancy between MOD17A2H and in situ peak GPP timing over forest sites was the highest (Figure 9b), with R and RMSE of 0.52 and 11.9 days, respectively, which demonstrates that the characteristics of the forests (mostly evergreen needleleaf) GPP, derived from MOD17A2H, are difficult to capture compared with the characteristics of other land cover types. This might be due to two reasons: (1) forests are composed of evergreen needleleaf forests, evergreen broadleaf forests, deciduous needleleaf forests, deciduous broadleaf forests, and mixed forests; therefore, the diversity of forests makes phenological characteristics difficult to capture by satellite observations; and (2) evergreen needleleaf forest is one of the major forest types in the Arctic, which is less sensitive to climate changes according to [46], thus the phenological characteristics of evergreen needleleaf forests are hard to capture. Figure 10 shows the annual maximum ( Figure 10a) and annual-averaged GPP (Figure 10b), derived from MOD17A2H, in the Arctic on a pixel basis. The two metrics generally present similar spatial distribution patterns. GPP is relatively low in the northeast of Canada and the regions surrounding Greenland, with an annual-maximum range of 0 to 1.200 g C m −2 d −1 and an annual-averaged range of 0 to 0.900 g C m −2 d −1 . The low GPP over these areas can be explained by the fact that these areas are almost completely covered by grassland and barren land, which have lower GPP (Figure 11c). Therefore, it can be inferred that the spatial distribution of GPP is related to land cover types. Furthermore, the GPP shows a decreasing trend as the latitude increases. This is especially true over the eastern hemisphere of the Arctic, where the annual maximum of GPP drops from 3.300 to 0.300 g C m −2 d −1 and the annual-averaged of GPP decreases from 2.400 to 0.300 g C m −2 d −1 .

Spatial Distribution of Annual-Averaged GPP
other land cover types. This might be due to two reasons: (1) forests are composed of evergreen needleleaf forests, evergreen broadleaf forests, deciduous needleleaf forests, deciduous broadleaf forests, and mixed forests; therefore, the diversity of forests makes phenological characteristics difficult to capture by satellite observations; and (2) evergreen needleleaf forest is one of the major forest types in the Arctic, which is less sensitive to climate changes according to [46], thus the phenological characteristics of evergreen needleleaf forests are hard to capture.  (a-e) Represent permanent wetlands, forests, grasslands, shrublands, and all sites combined, respectively. The gray belt refers to the confidence interval of 95%. Figure 11a,b display the distribution of GPP with latitude and elevation, respectively. It can be seen that GPP generally shows a decreasing trend with latitude, which is in line with the results of Gounand et al. [25]. Nevertheless, the sensitivity of GPP to latitude depends on the situation. The results define three groups of latitudes: (1) latitudes less than 62 • N and more than 80 • N; (2) latitudes higher than 62 • N but lower than 66 • N; and (3) latitudes between 66 • N and 80 • N. In the first group, the GPP distribution is sparse because the region located within this scope is quite limited. Thus, GPP shows irregular variation patterns with latitude. In the second group, GPP is relatively stable, indicating the insensitivity of GPP to latitude within this scope. Nevertheless, GPP presents a clear decreasing trend in the third group, demonstrating that GPP is most sensitive to latitude from 66 • N to 80 • N. Similar to latitude, GPP generally shows a decreasing trend with increasing elevation (Figure 11b). However, the decreasing rates are different depending on the elevation, which is smaller for elevations lower than 700 m, but larger for elevations higher than 700 m.
of Canada and the regions surrounding Greenland, with an annual-maximum range of 0 to 1.200 g C m −2 d −1 and an annual-averaged range of 0 to 0.900 g C m −2 d −1 . The low GPP over these areas can be explained by the fact that these areas are almost completely covered by grassland and barren land, which have lower GPP (Figure 11c). Therefore, it can be inferred that the spatial distribution of GPP is related to land cover types. Furthermore, the GPP shows a decreasing trend as the latitude increases. This is especially true over the eastern hemisphere of the Arctic, where the annual maximum of GPP drops from 3.300 to 0.300 g C m −2 d −1 and the annual-averaged of GPP decreases from 2.400 to 0.300 g C m −2 d −1 .  Figure 11a,b display the distribution of GPP with latitude and elevation, respectively. It can be seen that GPP generally shows a decreasing trend with latitude, which is in line with the results of Gounand et al. [25]. Nevertheless, the sensitivity of GPP to latitude depends on the situation. The results define three groups of latitudes: (1) latitudes less than 62° N and more than 80° N; (2) latitudes higher than 62° N but lower than 66° N; and (3) latitudes between 66° N and 80° N. In the first group, the GPP distribution is sparse because the region located within this scope is quite limited. Thus, GPP shows irregular variation patterns with latitude. In the second group, GPP is relatively stable, indicating the insensitivity of GPP to latitude within this scope. Nevertheless, GPP presents a clear decreasing trend in the third group, demonstrating that GPP is most sensitive to latitude from 66° N to 80° N. Similar to latitude, GPP generally shows a decreasing trend with increasing elevation (Figure 11b). However, the decreasing rates are different depending on the elevation, which is smaller for elevations lower than 700 m, but larger for elevations higher than 700 m. Figure 11c presents the annual-averaged GPP distribution over different land cover types for the entire study period. It can be seen that all land cover types show consider-    Figure 11c presents the annual-averaged GPP distribution over different land cover types for the entire study period. It can be seen that all land cover types show considerable interannual variation except for permanent snow and ice as well as barren land (as indicated by the wide distribution of boxplots). Figure 12 displays the spatial pattern of multiyear (2001-2019) averaged monthly GPP. The phenological cycle of vegetation is clearly shown in the figure. GPP is very low from November to March, with values close to 0. This is because the photosynthesis of vegetation was restricted due to the extremely harsh environment and limited lighting hours. From April to July, with rising temperatures, melting snow and sea ice, and increasing hours of light, the carbon fixation ability of vegetation is stronger, which contributes to the gradual expansion of GPP from the northwest of Canada and the low latitudes of Russia to the entire Arctic. In addition, the mean GPP over the whole Arctic shows a rapid increase during this period, with the values increasing from 0.0859 g C m −2 d −1 in April to 3.739 g C m −2 d −1 in July. However, the GPP gradually decreases from August to October as the mean values decrease from 2.215 to 0.0453 g C m −2 d −1 . This may be attributed to the decrease in temperature, the shortening of day-length, and the senescence of vegetation during this period. The northwestern region of Canada, a small low latitude portion of Russia, and the regions surrounding Iceland have the longest growth period since they are the first to begin and the last to stop photosynthesis. Figure 12 demonstrates that spatial heterogeneity is small during the dormant months (from November to March) and the mid and late summer months (July and August) when GPP is consistently small or large, but large during the transition months (April, May, September, and October) and early summer (June). There are two main reasons for this: the first is the spatial difference in land cover types, the second is the temporal difference in climate conditions [47].

Variation of Monthly GPP
The GPP in the Arctic has a distinct seasonality with the greatest values in July ( Figure 13). Throughout the year, most land cover types follow the general seasonality of GPP; they are lowest from January to March, begin to increase in April and reach a maximum in July, decrease from then on and fall back to the lowest values in November. Forest and savannas present the largest GPP from April to August, followed by shrublands, water bodies, permanent wetlands, and grasslands ( Figure 13). Permanent snow and ice as well as barren land ranks last. Nevertheless, from September to October, water bodies show slightly larger GPP than other land cover types. These results demonstrate that forests and savannas have a rather high carbon storage capacity during the growing season (from April to August) but water bodies are the biggest contributors to carbon fixation from September to October.
The results of Figure 13 are certainly not anticipated because land cover types such as permanent snow and ice, barren land, and water bodies, which cannot carry out photosynthesis, show considerable GPP. This can be explained by the definition of land cover type of IGBP, which is determined by the dominant land cover of a pixel. Water bodies refer to those pixels that are at least 60% covered by permanent water bodies; barren denotes that at least 60% of the pixel is non-vegetated barren (sand, rock, and soil) areas with less than 10% vegetation; permanent Snow and ice means that at least 60% of the pixel is covered by snow and ice for at least 10 months of the year. Therefore, the misclassified part of a pixel is the source of GPP for these three land cover types. From the results above, it is inferred that only extracting the vegetated pixels, as classified by the land cover, will lead to errors in calculating carbon storage in the Arctic. Figure 13 also shows that the monthly GPP has the widest distribution in June, indicating that the interannual variation of GPP is the most significant in June. This is understandable since vegetation grew at the fastest speed from May to June (indicated by the largest slope) and thus is more affected by climate change. Remote Sens. 2021, 13, x FOR PEER REVIEW 15 of 22 The GPP in the Arctic has a distinct seasonality with the greatest values in July ( Figure 13). Throughout the year, most land cover types follow the general seasonality of GPP; they are lowest from January to March, begin to increase in April and reach a maximum in July, decrease from then on and fall back to the lowest values in November. Forest and savannas present the largest GPP from April to August, followed by shrublands, water bodies, permanent wetlands, and grasslands ( Figure 13). Permanent snow and ice as well as barren land ranks last. Nevertheless, from September to October, water bodies show slightly larger GPP than other land cover types. These results demonstrate that forests and savannas have a rather high carbon storage capacity during the growing classified part of a pixel is the source of GPP for these three land cover types. From the results above, it is inferred that only extracting the vegetated pixels, as classified by the land cover, will lead to errors in calculating carbon storage in the Arctic. Figure 13 also shows that the monthly GPP has the widest distribution in June, indicating that the interannual variation of GPP is the most significant in June. This is understandable since vegetation grew at the fastest speed from May to June (indicated by the largest slope) and thus is more affected by climate change. Figure 13. Multiyear averaged monthly GPP distribution over the entire study period for different land cover types (denoted by colored lines) over the Arctic. The black boxplots denote the multiyear monthly GPP distribution of all land cover types.

Trend Estimates of GPP
The interannual variation trend of GPP is shown on a pixel basis ( Figure 14); almost half of the Arctic presents significant positive trends, but the magnitude of trends shows distinct spatial variation. In the northwest of Canada and the latitude lower than 70° N of Russia, the trend of GPP varies significantly, ranging from 0.005 to 0.08 g C m −2 year −1 . By contrast, in the northeast of Canada and the regions surrounding Greenland, the GPP trend varies slightly, ranging from 0 to 0.01 g C m −2 year −1 . The small range of the GPP trend is likely to be associated with grassland and barren land. The former has a relatively small spatial variation, as indicated by the centralized distribution of the interannual trends in GPP. The latter has a very small GPP, which is almost constant over time (Figure 11c).

Trend Estimates of GPP
The interannual variation trend of GPP is shown on a pixel basis ( Figure 14); almost half of the Arctic presents significant positive trends, but the magnitude of trends shows distinct spatial variation. In the northwest of Canada and the latitude lower than 70 • N of Russia, the trend of GPP varies significantly, ranging from 0.005 to 0.08 g C m −2 year −1 . By contrast, in the northeast of Canada and the regions surrounding Greenland, the GPP trend varies slightly, ranging from 0 to 0.01 g C m −2 year −1 . The small range of the GPP trend is likely to be associated with grassland and barren land. The former has a relatively small spatial variation, as indicated by the centralized distribution of the interannual trends in GPP. The latter has a very small GPP, which is almost constant over time (Figure 11c). As seen from Figure 12, there is almost no vegetation productivity in the Arctic from November to March. Thus, Figure 15 only presents the interannual trends from April to October. It is clear that the interannual trends show significant differences between different months. From April to June, the interannual trends in northwestern Canada gradually increase and reach a maximum in June, with an overall trend of 0.068 g C m −2 year −1 ( Figure  15c). From then on, the interannual trend decreases gradually until October (Figure 15d-g). As seen from Figure 12, there is almost no vegetation productivity in the Arctic from November to March. Thus, Figure 15 only presents the interannual trends from April to October. It is clear that the interannual trends show significant differences between different months. From April to June, the interannual trends in northwestern Canada gradually increase and reach a maximum in June, with an overall trend of 0.068 g C m −2 year −1 (Figure 15c). From then on, the interannual trend decreases gradually until October (Figure 15d-g). Furthermore, we also find that the interannual trends in July and August are more spatially heterogeneous than in other months. These results demonstrate that the response of vegetation to climate change is not consistent between different months and over different areas. It is interesting to find that the interannual trend is the most significant in June (Figure 15c) and shows a similar spatial pattern to the overall interannual trend ( Figure 14). This demonstrates that the interannual variations of GPP in the Arctic may be dominated by the change of vegetation productivity in June. The distributions of the interannual trends in GPP with latitude and elevation are shown in Figure 16a,b. The interannual trends first increase at the latitude of 51° N and reach a maximum of 0.018 g C m −2 year −1 at 57° N. From then on, the interannual trend decreases significantly until 62° N. However, the interannual trend seems to be independent of latitude from 62° N to 66° N. Then a clear decreasing trend can be observed from 66° N until 80° N. Therefore, we can conclude that the interannual trend of GPP is sensitive to the latitude, except in the regions located between 62° N and 66° N, which is similar to the distribution of annual-averaged GPP with the variation of latitude. This is mainly because the regions located between (62° N, 66° N) are in the same climatic zone, namely in the north temperate zone. Figure 16b demonstrates that the interannual trend of GPP generally shows a decreasing trend with increasing elevations. However, the sensitivity is relatively weak at low elevations (<700 m) and significant at larger eleva- Figure 15. Spatial distribution of the interannual variation trend (g C m −2 year −1 ) of monthly GPP over the Arctic, (a-g) refers to April to October.
The distributions of the interannual trends in GPP with latitude and elevation are shown in Figure 16a,b. The interannual trends first increase at the latitude of 51 • N and reach a maximum of 0.018 g C m −2 year −1 at 57 • N. From then on, the interannual trend decreases significantly until 62 • N. However, the interannual trend seems to be independent of latitude from 62 • N to 66 • N. Then a clear decreasing trend can be observed from 66 • N until 80 • N. Therefore, we can conclude that the interannual trend of GPP is sensitive to the latitude, except in the regions located between 62 • N and 66 • N, which is similar to the distribution of annual-averaged GPP with the variation of latitude. This is mainly because the regions located between (62 • N, 66 • N) are in the same climatic zone, namely in the north temperate zone. Figure 16b demonstrates that the interannual trend of GPP generally shows a decreasing trend with increasing elevations. However, the sensitivity is relatively weak at low elevations (<700 m) and significant at larger elevations (>700 m). The low GPP and interannual trend in the regions with relatively higher altitude and latitude Figure 11a,b and Figure 16a,b mainly due to the low precipitation and temperature, which are not conducive to plant growth [48,49]. These results could help us understand the latitude or elevation range in which the change of GPP mainly occurred, providing a basis for understanding the changes of the Arctic GPP.
Remote Sens. 2021, 13, x FOR PEER REVIEW 19 of 22 Figure 16. Interannual variation trend of GPP (g C m −2 year −1 ) distribution with latitude (a), elevation (m) (b) over the Arctic. The gray belt refers to the confidence interval of 95%, the blue line refers to the fit line, and the boxplots denote the distribution of the annual-averaged GPP with the latitude. (c) Land cover-dependent interannual variation trend of GPP (g C m −2 year −1 ) over the Arctic. FOR, SHR, SAV, GRA, WET, SNO, and BAR denote the forests, shrublands, savannas, grasslands, permanent wetlands, permanent snow and ice, barren, and water bodies, respectively. Figure 16c displays the boxplots of the interannual trends by combining the pixels for each land cover type. It can be seen that almost each land cover type presents positive trends. However, their magnitude, as well as their spatial variation, depend on the land cover types. Forests have the largest interannual trend, followed by savannas and shrublands. The interannual trend of grassland is not as large as we expected and is even slightly smaller than the trend of permanent wetlands. Savannas and forests show the largest spatial variations in interannual trends, as indicated by the most widespread distribution of boxplots. Permanent wetlands and grasslands present the smallest spatial variations when excluding permanent snow and ice as well as barren land. These results demonstrate that the interannual trends of savannas and forests are more influenced by other factors, while those of permanent wetlands and grasslands are less sensitive to other factors.

Conclusions
Arctic ecosystems have undergone great changes in the context of climate change. GPP is one of the most crucial indicators of the response of ecosystems to climate change. However, few efforts have been devoted to exploring the spatial variation and phenological characteristics of GPP in the Arctic. In response to this challenge, this study investigated the spatial distribution as well as the seasonal (phenological) and interannual variations of GPP in the Arctic using MOD17A2H. Furthermore, the GPP variation trends with land cover types, latitude, and elevation were also explored. In order to en- Figure 16. Interannual variation trend of GPP (g C m −2 year −1 ) distribution with latitude (a), elevation (m) (b) over the Arctic. The gray belt refers to the confidence interval of 95%, the blue line refers to the fit line, and the boxplots denote the distribution of the annual-averaged GPP with the latitude. (c) Land cover-dependent interannual variation trend of GPP (g C m −2 year −1 ) over the Arctic. FOR, SHR, SAV, GRA, WET, SNO, and BAR denote the forests, shrublands, savannas, grasslands, permanent wetlands, permanent snow and ice, barren, and water bodies, respectively. Figure 16c displays the boxplots of the interannual trends by combining the pixels for each land cover type. It can be seen that almost each land cover type presents positive trends. However, their magnitude, as well as their spatial variation, depend on the land cover types. Forests have the largest interannual trend, followed by savannas and shrublands. The interannual trend of grassland is not as large as we expected and is even slightly smaller than the trend of permanent wetlands. Savannas and forests show the largest spatial variations in interannual trends, as indicated by the most widespread distribution of boxplots. Permanent wetlands and grasslands present the smallest spatial variations when excluding permanent snow and ice as well as barren land. These results demonstrate that the interannual trends of savannas and forests are more influenced by other factors, while those of permanent wetlands and grasslands are less sensitive to other factors.

Conclusions
Arctic ecosystems have undergone great changes in the context of climate change. GPP is one of the most crucial indicators of the response of ecosystems to climate change. However, few efforts have been devoted to exploring the spatial variation and phenological characteristics of GPP in the Arctic. In response to this challenge, this study investigated the spatial distribution as well as the seasonal (phenological) and interannual variations of GPP in the Arctic using MOD17A2H. Furthermore, the GPP variation trends with land cover types, latitude, and elevation were also explored. In order to ensure that the results were reliable, the accuracy of MOD17A2H was first evaluated using in situ measurements from FLUXNET.
This study found that MOD17A2H generally underestimates GPP over the land cover types investigated in this study, and its accuracy tends to be higher over low latitude but lower over high latitude. However, the overall accuracy suggests that MOD17A2H is consistent with the FLUXNET 2015 dataset (RMSE = 1.26 g C m −2 d −1 , R = 0.8, Bias = −0.32), and MOD17A2H can represent the phenological characteristics of GPP (RMSE = 8.9 days, R = 0.62). Based on MOD17A2H, it was demonstrated that the maximum GPP occurred in July. In addition, the spatial distribution of GPP is related to land cover types; for example, forests and savannas have relatively high carbon storage capacity from April to August. By comparing the GPP variation with latitude and elevation, it was shown that GPP generally decreases as the latitude and elevation increase. However, the phenomenon is not evident for latitudes in the range (62 • N, 66 • N) and elevation lower than 700 m. The overall trend of GPP in the Arctic is greater than zero and is dominated by the variation of vegetation productivity in June. Furthermore, the response to climate change is different across these land cover types; for example, forests are most sensitive to climate warming. The distribution of the interannual trend in GPP across latitudes and elevations is consistent with the changes in GPP as a function of latitude and elevation.
This study is helpful for understanding the spatiotemporal distribution characteristics of GPP over the Arctic as well as the response of ecosystems to climate change. Nevertheless, the results need to be validated with different satellite products. The number of sites used for validation is limited. Therefore, the presented conclusion about the accuracy of MOD17A2H may not be transferable to other regions. Another limitation is that only individual factors such as land cover type, latitude, and elevation were considered in this paper. Other factors, such as air temperature, precipitation, and snow, that are related to vegetation growth status need to be explored further. This will also be our focus in the future.

Data Availability Statement:
Publicly available datasets were analyzed in this study. This data can be found in Google earth engine platform (https://earthengine.google.com/, accessed on 7 February 2021).