MODIS normalized difference vegetation index (NDVI) and vegetation phenology dynamics in the Inner Mongolia grassland

The Inner Mongolia grassland, one of the most important grazing regions in China, has long been threatened by land degradation and desertification, mainly due to overgrazing. To understand vegetation responses over the last decade, this study evaluated trends in vegetation cover and phenology dynamics in the Inner Mongolia grassland by applying a normalized difference vegetation index (NDVI) time series obtained by the Terra Moderate Resolution Imaging Spectroradiometer (MODIS) during 2002–2014. The results showed that the cumulative annual NDVI increased to over 77.10 % in the permanent grassland region (2002–2014). The mean value of the total change showed that the start of season (SOS) date and the peak vegetation productivity date of the season (POS) had advanced by 5.79 and 2.43 days, respectively. The end of season (EOS) was delayed by 5.07 days. These changes lengthened the season by 10.86 days. Our results also confirmed that grassland changes are closely related to spring precipitation and increasing temperature at the early growing period because of global warming. Overall, productivity in the Inner Mongolia Autonomous Region tends to increase, but in some grassland areas with grazing, land degradation is ongoing.


Introduction
Land degradation in arid and semiarid regions, such as the Sahel in Africa and temperate grasslands in Australia, has become a critical threat (Gisladottir and Stocking, 2005;Prober and Thiele, 2005;Sop and Oldeland, 2013).China's vast grassland has also suffered from land degradation, mainly in the northern and western cold areas over long periods (Chen and Tang, 2005;Li et al., 2014;Wang et al., 2013).The Inner Mongolia Autonomous Region (IMAR), which is located along the northern border of China, is 68 % grassland (Kawamura et al., 2005a).Typical and meadow steppes, which are mainly used for grazing and animal husbandry, are the primary grassland ecosystems found in the IMAR (Kang et al., 2007).In recent decades, the Chinese government has implemented a "return-farmland-to-grassland" strategy to reverse the land degradation in pastures (Shang et al., 2014;Wang et al., 2010).Vegetation restoration has been used to protect diverse degraded landscapes.A considerable number of studies suggest that vegetation cover helps mitigate soil erosion by stabilizing soil aggregation, reducing wind speeds and water erosion, improving soil porosity and increasing biological ac-Z.Gong et al.: MODIS NDVI and vegetation phenology dynamics in the Inner Mongolia grassland tivities in the soil (Fattet et al., 2011;Florentine et al., 2013;Lee et al., 2002;Lieskovský and Kenderessy, 2014).
To evaluate the vegetation restoration effect, anthropogenic and climatic impacts should be considered.Vegetation cover change represents the most direct response of vegetation to climate changes and human activities (Zhao et al., 2012).Akiyama and Kawamura (2003) analyzed the land cover change over 1979-1997 and indicated that the areas with productive grasslands decreased while low-productivity grasslands increased.The seasonal change in vegetation (phenology) is a key parameter for studying and analyzing climate change and vegetation responses (the feedback between the land surface and the atmosphere), which can improve the simulation quality of carbon, water, and energy exchanges between the atmosphere and the land surface (Ma et al., 2013).The vegetation productivity and phenology in the temperate region of China has already changed in response to the dramatic climatic changes (Jeong et al., 2011;Piao et al., 2006Piao et al., , 2010;;Peng et al., 2011).Earlier studies indicated that the recovery of vegetation from long-term degradation is related to the increase in precipitation (Eklundh and Olsson, 2003;Sop and Oldeland, 2013).The growing body of evidence suggests that climate warming has advanced the biological spring in temperate China (Chen et al., 2005;Piao et al., 2006;Zheng et al., 2002).Additionally, longer growing seasons, particularly earlier spring vegetation green-up, may significantly enhance the vegetation productivity in temperate and boreal regions (Cong et al., 2013;Hu et al., 2010;Kimball et al., 2004).
As traditional fieldwork is time-consuming and costly, remote sensing methods have been utilized as cost-effective approaches to detect vegetation changes at large spatial scales.Monitoring landscapes through satellite-derived vegetation indices (VIs), such as the normalized difference vegetation index (NDVI) and enhanced vegetation index (EVI), is a successful method for assessing vegetation conditions and phenology (Glenn et al., 2008;Zucca et al., 2015).Previous studies suggested that time series satellite data can reliably detect the phenology, forage quantity, and quality of grassland areas using the VIs derived from Advanced Very High Resolution Radiometer (AVHRR) data and Moderate Resolution Imaging Spectroradiometer (MODIS) satellite images (Kawamura et al., 2003(Kawamura et al., , 2005a, b, c), b, c).Significant delays in the vegetation green-up during 1982-1991 in the desert steppe and in part of the typical steppe of Inner Mongolia were detected using AVHRR NDVI (Yu et al., 2003).Moreover, although the signs of the trends in the vegetation green-up dates detected by various methods were broadly consistent spatially and for different vegetation types, large differences occurred in the magnitudes of the observed trends.The large variance obtained using different methods is particularly apparent for arid and semiarid vegetation types (Cong et al., 2012(Cong et al., , 2013;;Zhao et al., 2012).
Because of the complex vegetative species, diverse vegetation coverage and vast observational areas in semiarid con-tinental areas, limited research has been conducted in the IMAR to analyze the phenology shift using Terra MODIS NDVI time series.In this study, the objectives are to (1) evaluate the changes in land use and land cover, including vegetation cover and NDVI trends; (2) detect the vegetation phenology (the start/end dates, the maximum vegetation productivity dates and the length of the growing season) using Terra MODIS NDVI data; and (3) relate the trends of the NDVI and phenology to changes in the climate (precipitation and temperature) at 18 meteorological stations throughout the permanent grassland between 2002 and 2014.

Study area
The IMAR (Fig. 1), the third-largest province and autonomous region with plateaus in the country, occupies a total area of 1.18 million km 2 , or 12.3 % of the country's territory (Han et al., 2009).This area has a yearly average temperature of 0-8 • C, while the mean annual precipitation, which is concentrated in June and July, varies from less than 100 mm in the western Gobi Desert to more than 400 mm in the northeast forests and meadow steppes (Han et al., 2009;John et al., 2008;Yu et al., 2003;Zhou et al., 2006).The growing season for perennial species typically starts in April and ends in September.Annuals germinate from April to July due to the available soil moisture and subsequent rain events (Bai et al., 2004;John et al., 2008).Typical and meadow steppes are the major types of grassland ecosystems in the IMAR; these areas are often used for grazing and animal husbandry (John et al., 2008;Kang et al., 2007).Within the semiarid typical steppe, which has an annual precipitation of less than 350 mm, the dominant vegetation species are Stipa grandis, Leymus chinensis, and multiple species of Artemisia spp.and Festuca spp.(Kang et al., 2007;Yu et al., 2003).In comparison, meadow steppes develop in areas with a higher productivity, with an annual precipitation of approximately 450 mm; the primary plants are S. baicalensis, L. chinensis, and Cleistogenes mucronata (John et al., 2008;Kang et al., 2007).Moreover, the desert steppe is in an arid ecosystem, in which the annual precipitation is 150-200 mm; this steppe has the lowest biomass production, but it is also important to the grassland area in the IMAR.S. krylovii, S. bungeana and A. ordosica are dominant (Cheng et al., 2001;John et al., 2008;Kang et al., 2007;Yu et al., 2003).In the present study, 18 sampling sites throughout the permanent grassland area, according to the MODIS Land Cover Type Yearly product (MCD12Q1, 500 m) for 2002-2014, were selected to represent the major steppe areas in the IMAR., 6, 1185-1194, 2015 www.solid-earth.net/6/1185/2015/

Data
The meteorological data, including monthly precipitation and temperature, were acquired from the China Meteorological Data Sharing Service System (http://www.escience.gov.cn/metdata/page/index.html).Eighteen meteorological stations (Fig. 1 and Table 1) throughout the permanent grassland were selected to represent the diverse steppes in the IMAR.Two types of products derived from Terra MODIS satellite data were employed.The 500 m spatial resolution MCD12Q1 product (MODIS Land Cover Type Yearly L3 Global 500 m, Version 5 -which is generated by the International Geosphere-Biosphere Programme (IGBP) global veg-etation classification scheme from 2002 to the latest available year, 2012 -was used to detect the land use/cover change and to extract the permanent grassland area in the IMAR.MOD13Q1 (Vegetation Indices 16-Day L3 Global 250 m, Version 5) data from 2002 to 2014 were used to extract the NDVI for estimating trends in vegetation cover and phenology changes over the 11 years.The NDVI, which is a nonlinear combination of red and near-infrared (NIR) spectral radiances (NIR − red)/(NIR + red), exhibits a strong relationship with vegetation activity and green biomass.The index is usually employed for predicting vegetation phenology from space (Glenn et al., 2008;Hmimina et al., 2013;Karlsen et al., 2008).The NDVI time series were calculated as mean values over 6 × 6 pixels (1.5 km × 1.5 km) centered on each of the 18 meteorological stations to avoid land cover heterogeneity (Karlsen et al., 2008).

Reconstruction of NDVI time series data
After the preprocessing procedures (mosaic and reprojection), the data quality of the NDVI data was assessed using the corresponding quality assessment (QA) information that describes the utility of the VI values.When the VI usefulness = 1101, the quality of the pixel is low and not useful.If the VI usefulness = 1110, the data also cannot be used due to the faultiness in the original level data.And VI usefulness = 1111 means useless for other reasons.The invalid data (VI usefulness = 1101, 1110, 1111) were eliminated and interpolated linearly.
Then the NDVI time series were temporally smoothed by a Savitzky-Golay filter (Chen et al., 2004), which provides a simplified least-squares-fit convolution for smoothing and computing derivatives of a set of consecutive values (such as a spectrum) (Fig. 2).The S-G filter performs best in most situations when smoothing different vegetation types using various satellite data (Geng et al., 2014).The Savitzky-Golay filter is computed as follows: where Y is the original NDVI value, Y * is the resultant NDVI value, and C i is the coefficient for the ith NDVI of the smoothing window.N is the number of convoluting integers and is equal to the smoothing window size (2 m + 1).j represents the running index of the ordinate data in the original data table, and m represents the half width of the smoothing window (Savitzky and Golay, 1964).Two assumptions must be confirmed: the time series follows an annual cycle of growth and decline, and clouds and poor atmospheric conditions decrease the NDVI values.Therefore, sudden decreases in the NDVI, which are not compatible with the gradual process of vegetation change, should be regarded as noise and removed.Two key parameters, i.e., the half-size of the filtering window (m) and the degree of the filtering polynomial (d), determine the final effect.In this study, the m was set at 4 and the d was set at 2 to achieve the best-fitting result.Through using the smoothed NDVI time series data (299 values/pixel from 2002 through 2014), linear models were developed.The trends in the NDVI changes were quantified by the slope of the regression line derived from the linear models of the cumulated NDVI time series against time.The slopes were tested for significance with a significance level of 0.05.

Detecting phenological stages from the NDVI time series
The start and end dates of a growing season are usually selected as the indicators of phenology shifts.The globally constant NDVI threshold may be suitable for forested ecosystems but may also surpass the peak value of the NDVI in semiarid grasslands (White et al., 2003;White and Nemani, 2006).Thus, we determined the start/end dates (SOS/EOS), peak vegetation productivity date of growing season (POS), and the length of growing season (LOS) from the MODIS NDVI time series data by modifying the method applied by Butt et al. ( 2011): (i) the inflection point (the maxima of the second derivative) during the spring (from March to May) was identified as the SOS, while another inflection point (the maxima of the second derivative) during the autumn (from August to October) was identified as the EOS; (ii) the POS was defined as the date of the maximum NDVI during the growing season (Ma et al., 2013); and (iii) the LOS was defined as the difference between the EOS and SOS (Fig. 3).
The trends of the phenological stages were estimated by regressing the SOS, POS, EOS and LOS against years over the study period.

Correlation between meteorology and the phenology
Previous research indicated that only the precipitation and temperature in the spring and early summer determine vegetation growth (Bai et al., 2004).In the present study, the relationships between climate data (monthly and cumulative precipitation and mean monthly temperature) and the phenology dates (SOS, EOS, POS, LOS) were assessed using the population regression function.

Vegetation cover change
The primary land use and land cover change from 2002 to 2012 according to the latest MODIS land cover data are shown in Table 2. Furthermore, the standard deviation (SD) and coefficient of variance (CV = SD/mean × 100 %) were calculated to monitor the fluctuations of the land cover changes over these years.
Grassland was the most extensive land cover type in the IMAR, occupying 49.67 × 10 4 km 2 (53.36 %) in 2002 and 50.24 × 10 4 km 2 (53.97 %) in 2012.These values are smaller than the previous data, but a slight increase occurred during the study period.The forest area also expanded by 2.38 × 10 4 km 2 within this period.In contrast, the barren land, artificial land (crop land and urban area) and water bodies declined.Following grassland, the second-most-extensive cover of barren land decreased from 23.08 × 10 4 km 2 (24.79 %) to 21.17The permanent grassland (2002-2012) was then extracted, yielding an area of 37.57 × 10 4 km 2 (40.36 % of the total IMAR area) (Fig. 1).The cumulative annual NDVI trend was analyzed within the study area (Fig. 4).A positive trend occurred over 28.97 × 10 4 km 2 (77.10 % of the permanent grassland), whereas a negative trend was observed over 8.60 × 10 4 km 2 (22.90 %).A significant heterogeneous spatial distribution was observed in the cumulative annual NDVI change.The negatively changing area was mainly located in the center and east of the IMAR, where typical steppes are dominant.

Yearly change in precipitation and temperature
Yearly changes in the annual precipitation and mean temperature at 18 stations are shown in Fig. 5.The mean annual precipitation ranged from 221.90 to 404.38 mm year −1 , with a minimum value in 2007 and a peak value in 2012.Among the 18 stations, large variances of the precipitation SD were observed (56.47-143.67 mm year −1 ).The annual mean temperature ranged between 3.58 and 6.10 • C, with a minimum value in 2012 and a peak value in 2007.Meanwhile, the SD of the monthly temperature at the stations ranged from 2.23 to 2.67 • C over the study period.An upward trend in the mean annual precipitation and a slight downward trend in the mean monthly temperature occurred at the selected stations.

Phenology trends
The statistical results of the average phenology dates (SOS, POS, EOS, LOS) at all 18 meteorological stations in various years are shown in Table 3. Overall, the SOS at the 18 stations varied from mid-March to May, with a mean annual value of 122.21 ± 20.81 day of year (DOY), where 20.81 was the annual SD of the 18 stations (SD stations ).The POS occurred between July and August, with a mean annual value of 214.54 ± 9.33 DOY.The EOS varied between late August and October, with a mean annual value of 285.24 ± 17.67 DOY.The mean annual LOS was 164.03 ± 26.63 days.The SD stations in the vast region indicated significant spatial heterogeneity of the phenology across the IMAR grassland.The SD of the mean annual phenology dates was also computed as SD year .The SD year for the SOS was 5.89 days, while the SD year for the EOS and LOS equaled 4.47 and 9.03 days, respectively.The SD year for the POS (3.48 days) indicated the smallest yearly fluctuation of all of the phenological parameters.
After applying a linear analysis procedure, the average SOS and POS of the 18 stations advanced by 5.79 and 2.43 days, respectively, from 2002 to 2014.The average EOS of the stations was delayed by 5.07 days.Finally, the average LOS of the stations was extended by 10.86 days from 2002 to 2014 (Fig. 6).

Linkages between the phenology and meteorological data
The selected regression models at different phenological stages (SOS, POS, EOS, and LOS) between precipitation (monthly and accumulated value during different periods) and temperature (monthly and mean value during different periods) are presented in Table 4.The delayed effect from climate in SOS was obviously detected.Generally, the SOS negatively correlated with the cumulative precipitation, especially during the growing season in the last year (March- September, R 2 = 0.95, P < 0.001).Furthermore, the temperature in the last May also negatively correlated with SOS strongly (R 2 = 0.73, P < 0.05).Therefore, the increasing cumulated precipitation and the temperature in May in the last year could advance the SOS.The POS negatively correlated with the cumulative precipitation from May to June (R 2 = 0.62, P < 0.001) and positively correlated with the mean monthly temperature from June to July (R 2 = 0.52, P < 0.001).The increasing precipitation and colder weather can advance the peak vegetation activity date.The EOS positively correlated with the precipitation in the last August (R 2 = 0.68, P < 0.05) but correlated with the temperature in March with lower significance (R 2 = 64, P = 0.06).Thus the delay in the senescence of vegetation was considered because of the wetter autumn in the last year and the warmer spring in the current year.Overall, the LOS positively correlated with the cumulative precipitation from April to May (R 2 = 0.49, P < 0.05) and mean temperature from January to March (R 2 = 0.72 P < 0.001), indicating that the wetter and warmer weather during the early vegetation growing period could extend the LOS.

Discussion
This study investigated the change in the cumulative annual NDVI and phenology during the growing season of the permanent grassland in the IMAR.The phenological dynamics were correlated with the local meteorological variations.
Our results indicated that the cumulative annual NDVI had a positive trend mainly in the northern and western regions (Fig. 4).In the west, desert steppes are dominant.Previous research has reported a close relationship between vegetation changes and climate factors, particularly precipitation and temperature (Cao et al., 2013).The upward trend in the annual precipitation is considered the main meteorological factor that led to the cumulative annual NDVI increase in the Solid Earth, 6, 1185-1194, 2015 Regarding the phenology changes (Table 3) and trends (Fig. 6), the extension of the growing season was detected, as also revealed by other recent research in temperate China (Bai et al., 2004;Cong et al., 2013).An earlier onset of the start is most prominent (Linderholm, 2006).Referring to previous studies, the SOS was obviously earlier in the present study period (from mid-March to May).As Lee et al. (2002) reported that the onset of green-up for typical and desert steppes occurred from early May to early June (1982 to 1990), the SOS since 2002 had advanced by approximately 1 month compared with the previous 20 years.This finding was consistent with the previous suggestion that the SOS advanced between the 1980s and the 2000s in major biomes in China, according to various approaches (Cong et al., 2013;Ma and Zhou, 2012).However, limited research had investigated the senescence date of the growing season.Our results have provided the valuable finding that, besides of the advance in SOS, the delay in EOS (5.07 days) also contributed the extension of the growing season.
Some researchers have indicated that the phenology could be influenced by the climate several months before (Estrella and Menzel, 2006;Miller-Rushing and Primack, 2008).In our results, this delayed effect has been found.Our results were in agreement: global warming could promote the vegetation growth and extend the growing season (Linderholm, 2006).The temperature has increased significantly, particularly since the 1980s (Ding and Chen, 2008;Gao et al., 2009).However, from 2002 to 2014, the IMAR grassland tended to be slightly colder in the spring (from January to May).We speculated that the increasing precipitation might be the main driving factor of the advance in SOS and the delay in EOS.Nevertheless, previous work revealed that precipitation decreased slightly over the last 50 years compared with the obvious interannual change.Thus, the precipitation appeared to increase over the recent decade.Rather than the change in temperature, the wetter weather conditions were considered the main reason for the phenology change in the IMAR.
Our results indicated that plant productivity in the IMAR increased, but in some areas with grazing, land degradation is ongoing.As different vegetation types have various response to the climate change (Yuan et al., 2007), the spatial heterogeneity of the phenology dynamics needs to be detected.Meanwhile, the present study period only covered the last 13 years, limiting its conclusiveness regarding the change in the IMAR grassland ecosystem during a longer period.Thus further research should be conducted to identify the correlation between changes in phenology and meteorology.

Conclusions
In this study, we examined recent trends in the NDVI and phenology changes in the Inner Mongolia grassland using Terra MODIS time series data.The relationships between phenology change (SOS, EOS, POS and LOS) and climate data (precipitation and temperature) were also evaluated.The following conclusions can be drawn from the study: 1.The positive trends of the cumulative annual NDVI (77.10 %) could be interpreted as an increase in plant productivity in the Inner Mongolia permanent grassland.
2. The advance in the SOS and the delay in the EOS extended the LOS in Inner Mongolia between 2002 and 2014.
3. The increase in precipitation is the main factor for the extension of the growing season.
Overall, the results reveal recent trends in the Inner Mongolia grassland and their correlation with climate data.Further analysis using long-term satellite data and climate data should be conducted.
× 10 4 km 2 (22.74 %) by 2012.Based on the SD and CV values during 2002-2012, the artificially controlled area (crop land and urban area) decreased the most (SD = 1.31 × 10 4 km 2 ; CV = 1.30%), reflecting the effect of the return-farmland-to-grassland strategy of the Chinese government.Grassland and barren land had the smallest fluctuations.Water bodies were relatively stable among the land cover types.The changes in the land use dynamics indicated a positive trend in the total vegetation coverage.

Figure 4 .
Figure 4.The cumulated NDVI trend in permanent grassland from 2002 to 2014.The green parts represent the increasing area, while the red parts represent the decreasing area.

Figure 5 .
Figure 5.The annual total precipitation and mean temperature of 18 meteorology stations.

Figure 6 .
Figure 6.The trend of phenology dates from 2002 to 2014.SOS is the start of growing season; POS is the maximum NDVI date during the growing season; EOS is the date of the end of season; LOS is the length of growing season.

Table 1 .
Information on the meteorological stations.
Figure 1.The Inner Mongolia Autonomous Region (IMAR).The gray part represents the permanent grassland area, whereas the black dots represent the selected meteorology stations.
Figure 3.The extraction of phenology dates from NDVI.SOS is the start of growing season; POS is the maximum NDVI date during the growing season; EOS is the date of the end of season; LOS is the length of growing season.

Table 3 .
Annual mean and standard deviation (SD) of phenological stages at 18 meteorological stations during 2002-2014.The units are day of year (DOY); SOS, POS, EOS and LOS are the start, peak, end and length of the growing season; Mean is the mean value of the phenology dates; SDy is the standard deviation of the mean annual phenology dates.

Table 4 .
The climate variables most strongly correlated with phenology, and the corresponding parameters of its linear model.SOS is the start of growing season; POS is the maximum NDVI date during the growing season; EOS is the date of the end of season; LOS is the length of growing season.Prec and Temp represent the precipitation and temperature, respectively.