Arctic Greening Trends: Change Points in Satellite-Derived Normalized Difference Vegetation Indexes and Their Correlation with Climate Variables over the Last Two Decades

Abstract

In this study, we utilized NDVI data from the moderate resolution imaging spectroradiometer (MODIS) alongside climatic variables obtained from a reanalyzed dataset to analyze Arctic greening during the summer months (June–September) of the last two decades. This investigation entailed a detailed analysis of these changes across various temporal scales. The data indicated a continuous trend of Arctic greening, evidenced by a 1.8% per decade increment in the NDVI. Notably, significant change points were identified in June 2012 and September 2013. A comparative assessment of NDVI pre- and post-these inflection points revealed an elongation of the Arctic greening trend. Furthermore, an anomalous increase in NDVI of 2% per decade was observed, suggesting an acceleration in greening. A comprehensive analysis was conducted to decipher the correlation between NDVI, temperature, and energy budget parameters to elucidate the underlying causes of these change points. Although the correlation between these variables was relatively low throughout the summer months, a distinct pattern emerged when these periods were dissected and examined in the context of the identified change points. Preceding the change point, a strong correlation (approximately 0.6) was observed between all variables; however, this correlation significantly diminished after the change point, dropping to less than half. This shift implies an introduction of additional external factors influencing the Arctic greening trend after the change point. Our findings provide foundational data for estimating the tipping point in Arctic terrestrial ecosystems. This is achieved by integrating the observed NDVI change points with their relationship with climatic variables, which are essential in comprehensively understanding the dynamics of Arctic climate change, particularly with alterations in tundra vegetation.

1. Introduction

The Arctic region has experienced rapid climate change in recent decades, most notably a pronounced temperature increase compared to mid-latitude regions, a phenomenon attributed to Arctic amplification [1].

Consequently, various climate change events have occurred in the Arctic region. Vegetation, a key parameter in climate feedback, significantly influences temperature, surface energy budget, and hydrological balance [2,3,4]. Tundra vegetation dynamics profoundly influence permafrost, mitigating thawing in the summer and increasing soil temperature in snow-related winter [5]. Additionally, variations in atmospheric CO₂ concentrations drive vegetation changes, which impact climate change, primarily through alterations in albedo, evapotranspiration, and carbon stocks. These changes trigger various feedback mechanisms, for instance, the ice–albedo feedback where increased vegetation and increasing absorption of solar radiation accelerate Arctic warming [6]. Regional climate change is represented by different trends, such as cooling and warming, that significantly change the Arctic vegetation [7,8,9,10,11]. Therefore, a comprehensive understanding of Arctic vegetation is necessary to determine and respond to the large-scale significance of Arctic climate change [12].

Global warming has led to notable shifts in the density and distribution of tundra land cover [13]. The Arctic tundra is undergoing greening [8,14,15,16,17,18,19,20]. These vegetative changes are complicatedly linked to global warming, exerting substantial influence on the Arctic ecosystem. Understanding the trends and patterns associated with greening is crucial for deciphering its ecological implications and the direction of the Arctic climate. Previous satellite-based studies of Arctic vegetation have focused on overall greening trends [8,14,15,18,21,22,23]; however, the analyses of detailed growth characteristics, such as seasonal variations, annual variations, and change points, have been insufficient. Typically, these studies have leveraged various normalized difference vegetation index (NDVI) datasets such as the global inventory modeling and mapping studies (GIMMS) NDVI and the moderate resolution imaging spectroradiometer (MODIS) NDVI, with the latter being particularly advantageous for long-term, multi-scale trend analyses [21]. Different NDVIs, including annual peak NDVI (Max-NDVI) and time-integrated NDVI (TI-NDVI), have been employed to assess vegetation dynamics [8,14,21]. These NDVI datasets offer certain advantages, like enhanced outlier filtration, but also have limitations.

