Climate-Associated Seasonal Asymmetry in Forest Greening of Changbai Mountain

1. Introduction

Mountain forests are particularly sensitive to climate variability because topography compresses gradients in temperature, moisture, snow conditions, and growing-season length over short distances [1,2,3]. Satellite observations have documented widespread vegetation greening [4], but regional-scale trends provide limited insight into how forest greenness varies among seasons and elevation zones [5,6]. Along mountain slopes, the same climatic anomaly may occur at different phenological stages, resulting in contrasting greenness responses across elevations [3,5,7].

Growing-season averages are widely used for vegetation monitoring, but they can conceal marked seasonal differences [6]. Spring canopy development is commonly associated with heat accumulation, snowmelt [2,3,8], and soil thaw, whereas autumn greenness may reflect delayed senescence or prolonged late-season canopy activity [2,6]. Distinguishing these processes therefore requires a seasonal perspective.

Elevation introduces further spatial heterogeneity by influencing local climate, snow persistence, phenology, forest composition, and human accessibility. Elevation-dependent vegetation responses have been reported across mountain systems [5]. Changbai Mountain is characterised by pronounced vertical forest zonation [9] and climate-sensitive forest dynamics [10,11,12], making elevation-stratified analysis essential for detecting spatial differences in forest greenness and its climatic associations.

Climate–NDVI relationships may also change over time. Correlations calculated over an entire study period describe the average association between variables but may obscure shorter-term weakening, strengthening, or reversals under changing climatic conditions [13]. Sliding-window correlations can reveal such temporal variation [13], although estimates based on short, overlapping windows must be interpreted cautiously.

Nighttime light (NTL) data provide a spatial indicator of illuminated human and socioeconomic activity [14,15,16]. In forested mountain regions, NTL may therefore help identify human-activity signals concentrated at lower elevations. However, raw NTL and NDVI time series may appear correlated simply because both exhibit long-term trends [17]. Detrending is therefore needed to determine whether NTL–NDVI covariation persists beyond their shared temporal trends.

Changbai Mountain combines cold-temperate forests, strong elevation gradients, seasonal snow processes, and varying lower-elevation accessibility. Previous studies have examined seasonal vegetation dynamics across Northeast China [6]and vegetation dynamics, forest phenology, and climate associations specifically in Changbai Mountain [3,8]. Building on this work, we combine a fixed forest mask derived from the 2024 Annual China Land Cover Dataset (CLCD) with seasonal and elevation-stratified analyses, sliding-window spring climate correlations, and detrended NTL–NDVI regressions for 2000–2024. This integrated framework is used to compare seasonal, elevational, and temporal patterns in forest greenness rather than to attribute them to specific causal mechanisms.

Specifically, we address four questions: (1) How did forest NDVI change during the growing season, spring, summer, and autumn from 2000 to 2024? (2) How did the associations of NDVI with temperature and precipitation vary among seasons and elevation zones? (3) Did spring climate–NDVI associations vary over time? (4) Did NTL–NDVI covariation persist after shared long-term trends were removed?

2. Materials and Methods

2.1. Study Area

The study area is located within the Changbai Mountain Protection and Development Zone in northeastern China and covers parts of Antu County, Fusong County, and Changbai Korean Autonomous County (Figure 1). It covers approximately 3500 km² and extends from 127°28′ to 128°16′ E and from 41°42′ to 42°25′ N. Elevation ranges from approximately 600 to 2700 m, creating pronounced gradients in temperature, moisture, snow conditions, and phenology over short distances.

The region experiences a temperate monsoon climate, with most annual precipitation occurring from June to September. Forest vegetation exhibits clear vertical zonation. Korean pine–broad-leaved mixed forest predominates at lower elevations, dark coniferous forest occupies the middle-elevation belt, and Erman’s birch forest occurs near the upper forest limit, with alpine tundra above the forest belt. The thresholds of 1100 and 1700 m broadly correspond to the major transitions among the Korean pine–broad-leaved mixed forest, dark coniferous forest, and Erman’s birch forest belts [9]. These elevational gradients in climate and vegetation enable examination of spatial variation in forest greenness and its climatic associations. For the elevation-stratified analyses, forest pixels were classified into three zones: low (<1100 m), middle (1100–1700 m), and high (>1700 m).

2.2. Data Sources

We compiled vegetation index, climate, land-cover, nighttime light (NTL), and elevation data for the period 2000–2024 (Table 1).

Vegetation data were obtained from the MODIS/Terra Vegetation Indices Monthly L3 Global 1 km product (MOD13A3, Version 6.1) via Google Earth Engine [18]. This product provides monthly normalized difference vegetation index (NDVI) observations at 1 km spatial resolution.

Monthly mean temperature and precipitation data for 2000–2024 were obtained from the 1-km Monthly Mean Temperature Dataset for China and the 1-km Monthly Precipitation Dataset for China, respectively. Both datasets were provided by the National Tibetan Plateau/Third Pole Environment Data Center [19,20].

