The Response of Vegetation Phenology and Productivity to Drought in Semi-Arid Regions of Northern China

A major disturbance in nature, drought, has a significant impact on the vulnerability and resilience of semi-arid ecosystems by shifting phenology and productivity. However, due to the various disturbance mechanisms, phenology and primary productivity have remained largely ambiguous until now. This paper evaluated the spatio-temporal changes of phenology and productivity based on GIMMS NDVI3g time series data, and demonstrated the responses of vegetation phenology and productivity to drought disturbances with the standardized precipitation evapotranspiration index (SPEI) in semi-arid ecosystems of northern China. The results showed that (1): vegetation phenology exhibited dramatic spatial heterogeneity with different rates, mostly presented in the regions with high chances of land cover type variation. The delayed onset of growing season (SOS) and advanced end of growing season (EOS) occurred in Horqin Sandy Land and the eastern Ordos Plateau with a one to three days/decade (p < 0.05) rate and in the middle and east of Inner Mongolia with a two days/decade rate, respectively. Vegetation productivity presented a clear pattern: south increased and north decreased. (2) Spring drought delayed SOS in grassland, barren/sparsely vegetated land, and cropland, while autumn drought significantly advanced EOS in grassland and barren/sparsely vegetated lands. Annual drought reduced vegetation productivity and the sensitivity of productivity regarding drought disturbance was higher than that of phenology.


Introduction
Drought has affected most regions of East Asia in the early 21st century [1,2].The frequency of extreme climatic events will greatly increase for the next 40 years in Northern China, most notably in semi-arid regions [3,4].There will probably be an enhancement of the drought sensitivity of ecosystems [5], reduction in productivity and biodiversity, and change of vegetation species composition [6][7][8][9].Furthermore, ecosystem stability has been dramatically disturbed by drought [10].Therefore, it is vital to clearly understand the disturbance mechanisms of drought on ecosystems in order to cope with extreme drought events in the future.
Carbon and hydrology processes of terrestrial ecosystems are the main pathways disturbed by drought [11], particularly in a water-limited semi-arid ecosystem.Generally, moderate drought reduces the photosynthesis rate via stomatal control and improves water use efficiency with physiological adaption of vegetation at an individual level [7,11].At an ecosystem level, productivity and biodiversity significantly decrease with reduced photosynthesis and evapotranspiration rates [6,11,12].However, due to different response mechanisms to drought among vegetation species, the impacts distributed in these regions, which intensifies the vulnerability of the semi-arid ecosystem.Grassland, cropland, and barren or sparsely vegetated land are the main land cover types in the relatively stable area that occupy 69.49%, 0.97%, and 7.26% of the total region, respectively.In addition, a few shrubs and forest areas are scattered across the whole area (Figure 1b) based on 2001-2012 MODIS land cover data (MCD12Q1).A mean annual precipitation increase from less than 200 mm in the west arid grassland to more than 500 mm in the east cropland was recorded from 1982-2012 (Figure 1c).Mean annual temperature was found to have increased from −8.36 °C to 13.83 °C from north to south (Figure 1d) by analyzing Chinese Meteorological Data in this study.