For instance, the variable dating of Max-NDVI complicates the analysis of time-series features at fixed intervals. At the same time, the smoothed values in TI-NDVI pose challenges in discerning temporal anomalies. Although it is possible to analyze smoothed NDVI data after removing outliers, it is difficult to analyze change points, the strength of the NDVIs, or other factors for specific periods without considering the aforementioned feedback of climate change. In summary, while existing literature elucidates the progression of Arctic greening using satellite data, it should better explain the temporal dynamics of specific phenomena, such as large-scale tipping points. Analyzing the change points in long-term time series data can help explain the dynamics of changes in Arctic terrestrial ecosystems, including the acceleration in variability, the cumulative change, and the climate variables that have more impact relative to the change point.

It is imperative to discern the climatic drivers influencing Arctic greening trends. Generally, a strong correlation exists between changes in Arctic vegetation and temperature [23]. Various research has investigated the relationship between tundra vegetation changes and climatic variables, including the summer warmth index (SWI), snow cover, sea ice extent, and snow thickness [2,14,21,24,25]. However, these studies often vary in approach and have regional, seasonal, and methodological differences in vegetation data. Notably, Arctic vegetation has shown a weak correlation with temperature over time [26]. The tipping point of the tundra was explained using the ice–albedo feedback [6,23,27,28]. Beyond temperature, parameters related to the energy budget are intricately linked to Arctic greening. Most studies of these parameters have been grounded on in situ or modeling approaches [13,24,29,30]. While substantial research has focused on temperature and sea ice, there is a growing need to examine the impacts of energy budget variables and indices on a larger scale. This necessitates a comprehensive analysis integrating temperature and energy budget changes with Arctic vegetation dynamics.

This study aimed to investigate the temporal dynamics and driving forces of Arctic greening trends over the past two decades by employing MODIS NDVI data. We utilized monthly NDVI data, which are less affected by the transformation of peak periods and smoothing, to detect the change point in Arctic vegetation trends using Bayesian statistics. Based on this change point, we analyzed how the characteristics of Arctic greening change across various temporal scales. In addition, we analyzed the correlation between Arctic greening and both temperature and energy budget parameters (e.g., net radiation and turbulent heat flux) in the periods before and after the identified change point to determine the climatic factors governing these patterns across different temporal scales.

2. Materials

2.1. NDVI

NDVI is a vegetation index that uses the spectral properties of vegetation to evaluate their approximate health. NDVI is calculated using the reflectance difference between the near-infrared and visible spectra in the following way [31]:

NDVI = (NIR − Red)/(Red + NIR)

where NIR is the reflectance in the near-infrared region and Red is the reflectance in the red region. As most plants absorb solar radiation in the photosynthetically active radiation (PAR) spectral region, which ranges from 400 to 700 nm (visible), and reflect/re-emit energy at longer wavelengths (near-infrared), higher NDVI values represent healthier vegetation. NDVI has been widely used in many studies of vegetation dynamics [32,33]. Variations in NDVI data are driven by weather and ecosystem components; therefore, a vegetation condition index was developed to estimate the impact of weather on vegetation [34].

The MODIS NDVI dataset (MOD13C2V006) serviced by the National Aeronautics and Space Administration (NASA) was used for this study; the data covered the period 2000–2020. The data are formatted in the climate modeling grid (GMC) of 0.05-degree pixels that provide monthly averaged observations. According to the MODIS user guide [35], the valid range of the NDVI data is from −0.2 to 1.0; therefore, the pixels outside this range were first masked in our study.

2.2. Climate Parameters

Various variables have been collected to analyze the impact of climatic variables’ impact on changes to tundra vegetation [2,14,21,24,25]. This study focuses specifically on tundra vegetation and the complicated relationships between these variables that directly influence vegetation feedback.

Typically, vegetations thrive at high temperatures and under humid conditions; however, the Arctic region is generally characterized by its cold and dry climate. Due to climate change, the Arctic has become warmer and moister, and previous studies focused predominantly on temperature and atmospheric conditions. The previous studies used the summer warmth index (SWI) to represent temperature conditions [2,14,21,25], while this study used surface air temperature data from the same period as the NDVI data. Additionally, land surface temperature (LST) was used as a factor to analyze the impact of temperature on the ground surface where vegetation was present. Although albedo is often used as an energy budget related to surface conditions, this study examined net radiation, surface sensible heat flux (SSHF), and latent heat flux (SLHF) instead to provide a comprehensive interpretation of the surface radiation and heat budgets. The sign conventions of net radiation directions are as follows: positive for downward (atmosphere to surface) and negative for upward (surface to atmosphere). However, SSHF and SLHF use the opposite directions—positive for upward (surface to atmosphere) and negative for downward (atmosphere to surface).

