Vegetation Phenological Changes in Multiple Landforms and Responses to Climate Change

Abstract

Vegetation phenology is highly sensitive to climate change, and the phenological responses of vegetation to climate factors vary over time and space. Research on the vegetation phenology in different climatic regimes will help clarify the key factors affecting vegetation changes. In this paper, based on a time-series reconstruction of Moderate-Resolution Imaging Spectroradiometer (MODIS) normalized difference vegetation index (NDVI) data using the Savitzky–Golay filtering method, the phenology parameters of vegetation were extracted, and the Spatio-temporal changes from 2001 to 2016 were analyzed. Moreover, the response characteristics of the vegetation phenology to climate changes, such as changes in temperature, precipitation, and sunshine hours, were discussed. The results showed that the responses of vegetation phenology to climatic factors varied within different climatic regimes and that the Spatio-temporal responses were primarily controlled by the local climatic and topographic conditions. The following were the three key findings. (1) The start of the growing season (SOS) has a regular variation with the latitude, and that in the north is later than that in the south. (2) In arid areas in the north, the SOS is mainly affected by the temperature, and the end of the growing season (EOS) is affected by precipitation, while in humid areas in the south, the SOS is mainly affected by precipitation, and the EOS is affected by the temperature. (3) Human activities play an important role in vegetation phenology changes. These findings would help predict and evaluate the stability of different ecosystems.

1. Introduction

Phenology is the study of seasonally recurring plant and animal life cycle stages or phenophases, such as the leafing and flowering of plants, the maturation of agricultural crops, the emergence of insects, and the migration of birds [1]. Phenology is considered the pulse of our planet (www.usanpn.org). The International Biological Program (IBP) defined phenology as a recurring event that our ancestors used to guide agricultural production for thousands of years [2]. Phenology can reflect the climate and environmental changes and is an important indicator of changes in the natural environment [3,4]. Moreover, plants are an important part of the ecosystem, providing the main portions of matter and energy [5]. Phenology plays a crucial role in the carbon balance of terrestrial ecosystems [6]. Therefore, it is very important to study vegetation phenology. A change in vegetation phenology has many consequences for ecological processes, agriculture, forestry, human health, and the global economy [7]. Elucidating the trends of vegetation phenology can improve our understanding of the impacts of climate change on ecosystem productivity and carbon cycling and the associated feedback to the climate [8].

In recent decades, warming increasing trends and frequent extreme events due to climate changes have produced significant impacts on many ecosystems and changed the vegetation phenology [9,10]. Vegetation phenology is very sensitive to climate change, and it is regarded as an important indicator of climate change [11,12,13,14,15]. Various biotic and abiotic factors influence vegetation phenology, including the temperature, photoperiod, precipitation/moisture, vegetation composition, and climate changes [16,17,18,19,20]. In a study conducted in 21 European countries, Menzel et al. found that leaf-on conditions arrived 2.5 days earlier and that leaf-off conditions were delayed by 1 day for every 1°C increase in the temperature, respectively [21]. In the northern hemisphere, spring temperatures were found to have the highest correlation with the start of the growing season (SOS) [22,23,24]. Evidence of the effects of precipitation, nutrients, and soil physical properties on the spring phenology is scarce, and where effects have been found, the effects are small relative to the effect of the temperature in England [25]. However, as research develops, the temperature has been found to not be the only limiting factor of vegetation phenology [26]. Precipitation is also an important factor affecting changes in vegetation phenology [27]. Precipitation during winter and spring has also been found to affect the spring phenology in a complex manner in middle and high latitudes in the Northern Hemisphere [28,29,30]. Solar radiation (the total insolation, peak irradiance, and photoperiod) has also been implicated as a major cue, affecting the phenology in both seasonal and aseasonal periods [31,32]. In some temperate tree species, the photoperiod and winter chilling requirements are also known to play a role in the spring phenology [11]. In addition, studies have suggested that elevated CO₂ levels can enhance leaf longevity and impact the vegetation phenology in autumn [33].

