Spatiotemporal Patterns of Land Surface Phenology from 2001 to 2021 in the Agricultural Pastoral Ecotone of Northern China

Abstract

Vegetation phenology is one of the most sensitive indicators to understanding terrestrial ecosystem status and change. However, few studies have been conducted to reveal vegetation phenology variation characteristics over the past two decades, especially under the background of the global warming hiatus since 1998. The agricultural pastoral ecotone of northern China (APENC) is an ideal place to analyze land surface phenology (LSP) variation. Therefore, the spatiotemporal patterns of LSP were quantitatively analyzed at regional, basin and pixel scales based on time-series MODIS NDVI data (2001–2021). Results showed that: (i) The start of the growing season (SOS) occurred in 105–141 Julian days, the end of the growing season (EOS) was between 257 and 281 Julian days and the length of the growing season (LOS) varied from 130 to 172 days. The later SOS was mainly distributed in croplands and typical grassland areas, while the early SOS was observed in forests and sandy vegetation coverage areas. The early EOS occurred in typical grasslands, and the later EOS was concentrated in the southeast boundary. The magnitude of the SOS and LOS fluctuation was less than EOS. (ii) The SOS and EOS exhibited overall insignificant advanced and delayed trends at a rate of −0.09 days·yr⁻¹ and 0.12 days·yr⁻¹, respectively, and the LOS displayed an insignificant extended trend at a rate of 0.26 days·yr⁻¹ at a regional scale. The trends of phenological metrics were consistent with the APENC in the Yellow River and Haihe River Basins. The shortened trend of LOS occurred due to the delayed SOS and advanced EOS in the Songliao River and Continental Basins. (iii) The SOS variation gradually changed from an advanced trend to a delayed trend from a southwest to northeast direction in cropland and grassland ecosystems, whereas an opposite trend was found for EOS. The LOS exhibited a significant extended trend due to the significant advanced and delayed trend of SOS and EOS at p < 0.01 in forest ecosystems. This work provides a critical reference for the vegetation phenology dynamic research of semi-arid and semi-humid regions.

1. Introduction

Climate change has become the largest non-traditional security challenge facing human development [1]. The Intergovernmental Panel on Climate Change Sixth Assessment Report pointed out that the global average temperature has increased by 1.09 °C compared with pre-industrial levels [2]. Climate warming has had a profound influence on the structure and function of terrestrial ecosystems and has become the focus of scientists and the public [3]. Vegetation phenology refers to the periodic life activities of plants driven by climate and other environmental factors [4]. It is one of the most sensitive indicators to understanding terrestrial ecosystem status and change, and its change is also the primary manifestation of plant adaptation to climate change [5]. Many studies have shown that the occurrence and duration of phenological events has changed dramatically because of global warming, leading to a significant influence on global carbon cycle, water cycle, energy flows and even human health [6,7,8]. Therefore, obtaining vegetation phenological dynamics can help us to understand the evolution of plant communities and the behavioral responses of the Earth’s ecosystems in the face of climate change [9].