Five climatic variables (

Table 1. Summary of climatic variables used in the study (all data were collected at monthly intervals).

Variables	Spatial Resolution	Data Source	Reference
Surface Air Temperature (SAT)	0.25 degree	ERA5 reanalysis product	[36]
Land Surface Temperature (LST)	0.05 degree	NASA MODIS product	[37]
Net Radiation (Rn)	1 degree	NASA CERES EBAF product	[38]
Surface Latent Heat Flux (SLHF)	0.25 degree	ERA5 reanalysis product	[36]
Surface Sensible Heat Flux (SSHF)	0.25 degree	ERA5 reanalysis product	[36]

) were collected and all the data were extracted from satellite observation (MODIS) and climate reanalysis (ERA5) datasets. Snow and ice data were not included as the study was conducted for the summer season.

3. Methods

3.1. Study Area and Period

The focus of this research was the terrestrial Arctic region [39,40], a diverse landscape comprising barren lands, graminoid, prostrate shrub, and erect shrub vegetations, wetlands, and ice, as categorized by the circumpolar Arctic vegetation map (CAVM) [41]. The study area is poleward of the tree line, with low summer temperatures (an average of under 12 °C in July [41]), limited annual precipitation, and a short growing season [5]. The vegetation comprises mosses, lichens, dwarf shrubs, sedges, grasses, and rushes [5,42]. Notably, the vascular plant cover, including graminoid cover, shrubs, and trees, shows an increasing trend [5].

Due to the reliance of NDVI measurements on solar radiation, data acquisition is constrained during polar nights. NDVI data is obtained from the red and NIR bands, the primary source of which is solar radiation. Hence, a preliminary assessment of the spatiotemporal validity of the dataset is critical. This involves analyzing the observation frequency at each pixel location, guided by the quality control annotations for monthly pixel reliability provided in the MOD13C2V006 dataset. For a comprehensive understanding of the quality control methodologies employed, consultation of the MODIS user guide [35] is recommended.

Spatial distribution maps of valid NDVI points (Figure 1) reveal a significant presence of valid observations primarily between June and September. Outside this temporal window, the proportion of valid data points dropped below 20% during the study period. A notable reduction in valid data points was observed at higher latitudes, particularly in the vicinity of the Laptev and Kara Seas in the Russian Arctic regions, where the percentage of valid observations was markedly low. Additionally, a temporal analysis based on the monthly time series of valid data points indicated an increase in missing values, particularly in June of the early 2000s.

Therefore, the research period spanned from June to September, based on the spatial and temporal distribution of valid NDVI observations. Data were utilized in two approaches: Monthly (June–September) and Summer (average of data from June to August).

3.2. Characteristics of NDVI Time Series Data

This study aimed to characterize the trends in Arctic NDVI data from the last 20 years and assess the quantitative relationship between NDVI and climate variables. To achieve these objectives, we conducted various statistical evaluations and correlation analyses as follows.

We explored the spatiotemporal dynamics within the study period using multiple statistical methods. Arctic greening and browning can be described from a remote sensing perspective [26]. Greening refers to an increase in vegetation indices, while browning represents the reverse [18]. Selecting a robust statistical approach is crucial for estimating Arctic greening, given the pronounced seasonality in vegetation changes [15,18,21].

For trend detection, we applied the Mann–Kendall test [43,44] and Sen’s slope estimator [45]. The Mann–Kendall test is a widely recognized method for trend analysis in meteorological and hydrological studies [46,47], with foundational details provided by Mann [43] and Kendall [44]. In the research, the significance level was set at a p-value of 0.05 in the Mann–Kendall test. Although the Mann–Kendall test is commonly used to identify annual trends in various research domains, caution is advised when interpreting these trends, especially when they are nonlinear and exhibit change points. Sen’s method offers a nonparametric approach for estimating trends in time series data [45].