Annual NTL data for 2000–2024 were obtained from the Extended VIIRS-Like Artificial Nighttime Light Dataset of China, provided by the National Tibetan Plateau/Third Pole Environment Data Center [21]. The annual NTL rasters were aligned with the common analysis grid before mean values were calculated for each elevation zone.

Land-cover information was obtained from the 2024 Annual China Land Cover Dataset (CLCD) v01, Version 1.0.4, developed by Jie Yang and Xin Huang at Wuhan University and distributed through Zenodo [22]. The CLCD was derived from Landsat surface reflectance using random forest classification, spatiotemporal filtering, and logical post-processing. Forest was represented by class code 2, and the 2024 raster covering the study area was used to construct the forest-cover mask.

Elevation data were obtained from the ASTER Global Digital Elevation Model Version 3 (ASTER GDEM V003), provided by the Geospatial Data Cloud [23]. The dataset has a spatial resolution of 30 m and was used to derive elevation information.

2.3. Forest Mask and Spatial Harmonisation

The 2024 CLCD raster used to construct the forest mask was supplied in an Albers Equal-Area Conic projection based on the WGS 1984 datum.

Forest cover was quantified from source-pixel areas rather than by assigning a single land-cover class to each MODIS pixel through majority-rule or nearest-neighbour resampling. Forest pixels were first identified at the native 30 m CLCD resolution. For each source pixel, forest area and valid classified area were calculated separately and then aggregated independently to the MODIS NDVI grid by area summation. The forest-cover proportion of each MODIS pixel was calculated as

$$F_i={\displaystyle \frac{\sum _{j\in S_i}{I}_j{A}_j}{\sum _{j\in S_i}{V}_j{A}_j}},$$

(1)

where $F_i$is the forest-cover proportion of MODIS pixel $i$; $S_i$is the set of 30 m source pixels assigned to MODIS pixel $i$during spatial aggregation; $I_j$is the forest indicator for source pixel $j$, with forest pixels assigned a value of 1 and all other valid land-cover classes assigned a value of 0; $V_j$is the valid-classification indicator, with valid land-cover pixels assigned a value of 1 and NoData pixels assigned a value of 0; and $A_j$is the actual area of source pixel $j$.

The forest mask was defined as

$${\mathrm{M}}_{\mathrm{i}}=\{\begin{array}{cc}1,& {\mathrm{F}}_{\mathrm{i}}\ge 0.70,\\ \mathrm{N}\mathrm{A},& {\mathrm{F}}_{\mathrm{i}}<0.70,\end{array}$$

(2)

where ${\mathrm{M}}_{\mathrm{i}}=1$denotes a retained forest pixel and ${\mathrm{M}}_{\mathrm{i}}=\mathrm{N}\mathrm{A}$denotes a pixel excluded from subsequent analyses. The final mask therefore retained MODIS pixels with at least 70% forest cover according to the 2024 CLCD. This fixed mask was applied consistently to all NDVI, climate, and NTL analyses from 2000 to 2024, thereby maintaining a constant spatial sample throughout the study period.

The 70% threshold was selected by comparing candidate forest masks within the study area. A 50% threshold retained 4,627 MODIS pixels and 3,433.53 km² of CLCD-derived forest area, whereas the 70% threshold retained 4,568 pixels and 3,389.78 km². Increasing the threshold removed 59 pixels and 43.75 km² of forest area while retaining 98.73% of the pixels identified under the 50% criterion. The 70% threshold therefore increased the minimum forest proportion of retained pixels with limited loss of candidate forest area. It was treated as a study-specific screening criterion rather than a universal remote-sensing standard or statistical threshold.

After assignment to the elevation zones, the final elevation-stratified sample comprised 4,344 forest pixels and 3,223.62 km² of CLCD-derived forest area (Table 2). Low- and middle-elevation zones accounted for 47.64% and 44.54% of the masked forest area, respectively, whereas the high-elevation zone accounted for 7.82%.

MODIS NDVI defined the common analysis grid and was not resampled. The forest mask and DEM-derived elevation zones, both categorical layers, were aligned to the MODIS grid using nearest-neighbour assignment. Temperature and precipitation, as continuous variables, were aligned using bilinear interpolation. Annual NTL rasters were aligned to the same analysis grid before zonal means were calculated.

Area-weighted means were calculated for the entire masked forest area and for each elevation zone as

$${\bar{X}}_{y,z}={\displaystyle \frac{\sum _{i\in Z_z}{X}_{i,y}{A}_i}{\sum _{i\in Z_z}{A}_i}},$$

(3)

where ${\bar{X}}_{y,z}$is the area-weighted mean for year $y$and spatial zone $z$; $Z_z$is the set of valid forest pixels within zone $z$; $X_{i,y}$is the temperature, precipitation, or NTL value of pixel $i$in year $y$; and $A_i$is the actual area of pixel $i$. The spatial zone $z$represents either the entire masked forest area or one of the three elevation zones. Area weighting was used to account for differences in the actual area represented by geographic-grid pixels at different latitudes.

2.4. Seasonal Aggregation and Data Completeness

Seasonal periods were defined before statistical analysis as follows: growing season, May–September; spring, April–May; summer, June–August; and autumn, September–October. The growing-season metric was analysed separately from the three seasonal periods; therefore, these temporal windows were not mutually exclusive.

For year $y$and seasonal period $s$, seasonal mean NDVI was calculated as

$$\left({\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}}_{y,s}={\displaystyle \frac{1}{n_s}}{\sum }_{m\in s}{\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}}_{y,m}\right)$$

