Linking Vegetation Phenology to Net Ecosystem Productivity: Climate Change Impacts in the Northern Hemisphere Using Satellite Data

Abstract

With global climate change, linking vegetation phenology with net ecosystem productivity (NEP) is crucial for assessing vegetation carbon storage capacity and predicting terrestrial ecosystem changes. However, there have been few studies investigating the relationship between vegetation phenology and NEP in the middle and high latitudes of the Northern Hemisphere. This study comprehensively analyzed vegetation phenological changes and their climate drivers using satellite data. It also investigated the spatial distribution and climate drivers of NEP and further analyzed the sensitivity of NEP to vegetation phenology. The results indicated that the average land surface phenology (LSP) was dominated by a monotonic trend in the study area. LSP derived from different satellite products and retrieval methods exhibited relatively consistent responses to climate. The average SOS and POS for different retrieval methods showed a higher negative correlation with nighttime temperatures compared to daytime temperatures. The average EOS exhibited a higher negative correlation with daytime temperatures than a positive correlation. The correlations between VPD and the average SOS, POS, and EOS showed that the proportion of negative correlations was higher than that of positive correlations. The average annual NEP ranged from 0 to 1000 gC·m⁻². The cumulative trends of NEP were mainly monotonically increasing, accounting for 61.04%, followed by monotonically decreasing trends, which accounted for 17.95%. In high-latitude regions, the proportion of positive correlation between VPD and NEP was predominant, while the proportion of negative correlation was predominant in middle-latitude regions. The positive and negative correlations between soil moisture and NEP (48.08% vs. 51.92%) were basically consistent in the study area. The correlation between SOS and POS with NEP was predominantly negative. The correlation between EOS and NEP was overall characterized by a greater proportion of negative correlations than positive correlations. The correlation between LOS and NEP exhibited a positive relationship in most areas. The sensitivity of NEP to vegetation phenological parameters (SOS, POS, and EOS) was negative, while the sensitivity of NEP to LOS was positive (0.75 gC·m⁻²/d for EVI vs. 0.63 gC·m⁻²/d for LAI vs. 0.30 gC·m⁻²/d for SIF). This study provides new insights and a theoretical basis for exploring the relationship between vegetation phenology and NEP under global climate change.

1. Introduction

Vegetation is a crucial component of terrestrial ecosystems and is highly responsive to climate change [1,2,3]. Vegetation phenology is considered a sensitive and precise indicator for assessing how terrestrial vegetation responds to climate change [4,5,6,7,8]. Terrestrial vegetation adapts to climate change by altering its cyclical phenological patterns (e.g., budburst, leaf unfolding, leaf greening, flowering, fruiting, leaf browning, leaf shedding), such as early leaf greening and delayed leaf browning [9,10,11,12,13,14,15,16]. Vegetation phenology signifies the timing of canopy photosynthesis, playing a pivotal role in vegetation’s absorption of carbon dioxide from the atmosphere and carbon accumulation [17,18]. Furthermore, the extension of the vegetation growth period will significantly change the vegetation transpiration and surface albedo patterns, thereby affecting the water cycle and surface energy flow processes [2,19,20,21,22]. As global temperature rise becomes a general consensus, there are forecasts indicating that the global temperature may increase by 2.0 °C to 4.9 °C by the early 22nd century, which will significantly alter the patterns of vegetation phenology [23]. Therefore, further enriching the understanding of vegetation phenology is of great significance for addressing the challenge of climate change.

The traditional method of vegetation phenology monitoring relies on manual visual assessment at specific observation sites, recording the phenophases of individual or multiple plants, such as the dates of the first visible leaf stalk, the first appearance of flower buds, the first ripe fruits, and the date of 50% autumnal leaf fall [24,25]. Although manual observation can provide reliable phenological records, it is difficult to obtain complete data on large-scale plant communities and to conduct a comprehensive analysis of the spatial pattern of vegetation phenology in large areas due to the limitations of observation sites, time, and labor costs [26]. In recent decades, with the rapid development of remote sensing technology, the monitoring scope of vegetation phenology has expanded from single-site scale to regional, hemispheric, and even global scales [27,28,29]. The remote sensing of land surface phenology (LSP) utilizes high temporal and spatial resolution satellite datasets to retrieve LSP based on developed algorithms, with specific phenological indicators including the start of the growing season (SOS), the peak of the growing season (POS), the end of the growing season (EOS), and the length of the growing season (LOS) [13,30]. Remote sensing vegetation indices such as the Normalized Difference Vegetation Index (NDVI) and the Enhanced Vegetation Index (EVI) are considered effective proxies for reflecting vegetation greenness [31,32,33] and are extensively used in vegetation phenology studies [34,35,36]. Compared with EVI, NDVI is oversaturated in high biomass areas and is more sensitive to the influence of the atmosphere and soil [37]. Biophysical variables, such as leaf area index (LAI) and solar-induced chlorophyll fluorescence (SIF), can more directly represent vegetation structural information and approximate photosynthetic activity, respectively, showing stronger sensitivity to high-density vegetation and the ability to capture vegetation productivity [38,39].

Research has indicated that compared to vegetation indices, measurements of phenology by LAI and SIF are more robust and closer to ground-based phenological observations [40,41]. Assessing the differences in LSP retrieved by EVI, LAI, and SIF helps deepen the understanding of vegetation phenological changes and their connections and distinctions in the middle and high latitudes of the Northern Hemisphere. A comprehensive comparative analysis of vegetation phenology based on EVI, LAI, and SIF can provide a deeper investigation into the mechanisms through which climate change impacts vegetation phenology. Studies have shown that the NDVI and SIF have similar climatic controls on the SOS. However, compared to NDVI, warming has not significantly affected the EOS derived from SIF in high-latitude regions [42]. Some studies suggested that the NDVI had a longer LOS compared to SIF [43], which can impact the carbon sink of terrestrial ecosystems.

Vegetation fixes the carbon dioxide released into the atmosphere through photosynthesis, establishing a vast terrestrial carbon storage reservoir, which plays a crucial role in maintaining the balance of carbon dioxide in the atmosphere and alleviating the global warming trend [44,45,46]. NEP indicates the ecosystem’s ability to convert carbon sources into sinks, representing the net effect of ecosystem cycling [47,48,49,50]. Exploring the response of NEP to LSP derived from EVI, LAI and SIF will help us better understand the relationships among carbon sequestration, vegetation greenness, vegetation structure information, and vegetation photosynthesis. At present, studies are mainly focusing on exploring the relationship between LSP and NPP, which is directly linked to NEP. For instance, in Central Asia, NPP was positively correlated with delayed SOS and negatively correlated with advanced EOS [10]. In the Qinghai–Tibet Plateau, advanced SOS had a significant negative correlation with spring NPP (R = −0.808, p < 0.01), while delayed EOS had a positive correlation with autumn NPP (R = 0.587, p = 0.058) [51]. Research also revealed that with the prolonged LOS, the NPP of alpine meadows increased by 3.4 gCm⁻²/d, whereas the NPP of alpine grasslands decreased by 0.34 gCm⁻²/d in the Qinghai–Tibet Plateau [52]. Vegetation phenological changes also impact the carbon sequestration capacity of vegetation. Research utilizing the GIMMS LAI3g remote sensing dataset and simulated gridded NEP dataset has revealed that the early greening of vegetation enhanced the NEP of vegetation during summer in the middle and high latitudes of the Northern Hemisphere [53]. However, there are few studies that quantitatively assess the relationship between NEP and the LSP derived from EVI, LAI, and SIF. Further studies can enhance our understanding of the interplay among LSP, climate change, and NEP.

With global warming, the middle and high latitudes of the Northern Hemisphere have undergone significant vegetation phenological changes, which have also impacted the regional vegetation carbon sink capacity. This study utilized three satellite product datasets (EVI, LAI, and SIF) from 2001 to 2021, along with estimated NEP, to assess the responses of LSP and NEP to climate change, as well as the relationship between LSP and NEP. The specific works are as follows: (1) to quantitatively analyze the spatiotemporal changes in LSP and climate forces; (2) to investigate the spatiotemporal changes and climate forces of NEP in regional scale and site scale; (3) to examine the correlation and sensitivity between NEP and LSP; (4) to explore the preseason climate controls and climate sensitivity of LSP; (5) to analyze the climate sensitivity of NEP and its relationship with LSP. This study offers a scientific basis for enhancing the understanding of the relationship between LSP and NEP based on different satellite data under climate change patterns.

2. Materials and Methods

2.1. Study Area

The study area is the north of 30° in the Northern Hemisphere, including Eurasia and North America and parts of northern Africa (Figure 1), with a total land area of approximately 71.36 × 10⁶ km². The vegetation cover type is dominated by Grassland (33.30%), followed by Open Shrublands (17.60%), Woody Savannas (16.60%), Savannas (14.22%), Mixed Forests (9.43%), Evergreen Needleleaf Forests (4.85%), Deciduous Broadleaf Forests (3.34%), Deciduous Needleleaf Forests (0.45%), Evergreen Broadleaf Forests (0.14%) and Closed Shrublands (0.07%). The climate types are complex and diverse in the study area. According to the latest Köppen–Geiger climate classification system [54], there are a total of 26 climate types in this region. In areas with vegetation cover, there are 11 climate types that account for more than 1% of the area, encompassing four major climate categories: arid, temperate, cold, and polar climates (Figure S1). The highest region is located in the Tibetan Plateau, with an elevation of about 7596 m, while the lowest area is at Lake Baikal, with a depth of 1529 m. The area of plains below 200 m accounts for 32.52%. The climate of the study area is characterized by an average temperature of −15.39 °C in January, an average temperature of 16.10 °C in July, and an average annual precipitation of 41.72 mm. The region spans a wide range of latitudes and longitudes, and also exhibits distinct vertical zonation. The diverse and complex geographical environments and ecosystem structures are of significant importance for understanding vegetation phenology changes and carbon storage capacity [55,56].

2.2. Data Sources

2.2.1. Satellite Dataset