The change points in NDVI data can be used as reference points for land-based changes in the analysis of climate change. To detect change points in NDVI time series data, we employed the Bayesian model averaging time series decomposition algorithm called Bayesian Estimator of Abrupt change, Seasonality, and Trend (BEAST) [48]. BEAST is a widely used method in various fields, including phenology, vegetation changes, and dam displacement [49]. BEAST decomposes the time series into seasonal and trend changes to quantify change probabilities using a Bayesian model. The following equation includes the four decomposed components: trend, seasonal variation, abrupt change (i.e., change points), and noise (i.e., residual).

$$y\left(t\right)=S\left({\theta }_s\right)+T\left({\theta }_T\right)+\epsilon$$

where y(t) is the time series (NDVI in this study), S is the seasonality, T is the trend, ${\theta }_s$ is the change point of the seasonal signal, ${\theta }_T$ is the change point of the trend signal, and ε is the gaussian random error [48]. We used annual NDVI data, which are time series data with seasonality removed, as only one value is used for each year. Therefore, the seasonality-related component can be ignored, simplifying the time series decomposed in this study as follows [48]:

$$y\left(t\right)=T\left({\theta }_T\right)+\epsilon$$

In our study, we performed a time series analysis of NDVI data, dissecting them into trends, change points, and residuals. This allowed us to closely examine NDVI trends occurring before and after each identified change point. In the results of change points detected using BEAST, two crucial factors, namely number of change points and ranked by probability of occurrence, are essential for obtaining a more accurate interpretation of the change points detected.

The term number of change points refers to the change points observed in the analyzed time series data and their respective probabilities of occurrence. In this study, the change point occurrence probability was calculated for up to four change points, and subsequently, the case with the highest probability was determined. The mode denotes the number of change points with the highest probable occurrence.

The term ranked by probability of occurrence refers to the list of probable trend change points ranked by probability of occurrence. Therefore, it detects the points where change points can occur and provides information about their probabilities. The rank means priority and the trend change point means the time when a change point occurs. The rank is determined by the probability of change point occurrence. Notably, a change point does not necessarily occur when the detection process is performed. To deem the candidate a true change point, it is essential to verify primarily high probability and to consider the proportion of probabilities among the candidates. The number of ranks is not fixed because rank represents all periods within the significance level of the probability of occurrence.

Anomaly analysis was employed to differentiate between ongoing and accelerated NDVI changes. Generally, anomalies are calculated using differences in climatology. However, we conducted an anomaly analysis to determine whether the greening of the Arctic was accelerated based on the start year, as follows:

Anomaly = X_t − X₂₀₀₀

where X_t is the value of the time series data in the current year and X₂₀₀₀ is the value of the time series data in the starting year, 2000. These anomalies were computed relative to the baseline year in the time series. This method allowed us to assess the degree of change from the start year, with ‘anomaly’ in this context being distinctly defined by the ‘residual’ in the BEAST method.

3.3. Relationship with Climate Variables

In this study, we conducted a correlation analysis to examine the relationship between NDVI and various climatic factors over the study period. Our approach involved using a consistent NDVI dataset while aiming to establish a more comprehensive correlation with temperature and energy budget parameters. Furthermore, we distinguished correlations at both annual and monthly scales, specifically analyzing summer versus month-to-month variations. To match the spatial resolutions of different climatic variables, we employed the great circle distance (GCD) method.

We initiated our analysis using Pearson’s correlation to investigate the relationship between NDVI and five key climatic variables, as detailed in Table 1. This allowed us to compute correlation coefficients for each variable and to discern variations in the intensity of relationships before and after identified change points. Additionally, we analyzed the time series patterns of climatic variables through anomaly, differences between atmospheric and surface temperature, and the Bowen ratio (ratio of sensible heat flux to latent heat flux) to understand heat distribution dynamics.