A direct phenological record of individual species or a dominant community life cycle, based on ground observation and phenological cameras, is a major means to analyzing vegetation phenology changes at specie and plot scales [8,10,11,12]. Ground-based manual observations provide very accurate long-term phenology data for individual plant species, but have limitations in terms of spatial coverage and time resolution [13,14]. Phenological cameras monitor continuous vegetation changes at canopy and specie scales by configuring the visible, multispectral and hyperspectral automated cameras [15,16,17]. This approach has become a powerful tool for phenology research because of the less costly, easy to set up and high frequent real-time observation. However, the spatial distribution of ground sites is very rare in China. The Chinese Phenological Observation Network shows that the phenological information of thirty sites have been recorded since 1963 (http://www.cpon.ac.cn/, accessed on 2 December 2022). Phenological camera sites are mainly distributed in America, Europe, Japan and Australia, and the phenological data of only two sites (the Jurong and Shangqiu sites) can be accessed within China (https://phenocam.nau.edu/, accessed on 10 January 2023). Thus, it is difficult to reveal the integrative phenology patterns in continuous space due to the spatial discontinuity of the ground observation.

Satellite remote sensing has been commonly utilized to monitor the dynamic change of vegetation phenology since it has the advantage of obtaining spatial continuous information. Land surface phenology (LSP) refers to vegetation phenology derived from satellite data. Long time-series satellite data mainly come from sensors such as NOAA-AVHRR, Terra/Aqua-MODIS and SPOT-VGT, since 1980 [4,18,19]. These datasets are used to analyze vegetation phenology dynamics at regional and global scales [20,21]. Landsat and Sentinel data are used for LSP research at a local scale [22]. Phenological metrics, including the start of the growing season (SOS), end of the growing season (EOS) and length of growing season (LOS), can be identified based on LSP extraction methods and various time-series remote sensing data [23].

During the past few decades, many phenological metrics extraction methods, including thresholding, autoregressive moving average, derivative, machine learning and shape model fitting, have been developed based on optical vegetation indices, microwave vegetation opacity depth and solar induced chlorophyll fluorescence [24,25,26,27]. Among the vegetation indices, AVHRR GIMMS NDVI and MODIS NDVI are the most used data sources for large-scale phenological metrics simulation. The AVHRR, with a coarse spatial resolution (8 km) and 10-day composites, carried out LSP studies at global, hemisphere or continental scales since 1982. However, the LSP parameters may not be accurately estimated in the high heterogenous vegetation coverage areas based on AVHRR data due to the coarse resolution, data quality and data missing problems [28]. The MODIS, available from 2000 to the present, and at a 250 m–5.6 km spatial resolution with 16-day composites, is considered more appropriate for estimating phenological metrics and revealing the spatiotemporal patterns of LSP in the past two decades due to its moderate spatial resolution, high observation frequency and data quality. Previous LSP studies are mainly focused on the temperature sensitivity of spring phenology and the response of forest phenology in middle and high latitudes to climate change [3,8,13,29,30]. However, the research on EOS, LOS and the phenological changes of grassland and cropland ecosystems in middle latitude areas is not sufficient. Meanwhile, there is a lack of consideration of land use change in most LSP studies, which may not accurately reflect the spatiotemporal evolution characteristics of phenology in different ecosystems. In addition, previous LSP research mainly focuses on the period from 1982 to 2010, based on the AVHRR GIMMS data due to the limitation of data availability. How LSP changed in the last 20 years should be addressed, especially in the context of the global warming hiatus since 1998 [31].

The agricultural pastural ecotone of northern China (APENC) is an ecological transition zone connecting agricultural cultivation and animal husbandry [32]. To protect the north China plain and Beijing-Tianjin-Hebei region against a sandstorm invasion, a great deal of ecological restoration projects have been performed by the Chinese government since 1999 [33]. The land use has changed dramatically, vegetation productivity has been increased significantly and vegetation has shown a greening trend under the influence of man-made protection and a warming-drying climate trend over the APENC [34]. Against this background, there is an urgent need to reveal the spatiotemporal patterns and trend characteristics of LSP in the past two decades. This work will provide a critical reference for the vegetation phenology dynamic research of semi-arid and semi-humid regions. The purposes of this study are to: (i) analyze the spatiotemporal patterns and trend characteristics of LSP quantitatively at regional, basin and pixel scales from 2001 to 2021, and (ii) reveal the phenological evolution process of different ecosystems.

2. Materials and Methods

2.1. Study Area

APENC is situated in the southeast edge of the Mongolia Plateau and northern part of the Loess Plateau, between 33.5°–48.6° N and 101°–126.5° E (Figure 1 [35]). The APENC covers an area of 8.35 × 10⁵ km², and more than 50% of the area is grassland. The mean annual precipitation ranges from 300 to 450 mm, and about 60–70% of precipitation is mainly concentrated in summer. Precipitation variability belongs to 15–30%. The mean annual evaporation ranges from 1600 to 2500 mm, and the frost-free period lasts 130–160 days. The land use pattern is characterized by the interlacing of agriculture and animal husbandry affected by the monsoon climate [36,37].

2.2. Datasets and Pre-Processing

2.2.1. MODIS Data

We used the 16-day composite data of the MODIS NDVI product (MOD13Q1 V006) to estimate the key LSP parameters, with the spatial resolution of 231.7 m from 2001 to 2021. The preprocessing such as image mosaic and clipping, projection and data format conversion, etc., were carried out using the MODIS Reprojection Tool. However, there are still many singular values, although the pre-composited eight-day surface reflectance data were used to generate the MOD13Q1 dataset. In this study, the harmonic analysis of time series approach was utilized to reconstruct cloud-free time-series MODIS NDVI data [38]. The optimal harmonic parameter combination suited for reconstructing high-quality MODIS NDVI data over the APENC has been determined by Wei et al. (2020) [39]. Furthermore, areas such as bare rock, Gobi, desert and waterbody were masked based on the threshold method due to no or sparse vegetation coverage, and 98% of study area was finally considered.

2.2.2. Basin and Land Use Data

The basin data was derived from the data center for resource and environmental sciences of Chinese Academy of Sciences (https://www.resdc.cn/, accessed on 15 January 2023), and the APENC was divided into four basins: Yellow River Basin, Haihe River Basin. Songliao River Basin and Continental Basin (Figure 1). The GlobeLand30 data, with a spatial resolution of 30 m in 2000, 2010 and 2020, were used to reflect the ecosystem situation [40]. The main land use types are grasslands, croplands and forests in AEPNC, accounting for more than 96% of the study area [35]. Thus, this study mainly focused on the change characteristics of the vegetation phenology in these three ecosystems. To avoid the impact of land use change on analysis results, the Globeland30 in 2000, 2010 and 2020 were firstly upscaled to 231.7 m using the preponderant area method, and then the unchanged land use type pixels were chosen based on the spatial overlay analysis method, and the corresponding phenological metrics of grasslands, croplands and forests were finally extracted.

2.2.3. Plant Phenological Observation Dataset

The vegetation phenology observation dataset was collected from the Chinese Ecosystem Research Network Data Center (CERN, http://www.nesdc.org.cn/, accessed on 20 December 2022) and was used to evaluate the accuracy of the SOS and EOS parameters. This dataset records the phenological observation information of more than 660 plant species at 21 ecological sites during the period of 2003–2015 [41]. The dataset was divided into a woody plant subset and an herbaceous plant subset because of the distinct differences between the phenological stages of woody and herbaceous plants. The woody plant subset reflects the leaf bud breaking phase, leaf unfolding phase, first bloom phase, full flowering phase, fruit or seed ripening phase, leaf turning to autumn color phase and leaf falling phase. The herbaceous plant subset reflects the germination or turning green phase, flowering phase, fruit or seed ripening phase, seed dispersal phase and autumn wilting phase. Four ecological sites that are distributed within the boundary of the study area were selected, including Inner Mongolia, Naiman, Ordos and Beijing sites (Figure 1 and Table 1). In this study, the phenological information of dominant species were used to reflect the basic status of plant phenology in four sites.

2.3. Phenological Metrics Extraction and Validation

2.3.1. Phenological Metrics Extraction

The cumulative NDVI-based logistic regression curves, as one of the derivative methods, was proven to have a higher accuracy for extracting the grassland phenology in the Mongolia Plateau and Northern China [42,43]. This approach was put forward on the basis of the piecewise logistic function by Zhang et al. (2003) [24], and it has an advantage since it does not identify the growth or senescence phase modelling the NDVI curve using a logistic function [42]. Therefore, the cumulative NDVI-based logistic regression curves approach was used to identify the LSP parameters.

The cumulative NDVI curve is calculated through a logistic regression function:

(1)${{\displaystyle y}}_t=\frac{c}{1+{{\displaystyle e}}^{a+bt}}+d$

where t indicates the Julian day, y_t indicates the fitted NDVI value at time t, d indicates the initial background value of NDVI, and is replaced by the minimum NDVI value in a year in this study, c + d indicates the maximum value of the cumulative NDVI and a and b are fitting coefficients and can be determined based on the Levenberg-Marquardt method. To realize the matching between simulated phenology and ground observation data, the 16-day NDVI data was resampled to daily NDVI data using a linear interpolation algorithm.

The first curvature for Equation (1) at time t is calculated as follows:

(2)$K=\frac{b^{2}cz(1-z)({{\displaystyle 1+z)}}^3}{{{\displaystyle [{(1+z)}^4+{(bcz)}^2]}}^{1.5}}$

(3)$z=e^{a+bt}$

where K is the first curvature, and the local maximum and minimum value of first curvature correspond to the date of SOS and EOS, respectively.

2.3.2. Accuracy validation

Three indicators, including the root mean squared error (RMSE), determination coefficient (R²) and relative root mean squared error (RRMSE), were computed to validate the accuracy of simulated phenology [44]. The formulas are shown:

(4)${\mathrm{R}}^2=1-\frac{{\displaystyle \sum _{i=1}^n(Pob{e}_i-Pex{t}_i}{)}^2}{{\displaystyle \sum _{i=1}^n{(Pob{e}_i-\overline{Pobe})}^2}}$

(5)$\mathrm{RMSE}=\sqrt{{\displaystyle \sum _{i=1}^n\frac{{(Pex{t}_i-Pob{e}_i)}^2}{n}}}$

(6)$\mathrm{RRMSE}=(1-\frac{RMSE}{\overline{Pobe}})\times 100\%$

where n indicates the number of samples, Pobe_i and Pext_i indicate the observed and simulated values, respectively, and $\overline{Pobe}$ indicates the average of all the observation samples.

2.4. Data Analyses

The multiyear average and standard deviation of phenological metrics were computed at pixel level to reflect the spatial distribution and fluctuation degree of vegetation phenology. The temporal, trend characteristics and phenological dynamics of different ecosystems were analyzed at regional, basin and pixel scales based on two non-parametric methods—Sen’s slope estimator and the Mann-Kendall trend test. These two methods are neither affected by outliers nor in need of data to obey a certain distribution, and they have a good avoidance ability for measurement errors or discrete data [45].

2.4.1. Sen’s Slope Estimator

Sen’s slope estimator is a non-parametric procedure that involves choosing the median of slopes, which was developed by Sen (1968) [46]. The trend detection of vegetation phenology was completed using Sen’s slope. The formula is as follows:

(7)${{\displaystyle \beta }}_i=Median(\frac{{{\displaystyle x}}_j-{{\displaystyle x}}_k}{j-k}) for i=1,2, \cdots \cdots N$

where β_i indicates the Sen’s slope and x_j and x_k indicate the NDVI data values at time j and k (j > k), respectively.

2.4.2. Mann-Kendall Trend Test

The Mann-Kendall test (M-K test) was used to identify the trend significance of vegetation phenological variation. The M-K test statistic S is computed as:

(8)$S={\displaystyle \sum _{i=1}^{n-1}{\displaystyle \sum _{j=i+1}^{n}sgn({{\displaystyle x}}_j-{{\displaystyle x}}_i)}}$

where x_i and x_j stand for NDVI values at time i and j (j > i), respectively, n indicates the number of NDVI values and sgn(x_j − x_i) is the sign function:

(9)$sgn({{\displaystyle x}}_j-{{\displaystyle x}}_i)=\{\begin{array}{ll}+1& if{{\displaystyle x}}_j>{{\displaystyle x}}_i\\ 0& if{{\displaystyle x}}_j={{\displaystyle x}}_i\\ -1& if{{\displaystyle x}}_j<{{\displaystyle x}}_i\end{array}$

The variance Var(s) is calculated as:

(10)$Var(S)=\frac{n(n-1)(2n+5)-{\displaystyle \sum _{p=1}^g{{\displaystyle t}}_p}({{\displaystyle t}}_p-1)({{\displaystyle 2t}}_p+5)}{18}$

where n indicates the number of NDVI values, g indicates the number of tied groups and t_p indicates the number of ties of extent p. The test statistic Z is calculated as:

(11)$Z=\{\begin{array}{ll}\frac{S-1}{\sqrt{Var(S)}}& if S>0\\ \hfill 0\hfill & if S=0\\ \frac{S+1}{\sqrt{Var(S)}}& if S<0\end{array}$

The Z value with a positive sign indicates an increasing trend, and a negative sign indicates decreasing trends. In this study, testing trends were conducted when the significance level was set to 0.05 and 0.01, respectively. The significant trend was accepted if |Z| > 1.96 at 0.05 significance level and accepted if |Z| > 2.576 at 0.01 significance level. All analyses were performed in ArcGIS 10.2, ENVI/IDL 5.3 and MATLAB R2017a software.

3. Results

3.1. Validation of Land Surface Phenology

Overall, the extracted results of phenological metrics were highly consistent with ground observation using the cumulative NDVI-based logistic regression curves method (Figure 2). Results showed that the determination coefficient (R²), RMSE and RRMSE between the simulated SOS and ground observation value were 0.678, 8.7 days and 92.37%, respectively (Figure 2a); while the three indicators corresponding to EOS were 0.695, 5.8 days and 97.84% (Figure 2b). This meant that the extracted phenological metrics could be used to analyze the spatiotemporal pattern of vegetation phenology over the APENC. Meanwhile, we found that the simulation accuracy of EOS was obviously higher than that of SOS, although the SOS, EOS and ground observed phenology showed significant correlation at p < 0.01.

The extraction precision of phenological metrics in different ecosystems was further assessed (Figure 3). We could see that the grassland-desert ecosystem had the highest accuracy, followed by the grassland ecosystem, and the forest ecosystem was relatively inferior. There were significant correlations between SOS, EOS and ground observation at p < 0.01 and p < 0.05 in the grassland-desert ecosystem (Naiman and Ordos sites; Figure 3a,b). The values of RMSE and RRMSE were 6.1 days (SOS), 6.2 days (EOS), 94.72% (SOS) and 97.72% (EOS). No significant correlations (p > 0.05) were identified between SOS, EOS and ground observed phenology in the grassland ecosystem (Inner Mongolia site) with a RMSE and RRMSE of 9.1 days and 93.1% for SOS and 5.8 days and 97.76% for EOS, respectively (Figure 3c,d). The RMSE and RRMSE were 11.6 days and 88.76% for SOS and 5.3 days and 98.06% for EOS in the forest ecosystem (Beijing site), and the simulated values for SOS and EOS showed significant correlations at p < 0.05 and p < 0.01 with measured phenology, respectively (Figure 3e,f).

3.2. Spatiotemporal and Trend Variation of Land Surface Phenology

3.2.1. Spatial Pattern of Land Surface Phenology

Figure 4 shows the distribution of the average and standard deviation of LSP parameters from 2001 to 2021 over the APENC. It was found that about 91.5% of the average SOS primarily occurred between 105 to 141 Julian days, that is, from April 15 to May 20, and the SOS variation showed a significant spatial heterogeneity (Figure 4a). The later SOS (>129 Julian days) was mainly distributed in the typical grasslands to the west of Greater Khingan and croplands areas, including the Songliao Plain in the northeast, Xinzhou, Datong and Xuanhua Basins, Tumochuan Plain and the hinterland of the Loess Plateau to the west of Lvliang Mountain in the middle, along Qingshui river in the southwest of the study area. The relatively early SOS (<111 Julian days) was found in the forest coverage areas such as the Greater Khingan and Yanshan Mountains in the northeast, Lvliang Mountain in the middle, Liupan and Lianhua Mountains in the southwest of APENC and the sparse vegetation coverage areas located in Horqin and Mu Us Sandy Lands. The spatial distribution of the standard deviation showed that the fluctuation degree of SOS was in the range of 3–15 days, accounting for 96.5% of the study area in the recent 20 years (Figure 4b). The magnitude of SOS fluctuation exceeded 15 days in part of the downstream of Huolin River and the northern margin of the Yinshan Mountain regions.

The average EOS occurred between 257 to 281 Julian days (from September 14 to October 8; Figure 4c) with the counting of pixels for 95.4%. The early EOS (<266 Julian days) mainly occurred in the typical grassland to the west of Greater Khingan and the north of the Yanshan Mountain. The later EOS (>278 Julian days) was mainly distributed in the southeast boundary of APENC, such as forest coverage areas in Yanshan Mountain, croplands and grassland coverage areas in the Yellow River Basin. The magnitude of EOS fluctuation was relatively small compared with SOS, and the fluctuation degree with 89.5% of the research area ranged from 3 to 9 days (Figure 4d).

The average LOS varied from 130 to 172 days, accounting for 93.6% of APENC (Figure 4e). The values of LOS were controlled by SOS and EOS. The LOS was less than 130 days in the typical grassland areas with later SOS and early EOS, and the LOS was primarily concentrated within 130 to 151 days in the croplands area with later SOS. The LOS exceeded 169 days in forest coverage areas with early SOS. The fluctuation degree of LOS ranged from 3 to 15 days with 95% of study area, and the LOS fluctuated more than 15 days in desert steppe in the north of Yinshan Mountain (Figure 4f).

3.2.2. Temporal Trend of Land Surface Phenology at Regional and Basin Scales

Figure 5 shows the temporal trend variation of phenological metrics at regional and basin scales from 2001 to 2021. We could clearly see that the SOS exhibited an advanced trend (slope = −0.09 days·yr⁻¹), and a delay trend occurred in EOS (slope = 0.12 days·yr⁻¹), leading to an extended trend for LOS (slope = 0.26 days·yr⁻¹) over the APENC (Figure 5e). However, the trend variation of phenological metrics was not significant (p > 0.05) for the entire whole period. The temporal trends of SOS, EOS and LOS in the Yellow River and Haihe River Basins were consistent with that of the whole study area, and the significant delayed EOS and extended LOS trends were observed with a rate of 0.11 days·yr⁻¹ (p < 0.05) and 0.39 days·yr⁻¹ (p < 0.01) in the Yellow River Basin, and 0.28 days·yr⁻¹ (p < 0.01) and 0.63 days·yr⁻¹ (p < 0.01) in the Haihe River Basin, respectively (Figure 5a,b). However, the change of LSP parameters in the Songliao River and Continental Basins was opposite to that of the whole study area (Figure 5c,d). The LOS exhibited a shortened trend because of the delayed SOS and advanced EOS. Nevertheless, this trend did not pass the significance test at a 0.05 level.

3.2.3. Trend Change Characteristic of Land Surface Phenology at Pixel Level

The trend characteristic and significance of phenological metrics were calculated at pixel level based on Sen’s slope and M-K test methods (Figure 6). From 2001 to 2021, the change rate of SOS ranged from −2.18 days·yr⁻¹ to 1.99 days·yr⁻¹. The mean value was −0.14 days·yr⁻¹, and it generally showed an advanced trend. The area of advanced SOS accounted for 60.47% of APENC. Among them, the SOS exhibited an extremely significant and significant advancing trend with the area of 13.02% and 8.1%, respectively. These regions were mainly distributed in the southeast boundary of APENC. The advanced trend of SOS reached −0.85 days·yr⁻¹ (p < 0.01) and −0.60 days·yr⁻¹ (p < 0.05). Meanwhile, the SOS with an area of 2.83% and 1.9% displayed a significant and extremely significant delaying trend, and the delayed rates reached 0.67 days·yr⁻¹ (p < 0.05) and 0.88 days·yr⁻¹ (p < 0.01), respectively. These regions mainly occurred in Horqin Sandy Land to the east of the Great Khingan, Mu Us Sandy Land to the south of the Yellow River and Hulu River Basin to the west of Liupan Mountain.

The change rate of EOS ranged from −1.25 days·yr⁻¹ to 1.93 days·yr⁻¹, and the average was 0.15 days·yr⁻¹. The EOS with an area of 68.54% was greater than 0, which indicated that the EOS generally showed a delayed trend over the APENC. The EOS with areas of 7.52% and 8.29% displayed a significant and extremely significant delayed trend, and the delayed rates were 0.43 d·yr⁻¹ (p < 0.05) and 0.62 d·yr⁻¹ (p < 0.01), respectively. From the perspective of spatial distribution, these areas were mainly concentrated in the southeast boundary, that is, the area extending from northeast to southwest to the east of Lvliang Mountain, the hinterland of Mu Us Sandy Land and the Qingshui, Zhuanglang and Taohe River Basins in the southwest. In addition, the EOS had a significant advanced trend at a rate of −0.47 days·yr⁻¹ (p < 0.05) with the area of 1.37%, which was mainly scattered in Horqin Sandy Land and the Songhua River Basin in the northeast of APENC.

The LOS tended to an extended trend due to the advanced SOS and delayed EOS. The change rate of LOS ranged from −2.01 day·yr⁻¹ and 2.99 days·yr⁻¹, and the average was 0.27 days·yr⁻¹. The area of extended LOS accounted for 64.86% of the study area, and the LOS with areas of 8.46% and 22.57% displayed a significant and extremely significant extended trend. The extended rates were 0.67 days·yr⁻¹ (p < 0.05) and 1.04 days·yr⁻¹ (p < 0.01), which mainly occurred in the Great Khingan, along southeast boundary, Tumochuan Plain and part of a region southwest of the study area. The significant advanced LOS was basically consistent with the regions of significant delayed SOS. These areas accounted for 2.28% and 1.52% of the study area, and the advanced rate of LOS were −0.78 days·yr⁻¹ and −0.99 days·yr⁻¹.

3.3. Land Surface Phenology Variation in Relation to Different Ecosystems

The change characteristics of phenological metrics were computed in different ecosystems at regional and basin scales in the recent 20 years (Table 2). The SOS of the cropland ecosystem displayed a no-significant advanced trend at a regional scale with a rate of −0.053 days·yr⁻¹ (p > 0.05). At the basin scale, we found that the SOS of croplands changed from an advanced trend to a delayed trend from a southwest to northeast direction (from the Yellow River Basin to the Songliao River Basin). The SOS was advanced with a rate of −0.182 days·yr⁻¹ (p < 0.05) and −0.16 days·yr⁻¹ (p > 0.05) in the Yellow River and Haihe River Basins, respectively. However, the opposite trend occurred in the Songliao River and Continental Basins with a rate of 0.059 days·yr⁻¹ and 0.165 days·yr⁻¹, respectively. The EOS of croplands showed a no-significant delayed trend with a rate of 0.065 days·yr⁻¹ (p > 0.05) at a regional scale. At the basin scale, the EOS gradually changed from an extremely significant delayed trend of 0.196 days·yr⁻¹ (p < 0.01) to a no-significant advanced trend of −0.014 days·yr⁻¹ from a southwest to northeast direction. The LOS tended to an extended trend because of advanced SOS and delayed EOS in APENC, Yellow River and Haihe River Basins. The LOS tended to a shortened trend due to the delayed SOS and advanced EOS in the Songliao River. The LOS also tended to a shortened trend since the delayed rate of SOS was faster than the delayed rate of EOS in the Continental Basin.

The SOS variations of the grassland ecosystem were consistent with the cropland ecosystem at regional and basin scales. The EOS exhibited an advanced trend in APENC, Songliao and Continental Basins while the delayed trend was observed in the Yellow River and Haihe River Basins. Meanwhile, we noticed that the EOS of the grassland ecosystem changed from a delayed trend with a rate of 0.07 days·yr⁻¹ to an advanced trend with a rate of −0.01 days·yr⁻¹ from a southwest to northeast direction, which was consistent with the EOS variation of croplands. The change of LOS tended to a shortened trend in the Songliao River and Continental Basin while a delayed trend occurred in APENC, Yellow River and Haihe River Basins.

The SOS of the forest ecosystem exhibited an advanced trend, and the EOS exhibited a delayed trend, leading to an extended trend for LOS at regional and basin scales. The trends passed the significant test at the 0.01 level. The forests were mainly concentrated in the Lianhua, Liupan, Yinshan and Lvliang Mountains in the Yellow River Basin with the SOS, EOS and LOS of −0.77 days·yr⁻¹, 0.32 days·yr⁻¹ and 1.1 days·yr⁻¹, respectively. The change rate of the SOS, EOS and LOS were −0.58 days·yr⁻¹, 0.41 days·yr⁻¹ and 0.98 days·yr⁻¹, which was observed in the Wutai, Hengshan and Yanshan Mountains in the Haihe River Basin. The forests were mainly distributed in the Great Khingan, Songling, Nuluerhu and Yiwulv Mountains with the SOS, EOS and LOS of −0.63 days·yr⁻¹, 0.26 days·yr⁻¹ and 0.87 days·yr⁻¹, respectively.

4. Discussion

The vegetation phenology is known as “canary in a coal mine”, and it can directly reflect the characteristics of climate change and is a comprehensive indicator to measure environmental change [47]. The vegetation coverage of APENC has changed dramatically under the influence of human activities and the warming-drying climate trend [34,48]. In this context, how the vegetation phenology changed in spatiotemporal and trend characteristics in the recent 20 years need to be further explored, especially since the global warming hiatus in 1998.

The simulated phenology was verified using ground observation data, and it was seen that the cumulative NDVI-based logistic regression curves approach could achieve the extraction of vegetation phenology information with high precision over the APENC. The simulation of EOS was better than that of SOS. This result was consistent with the finding of the analysis of the grassland phenology in northern China based on SPOT VEGETATIONT NDVI data, conducted by Hou et al. (2014) [42]. Climate change, especially temperature change, has a significant influence on vegetation phenology [5,18,44]. The climate change characteristic of the study area is that the air temperature fluctuates greatly in spring, and the SOS varied considerably in different years, affecting the simulation results of remote sensing. The fluctuation of air temperature is small in autumn, and the climate conditions are relatively stable. So, the EOS could be obtained with higher accuracy [49,50]. Furthermore, the 16-day time-series MODIS NDVI data were resampled to the daily data through the linear interpolation method in order to realize the matching between simulated and ground observation data. The vegetation change, however, is a nonlinear process, which might increase the simulation error of the SOS to some extent [51,52].

The areas with later SOS were mainly distributed in croplands and typical grasslands. The phenological variations of croplands are more complex than that of natural vegetation [53]. The majority of croplands belong to a typical rain-fed agriculture over the APENC. The crops’ sowing date was usually chosen before the rainy season (from late April to early May), and was also affected by human farming activities such as planting structure and crop layout, resulting in the later SOS [35]. In typical grassland areas, we noted that SOS appeared later while EOS appeared early, making the LOS less than 130 days. This is because the region is affected by the southeast movement of cyclones and anticyclones in the mid-latitude westerly belt of the East Asian continent, and blocked by the Great Khingan and Yinshan Mountains, which brings abundant snowfall in winter and early spring, accompanied by the appearance of strong wind and cold wave weather [54,55]. At the same time, the frost weather comes earlier in autumn, which makes the temperature stay at a low level for a long time, leading to such phenological changes [56].

The early SOS and later EOS were observed in forest coverage areas. The forests were mainly concentrated in the mountains at the southeast boundary of APENC. The distribution of surface heat is sufficient in these areas, and the land surface temperature has increased significantly over the past two decades [57]. The snow cover in winter provides sufficient water conditions for spring vegetation green-up, leading to the early SOS. The East Asian monsoon climate has brought abundant precipitation in summer and autumn, which might be the reason for later EOS [58]. Meanwhile, the early SOS also occurred in sandy vegetation coverage areas. The reason was that the sandy vegetation had a higher sensitivity to precipitation, and a small amount of precipitation could promote the rapid green-up of vegetation [19].

The estimated key LSP parameters indicated that the LOS was extended due to the advanced SOS and delayed EOS at a regional scale. The results were consistent with the conclusion analysis of the vegetation phenology change at global and continental scales based on ARHRR GIMMS NDVI data [28,59,60]. However, the trends of phenological metrics did not pass the significant test at a 0.05 level. Wei et al. (2021) [57] analyzed the variation trend of land surface temperature based on MYD11A2 data product over the APENC during the period of 2003–2020 and found that the land surface temperature exhibited a no-significant warming trend with the rate of 0.007 °C·yr⁻¹. This might be the reason why the change in the vegetation phenology trend was not significant in this area.

The variation of vegetation phenology in different ecosystems varied considerably at the basin scale. The SOS gradually changed from an advanced trend to a delayed trend from southwest to northeast in the cropland and grassland ecosystems, whereas a contrary trend was found for EOS. These trends were also observed by Tao et al. (2011) [61] when analyzing the phenological dynamics of croplands in central and eastern Inner Mongolia based on MODIS EVI data. However, the driving mechanism is still unclear.

The phenological variation of the forest ecosystem was consistent at regional and basin scales, and the trends passed the significant test at a 0.01 level. The results were consistent with the conclusion analysis of forests phenology in northern China based on AVHRR GIMMS and MODIS NDVI data, respectively, conducted by Zhao et al. (2023) [62] and Zheng et al. (2022) [44]. For example, Zhao et al. (2023) [62] argued that the SOS and EOS of the forest ecosystem were significantly advanced and delayed by −0.18 days·yr⁻¹ and 0.14 days·yr⁻¹ at p < 0.01, respectively, resulting in a significant extended (p < 0.01) for LOS with a rate of 0.32 days·yr⁻¹ during 1982–2015. In addition, we noted that the magnitude of forest phenological metrics change was obviously higher than that of cropland and grassland ecosystems. This partly reflects that forest ecosystems might be more sensitive to climate change.

In this study, we systematically revealed the space, time and trend characteristics of vegetation phenology and its variations in different ecosystems over the APENC during the period of 2001–2021. However, the impact of climate change (e.g., temperature, rainfall and soil moisture, etc.) and human activities on LSP variations cannot be analyzed, which could limit the interpretability and relevance of the results. Therefore, it is necessary to make further efforts to discuss the contribution of climate change and human activities on LSP variations in future studies. In addition, the potential impacts of LSP variations on ecosystem services, such as carbon sequestration, water regulation and biodiversity conservation, will be considered in the future.

5. Conclusions

We systematically analyzed the spatiotemporal patterns and trend characteristics of LSP based on time-series Terra/MODIS NDVI data at regional, basin and pixel scales over the past two decades (2001–2021). The major findings are as follows:

(1) There was high consistency between simulated phenology and ground observation. The determination coefficient, RMSE and RRMSE were 0.678, 8.7 days and 92.37% for SOS and 0.697, 5.8 days and 97.84% for EOS, respectively. The simulation accuracy of EOS was obviously higher than that of SOS.

(2) The SOS of vegetation occurred in 105–141 Julian days, the EOS was mainly between 257 and 281 Julian days and the LOS varied from 130 to 172 days. The later SOS was distributed in croplands and typical grassland areas, while the early SOS was observed in forests and sandy vegetation coverage areas. The early EOS occurred in typical grasslands, and the later EOS was mainly concentrated in the southeast boundary of APENC. The magnitude of SOS and LOS fluctuation ranged from 3 to 15 days, and the fluctuation degree of EOS was relatively small, ranging from 3 to 9 days.

(3) The SOS and EOS exhibited overall insignificant advanced and delayed trends at a rate of −0.09 days·yr⁻¹ and 0.12 days·yr⁻¹ at a regional scale, and the LOS showed an insignificant extended trend at a rate of 0.26 days·yr⁻¹. The trends of phenological metrics were consistent with the APENC in the Yellow River and Haihe River Basins. While an opposite trend was observed in the Songliao River and Continental Basins, the shortened trend of LOS occurred due to the delayed SOS and advanced EOS.

(4) The SOS variation gradually changed from an advanced trend to a delayed trend from southwest to northeast in cropland and grassland ecosystems, whereas an opposite trend was found for EOS. The LOS exhibited a significant extended trend at p < 0.01 due to the significant advanced and delayed trend of SOS and EOS in the forest ecosystem.

Author Contributions

Conceptualization, Z.Y.; Methodology, B.W. and Z.Y.; Software, X.J.; Validation, J.W.; Resources, J.W.; Writing—original draft, B.W.; Writing—review & editing, B.W., X.J., Z.Y., S.Y. (Shan Yu) and S.Y. (Shan Yin); Funding acquisition, B.W. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

No new data were created or analyzed in this study.

Conflicts of Interest

The authors declare no conflict 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 and Tables

Figure 1. Location of APENC and spatial distribution of ecological sites.

Figure 2. Relationship between observed and simulated values (a) SOS and (b) EOS.

Figure 3. Relationship between estimated and measured value in grassland-desert ecosystem (a,b), grassland ecosystem (c,d) and forest ecosystem (e,f).

Figure 3. Relationship between estimated and measured value in grassland-desert ecosystem (a,b), grassland ecosystem (c,d) and forest ecosystem (e,f).

Figure 4. Spatial pattern of average and standard deviation from 2001 to 2021 (a,b) SOS; (c,d) EOS and (e,f) LOS.

Figure 4. Spatial pattern of average and standard deviation from 2001 to 2021 (a,b) SOS; (c,d) EOS and (e,f) LOS.

Figure 5. Temporal trend of LSP (a) Yellow River Basin, (b) Haihe River Basin, (c) Songliao River Basin, (d) Continental Basin and (e) APENC.

Figure 5. Temporal trend of LSP (a) Yellow River Basin, (b) Haihe River Basin, (c) Songliao River Basin, (d) Continental Basin and (e) APENC.

Figure 6. Trend variation and significant of LSP from 2001 to 2021 (a,b) SOS; (c,d) EOS and (e,f) LOS.

Figure 6. Trend variation and significant of LSP from 2001 to 2021 (a,b) SOS; (c,d) EOS and (e,f) LOS.

Table 1

The basic information of four ecological sites in the APENC.

Scheme	Location	Ecosystem Type	Dominant Species
Inner Mongolia	43.55° N, 116.68° E	Grassland	Leymus chinensis, Stipa grandis, Agropyron cristatum and Cleistogenes squarrosa
Naiman	42.93° N, 120.70° E	Grassland and desert	Caragana microphylla, Artemisia halodendron and Setaria viridis
Ordos	39.49° N, 110.19° E	Grassland and desert	Artemisia ordosica, Hedysarum leave, Salix psammophila, Dysphania aristate, Leymus secalinus and Cleistogenes squarrosa
Beijing	39.96° N, 115.43° E	Forest	Pinus tabuliformis, Larix principis-rupprechtii and Quercus wutaishanica

Table 2

LSP variation in different ecosystems from 2001 to 2021 at regional and basin scales.

		APENC	YB	HB	SLB	CB	APENC	YB	HB	SLB	CB	APENC	YB	HB	SLB	CB
		SOS					EOS					LOS
Croplands	Slope	−0.053	−0.182 *	−0.160	0.059	0.165	0.065	0.196 **	0.149	−0.014	0.016	0.176	0.436 **	0.297 *	−0.097	−0.184
	Z value	−0.997	−1.962	−1.057	0.513	0.936	0.876	2.990	1.782	−0.091	0.151	1.842	3.835	1.963	−0.393	−0.815
	p value	0.319	0.047	0.291	0.608	0.349	0.381	0.003	0.075	0.928	0.880	0.066	0.000	0.049	0.695	0.415
Grasslands	Slope	−0.045	−0.250 *	−0.128	0.008	0.255	−0.008	0.069	0.101	−0.005	−0.160	0.106	0.388 **	0.145	−0.078	−0.177
	Z value	−0.453	−2.084	−0.634	0.211	1.299	−0.030	1.117	0.936	−0.030	−0.936	1.117	3.050	1.178	−0.091	−1.117
	p value	0.651	0.037	0.526	0.833	0.194	0.976	0.264	0.349	0.976	0.349	0.264	0.002	0.239	0.928	0.264
Forests	Slope	−0.665 **	−0.775 **	−0.578 **	−0.627 **	-	0.334 **	0.321 **	0.413 **	0.261 *	-	0.964 **	1.102 **	0.982 **	0.867 **	-
	Z value	−3.533	−4.439	−3.231	−2.808	-	3.895	3.110	3.775	2.386	-	4.318	5.164	4.801	3.352	-
	p value	0.000	0.000	0.001	0.005	-	0.000	0.002	0.000	0.017	-	0.000	0.000	0.000	0.001	-

Notes: YB, Yellow River Basin; HB, Haihe River Basin; SLB, Songliao River Basin; CB Continental Basin. The characters (‘*’ and ‘**’) indicate the significant (trends at p < 0.05 and p < 0.01, respectively). The sample sizes of cropland pixels are 4,755,492, 1,559,074, 570,867, 2,282,696 and 342,855 in APENC, YB, HB, SLB and CB, respectively; the grasslands are 6,261,773, 2,162,689, 504,429, 1,961,625 and 1,633,030, and the forests are 586,026, 116,741, 298,032, 156,765 and 14,488.

References

1. Hu, A. China’s goal of achieving carbon peak by 2030 and its main approaches. J. Beijing Univ. Technol.; 2021; 21, pp. 1-15.

2. IPCC. Summary for Policymakers of Climate Change 2021: The Physical Science Basis; Cambridge University Press: Cambridge, UK, 2021.

3. Cheng, W.; Li, Z.; Yan, L. Uniforming spring phenology under non-uniform climate warming across latitude in China. Sci. Total Environ.; 2021; 762, 143177. [DOI: https://dx.doi.org/10.1016/j.scitotenv.2020.143177]

4. Xie, Q.; Cleverly, J.; Moore, C.E.; Ding, Y.; Hall, C.C.; Ma, X.; Brown, L.A.; Wang, C.; Beringer, J.; Prober, S.M. et al. Land surface phenology retrievals for arid and semi-arid ecosystems. ISPRS J. Photogramm. Remote Sens.; 2022; 185, pp. 129-145. [DOI: https://dx.doi.org/10.1016/j.isprsjprs.2022.01.017]

5. Li, X.; Jiang, L.; Meng, F.; Wang, S.; Niu, H.; Iler, A.M.; Duan, J.; Zhang, Z.; Luo, C.; Cui, S. et al. Responses of sequential and hierarchical phenological events to warming and cooling in alpine meadows. Nat. Commun.; 2016; 7, 12489. [DOI: https://dx.doi.org/10.1038/ncomms12489] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/27535205]

6. Piao, S.; Ciais, P.; Friedlingstein, P.; Peylin, P.; Reichstein, M.; Luyssaert, S.; Margolis, H.; Fang, J.; Barr, A.; Chen, A. et al. Net carbon dioxide losses of northern ecosystems in response to autumn warming. Nature; 2008; 451, pp. 49-52. [DOI: https://dx.doi.org/10.1038/nature06444] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/18172494]

7. Khwarahm, N.R.; Dash, J.; Skjoth, C.A.; Newnham, R.M.; Adams-Groom, B.; Head, K.; Caulton, E.; Atkinson, P.M. Mapping the birch and grass pollen seasons in the uk using satellite sensor time-series. Sci. Total Environ.; 2017; 578, pp. 586-600. [DOI: https://dx.doi.org/10.1016/j.scitotenv.2016.11.004]

8. Chen, M.; Melaas, E.K.; Gray, J.M.; Friedl, M.A.; Richardson, A.D. A new seasonal-deciduous spring phenology submodel in the community land model 4.5: Impacts on carbon and water cycling under future climate scenarios. Glob. Chang. Biol.; 2016; 22, pp. 3675-3688. [DOI: https://dx.doi.org/10.1111/gcb.13326]

9. Qader, S.H.; Atkinson, P.M.; Dash, J. Spatiotemporal variation in the terrestrial vegetation phenology of iraq and its relation with elevation. Int. J. Appl. Earth Obs. Geoinf.; 2015; 41, pp. 107-117. [DOI: https://dx.doi.org/10.1016/j.jag.2015.04.021]

10. Bórnez, K.; Richardson, A.D.; Verger, A.; Descals, A.; Peñuelas, J. Evaluation of vegetation and proba-v phenology using phenocam and eddy covariance data. Remote Sens.; 2020; 12, 3077. [DOI: https://dx.doi.org/10.3390/rs12183077]

11. Ge, Q.; Dai, J.; Zheng, J.; Bai, J.; Zhong, S.; Wang, H.; Wang, W.-C. Advances in first bloom dates and increased occurrences of yearly second blooms in eastern china since the 1960s: Further phenological evidence of climate warming. Ecol. Res.; 2011; 26, pp. 713-723. [DOI: https://dx.doi.org/10.1007/s11284-011-0830-7]

12. Pertille, R.H.; Citadin, I.; Oliveira, L.d.S.d.; Broch, J.d.C.; Kvitschal, M.V.; Araujo, L. The influence of temperature on the phenology of apple trees grown in mild winter regions of Brazil, based on long-term records. Sci. Hortic.; 2022; 305, 111354. [DOI: https://dx.doi.org/10.1016/j.scienta.2022.111354]

13. Studer, S.; Stockli, R.; Appenzeller, C.; Vidale, P.L. A comparative study of satellite and ground-based phenology. Int. J. Biometeorol.; 2007; 51, pp. 405-414. [DOI: https://dx.doi.org/10.1007/s00484-006-0080-5] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/17235537]

14. Rodriguez-Galiano, V.F.; Dash, J.; Atkinson, P.M. Intercomparison of satellite sensor land surface phenology and ground phenology in europe. Geophys. Res. Lett.; 2015; 42, pp. 2253-2260. [DOI: https://dx.doi.org/10.1002/2015GL063586]

15. Alberton, B.; Torres, R.D.S.; Cancian, L.F.; Borges, B.D.; Almeida, J.; Mariano, G.C.; dos Santos, J.; Morellato, L.P.C. Introducing digital cameras to monitor plant phenology in the tropics: Applications for conservation. Perspect. Ecol. Conserv.; 2017; 15, pp. 82-90. [DOI: https://dx.doi.org/10.1016/j.pecon.2017.06.004]

16. Weil, G.; Lensky, I.M.; Levin, N. Using ground observations of a digital camera in the vis-nir range for quantifying the phenology of Mediterranean woody species. Int. J. Appl. Earth Obs. Geoinf.; 2017; 62, pp. 88-101. [DOI: https://dx.doi.org/10.1016/j.jag.2017.05.016]

17. Richardson, A.D.; Hufkens, K.; Li, X.; Ault, T.R. Testing hopkins’ bioclimatic law with phenocam data. Appl. Plant Sci.; 2019; 7, e01228. [DOI: https://dx.doi.org/10.1002/aps3.1228]

18. Zhou, J.; Cai, W.; Qin, Y.; Lai, L.; Guan, T.; Zhang, X.; Jiang, L.; Du, H.; Yang, D.; Cong, Z. et al. Alpine vegetation phenology dynamic over 16years and its covariation with climate in a semi-arid region of China. Sci. Total Environ.; 2016; 572, pp. 119-128. [DOI: https://dx.doi.org/10.1016/j.scitotenv.2016.07.206]

19. Bohovic, R.; Dobrovolny, P.; Klein, D. The spatial and temporal dynamics of remotely-sensed vegetation phenology in central asia in the 1982–2011 period. Eur. J. Remote Sens.; 2017; 49, pp. 279-299. [DOI: https://dx.doi.org/10.5721/EuJRS20164916]

20. 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]

21. Wu, W.; Sun, Y.; Xiao, K.; Xin, Q. Development of a global annual land surface phenology dataset for 1982–2018 from the avhrr data by implementing multiple phenology retrieving methods. Int. J. Appl. Earth Obs. Geoinf.; 2021; 103, 102487. [DOI: https://dx.doi.org/10.1016/j.jag.2021.102487]

22. Qader, S.H.; Priyatikanto, R.; Khwarahm, N.R.; Tatem, A.J.; Dash, J. Characterising the land surface phenology of middle eastern countries using moderate resolution landsat data. Remote. Sens.; 2022; 14, 2136. [DOI: https://dx.doi.org/10.3390/rs14092136]

23. 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]

24. Zhang, X.Z.; Friedla, M.A.; Schaaf, C.B.; Strahler, A.H.; HodgesJohn, C.F.; Feng, G.; Reed, B.C.; Huete, A. Monitoring vegetation phenology using modis. Remote Sens. Environ.; 2003; 84, pp. 471-475. [DOI: https://dx.doi.org/10.1016/S0034-4257(02)00135-9]

25. Tong, X.; Tian, F.; Brandt, M.; Liu, Y.; Fensholt, R. Trends of land surface phenology derived from passive microwave and optical remote sensing systems and associated drivers across the dry tropics 1992–2012. Remote Sens. Environ.; 2019; 232, 111307. [DOI: https://dx.doi.org/10.1016/j.rse.2019.111307]

26. Jeong, S.-J.; Schimel, D.; Frankenberg, C.; Drewry, D.T.; Fisher, J.B.; Verma, M.; Berry, J.A.; Lee, J.-E.; Joiner, J. Application of satellite solar-induced chlorophyll fluorescence to understanding large-scale variations in vegetation phenology and function over northern high latitude forests. Remote Sens. Environ.; 2017; 190, pp. 178-187. [DOI: https://dx.doi.org/10.1016/j.rse.2016.11.021]

27. Xin, Q.; Li, J.; Li, Z.; Li, Y.; Zhou, X. Evaluations and comparisons of rule-based and machine-learning-based methods to retrieve satellite-based vegetation phenology using modis and USA national phenology network data. Int. J. Appl. Earth Obs. Geoinf.; 2020; 93, 102189. [DOI: https://dx.doi.org/10.1016/j.jag.2020.102189]

28. Luo, Z.; Yu, S. Spatiotemporal variability of land surface phenology in China from 2001–2014. Remote Sens.; 2017; 9, 65. [DOI: https://dx.doi.org/10.3390/rs9010065]

29. Wang, L.; De Boeck, H.J.; Chen, L.; Song, C.; Chen, Z.; McNulty, S.; Zhang, Z. Urban warming increases the temperature sensitivity of spring vegetation phenology at 292 cities across China. Sci. Total Environ.; 2022; 834, 155154. [DOI: https://dx.doi.org/10.1016/j.scitotenv.2022.155154]

30. Piao, S.; Cui, M.; Chen, A.; Wang, X.; Ciais, P.; Liu, J.; Tang, Y. Altitude and temperature dependence of change in the spring vegetation green-up date from 1982 to 2006 in the Qinghai-Xizang plateau. Agric. For. Meteorol.; 2011; 151, pp. 1599-1608. [DOI: https://dx.doi.org/10.1016/j.agrformet.2011.06.016]

31. 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]

32. Zhao, S.Q. Ximeng-an economic geography investigation of farming-pastoral region of northern China. Acta Geogr. Sin.; 1953; 19, pp. 43-60.

33. Gao, J.; Lv, S.; Zheng, Z.; Liu, J. Typical ecotones in China. J. Resour. Ecol.; 2012; 3, pp. 297-307.

34. Fang, Z.H.; He, C.Y.; Liu, Z.F.; Zhao, Y.Y.; Yang, Y.J. Climate change and future trends in the agro-pastoral transitional zone in northern china: The comprehensive analysis with the historical observation and the model simulation. J. Nat. Resour.; 2020; 35, pp. 358-370.

35. Wei, B.; Xie, Y.; Jia, X.; Wang, X.; He, H.; Xue, X. Land use/land cover change and it’s impacts on diurnal temperature range over the agricultural pastoral ecotone of northern China. Land Degrad. Dev.; 2018; 29, pp. 3009-3020. [DOI: https://dx.doi.org/10.1002/ldr.3052]

36. Zhao, H.L.; Zhao, X.Y.; Zhang, T.H.; Zhou, R.L. Boundary line on agro-pasture zigzag zone in north china and its problems on eco-environment. Adv. Earth Sci.; 2002; 17, pp. 739-747.

37. Liu, J.H.; Gao, J.X. Changes of land use and landscape pattern in the boundary change areas in farming-pastoral ecotone of northern China. Trans. Chin. Soc. Agric. Eng.; 2008; 24, pp. 76-82.

38. Jakubauskas, M.E. Harmonic analysis of time-series avhrr ndvi data. Photogramm. Eng. Remote Sens.; 2001; 67, pp. 461-470.

39. Wei, B.; Xie, Y.; Wang, X.; Jiao, J.; He, S.; Bie, Q.; Jia, X.; Xue, X.; Duan, H. Land cover mapping based on time series modis-ndvi using a dynamic time warping approach: A case study of the agricultural pastoral ecotone of northern China. Land Degrad. Dev.; 2020; 31, pp. 1050-1068. [DOI: https://dx.doi.org/10.1002/ldr.3502]

40. Chen, J.; Chen, J.; Liao, A.; Cao, X.; Chen, L.; Chen, X.; He, C.; Han, G.; Peng, S.; Lu, M. et al. Global land cover mapping at 30m resolution: A pok-based operational approach. ISPRS J. Photogramm. Remote Sens.; 2015; 103, pp. 7-27. [DOI: https://dx.doi.org/10.1016/j.isprsjprs.2014.09.002]

41. CERN. Plant Phenological Observation Dataset of the Chinese Ecosystem Research Network (2003–2015); Science Data Bank CERN: Meyrin, Switzerland, 2017.

42. Hou, X.; Gao, S.; Niu, Z.; Xu, Z. Extracting grassland vegetation phenology in north china based on cumulative spot-vegetation ndvi data. Int. J. Remote Sens.; 2014; 35, pp. 3316-3330. [DOI: https://dx.doi.org/10.1080/01431161.2014.903437]

43. Bao, G.; Bao, Y.; Alateng, T.; Bao, Y.; Qin, Z.; Wang, M.; Zhou, Y. Spatio-temporal dynamics of vegetation phenology in the mongolia plateau duing 1982–2011. Remote Sens. Technol. Appl.; 2017; 32, pp. 866-874.

44. Zheng, W.; Liu, Y.; Yang, X.; Fan, W. Spatiotemporal variations of forest vegetation phenology and its response to climate change in northeast China. Remote Sens.; 2022; 14, 2909. [DOI: https://dx.doi.org/10.3390/rs14122909]

45. Gocic, M.; Trajkovic, S. Analysis of changes in meteorological variables using mann-kendall and sen’s slope estimator statistical tests in Serbia. Glob. Planet. Chang.; 2013; 100, pp. 172-182. [DOI: https://dx.doi.org/10.1016/j.gloplacha.2012.10.014]

46. Sen, P.K. Estimates of the regression coefficient based on kendall’s tau. J. Am. Stat. Assoc.; 1968; 63, pp. 1379-1389. [DOI: https://dx.doi.org/10.1080/01621459.1968.10480934]

47. Brenskelle, L.; Stucky, B.J.; Deck, J.; Walls, R.; Guralnick, R.P. Integrating herbarium specimen observations into global phenology data systems. Appl. Plant Sci.; 2019; 7, e01231. [DOI: https://dx.doi.org/10.1002/aps3.1231]

48. Liu, Z.; Liu, Y.; Li, Y. Anthropogenic contributions dominate trends of vegetation cover change over the farming-pastoral ecotone of northern China. Ecol. Indic.; 2018; 95, pp. 370-378. [DOI: https://dx.doi.org/10.1016/j.ecolind.2018.07.063]

49. Dong, M.Y.; Jiang, Y.; Ren, P.P.; Wu, Z.F. Variation trend and catastrophe change of air temperature in the farming pastoral ecotone of northern china during recent 50 years. J. Desert Res.; 2010; 30, pp. 926-932.

50. Zhao, W.; Wei, Z.; Zheng, Z.; Dong, W. Surface temperature and precipitation variation of pastoral transitional zone in northern china during 1964–2013. Plateau Meteorol.; 2016; 35, pp. 979-988.

51. Kewei, J.; Jiangbo, G.; Shaohong, W.; Wenjuan, H. Research progress on the response processes of vegetation activity to climate change. Acta Ecol. Sin.; 2018; 38, pp. 2229-2238.

52. Zhang, M.; Liu, H.; Wang, K.; Chen, Y.; Ren, Y.; Yue, Y.; Deng, Z.; Zhang, C. Nonlinear trends of vegetation changes in different geomorphologic zones and land use types of the yangtze river basin, China. Land Degrad. Dev.; 2023; 34, pp. 1-12. [DOI: https://dx.doi.org/10.1002/ldr.4627]

53. Duncan JM, A.; Dash, J.; Atkinson, P.M. Spatio-temporal dynamics in the phenology of croplands across the indo-gangetic plains. Adv. Space Res.; 2014; 54, pp. 710-725. [DOI: https://dx.doi.org/10.1016/j.asr.2014.05.003]

54. Sun, A.J.; Fan, J.H. Characteristics of the distribution of temperature variabiity in China. Acta Geogr. Sin.; 1985; 40, pp. 11-19.

55. Sa, C.; Meng, F.; Luo, M.; Li, C.; Wang, M.; Adiya, S.; Bao, Y. Spatiotemporal variation in snow cover and its effects on grassland phenology on the Mongolian plateau. J. Arid. Land; 2021; 13, pp. 332-349. [DOI: https://dx.doi.org/10.1007/s40333-021-0056-7]

56. Tong, S.; Li, X.; Zhang, J.; Bao, Y.; Bao, Y.; Na, L.; Si, A. Spatial and temporal variability in extreme temperature and precipitation events in inner Mongolia (China) during 1960–2017. Sci. Total Environ.; 2019; 649, pp. 75-89. [DOI: https://dx.doi.org/10.1016/j.scitotenv.2018.08.262] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/30172136]

57. Wei, B.; Bao, Y.; Yu, S.; Yin, S.; Zhang, Y. Analysis of land surface temperature variation based on modis data a case study of the agricultural pastural ecotone of northern China. Int. J. Appl. Earth Obs. Geoinf.; 2021; 100, 102342. [DOI: https://dx.doi.org/10.1016/j.jag.2021.102342]

58. Hong, Y.; Zhang, H.; Zhao, J.; Shan, Y.; Zhang, Z.; Guo, X.; Wu, R.; Deng, G. Effects of spring and summer extreme climate events on the autumn phenology of different vegetation types of inner Mongolia, China, from 1982 to 2015. Ecol. Indic.; 2020; 111, 105974.

59. Julien, Y.; Sobrino, J.A. Global land surface phenology trends from gimms database. Int. J. Remote Sens.; 2009; 30, pp. 3495-3513. [DOI: https://dx.doi.org/10.1080/01431160802562255]

60. Park, T.; Ganguly, S.; Tømmervik, H.; Euskirchen, E.S.; Høgda, K.-A.; Karlsen, S.R.; Brovkin, V.; Nemani, R.R.; Myneni, R.B. Changes in growing season duration and productivity of northern vegetation inferred from long-term remote sensing data. Environ. Res. Lett.; 2016; 11, 084001. [DOI: https://dx.doi.org/10.1088/1748-9326/11/8/084001]

61. Tao, J.; Zhang, G.; Wang, J.; Dong, J. Variation of cropland phenology in mid-eastern inner Mongolia. Sci. Agric. Sin.; 2011; 44, pp. 4583-4592.

62. Zhao, X.; Liu, J.; Yang, S.; Zhang, Q.; Gao, F.; Liu, Y. Spatiotemporal variations of typical woodland and grassland phenology and its response to meteoroligical factors in northern China. Acta Ecol. Sin.; 2023; 43, pp. 1-12.