(4)

where ${\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}}_{y,s}$is the mean NDVI for seasonal period $s$in year $y$, ${\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}}_{y,m}$is the monthly NDVI in month $m$, and $n_s$is the number of months included in seasonal period $s$.

Seasonal mean temperature was calculated as

$$\left(T_{y,s}={\displaystyle \frac{1}{n_s}}\sum _{m\in s}{T}_{y,m}\right),$$

(5)

where $T_{y,s}$is the mean temperature for seasonal period $s$in year $y$, and $T_{y,m}$is the monthly mean temperature in month $m$.

Seasonal precipitation was calculated as the sum of monthly precipitation:

$$P_{y,s}=\sum _{m\in s}{P}_{y,m},$$

(6)

where $P_{y,s}$is the total precipitation for seasonal period $s$in year $y$, and $P_{y,m}$is the monthly total precipitation in month $m$. Annual NTL data required no temporal aggregation.

The fixed 2024 CLCD-derived forest mask and the same elevation zones were applied throughout 2000–2024. Seasonal values were calculated only when observations were available for all constituent months. The final dataset contained 25 complete annual observations for each seasonal period and spatial unit. Spatial means for the entire masked forest area and for each elevation zone were calculated using the area-weighting procedure described in Section 2.3.

2.5. Statistical Analyses

2.5.1. Seasonal NDVI Trends

Long-term trends in seasonal forest NDVI from 2000 to 2024 were estimated using ordinary least-squares regression:

$${\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}}_{y,s}={\alpha }_s+{\beta }_s{t}_y+{\epsilon }_{y,s},$$

(7)

where ${\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}}_{y,s}$is the mean NDVI for seasonal period $s$in year $y$, $t_y$is the corresponding calendar year, ${\alpha }_s$is the intercept, ${\beta }_s$is the annual rate of NDVI change, and ${\epsilon }_{y,s}$is the residual error. Trend direction and magnitude were evaluated using the regression slope, coefficient of determination ($R^2$), and $P$-value for ${\beta }_s$.

2.5.2. NDVI–Climate Correlations and Standardised Regression

Pearson correlation coefficients were calculated between seasonal NDVI and concurrent seasonal mean temperature and between seasonal NDVI and concurrent seasonal total precipitation. Correlations were first calculated for the entire masked forest area and then separately for the low-, middle-, and high-elevation zones.

For two annual variables $x$and $y$, the Pearson correlation coefficient was calculated as

$$r={\displaystyle \frac{{\displaystyle {\sum }_{i=1}^n}\left(x_i-\bar{x}\right)\left(y_i\right|\bar{y})}{\sqrt{{\displaystyle {\sum }_{i=1}^n}{\left(x_i\right|\bar{x})}^2}\sqrt{{\displaystyle {\sum }_{i=1}^n}{\left(y_i\right|\bar{y})}^2}}},$$

(8)

where $x_i$and $y_i$are the paired NDVI and climate observations in year $i$, $\bar{x}$and $\bar{y}$are their respective means, and $n$is the number of paired observations. Each full-period correlation was based on 25 annual observations. Correlation coefficients, two-sided $P$-values, and sample sizes were reported.

To compare the relative magnitudes of the temperature and precipitation coefficients within each seasonal period and elevation zone, NDVI, temperature, and precipitation were standardised as

$$z\left(X\right)={\displaystyle \frac{X-\bar{X}}{\mathrm{S}\mathrm{D}\left(X\right)}},$$

(9)

where $X$denotes NDVI, temperature, or precipitation, $\bar{X}$is its mean, and $\mathrm{S}\mathrm{D}\left(X\right)$is its standard deviation.

A standardised multiple linear regression model was then fitted:

$$z\left(\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}\right)=\alpha +{\beta }_{T}z\left(T\right)+{\beta }_{P}z\left(P\right)+\epsilon ,$$

(10)