Climate changes may accelerate in the future, so there is an urgent need to understand the Spatio-temporal phenological driving forces behind these changes [34]. Advances in remote-sensing technology have greatly expanded the capability of phenology observation [5]. The greenness index has been widely used in vegetation phenology research around the world, and it is retrieved from various satellite remote-sensing data, such as the normalized difference vegetation index (NDVI) [35]. Piao et al. used the method of the greatest increase in the NDVI to extract vegetation phenological parameters in temperate regions of China [36]. The phenological changes in the northern hemisphere over the past 30 years have been discussed by the maximum ratio method [3]. The method of the maximum curvature change has been used to estimate the autumn phenology, which is suitable for the comparison of ground observation data with Moderate-Resolution Imaging Spectroradiometer (MODIS) NDVI data [37]. However, at present, most phenological studies are focused on plateaus, watersheds, forests, and other large-scale spaces [38,39,40,41,42]. There are few studies on the response mechanisms of the vegetation phenology to climate changes under different climatic conditions at medium and small scales. Furthermore, the phenological sensitivity to climate changes across all taxa and scales is different, and there is potential ecological asynchrony [43]. Due to differences in topography and geomorphology, unique regional microclimates are formed, and local models of vegetation phenology arise [37]. We expect that even in areas with similar geographical locations, different local topographical conditions would have a significant impact on the vegetation phenology. The main purpose of this study was to address the following questions:

(1)

What is the Spatio-temporal variation pattern of the vegetation phenology?

(2)

How do climatic factors, such as temperature, precipitation, and light, affect phenological changes?

Clarifying these issues would help predict and evaluate the stability of different ecological communities and provide an important scientific basis for local sustainable development.

2. Data and Methods

2.1. Study Area

Shaanxi Province (Figure 1) is located at the junction of the eastern humid area and western arid area in China, and its area is approximately 205,600 km². The climate differs from north to south, and the landforms are composed of plateaus, mountains, plains, and basins. According to the landform types, the province can be divided into six parts: the Mu Us Desert (MUD), Loess Plateau (LP), Guanzhong Plain (GZP), Qinling Mountains (QLM), Hanjiang Basin (HJB), and Daba Mountains (DBM) [44]. The only desert landform is a MUD, where the climate is dry and cold, and the vegetation types are mainly desert vegetation, such as Artemisia and Salix. In recent years, the ecological environment has changed due to the Grain for Green Program. The LP is the frailest ecological area, where Yellow River is divided into numerous small tributaries, and the vegetation types are mainly grassland and shrub. The GZP is the largest plain and grain-producing zone, and its vegetation types are mixed deciduous forests and crops. The highest-elevation zone is the QLM, which is predominantly forested, with abundant animal and plant species. The HJB is the source of the HanJiang River, where the vegetation types are broad-leaved forests and crops. There are many mountain ravines crisscrossing the DBM, where the district basically consists of evergreen broad-leaved forests and evergreen deciduous forests.

2.2. Data Source and Preprocessing

2.2.1. NDVI