It is also crucial to consider potential time lags when analyzing vegetation–climate relationships, as reactions may not be immediate. For NDVIs, valid pixels were primarily observed from June to September. To accommodate this, we explored time-lagged correlations by adjusting the periods of climatic variables in one-month steps relative to the fixed NDVI period. For example, we calculated correlation coefficients for the NDVI data from June to September against climatic variables from May to August (Month − 1), from June to September (Month + 0), and from July to October (Month + 1), effectively analyzing NDVI correlation with climatic variables over a span of three months both prior to and following the fixed period.

4. Results

4.1. Characterization of NDVI Time Series Data

4.1.1. Long-Term Spatiotemporal Trends in NDVI Data

The summer NDVI changes in tundra over the last two decades were estimated to understand large-scale vegetation changes. The annual trend was determined based on the summer NDVI using Sen’s slope, and the Mann–Kendall test was employed to extract the greening area. Figure 2 shows the resulting spatial distributions of the annual trends and the Arctic greening area.

The annual NDVI trends turned out to be positive throughout the tundra, regardless of the region. The spatial distribution of trends in Figure 2a is in line with those reported in previous studies. Annual changes in Max-NDVI and TI-NDVI also showed an increasing trend over the tundra based on MODIS or Landsat products [17,18,21]. Summer NDVI, Max-NDVI, and TI-NDVI commonly indicated significant Arctic greening. Moreover, the greening persisted for two decades from 2000 to 2020. Almost all the tundra areas undergo Arctic greening annually, as shown in Figure 2b. However, some areas are presented in shaded gray (indicating no trend) because they had positive trends but did not meet the significance level.

We conducted another analysis to assess large-scale changes and comprehensively evaluate how the entire Arctic tundra had changed over the given period. The summer NDVI data were spatially averaged across the entire Arctic tundra and the year-to-year changes were estimated, as depicted in Figure 2c. The time series data were decomposed into trends and residuals that describe the temporal evolution of NDVI. The black dotted line, red line, and gray shade represent the summer NDVI data, the trend calculated using BEAST, and the standard deviation of the trend, respectively. The trend indicates an increase of 1.8% per decade at a significant level (p < 0.001). The bar plot in Figure 2c displays the residuals, which are the differences between the NDVI and its trend. The time series residuals indicate that there were dynamic changes in 2013 and 2015; higher NDVIs exceeding 0.02 were observed compared to other periods.

4.1.2. Change Points in Monthly and Summer NDVI Data

The characteristics of the Arctic vegetation time series data were further analyzed using BEAST. Since residuals at a specific point in Figure 2c show a significant increase compared to other periods, it suggests the presence of a singularity during that specific period rather than a typical linear trend. Change point detection in both summer and monthly data was carried out using BEAST; the results are summarized in

Table 2. Summary of change points detected using BEAST.

	Number of Change Points		Ranked by Probability of Occurrence
	Mode	Probability (%)	Rank	Trend Change Points	Probability (%)
Summer	0	52.1	1	2012	32.8
			2	2016	10.2
June	1	48.2	1	2012	39.7
			2	2004	6.7
July	0	58.5	1	2012	18.7
			2	2004	7.7
August	0	47.7	1	2010	16.0
			2	2007	15.3
			3	2016	15.2
September	1	67.8	1	2013	73.6
			2	2004	4.4

. They are divided into ‘number of changepoints’ and ‘ranked by probability of occurrence’ as described in the ‘Method’ section.

First, we had to check the ‘Mode’ of the number of change points. Mode indicates the number of change points in the time series data. For summer (averaged), the mode was zero, indicating that there is likely no change point. In contrast, for June, the mode was one, implying that there is a possibility of one change point during the period. Based on this outcome, it was possible to infer the number of change points in the NDVI data during the particular period.

Second, we focused on ‘Rank’ and ‘Trend Changepoints’ in ranked by probability of occurrence. Rank is a ranking that occurs at a change point in a time series, and Trend Change is the corresponding period. Therefore, Mode and Rank must be considered when detecting change points. In other words, if the mode is 1, the change point is up to the first rank, and if the mode is 2, the change points are up to the second rank.