where $T$and $P$denote temperature and precipitation, respectively; $\alpha$is the intercept; ${\beta }_T$and ${\beta }_P$are the standardised regression coefficients for temperature and precipitation; and $\epsilon$is the residual error.

The relative coefficient share of temperature was calculated as

$$C_T={\displaystyle \frac{\mid {\beta }_T\mid }{\mid {\beta }_T\mid +\mid {\beta }_P\mid }}\times 100\mathrm{\%},$$

(11)

whereas the relative coefficient share of precipitation was calculated as

$$C_P={\displaystyle \frac{\mid {\beta }_P\mid }{\mid {\beta }_T\mid +\mid {\beta }_P\mid }}\times 100\mathrm{\%}.$$

(12)

The signs of ${\beta }_T$and ${\beta }_P$were retained when interpreting the direction of each relationship, whereas their absolute values were used only to calculate the relative coefficient shares. These percentages describe the relative magnitudes of the two standardised coefficients within the same model and do not represent independent variance decomposition or causal contributions.

Multicollinearity between temperature and precipitation was assessed using the variance inflation factor:

$$\mathrm{V}\mathrm{I}\mathrm{F}={\displaystyle \frac{1}{1-r_{TP}^2}},$$

(13)

where $r_{TP}$is the Pearson correlation coefficient between temperature and precipitation. VIF values below 5 were considered to indicate no problematic multicollinearity[24].

2.5.3. Sliding-Window Correlations

Temporal variation in NDVI–climate correlations was examined using 5-year sliding windows for each seasonal period and elevation zone. For a window beginning in year $t$, the included annual observations were defined as

$$W_t=\left\{t,t+1,t+2,t+3,t+4\right\},$$

(14)

and the corresponding centre year was calculated as

$$c_t=t+2.$$

(15)

The first window covered 2000–2004 and was represented by the centre year 2002, whereas the final window covered 2020–2024 and was represented by 2022. This procedure produced 21 overlapping windows.

Within each window, Pearson correlation coefficients were calculated between NDVI and temperature and between NDVI and precipitation. A correlation coefficient was calculated only when all five paired annual observations were available.

Because each window contained only five observations and adjacent windows shared four years, window-level $P$-values were interpreted as exploratory. The sliding-window analysis was therefore used to describe temporal variation in correlation strength and direction rather than to identify statistically independent transitions or precise ecological thresholds.

2.5.4. NTL–NDVI Regression and Detrending

The NTL analysis was exploratory and was designed to assess whether raw NTL–NDVI covariation persisted after common linear temporal trends had been removed. Annual mean NTL and growing-season NDVI were summarised within the same forest mask and elevation zones. For each elevation zone, the raw NTL–NDVI relationship was evaluated using ordinary least-squares regression:

$${\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}}_y=\alpha +\gamma {\mathrm{N}\mathrm{T}\mathrm{L}}_y+{\epsilon }_y,$$

(16)

where ${\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}}_y$is growing-season mean NDVI in year $y$, ${\mathrm{N}\mathrm{T}\mathrm{L}}_y$is annual mean nighttime light intensity in the same year and elevation zone, $\alpha$is the intercept, $\gamma$is the regression coefficient, and ${\epsilon }_y$is the residual error. Each elevation-specific regression was based on 25 annual observations.

To compare interannual variation in NDVI and NTL on a common scale, both variables were separately transformed to $z$-scores within each elevation zone using Equation (9). The resulting standardised values were used to construct annual time series.

Because both NDVI and NTL could contain long-term temporal trends, detrended analyses were also conducted. A linear temporal trend was first fitted to NDVI:

$${\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}}_y=a+{\beta }_{\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}}{t}_y+e_{\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I},y},$$

(17)

where $a$is the intercept, ${\beta }_{\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I}}$is the temporal trend coefficient, and $e_{\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I},y}$is the residual after removal of the linear NDVI trend.

A corresponding temporal trend was fitted to NTL:

$${\mathrm{N}\mathrm{T}\mathrm{L}}_y=c+{\beta }_{\mathrm{N}\mathrm{T}\mathrm{L}}{t}_y+e_{\mathrm{N}\mathrm{T}\mathrm{L},y},$$

(18)

where $c$is the intercept, ${\beta }_{\mathrm{N}\mathrm{T}\mathrm{L}}$is the temporal trend coefficient, and $e_{\mathrm{N}\mathrm{T}\mathrm{L},y}$is the residual after removal of the linear NTL trend.

The relationship between the two detrended residual series was subsequently evaluated using

$$e_{\mathrm{N}\mathrm{D}\mathrm{V}\mathrm{I},y}=\alpha +\delta e_{\mathrm{N}\mathrm{T}\mathrm{L},y}+{\epsilon }_y,$$

(19)

where $\delta$is the regression coefficient between detrended NTL and detrended NDVI. Comparison of the raw and detrended regressions was used to assess whether the NTL–NDVI relationship persisted after linear temporal trends had been removed.