This paper used MODIS normalized difference vegetation index (NDVI) data as the basic data for the production of remote-sensing vegetation phenology analyses. The NDVI with a 16-day temporal resolution and a 250-m spatial resolution (MOD13Q1, collection v006) was obtained from the Level-1 and Atmosphere Archive & Distribution System (LAADS) Distributed Active Archive Center (DAAC) website (https://ladsweb.modaps.eosdis.nasa.gov/), and two tiles (h26v05 and h27v05) covered the study area. The sixteen-year period between 2001 and 2016 was used. After we obtained the data, we extracted, mosaicked, reformatted, and reprojected the NDVI and pixel-reliability data layer using the MODIS projection tool (MRT) provided by the Land Processes Distributed Active Archive Center (LP DAAC). Subsequently, according to the methods in (SONG Chun-qiao 2011) [45], the pixel-reliability was used to process the NDVI data to improve reliability. In detail, for the same one-pixel position, NDVI values with pixel-reliability values of 0 and 1 were weighted with 1 and 0.8, respectively, and pixels with pixel-reliability data of 2 and 3 were weighted with 0.2. We performed this process in Matlab and marked the NDVI processed according to the pixel-reliability data as the pNDVI (Figure 2).

2.2.2. Meteorological Data

The meteorological data used in this study included the monthly mean temperature, monthly precipitation, and monthly sunshine hours. These data were calculated from the daily values, which were obtained from the Chinese ground meteorological data daily dataset (V3.0) downloaded from the China Meteorological Data website (http://cdc.cma.gov.cn/). The meteorological data were averaged monthly or summed after spatial interpolation using the inverse distance-weighted (IDW) method in ArcGIS, and Spatio-temporal meteorological data at the monthly time scale with a 250-m resolution were obtained covering 2001 to 2016.

2.2.3. Topographic Data

Because there are many landform types in the study area (Figure 1) and due to the altitudinal differences between different landform types, the relationship between the altitude and phenology was analyzed using digital elevation model (DEM) data with a 30-m spatial resolution to explore the influence of the altitude on phenological parameters. The DEM data were derived from the ASTER GDEM V2 dataset, which was downloaded from the Geospatial Data Cloud (http://www.gscloud.cn) of the Computer Network Information Center, Chinese Academy of Sciences. In ArcMap, the DEM data were obtained in separate tiles, which were subsequently mosaicked, resampled to 250 m, and then subset to Shaanxi Province (study area).

Land use data with a 30-m spatial resolution from the Resources and Environment Scientific Data Center, Chinese Academy of Sciences, were used to extract the NDVI in areas with vegetation coverage in the study area by resampling and masking in ArcMap. Pixels with NDVI<0.05 were eliminated to ensure that each NDVI pixel had vegetation cover.

2.3. Remote-Sensing Vegetation Phenology Extraction

2.3.1. Reconstruction of the Vegetation Index Time Series Curve

The time series curve of the vegetation index (Figure 2) reflects the growth of vegetation and is the basis for extracting the parameters necessary to assess the vegetation phenology via remote sensing. The NDVI values from the same year were layer-stacked in ArcGIS to obtain multi-band images of 23 periods in that year, with each band representing the NDVI for a certain period, and then the multi-band images were read as a multi-dimensional matrix in Matlab, pixel by pixel. Subsequently, the multi-dimensional matrix was converted into a 1-dimensional sequence to obtain the NDVI time series of the pixel, as shown in Figure 2 NDVI. However, the original NDVI time series curve (Figure 2, NDVI) is uneven and abrupt due to interference from poor atmospheric conditions, clouds, and other unfavorable factors [46], and it is not suitable for use in the direct extraction of vegetation phenological parameters. Therefore, the Savitzky–Golay (S-G) filtering method was used to reconstruct the time series curve before extracting the vegetation phenological parameters.

The S-G filtering method was proposed by Savitzky and Golay in 1964 and is a weighted fitting algorithm for moving windows. The weight coefficients are derived by the least-squares fitting of a given high-order polynomial within the moving window [47]. Therefore, during the process of smoothing and denoising the NDVI time-series data with the S-G filter, the establishment of two parameters is crucial. The first parameter is m, which is the half-width of the smoothing window and is set to range from 3 to 7 [47]. Usually, a large m value produces a smooth result at the expense of flattening sharp peaks. The second parameter is an integer (d) that specifies the degree of smoothing of the polynomial, which is typically set to range from 2 to 4 [47]. A small d value will produce a smooth result but may introduce bias, while a high d value will reduce the filter bias but may “overfit” the data and yield a noisy result. Therefore, after comparing the results of different (m, d) combinations of parameters, the moving window (m) value was set to 5, the polynomial order (d) was set to 2, and the S-G filter was employed to reconstruct the NDVI time series (Figure 2, sgNDVI).

2.3.2. Extraction of Phenological Parameters

The extraction of vegetation phenology assumes that there is only one growing season in a year (for areas with two or more growing seasons in a year, this paper extracted the growing season with the largest NDVI value during the vegetation growing seasons in the area). Three commonly used metrics for the assessment of vegetation growth, namely, the SOS, end of the growing season (EOS), and length of the growing season (LOS), were extracted. The SOS refers to the date of onset of photosynthetic activity, and the EOS refers to the date at which the photosynthetic activity and green leaf area begin to rapidly decrease [48]. The LOS is calculated as the difference between the EOS and SOS.

Field-based ecological studies have demonstrated that after leaf emergence, a plant will grow rapidly for a period of time until it reaches a growth peak at a relatively stable stage and will then enter a similar but opposite pattern, namely, the rapid decline [48]. When this pattern is reflected in the growth curve of plants, the growth state of plants can be determined by observing the shape of the curve. There are many methods for extracting vegetation phenology based on remote-sensing vegetation index time series. The commonly used methods are the threshold method, the maximum slope method, the maximum ratio method, and the cumulative frequency method [36,49,50,51]. Considering the computational efficiency and operability characteristics of the methods, the maximum ratio method was used to extract the vegetation phenology in the study area. The maximum ratio describes the maximum change rate of the vegetation index, which is calculated as [3]:

$${\mathrm{NDVI}}_{\mathrm{ratio}}\left(i\right)=\frac{{\mathrm{NDVI}}_{i+1}-{\mathrm{NDVI}}_i}{{\mathrm{NDVI}}_i}$$

(1)

where i is the time sequence number. We determined the i value with the maximum NDVI_ratio and then used the corresponding Julian day as the SOS. Likewise, we determined the i value with the minimum NDVI_ratio and used the corresponding Julian day at i + 1 as the EOS. We extracted the phenological parameters of the vegetation in the study area pixel by pixel and saved the results as a TIF-format image in Matlab. Finally, we made a thematic map in ArcGIS, as shown in Figure 3.

2.4. Trend and Correlation Analysis

In the investigation of the temporal and spatial changes in phenology over the study area, we performed a linear regression of the phenological parameters against the SOS, EOS, and LOS and used the resulting regression coefficients to analyze the spatial-temporal trends of the phenology. The mathematical formula can be expressed as follows:

$$\mathrm{trend}=\frac{n\times {{\displaystyle \sum }}_{t=1}^{n}t\times {\mathrm{y}}_t-{{\displaystyle \sum }}_{t=1}^{n}t{{\displaystyle \sum }}_{t=1}^n{\mathrm{y}}_t}{n\times {{\displaystyle \sum }}_{t=1}^n{t}^2-{\left({{\displaystyle \sum }}_{t=1}^{n}t\right)}^2}$$

(2)

where n is the years during the monitoring period, t is a certain year, and y_t is the vegetation phenological parameter value for year t. trend > 0 indicates a delayed trend, while trend < 0 indicates an advancement trend.

In order to explore the response of the phenology to the temperature, precipitation, and sunshine hours, we used Pearson’s correlation analysis. Specifically, (1) the correlations of SOS with the temperature and precipitation from February to April were evaluated; (2) the correlations between the EOS and the temperature and precipitation from September to November and the correlation between the EOS and the sunshine hours from May to July were evaluated; and (3) the correlations between the LOS and the temperature and precipitation from April to September and the correlation between the LOS and the monthly sunshine hours from May to July were evaluated. The R-value was calculated as follows:

$$\mathrm{r}=\frac{{{\displaystyle \sum }}_{t=1}^n\left({\mathrm{x}}_t-\overline{\mathrm{x}}\right)\left({\mathrm{y}}_t-\overline{\mathrm{y}}\right)}{\sqrt{{{\displaystyle \sum }}_{t=1}^n{\left({\mathrm{x}}_t-\overline{\mathrm{x}}\right)}^2}\sqrt{{{\displaystyle \sum }}_{t=1}^n{\left({\mathrm{y}}_t-\overline{\mathrm{y}}\right)}^2}}$$

(3)

where x_t is the meteorological data for year t. y_t is the phenological parameter value for year t. We performed Pearson’s correlation analysis and saved the results as a TIF-format image in Matlab, which were used to make a thematic map in ArcGIS.

3. Results

3.1. Spatial Distribution Patterns of the Multi-Annual Average Vegetation Phenology

Figure 3 shows the spatial distributions of the SOS, EOS, and LOS in different geomorphic regions of the study area during 2001–2016. As shown in the figure, the vegetation phenology throughout the study area varied with the latitude from north to south, but due to the differences in topography and climate in each geomorphic zone, the vegetation phenology exhibited spatial heterogeneity. Compared with the southern area, the northern area had a later SOS and EOS and a shorter LOS. The difference in SOS between the north and south was the most obvious, with the SOS occurring early in the south and late in the north.

Table 1. Multi-annual mean vegetation phenology in different regions (day).

	SOS		EOS		LOS
	mean	std	mean	std	mean	std
Shaanxi	120.68	28.36	280.04	18.47	160.36	29.33
Mu Us Desert	132.58	36.80	291.07	12.59	158.49	46.64
Loess Plateau	134.16	27.72	289.65	11.97	155.50	28.21
Guanzhong Plain	101.02	27.31	274.90	25.70	173.88	26.66
Qinling Mountain	113.43	16.87	270.24	16.27	156.81	25.48
Hanjiang Basin	88.10	15.21	271.81	15.29	183.72	23.14
Daba Mountain	109.51	17.96	267.08	17.13	157.57	25.42

SOS: start of the growing season; EOS: end of the growing season; LOS: length of the growing season.

presents a statistical analysis of the spatial distribution of the vegetation phenology illustrated in Figure 3. In Shaanxi, the annual average SOS was 120 d, the annual average EOS was 280 d, and the annual average LOS was 160 d. In the northern LP, the annual average SOS was 134 d, which was greater than that in other geomorphic regions. The SOS first appeared in the HJB, where the climate is humid all year, and the altitude is relatively low, with an annual average of 88 d. The EOS was the latest in the MUD, with an annual average of 291 d, and the earliest in the DBM, at 267 d. The HJB had the longest LOS of 183 d. In the other geomorphic regions, the LOS exhibited little difference, ranging from 155 d to 173 d.

There was a close relationship between the vegetation phenology and altitude in Shaanxi. Pearson’s correlation analysis indicated that the relationship was significant. The characteristics of the changes in vegetation phenology with altitude are shown in Figure 3. The SOS was significantly delayed with increasing altitudes, and the extent of the delay was approximately 1.4 d/100 m (r = 0.92, P < 0.01). The LOS was significantly shortened at a rate of approximately 2 d/100 m (r = −0.97, P < 0.01). The EOS was both postponed and advanced with increasing altitudes, but it showed a significant advancing trend overall, with the rate of advancement of approximately 0.6 d/100 m (r = −0.80, P < 0.01).

3.2. Spatial Distribution Patterns of Interannual Phenological Changes

As shown as Figure 4, throughout the study area, over 66.39% of the pixels showed that the SOS had an advancing trend (0.89 d/years), and the trends of 6.89% of the pixels were significant. In the LP, the percentage of pixels with advancing trends in this area was 78.39%, and the trends were significant in 11.47% of the pixels, with a rate of advancement of 1.51 d/year. However, in the northern MUD, delayed-SOS pixels accounted for the majority, covering 71.97% of the pixels. In some geomorphic regions, the advancement of the SOS gradually weakened from north to south—for example, in the LP, the trend in the north was larger than that in the south. There were no significant differences in the EOS in terms of advancements and delays, and 53.81% of the pixels exhibited delays, while 2.48% of the pixels exhibited significant delays, with the average rate of delay of 0.18 d/year. These pixels were mainly distributed in the QLM, HJB, and LP. The highest rate of delay of 64.54% occurred in the QLM, with a rate of delay of 0.83 d/year. The areas with advancements were mainly distributed in the MUD, DBM, and GZP, with a rate of advancement of 0.34 d/year. The LOS showed an overall trend of extension, covering 65.05% of the pixels, of which 5.57% showed significant extension, with the rate of extension of 1.07 d/year, which was greater than the changes in SOS and EOS. The prolongation trends of the area were mostly in the LP and QLM, where 73.03% and 72.32% of the pixels exhibited delays. In the LP, 7.32% of the pixels exhibited significant prolongation, and the rate of prolongation in the QLM was significantly higher than the average, reaching 1.88 d/year. The shortening of the LOS was mainly concentrated in the MUD, where the shortening ratio reached 72.03% and significantly decreased to 9.76%, and the shortening time was 2.52 d/year.

3.3. The Response of the Vegetation Phenology to Climatic Factors

3.3.1. The Response of the SOS to Precipitation and Temperature

We can see from Figure 5, for the entire study area, the average precipitation from February to April had a negative correlation with the SOS, with negative correlations occurring in 61.88% of the pixels, indicating that with the increase of precipitation, the SOS advanced. The greatest influence of precipitation on the SOS occurred in March when 7.5% of the pixels exhibited significant negative correlations. The average temperature from February to April also had a negative correlation with the SOS in 58.14% of the pixels, where the SOS advanced with the increase of the temperature. However, the greatest influence of the temperature on the SOS occurred in February, with the negative correlation being significant in only 3.17% of the pixels.

The influences of the precipitation and temperature on the SOS in different geomorphic regions of the study area differed from north to south. The MUD and LP in the north were mainly affected by the temperature. In February, the temperature had the greatest influence in the MUD and LP, where the proportion of significant pixels reached 14.42% and 11.83%, respectively, which mainly exhibited positive correlations. Precipitation had a mainly negative correlation with the SOS in the MUD and LP, but fewer pixels were affected by precipitation than by the temperature. The GZP and QLM in the central region were affected by temperature and precipitation. Precipitation and temperature had mainly negative correlations in the QLM. In the GZP, precipitation and temperature had mainly positive correlations in February and March and negative correlations in April. The HJB and DBM in the south were strongly affected by precipitation. Precipitation had the greatest influence on the SOS in March in the DBM, with significant pixels accounting for 12.16% of the area, in which the correlations were positive. The delays in these regions might be due to the increase in precipitation, which might have resulted in a decrease in the temperature and a delay of the SOS. However, the response to temperature was weak in the DBM, and the largest proportion of significant pixels in March was only 2.74%, which was far less than those in other areas.

3.3.2. The Response of the EOS to Precipitation, Temperature, and Sunshine Hours

As shown as Figure 6, compared with the SOS, the impacts of the Spatio-temporal variations in the climate on the EOS were more complex. There were limited positive and negative correlations between the EOS and precipitation in September and October, and there were few significant pixels. Precipitation had the greatest influence on the EOS in November. positive correlations covered 64.88% of the area, and the influence was significant in 10.73% of the pixels, indicating that the EOS was delayed in most regions with the increase of precipitation. There were limited positive and negative correlations between the EOS and temperature from September to November. The greatest influence of the temperature on the EOS occurred in September, and 10.13% of the pixels exhibited significant positive correlations. The influence of the sunshine hours on the EOS was positive from May to July. This result indicated that the EOS was delayed as the sunshine hours increased. In May and July, there were no significant differences in the effects of the sunshine hours on the EOS, and the influences were significant in 8.75% and 9.63% of the pixels, respectively, compared with 4.58% in June.

Similar to the SOS, the response of the EOS to the climate also showed differences between the north and south. However, in contrast to the SOS, the temperature had a mainly positive correlation in the central and southern areas of the study area to the south of the GZP. The temperature in September had the greatest influence in the GZP and QLM, at 15.97%, while 12.61% of the pixels exhibited significant positive correlations. In contrast, in the DBM and HJB, the temperature had the greatest influence in October, when the proportions of pixels with significant positive correlations reached 12.97% and 10.82%, respectively. Precipitation mainly affected the MUD and LP in the north of the study area. The greatest impact from precipitation on the EOS in the MUD and LP occurred in November, when significant positive correlations occurred in 17.85% and 16.11% of the pixels, respectively. This finding indicated that the EOS was delayed with the increase of precipitation. In summer, most of the plants are luxuriant. The sunshine hours affect changes in the EOS and LOS by affecting photosynthesis [52]. The greatest effects of the sunshine hours on the EOS occurred in the MUD and LP, and the correlations were mainly positive. This finding indicated that the EOS was delayed as the sunshine intensified. The highest proportion of significant pixels was 16.61%, which occurred in the MUD in May, while 13.61% were significant in the LP in July. In other areas, the effects of the sunshine hours on the EOS were mainly positive, but the number of significant pixels was low.

3.3.3. The Response of the LOS to the Precipitation, Temperature, and Sunshine

Figure 7 shows the direction of the correlation coefficient between the precipitation and LOS changed from April to September. Positive correlations covered 70.12% of the area from April to June, and negative correlations covered 74.18% of the area in July and August. In September, about 57.6% of the area showed a positive correlation between the LOS and precipitation. The greatest influence of precipitation on the LOS reached 12.62% in July, which was a significant negative correlation. This result showed that the LOS in most areas would be shortened with the increase of precipitation in July. The proportions of positive and negative correlations of the LOS with the temperature from April to September were 45.03% and 54.97%, respectively. The greatest influence of the temperature on the LOS reached 72.41% in May when the correlation was negative, whereas 12.03% of the pixels exhibited positive correlations. The influence of the sunshine hours on the LOS was less significant than the influences of precipitation and temperature. The effect of the sunshine hours on the LOS was the most significant in May when the proportion of significant pixels was 6.85%. There were positive correlations in May and July and negative correlations in June.

In contrast to the SOS and EOS, the influence of the climate on the LOS exhibited Spatio-temporal heterogeneity and no regularity. Precipitation had the greatest effect on the LOS in most regions in July, and the correlations were mainly negative. The largest correlation occurred in the HJB, where 17.41% of the pixels showed significant correlations. In the LP, the greatest effects of precipitation on the LOS occurred in August, when 10.19% of the pixels showed significant correlations. In April, the temperature had the greatest effect in the DBM, and the significant pixel ratio was 8.55%, with most being positively correlated. In May, the temperature had the greatest effect in the LP and QLM, and the proportions of significant pixels were 13.66% and 14.23%, respectively, which mainly exhibited negative correlations. The GZP was most affected by the temperature in September when the proportion of significant pixels was 8.98%. The LOS in the MUD was not affected by the temperature, and the average ratio of significant pixels was only 2.29%, which was much lower than those in other areas. The influence of the sunshine hours on each geomorphic area was positive, and the greatest effect was observed in May. The highest proportion of significant pixels was 9.93% in the DBM. The influence of the sunshine hours in the other regions was not significantly different, and the proportion of significant pixels was less than 6% from May to July.

4. Discussion

The extent to which climatic factors influence the vegetation phenology is dependent on the different developmental stages of the phenological events [5]. Moreover, the responses of the vegetation phenology to environmental changes are non-uniform in different regions [52], in which the climate and vegetation type can differ greatly. Under normal conditions, the northern arid area experiences less precipitation and lower temperatures than does the south in autumn, so the EOS should be earlier. However, interestingly, the EOS in the northern area was found to be significantly later than that in the southern area in our study. This difference might arise because the vegetation in the north is mostly Artemisia and herbaceous plants, for which the EOS occurs around the end of October or the middle of November, which is later than the EOS for corn, wheat, and other crops in the south. This finding showed that the responses of the EOS to climate changes were more complex than those of the SOS, and the two trends are asymmetric [3,11].

Vegetation phenology was influenced by human activities. Interannual variations in the phenology indicated that the SOS was undergoing an advancement trend, while the LOS was being prolonged. However, in the two adjacent climatic regions in the north, the results were the opposite. In the LP, the SOS underwent the most significant advancement, while the LOS was extended. For the MUD, the SOS was delayed, and the LOS was shortened. This trend occurred mainly because the main vegetation comprises shrubs and grassland in the LP, which are more sensitive to temperature changes than the other types of vegetation. Because trees have been planted to control desertification in the MUD over the past decade, the ecological environment has improved, and the vegetation coverage has increased from 9.9% in 2000 to 38.03% in 2018 (http://www.shaanxi.gov.cn). The vegetation type also changed from Salix and other wind-break plants to crops, such as potatoes and corn (sprouting in May and harvested in October), which have a later SOS and a shorter LOS than those of forest and grass plants.

Most studies have found temperature to be the main factor that controls the vegetation phenology [53,54]. In this study, precipitation was shown to play a role as critical as that of the temperature. For the SOS, most areas of the MUD and LP were mainly affected by the temperature in the north. The central GZP and QLM were affected by both the temperature and precipitation, and the precipitation in the HJB and DBM was the key factor in the south. However, for the EOS, precipitation mainly affected most of the northern region, while the central and southern regions were mainly affected by the temperature. These results showed that the effects of climate changes on the vegetation phenology were complex. These impacts differed by region, and the climate factors that control the vegetation phenology were different in different geomorphic areas. Similarly, we found that the sunshine hours were positively correlated with the EOS from May to July. This finding was consistent with the results of previous studies [55]. When summer transitions into autumn, the day gradually shortens (the sunshine hours decrease), which induces bud setting and leaf senescence [56].

Vegetation phenology was affected by different climate factors. In the northern MUD, the SOS was most significantly affected by the temperature. In contrast, in the southern DBM, precipitation had the largest influence on the SOS, whereas the temperature had the least influence. This phenomenon might be due to the lower temperature in spring in the north, and the increase of the temperature was key to breaking the dormancy of plants. For the warm south, the effect of precipitation was more important than that of temperature for the growth of plants. Furthermore, in different geomorphic areas, the response relationships and correlation coefficients were different between the vegetation phenology and climate changes. In February, the temperatures in the LP and MUD were positively correlated with the SOS, while in the HJB and DBM, they were negatively correlated. The precipitation in October was negatively correlated with the EOS in the QLM, HJB, DBM, and GZP, while in the LP and MUD, they were positively correlated in November. In some northern arid areas, the increase of the temperature in spring might lead to the decrease of the soil moisture and have a negative impact on plant growth, while in southern humid areas, the increase of precipitation in autumn would lead to the decrease of the temperature, and plant growth would stop early.

5. Conclusions

In the context of global warming, quantifying the responses of phenology to climate drivers is essential for the sustainable development of ecosystems. This study analyzed the spatial differentiation of the phenology and the spatial-temporal variation of the phenology in six different geomorphological landscapes in Shaanxi, China, and discussed the responses of the phenology to climate changes. The main conclusions are as follows:

(1)

The phenology varied with latitude. The SOS in the south was earlier, whereas that in the north was later.

(2)

In the whole region, the SOS was gradually being advanced, while the LOS was being prolonged, but in the MUD, the SOS was being delayed, and the LOS was being shortened. The variation in the phenology had spatial heterogeneity among different geomorphological areas.

(3)

The responses of the phenology to climatic factors were very complex. In the north, the SOS was mainly affected by the temperature in February, and the EOS was mainly affected by the precipitation in November. In the south, the SOS was mainly affected by the precipitation in March, and the EOS was affected by the temperature in October.

Our results showed that the phenology was very sensitive to the climate, and there might be different types of phenology due to different topography and vegetation types. These findings might help reveal local phenological patterns related to microclimates and in understanding phenology at different scales.