Comparative evaluation of vegetation greenness trends over circumpolar Arctic tundra using multi-sensors satellite datasets

ABSTRACT

The circumpolar arctic tundra, located at Earth’s highest latitudes, is extremely sensitive to climate warming. Studies on arctic greening, based on satellite data and field measurements, show discrepancies due to differences in spatial resolution across datasets (e.g., Landsat 30-m, MODIS 250-m, and AVHRR GIMMS 8 km). Research on scale effects has been limited, mostly focusing on small areas rather than the entire 7.11 million km² arctic tundra. Our study addresses this by mapping scale effects across the entire tundra using Normalized Difference Vegetation Index (NDVI) measurements. Findings reveal: (1) Landsat data provides detailed spatial trends, identifying 18.7% of the area as significantly greening, whereas GIMMS data detects more browning due to spectral mixing; (2) GIMMS underestimates the greening to browning ratio at 2.2:1, compared to Landsat and MODIS ratios of 14.1:1 and 15.1:1, respectively; (3) Over 93% agreement exists between Landsat and MODIS or GIMMS trends, with discrepancies in limited areas. This highlights the importance of high-resolution data and field studies for accurately understanding vegetation trends across the arctic tundra.

KEYWORDS:

• Arctic greening
• Landsat
• MODIS
• AVHRR GIMMS_3g

1. Introduction

The Arctic has warmed approximately four times more quickly than the rest of the world since 1979 (Rantanen et al. 2022). Arctic warming influences how the tundra ecosystem functions due to its interactions with flora, wildlife, and populated areas. Arctic warming has the potential to boost ecosystem respiration and increase greenhouse gas emissions into the atmosphere due to degradation of permafrost. In the meantime, rising temperatures encourage plant growth and boost carbon uptake (Li et al. 2021). Researchers have shown that the Arctic was once lush and green, according to a study of plant DNA found in lake sediments, and if the area continues to warm due to climate change, it may someday revert to this state (Crump et al. 2021).

Satellite observations over tundra ecosystems have fully documented increases in vegetation productivity (Berner et al. 2020; Ju and Masek 2016; Liu, Huang, and Sun 2021). Some of present changes have been connected to shrub expansion and rapid yearly growth at locations throughout the Arctic (Berner et al. 2020; Berteaux et al. 2017; Forbes, Macias Fauria, and Zetterberg 2010; Loranty et al. 2018; Wheeler, Høye, and Svenning 2018). However, the Arctic shrub expansion is spatially unequal. Studies have discovered that willow shrubs are growing more frequently on outcrops of higher salinized sea rocks (Ukraintseva et al. 2014). A considerable contribution to the greening of the tundra can also be made by the emergence of highly productive herbaceous plant communities in the basins of drained thermokarst lakes (Chen, Liu, and Cheng 2021; Loiko et al. 2020; Veremeeva et al. 2021).

At present, we still don't fully comprehend this prevalent and significant ‘Arctic greening’ process. This is caused, in part, by the fact that going green is not a consistent, all-encompassing reaction, but rather is heterogeneous. Throughout the Arctic, the rate of ‘greening’ varies both spatially and temporally, while some areas show ‘browning’. Accurate predictions of future carbon, nitrogen, and water cycling in the Arctic are significantly hampered by these complexity and substantial uncertainty on future vegetation change in the region (Berner et al. 2020; Phoenix and Treharne 2022).

To assess the Arctic vegetation response to climate change or anthropogenic influences, it is common practice to assess greening by examining maps showing statistically significant greening or browning trends (every grid cell with a p-value below the nominal value is statistically significant) (Chen et al. 2019; Ju and Masek 2016; Liu, Huang, and Sun 2021). Satellite-derived vegetation indices have been used to assess these changes at high northern latitudes over the past four decades (Berner et al. 2020; Elsakov 2017; Frost et al. 2020; Ju and Masek 2016). The most used vegetation indicator for measuring the greenness of aboveground vegetation is the normalized difference vegetation index (NDVI), which has a strong relationship with vegetation health and productivity which can demonstrate various alterations in the features of the land surface properties. Although the NDVI is not sensitive to changes in land surface attribute at high values (Pastick et al. 2019), most of the maximum NDVI values over circumpolar tundra are in low value ranges that are less than 0.6 (Walker et al. 2002). Therefore, NDVI was selected for the calculation of greenness trend in this study. The three most often used datasets come from the USGS Landsat (1972- present), the Moderate Resolution Imaging Spectroradiometer (MODIS; 2001- present), and the Advanced Very High-Resolution Radiometer (AVHRR; 1982–2015). The NASA Goddard Space Flight Centre’s Global Inventory Modelling and Mapping Studies (GIMMS) NDVI was discovered to be the best for AVHRR vegetation monitoring. Index-3rd generation is the most recent iteration of the GIMMS dataset. Although AVHRR is well designed for circumpolar-scale trend detection, the coarse spatial resolution of AVHRR NDVI data (typically 8 km) fails to satisfy the requirements for straight linking sensor signal with ecological change in heterogeneous tundra landscapes, which further restricts the ability to conduct an attribution analysis on greenness trends, which includes permafrost thaw (Verdonen et al. 2020), wildfires (Jenkins et al. 2014; Zhou et al. 2019), and human disturbances (Berner et al. 2020; Liu, Huang, and Sun 2021). Furthermore, it is challenging to correlate trends with field measurements due to the coarse spatial resolution (Berner et al. 2020).