All statistical analyses were conducted in R version [version number]. Unless otherwise stated, statistical tests were two-sided, and $P<0.05$was considered statistically significant.

3. Results

3.1. Seasonal NDVI Trends from 2000 to 2024

From 2000 to 2024, forest NDVI showed positive linear trends in all four seasonal periods, but the magnitude and statistical significance of these trends differed substantially among seasons (Figure 2). Autumn exhibited the strongest increase, with both the steepest slope and the highest $R^2$among the four periods (slope = 0.0019 yr⁻¹, $R^2$= 0.280, $P$< 0.01). Significant but weaker increases were observed during the growing season (slope = 0.0006 yr⁻¹, $R^2$= 0.166, $P$< 0.05) and summer (slope = 0.0006 yr⁻¹, $R^2$= 0.188, $P$< 0.05). Spring NDVI showed a positive but non-significant trend (slope = 0.0011 yr⁻¹, $R^2$= 0.087, $P$= 0.152).

Spring NDVI displayed substantial interannual variability, as reflected by the low $R^2$of its linear trend. Growing-season and summer NDVI generally increased over the study period, although both declined markedly in 2024. Autumn NDVI also fluctuated among years, with a notable peak in 2010, but exhibited a more pronounced long-term increase than the other seasonal series. Overall, forest greening was strongly season-dependent, with the strongest long-term increase in autumn and a weak, non-significant trend in spring.

3.2. Seasonal NDVI Correlations with Temperature and Precipitation

NDVI correlations with temperature and precipitation differed substantially among seasonal periods (Figure 3). The strongest correlation was observed in spring, when NDVI was strongly and positively correlated with temperature ($r=0.760$, $P=1.04\times {10}^{-5}$, $n=25$). In contrast, the correlation between spring NDVI and precipitation was negative but not statistically significant ($r=-0.299$, $P=0.146$). Growing-season NDVI was also positively correlated with temperature ($r=0.446$, $P=0.025$), whereas its correlation with precipitation was weak and non-significant ($r=0.186$, $P=0.375$).

NDVI correlations with both climate variables were generally weak during summer and autumn. In summer, the correlation with temperature was close to zero ($r=-0.006$, $P=0.978$), while the correlation with precipitation was weak and positive ($r=0.074$, $P=0.725$). In autumn, weak positive correlations were observed between NDVI and precipitation ($r=0.192$, $P=0.358$) and between NDVI and temperature ($r=0.182$, $P=0.383$); neither correlation was statistically significant. Overall, significant study-area-wide correlations were detected only between NDVI and temperature during spring and the growing season, whereas precipitation was not significantly correlated with NDVI in any seasonal period.

3.3. Elevation-Dependent Seasonal NDVI–Climate Correlations and Standardised Coefficient Shares

Elevation-stratified analyses revealed substantial seasonal variation in both NDVI–climate correlations and standardised regression coefficients (Figure 4 and Figure 5). In spring, NDVI was positively correlated with temperature across all three elevation zones. The correlations were significant at low elevation ($r=0.543$, $P=0.005$), middle elevation ($r=0.846$, $P=9.98\times {10}^{-8}$), and high elevation ($r=0.769$, $P=6.95\times {10}^{-6}$), with the strongest correlation at middle elevation. In contrast, spring NDVI was negatively correlated with precipitation at low ($r=-0.188$, $P=0.369$), middle ($r=-0.343$, $P=0.094$), and high elevations ($r=-0.390$, $P=0.054$), although none of these correlations was statistically significant.

During the growing season, NDVI correlations with temperature were positive but weakened with increasing elevation. The correlation was significant at low elevation ($r=0.458$, $P=0.021$), but did not reach statistical significance at middle ($r=0.384$, $P=0.058$) or high elevation ($r=0.153$, $P=0.465$). Growing-season correlations with precipitation were weak and positive across all elevation zones ($r=0.047$–0.237; all $P>0.05$). No significant correlations with either temperature or precipitation were detected during summer or autumn. In summer, temperature correlations ranged from weakly negative to weakly positive ($r=-0.086$–0.057), whereas precipitation correlations were weakly positive ($r=0.043$–0.129). In autumn, weak positive correlations were observed for temperature ($r=0.096$–0.235) and precipitation ($r=0.162$–0.187), with all $P>0.05$.

The standardised multiple regression models showed little multicollinearity between temperature and precipitation (all VIFs ≤ 1.22). In spring, temperature represented 97.8%, 93.6%, and 84.1% of the combined absolute standardised coefficient magnitude at low, middle, and high elevations, respectively, whereas the corresponding precipitation shares were 2.2%, 6.4%, and 15.9% (Figure 5). The spring regression models were statistically significant at low ($R^2=0.295$, $P=0.021$), middle ($R^2=0.719$, $P=8.74\times {10}^{-7}$), and high elevations ($R^2=0.608$, $P=3.32\times {10}^{-5}$).