GIMMS NDVI 3g Dataset
The Global Inventory Modeling and Mapping Studies (GIMMS) 15-day composite NDVI3g dataset with an 8 km spatial resolution applied here has been shown to be more accurate than the GIMMS NDVI for monitoring vegetation activity and phenological change [33].Data were derived from the AVHRR instrument on board the NOAA satellite series (7,9,11,14,(16)(17)(18)(19) for the time period of July 1981 to December 2013.These data had been corrected for calibration, solar geometry, heavy aerosols, clouds, and other effects not related to vegetation change with flag information [34].If the flag shows 1, 2, or 7, it represents a good value or missing data, respectively.Other flag values mean that NDVI was retrieved from different methods [35].Although GIMMS NDVI3g data has limitations, it is the longest time series of NDVI for monitoring the characterization and variability of vegetation.The dataset was directly downloaded from the NASA website [35].For each pixel, the GIMMS NDVI3g time series dataset was pre-processed using the Savitzky-Golay Filter to smooth NDVI with erroneous spikes due to clouds [36].

GIMMS NDVI 3g Dataset
The Global Inventory Modeling and Mapping Studies (GIMMS) 15-day composite NDVI3g dataset with an 8 km spatial resolution applied here has been shown to be more accurate than the GIMMS NDVI for monitoring vegetation activity and phenological change [33].Data were derived from the AVHRR instrument on board the NOAA satellite series (7,9,11,14,(16)(17)(18)(19) for the time period of July 1981 to December 2013.These data had been corrected for calibration, solar geometry, heavy aerosols, clouds, and other effects not related to vegetation change with flag information [34].If the flag shows 1, 2, or 7, it represents a good value or missing data, respectively.Other flag values mean that NDVI was retrieved from different methods [35].Although GIMMS NDVI3g data has limitations, it is the longest time series of NDVI for monitoring the characterization and variability of vegetation.The dataset was directly downloaded from the NASA website [35].For each pixel, the GIMMS NDVI3g time series dataset was pre-processed using the Savitzky-Golay Filter to smooth NDVI with erroneous spikes due to clouds [36].

Drought Index and Meteorological Data
As a site-specific drought indicator, SPEI is particularly suited to detecting, monitoring, and exploring the consequences of global warming on drought conditions [37].It is formulated based on precipitation and potential evapotranspiration (PET) [31,32].Different SPEIs are obtained for different time-scales representing the cumulative water balance over the previous several months [15].To calculate SPEI at different scales, we collected mean monthly air temperature and total monthly precipitation data during the period 1982-2012 from the China Meteorological Data Sharing Service System [38].Some missing data in the specific stations were replaced by linear regression estimation according to the closest series.
Here, three-month scale and 12-month scale SPEI values were calculated at 116 weather stations via the SPEI Calculator (SPEI calc) [39,40]).We chose the three-month scale SPEI from May as the spring drought index (SPEI_May) and from October as the autumn drought index (SPEI_Oct), and the 12-month scale SPEI value from December as the annual drought index (SPEI_year).Based on the categorization of dryness/wetness grade by SPEI [41], drought is defined as when SPEI was less than −1.0, and wet was defined as greater than 1.0.The normal period was defined as −1.0 < SPEI < 1.0.
All meteorological monthly data were interpolated 8 km spatial resolution grid data based on Partial Thin Plate Smoothing Splines (PTPSS) Models via ANUSPLIN software [42,43].The interpolation method accuracy of ANUSPLIN is higher than that based on the inverse distance weight method and ordinary Kriging method in China [44].According to the interpolation, the average signal to noise ratio (SNR) of precipitation, temperature, and SPEI is 0.340, 0.372, and 0.371, respectively, which is smaller than 1.The mean square for error (MSE) of the meteorological elements (precipitation, 0.0012; temperature, 0.0018; SPEI, 0.0065) is less than the generalized cross validation (GCV) value (precipitation, 0.0052; temperature, 0.0055; SPEI, 0.0336).These assessment indicators showed that the PTPSS Models and simulated surface are good for interpolating meteorological elements.

Other Auxiliary Data
MODIS IGBP land cover data (MCD12Q1) with a 500m spatial resolution and 75% accuracy [45] from 2001 to 2012 were analyzed to trace land cover changes and a relatively stable area (Figure S1), where change times are non-zero and zero values, respectively.Then, the stable areas were reclassified into five main classes, including forest, grassland, cropland, shrub, and barren/sparsely vegetated lands (Figure 1b).We extracted the percentage for the five land cover types at 8 km.The largest percentage type was defined as a new type for the 8-km pixel.Additionally, elevation data were derived from the "China Western Environment and Ecology Science Data Center" [46] with a spatial resolution of 1 km, and were resampled into 8 km to interpolate meteorological data via ANUSPLIN.

Extraction of Phenological Metrics and Vegetation Productivity
Here, the dynamic threshold method and derivative method were selected to extract phenological metrics (e.g., SOS, EOS, LOS) with smoothed GIMMS NDVI3g data, based on previous research and characteristics of the study area [47][48][49][50][51].

Dynamic Threshold Method (m1)
Taking 50% of the maximum NDVI value as the threshold to extract SOS and EOS (Figure 2a) could be reliable for different land cover types in semi-arid regions of China referencing previous studies [48].LOS was calculated as the difference between SOS and EOS.Finally, the integral value of smoothed NDVI during the growing season (iNDVI) was taken as a proxy of NPP (light blue part in Figure 2a).

Derivative Method Based on the Fitted Double Logistic Function (m2)
The theory of the derivative method is based on the NDVI time series curve and its derivation function to retrieve SOS and EOS.The date of maximum and minimum derivation is SOS and EOS, respectively (Figure 2b) [51].The Double Logistic function was used for fitting the NDVI time series curve.The fitted equitation is as follows [52,53]: where NDVI min and NDVI max are the minimum and the maximum NDVI during the growing year season, respectively; S and A are two inflection points, one as the curve rises (S) and one as it drops (A); and r i and r m are the rate of increase or decrease of the curve at the inflection points, respectively.

Trend Analysis
Temporal trends in the datasets were examined by applying the Theil-Sen (TS) median slope trend analysis, which is a robust trend detection method.Functionally similar to linear least squares regression, it operates on non-parametric statistics and is independent of the assumptions of linear regression.Trends are estimated using the median values and are therefore less susceptible to noise and outliers [54,55].
where, x i and x j are the value of the ith or jth year, respectively.A non-parametric test (Mann-Kendall significance test) was used to test the significance of the TS median slope value, producing z scores and providing information on both the significance and direction of the trend.A positive slope (z ≥ 1.96) represents a significant increase and a negative slope (z ≤ −1.96) indicates a significant decrease (p < 0.05) over time [56].

Correlation Analysis
The anomalies of SOS, EOS, LOS, and iNDVI were calculated to assess the response of phenology and productivity to drought by correlation analysis.Standardized anomalies were calculated by dividing anomalies by the standard deviation: where, x sd is the standard anomalies of SOS, EOS, LOS, and iNDVI; x is the SOS, EOS, LOS, and iNDVI at any given year; and µ and σ are the mean and standard deviation of SOS, EOS, LOS, and iNDVI during the 1982-2012 time period, respectively.The slope, p-value and correlation coefficients of the linear regression between anomalies and SPEI were analyzed.All the data processing was finished in IDL 8.3 and ArcGIS 10.2.

Climatic Variation Analysis in Semi-Arid Regions of Northern China
The climatic variation or dry-wet conditions in the study area were characterized by precipitation, temperature, and SPEI values.Four main land cover types, namely: grassland, cropland, barren/sparsely vegetated lands, and shrub, showed a significant warming trend from 1982 to 2012 (Figure S2).Although precipitation has seen no significant trend during the last three decades, the amplitudes of inter-annual variation for grassland and shrub regions were higher than those for cropland and barren/sparsely vegetated lands (Figure S2).For example, precipitation increased abnormally in 1998 for grassland and decreased abnormally in 1991 for shrubland.From the SPEI perspective, grassland (Figure S2b), cropland, and barren or sparsely vegetated lands (Figure S2d) presented a wetter trend during the first 15 years and a drier one for the last 15 years, except for 2010.Annual drought was more severe in 1999, 2000, 2005, and 2009, spring drought in 2000 and 2009, and autumn drought in 1999.

Spatial-Temporal Patterns of Vegetation Phenology and Productivity
In semi-arid regions, phenology showed different mean values via two different methods.SOS generally falls in the range between DOY 122 and 160 (Figure 3(a-1)) based on the threshold method, whereas it advanced for 10-20 days based on the derivative method (Figure 3(a-5)).Similarly, EOS also advanced for around 10 days via the derivative method, with a range of DOY 260-290 (Figure 3(a-2),(a-6)).Consequently, the retreat of LOS for the derivative method is longer than that for the threshold method in Mu Us Sandy Land and north of Yinshan Mountains due to advanced SOS.LOS ranged from 121 to 180 days in most semi-arid regions (Figure 3(a-3),(a-7)).
Spatial distributions of vegetation phenology, however, are fairly similarly derived from the threshold and derivative method.For instance, late SOS (e.g., after DOY 160) regions mainly surrounded the Horqin Sandy Land.Early EOS (260 < DOY < 280) regions mainly included the whole Horqin Sandy Land, the Da Hinggan Ling Mountains, the Yinshan Mountains, and the middle of the Mu Us Sandy Land, while late EOS has a scattered distribution in the study area (Figure 3(a-2),(a-6)).The iNDVI showed the same spatial patterns over this region.It decreased from east to west, except along the Da Hinggan Ling Mountains (Figure 3(a-4),(a-8)).In contrast, SOS, EOS, LOS, and iNDVI took on more obvious latitude and longitude trends via the threshold method than via the derivation method.
From the perspective of phenology dynamic trends, the results of the two methods were almost consistent, even though the rates of variation show little difference.SOS caused delays of one to five days per decade (1-3 days for m1, 3-5 days for m2, p < 0.05) in the Horqin Sandy Land and the eastern Ordos Plateau.Advanced SOS regions (Figure 3(b-1),(b-4)) were found along the Da Hinggan Ling Mountains with the rate of three to five days per decade (p < 0.05).EOS caused an advance of two days per decade (p < 0.05) in the middle and east of Inner Mongolia, and a delay of two days per decade (p <0.05) around the junction regions of Inner Mongolia and Shanxi province during 1982-2012 (Figure 3(b-2),(b-4)).The shortened LOS surrounded the Horqin Sandy Land and local places of the Yinshan Mountains.Other regions showed a decreased trend from east to west (Figure 3(b-3),(b-7)), that was similar to the delayed SOS regions (Figure 3(b-1),(b-3)).A significantly decreased area of iNDVI was mostly distributed in the north of semi-arid regions, while in the south of semi-arid regions, such as the Mu Us sandy land and Loess Plateau, an increased trend of iNDVI was found to be significant via the threshold method, and insignificant via the derivation method.

Response of Phenology and Productivity to Drought
In semi-arid regions, the significant correlations between SOS anomalies to spring drought were almost negative (8.28% for m1, 11.04% for m2, Table 1), and were located in the north of semi-arid regions (Figure 4).These regions, however, were not completely consistent with delayed SOS regions (Figure 3b).The significantly positive responses of EOS to autumn drought were in the southwest and middle-east of the study area and two-sides of the agro-pastoral transition zone, similarly to the EOS delayed or advanced regions (Figure 3b).Compared with EOS, an opposite result was found for LOS.Some shortened growing season length regions did not show a significant correlation with SPEI, e.g., in the region of the Horqin Sandy Land (Figures 3b and 4).Similar to EOS, iNDVI presented a significant positive correlation with SPEI in the north of semi-arid regions, where productivity had an obvious negative trend (Figures 3b and 4).It should be noted that the significantly correlated area of m1 is smaller than that of m2 for SOS, and larger for EOS and iNDVI.The difference between the two methods for correlated patterns of LOS and SPEI is obvious except in the Hulun Buir Sandy Land.However, the sensitivity of phenology and productivity to drought disturbance was consistent based on these two methods.The following relation holds: SOS < LOS < EOS < iNDVI.

The Phenology and Productivity among Different Land Cover Types
Since there are more obvious latitude and longitude trends of phenology, a more significant iNDVI dynamics area, and a larger correlated area than that of derivation method, the threshold method was chosen to analyze the dynamics of phenology and productivity among different land cover types.

The Spatio-Temporal Trends of Phenology among Different Land Cover Types
At land cover type stable regions, grassland, shrub, and barren or sparsely vegetated land did not show a significant trend of phenology and iNDVI, except in the case of cropland (Figure 5).Combined with climate analysis, EOS advanced in most autumn drought years at shrub, cropland and barren or sparsely vegetated land, such as in 1999.Grassland and cropland showed that SOS caused a delay in the spring drought years (e.g., in 1994 for grassland, 2001-2004 for cropland) and LOS was shortened in the 1999 drought year (Figure 5, Figure S2).Additionally, most significant changes in the regions of phenology have high chances of land cover variation (Figure 3b, Figure S1).

The Response of Phenology and Productivity among Different Land Cover Types
Among the different land cover types, dramatic heterogeneity phenology parameters (SOS and EOS) with dry-wet conditions were found (Table 2).An insignificant negative correlation between SOS anomaly and SPEI_May presented in grassland, barren/sparsely vegetation, and cropland.EOS anomalies have a significant positive correlation with SPEI_Oct in grassland and barren/sparsely vegetated lands (p < 0.05); LOS has a significant positive correlation with SPEI_ year for grassland (p < 0.05).Similarly, iNDVI also showed a positive correlation with SPEI_year at grassland, shrub, and barren/sparsely vegetated land, and the correlation coefficients are higher than those of LOS (Tables 1  and 2; Figure 5).This revealed that vegetation productivity presented much more sensitivity than phenology metrics.

Discussions
In this study, we analyzed the spatio-temporal changes of vegetation phenology and productivity, and investigated the responses of phenology and productivities to drought disturbance among different land cover types.The results showed that vegetation phenology and productivity have different trends and responses to drought disturbance.

Spatio-Temporal Trends of Vegetation Phenology and Productivity
Vegetation phenological metrics showed dramatic spatial heterogeneity over the study area, indicating the complicated effects of climate, vegetation, and soil conditions on phenology [17].Certainly, different methods increase the heterogeneity and uncertainty of phenology.In this study, the trends have fairly similar patterns based on two different methods, even though the rates are slightly different.For instance, the majority of advanced SOS trends occurred along the Da Hinggan Ling Mountains, with a minority occurring along the Yinshan Mountains.Delayed EOS is mainly distributed in the junction regions of Inner Mongolia and Shanxi province (Figure 3b).These results are supported by the conclusions of Wu et al. [17], Cong et al. [57], and Liu et al. [23].On account of different vegetation types, climatic zones, and water-heat conditions, the positive and negative change trends of phenology have different rates: SOS advanced by three days per decade (p < 0.05), and EOS delayed by two days per decade (p < 0.05) in semi-arid regions of Northern China.These results strongly support previous findings of phenology trends in the Northern Hemisphere [58] and China [59], but the results are lower than in the Inner Mongolia grassland regions [44].For grassland, advanced EOS was confirmed by Yang et al. [59], whereas phenology variation in shrubland was slightly different than that described by Liu et al. [23].The positive and negative trend of iNDVI during the growing season based on the threshold method is consistent with the results of Gong et al. [19] and Cao et al. [60].Hence, compared with the derivative method, the threshold method may be relatively effective for phenology analysis in the semi-arid region of northern China.

The Responses of Vegetation Phenology and Productivity to Drought
Commonly, temperature and precipitation are widely used to reveal the impacts of climate change on vegetation phenology and productivity [15][16][17][18][19]. Here, focusing on the responses of phenology and productivity to dry-wet conditions based on SPEI, several new insights were obtained by our results.First, drought could explain the anomalies of phenology and productivity (Figure 4).Spring drought delayed SOS and autumn drought significantly advanced EOS, which confirms our assumption and lies in the dynamics of hydrothermal conditions [48].As we mentioned in introduction, drought shifts the phenology by carbon and water processes [7,11,12].Due to spring drought, there is not enough available soil water for carbon synthesis during the photosynthesis process, so vegetation grows much more slowly than the normal situation [61,62].In contrast, on the one hand, autumn drought prompts leaves to close the stoma to decrease the transpiration and photosynthesis rates.On the other hand, respiration maintain higher rates to accelerate carbon degradation, and vegetation starts to stop growing earlier than normal years [62,63].Similar results also occurred in the cropland/grassland mosaic of North East China Transect [30].
Second, the sensitivity of productivity on drought disturbance is larger than vegetation phenology (Tables 1 and 2 and Figure 4).Productivity is the direct response to drought, while phenology is the indirect response to drought.It should be noted that the uncertainty of productivity assessment based on NDVI and phenology parameters might overestimate the sensitivity.Additionally, some researches realize the positive correlation between SOS and EOS [64,65].Nonetheless, it should be noted that the influence of earlier SOS on the determination of EOS was weaker than climatic variables [64,65].Although we did not take the relationship into consideration for phenology analysis, it has been noted that negative correlations were found in barren/sparsely vegetated land during 1982-2012 (Figure 5).

Analysis at Different Land Cover Types
The most significant regions of change in phenology occurred not in stable land cover type regions, but in the regions with high chances of land cover variation (Figure 3b, Figure S1).This verifies that land cover is an important factor for phenology change [51].Significant phenology dynamics in stable cropland may result from human disturbance.Different land cover types showed different responses to drought (Table 2).Because water is the main limiting factor in the grassland of Northern China's semi-arid regions, grassland showed a much more significant response to drought.In contrast, irrigation could help the phenology of cropland to maintain stable changes.

Uncertainty Analysis
One of the major uncertainties associated with satellite-derived phenology lies in the different methods.As the results showed, the phenology differences of two methods are obvious (Figure S3).Both mean SOS and EOS advanced by 10 days via the derivative method, compared with the dynamic threshold method (Figure 3a).The SOS values derived by the threshold method and EOS by the derivative method agree with that of Wu et al. [20].The EOS derived by the threshold method is similar to that of the Double Logistic methods by Liu et al. [57] and that of the Polyfit method by Yang et al. [59].Due to lack of field investigation data, it is difficult to decide which method is best for different land cover types and different study areas.Fortunately, several researches focus on productivity dynamics.The consistent conclusion is that the vegetation productivity of Mu Us sandy land has significant increased trends [19,66], which agree with our results based on the threshold method.Furthermore, more obvious latitude and longitude trends of the mean values of phenology via the threshold method are related to thermal and hydrological conditions.Therefore, relatively speaking, the dynamic threshold method is much more reliable in the semi-arid region of northern China.
Another uncertainty results from satellite data sources.GIMMS NDVI3g, SPOT-VGT, and MODIS NDVI data are the most commonly used for monitoring phenology [51].Due to the characteristics of resolution, different NDVI data sources have different accuracies of phenology.Some showed large discrepancies for SOS, comparing GIMMS NDVI 3g NDVI and MODIS NDVI data in Europe [67], equatorial areas, Artic areas, and arid areas [68].Other research considered that long-term GIMMS data was more accurate than SPOT-4 VGT NDVI [69], showing consistency and strong correlations between GIMMS NDVI3g and MODIS NDVI in the semi-arid West Africa zone [70,71] and China [68].Furthermore, the coarse time resolution of GIMMS NDVI3g may underestimate the magnitude of phenology variation.Daily data, however, are not sufficient to detect the change accurately enough to match the weather data.Integrated GIMMS NDVI and MODIS NDVI might show spatial heterogeneity for different regions and vegetation types [72], enabling phenology analysis.

Figure 2 .Figure 2 .
Figure 2. Diagram of algorithm for deriving phenological metrics and Net Primary Productivity (NPP) proxy.(a) is the algorithm of Dynamic Threshold Method,(b) is the algorithm of Derivative Method Figure 2. Diagram of algorithm for deriving phenological metrics and Net Primary Productivity (NPP) proxy.(a) is the algorithm of Dynamic Threshold Method, (b) is the algorithm of Derivative Method.

Figure 3 .Figure 3 .
Figure 3. Phenology and productivity average values and spatio-temporal trends based on two different methods.Note, m1 is threshold method, m2 is derivation methods, (a) is the mean value of phenology and iNDVI during 1982-2012, (b) is the temporal trends of phenology and iNDVI during 1982-2012.The unit for phenology mean values and trends are day or days per decade, respectively Figure 3. Phenology and productivity average values and spatio-temporal trends based on two different methods.Note, m1 is threshold method, m2 is derivation methods, (a) is the mean value of phenology and iNDVI during 1982-2012, (b) is the temporal trends of phenology and iNDVI during 1982-2012.The unit for phenology mean values and trends are day or days per decade, respectively.

Figure 4 .
Figure 4.The correlation coefficients and p-values between phenological metrics, productivity, and SPEI, (a,b) is based on m1 and m2, respectively.

Figure 4 .
Figure 4.The correlation coefficients and p-values between phenological metrics, productivity, and SPEI, (a,b) is based on m1 and m2, respectively.

Figure 5 .
Figure 5. Phenology and productivity analysis among different land cover types based on threshold method.(a), (b), (c) and (d) represent shrub, grassland, cropland and barren/sparsely vegetated land.

Figure 5 .
Figure 5. Phenology and productivity analysis among different land cover types based on threshold method.(a-d) represent shrub, grassland, cropland and barren/sparsely vegetated land.

Table 1 .
Statistical results of significantly correlated area between phenological metrics, productivity, and SPEI based on threshold and derivation methods.

Table 2 .
Correlation coefficient between anomalies of SOS, EOS, LOS, iNDVI, and SPEI among different land cover types.