Satellites with higher resolution, like Landsat, can lessen but not always eliminate the spectrum mixing. The Landsat satellites have been utilized to assess regional or even local trends in tundra greenness with vegetation index like NDVI or EVI. The prior greenness trend studies over the arctic tundra with Landsat satellites mainly focused on reginal sites in Alaska, Canada, and the Russian tundra for two or three decades (Ju and Masek 2016; Liu, Huang, and Sun 2021; Miles and Esau 2016; Pastick et al. 2019). Even though some researchers (Berner et al. 2020) extracted summer Landsat data for 50,000 randomly selected sampling sites (each with a radius of 50 m) in the circumpolar tundra and discovered that greening (tundra greenness) rose in about 37.3% of sampling sites and browning (tundra browning) increased in about 4.7% of sampling sites from 1985 to 2016, they did not show those changed magnitude pixel by pixel and it is still unclear how they are distributed geographically. In order to improve our ability to predict how the tundra will change in the future, it is imperative that we understand how the tundra as a whole is changing over time and space, as well as how these changes may impact the processes that shape the landscape.

Although it is difficult to quantify how much the scale of measurement affects the sometimes-contradictory greening and browning signals reported in various spectrum datasets (Berner et al. 2020), the variations in Arctic greenness trends can be examined using various data resolutions. In this work, we reanalyse circumpolar arctic tundra greening trends pixel by pixel based on Landsat NDVI datasets between 2000 and 2020. By doing this, we demonstrate the ability of the Landsat archive to identify vegetation trends throughout the whole circumpolar arctic tundra and offer the highest resolution spatially complete trend maps of the region to date. We used the same trend analysis method with rigorous significance testing to compare Landsat-based trends with MODIS-based and AVHRR-based trends in order to better understand the regional variability of the reported greening/browning trend. Our goals are to (1) generate a spatially explicit map of vegetation dynamics across the circumpolar arctic tundra at 30 m resolution between 2000 and 2020, (2) compare the fine and coarse data-driven results to find out whether satellite data at coarse resolution misleading the vegetation greening due to the ignorance of spatial heterogeneous at a fine scale, and (3) what are the differences in arctic greening and browning across bioclimatic zone and country level?

2. Data materials and method

2.1. Study areas

The study region was defined in the terrestrial Arctic. The terrestrial Arctic is defined as the northernmost part of the Earth characterized by tundra vegetation, an arctic climate and arctic flora, with the tree line and continental coastlines jointly determining the extent of borders (Figure 1). Spatially, the present study covers an area of approximately 7.11 million km², overlapping with parts of six countries including Canada (CA), Denmark (Greenland, GR), Iceland (IC), Norway (NO), Russia (RU), and the United States (Alaska, AK). The tundra regions of these nations have all experienced different levels of anthropogenic disturbance and Russian Arctic has the most densely populated clusters of impervious surfaces (Liu, Huang, and Sun 2021; Liu et al. 2022a; Liu et al. 2022b; Xu et al. 2022).

Figure 1. The distribution of circumpolar arctic tundra which lives north of arctic treeline.