During the growing season, temperature had the larger absolute standardised coefficient share at all three elevations, with shares of 89.5%, 85.0%, and 92.3% at low, middle, and high elevations, respectively. In summer, precipitation had the larger coefficient share at low elevation (57.1%) and high elevation (59.1%), whereas temperature had a slightly larger share at middle elevation (53.6%). Autumn coefficient shares were more balanced. Temperature shares were 53.0%, 42.6%, and 52.6% at low, middle, and high elevations, respectively, and the corresponding precipitation shares were 47.0%, 57.4%, and 47.4%.

Together, the correlation and regression results showed that the spring NDVI–temperature relationship was consistently strong across elevation zones. The growing-season temperature relationship was weaker and varied more strongly with elevation, whereas summer and autumn correlations remained weak despite seasonal shifts in relative magnitudes of the standardised temperature and precipitation coefficients. The reported percentages represent relative shares of the absolute standardised coefficients within each model and should not be interpreted as independent contributions to explained variance or as causal effects.

3.4. Temporal Variation in Spring NDVI Correlations with Temperature and Precipitation

The 5-year sliding-window analysis showed temporal variation in spring NDVI correlations with temperature and precipitation across elevation zones (Figure 6). Spring NDVI–temperature correlations were predominantly positive, but their magnitude varied among elevations and windows. Positive correlations were generally stronger at middle and high elevations than at low elevation. At high elevation, the temperature correlation weakened in several windows centred around the mid-2010s.

Spring NDVI–precipitation correlations varied more strongly in both direction and magnitude. At low and middle elevations, correlations were mainly negative in earlier windows but moved towards neutral or positive values later in the record. At high elevation, precipitation correlations were predominantly negative during the early and middle portions of the record and became more positive towards the end of the series.

Overall, temperature correlations remained more consistently positive than precipitation correlations, particularly at middle and high elevations. Precipitation correlations showed more frequent changes in direction. Because each window contained only five paired annual observations and adjacent windows shared four years, the window-level P-values were treated as nominal and descriptive. The patterns are therefore interpreted as exploratory temporal variation rather than as statistically independent transitions or precisely dated ecological changes.

3.5. Elevation-Dependent Relationships between Nighttime Light Intensity and Growing-Season NDVI

The raw relationship between annual mean nighttime light intensity and growing-season mean NDVI differed among elevation zones (Figure 7). At low elevation, the relationship was positive and statistically significant, with the highest $R^2$among the three zones (slope = 0.159, $R^2=0.198$, $P=0.026$, $n=25$). At middle elevation, the relationship was also positive but not statistically significant (slope = 0.438, $R^2=0.049$, $P=0.288$). At high elevation, the relationship was negative and non-significant, with negligible variance explained (slope = −0.514, $R^2=0.002$, $P=0.853$).

The standardised time series showed an overall increase in nighttime light intensity at low elevation, particularly during the latter part of the study period, accompanied by a weaker increase in growing-season NDVI. At middle elevation, both variables showed more variable interannual patterns, whereas nighttime light intensity remained low at high elevation. After linear temporal trends were removed from both variables, the low-elevation residual relationship was no longer statistically significant (slope = 0.154, $R^2=0.042$, $P=0.324$; Figure 8). Detrended residual regressions were also non-significant at middle elevation (slope = −0.181, $R^2=0.004$, $P=0.774$) and high elevation (slope = −1.028, $R^2=0.007$, $P=0.697$).

Overall, a significant raw low-elevation NTL–NDVI relationship was detected, but it did not persist after detrending, and no significant residual relationship was found in any elevation zone.

4. Discussion

4.1. Autumn-Dominated Seasonal Asymmetry in Forest Greening

Autumn is a phenologically sensitive period in which relatively small changes in the timing or rate of canopy decline can produce clear NDVI signals. Summer canopies in Changbai Mountain are already close to their seasonal maximum, and NDVI may become less responsive to further increases in leaf area or canopy density. During autumn, chlorophyll degradation, leaf loss, and declining green leaf area create a steeper seasonal trajectory. Slower canopy decline or prolonged leaf retention would therefore be more readily detected in autumn than during peak summer conditions. The strong autumn trend observed here is consistent with a longer period of canopy greenness late in the growing season.

Earlier studies provide regional support for this result. Hou et al. [10] documented increasing vegetation greenness during 2000–2009 and reported closer NDVI relationships with temperature than with precipitation, particularly in spring and autumn. Guo et al. [11] also identified widespread regional greening, including a significant autumn increase and a positive spring temperature–NDVI relationship. Extending the record to 2024 indicates that the autumn signal persisted beyond the shorter periods examined previously. The consistent 70% forest-cover mask also narrows the analysis to forest-dominated pixels and reduces interference from agricultural land, grassland, settlements, and other non-forest surfaces.

