*3.1. 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–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:

$$R = \sum\_{l=1}^{h} \left( P\_l - \overline{P} \right) \left( O\_l - \overline{O} \right) / \sqrt{\sum\_{l=1}^{h} \left( P\_l - \overline{P} \right)^2 \sum\_{l=1}^{h} \left( O\_l - \overline{O} \right)^2},\tag{2}$$

$$\text{RMSE} = \sqrt{\sum\_{l=1}^{h} \left(P\_l - O\_l\right)^2 / h} \,\text{/}h\tag{3}$$

$$\text{Bias} = \sum\_{l=1}^{h} (P\_l - O\_l) / h \tag{4}$$

where *Pl* and *Ol* 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.

#### *3.2. 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.

#### *3.3. 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.

#### **4. Results and Discussion**

*4.1. 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 agreemen<sup>t</sup> between them is significantly different from site to site.

**Figure 2.** Time series comparison between site-based GPP and MOD17A2H over wetland sites, (**<sup>a</sup>**–**g**) refers to FI\_Lom, GL\_NuF, US\_Atq, RU\_Che, GL\_ZaF, US\_Ivo and SJ\_Adv, respectively.

**Figure 3.** Scatter plots of site-based GPP against MOD17A2H over different wetland sites (**<sup>a</sup>**–**g**) and the overall scatter plots by combining all wetland sites (**h**). Flux tower GPP means in situ GPP. The gray belt refers to the confidence interval of 95%.

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.69gCm−<sup>2</sup> d−<sup>1</sup> and R of 0.9 and 0.86, respectively (Figure 3a,c). For GL\_ZaF and SJ\_Adv, the agreemen<sup>t</sup> is worse, as the RMSEs are 2.01 and 1.79 g C m<sup>−</sup><sup>2</sup> 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<sup>−</sup><sup>2</sup> 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 agreemen<sup>t</sup> 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<sup>−</sup><sup>2</sup> d−<sup>1</sup> 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.

Figures 4 and 5 show the evaluation results of MOD17A2H over forest sites (i.e., FI\_Sod and US\_Prr). MOD17A2H underestimates GPP at FI\_Sod but overestimates GPP at US\_Prr. The extent of misestimation differs between the two sites and is weak for FI\_Sod but strong for US\_Prr. Over FI\_Sod, MOD17A2H agrees well with in situ measurements (Figure 4a), with the RMSE of 1.33 g C m<sup>−</sup><sup>2</sup> d−<sup>1</sup> and bias of −0.77 (Figure 5a). However, for US\_Prr, their agreemen<sup>t</sup> is worse (Figure 4b), with the RMSE of 2.05gCm−<sup>2</sup> d−<sup>1</sup> and bias of 1.19 (Figure 5b). It is important to note that the extent of misestimating MOD17A2H varies from year to year. This is especially true over US\_Prr, where the overestimation of MOD17A2H is more significant from 2012 to 2014 (Figure 4b). Although the two sites feature forests, their locations are different (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<sup>−</sup><sup>2</sup> d−<sup>1</sup> and 0.79, respectively.

**Figure 4.** Time series comparison between site-based GPP and MOD17A2H over forest sites, (**a**) and (**b**) refers to FI\_Sod and US\_Prr.

**Figure 5.** Scatter plots of site-based GPP against MOD17A2H over different forest sites (**<sup>a</sup>**,**b**), and the overall scatter plots by combining all forest sites (**c**). Flux tower GPP means in situ GPP. Gray belt refers to the confidence interval of 95%.

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.02gCm−<sup>2</sup> d−1, Bias = −0.44) but strong for GL\_ZaH and SJ\_Blv (with the RMSEs of 0.61 and 0.19gCm−<sup>2</sup> d−<sup>1</sup> 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.26gCm−<sup>2</sup> d−<sup>1</sup> 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 sugges<sup>t</sup> that the performance of MOD17A2H is better over shrublands and wetlands than it is over forests.

**Figure 6.** Time series comparison between site-based GPP and MOD17A2H over grasslands (**a**), shrublands (**b**), and permanent snow and ice (**c**) sites.

**Figure 7.** Scatter plots of site-based GPP against MOD17A2H over grasslands (**a**), shrublands (**b**), permanent snow and ice (**c**) sites. Flux tower GPP means in situ GPP. Gray belt refers to the confidence interval of 95%.

**Figure 8.** The overall scatterplots of in situ based GPP against MOD17A2H. Flux tower GPP means in situ GPP. Gray belt refers to the confidence interval of 95%.

#### *4.2. 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.

#### *4.3. Spatial Distribution Characteristics*

#### 4.3.1. Spatial Distribution of Annual-Averaged GPP

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.200gCm−<sup>2</sup> d−<sup>1</sup> and an annual-averaged range of 0 to 0.900 g C m<sup>−</sup><sup>2</sup> 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<sup>−</sup><sup>2</sup> d−<sup>1</sup> and the annual-averaged of GPP decreases from 2.400 to 0.300 g C m<sup>−</sup><sup>2</sup> d−1.

**Figure 9.** The DOY comparison of the timing of the maximum GPP between the sites and MOD17A2H. (**<sup>a</sup>**–**<sup>e</sup>**) 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.

**Figure 10.** The spatial distribution characteristics of annual-maximum (**a**) and annual-averaged GPP (**b**) over the Arctic. Units of GPP aregCm−<sup>2</sup> d−1.

**Figure 11.** The annual-averaged GPP distribution with latitude (**a**) and elevation (m) (**b**) for the whole 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**) The annual-averaged GPP distribution for the entire study period 2001–2019 over different land cover types. FOR, SHR, SAV, GRA, WET, SNO, BAR, and WAT denote the forests, shrublands, savannas, grasslands, permanent wetlands, permanent snow and ice, barren, and water bodies, respectively.

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).

#### 4.3.2. Variation of Monthly GPP

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<sup>−</sup><sup>2</sup> d−<sup>1</sup> in April to 3.739 g C m<sup>−</sup><sup>2</sup> d−<sup>1</sup> in July. However, the GPP gradually decreases from August to October as the mean values decrease from 2.215 to 0.0453gCm−<sup>2</sup> 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].

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.

**Figure 12.** Spatial distribution characteristics of the multiyear averaged monthly GPP over the Arctic, (**<sup>a</sup>**–**l**) refers to January to December. The mean values of the whole Arctic are also shown in the figure for each month. Units of GPP aregCm−<sup>2</sup> d−1.

**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.

#### *4.4. 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<sup>−</sup><sup>2</sup> year<sup>−</sup>1. By contrast, in the northeast of Canada and the regions surrounding Greenland, the GPP trend varies slightly, ranging from 0 to 0.01gCm−<sup>2</sup> year<sup>−</sup>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).

**Figure 14.** The spatial distribution of interannual trend of GPP (g C m<sup>−</sup><sup>2</sup> year<sup>−</sup>1) over the Arctic (at 97.5% confidence level based on MK test).

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<sup>−</sup><sup>2</sup> year<sup>−</sup><sup>1</sup> (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.

**Figure 15.** Spatial distribution of the interannual variation trend (g C m<sup>−</sup><sup>2</sup> year<sup>−</sup>1) of monthly GPP over the Arctic, (**<sup>a</sup>**–**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<sup>−</sup><sup>2</sup> year<sup>−</sup><sup>1</sup> 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.

**Figure 16.** Interannual variation trend of GPP (g C m<sup>−</sup><sup>2</sup> year<sup>−</sup>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<sup>−</sup><sup>2</sup> year<sup>−</sup>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.