Based on the lowest to highest mean July temperatures, the Circumpolar Arctic bioclimate is classified into five subzones, A to E (Walker et al. 2002) (Figure 1). The coldest subzone is A, where mean July temperatures range from 0 to 3°C, while the hottest subzone is E, where mean July temperatures range from 9 to 12°C. The mean July temperature at the most southern boundary of subzone B, C and D are about 5°C, 7°C and 9°C, respectively. Sub-area A is mostly composed of bare rock, has short vegetation heights, 40% lichen and moss cover, and has a total biomass of less than 3 t/ha; sub-area B has a two-layered vertical structure and has a herbaceous layer up to 10 cm tall, 60% cryptogamous plant cover, and a total biomass of 5–20 t/ha; The vascular plant cover is 5–50% and the height of the creeping and semi-prostrate dwarf shrubs in sub-area C can reach 15 cm, forming an open patchy vegetation zone. In sub-area D, the herbaceous and dwarf shrub layer grows to a height of 40 cm, and vascular plants cover 50–80 percent of the subarea, creating an intermittent and closed canopy zone. The vegetation in the E zone is more densely populated and has two to three layers of vertical structure. The shrub layer can occasionally reach a height of 80 cm. Elevation has a significant impact on the patterns of tundra plant communities and soil moisture. Most of the study area are low-elevation plains less than 100 metres above sea level except some mountains and hills around 1000 metres above sea level located in Alaska and eastern Siberian arctic tundra. Erect shrublands make up around 26% of the vegetated area, followed by peaty graminoid tundras (18%), mountain complexes (13%), barrens (12%), mineral graminoid tundras (11%), prostrate-shrub tundras (11%), and marshes (7%). Russia has the greatest area in the Low Arctic, largely low-shrub tundra, whereas Canada has by far the biggest terrain in the High Arctic, mostly linked with plentiful barren kinds and prostrate dwarf-shrub tundra (Walker et al. 2002).

2.2. Landsatbased NDVI datasets

This study utilizes Landsat 5 TM, Landsat 7 ETM+, and Landsat 8 OLI orthorectified top-of-atmosphere (TOA) reflectance data (L1 T) Tier 1. The TOA data were utilized instead of surface reflectance data because they are compatible with previous research due to known errors associated with Landsat SR data over high latitudes (above 65°N) and the lack of MODIS data required to constrain atmospheric correction models for a section of the time series (1984–2000) (Arndt et al. 2019; Ju and Masek 2016; Pastick et al. 2019). All the employed Landsat pictures have a 30 m spatial resolution and were collected between 2000 and 2020 during the melt season (from 1 July to 30 August). Figure 2 displays the number of Landsat observations made in the Arctic between 1984 and 2020 which shows Landsat archives had numerous (more than 20%) data gaps during the 1980s and 1990s, but could cover the entire region after the 2000s. The cloud cover in the TOA reflectance images was assessed using the SimpleCloudScore technique implemented in the Google Earth Engine (GEE). According to our visual interpretation, the threshold score at 20 was set to remove pixels with high cloud. For each land pixel after masking cloud, shadow, and snow with band ‘'BQA'’ in the TOA dataset, the NDVI time series data were calculated. Compared to the previous Landsat-7 Enhanced Thematic Mapper (ETM +), the Landsat-8 Operational Land Imager (OLI) offers better calibration, signal-to-noise characteristics, a greater 12-bit radiometric resolution, and narrower spectral bands, so cross calibration between the two sensors are required before time series analysis. The data from Landsat 7, and 8 were cross-calibrated using the equations from already existing studies (Roy et al. 2016) in which statistical functions have been given using ordinary least squares (OLS).

Figure 2. Number of valid Landsat observations during two periods (a) 1984–1999 and (b) 2000-2020.

2.3. MODIS based NDVI datasets

The MOD13Q1 product of the Moderate Resolution Imaging Spectroradiometer (MODIS) provides global vegetation index observations at a spatial resolution of 250 m every 16 days. The version number of MOD13Q1 vegetation index product used in this study is 006. The product's V006 iteration delivers more spatial and temporal consistency and lessens systematic inaccuracy caused by sensor degradation. The MOD13Q1 package includes 12 bands, including red, green, blue, and mid-infrared reflectance, as well as quality control (QA) and observation azimuth bands, in addition to the two commonly used vegetation indices, NDVI and EVI. Only high-quality pixel observations were chosen for further investigation. Each image pixel had its unique QA band utilized to measure the observation quality level. In order to obtain the highest NDVI sequence during the growing season, the best inter-annual synthesis was performed based on this information using all valid observation sequences during melt season from 2000 to 2020.

2.4. AVHRR based NDVI datasets