The pronounced autumn trend was not accompanied by significant correlations with concurrent autumn temperature or precipitation. Year-to-year variation in September–October climate did not closely track the long-term rise in NDVI. Autumn canopy greenness may integrate antecedent climate, frost and snowfall timing, and gradual changes in forest structure. Stand development and forest recovery could also increase late-season NDVI without generating strong correlations with annual autumn climate. Concurrent seasonal means alone are unlikely to capture these cumulative, lagged, and structural influences.

Autumn NDVI represents late-season spectral greenness, not an equivalent change in ecosystem productivity or physiological senescence[25]. NDVI responds to green leaf area, pigment status, and canopy structure, but it does not necessarily track photosynthetic carbon uptake. The pattern agrees with reports of prolonged canopy greenness and delayed canopy decline in other temperate forests [26,27], although broad-band NDVI cannot resolve the underlying physiological process[28]. The autumn signal nevertheless shows that late-season greenness contributed substantially to the overall seasonal pattern in Changbai Mountain.

4.2. Spring NDVI–Temperature Correlations and Early-Season Canopy Development

The close spring NDVI–temperature correlation fits the phenological setting of Changbai Mountain, where early canopy development is linked to heat accumulation, snowmelt, soil thaw, and leaf emergence[29]. Despite this strong correlation, spring NDVI did not increase significantly over the full study period. Temperature was closely related to year-to-year variation in spring greenness, while the direction of long-term change remained weak[30]. Late frost[31], variable snow cover, cloudiness, and differences in soil thaw may have interrupted a steady increase.

Previous work in Changbai Mountain supports this interpretation. Hou et al. [10] found that spring NDVI was more closely related to temperature than to precipitation. Local phenological studies have also linked growing-season onset to spring temperature, snow-cover duration, and snowmelt timing [3,7,8]. Green-up generally occurs later, and the growing season becomes shorter, at higher elevations [5]. In the present study, the strongest temperature correlation occurred at middle elevation, followed by high and low elevations. The middle-elevation zone overlaps the dark coniferous forest belt, but our analysis does not allow us to isolate forest composition from snow persistence, terrain exposure, or local climate. We therefore interpret the stronger middle-elevation correlation as an elevation-zone pattern, not as a forest-type response. Elevation alone therefore does not determine the strength of the spring relationship[32]. Forest composition, snow persistence, terrain exposure, and local moisture conditions may also influence canopy development.

Spring precipitation was negatively correlated with NDVI in all three elevation zones, although none of these correlations was statistically significant. The negative direction may partly reflect the form and timing of precipitation. At higher elevations, some precipitation may fall as snow and prolong snow cover. Wet spring conditions may also coincide with cooler temperatures, lower radiation, and later soil thaw. Seasonal precipitation totals cannot distinguish rainfall from snowfall or describe the timing of snowmelt. The negative correlations should therefore not be taken as evidence that precipitation suppressed spring forest greenness.

4.3. Seasonal Shifts in Temperature- and Precipitation-Related NDVI Variation

The regression pattern changed after spring canopy development. In summer, temperature correlations approached zero, while precipitation had the larger standardised coefficient at low and high elevations. Temperature retained a slightly larger coefficient at middle elevation. By autumn, the two variables showed more similar coefficient magnitudes, although neither was significantly correlated with NDVI.

Greater summer water demand offers one possible explanation. Peak canopy development is accompanied by high transpiration, and rainfall can help maintain soil moisture during warm periods. This may be especially relevant at low elevations, where summer temperatures are higher, and at exposed high-elevation sites, where shallow soils may retain less water. Summer precipitation, however, showed little correlation with annual NDVI variation. The larger precipitation coefficients are therefore compatible with greater summer water relevance, but they do not demonstrate water limitation.

Delayed precipitation effects may also contribute to the seasonal contrast. Moisture supplied during spring may influence soil-water availability and canopy conditions later in summer. Previous studies have reported NDVI responses to precipitation after delays of several weeks or months [33,34]. Because the models paired NDVI with climate from the same seasonal period, they could not test this cross-seasonal pathway.

The reported percentages represent the relative magnitudes of the absolute standardised coefficients for temperature and precipitation. They are not independent shares of explained variance or measures of causal climatic effects. The summer results indicate a change in the relative regression pattern, not a confirmed transition from temperature limitation to water limitation.

4.4. Temporal Variation in Spring NDVI–Climate Correlations

The sliding-window analysis extends the full-period correlations by revealing how spring NDVI–climate relationships varied over time. Spring NDVI remained positively correlated with temperature in most windows, especially at middle and high elevations. The relationship weakened during several windows centred in the mid-2010s, most clearly at high elevation. Precipitation correlations changed more markedly in both direction and magnitude. They were mainly negative in earlier windows but moved towards neutral or positive values later in the record. Positive precipitation correlations appeared in later windows at low and middle elevations and in the final window at high elevation.