The LAI dataset was obtained from the global land surface reflectance satellite product (http://www.glass.umd.edu/LAI/MODIS, accessed on 1 October 2024). It was generated by integrating existing global LAI products such as MODIS C6, GLASS V5, and PROBA-V V1, and utilizing the bidirectional long short-term memory (Bi-LSTM) model. Compared to alternative products, this dataset demonstrated superior performance in resisting cloud and snow contamination, as well as higher accuracy [57]. In our study, we selected the LAI dataset from 2001 to 2021, with a spatial resolution of 0.05°and 500 m, and a temporal resolution of 8-day intervals.

The MODIS Seamless Data Cube (SDC) with a spatial resolution of 500 m and 0.05° is sourced from the Pengcheng Star Cloud platform (https://data-starcloud.pcl.ac.cn/resource/28, accessed on 1 October 2024). Based on the daily MOD09GA 500 m dataset, the MODIS SDC product used time series algorithms to fill in cloud-contaminated observations, eliminated anomalous fluctuation information, and better retained the characteristics of surface dynamic changes, with good spatio-temporal continuity [58]. The MODIS SDC product has been utilized for phenological studies in high latitude and high altitude regions of the Northern Hemisphere, demonstrating a good correlation with phenological records from observation sites [55]. To match the temporal resolution of the LAI dataset, we selected the MODIS SDC reflectance dataset with an 8-day interval for further research from 2001 to 2021.

SIF is an important indicator to characterize photosynthetic rate [59]. Due to its strong correlation with gross primary productivity (GPP), SIF can be used as an approximation for GPP [60]. To maintain consistency with the temporal resolution of LAI and MODIS SDC products, a continuous SIF dataset with an 8-day interval was selected, spanning from 2001 to 2021, accompanied by a spatial resolution of 0.05° (https://data.tpdc.ac.cn/home, accessed on 1 October 2024) [61]. The dataset has been widely used for research on vegetation phenology changes at the regional and hemispheric scales [42,43,62,63,64].

The selected annual NPP tile dataset is the MOD17A3HGF V6.1 product, which has a resolution of 500 m and was obtained from the NASA website (https://lpdaac.usgs.gov/products/mod17a3hgfv061, accessed on 1 October 2024). After a series of post-processing steps on the dataset, including mosaicking, resampling, removing fill values, and clipping, an annual NPP dataset with a spatial resolution of 0.05° was ultimately obtained for the years 2001 to 2021.

2.2.2. Climate and Other Auxiliary Datasets

The monthly daytime and nighttime temperature datasets are sourced from the MODIS MOD11C3 product (https://lpdaac.usgs.gov/, accessed on 1 October 2024), with a resolution of 0.05°. The monthly soil moisture, VPD, maximum temperature, and minimum temperature datasets were obtained from the TerraClimate dataset (https://www.climatologylab.org/, accessed on 1 October 2024), with a resolution of 4.6 km. The time span of the entire climate dataset ranges from 2001 to 2021. To maintain consistency with the spatial resolution of the satellite dataset, the soil moisture, VPD, maximum temperature, and minimum temperature datasets derived from the TerraClimate products were resampled to 0.05°.

The MODIS land cover (MCD12Q1) product, accompanied by the International Geosphere-Biosphere Programme (IGBP) classification system, was adopted to study vegetation phenological changes within specific vegetation functional types. Due to the significant impact of human activities on cultivated land, including planting and irrigation, we have excluded the cultivated land areas. Eventually, ten vegetation types were discussed in our study, including Evergreen Needleleaf Forests (ENF), Evergreen Broadleaf Forests (EBF), Deciduous Needleleaf Forests (DNF), Deciduous Broadleaf Forests (DBF), Mixed Forests (MF), Closed Shrublands (CS), Open Shrublands (OS), Woody Savannas (WS), Savannas (SA) and Grasslands (GRA).

The net ecosystem carbon exchange (NEE) dataset was used in our study from 2001 to 2019. The dataset integrated eddy covariance (EC) flux tower data and remote sensing data and employed a random forest model along with an adaptive genetic optimization algorithm to estimate NEE at 8-day intervals [65]. The estimated NEE site data includes the IGBP land cover classification codes (IGBP_encode), prediction accuracy (pred_r2), and years. We excluded sites with a classification code of cropland, sites with a prediction accuracy below 0.5, and sites with a time span covering fewer than 10 years. Ultimately, 1289 sites were retained for further research. The statistics of sites with different vegetation type functions can be found in Table S1.

The surface elevation data originates from the ETOPO 2022 product released by the NOAA National Centers for Environmental Information, with a resolution of 60 arc-seconds (https://www.ncei.noaa.gov/products/etopo-global-relief-model, accessed on 1 October 2024). The ETOPO 2022 product integrated various topographic data sources, such as airborne lidar, satellite-derived topography, and other terrain survey data, and employed the latest algorithms to correct issues with data stitching, point artifacts, and elevation discrepancies, resulting in higher accuracy compared to the older ETOPO products. To ensure consistency in spatial resolution with other datasets, the ETOPO 2022 product was resampled to 0.05° for further analysis.

2.3. Methods

2.3.1. Retrieval of LSP

In the middle and high latitudes of the Northern Hemisphere, snow cover greatly affects the uncertainty of LSP retrieval. To mitigate the interference of snow on EVI, LAI, and SIF, this study referenced previous research and employed the Normalized Difference Snow Index (NDSI) to identify pixels affected by snow contamination [55]. When NDSI is greater than 0, it indicates the presence of snow cover in the pixel, while NDSI less than 0 indicates the absence of snow contamination. This study first calculated the 8-day interval NDSI and EVI using the MODIS SDC reflectance dataset (Table S2). Then, the NDSI was used to filter snow-affected pixels from the time series of EVI, LAI, and SIF, with background values filled using the 5th percentile of snow-free observations. To eliminate the interference of bare soil, gravel, and sparse vegetation, pixels with annual average EVI < 0.08, annual average LAI < 0.1, or annual average SIF < 0.03 were removed. In order to minimize the impact of different phenological retrieval methods on the phenological results, we evaluated the LSP extracted by three methods and averaged the corresponding phenological parameters for further research.

The first method used cubic spline and Savitzky–Golay (SG) filtering to reconstruct the time series curve. Firstly, the SG filtering was used to smooth the EVI, LAI, and SIF values after removing snow contamination. Subsequently, the daily EVI, LAI, and SIF values were reconstructed using the cubic spline function (Equation (1)) [66]. Finally, the LSP was retrieved by the method of amplitude threshold (Equation (2)) [67]. SOS was defined as the date when the values of EVI, LAI, and SIF ascend to an amplitude of 0.5 during the vegetation development period. EOS was defined as the date when the values of EVI, LAI, and SIF descend to an amplitude of 0.5 during the vegetation dormancy period. POS was defined as the date corresponding to the peak value in the fitted time series.

(1)$S_i(x)=a_i+b_i(x-x_i)+c_i{(x-x_i)}^2+d_i{(x-x_i)}^3$

where $S_i(x)$ represents the fitted values of EVI, LAI, and SIF. $x$ and $x_i$ are two adjacent points (the time step for EVI, LAI and SIF is 8 days, a total of 46 values in a year). $a_i$, $b_i$, $c_i$ and $d_i$ are the fitting coefficients.

(2)$Threshol{d}_{SOS(EOS)}=Valu{e}_{\min ,up(\min ,down)}+C\times (Valu{e}_{\max }-Valu{e}_{\min ,up(\min ,down)})$

where $Valu{e}_{\min ,up}$ represents the minimum value of EVI, LAI, and SIF during the ascending phase of the year. $Valu{e}_{\min ,down}$ represents the minimum value of EVI, LAI, and SIF during the descending phase of the year. $Valu{e}_{\max }$ represents the maximum value of EVI, LAI, and SIF in a year. $C$ is the threshold coefficient, which is 0.5 in this paper [5,27,68,69].

The second method utilized the harmonic analysis of time series (HANTS) algorithm (Equations (3) and (4)) to reconstruct daily time series of EVI, LAI, and SIF, which has been widely employed for reconstructing time series of NDVI, EVI, LAI, and land surface temperature (LST) [70,71,72,73]. Subsequently, the method of relative changing rate (Equation (5)) was applied to retrieve LSP [74]. SOS was defined as the date of the maximum relative changing rate during the vegetation development period, while EOS was defined as the date of the minimum relative changing rate during the vegetation dormancy period. POS was defined as the date of the maximum value of the fitted time series.

(3)${\tilde{y}}^{ij}=a_0+{\displaystyle \sum _{n-1}^N{a}_n\cos (2\pi f_n{t}^{ij})}+b_n\sin (2\pi f_n{t}^{ij})$

(4)$y^{ij}={\tilde{y}}^{ij}+{\epsilon }^{ij}$

where $y^{ij}$, ${\tilde{y}}^{ij}$ and ${\epsilon }^{ij}$ represent the snow-free EVI, LAI, and SIF, the fitted EVI, LAI, and SIF, and the error at pixel position $(i,j)$. $t$ denotes the date corresponding to EVI, LAI, and SIF. $N$ is the maximum number of harmonic components associated with the frequencies ($f_n$). $a_n$ and $b_n$ are the coefficients of cosine and sine components with $f_n$. $a_0$ can be considered as the coefficient at zero frequency.

(5)$f_{ratio(t)}=\left[y_{(\mathrm{t}+1)}-{\mathrm{y}}_{(t)}\right]/{\mathrm{y}}_{(t)}$

where $y_{(t)}$ means the EVI, LAI, and SIF at the date of $t$, and $y_{(t+1)}$ denotes the EVI, LAI, and SIF at the date of $t+1$. $f_{ratio(t)}$ represents the relative changing rate at the date of $t$.

The third method employed the double logistic function (Equation (6)) to reconstruct a daily series of EVI, LAI, and SIF [75]. SOS was defined as the date of the local maximum of the first derivative, while EOS was defined as the date of the local minimum of the first derivative. Similarly, POS was defined as the date of the maximum value of the fitted time series using the double logistic function.

(6)$y_{(t)}=a_1+(a_2-a_1)\times ({\displaystyle \frac{1}{1+e^{-{\mathrm{a}}_3\times (t-a_4)}}}+{\displaystyle \frac{1}{1+e^{a_5\times (t-a_6)}}})$

where $y_{(t)}$ denotes the fitted EVI, LAI, and SIF at the date of $t$. $a_1$ and $a_2$ denote minimum and maximum value of EVI, LAI and SIF for each year. $a_3$ and $a_5$ are the rates of increase and decrease at the inflection point. $a_4$ and $a_6$ are the initial estimates of the inflection points.

2.3.2. NEP Estimation Method in Regional Scale and Site Scale

On a regional scale, NEP is defined as the difference between NPP and soil microbial respiration [76], and it has been widely used to estimate vegetation carbon sources/sinks. When NEP > 0, it indicates that the carbon fixed by the ecosystem is greater than the carbon released into the atmosphere, which is called a carbon sink. When NEP < 0, it signifies that the carbon fixed by the ecosystem is less than the carbon released into the atmosphere, which is called a carbon source [77]. The formula for NEP is as follows:

(7)$NEP=NPP-R_H$

where $R_H$ means soil microbial respiration, which can be estimated through an empirical model [78], which has been used by many scholars [79,80,81,82]. The formula for $R_H$ is as follows:

(8)$R_H=0.22\times \left({\mathrm{e}}^{0.0912T}+{\ln }^{0.3145P+1}\right)\times 30\times 46.5\%$

where $T$ represents the average annual temperature (°C), and $R$ represents the average annual precipitation (mm), and 46.5% is the proportion of carbon emissions released by soil microbial respiration.

In site scale, NEE site data was converted to NEP (NEP = −NEE) for site-scale research [83,84,85].

2.3.3. Trend Analysis Using Ensemble Empirical Mode Decomposition (EEMD) Method

The EEMD method enhances signal analysis by decomposing the original time series into multiple intrinsic mode functions and a secular trend function. This improvement in decomposition capabilities addresses issues such as mode mixing and spurious modes present in the EMD method, making it more effective for analyzing complex time series [86]. The original time series $\chi (t)$ can be expressed as follows:

(9)$\chi (t)={\displaystyle \sum _{i=1}^{n}IM{F}_i(t)+R(t)}$

where $n$ represents the number of intrinsic mode functions (IMFs), $i$ denotes the $i$-th IMF, and $R(t)$ represents the secular trend.

Considering computational efficiency and the robustness of the decomposition, we repeated the above steps 500 times, each time inputting different white noise to obtain varying decomposition results. Finally, we averaged these results to obtain the mean set of IMFs and the secular trend. The cumulative trend (Equation (10)) was defined as the overall change between the initial time ($t_0$) and the ending time ($t_{end}$) [87]. The average changing rate of the instantaneous trend of EEMD (Equation (11)) was defined as the cumulative trend divided by the time interval, representing the average EEMD trend [88]. According to the EEMD method, the types of cumulative trends can be classified into four categories: monotonically increasing, monotonically decreasing, increasing to decreasing, and decreasing to increasing [89,90,91].

(10)$Tren{d}_{accu}=R(t_{end})-R(t_0)$

(11)$Tren{d}_{average}=Tren{d}_{accu}/(t_{end}-t_0)$

where $R(t_0)$ and $R(t_{end})$ represent the initial and final values of the observed time series.

2.3.4. Correlation Analysis

A Pearson correlation analysis (Equation (12)) was conducted to detect the spatial linear correlation between NEP and LSP. When the correlation coefficient is greater than 0, it indicates a positive correlation, while a coefficient less than 0 indicates a negative correlation. The larger the absolute correlation coefficient, the stronger their correlation. When the p-value is less than 0.05, it indicates that the correlation coefficient has passed the significance test, suggesting that there is a significant linear relationship between the two variables.

(12)$r={\displaystyle \frac{{\displaystyle \sum _{i=1}^n(x_i-\overline{x})(y_i-\overline{y})}}{\sqrt{{\displaystyle \sum _{i=1}^n{(x_i-\overline{x})}^2{\displaystyle \sum _{i=1}^n{(y_i-\overline{y})}^2}}}}}$

where $r$ represents the Pearson correlation coefficient. $x_i$ and $y_i$ represent the LSP and NEP for the $i$-th year. $n$ represents the number of years. $\overline{x}$ and $\overline{y}$ represent the average values of multi-year LSP and NEP, respectively.

2.3.5. Partial Correlation Analysis

Partial correlation analysis was used to detect the partial correlations between individual preseason climate factors (soil moisture, nighttime temperature, daytime temperature, and VPD) and LSP, as well as between annual climate factors (minimum temperature, maximum temperature, soil moisture, and VPD) and NEP, while controlling other preseason and annual climate variables. The formula for partial correlation is as follows:

(13)$R_{xy\times z}={\displaystyle \frac{R_{xy}-R_{xz}{R}_{yz}}{\sqrt{(1-R_{xz}^2)(1-R_{yz}^2)}}}$

where $R_{xy}$, $R_{xz}$ and $R_{yz}$ refer to the correlations between $x$ and $y$, $x$ and $z$, and $y$ and $z$, respectively, when controlling for $z$, $y$ and $x$, respectively. Due to the lag and cumulative effects of climate responses on vegetation phenology, the preseason length was set to 0–4 months, and the maximum partial correlation coefficient was used as the correlation value between vegetation phenology and climate factors, with the preseason length defined as the optimal preseason length [74,92,93]. The annual NEP and annual climate factors were directly input into the model without considering their preseason length.

2.3.6. The Analysis of Sensitivity and Three-Dimensional Change Patterns

In this study, the slope of the linear regression model was used to characterize the sensitivity of vegetation phenology and NEP to preseason climate factors and annual climate factors, as well as the sensitivity of NEP to vegetation phenology (Equation (14)). The preseason climate data was derived from pixel values at the optimal preseason length for each climate factor.

(14)$y=S\times x+\epsilon$

where $S$ represents the slope of the linear regression, that is, the sensitivity. For the sensitivity of LSP to preseason climate, $y$ represents the LSP, and $x$ represents the preseason climate factor. For the sensitivity of NEP to annual climate factors, $y$ is NEP, and $x$ is the annual climate factor. For the sensitivity of NEP to LSP, $y$ is NEP, and $x$ is the LSP.

In order to explore the variation patterns of global NEP with altitude, latitude, and longitude, we adopted a research method similar to that of Gao et al. [56], and applied a multiple linear regression model (Equation (15)) combined with a 2.5° × 2.5° sliding window for analysis.

(15)$y=a_1\times Alt+a_2\times Lat+a_3\times Lon+\epsilon$

where $y$ corresponds to the NEP of each pixel. $Alt$, $Lat$ and $Lon$ correspond to the altitude, latitude, and longitude of each pixel, respectively. $a_1$, $a_2$ and $a_3$ are the regression coefficients of multiple regression, characterizing the gradient of NEP with altitude, latitude, and longitude. $\epsilon$ represents random error. According to Gao’s perspective, multiple regression analysis was conducted in the window only when the number of valid values exceeds half and the elevation is above 50 m.

3. Results

3.1. The Spatiotemporal Variation and Climate Driving of LSP

Based on three phenological retrieval methods, the overall spatial distribution of LSP for the three satellite products was similar (Figures S2–S4). However, in terms of specific LSP dates, the LSP retrieved by the HANTS method showed significantly larger differences compared to the LSP retrieved by cubic spline method and double logistic method (Table S3). For instance, the maximum difference in LOS among the various satellite products was 10.26 days and 11.35 days using the cubic spline method and double logistic method, respectively. However, for the HANTS method, the maximum difference in LOS among the different satellite products was 47.06 days. After averaging the three phenological retrieval methods, the maximum difference in LOS across different satellite products was reduced to 6.5 days.

The average LSP derived from three phenological retrieval methods showed that for SOS, EVI_SOS primarily concentrated between 140 and 180 days, while LAI_SOS and SIF_SOS mainly concentrated between 120 and 160 days (Figure 2). For POS, EVI_POS mainly concentrated in the range of 190–230 days, while LAI_POS and SIF_POS primarily focused on the range of 170–210 days. For EOS, EVI_EOS and SIF_EOS mainly concentrated in the range of 240–280 days, while LAI_EOS primarily focused on the range of 220–280 days. For LOS, EVI_LOS, LAI_LOS, and SIF_LOS primarily concentrated in the range of 90–110 days. Based on the cumulative trends of LSP, it was found that the spatial distribution of the cumulative trends was similar across different retrieval methods and satellite products, and the proportion of the four trend types was consistent (Figures S6–S11). The cumulative trends for average SOS and POS were primarily advanced by 0 to 5 days, while the cumulative trend for average EOS was mainly concentrated between −5 and +5 days (Figure 3). The cumulative trend of average LOS showed a predominant extension trend in the middle and high latitudes of the Northern Hemisphere. The trend types for average SOS and POS across different satellite products were predominantly characterized by a monotonic decrease, with proportions generally ranging from 50% to 60% (Figure S5). The trend types for average EOS were predominantly characterized by a monotonic increase, with proportions ranging from 36% to 46%. Across different satellite products, the predominant trend type for average LOS was a monotonic increase, accounting for proportions between 39% and 49%, followed by monotonic decreasing trends with proportions ranging from 24% to 33%.

From the interannual variation analysis (Figure 4), it could be found that the average EEMD trend of SOS, POS, and LOS maintained a consistent direction across the three satellite products. However, the average EEMD trend of EOS showed different trend directions among the three satellite products, with EVI_EOS and SIF_EOS advancing by 0.03 days/yr and 0.01 days/yr, respectively, while LAI_EOS extended by 0.02 days/yr. For different vegetation types, the average SOS and POS (Table S4) showed a consistent trend of advancement across various satellite products in terms of the average EEMD trend. However, there were differences in the direction of the average EEMD trend for EOS and LOS across most vegetation types. For instance, for the LOS of ENF, EVI_LOS and SIF_LOS showed advancements of 0.13 days/yr and 0.03 days/yr, while LAI_LOS showed an extension of 0.15 days/yr. For the LOS of MF, EVI_LOS and SIF_LOS showed advancements of 0.03 days/yr and 0.05 days/yr, while LAI_LOS showed an extension of 0.14 days/yr.

Daytime temperature, nighttime temperature, soil moisture, and VPD were selected as climate factors, and partial correlation analysis was applied to quantify the response of LSP to climate based on different satellite products and phenological retrieval methods (Figure 5 and Figures S12–S21). The response of LSP to climate, derived from different satellite products and various retrieval methods, exhibited a relatively consistent characteristic. Soil moisture showed a higher proportion of positive correlation with SOS and POS than negative correlation, while the proportion of positive correlation with EOS was close to that of negative correlation. The negative correlation between SOS and POS with nighttime temperatures was higher than that with daytime temperatures. The proportion of negative correlation between EOS and daytime temperature exceeded the proportion of negative correlation between EOS and nighttime temperature. The partial correlation of VPD with the phenological parameters (SOS, POS, and EOS) showed a higher proportion of negative correlations than positive correlations, indicating that as the saturation vapor pressure increased, the dates of SOS, POS, and EOS tended to occur earlier in most areas of the middle and high latitudes of the Northern Hemisphere.

3.2. The Spatiotemporal Variation and Climate Driving of NEP

In the spatial distribution of the multi-year average NEP in the middle and high latitudes of the Northern Hemisphere (Figure 6a), regions with NEP greater than 700 gC·m⁻² were primarily found in the eastern part of North America, southern Europe, northeastern China, the Korean Peninsula, and southern Japan. In certain areas of the Appalachian Mountains in North America, the coastal regions of southern Europe, and parts of southern Japan, the NEP could reach 900–1000 gC·m⁻². The NEP below 200 gC·m⁻² was primarily distributed in the vicinity of the Arctic Circle, the arid regions of Central Asia, and the southwestern areas of North America. In the spatial distribution of the NEP cumulative trend (Figure 6b), NEP showed a decreasing trend in some areas, accounting for 24.93%, mainly distributed in the northeastern part of North America. The cumulative trend of NEP in the range of 0–20 gC·m⁻² accounted for the highest proportion, approximately 45%.

For the trend type of NEP (Figure S22), the monotonically increasing trend accounted for the largest proportion, which was 61.04%, followed by the monotonically decreasing trend, which accounted for 17.95%. The trends of increasing to decreasing and decreasing to increasing accounted for a relatively small proportion, which were 12.50% and 8.51%, respectively. Among different vegetation types (Figure S23), EBF had the highest NEP (815.47 gC·m⁻²), followed by DBF (649.10 gC·m⁻²), and OS had the lowest NEP (166.04 gC·m⁻²). The average EEMD trend of NEP among various vegetation types showed that the DBF had the fastest NEP growth rate (1.34 gC·m⁻²a⁻¹), followed by Grasslands (0.74 gC·m⁻²a⁻¹), while DNF exhibited the slowest increase rate of 0.42 gC·m⁻²a⁻¹.

The spatial distribution of partial correlations between NEP and annual average climatic factors (Figure 7) showed that the proportions of positive and negative correlations of soil moisture with NEP (48.08% vs. 51.92%) were essentially consistent across the study area. In the middle latitude region (30°–60°N), soil moisture was predominantly positively correlated with NEP, while soil moisture was mainly negatively correlated with NEP in high-latitude areas (>60°N). In high-latitude regions, NEP was predominantly positively correlated with maximum temperature, while it was mainly negatively correlated with minimum temperature. In high-latitude regions, VPD and NEP primarily exhibited a positive correlation, while in middle-latitude regions, they mainly showed a negative correlation. Across different vegetation types (Figure 8), NEP responded significantly differently to various climate factors, showing a predominant positive correlation with maximum temperature and a negative correlation with minimum temperature, except for EBF and GRA. In particular, the NEP of DNF exhibited a 77.12% positive response to maximum temperature (6.25% significant at p < 0.05 level) and a 76.31% negative response to minimum temperature (11.79% significant). The vegetation types where NEP responded predominantly positively to soil moisture included EBF (59.34%), DNF (69.05%), DBF (72.13%), and GRA (54.52%). The response of NEP to VPD was predominantly positive, except for CS (49.48%) and GRA (44.50%).

Based on site-scale data (1289 sites) in the middle and high latitudes of the Northern Hemisphere, we found that the average NEP was 0.52 gC·m⁻² d⁻¹ during the growing season. The distribution area of the maximum NEP derived from the site was relatively consistent with the distribution area of the maximum NEP obtained from NEP estimation methods based on regional scale (Figure S24a). The average cumulative trend change was 0.02 gC·m⁻² during the vegetation growth period, and the sites with the fastest growth trends were mainly distributed in Europe and North America (Figure S24b). In the NEP trend types based on sites, the proportion of monotonically increasing trends was the highest (46.86%), followed by monotonically decreasing (28.47%), increasing to decreasing (9.31%), and decreasing to increasing (15.36%) (Figure S25). In the site-scale partial correlation of NEP and climatic factors (Figures S26 and S27), the response of NEP to daytime temperature and nighttime temperature was dominated by negative correlation (72.85%) and positive correlation (67.11%), respectively, which was consistent with the response of the NEP estimated at the regional scale to maximum temperature and minimum temperature in the middle latitude region. For the response of NEP to climate across different vegetation types (Figure 9), due to the limited number of sites for certain vegetation types, there were discrepancies between site scale and regional scale. For instance, in ENF (15 sites) and MF (6 sites), the response of NEP to nighttime temperature was predominantly positive at the site scale, whereas the response of NEP to minimum temperature based on the regional scale method was predominantly negative.

3.3. Analysis of the Correlation and Sensitivity Between LSP and NEP

The spatial distribution of the correlation between LSP retrieved from different satellite products and NEP was relatively consistent (Figure 10). Overall, the correlation between SOS and POS and NEP was predominantly negative (SOS: 80.08% for EVI, 76.56% for LAI, 74.35% for SIF; POS: 77.09% for EVI, 74.12% for LAI, 79.09% for SIF). In the study area, the correlation between EOS and NEP showed a higher proportion of negative correlations compared to positive correlations (negative ratio vs. positive ratio: 63.50% vs. 36.50% for EVI; 55.26% vs. 44.74% for LAI; 57.16% vs. 42.84% for SIF). In most regions, the correlation between LOS and NEP showed a positive relationship (73.67% for EVI, 70.80% for LAI, 66.06% for SIF), indicating that the longer LOS, the higher NEP. The sensitivity of NEP to phenology was negative for SOS (−1.38 gC·m⁻²/d for EVI vs. −1.16 gC·m⁻²/d for LAI vs. −0.69 gC·m⁻²/d for SIF), POS (−1.32 gC·m⁻²/d for EVI vs. −1.33 gC·m⁻²/d for LAI vs. −1.52 gC·m⁻²/d for SIF), and EOS (−0.50 gC·m⁻²/d for EVI vs. −0.18 gC·m⁻²/d for LAI vs. −0.19 gC·m⁻²/d for SIF) (Figure 11). The sensitivity of NEP to LOS was positive (0.75 gC·m⁻²/d for EVI vs. 0.63 gC·m⁻²/d for LAI vs. 0.30 gC·m⁻²/d for SIF) (Figure 11).

Based on the NEP sites and the corresponding LSP sites derived from 500 m EVI and LAI retrievals, the correlation between EVI and LAI with NEP showed spatial similarity (Figures S28–S31). The proportion of positive and negative correlations between LSP and NEP sites ranged between 45% and 55% (approximately 10% significant at p < 0.05 level). The sensitivity of NEP sites to phenology was negative for SOS (−0.06 × 10⁻² gC·m⁻²/d for EVI vs. −0.03 × 10⁻² gC·m⁻²/d for LAI), POS (−0.02 × 10⁻² gC·m⁻²/d for EVI and LAI), and EOS (−0.01 × 10⁻² gC·m⁻²/d for EVI and LAI). The sensitivity of NEP sites to LOS was positive (0.04 × 10⁻² gC·m⁻²/d for EVI vs. 0.01 × 10⁻² gC·m⁻²/d for LAI) (Figures S32–S35).

4. Discussion

4.1. The Preseason Climate Controls and Climate Sensitivity of LSP

The preseason climate factors of LSP across different satellite products were analyzed using the EEMD cumulative trend method (Figure 12 and Figures S36 and S37). The spatial distribution of preseason climatic change patterns in LSP from different products was similar. The cumulative trend of preseason daytime temperatures for SOS, POS, and EOS was predominantly characterized by a positive growth trend, with average proportions of 71.93%, 63.80%, and 64.58%, respectively. In the middle and high latitudes of the Northern Hemisphere, the preseason daytime temperatures of SOS, POS, and EOS showed an average annual increase of 0.03 °C, 0.017 °C, and 0.014 °C, respectively, based on EEMD trends (Tables S5 and S6). The cumulative trend of preseason nighttime temperatures showed a higher positive growth proportion compared to daytime temperatures for SOS, POS, and EOS, with average proportions of 78.98%, 76.31%, and 75.58%, respectively. Additionally, the average warming rate of preseason nighttime temperatures exceeded that of daytime temperatures, recorded at 0.033 °C yr⁻¹, 0.02 °C yr⁻¹, and 0.02 °C yr⁻¹ based on the EEMD trend, aligning with global warming patterns [94].

The cumulative trend of preseason soil moisture for SOS, POS, and EOS was predominantly characterized by negative growth, averaging 64.45%, 56.48%, and 55.21%, respectively. The cumulative trend of preseason VPD for SOS, POS, and EOS was predominantly characterized by positive growth, averaging 83.65%, 80.32%, and 75.57%, respectively, indicating a decreasing trend in preseason soil moisture and an increasing trend in air dryness, which have profoundly affected the productivity of terrestrial vegetation ecosystems in the middle and high latitudes of the Northern Hemisphere. In the preseason climate trend types, SOS, POS, and EOS showed a monotonic increasing trend in daytime temperature, nighttime temperature, and VPD with proportions of 45–54%, 62–65%, and 57–68%, respectively, while the proportion of a monotonic decreasing trend in preseason soil moisture was 38–43%. For SOS, POS, and EOS, the proportion of trend types showing a pattern of increasing to decreasing or decreasing to increasing ranged between 7% and 16% (Figures S38–S40).

For different satellite datasets, the average climate sensitivity of LSP was negative for daytime temperature (−1.71 d °C⁻¹ for SOS; −2.31 d °C⁻¹ for POS; −1.17 d °C⁻¹ for EOS), nighttime temperature (−1.88 d °C⁻¹ for SOS; −2.80 d °C⁻¹ for POS; −1.13 d °C⁻¹ for EOS), and VPD (−31.94 d kPa⁻¹ for SOS; −26.72 d kPa⁻¹ for POS; −14.72 d kPa⁻¹ for EOS). Similarly, the average soil moisture sensitivity of LSP was positive (0.38 d mm⁻¹ for SOS; 0.48 d mm⁻¹ for POS; 0.22 d mm⁻¹ for EOS) (Figures S41–S43, and Table S7).

Temperature is the primary factor influencing the spring greening of vegetation in the middle and high latitudes of the Northern Hemisphere [95,96]. With increasing temperatures, the critical threshold for vegetation to green up earlier is met, thereby promoting the early greening process of vegetation. Studies have shown that compared to nighttime warming, the increase in daytime temperature has a greater impact on the slowdown of vegetation greening rates [97]. Asymmetrical warming of daytime and nighttime temperatures may alter the photosynthetic and respiratory rhythms of vegetation, thereby disrupting the balance of terrestrial carbon storage patterns [98]. In high-latitude regions, the negative correlation intensity between SOS_SIF and daytime temperature was lower than that of EVI_SOS and LAI_SOS, which may be due to cold-adapted vegetation having less pigment content and lower photosynthetic activity during the early growth stage, resulting in a weaker response of SIF to increases in daytime temperature compared to EVI and LAI. In high-latitude regions, daytime temperature and nighttime temperature showed positive and negative correlations with the SOS, respectively. This opposing trend may stem from the fact that in high-latitude regions, the increase in daytime temperature promotes vegetation photosynthesis and physiological activities, leading to an earlier SOS. However, despite a warming trend in nighttime temperatures in these regions, they still struggle to meet the temperature thresholds required for vegetation growth, resulting in a delay in the SOS. The increase in daytime and nighttime temperatures promoted the advancement of POS, especially nighttime temperatures (Table S7). The increase in diurnal temperature may have enhanced enzyme activity in vegetation, promoted metabolic processes, and accelerated plant growth and development, leading to an earlier peak in photosynthesis [99]. Compared to daytime temperatures, the increase in nighttime temperatures has a more widespread impact on the extension of EOS. Higher night temperatures prevent the oxidation and destruction of membrane fatty acids and the reduction of chlorophyll due to night frost, thereby delaying the color change of vegetation leaves [100]. Soil moisture is another important factor that limits the growth and development of vegetation. During the greening period of vegetation, excessive preseason soil moisture reduces the heat accumulation necessary for vegetation greening, which adversely affects vegetation. The increase in preseason soil moisture simultaneously weakens the threat of water scarcity that vegetation faces during the growth process, overall extending the peak and dormancy dates of vegetation in the middle and high latitudes of the Northern Hemisphere. When certain vegetation enters the dormancy period, excessive preseason soil moisture, along with low-temperature effects and cloud cover, leads to reduced photosynthesis, accelerating the accumulation of leaf abscission enzymes and leaf aging, thereby causing vegetation to enter dormancy prematurely [68]. In the middle and high latitudes of the Northern Hemisphere, early vegetation greening and extended growing seasons reduce soil moisture, hampering surface temperature regulation by evapotranspiration, thereby exacerbating droughts and extreme weather events, potentially leading to premature leaf aging in certain regions [1,101,102,103]. The increase in preseason VPD is beneficial for vegetation to accumulate the necessary heat for greening, thereby promoting early vegetation greening. However, when vegetation enters the dormant period, the increase in VPD can accelerate the rate of transpiration, increasing the risk of dehydration for the vegetation and potentially leading to an earlier EOS.

4.2. The Gradient Changes in NEP Along Longitude, Latitude, and Elevation, and the Climate Sensitivity of NEP

A global NEP dataset was generated based on regional-scale NEP estimation models from 2001 to 2021. From the spatial distribution of the average NEP over the years, we found that the largest NEP was relatively concentrated, primarily located in tropical rainforest regions, while the smallest NEP was more dispersed, mainly found in alpine tundra and desert grassland areas (Figure 13a). In the high-latitude regions of the Northern Hemisphere, the altitude gradient of NEP tended to be negative, with the gradient changes mainly concentrated between −0.2 and 0 gC·m⁻²/m. In most areas of Africa, South America, and Australia, the altitude gradient of NEP was primarily concentrated between 0 and 0.7 gC·m⁻²/m. In the middle and high latitudes of the Northern Hemisphere and Australia, the latitudinal gradient control of NEP was evident, with the gradient changing from positive to negative as the latitude increased. The latitudinal gradient of NEP exhibited a fragmented pattern in the Northern Hemisphere, mainly concentrated in the range of −50 to 80 gC·m⁻²°⁻¹.

The cumulative changes and trend types for annual average climate variables from 2001 to 2021 were studied (Figure 14 and Figure S44). The spatial distribution of the global average maximum and minimum temperature cumulative change trends was relatively consistent, with the most obvious increase in the northern part of Russia and the Alaska Peninsula, with increases of more than 1.1 °C and 0.98 °C, respectively. The cumulative trend of the average maximum and minimum temperatures was negative, mainly concentrated in the eastern highland region of Africa. The proportion of positive cumulative trend of global soil moisture was 43.02%. The global cumulative trend of VPD with a positive percentage was 82.30%, with the largest increases observed in the central and eastern regions of Africa and South America, exceeding approximately 20.4 kPa. The trends of annual average maximum temperature, minimum temperature, and VPD showed that approximately 70% to 74% of the pixels exhibited a monotonically increasing trend. In contrast, the trend of annual average soil moisture indicated that 44.35% of the pixels had a monotonically decreasing trend.

The sensitivity of global NEP to maximum temperature, minimum temperature and VPD was positive, which were 19.05 gC·m⁻²/°C, 19.51 gC·m⁻²/°C and 302.36 gC·m⁻²/kPa, respectively. However, the sensitivity of NEP to soil moisture was negative, which was −1.19 gC·m⁻²/mm. The spatial distribution of the sensitivity of NEP to the maximum and minimum temperatures was similar (Figure 15). For instance, in high-latitude regions of the Northern Hemisphere and high-altitude areas like the Tibetan Plateau, the sensitivity of NEP to maximum and minimum temperatures was positive. In contrast, in southern Africa, parts of South America, and Australia, the sensitivity of NEP to maximum and minimum temperatures was negative. The sensitivity of NEP to soil moisture and VPD was negative and positive in high-latitude regions of the Northern Hemisphere and positive and negative in some arid regions, such as Central Asia, western North America, western Australia, and southern Africa.

In high-latitude and high-altitude regions of the Northern Hemisphere, temperature is the dominant factor controlling the productivity of cold-adapted vegetation. The increase in temperature is beneficial for the growth and development of vegetation. However, the melting of permafrost and the increase in soil moisture can consume heat and lower soil temperatures. Additionally, high soil humidity is not conducive to the growth of plant roots, which ultimately affects the productivity of vegetation [104]. In arid regions, due to the inherently high temperatures and significant evaporation, soil moisture becomes the most critical factor limiting vegetation growth [47,105,106]. In arid areas, the increase in VPD may accelerate the rate of transpiration of vegetation, which may lead to a more serious risk of dehydration of vegetation. To mitigate water loss, plants may close stomata to reduce water evaporation, which will affect the photosynthesis rate of vegetation [107]. Due to the closure of plant stomata, the absorption of carbon dioxide by plants decreases, inhibiting photosynthesis and thereby affecting the productivity of the entire ecosystem [108]. In high-latitude regions, elevated VPD, typically associated with warmer climates, may extend the plant growing season, offering more temperature and light resources that vegetation could more effectively utilize, potentially boosting photosynthesis rates, fostering plant growth and development, and increasing vegetation productivity [109].

4.3. Relationship Between LSP Derived from Different Satellite Products and NEP

Remote sensing phenology can effectively reflect the seasonal cycles of photosynthetic carbon assimilation in vegetation other than tropical plants, making it a powerful tool for studying the carbon cycle in the context of climate change [42]. Many vegetation indices, such as NDVI and EVI, are commonly used to characterize vegetation coverage and growth conditions. The LAI can measure the ability of the vegetation canopy to intercept light energy. A higher LAI indicates that more light energy can be captured. The SIF directly reflects the photosynthetic activity of vegetation. In phenological studies, these indicators are widely used as alternative metrics for vegetation photosynthesis. This article utilizes various remote sensing products with different characteristics to investigate LSP, aiming to deepen the understanding of how LSP promotes NEP dynamics.

Studies have shown that the relationship between SOS and NEP was predominantly negative, implying that an earlier SOS extends the duration of spring vegetation photosynthesis, thereby promoting NEP accumulation. The strongest negative correlation was found between SOS_EVI and NEP, with NEP showing the highest sensitivity to SOS_EVI (Figure 11 and Figure 12). In contrast, the smallest proportion of negative correlation was observed between SOS_SIF and NEP, with NEP exhibiting the lowest sensitivity to SOS_SIF. This may be because, during the early stages of vegetation growth, EVI better reflects the impact of vegetation growth status on photosynthesis, while SIF focuses more on photosynthetic efficiency rather than SOS. The earlier onset of POS indicates that the vegetation canopy achieves maximum efficiency in intercepting and utilizing light energy, which explains the negative correlation between the increase in NEP and the advancement of POS. The negative correlation between POS_SIF and NEP is the strongest, and the sensitivity of NEP to POS_SIF is also the highest, which further supports this viewpoint. In the later stages of vegetation growth, as leaves age and exhibit yellowing or shedding, heterotrophic respiration in the soil increases. Simultaneously, environmental conditions such as light intensity, temperature, and soil moisture deteriorate, leading to a decline in vegetation photosynthetic efficiency. These factors are all detrimental to the accumulation of NEP. Due to the influence of the pigments in vegetation leaves and their growth conditions, the sensitivity of NEP to EVI_EOS was enhanced again in the later stages of the growing season. Extended LOS is conducive to the accumulation of NEP. However, some argue that extending the LOS and EOS may not increase vegetation productivity accumulation for short-lived vegetation due to decreased photosynthesis and respiration exceeding photosynthesis in the later stages of growth [110].

In this study, the MODIS NPP dataset and the empirical soil microbial respiration model were used to estimate the NEP at the middle and high latitudes of the Northern Hemisphere and at the global scale, respectively. However, factors such as cloud cover, aerosols, atmospheric environmental changes, and sensor performance pose challenges to the quality of MODIS NPP products [111]. As the soil microbial empirical model used was widely applied in northern China [79,80,81], its regional adaptability needed to be assessed before being extended to hemispheric and global scales. Therefore, future research needs to employ advanced algorithms to further improve the NPP model, enhance the collection and sharing of soil microbial data on a global scale, develop more universally applicable soil microbial models, and conduct multi-scale and multi-factor model validation and assessment. Understanding the relationship between LSP and NEP can indirectly forecast the impact of climate change on vegetation carbon sequestration. For instance, sustained early spring warming may induce vegetation in the middle and high latitudes of the Northern Hemisphere to bud and reach peak growth earlier, enabling plants to more effectively utilize photosynthetic energy, thus storing more carbon, provided that water and other resources are not limited. However, in the later stages of vegetation growth, climate change may affect the timing of leaf drop, which in turn could lead to a reduction in carbon sequestration. Furthermore, incorporating different phenological parameters into global carbon uptake estimation models can more accurately assess carbon absorption and release at different stages globally. By utilizing knowledge from LSP, we can selectively plant vegetation types with longer growing periods, significantly enhancing their carbon sink function and effectively mitigating the challenges posed by climate change. This research has a profound impact on promoting sustainable development strategies and building a green, low-carbon future.

5. Conclusions

In this study, three different satellite products (EVI, LAI, and SIF) were applied from 2001 to 2021 to explore vegetation phenological changes and climate drivers in the middle and high latitudes of the Northern Hemisphere based on various phenological retrieval methods. Furthermore, the spatial distribution, changes, and climate response of NEP at both regional and site scales were quantified. Additionally, the sensitivity and correlation of NEP with LSP were explored. In the LSP retrieved from different remote sensing products, the cumulative trends of SOS and POS showed a primary advancement of around 5 days from 2001 to 2021. The cumulative trend of average LOS showed a predominant extension trend. The average LSP exhibited a predominantly monotonic trend, with increasing to decreasing and decreasing to increasing trends accounting for less than 15%. The response of LSP to climate under different satellite products and inversion methods showed a consistent overall trend. The proportion of positive correlations between soil moisture, SOS and POS was greater than that of negative correlations. Daytime and nighttime temperatures were generally negatively correlated with SOS and POS, with nighttime temperatures showing a higher negative correlation with SOS and POS compared to daytime temperatures. The negative correlation between VPD and vegetation phenological parameters (SOS, POS, EOS) was greater than the positive correlation.

From 2001 to 2021, the trend types of NEP showed that the largest proportion was the monotonically increasing trend (61.04%), followed by the monotonically decreasing (17.95%), the increasing to decreasing trend (12.50%), and the decreasing to increasing trend (8.51%) in the middle and high latitudes of the Northern Hemisphere. The proportion of positive and negative correlations between soil moisture and NEP (48.08% vs. 51.92%) was roughly consistent in the study area. In high-latitude regions, NEP was primarily positively correlated with maximum temperature and negatively correlated with minimum temperature. For SOS, the proportion of negative correlation between SOS_EVI and NEP was the largest, and NEP exhibited the strongest sensitivity to SOS_EVI. As vegetation growth reached its peak, the vegetation canopy gradually attained its maximum, resulting in higher light energy interception and utilization efficiency. The negative correlation between POS_SIF and NEP was the strongest, and NEP showed the highest sensitivity to POS_SIF. As leaf color and vegetation growth conditions change, the sensitivity of NEP to EVI_EOS increases once again. The study emphasized that, while prioritizing biodiversity conservation, we could leverage remote sensing phenology to promote the cultivation of vegetation with longer growing seasons, thereby enhancing the carbon sequestration capacity of ecosystems. The study provided a scientific basis for climate change mitigation and carbon management.

Author Contributions

Conceptualization, X.M. and H.Y. (Hanmin Yin); methodology, H.Y. (Hanmin Yin) and Q.L.; software, H.Y. (Hanmin Yin), Q.L., W.Y., Y.L., J.W. and J.Y.; resources, Q.L.; data curation, H.Y. (Hanmin Yin); writing—original draft preparation, H.Y. (Hanmin Yin); writing—review and editing, H.Y. (Hanmin Yin) and Q.L.; supervision, X.L. and H.Y. (Huping Ye); project administration, Q.L.; funding acquisition, Q.L. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

The MODIS Seamless Data Cube (SDC) is obtained from the Pengcheng Star Cloud platform (http://data.starcloud.pcl.ac.cn/resource/28, accessed on 1 October 2024).

Acknowledgments

We would like to express our gratitude to the Land Processes Distributed Active Archive Center (LPDAAC) of NASA for providing the MOD11C3, MCD12C1, and MOD17A3HGF V6.1 datasets, as well as to TerraClimate for the monthly maximum temperature, minimum temperature, soil moisture, and vapor pressure deficit data. We also thank Zhang Wenqiang for sharing the global net ecosystem carbon exchange dataset and Zhang Yao for sharing the global continuous satellite-induced fluorescence (CSIF) dataset. Additionally, we appreciate the support from the GLASS LAI product of MODIS surface reflectance data. Their contributions have been crucial to the success of this project, and we sincerely thank them for their support.

Conflicts of Interest

The authors declare no conflicts of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures

Figure 1. The spatial distribution of vegetation types (a) and elevation (b) in the middle and high latitudes of the Northern Hemisphere. The left histogram (a) represents the proportion of the corresponding land cover type, and the right histogram (b) represents the proportion of the corresponding elevation. ENF: Evergreen Needleleaf Forests; EBF: Evergreen Broadleaf Forests; DNF: Deciduous Needleleaf Forests; DBF: Deciduous Broadleaf Forests; MF: Mixed Forests; CS: Closed Shrublands; OS: Open Shrublands; WS: Woody Savannas; SA: Savannas; GRA: Grasslands.

Figure 2. The spatial distribution of multi-year average LSP retrieved by three different methods and different satellite products over the period of 2001–2021. Each histogram in the bottom left corner represents the frequency distribution of the corresponding phenological date.

Figure 3. The spatial distribution of the cumulative trend of average LSP using the EEMD method based on different satellite products and retrieval methods from 2001 to 2021. The histogram inserted in the bottom left corner of each panel describes the frequency distribution of the corresponding cumulative trend.

Figure 4. Interannual variations of average LSP from 2001 to 2021. Panels (a–d) represent the interannual variations of average SOS, POS, EOS, and LOS based on different remote sensing satellite datasets and retrieval methods. The “linear trend” refers to the trend of phenological parameters using a linear regression model. The “EEMD trend” signifies the average rate of instantaneous trend change, representing the average EEMD trend.

Figure 5. Spatial distribution of partial correlation coefficients between average phenological parameters (SOS and EOS) and climate factors. The histograms represent the frequency distribution of the corresponding partial correlation coefficients. The abbreviations SM, DT, NT, and VPD represent soil moisture, daytime temperature, nighttime temperature, and vapor pressure deficit, respectively. The black dots represent pixels that have passed the significance test (p [less than] 0.05). The spatial distribution of the partial correlation coefficients between average POS and climatic factors can be found in Figure S12.

Figure 6. Spatial distribution of multi-year average NEP (a) and cumulative trend changes (b) based on the NEP estimation model in regional scales from 2001 to 2021.

Figure 7. Spatial distribution of partial correlation coefficients between NEP and annual average climate factors. The histograms represent the frequency distribution of the corresponding partial correlation coefficients. The abbreviations SM, TEMMAX, TEMMIN, and VPD represent soil moisture, maximum temperature, minimum temperature, and vapor pressure deficit, respectively. The black dots represent pixels that have passed the significance test (p [less than] 0.05).

Figure 8. Frequency of partial correlation between NEP estimated at the regional scale and climate driving factors for different vegetation functional types. The upward and downward bars represent the percentages of positive and negative correlations, respectively. The white bar indicates the proportion of significant correlation (p [less than] 0.05). The definitions of the abbreviations on the x-axis are the same as in Figure 1.

Figure 9. Frequency of partial correlation between NEP sites and climate driving factors for different vegetation functional types. The upward and downward bars represent the percentages of positive and negative correlations, respectively. The white bar indicates the proportion of significant correlation (p [less than] 0.05).

Figure 10. The spatial distribution of the correlation between LSP and NEP based on different satellite products. The black dots represent pixels that have passed the significance test (p [less than] 0.05).

Figure 11. The spatial distribution of the sensitivity of NEP to LSP based on different satellite products. The black dots represent pixels that have passed the significance test (p [less than] 0.05).

Figure 12. The spatial distribution of the cumulative trend changes in preseason climate factors based on EVI data sources.

Figure 13. The spatial distribution of multi-year average NEP (a) and its three-dimensional variations (b–d), along with their average latitudinal distribution (e–h).

Figure 14. The spatial distribution (a–d) and latitude distribution (e–h) of the cumulative trends of climate variables from 2001 to 2021.

Figure 15. The spatial distribution (a–d) and latitude distribution (e–h) of the climate sensitivity of NEP.

Supplementary Materials

The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/rs16214101/s1.

References

1. Piao, S.; Liu, Q.; Chen, A.; Janssens, I.A.; Fu, Y.; Dai, J.; Liu, L.; Lian, X.; Shen, M.; Zhu, X. Plant phenology and global climate change: Current progresses and challenges. Glob. Chang. Biol.; 2019; 25, pp. 1922-1940. [DOI: https://dx.doi.org/10.1111/gcb.14619] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/30884039]

2. Moon, M.; Zhang, X.; Henebry, G.M.; Liu, L.; Gray, J.M.; Melaas, E.K.; Friedl, M.A. Long-term continuity in land surface phenology measurements: A comparative assessment of the MODIS land cover dynamics and VIIRS land surface phenology products. Remote Sens. Environ.; 2019; 226, pp. 74-92. [DOI: https://dx.doi.org/10.1016/j.rse.2019.03.034]

3. Peters-Lidard, C.D.; Zion, M.S.; Wood, E.F. A soil-vegetation-atmosphere transfer scheme for modeling spatially variable water and energy balance processes. J. Geophys. Res. Atmos.; 1997; 102, pp. 4303-4324. [DOI: https://dx.doi.org/10.1029/96JD02948]

4. Li, B.; Wang, R.; Chen, J.M. Responses of phenology to preseason drought and soil temperature for different land cover types on the Mongolian Plateau. Sci. Total Environ.; 2024; 926, 171895. [DOI: https://dx.doi.org/10.1016/j.scitotenv.2024.171895] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/38531448]

5. Dong, L.; Wu, C.; Wang, X.; Zhao, N. Satellite observed delaying effects of increased winds on spring green-up dates. Remote Sens. Environ.; 2023; 284, 113363. [DOI: https://dx.doi.org/10.1016/j.rse.2022.113363]

6. Li, N.; Zhan, P.; Pan, Y.; Zhu, X.; Li, M.; Zhang, D. Comparison of Remote Sensing Time-Series Smoothing Methods for Grassland Spring Phenology Extraction on the Qinghai–Tibetan Plateau. Remote Sens.; 2020; 12, 3383. [DOI: https://dx.doi.org/10.3390/rs12203383]

7. Wang, X.; Xiao, J.; Li, X.; Cheng, G.; Ma, M.; Zhu, G.; Altaf Arain, M.; Andrew Black, T.; Jassal, R.S. No trends in spring and autumn phenology during the global warming hiatus. Nat. Commun.; 2019; 10, 2389. [DOI: https://dx.doi.org/10.1038/s41467-019-10235-8]

8. Linderholm, H.W. Growing season changes in the last century. Agric. For. Meteorol.; 2006; 137, pp. 1-14. [DOI: https://dx.doi.org/10.1016/j.agrformet.2006.03.006]

9. Gao, X.; Zhao, D. Impacts of climate change on vegetation phenology over the Great Lakes Region of Central Asia from 1982 to 2014. Sci. Total Environ.; 2022; 845, 157227. [DOI: https://dx.doi.org/10.1016/j.scitotenv.2022.157227]

10. Wu, L.; Ma, X.; Dou, X.; Zhu, J.; Zhao, C. Impacts of climate change on vegetation phenology and net primary productivity in arid Central Asia. Sci. Total Environ.; 2021; 796, 149055. [DOI: https://dx.doi.org/10.1016/j.scitotenv.2021.149055]

11. Ren, P.; Liu, Z.; Zhou, X.; Peng, C.; Xiao, J.; Wang, S.; Li, X.; Li, P. Strong controls of daily minimum temperature on the autumn photosynthetic phenology of subtropical vegetation in China. For. Ecosyst.; 2021; 8, 31. [DOI: https://dx.doi.org/10.1186/s40663-021-00309-9] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/34721934]

12. Qi, Y.; Wang, H.; Ma, X.; Zhang, J.; Yang, R. Relationship between vegetation phenology and snow cover changes during 2001–2018 in the Qilian Mountains. Ecol. Indic.; 2021; 133, 108351. [DOI: https://dx.doi.org/10.1016/j.ecolind.2021.108351]

13. Caparros-Santiago, J.A.; Rodriguez-Galiano, V.; Dash, J. Land surface phenology as indicator of global terrestrial ecosystem dynamics: A systematic review. ISPRS J. Photogramm. Remote Sens.; 2021; 171, pp. 330-347. [DOI: https://dx.doi.org/10.1016/j.isprsjprs.2020.11.019]

14. Zhang, J.; Zhao, J.; Wang, Y.; Zhang, H.; Zhang, Z.; Guo, X. Comparison of land surface phenology in the Northern Hemisphere based on AVHRR GIMMS3g and MODIS datasets. ISPRS J. Photogramm. Remote Sens.; 2020; 169, pp. 1-16. [DOI: https://dx.doi.org/10.1016/j.isprsjprs.2020.08.020]

15. Chen, X.; Yang, Y. Observed earlier start of the growing season from middle to high latitudes across the Northern Hemisphere snow-covered landmass for the period 2001–2014. Environ. Res. Lett.; 2020; 15, 034042. [DOI: https://dx.doi.org/10.1088/1748-9326/ab6d39]

16. Garonna, I.; de Jong, R.; Schaepman, M.E. Variability and evolution of global land surface phenology over the past three decades (1982–2012). Glob. Chang. Biol.; 2016; 22, pp. 1456-1468. [DOI: https://dx.doi.org/10.1111/gcb.13168]

17. Keenan, T.F.; Gray, J.; Friedl, M.A.; Toomey, M.; Bohrer, G.; Hollinger, D.Y.; Munger, J.W.; O’Keefe, J.; Schmid, H.P.; Wing, I.S. et al. Net carbon uptake has increased through warming-induced changes in temperate forest phenology. Nat. Clim. Chang.; 2014; 4, pp. 598-604. [DOI: https://dx.doi.org/10.1038/nclimate2253]

18. Barichivich, J.; Briffa, K.; Osborn, T.; Melvin, T.; Caesar, J. Thermal growing season and timing of biospheric carbon uptake across the Northern Hemisphere. Glob. Biogeochem. Cycles; 2012; 26, [DOI: https://dx.doi.org/10.1029/2012GB004312]

19. Li, P.; Peng, C.; Wang, M.; Luo, Y.; Li, M.; Zhang, K.; Zhang, D.; Zhu, Q. Dynamics of vegetation autumn phenology and its response to multiple environmental factors from 1982 to 2012 on Qinghai-Tibetan Plateau in China. Sci. Total Environ.; 2018; 637–638, pp. 855-864. [DOI: https://dx.doi.org/10.1016/j.scitotenv.2018.05.031]

20. Richardson, A.D.; Keenan, T.F.; Migliavacca, M.; Ryu, Y.; Sonnentag, O.; Toomey, M. Climate change, phenology, and phenological control of vegetation feedbacks to the climate system. Agric. For. Meteorol.; 2013; 169, pp. 156-173. [DOI: https://dx.doi.org/10.1016/j.agrformet.2012.09.012]

21. Peñuelas, J.; Rutishauser, T.; Filella, I. Phenology feedbacks on climate change. Science; 2009; 324, pp. 887-888. [DOI: https://dx.doi.org/10.1126/science.1173004] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/19443770]

22. Huntington, T.G. Climate change, growing season length, and transpiration: Plant response could alter hydrologic regime. Plant Biol.; 2004; 6, pp. 651-653. [DOI: https://dx.doi.org/10.1055/s-2004-830353] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/15570468]

23. Raftery, A.E.; Zimmer, A.; Frierson, D.M.W.; Startz, R.; Liu, P. Less than 2 °C warming by 2100 unlikely. Nat. Clim. Chang.; 2017; 7, pp. 637-641. [DOI: https://dx.doi.org/10.1038/nclimate3352] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/30079118]

24. Templ, B.; Koch, E.; Bolmgren, K.; Ungersböck, M.; Paul, A.; Scheifinger, H.; Rutishauser, T.; Busto, M.; Chmielewski, F.-M.; Hájková, L. Pan European Phenological database (PEP725): A single point of access for European data. Int. J. Biometeorol.; 2018; 62, pp. 1109-1113. [DOI: https://dx.doi.org/10.1007/s00484-018-1512-8] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/29455297]

25. Aono, Y.; Kazui, K. Phenological data series of cherry tree flowering in Kyoto, Japan, and its application to reconstruction of springtime temperatures since the 9th century. Int. J. Climatol.; 2008; 28, pp. 905-914. [DOI: https://dx.doi.org/10.1002/joc.1594]

26. Li, G.; Wu, C.; Chen, Y.; Huang, C.; Zhao, Y.; Wang, Y.; Ma, M.; Ding, Z.; Yu, P.; Tang, X. Increasing temperature regulates the advance of peak photosynthesis timing in the boreal ecosystem. Sci. Total Environ.; 2023; 882, 163587. [DOI: https://dx.doi.org/10.1016/j.scitotenv.2023.163587]

27. Guo, M.; Wu, C.; Peng, J.; Lu, L.; Li, S. Identifying contributions of climatic and atmospheric changes to autumn phenology over mid-high latitudes of Northern Hemisphere. Glob. Planet. Chang.; 2021; 197, 103396. [DOI: https://dx.doi.org/10.1016/j.gloplacha.2020.103396]

28. Zhou, X.; Geng, X.; Yin, G.; Hänninen, H.; Hao, F.; Zhang, X.; Fu, Y.H. Legacy effect of spring phenology on vegetation growth in temperate China. Agric. For. Meteorol.; 2020; 281, 107845. [DOI: https://dx.doi.org/10.1016/j.agrformet.2019.107845]

29. Zhang, X.; Friedl, M.A.; Schaaf, C.B. Global vegetation phenology from Moderate Resolution Imaging Spectroradiometer (MODIS): Evaluation of global patterns and comparison with in situ measurements. J. Geophys. Res. Biogeosci.; 2006; 111, [DOI: https://dx.doi.org/10.1029/2006JG000217]

30. Helman, D. Land surface phenology: What do we really ‘see’from space?. Sci. Total Environ.; 2018; 618, pp. 665-673. [DOI: https://dx.doi.org/10.1016/j.scitotenv.2017.07.237]

31. Yan, E.; Wang, G.; Lin, H.; Xia, C.; Sun, H. Phenology-based classification of vegetation cover types in Northeast China using MODIS NDVI and EVI time series. Int. J. Remote Sens.; 2015; 36, pp. 489-512. [DOI: https://dx.doi.org/10.1080/01431161.2014.999167]

32. Atkinson, P.M.; Dash, J.; Jeganathan, C. Amazon vegetation greenness as measured by satellite sensors over the last decade. Geophys. Res. Lett.; 2011; 38. [DOI: https://dx.doi.org/10.1029/2011GL049118]

33. Myneni, R.B.; Hall, F.G.; Sellers, P.J.; Marshak, A.L. The interpretation of spectral vegetation indexes. IEEE Trans. Geosci. Remote Sens.; 1995; 33, pp. 481-486. [DOI: https://dx.doi.org/10.1109/TGRS.1995.8746029]

34. Ge, C.; Sun, S.; Yao, R.; Sun, P.; Li, M.; Bian, Y. Long-term vegetation phenology changes and response to multi-scale meteorological drought on the Loess Plateau, China. J. Hydrol.; 2022; 614, 128605. [DOI: https://dx.doi.org/10.1016/j.jhydrol.2022.128605]

35. Zhao, J.; Zhang, H.; Zhang, Z.; Guo, X.; Li, X.; Chen, C. Spatial and Temporal Changes in Vegetation Phenology at Middle and High Latitudes of the Northern Hemisphere over the Past Three Decades. Remote Sens.; 2015; 7, pp. 10973-10995. [DOI: https://dx.doi.org/10.3390/rs70810973]

36. Wu, D.; Zhao, X.; Liang, S.; Zhou, T.; Huang, K.; Tang, B.; Zhao, W. Time-lag effects of global vegetation responses to climate change. Glob. Chang. Biol.; 2015; 21, pp. 3520-3531. [DOI: https://dx.doi.org/10.1111/gcb.12945]

37. Glenn, E.P.; Huete, A.R.; Nagler, P.L.; Nelson, S.G. Relationship between remotely-sensed vegetation indices, canopy attributes and plant physiological processes: What vegetation indices can and cannot tell us about the landscape. Sensors; 2008; 8, pp. 2136-2160. [DOI: https://dx.doi.org/10.3390/s8042136]

38. Chen, A.; Mao, J.; Ricciuto, D.; Xiao, J.; Frankenberg, C.; Li, X.; Thornton, P.E.; Gu, L.; Knapp, A.K. Moisture availability mediates the relationship between terrestrial gross primary production and solar-induced chlorophyll fluorescence: Insights from global-scale variations. Glob. Chang. Biol.; 2021; 27, pp. 1144-1156. [DOI: https://dx.doi.org/10.1111/gcb.15373]

39. Wang, C.; Li, J.; Liu, Q.; Zhong, B.; Wu, S.; Xia, C. Analysis of differences in phenology extracted from the enhanced vegetation index and the leaf area index. Sensors; 2017; 17, 1982. [DOI: https://dx.doi.org/10.3390/s17091982]

40. Fang, J.; Li, X.; Xiao, J.; Yan, X.; Li, B.; Liu, F. Vegetation photosynthetic phenology dataset in northern terrestrial ecosystems. Sci. Data; 2023; 10, 300. [DOI: https://dx.doi.org/10.1038/s41597-023-02224-w]

41. Verger, A.; Filella, I.; Baret, F.; Peñuelas, J. Vegetation baseline phenology from kilometric global LAI satellite products. Remote Sens. Environ.; 2016; 178, pp. 1-14. [DOI: https://dx.doi.org/10.1016/j.rse.2016.02.057]

42. Chen, A.; Meng, F.; Mao, J.; Ricciuto, D.; Knapp, A.K. Photosynthesis phenology, as defined by solar-induced chlorophyll fluorescence, is overestimated by vegetation indices in the extratropical Northern Hemisphere. Agric. For. Meteorol.; 2022; 323, 109027. [DOI: https://dx.doi.org/10.1016/j.agrformet.2022.109027]

43. Meng, F.; Huang, L.; Chen, A.; Zhang, Y.; Piao, S. Spring and autumn phenology across the Tibetan Plateau inferred from normalized difference vegetation index and solar-induced chlorophyll fluorescence. Big Earth Data; 2021; 5, pp. 182-200. [DOI: https://dx.doi.org/10.1080/20964471.2021.1920661]

44. Mu, Y.; Jia, X.; Ye, Z.; Zha, T.; Guo, X.; Black, T.A.; Zhang, Y.; Hao, S.; Han, C.; Gao, S. Dry-season length affects the annual ecosystem carbon balance of a temperate semi-arid shrubland. Sci. Total Environ.; 2024; 917, 170532. [DOI: https://dx.doi.org/10.1016/j.scitotenv.2024.170532] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/38296104]

45. Qiu, S.; Liang, L.; Wang, Q.; Geng, D.; Wu, J.; Wang, S.; Chen, B. Estimation of European terrestrial ecosystem NEP based on an improved CASA model. IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens.; 2022; 16, pp. 1244-1255. [DOI: https://dx.doi.org/10.1109/JSTARS.2022.3233128]

46. Poulter, B.; Frank, D.; Ciais, P.; Myneni, R.B.; Andela, N.; Bi, J.; Broquet, G.; Canadell, J.G.; Chevallier, F.; Liu, Y.Y. Contribution of semi-arid ecosystems to interannual variability of the global carbon cycle. Nature; 2014; 509, pp. 600-603. [DOI: https://dx.doi.org/10.1038/nature13376]

47. Hou, Q.; Yang, H.; Wu, J.; Yu, X. Carbon budget response to climate change varies with grassland type in Qilian Mountains, China. Glob. Ecol. Conserv.; 2023; 47, e02670. [DOI: https://dx.doi.org/10.1016/j.gecco.2023.e02670]

48. Cui, E.; Bian, C.; Luo, Y.; Niu, S.; Wang, Y.; Xia, J. Spatial variations in terrestrial net ecosystem productivity and its local indicators. Biogeosciences; 2020; 17, pp. 6237-6246. [DOI: https://dx.doi.org/10.5194/bg-17-6237-2020]

49. Fernández-Martínez, M.; Vicca, S.; Janssens, I.A.; Sardans, J.; Luyssaert, S.; Campioli, M.; Chapin Iii, F.S.; Ciais, P.; Malhi, Y.; Obersteiner, M. et al. Nutrient availability as the key regulator of global forest carbon balance. Nat. Clim. Chang.; 2014; 4, pp. 471-476. [DOI: https://dx.doi.org/10.1038/nclimate2177]

50. Hassan, Q.K.; Bourque, C.P.A.; Meng, F.-R. Estimation of daytime net ecosystem CO2 exchange over balsam fir forests in eastern Canada: Combining averaged tower-based flux measurements with remotely sensed MODIS data. Can. J. Remote Sens.; 2006; 32, pp. 405-416. [DOI: https://dx.doi.org/10.5589/m07-009]

51. Wang, S.; Zhang, B.; Yang, Q.; Chen, G.; Yang, B.; Lu, L.; Shen, M.; Peng, Y. Responses of net primary productivity to phenological dynamics in the Tibetan Plateau, China. Agric. For. Meteorol.; 2017; 232, pp. 235-246. [DOI: https://dx.doi.org/10.1016/j.agrformet.2016.08.020]

52. Sun, H.; Chen, Y.; Xiong, J.; Ye, C.; Yong, Z.; Wang, Y.; He, D.; Xu, S. Relationships between climate change, phenology, edaphic factors, and net primary productivity across the Tibetan Plateau. Int. J. Appl. Earth Obs. Geoinf.; 2022; 107, 102708. [DOI: https://dx.doi.org/10.1016/j.jag.2022.102708]

53. Ren, Y.; Qiu, J.; Zeng, Z.; Liu, X.; Sitch, S.; Pilegaard, K.; Yang, T.; Wang, S.; Yuan, W.; Jain, A.K. Earlier spring greening in Northern Hemisphere terrestrial biomes enhanced net ecosystem productivity in summer. Commun. Earth Environ.; 2024; 5, 122. [DOI: https://dx.doi.org/10.1038/s43247-024-01270-5]

54. Beck, H.E.; Zimmermann, N.E.; McVicar, T.R.; Vergopolan, N.; Berg, A.; Wood, E.F. Present and future Köppen-Geiger climate classification maps at 1-km resolution. Sci. Data; 2018; 5, 180214. [DOI: https://dx.doi.org/10.1038/sdata.2018.214] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/30375988]

55. Yin, H.; Liu, Q.; Liao, X.; Ye, H.; Li, Y.; Ma, X. Refined Analysis of Vegetation Phenology Changes and Driving Forces in High Latitude Altitude Regions of the Northern Hemisphere: Insights from High Temporal Resolution MODIS Products. Remote Sens.; 2024; 16, 1744. [DOI: https://dx.doi.org/10.3390/rs16101744]

56. Gao, M.; Wang, X.; Meng, F.; Liu, Q.; Li, X.; Zhang, Y.; Piao, S. Three-dimensional change in temperature sensitivity of northern vegetation phenology. Glob. Chang. Biol.; 2020; 26, pp. 5189-5201. [DOI: https://dx.doi.org/10.1111/gcb.15200]

57. Ma, H.; Liang, S. Development of the GLASS 250-m leaf area index product (version 6) from MODIS data using the bidirectional LSTM deep learning model. Remote Sens. Environ.; 2022; 273, 112985. [DOI: https://dx.doi.org/10.1016/j.rse.2022.112985]

58. Liang, X.; Liu, Q.; Wang, J.; Chen, S.; Gong, P. Global 500 m seamless dataset (2000–2022) of land surface reflectance generated from MODIS products. Earth Syst. Sci. Data Discuss.; 2024; 16, pp. 177-200. [DOI: https://dx.doi.org/10.5194/essd-16-177-2024]

59. Baker, N.R. Chlorophyll fluorescence: A probe of photosynthesis in vivo. Annu. Rev. Plant Biol.; 2008; 59, pp. 89-113. [DOI: https://dx.doi.org/10.1146/annurev.arplant.59.032607.092759]

60. Meng, F.; Hong, S.; Wang, J.; Chen, A.; Zhang, Y.; Zhang, Y.; Janssens, I.A.; Mao, J.; Myneni, R.B.; Peñuelas, J. Climate change increases carbon allocation to leaves in early leaf green-up. Ecol. Lett.; 2023; 26, pp. 816-826. [DOI: https://dx.doi.org/10.1111/ele.14205]

61. Zhang, Y.; Joiner, J.; Alemohammad, S.H.; Zhou, S.; Gentine, P. A global spatially contiguous solar-induced fluorescence (CSIF) dataset using neural networks. Biogeosciences; 2018; 15, pp. 5779-5800. [DOI: https://dx.doi.org/10.5194/bg-15-5779-2018]

62. Zhang, L.; Shen, M.; Jiang, N.; Lv, J.; Liu, L.; Zhang, L. Spatial variations in the response of spring onset of photosynthesis of evergreen vegetation to climate factors across the Tibetan Plateau: The roles of interactions between temperature, precipitation, and solar radiation. Agric. For. Meteorol.; 2023; 335, 109440. [DOI: https://dx.doi.org/10.1016/j.agrformet.2023.109440]

63. Huang, Z.; Zhou, L.; Chi, Y. Spring phenology rather than climate dominates the trends in peak of growing season in the Northern Hemisphere. Glob. Chang. Biol.; 2023; 29, pp. 4543-4555. [DOI: https://dx.doi.org/10.1111/gcb.16758] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/37198735]

64. Dang, C.; Shao, Z.; Huang, X.; Zhuang, Q.; Cheng, G.; Qian, J. Vegetation greenness and photosynthetic phenology in response to climatic determinants. Front. For. Glob. Chang.; 2023; 6, 1172220. [DOI: https://dx.doi.org/10.3389/ffgc.2023.1172220]

65. Zhang, W.; Luo, G.; Yuan, X.; Li, C.; Xie, M.; Wang, Y.; Ma, X.; Shi, H.; Hamdi, R.; Hellwich, O. New data-driven method for estimation of net ecosystem carbon exchange at meteorological stations effectively increases the global carbon flux data. Methods Ecol. Evol.; 2023; 14, pp. 2449-2463. [DOI: https://dx.doi.org/10.1111/2041-210X.14188]

66. Chen, J.M.; Deng, F.; Chen, M. Locally adjusted cubic-spline capping for reconstructing seasonal trajectories of a satellite-derived surface parameter. IEEE Trans. Geosci. Remote Sens.; 2006; 44, pp. 2230-2238. [DOI: https://dx.doi.org/10.1109/TGRS.2006.872089]

67. Fischer, A. A model for the seasonal variations of vegetation indices in coarse resolution data and its inversion to extract crop parameters. Remote Sens. Environ.; 1994; 48, pp. 220-230. [DOI: https://dx.doi.org/10.1016/0034-4257(94)90143-0]

68. Du, H.; Wang, M.; Liu, Y.; Guo, M.; Peng, C.; Li, P. Responses of autumn vegetation phenology to climate change and urbanization at northern middle and high latitudes. Int. J. Appl. Earth Obs. Geoinf.; 2022; 115, 103086. [DOI: https://dx.doi.org/10.1016/j.jag.2022.103086]

69. White, M.A.; Thornton, P.E.; Running, S.W. A continental phenology model for monitoring vegetation responses to interannual climatic variability. Glob. Biogeochem. Cycles; 1997; 11, pp. 217-234. [DOI: https://dx.doi.org/10.1029/97GB00330]

70. Priyadarshi, N.; Pathak, S.; Chakraborty, D.; Chowdary, V.M.; Srivastav, S.K.; Kamalakannan, C.; Chockalingam, J.; Bandyopadhyay, S. Estimation of natural vegetation phenology metrics using time series EVI over Jharkhand state, India. J. Spat. Sci.; 2024; 69, pp. 721-743. [DOI: https://dx.doi.org/10.1080/14498596.2023.2281926]

71. Song, Y.; Wang, J.; Yu, Q.; Huang, J. Using MODIS LAI data to monitor spatio-temporal changes of winter wheat phenology in response to climate warming. Remote Sens.; 2020; 12, 786. [DOI: https://dx.doi.org/10.3390/rs12050786]

72. Zhou, J.; Jia, L.; Menenti, M. Reconstruction of global MODIS NDVI time series: Performance of Harmonic ANalysis of Time Series (HANTS). Remote Sens. Environ.; 2015; 163, pp. 217-228. [DOI: https://dx.doi.org/10.1016/j.rse.2015.03.018]

73. Xu, Y.; Shen, Y. Reconstruction of the land surface temperature time series using harmonic analysis. Comput. Geosci.; 2013; 61, pp. 126-132. [DOI: https://dx.doi.org/10.1016/j.cageo.2013.08.009]

74. Piao, S.; Fang, J.; Zhou, L.; Ciais, P.; Zhu, B. Variations in satellite-derived phenology in China’s temperate vegetation. Glob. Chang. Biol.; 2006; 12, pp. 672-685. [DOI: https://dx.doi.org/10.1111/j.1365-2486.2006.01123.x]

75. Beck, P.S.; Atzberger, C.; Høgda, K.A.; Johansen, B.; Skidmore, A.K. Improved monitoring of vegetation dynamics at very high latitudes: A new method using MODIS NDVI. Remote Sens. Environ.; 2006; 100, pp. 321-334. [DOI: https://dx.doi.org/10.1016/j.rse.2005.10.021]

76. Woodwell, G.M.; Whittaker, R.; Reiners, W.; Likens, G.E.; Delwiche, C.; Botkin, D. The Biota and the World Carbon Budget: The terrestrial biomass appears to be a net source of carbon dioxide for the atmosphere. Science; 1978; 199, pp. 141-146. [DOI: https://dx.doi.org/10.1126/science.199.4325.141]

77. Zhang, J.; Hao, X.; Hao, H.; Fan, X.; Li, Y. Climate change decreased net ecosystem productivity in the arid region of central Asia. Remote Sens.; 2021; 13, 4449. [DOI: https://dx.doi.org/10.3390/rs13214449]

78. Pei, Z.Y.; Ouyang, H.; Zhou, C.P.; Xu, X.L. Carbon balance in an alpine steppe in the Qinghai-Tibet Plateau. J. Integr. Plant Biol.; 2009; 51, pp. 521-526. [DOI: https://dx.doi.org/10.1111/j.1744-7909.2009.00813.x]

79. Zhang, K.; Zhu, C.; Ma, X.; Zhang, X.; Yang, D.; Shao, Y. Spatiotemporal variation characteristics and dynamic persistence analysis of carbon sources/sinks in the Yellow River Basin. Remote Sens.; 2023; 15, 323. [DOI: https://dx.doi.org/10.3390/rs15020323]

80. Lyu, J.; Fu, X.; Lu, C.; Zhang, Y.; Luo, P.; Guo, P.; Huo, A.; Zhou, M. Quantitative assessment of spatiotemporal dynamics in vegetation NPP, NEP and carbon sink capacity in the Weihe River Basin from 2001 to 2020. J. Clean. Prod.; 2023; 428, 139384. [DOI: https://dx.doi.org/10.1016/j.jclepro.2023.139384]

81. Wei, X.; Yang, J.; Luo, P.; Lin, L.; Lin, K.; Guan, J. Assessment of the variation and influencing factors of vegetation NPP and carbon sink capacity under different natural conditions. Ecol. Indic.; 2022; 138, 108834. [DOI: https://dx.doi.org/10.1016/j.ecolind.2022.108834]

82. Tang, X.; Fan, S.; Qi, L.; Guan, F.; Du, M.; Zhang, H. Soil respiration and net ecosystem production in relation to intensive management in Moso bamboo forests. Catena; 2016; 137, pp. 219-228. [DOI: https://dx.doi.org/10.1016/j.catena.2015.09.008]

83. Yao, Y.; Li, Z.; Wang, T.; Chen, A.; Wang, X.; Du, M.; Jia, G.; Li, Y.; Li, H.; Luo, W. A new estimation of China’s net ecosystem productivity based on eddy covariance measurements and a model tree ensemble approach. Agric. For. Meteorol.; 2018; 253, pp. 84-93. [DOI: https://dx.doi.org/10.1016/j.agrformet.2018.02.007]

84. Ziemblińska, K.; Urbaniak, M.; Chojnicki, B.H.; Black, T.A.; Niu, S.; Olejnik, J. Net ecosystem productivity and its environmental controls in a mature Scots pine stand in north-western Poland. Agric. For. Meteorol.; 2016; 228–229, pp. 60-72. [DOI: https://dx.doi.org/10.1016/j.agrformet.2016.05.022]

85. Ohtsuka, T.; Mo, W.; Satomura, T.; Inatomi, M.; Koizumi, H. Biometric Based Carbon Flux Measurements and Net Ecosystem Production (NEP) in a Temperate Deciduous Broad-Leaved Forest Beneath a Flux Tower. Ecosystems; 2007; 10, pp. 324-334. [DOI: https://dx.doi.org/10.1007/s10021-007-9017-z]

86. Wu, Z.; Huang, N.E. Ensemble empirical mode decomposition: A noise-assisted data analysis method. Adv. Adapt. Data Anal.; 2009; 1, pp. 1-41. [DOI: https://dx.doi.org/10.1142/S1793536909000047]

87. Pan, N.; Feng, X.; Fu, B.; Wang, S.; Ji, F.; Pan, S. Increasing global vegetation browning hidden in overall vegetation greening: Insights from time-varying trends. Remote Sens. Environ.; 2018; 214, pp. 59-72. [DOI: https://dx.doi.org/10.1016/j.rse.2018.05.018]

88. Sun, R.; Chen, S.; Su, H. Trend analysis and driving factors of vegetation dynamics in northern China from 1982 to 2015. Remote Sens.; 2022; 14, 6163. [DOI: https://dx.doi.org/10.3390/rs14236163]

89. Xiong, T.; Du, S.; Zhang, H.; Zhang, X. Satellite observed reversal in trends of spring phenology in the middle-high latitudes of the Northern Hemisphere during the global warming hiatus. Glob. Chang. Biol.; 2023; 29, pp. 2227-2241. [DOI: https://dx.doi.org/10.1111/gcb.16580]

90. Jiao, F.; Liu, H.; Xu, X.; Gong, H.; Lin, Z. Trend evolution of vegetation phenology in China during the period of 1981–2016. Remote Sens.; 2020; 12, 572. [DOI: https://dx.doi.org/10.3390/rs12030572]

91. Wei, F.; Wang, S.; Fu, B.; Pan, N.; Feng, X.; Zhao, W.; Wang, C. Vegetation dynamic trends and the main drivers detected using the ensemble empirical mode decomposition method in East Africa. Land Degrad. Dev.; 2018; 29, pp. 2542-2553. [DOI: https://dx.doi.org/10.1002/ldr.3017]

92. Wu, C.; Wang, X.; Wang, H.; Ciais, P.; Peñuelas, J.; Myneni, R.B.; Desai, A.R.; Gough, C.M.; Gonsamo, A.; Black, A.T. Contrasting responses of autumn-leaf senescence to daytime and night-time warming. Nat. Clim. Chang.; 2018; 8, pp. 1092-1096. [DOI: https://dx.doi.org/10.1038/s41558-018-0346-z]

93. Jeong, S.J.; HO, C.H.; GIM, H.J.; Brown, M.E. Phenology shifts at start vs. end of growing season in temperate vegetation over the Northern Hemisphere for the period 1982–2008. Glob. Chang. Biol.; 2011; 17, pp. 2385-2399. [DOI: https://dx.doi.org/10.1111/j.1365-2486.2011.02397.x]

94. Xia, J.; Chen, J.; Piao, S.; Ciais, P.; Luo, Y.; Wan, S. Terrestrial carbon cycle affected by non-uniform climate warming. Nat. Geosci.; 2014; 7, pp. 173-180. [DOI: https://dx.doi.org/10.1038/ngeo2093]

95. Li, C.; Wang, R.; Cui, X.; Wu, F.; Yan, Y.; Peng, Q.; Qian, Z.; Xu, Y. Responses of vegetation spring phenology to climatic factors in Xinjiang, China. Ecol. Indic.; 2021; 124, 107286. [DOI: https://dx.doi.org/10.1016/j.ecolind.2020.107286]

96. Ren, S.; Li, Y.; Peichl, M. Diverse effects of climate at different times on grassland phenology in mid-latitude of the Northern Hemisphere. Ecol. Indic.; 2020; 113, 106260. [DOI: https://dx.doi.org/10.1016/j.ecolind.2020.106260]

97. Wang, J.; Xi, Z.; He, X.; Chen, S.; Rossi, S.; Smith, N.G.; Liu, J.; Chen, L. Contrasting temporal variations in responses of leaf unfolding to daytime and nighttime warming. Glob. Chang. Biol.; 2021; 27, pp. 5084-5093. [DOI: https://dx.doi.org/10.1111/gcb.15777] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/34263513]

98. Lobell, D.B. Changes in diurnal temperature range and national cereal yields. Agric. For. Meteorol.; 2007; 145, pp. 229-238. [DOI: https://dx.doi.org/10.1016/j.agrformet.2007.05.002]

99. Liu, Y.; Wu, C.; Wang, X.; Jassal, R.S.; Gonsamo, A. Impacts of global change on peak vegetation growth and its timing in terrestrial ecosystems of the continental US. Glob. Planet. Chang.; 2021; 207, 103657. [DOI: https://dx.doi.org/10.1016/j.gloplacha.2021.103657]

100. Hou, X.; Gao, S.; Li, X.; Wang, M.; Zhao, J.; Sui, X.; Liang, S. Responses of vegetation phenology to the asymmetric changes of temperature in daytime and night-time in the north of 20° N. Int. J. Climatol.; 2022; 42, pp. 8146-8161. [DOI: https://dx.doi.org/10.1002/joc.7699]

101. Li, Y.; Zhang, W.; Schwalm, C.R.; Gentine, P.; Smith, W.K.; Ciais, P.; Kimball, J.S.; Gazol, A.; Kannenberg, S.A.; Chen, A. Widespread spring phenology effects on drought recovery of Northern Hemisphere ecosystems. Nat. Clim. Chang.; 2023; 13, pp. 182-188. [DOI: https://dx.doi.org/10.1038/s41558-022-01584-2]

102. Fu, Y.H.; Li, X.; Chen, S.; Wu, Z.; Su, J.; Li, X.; Li, S.; Zhang, J.; Tang, J.; Xiao, J. Soil moisture regulates warming responses of autumn photosynthetic transition dates in subtropical forests. Glob. Chang. Biol.; 2022; 28, pp. 4935-4946. [DOI: https://dx.doi.org/10.1111/gcb.16227] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/35642473]

103. Zhang, Y.; Parazoo, N.C.; Williams, A.P.; Zhou, S.; Gentine, P. Large and projected strengthening moisture limitation on end-of-season photosynthesis. Proc. Natl. Acad. Sci. USA; 2020; 117, pp. 9216-9222. [DOI: https://dx.doi.org/10.1073/pnas.1914436117] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/32284402]

104. Fensholt, R.; Langanke, T.; Rasmussen, K.; Reenberg, A.; Prince, S.D.; Tucker, C.; Scholes, R.J.; Le, Q.B.; Bondeau, A.; Eastman, R. et al. Greenness in semi-arid areas across the globe 1981–2007—An Earth Observing Satellite based analysis of trends and drivers. Remote Sens. Environ.; 2012; 121, pp. 144-158. [DOI: https://dx.doi.org/10.1016/j.rse.2012.01.017]

105. Fang, X.; Zhang, C.; Wang, Q.; Chen, X.; Ding, J.; Karamage, F. Isolating and quantifying the effects of climate and CO₂ changes (1980–2014) on the net primary productivity in arid and semiarid China. Forests; 2017; 8, 60. [DOI: https://dx.doi.org/10.3390/f8030060]

106. Ma, X.; Huete, A.; Moran, S.; Ponce-Campos, G.; Eamus, D. Abrupt shifts in phenology and vegetation productivity under climate extremes. J. Geophys. Res. Biogeosci.; 2015; 120, pp. 2036-2052. [DOI: https://dx.doi.org/10.1002/2015JG003144]

107. Konings, A.G.; Williams, A.P.; Gentine, P. Sensitivity of grassland productivity to aridity controlled by stomatal and xylem regulation. Nat. Geosci.; 2017; 10, pp. 284-288. [DOI: https://dx.doi.org/10.1038/ngeo2903]

108. Grossiord, C.; Buckley, T.N.; Cernusak, L.A.; Novick, K.A.; Poulter, B.; Siegwolf, R.T.W.; Sperry, J.S.; McDowell, N.G. Plant responses to rising vapor pressure deficit. New Phytol.; 2020; 226, pp. 1550-1566. [DOI: https://dx.doi.org/10.1111/nph.16485]

109. Li, X.; Zhao, C.; Kang, M.; Ma, M. Responses of net primary productivity to phenological dynamics based on a data fusion algorithm in the northern Qinghai-Tibet Plateau. Ecol. Indic.; 2022; 142, 109239. [DOI: https://dx.doi.org/10.1016/j.ecolind.2022.109239]

110. Chen, X.; Zhang, Y. Impacts of climate, phenology, elevation and their interactions on the net primary productivity of vegetation in Yunnan, China under global warming. Ecol. Indic.; 2023; 154, 110533. [DOI: https://dx.doi.org/10.1016/j.ecolind.2023.110533]

111. Liu, Z.; Wang, T.; Qu, Y.; Liu, H.; Wu, X.; Wen, Y. Prediction of high-quality MODIS-NPP product data. Remote Sens.; 2019; 11, 1458. [DOI: https://dx.doi.org/10.3390/rs11121458]