For the summer, the probability of a change point trend occurring was highest in 2012 (i.e., rank 1); however, the number of change point modes detected was zero, because the difference with the second highest (i.e., rank 2) was not sufficiently large. For that signal, no change point was observed, implying that the imaginary change point in 2012 was not acceptable, considering the trend. In contrast, for June 2012, a trend in change point occurred (i.e., mode 1), with a larger difference in probabilities between rank 1 and rank 2. The 48% probability indicates that there was one change point in 2012.

Overall, the summer NDVI data had a 52.1% probability of having no change points, suggesting that no distinct change point was identified in the summer NDVI data. However, during time series decomposition, the residuals had a dynamic pattern in 2013 and 2016, as shown in Figure 2c. To have a closer look at this, we searched for change points in each of the monthly NDVI data, and two change points were identified in the June NDVI and September NDVI data.

For the June NDVI data, one change point occurred with a 48.2% probability. The highest probability of change point occurrence was in 2012 at 39.7%. Similarly, for the September NDVI data, there was a 67.8% probability of one change point occurring, with the highest probability of 73.6% occurring in 2013. The years 2012 and 2013 were notable periods for changes in sea ice extent. The year 2012 marked a record low in sea ice extent; however, the area recovered quickly and was near average in 2013 [50]. The change points in NDVI occurred in 2012 and 2013. In other words, change points occurred in both land and ocean climate variables simultaneously.

4.1.3. Analysis of NDVI Time Series Data Based on a Temporal Scale

The time series analysis applied to the monthly NDVI data resulted in a more detailed characterization of the dynamic changes observed in the summer NDVI data. Figure 3 displays the monthly NDVI trend, the difference between before and after change point occurrence, and the anomaly in monthly NDVI. In Figure 3a, the trends for July and August NDVI data showed a steady increase with no change points, while the existence of change points in Table 2 was confirmed, showing clearly different trends for June and September NDVIs. The overall upward changes are common to all the monthly NDVIs; however, the jumps between 2012 and 2013 are specific for the June and September NDVIs only.

We estimated the changes in NDVI in seasonal growing cycles over the study period by analyzing how they changed before and after the detected change points. Since the change point in 2012–2013 occurred with high probability, we selected 2012 as the division year. The year-to-year averaged monthly NDVIs for 2000–2012 were termed NDVI_Pre, and those for 2013–2020 were termed NDVI_Aft. NDVI_Aft was higher than NDVI_Pre for all the summer months, as shown in Figure 3b.

Besides the overall difference between NDVI_Pre and NDVI_Aft, larger increases were also observed in June and September. This indicates a strengthened pattern over time in a growing seasonal cycle with the same trend.

The anomaly in the acceleration of Arctic greening was calculated based on the 2000 start year. If the anomaly represents an increase over time and is a positive value, it indicates an acceleration of Arctic greening from the starting year. The anomaly showed a tendency to increase across all the given periods. The anomaly indicates an increasing trend of 0.2% per decade, with a p-value of < 0.001 (p-value < 0.002 in the September NDVI) based on the Mann–Kendall test. In addition, all periods had positive anomalies, except the September NDVI. Although negative anomalies were observed in the September NDVI until 2012, they subsequently turned positive. This is in line with the change point detection results where the September NDVI increased significantly after 2013.

4.2. Relationship between Arctic NDVI and Climate Variables

4.2.1. Relationships between All-Study Period NDVI Data

We analyzed the relationships between vegetation and climatic variables at the same spatial and temporal scales over the past two decades. This analysis assists with explaining changes in Arctic vegetation.

Table 3. Correlation coefficients between vegetation and climate variables in the study periods.

	Surface Air Temperature	Land Surface Temperature	Net Radiation	Surface Sensible Heat Flux	Surface Latent Heat Flux
Summer	0.29	0.38	0.27	−0.13	−0.42
June	0.25	0.42	0.27	−0.28	−0.41
July	0.33	0.34	−0.23	0.24	0.26
August	0.21	0.26	0.43	−0.02	−0.37
September	−0.17	−0.24	−0.42	0.05	−0.06