AVHRR GIMMS NDVI3 g v3.1 time series (hereafter GIMMS) (Pinzon and Tucker 2014) was also used and it was download from NCAR (https://climatedataguide.ucar.edu/climate-data/ndvi-normalized-difference-vegetation-index-3rd-generation-nasagfsc-gimms) to compare geographic and trend differences. The GIMMS data are two estimates per month with a spatial resolution of 1/12 degrees (∼8 km at the equator). A maximum NDVI composite in every year (2000-2015) in July and August were derived and then used for trend calculation.

2.5. NDVI trend analysis

Every land pixel, with the exception of non-water and non-snow places, was used with the NDVI time series data during melt seasons. Valid time series were considered as having at least two observations within the research period. A Theil-Sen slope estimator was used to determine linear trends for the NDVI time series. Theil-Sen slope estimator performed better than a conventional least squares regression analysis since it is insensitive to outliers with a breakdown point of roughly 29% by using median of slope estimates [16]. The trend analysis of Landsat NDVI time series in Alaska [11], Landsat Tasseled Cap (TC) time series in the Canadian North-West (Fraser et al. 2014), and the Arctic Lena Delta in Russia [17–19] have been effectively conducted using Theil-Sen slope estimator because accurate slope estimates can be produced by Theil-Sen slope estimator in the face of arbitrarily extreme observations. The slope between each data pair was obtained at every single pixel, and the median value was then computed as follows. (1) $T_{\mathrm{slope}}=\mathrm{median}({\displaystyle \frac{x_j-x_i}{j-i}),\mathrm{\forall j}>i}$(1) Where $T_{\mathrm{slope}}$ is the Theil-Sen slope, $x_i$ and $x_j$ are the NDVI values in the time of $i$ and $j$.

The annual NDVI trend is set as the Theil-Sen slope. When there are more valid data, the trend prediction accuracy might be elevated. A greening trend has a considerable positive slope, while a browning trend has a significant but negative slope. The remaining pixels other than greening and browning are considered to be stable regions where the greenness remains constant.

The Mann-Kendall significance test can be used to validate the significance of Theil-Sen slope estimator. The primary benefit is that it does not require the data to follow a normal distribution and is a widely used for trend test. The following is the relevant calculation formula: (2) $Z_{\mathrm{MK}}=\left\{\begin{array}{c}{\displaystyle \frac{S-1}{\sqrt{n(n-1)(2n+5)/18}},S>0}\\ 0,S=0\\ {\displaystyle \frac{S+1}{\sqrt{n(n-1)(2n+5)/18}},S<0}\end{array}\right\}$(2) (3) $S={\sum }_{k=1}^{n-1}{\sum }_{j=k+1}^n\mathrm{sgn}(x_j-x_k)$(3) (4) $\mathrm{sgn}(x_j-x_k)=\left\{\begin{array}{c}1,(x_j-x_k)>0\\ 0,(x_j-x_k)=0\\ -1,(x_j-x_k)<0\end{array}\right\}$(4) Where $S$ is the Mann-Kendall test statistics, $x_j$ and $x_k$ are the NDVI values at the time of $j$ and $k$. $n$ is the time length of the NDVI values. The threshold of significance level is $\pm 1.96$ (95% significance test). Only the altered pixels were included in the final results because trends with p < 0.05 were deemed to be statistically significant in this study. Figure 3 shows the workflow of this study.

2.6. Trend comparison

The comparison of three spatial resolution datasets (Landsat, MODIS, and GIMMS NDVI3g) involved assessing spatial consistency, where we summarize the agreement percentages between Landsat and MODIS or Landsat and GIMMS. Before comparison, bilinear interpolation was applied to resample Landsat results to match the spatial resolution of the other datasets. Additionally, we calculated and compared the green-brown ratio among the datasets. The comparison can help to facilitate a detailed examination of vegetation trends in the circumpolar arctic tundra, highlighting spatial agreement, resolution matching, and overall trends among the datasets.

Figure 3. Workflow of this study.

3. Results

3.1. Landsat-based NDVI trends over the past two decades

Over the past two decades, arctic tundra has witnessed widespread greening (Figure 4). 18.7% of the region experience greening trend at 95% significance level. Only 1.3% of the region were found has vegetation browning. When we take a close look at the spatial pattern, the most notable greening was detected over the northern and southern parts of Russia's Yamal Peninsula and the central Gydan Peninsula, the northeastern Chukchi Peninsula, northern Quebec tundra, and the North Slope of Alaska. Although browning was much less common than greening, it was more likely to occur at higher latitudes in Canada and lowlands in central and eastern Russia. The average greening trend value was 0.003 NDVI units per year.

Figure 4. Greenness trend (at p < 0.05) over circumpolar tundra acquired with Landsat images from 2000 to 2020.

The results show that the area with changes in vegetation greenness is 1.04 million km², accounting for 20% of the total area of the study area. Greening, or growing greenness, makes up 95% of the altered regions, while browning, or declining greenness, makes up 5%.

Greening statistics can be summarized on the basis of national scale and bioclimatic subzone scale as shown in Figure 5. On a national scale, Canada has the biggest amount of greening areas, accounting for 29% of its total area, compared to 9% in the US and 11% in Russia. The difference in browning regions was not statistically significant (1%). The yearly rate of change of NDVI in greening regions was mostly centred in the 0.002–0.004 value range, with the US and RU having a small proportion (5%) of fast greening regions (annual rate of change of NDVI >0.008). At the bioclimatic subzones scale, areas of rapid greening (annual rate of change of NDVI more than 0.008) existed in subzones A and B, accounting for 17% and 9% of the total, respectively. Subzones C、D and E had nearly the same percentage of the rate of change in the 0.002–0.004 interval, which equated to around 55%. As latitude declined, the amplitude of greening rose.

Figure 5. Landsat-based NDVI trend statistics in country and bioclimatic subzone scale. (a). Spatial distribution of greenness trend. (b). statistics of changed and unchanged area percentage over three main countries. (c). greening magnitude over three main countries. (d). greening magnitude over 5 bioclimatic subzones.

3.2. Landsat greening trends compared to MODIS and GIMMS

The trend findings of MODIS and GIMMS are given in Figure 6, and the results show that the browning areas in the GIMMS results are quite big, cantered in the western and central sections of Russia, whereas the extent of the greening areas in MODIS is similar to that in Landsat results. The higher resolution trend results from MODIS and GIMMS can be found in Figure S4 and S5.

Figure 6. Greenness trend (at p < 0.05) over circumpolar tundra acquired with MODIS and GIMMS NDVI time series. (a). MODIS (b). GIMMS.

Table 1 and Figure S6 provide a thorough summary of the three data sources used to detect greenness trends over the circumpolar arctic tundra. Landsat NDVI demonstrates a stability rate of 79.98%, indicating that almost 80% of the circumpolar arctic tundra exhibit consistent greenness over the study period. MODIS NDVI follows closely with a stability rate of 79.76%, while GIMMS NDVI3 g shows a slightly higher stability rate of 83.05%. These stability percentage provide valuable insights that the vegetation productivity remains unchanged. Landsat NDVI indicates a greening rate of 18.7%, suggesting that approximately 18.7% of the area exhibit a positive trend of vegetation growth. In contrast, MODIS NDVI shows a slightly higher greening rate of 19%, while GIMMS NDVI3 g exhibits a lower greening rate of 11.6%. Similarly, Landsat NDVI indicates a browning rate of 1.3%, implying a decreasing trend in vegetation growth in 1.3% of the areas. MODIS NDVI and GIMMS NDVI3 g exbibit lower browning rates of 1.3% and 5.3%, respectively.

Table 1. Comparison in the percentage of vegetation greening trend over the circumpolar tundra across three spatial resolution data.

Data source	Spatial resolution	Stable (%)	Greening (%)	Browning (%)
Landsat NDVI	30 m	79.98	18.69	1.33
MODIS NDVI	250 m	79.76	18.98	1.26
GIMMS NDVI3g	8 km	83.05	11.62	5.34

Figure 7 summarizes the stable, greening, and browning areas that were obtained from the three NDVI datasets with regard to the greening tendency in each bioclimate zone. Overall, the areas of greenness change grow as latitude decreases. In Subzones A and B, Landsat and GIMMS both found small areas of greening in addition to patches of browning, while MODIS record few browning spots. The results of the three datasets are more dissimilar in Subzones C, D, and E. For Subzones C, D, and E, the findings from Landsat and MODIS were mostly similar, with greening areas ranged from 9 to 27 times greater than browning areas. However, the GIMMS data showed a less pronounced size difference between areas that are greening and those that are browning.

Figure 7. Area percent of greening trend in each bioclimate zone.

There is a significant variation between the three datasets regarding the area of change in greenness (Figure 8). The biggest green-brown ratio, 15.1:1, was discovered by MODIS data, while 14.1:1 and 2.2:1 were computed by Landsat and GIMMS, respectively. The name of the bioclimate zones with the highest percentage of greening is different when it comes to the bioclimate zones of greening. Subzone D (41.251%) has the highest greening value according to Landsat data; however, subzone E has the highest greening value according to MODIS and GIMMS results. The name of the bioclimate zones with the highest level of browning in the browning zones is a distinct one. With percentages of 37.558% and 34.50%, respectively, subzone C has the largest browning proportion in the data from Landsat and GIMMS, whereas subzone E has the highest browning proportion in the results from MODIS (with a percentage of 45.54%). Overall, NDVI data at various resolutions demonstrate disparities in the bioclimate locations where changes in vegetation greenness occur.

Figure 8. Variation between the three datasets regarding the area of change in greenness over bioclimatic subzones. a. Landsat b. MODIS c.GIMMS.

The trend consistency between the Landsat, MODIS, and GIMMS products is distributed spatially in Figure 9. The results from Landsat and MODIS and GIMMS agree well spatially, with a 95.32% and 93.95% agreement, respectively. The disparity, however, is greater in limited locations. Before the three datasets were compared, bilinear interpolation was used to resample the Landsat result to 250 m and 8 km, respectively. The consistency of Landsat's greening pixels with GIMMS and MODIS is 15.8% and 34.61%, respectively. The reason for the low agreement is that, particularly in the northern Canada and eastern Russia regions where the divergence is greater, the probability that GIMMS and MODIS deem them to be stable pixels exceeds 80% and 65%, respectively. In pixels where Landsat is stable, GIMMS and MODIS exhibit high agreement, with 97.67% and 98.06%, respectively. This may explain why the coarse resolution data overestimated vegetation greenness, with areas of greater disagreement with GIMMS occurring in the Alaska region, portions of the southern Yamal and Gydan peninsulas of Russia, and central Russia and areas of greater disagreement with MODIS occurring in the northern Alaskan, southern Canadian, central, and eastern Russian tundra. All these regions are all considered to be greening by both GIMMS and MODIS. Additionally, MODIS observes some browning in western Russia. Despite the modest browning area in the tundra region, the three datasets are noticeably more different in terms of Landsat browning pixels, with inconsistency rates of 82.9% (GIMMS) and 88.31% (MODIS), respectively.

Figure 9. Comparison of Landsat trend’s consistency with MODIS and GIMMS trend. The pixel frequency distribution (%) of the agreement and disagreement is displayed as a bar plot in each panel. When both two trends are constant, the agreement is displayed in blue; if not, it is displayed in red or green. A circle denotes no change, a + indicates greening, and a - indicates browning.

4. Discussion

4.1. Advantages of greenness trend derived from Landsat NDVI

In order to provide a pixel-by-pixel analysis of the greenness trend in the circumpolar Arctic tundra with fine spatial resolution, Landsat NDVI data collected between 2000 and 2020 were used. Apart from a few locations (less than 1%), where data gaps existed, more than 100 valid observations per pixel were used for trend analysis from 2000 to 2020 (Figure 2). Landsat data on greenness trends indicated that approximately 18.69% and 1.33% of the Arctic were greening and browning, respectively, at a significance level of 0.05 when water bodies were eliminated in the first stage. The Yamal Peninsula, the Gydan Peninsula, and the northeastern Chukchi Peninsula in Russian tundra, southeast Canada where shrubs and grasses are growing thicker (Ju and Masek 2016), and northwestern Alaska all experienced significant greening, which is consistent with earlier studies that used sampled Landsat time series (Berner et al. 2020). The values of the NDVI/year greening trend are consistent with those obtained in other arctic regions (Cheţan et al. 2020; Fraser et al. 2014; Ju and Masek 2016; Nitze and Grosse 2016; Verdonen et al. 2020).

Furthermore, a number of examples of locations that are clearly getting greener were provided through a review of the literature and analysis of very high-resolution remote sensing images (Arndt et al. 2019; Davis et al. 2021; Edwards and Treitz 2017; Myers-Smith et al. 2019; Verdonen et al. 2020) as shown in Figure 10. These locations have all undergone a totally greening process except a small portion of the southern areas in Cape Bounty Arctic Watershed Observatory (CBAWO) in Figure 10(b). Although CBAWO showed a clear greening trend overall, especially in the dry and mesic vegetation types (Edwards and Treitz 2017), our Landsat trend demonstrates browning in the polar semi-desert types and GIMMS also shows a magnitude of −0.002 browning over the entire region. In other locations, Landsat and MODIS trend results show relatively high-level spatial consistency. However, GIMMS has no findings in Mordy Yakha in Russia and Barrow in Alaska because of coarse resolution as shown in Figure 10(d) and 10(e). These case studies allow us to conclude that Landsat's trend maps are more trustworthy than coarser dataset with finer resolution spatial details.

Figure 10. Greening cases retrieved from Landsat, MODIS and GIMMS. These greening cases are validated by very high-resolution data or field photos. These five spots are (a) Torngat Mountains National Park (TMNP) in Canada, centred at 63.544°W,58.781°N. (b) CBAWO in Canada, centred at 109.531°W,74.885°N. (c) Herschel Island in Canada, centred at 139.012°W, 69.591°N. (d) Mordy Yakha in Russia, centred at 68.481°E,70.167°N. (e) Barrow in Alaska, US centred at 156.613°W,71.295°N. Colour code of trend values are the same as those in Figures 4 and 6, respectively.

Disparities between the datasets may be seen in the browning and greening areas’ spatial distribution, especially in Subzones A, B, C, D, and E. In these subzones, Landsat and MODIS display greater greening areas than GIMMS. The disparities in trend results for the Cape Bounty Arctic Watershed Observatory (CBAWO) underscore the significance of dataset resolution, as Landsat offers more precise spatial information.

4.2. Possible reasons for the differences between the three datasets

Different stability rates are shown by Landsat, MODIS, and GIMMS in Figure 9, suggesting that greenness has persisted throughout the research period. The smoothing effect and poor resolution of GIMMS could be the reason for its greater stability. Interestingly, MODIS's greening rate is higher than Landsat's, indicating that they may be more or less sensitive to changes (Berner et al. 2020). The agreement percentages show that there is generally good spatial agreement between Landsat, MODIS, and GIMMS. Nonetheless, the reduced agreement in some areas, especially in eastern Russia and northern Canada, raises the possibility of constraints, maybe related to the coarse resolution of GIMMS and MODIS. The necessity of resolution matching for precise comparison is emphasized by the interpolation of Landsat data.

Before findings at various resolutions can be compared, the spatial resolution needs to be harmonized. Resampling ensures spatial consistency, allowing for the integration of datasets with different native resolutions. By using bilinear interpolation, which can increase the spatial distribution of continuous values to the maximum extent possible, the original 30 m Landsat NDVI trend could be aggregated to 250 m and 8 km, respectively, to match the resolution of the GIMMS NDVI and MODIS NDVI trend results, which can provide broader perspective on vegetation trends across the entire circumpolar tundra and is useful for understanding landscape-level changes and regional patterns. However, the smoothed trends may mask localized variations in vegetation health or changes, potentially affecting the interpretation of ecological process. In order to investigate the reasons behind the contradicting data, the discrepancies between Landsat vs. MODIS and Landsat vs. GIMMS are summarized in Figure 11. Although Landsat's results provide a better representation of spatial detail, there may be information loss during the downscaling process when using bilinear interpolation in resampling based on the weighted average distance between the centres of the four closest input pixels.

Figure 11. The discrepancies between Landsat vs. MODIS and Landsat vs. GIMMS due to resampling. Letter A denotes a typical greening area, and letter B denotes a typical browning area in the results of the three datasets for these two locations, which are listed in (c) and (d). Agreement means both has the same greenness trend (stable/greening/browning) and disagreement means neither of the two maps has the same greenness trend.

A potential reason for the significant browning regions observed in the GIMMS data could be the disparity in the time interval for acquiring data between MODIS and Landsat. Since the Landsat and MODIS dataset covers the years 2000–2020but the GIMMS dataset we used covers the years 2000–2015. During the year around 2014, browning in GIMMS result may be caused by a variety of factors. Several different types of unusual weather events were discovered that induced browning in the record-low productivity for north-west Scandinavia reported in 2012 (Bjerke et al. 2014). Although they happened further south than our research area, higher latitudes may act as an early warning system. While this might indicate that the Arctic is at risk of browning, vegetation also shows hints of rapid recovery from such episodes. Several studies (Jenkins et al. 2014; Zhou et al. 2019) have shown that plants can recover from tundra fire within 2–5 years as well as indicated in our study shown in Figure 12. Those fire events are considered as browning factors and occur at different times. Allman Lake fire started in 2009 (Figure 12(a)) and Anaktuvuk River fire (Figure 12(b)) started in 2007 (Zhou et al. 2019) and both of the vegetation is fully restored in 2020 and the vegetation is greener than it was before the fires in our Landsat and MODIS results. Nevertheless, Gechia Lake fire scar is browning in our trend results because it was a recent fire occurred in 2015 and has not recovered yet. The annual NDVI maps retrieved by Landsat data are listed in Figure S1–S3 and can confirm our trend results.

Figure 12. NDVI trend after browning event caused by tundra fires in Alaska. (a)-(c) are the trend results from Landsat, MODIS and GIMMS around Allman Lake fire scar. The black polygon is the maximum boundary of Allman Lake fire scar. (d)-(f) are the trend results from Landsat, MODIS and GIMMS around Anaktuvuk River fire scar. The black polygon is the maximum boundary of Anaktuvuk River fire scar. (g)-(i) are the trend results from Landsat, MODIS and GIMMS around Gechia Lake fire scar. The black polygon is the maximum boundary of Gechia Lake fire scar.

Another explanation for the disparity in the results of the three resolutions could be that small-scale biological changes in vegetation are overlooked by coarse resolution analysis since several of the areas that GIMMS perceive an area to be static when Landsat or MODIS deems it to be changing (greening or browning) shown in Figure 9. Figure 12 provides additional examples of how tundra fire has changed localized vegetation to highlight the issue. GIMMS missed the greening trend in Allman Lake fire and the browning trend in Gechia Lake fire scar. Ecological contributions to spectral greening signals may coexist with browning signals from disturbances like retrogressive thaw slumps or vegetation trampling by herbivores, for example, within this aggregated pixel (Berner et al. 2020). Large-scale patterns of spatial variability in pixel-level browning and greening are influenced by spatial resolution. Higher resolution satellites like Landsat can lessen but not entirely stop this spectrum mixing. The significant discrepancies observed in Landsat browning pixels, particularly when compared to GIMMS and MODIS, imply difficulties in precisely identifying areas experiencing vegetation browning (Frost et al. 2021; Phoenix and Treharne 2022). This aspect holds crucial importance in comprehending the dynamics of Arctic tundra vegetation, and the elevated rates of inconsistency may signal constraints in the ability of coarser datasets to capture subtle changes effectively.

4.3. Limitations of the Landsat greenness trend over the circumpolar Arctic

The limited field measurements that were made could not confirm the general trend in greenness that was observed throughout the entire circumpolar tundra due to the research area's inaccessibility. On the other hand, a 14.1:1 ratio of greening to browning was obtained from 2000 to 2020 using Landsat NDVI time series. These ratios are partially consistent with the ratios found in existing study (Berner et al. 2020), who examined 50,000 sampling sites (50-metre circle buffer) across the whole Arctic. To further support the trends found in the Landsat NDVI, a number of hotspots that have been identified by other researchers using high-resolution data or field measurements were found. (see Figure 10 and 12). There is a high demand for higher spatial resolution observations in upcoming research (Berner et al. 2020), as well as for gathering field measurement data over focal scales and even the tundra biome scale (Cunliffe et al. 2020; Myers-Smith et al. 2019; Shabanov et al. 2021). These observations will aid in improving our understanding of vegetation greening in circumpolar areas in various vegetation types and geographic gradients.

The piecewise trend analysis of vegetation greenness is highly needed for Arctic tundra. Time series analysis techniques, such as Breaks For Additive Seasonal and Trend (BFAST) (Liu et al. 2019; Verbesselt et al. 2010a; Verbesselt et al. 2010b), CCDC (Zhu and Woodcock 2014) or LandTrendr (Kennedy, Yang, and Cohen 2010), which have been integrated into the platform of Google Earth Engine (Arévalo et al. 2020; Hamunyela et al. 2020; Kennedy et al. 2018), can be used to identify the breakpoints in the time series data. In addition, widespread active layer detachments and retrogressive thaw slumps are very common in the low arctic tundra zone (Verdonen et al. 2020). After a landslide, vegetation cannot fully regrow and the regeneration to original levels of composition, cover, and productivity could take several more decades. Those piecewise time series analysis algorithms are required for this regeneration study.

5. Conclusion

In conclusion, our comprehensive evaluation of pan-Arctic tundra greenness over the past two decades, using high-resolution Landsat satellite observations alongside coarser datasets, provides key insights into the dynamics of plant productivity and biomass. While a general increase in greenness was observed across 18.7% of the Arctic territory, our findings highlight nuanced patterns and challenges:

Disparities in the spatial distribution of browning and greening areas, notably in Subzones A, B, C, D, and E, underscore the pivotal role of dataset resolution. Landsat's finer resolution emerges as crucial for capturing detailed trends, as evidenced in comparisons with MODIS and GIMMS. The finer granularity of Landsat data allows for a more precise understanding of vegetation dynamics, particularly in heterogeneous Arctic landscapes.

Divergent stability rates among Landsat, MODIS, and GIMMS reveal variations in their ability to capture greenness persistence. GIMMS has a higher stability rate attributed to its smoothing effect, while MODIS displays a higher greening rate, suggesting differences in sensitivity to changes. These nuances in stability rates emphasize the importance of considering dataset characteristics when interpreting long-term trends in Arctic vegetation.

Notable inconsistencies in Landsat browning pixels, especially when compared to GIMMS and MODIS, highlight challenges in precisely identifying areas experiencing vegetation browning. These discrepancies underscore potential limitations in coarser datasets, emphasizing the necessity of finer resolution datasets to accurately capture subtle changes in Arctic tundra vegetation.

Challenges in precisely identifying vegetation browning areas, particularly in GIMMS data, include disparities in acquisition time intervals between datasets, spectrum mixing at coarser resolutions affecting detection accuracy, and the complex interplay of ecological contributions to spectral signals. These challenges emphasize the need for a nuanced approach to browning detection and the importance of incorporating additional contextual information for accurate interpretation.

In summary, our study quantifies the changing greenness patterns in the Arctic tundra and highlights the impact of dataset resolution on the observed trends, providing valuable insights for future research in understanding the intricate dynamics of this vital ecosystem.

Acknowledgements

We acknowledge Compton J. Tucker and Jorge E. Pinzon of NASA GSFC for providing GIMMS NDVI3 g data, and NASA for providing MODIS data products. We are appreciative of the USGS and Google Earth Engine's (https://earthengine.google.com/, last access: Feb 2024) assistance with Landsat data. Grateful appreciation to the reviewers for invaluable feedback and contributions to the enhancement of this manuscript.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

The data that support the findings of this study can be obtained from the corresponding author upon reasonable request.

Additional information

Funding

Research was supported by the Open Research Program of the International Research Center of Big Data for Sustainable Development Goals [grant number CBAS2022ORP04].