Similar temporal changes in vegetation–climate relationships have been reported under varying climatic conditions [13]. Year-to-year changes in snow conditions, freeze–thaw timing, radiation, and frost occurrence may alter the strength of the spring temperature signal [3,8]. NDVI may track temperature closely in colder years, when canopy development remains strongly constrained by thermal and snow conditions. In milder years, radiation, moisture, frost timing, and local forest conditions may account for more of the observed variation. The available data do not allow us to separate these processes.

Changes in precipitation correlations may reflect differences in rainfall–snowfall partitioning and antecedent moisture. In warmer springs, a larger proportion of precipitation may fall as rain or enter the soil earlier. Water stored in snow or soil can also affect vegetation later through snowmelt and soil-moisture carry-over [8,33,34]. Same-season correlations cannot distinguish these immediate and delayed pathways.

Unusual climate years may further amplify variation among windows. Late frost, exceptionally warm springs, drought, and heavy precipitation can alter green-up and canopy greenness [35]. Their influence is especially large in a 5-year window: each estimate contains only five observations, and adjacent windows share four. A single anomalous year can affect several consecutive estimates. The window patterns are therefore treated as exploratory evidence of temporal variation rather than ecological thresholds.

4.5. Methodological Considerations and Future Research

The analysis was restricted to pixels with at least 70% forest cover according to the fixed 2024 CLCD mask. This maintained a consistent spatial sample across years and reduced interference from non-forest surfaces. It does not establish that every retained pixel remained forest throughout 2000–2024. The 1 km MOD13A3 grid also represents broad canopy conditions rather than individual forest stands. A single pixel may contain different forest types, terrain settings, and local disturbances. Higher-resolution imagery and annual land-cover records could help separate persistent forest areas from local vegetation changes.

Temperature and precipitation described the main seasonal climate conditions, but they did not cover all processes relevant to mountain forest phenology. Snow duration, rainfall–snowfall partitioning, radiation, frost, soil temperature, soil moisture, and atmospheric water demand may alter the timing and strength of NDVI–climate relationships. These processes are particularly relevant to spring canopy development, possible summer water effects, and autumn canopy decline. The mechanisms discussed in the preceding sections are supported by ecological reasoning and previous studies, but they were not tested directly in the present models.

Our statistical analyses quantify covariation, not causation. Pearson correlations describe synchronous year-to-year relationships, while the standardised regression coefficients indicate the relative magnitudes of temperature and precipitation within each fitted model. Their percentage shares are not independent partitions of explained variance. The 5-year sliding-window results require additional caution because each estimate contains only five observations and adjacent windows overlap by four years. Individual anomalous years may therefore influence several consecutive windows. Longer records and formal change-point or time-varying coefficient models would provide stronger evidence for persistent changes in climate sensitivity.

Interpreting nighttime light data also requires caution. NTL indicates the broad concentration of illuminated human activity, but it cannot distinguish roads, settlements, tourism, restoration, forest management, or changes in lighting technology. The significant low-elevation NTL–NDVI relationship was no longer significant after detrending, indicating that the raw result largely reflected parallel long-term change. Attributing the pattern to particular human processes would require spatially explicit data on land-use change, infrastructure, tourism activity, restoration, and management. These data should be analysed together with climate, forest structure, and disturbance history.

Elevation stratification revealed clear spatial contrasts, but elevation integrates several environmental and human gradients. Temperature, snow conditions, forest composition, terrain exposure, growing-season length, and accessibility all vary with altitude. Future analyses could refine the current elevation zones by incorporating forest type, slope, aspect, snow dynamics, and local disturbance. Nevertheless, separate assessment of spring canopy development, summer precipitation-related NDVI variation, and autumn canopy decline provides more informative insight into forest change than annual or growing-season averages alone.

5. Conclusions

Changbai Mountain forests showed pronounced seasonal asymmetry in greenness from 2000 to 2024. Autumn NDVI exhibited the strongest significant increase, whereas the spring trend was weak and non-significant. Despite the absence of a clear long-term spring trend, spring NDVI had the strongest correlation with temperature. This positive correlation was observed across all elevation zones and was strongest at middle elevation. By contrast, NDVI correlations with temperature and precipitation were generally weak during summer and autumn.

Spring NDVI–climate correlations also varied over time, although the 5-year sliding-window results remain exploratory because of the small and overlapping samples. At low elevation, the raw relationship between nighttime light intensity and growing-season NDVI was significant but no longer significant after detrending, indicating that the original relationship largely reflected parallel long-term trends. Overall, forest greening in Changbai Mountain cannot be represented adequately by a single growing-season trend or full-period climate correlation. Seasonal, elevational, and temporal differences should be considered when interpreting long-term changes in mountain forest greenness.