represents Pearson’s correlation coefficients between NDVI and the seven climatic variables listed in Table 1. The highest correlation coefficient for each category is indicated in bold. The correlations estimated in this study were slightly weaker than those reported in the previous studies [2,14,21,25] by 0.45–0.7. The decreased correlation over time suggests that, in addition to temperature, other climate variables have become more influential. This was supported by comparing the relative strengths of the correlation coefficients between the seven climatic variables. For example, the surface latent heat flux (SLHF) exhibited a stronger correlation with summer NDVI than temperature, even though it was previously widely recognized that temperature drives vegetation change.

4.2.2. Relationship between Summer NDVIs

We calculated correlation coefficients for the before- and after-periods based on the change point to better analyze the relationship between vegetation changes and climate variables. We compared them with all periods (Table 3 and

Table 4. Correlation coefficients between vegetation and climate variables based on the 2012 change point.

	Surface Air Temperature	Land Surface Temperature	Net Radiation	Surface Sensible Heat Flux	Surface Latent Heat Flux
Before	0.65	0.69	0.55	−0.33	−0.59
After	−0.23	−0.05	−0.03	−0.22	0.27

). Compared to the correlation coefficients for the summer in Table 3, the correlation is halved when the period after the change point is included. Compared to the previous study [26], the correlation between climate variables and tundra NDVI decreases when more recent periods are included.

The time series patterns of climate variables were analyzed to understand the reasons for the low correlation with NDVI after the change point. For most variables, the interannual differences were amplified around the NDVI change point, resulting in a pattern of high variability (Figure 4a,b). Regarding temperature, both SAT and LST showed a low–high peak change on a 4-year cycle, and the variation increased, with higher amplification after 2012 than before.

The lowest peak in the difference between air and surface temperature occurred in 2015, a period of anomalously high summer NDVI. Around this time, the difference between air and surface temperature decreased and then began to increase again, and the trend of the difference changed from negative to positive (Figure 4c). This means that the mixing between the atmosphere and the surface became more robust, and the blocking of the two boundaries became stronger again around this period.

The energy budget variables showed a pattern similar to that of temperature; anomaly variation increased after 2012. The net radiation showed an increasing trend, and SSHF and SLHF showed decreasing trends. These indicate an increase in the energy absorbed by the surface, suggesting an increase in surface thermal energy. In addition, the Bowen ratio of the time series data was analyzed to determine the heat distribution in the atmosphere. The Bowen ratio pattern peaked in 2011 (Figure 4d). At this time, the Bowen ratio pattern switched from increasing to decreasing. The increase in sensible heat was followed by a decrease in latent heat, indicating that the Arctic has been experiencing an increase in atmospheric evaporation in recent years. The variability in all the variables has increased since 2012; this characteristic reduced the correlation between NDVI and climate variables.

4.2.3. Relationships between Monthly NDVI and Climate Variables

Monthly correlations were analyzed to detect the relationship between climate variables and NDVI during the year. Cross-correlation was calculated to determine the before and after values of the variables (Figure 5). Since the NDVI data are only available from June to September, we analyzed a three-month cross-correlation. As a result, all climate variables (temperature and energy budget) and the vegetation index have a time gap of one month. The results showed that all climate variables exhibited their highest correlations in the previous month, suggesting that the changes in temperature [26] and energy budget variables had a significant impact on NDVI. This finding implies that the change in NDVI, and hence vegetation, could be used as an indicator to understand Arctic climate change better.

5. Discussion

Our study analyzed the spatiotemporal dynamics of tundra vegetation over the last two decades, examining the evolution of satellite-derived NDVI and its correlation with various climatic variables. Consistent with previous research [8,14,18,22], we observed a long-term increasing trend of 1.8% per decade in summer monthly NDVI data, providing strong evidence for Arctic greening. Beyond this, our study utilized a state-of-the-art Bayesian method, BEAST, to provide the characterizations of vegetation changes, including detecting critical change points across the Arctic tundra. The time series residuals indicate that there were dynamic changes in 2013 and 2015; higher NDVIs exceeding 0.02 were observed compared to other periods. We suspect that these periods included the vegetation change points. This large-scale time series singularity could be used as a critical point for detecting changes in the Arctic terrestrial area and analyzing Arctic climate change.

Two significant change points were identified in June and September of 2012 and 2013, similar to the characteristics of sea ice extent [50]. This finding underscores the interconnection between terrestrial and sea ice changes, resonating with the existing literature on the association between reduced sea ice and vegetation dynamics [14]. The change points suggest several possibilities. The NDVI trend strengthened and the range of values generally increased compared to the previous period. It is possible to change the composition of the vegetation type (e.g., barren to grass or grass to shrub) in the tundra region. It also suggests a possible lengthening of the growing season, as the change points in the NDVI data were observed in June and September. In other words, the change points are likely periods when various changes in tundra vegetation may have occurred.

Our analysis focused on the June–September period. This is crucial for understanding the changes in the length of the growing season of Arctic vegetation, which typically spans May–September [16]. These months are the periods associated with the onset and offset of growth, respectively. These two factors are essential when analyzing the length of the vegetation growth season. An extended growing season not only indicates a longer thaw period for the tundra but also impacts permafrost dynamics. Our study specifically investigated NDVI variations before and after the change point, uncovering notable differences, particularly in the June and September NDVI data for the critical growing season. An anomaly analysis revealed that Arctic greening is not merely continuing but accelerating at a rate of 2% per decade. Therefore, it simultaneously exhibits a change point. Overall, our results indicate accelerated Arctic greening over time. This is supported by the increasing annual and anomalous trends in the tundra NDVI.

All variables were highly correlated in the period before the change point, similar to the results of previous studies [2,14,21,25]. It has also been reported that the correlation between vegetation and temperature has decreased in recent years [26]. Our results indicate similar levels of correlations when compared with recent temperature correlations (1982–2019), implying the other correlation values may reflect the same changing conditions adequately. This explains why the correlations for the entire summer period were lower than in previous studies and suggests that there are external factors other than temperature and energy variables after the change point. Additionally, the summer Arctic Bowen ratio entered an unstable phase after 2012, contributing to the observed decline in the correlation between vegetation and climate variables (energy and temperature).

The analysis of time-lagged correlations revealed that all key climate variables exhibit their highest correlations with NDVI in the preceding month. This suggests that climatic variables influence vegetation, with NDVI responding accordingly. Therefore, surface atmosphere energy changes appear to significantly impact vegetation changes, highlighting this relationship as a critical parameter for unraveling the mechanisms driving Arctic vegetation dynamics.

6. Conclusions

This study specifically focused on tundra vegetation, aiming to elucidate its relationship with various climatic variables. Notably, a marked acceleration in Arctic greening was observed around 2012–2013, which can be identified as a potential tipping point. During this period, the correlation between vegetation, temperature, and energy diminished by half. This trend implies that factors beyond direct terrestrial influences are likely playing a significant role in driving tundra vegetation dynamics. Additionally, this timeframe coincides with heightened variability in vegetation and climate factors, suggesting a critical juncture for regional changes in the Arctic tundra. Although it is challenging to estimate the exact growing season length from monthly averaged data, our results indicate that the NDVIs of months including the start and end points increased more than in the other periods, implying an indirect possibility of a faster onset and later offset of growth. It was confirmed that the tundra vegetation is becoming greener, and the growing season is increasing.

However, this study’s reliance on observational data (satellite-derived NDVI) and reanalysis data (climatic variables) introduces complexities in establishing a direct causal link between climatic factors and vegetation changes. While our analysis presents an intuitive relationship between these variables, further research incorporating climate modeling is essential for unraveling the underlying mechanisms. A comprehensive understanding of the forces driving vegetation changes necessitates an analysis that extends beyond terrestrial variables to include a broader energy budget, which has emerged as a critical factor in vegetation dynamics. Additionally, examining the horizontal energy transfer between land and ocean and the vertical energy dynamics within the atmosphere, ocean, and sea ice is pivotal to fully comprehending these complex climate feedback processes.