The Impacts of Weather and Conservation Programs on Vegetation Dynamics in China ’ s Loess Plateau

We present an analysis of the impacts of weather change and large-scale vegetation conservation programs on the vegetation dynamics in China’s Loess Plateau from 2000 through 2009. We employed a multiple lines of evidence approach in which multi-scale data were used. We employed Normalized Difference Vegetation Index (NDVI) data, acquired by the Moderate Resolution Imaging Spectroradiometer (MODIS) at 500 m to identify significant vegetation increases in the Loess Plateau since 2000. We found increases in NDVI for 48% of the Loess Plateau between 2000 and 2009. We were able to attribute up to 37.5% of the observed vegetation increases to weather change, vegetation conservation activities and crop yield increases. We demonstrate that the impact of vegetation conservation programs on vegetation change in the Loess Plateau is twofold. On the one hand, vegetation conservation programs target marginal lands. Thus, significant vegetation increases due to cropland conversion and afforestation can be found in these regions. On the other hand, intensified agricultural production can be found in croplands with suitable topography and well-established irrigation systems, which were not enrolled in conservation programs to offset the agricultural production loss caused by vegetation conservation programs elsewhere.


Introduction
Forest transitions refer to the change from a net loss in forest cover to a net gain of forested areas within a specific region [1,2].The occurrence of forest transitions can improve ecosystem services by increases in forest cover or woody biomass [3].For instance, expanded forest cover can improve water quality by reducing sediment discharges resulting from soil erosion [3].Woody biomass increases can lead to increased carbon sequestration [3].There are two pathways identified to facilitate the development of forest transitions: the "economic development path" and the "forest scarcity path" [3].The "economic development path" describes forest regrowth on abandoned agricultural lands as a result of a decreasing agricultural population, which migrates to urban areas during the process of industrialization and urbanization [3,4].The "forest scarcity path" is a result of elevated prices of forest products due to the forest scarcity caused by previous agricultural expansion [3,4].
Previous studies on forest transitions in China reveal distinct characteristics in terms of the driving forces and the impacts of forest transitions.Rather than the loss of agricultural population and rising prices of forest products, large scale conservation programs, launched in response to the deteriorated natural environment, were found to be the most significant driving force of Chinese forest transitions [4].Historically, China has suffered from deforestation and widespread soil erosion [5].Severe soil erosion associated with degraded vegetation conditions was recognized as the major factor in causing the massive floods in China during the summer of 1998 [6,7].In order to mitigate the degraded environment by protecting natural forests and reducing soil erosion, the Chinese government launched two of the largest conservation programs in the world: the Natural Forest Conservation Program (NFCP) and Grain to Green Program (GTGP) [8].NFCP was initiated in 1998 and designed to be implemented in multiple stages between 1998 and 2050 [8].There are three key goals associated with NFCP: (1) restoring natural forest by issuing logging bans and mountain closure, (2) shifting the major source of timber production from natural forests to plantation forests, and (3) reducing soil and water loss by afforestation [8,9].The NFCP has been established in 18 provinces and autonomous regions in China since 2000 with a total investment of about 61 billion Yuan (~9.7 billion dollars), from 1998 to 2005 [8,9].GTGP was launched in 1999 in three pilot study provinces-Gansu, Shanxi, and Sichuan [8], and expanded to 25 provinces and autonomous regions by 2002 [8].The goal of GTGP is twofold.The primary goal is to mitigate soil erosion by increasing vegetation cover through converting steep slope agriculture (croplands on slopes ≥ 15° in Northwestern China and ≥25° in areas other than Northwestern China) to forest or grassland and to carry out afforestation in barren areas [8].GTGP also aims to fight poverty and improve the socioeconomic well-being of participating farmers [8].GTGP received an investment of more than 90 billion Yuan (~14.3 billion dollars) by 2005 and aimed to achieve an increase in vegetation cover of 320,000 km 2 by 2010 [8].
In this study, we focused on the Loess Plateau which is an area covered by thick and highly erodible loess in Central-Northern China [10] (Figure 1).Severe soil erosion affects about 70% of China's Loess Plateau [5].Intensive cultivation on marginal land (e.g., land with steep slopes) without effective soil loss protection was identified as one of the major factors contributing to the high erosion rates within this region [10].Severe soil erosion results in significant loss of fertile land resources and imposes food security threats to this region [11].The Loess Plateau was a first priority area for NFCP and two of the three pilot study provinces of GTGP (Gansu and Shanxi) [8,9] were located here.The Loess Plateau is frequently included in studies which investigate the effects of NFCP and GTGP on vegetation cover change [12,13], desertification control [14], and cost-effectiveness of the programs.Typically these studies are only conducted in a small portion of the plateau (e.g., a few counties) and there has been no evaluation of the impacts of NFCP and GTGP on vegetation development in the entire area.In addition, there are debates over the effectiveness of NFCP and GTGP in improving vegetation conditions in the Loess Plateau.For example, Zhou et al. [13], and Cao et al. [12], attributed the increase in vegetation cover in Shaanxi Province to GTGP.Afforestation failure by GTGP, however, was also found in some counties of Shaanxi Province due to inappropriately selected tree species [12].In this paper, we employed a multiple lines of evidence approach to investigate the effects of large-scale vegetation conservation programs on the vegetation development in China's entire Loess Plateau.We carried out five different analyses to determine the effect of these large-scale vegetation conservation programs on the vegetated land surface: (1).We employed satellite derived vegetation indices with a spatial resolution of 500 m to identify significant vegetation changes since 2000.(2).We examined monthly temperature and precipitation data to determine if significant weather changes were contributing to the observed vegetation increases.(3).We employed satellite derived land cover data and vegetation conservation data reported by statistical yearbooks to investigate if the observed vegetation increases can be attributed to land cover change caused by vegetation conservation activities.(4).We used anthropogenic biomes to evaluate how large-scale vegetation conservation programs choose target regions.(5).We investigated crop yield data to investigate the effects of vegetation conservation programs on croplands in areas with relatively shallow slopes.

Loess Plateau
The Loess Plateau, with an area of about 626,000 km 2 , is located in the middle reaches of the Yellow River.Elevation of the Loess Plateau ranges from 1,000 m to 1,600 m above sea level [16].The Loess Plateau extends from 101°E to 114°E and 34°N to 41°N, covering the entire area of Ningxia Hui Autonomous Region and Shanxi Province and parts of the Shaanxi, Gansu, Henan, and Qinghai Provinces and Inner Mongolia Autonomous Region (Figure 1).The Loess Plateau hosts a population of about 86 million mostly rural residents [17].The provincial capitals located in the Loess Plateau include Hohhot, Lanzhou, Taiyuan, Xi'an, Xining, and Yinchuan.The Loess Plateau features a northwest-southeast climate gradient.The mean annual precipitation increases from 300 mm in the northwest to 700 mm in the southeast [18].The mean annual temperature ranges from about 4 °C in the northwest to 14 °C in the southeast [18].

Study Regions for Conservation Program Analysis and Crop Yield Analysis
Besides studying the entire Loess Plateau, we also selected 20 counties in the Northern Loess Plateau to analyze the effects of vegetation conservation programs on vegetation dynamics in more detail.Of these 20 counties, four were in the region of Ordos, Inner Mongolia, nine counties were in the region of Yulin, Shaanxi Province and the remaining seven counties were located in the region of Yan'an, Shaanxi Province (Figure 1).In addition, we selected nine counties as case study areas to investigate the effects of GTGP on the crop yield of croplands that were not directly participating in GTGP (Figure 1).
Detailed explanations of the selection criteria for the case study counties for the vegetation conservation and crop yield analyses are available in Sections 3.4 and 3.6, respectively.

Data and Methods
Determining the effect of large conservation programs over large areas is not straightforward.There are several different factors that play a role in the location and the success of these programs.In this study we took a multiple lines of evidence approach to determine the effect of the large scale conservation programs on the vegetated land cover in the Loess Plateau.First, we determined where vegetation is changing according to remotely sensed vegetation index data (Section 3.1).Then we went through multiple steps to attribute the vegetation increases to one of the following potential causes: weather (Section 3.2), land cover change (Section 3.3), vegetation conservation activities reported in statistical yearbooks (Section 3.4), the effect of GTGP and NFCP on elevated areas and areas with steep slopes (Section 3.5), and, finally, we investigate how crop yield has changed as an indirect result of conservation programs (Section 3.6).An overview of the data used in the following seven sections can be found in Table 1.We employed a stepwise approach rather than a conventional multiple regression procedure to attribute observed vegetation increases to different factors because data availability was variable in space and time and, thus, a unified multiple regression approach could not be used.For example, for any location we were interested in, a multiple regression procedure would require us to have a complete set of response and predictor variables.However, the crop yield and vegetation conservation data only covered a limited portion of the area with observed vegetation increases.
The MODIS reflectance dataset is available as 10° × 10° latitude/longitude tiles and the Loess Plateau is covered by four tiles (h26v4, h25v5, h26v5, h27v5).The MODIS reflectance dataset is originally produced every eight days based on 16 days of acquisition.After calculating NDVI for each original reflectance composite, we generated 16-day NDVI composites by averaging the original consecutive 8-day rolling NDVI composites.With a temporal resolution of 16 days there are a total of 227 NDVI composites for the time period from 2000 to 2009 (20 NDVI composites for year 2000 and 23 composites for each year, from 2001 to 2009).
We applied the seasonal Mann-Kendall (SMK) trend test to estimate vegetation change from the 272 NDVI composites between 2000 and 2009.The SMK test statistic was calculated based on the rank of a season (in this case, a 16-day period) within the time series of the same season across all years [29].Two steps were carried out to calculate the SMK test statistic: (1) the Mann-Kendall (MK) test statistic is the sum of the number of times each season has a higher rank than the same season in any previous year; (2) The SMK test statistic for the entire time series is calculated by summing the MK test statistic for each season [29].SMK is then corrected for autocorrelation among seasons.The advantage of the seasonal Mann-Kendall trend test over simple linear regression in terms of estimating changes in NDVI time series was previously discussed by de Beurs and Henebry [29].The SMK trend test was applied to each 500 m pixel in the MODIS NDVI dataset.We interpreted the area having a positive SMK test statistic with a p-value lower than or equal to 0.01 and annual missing data rates lower than or equal to 40% as highly significant positive vegetation change.The area having a negative SMK test statistic with a p-value lower than or equal to 0.01 and missing data rates lower than or equal to 40% was interpreted as highly significant negative vegetation change.
Aside from applying SMK change analysis, we generated monthly NDVI maximum value composites to represent the overall vegetation condition for each month.We also calculated annual accumulated NDVI on a pixel basis by summing NDVI for the entire year using the 16-day NDVI composites.

Weather Data
Weather data used in this study included monthly temperature and precipitation records from 138 weather stations, spread over the Loess Plateau (Figure 1).The weather data were collected and compiled by the Chinese National Meteorological Information Center [21].The data were of high accuracy after eliminating errors by extreme value testing, temporal consistency testing and manual error checks [21].We applied kriging interpolation to generate monthly raster layers at a spatial resolution of 500 m for temperature and precipitation from 1999 to 2009.We then applied the SMK trend test to the time series of monthly temperature and precipitation data from 1999 to 2009 on a pixel by pixel base in order to examine whether there was significant weather change that could be related to the vegetation dynamics.We took the area with a p-value lower than or equal to 0.1 as a region showing significant weather change since 1999.
In order to determine whether there were significant impacts on vegetation change imposed by weather changes, we employed the Spearman's rank correlation coefficient to calculate the lagged correlation between the monthly NDVI maximum value composites and the precipitation in the previous month.We also computed the correlation coefficient between temperature and NDVI with the same monthly lag.
In addition, we employed the p-values associated with NDVI and weather trend analyses to assess the probability of error in terms of attributing NDVI increases to significant changes in temperature or precipitation.We calculated the error by estimating the probability that both changes in NDVI and weather were correctly estimated.

Land Cover Data
The MODIS global 500 m land cover product (MCD12Q1) was employed to derive the land cover information over the study period.As the quality of the data was higher for the years starting in 2003, MODIS land cover data from 2003 through 2008 were selected for this study.The MODIS land cover data provides multiple land cover classification schemes including: IGBP global vegetation, University of Maryland land cover classification, MODIS-derived LAI/fPAR, and the Plant Functional Type.We picked the IGBP classification scheme, which was accompanied by a separate data layer that provides the assessment of relative classification quality for each 500 m pixel.The classification quality ranges from 0 to 100 with a higher number representing a higher probability that the classification was correct.The IGBP classification scheme allowed an analysis of sixteen land cover categories including twelve vegetation classes and four classes of non-vegetated land.
Year-to-year land cover change inferred from direct comparison of the MODIS land cover data across years could contain substantial spurious changes due to inaccurate classifications resulting from potential mixed pixels at 500 m spatial resolution, the limited separability of certain land cover classes in the MODIS spectral-temporal space, phenological changes, and disturbances (e.g., fires) [20].In order to avoid characterizing the land cover change during the study period by direct comparison of the MODIS land cover data, we divided the six-year MODIS land cover datasets into two groups: one from 2003 to 2005, and the other from 2006 to 2008.For each group, we derived a three-year-summarized land cover map.Pixels that were classified as the same land cover type for at least two consecutive years were considered as stable.The land cover type of these pixels was assigned as the land cover type with the highest frequency during the three years.Other pixels were assigned to a new class labeled as random classification.In order to evaluate the efficiency of this grouping strategy, we compared the land cover assessment scores of the pixels that were considered as stable with those considered as random classifications.We employed a two-sample t test to determine if the mean land cover assessment scores of stable pixels were significantly higher than the scores of the random classification pixels.
We explored the land cover change for three broad land cover categories: natural vegetation, cropland, and barren or sparsely vegetated to identify vegetation dynamics that could be caused by large scale vegetation conservation programs.To create the natural vegetation class we combined the ten natural vegetation classes under the IGBP classification scheme.
We employed the relative classification quality scores to assess the potential errors in the land cover change analysis.We first calculated the average classification quality scores for each of the three land cover categories in the two time periods: from 2003 to 2005 and from 2006 to 2008.The error rate for any type of land cover change was calculated by estimating the probability that the land cover types in both time periods were correctly classified.
Besides investigating land cover changes according to MODIS data we also explored anthropogenic biomes (Anthromes).Anthromes provide a way of understanding the terrestrial biosphere in terms of form and intensity of human influence on biomes [15,30].Anthromes are organized according to population density with village anthromes supporting more than 100 people/km 2 and residential anthromes supporting between 10 and 100 people/km 2 .Populated (1-10) and remote anthromes (<1) support fewer people.Anthromes from the year 2000 were employed in this study to identify the preference of large-scale vegetation conservation programs in choosing target regions in the Loess Plateau based on population density and agricultural productivity.Anthromes are mapped at a spatial resolution of 0.083°.We examined the distributions of land cover changes that may be caused by vegetation conservation programs (e.g., croplands and barren lands converted to natural vegetation) within different anthromes.We compared the distribution of different types of land cover change (e.g., croplands converted to natural vegetation vs. croplands left unchanged) to understand how large-scale vegetation conservation programs choose target regions.

Vegetation Conservation Data from Statistical Yearbooks
In order to identify land cover change that occurred at scales finer than 500 m and thus may not have been identified in the land cover change analysis using MODIS derived land cover maps, we employed vegetation conservation data documented in regional statistical yearbooks.
Forestry and agricultural statistics documented in a regional statistical yearbook are published by the regional bureau of statistics using a methodology that integrates household survey and a bottom-up reporting system [22,31].The bottom-up reporting system works in such a way that the collection of the target statistics begins with the statistics reported by officials at the bottom level (e.g., the head official of a village) [31].The reported statistics are aggregated at each upper level before being integrated with data acquired using probability based survey methods and published by the regional bureau of statistics [22,31].
Inaccurate national agricultural statistics such as inflated crop yield resulting from underreported sown or cultivated area had been reported by previous studies using these data [32].The quality of the official statistics has been improved after probability-based survey methods were integrated with the bottom-up reporting system.Numerous peer-reviewed articles on the economic development of China that employ the official statistics have been published in professional journals [33].
We examined the effects of the vegetation conservation programs on vegetation increases in 20 counties located in northern Shaanxi Province and western Inner Mongolia.The 20 counties were chosen based on: (1) the availability of vegetation conservation data; and (2) observed NDVI increases not explained by land cover change in natural vegetation accounted for at least 80% of the total natural vegetation within the county.
We identified 20 counties in Shaanxi Province and western Inner Mongolia (Figure 1).SMK derived significant increases in precipitation coincided with more than 40% of the unexplained vegetation change in eight of the 20 counties.No more than 15% coincidence between SMK derived increases in precipitation and NDVI were found in the remaining 12 counties.Significant increases in temperature were not identified in any of the selected 20 counties.
Vegetation conservation data are available for the regions of Yan'an and Yulin in Shaanxi Province (2000 to 2008 for Yan'an; 2000 to 2007 for Yulin) and for the region of Ordos in Inner Mongolia from 2000 to 2007 [26][27][28].These data document the annual area of vegetation conservation activity including artificial planting of forest and grassland, aerial seeding, cropland converted to forest or grassland, and mountain closures (Table 2).We employed spearman rank partial correlation analysis to understand if vegetation conservation activities documented in statistical yearbooks in combination with weather variables (precipitation and temperature) can be used to explain the observed positive vegetation change.The dependent variable was defined as the spatially averaged annual maximum NDVI of natural vegetation pixels with 80% observed vegetation increases not explained by land cover change.Three independent variables were used in the analysis: (1) accumulated vegetation conservation activities from 2000 to the previous year (e.g., accumulated conservation activities for the year 2004 were the sum of conservation activities conducted from 2000 through 2003); (2) spatially averaged annual precipitation for each year; and (3) spatially averaged mean annual temperature for each year.The partial correlation analysis was carried out separately for each of the three independent variables while the other two independent variables were controlled.A two-tailed Student's test was carried out to test if the partial correlations were significant.We examined the partial correlation between these variables from 2003 through 2009.We did not include the years between 2000 and 2002 in the analysis as newly afforested/mountain closure/cropland conversion areas may not exhibit increases in NDVI immediately.We compared the results for the three independent variables to see which combination could explain most of the observed vegetation change.

Elevation and Slope
The goals of GTGP and NFCP have led us to hypothesize that the observed land cover change including conversion from cropland or sparsely vegetated area to natural vegetation occurring on steep slopes (≥15° in the Loess Plateau [8]) and in relatively high elevation regions (above 1,200 m) could be the result of the implementation of one of the two vegetation conservation programs.In order to determine whether the observed vegetation changes were indeed related to GTGP and NFCP, we compared the spatially averaged slope and elevation (derived in ArcGIS) on which the land-cover change occurred to that of the cropland and sparsely vegetated area experiencing no dynamics.We retrieved the elevation data from the ASTER Global Digital Elevation Map (GDEM) which has a spatial resolution of 1 arc-second (~30 m) and is comprised of 22,600 1° × 1° lat/lon tiles with a geographic projection [25].We generated a mosaic of 89 tiles to cover the Loess Plateau.We re-projected and resampled the ASTER DEM mosaic to Albers Equal Area Projection with 500 m spatial resolution to match the projection and spatial resolution of other datasets used in this study.

Crop Yield Data
We investigated the relationship between observed crop yield and observed NDVI changes for croplands that experienced no land cover change during the study period.We only investigated counties with available crop yield data and in which the area of croplands showing significant increases in NDVI accounted for more than 60% of the total cropland area of the county.We found nine qualified counties in the Loess Plateau.Five of these counties were in Shaanxi Province, three in Gansu Province and one in Henan Province.Significant increases in temperature and precipitation detected by the SMK analysis were not found in any of the nine counties.Dominant anthromes within the selected nine counties are rainfed villages and residential rainfed croplands according to the anthromes in 2000 (Figure 1).
Agricultural productivity data were obtained from officially published regional statistical yearbooks [22][23][24].Agricultural productivity data are available for the region of Xianyang in Shaanxi Province, the region of Qingyang in the Gansu province, and for the region of Sanmenxia in the Henan Province from 2000 to 2008 (except in 2002 and 2005 for Sanmenxia and Xianyang, respectively).Data for annual sown area, crop production and crop yield are available in the agricultural productivity dataset for major crops including: wheat, corn, bean, potato, rapeseed, tobacco, vegetables, and apple.
We calculated the sown area weighted crop yield for each county as the sum of crop yields of major crops, which were weighted by the sown area ratio.The sown area ratio of a specific crop was computed by dividing the sown area of the crop by the total sown area of all crops in a county.
We employed spearman rank partial correlation analysis to understand if crop yield change or weather change or the combination of the two factors were responsible for the observed positive vegetation changes that were not attributable to land cover change.We used spatially averaged annual accumulated NDVI of cropland pixels in a specific county as the dependent variable.Three independent variables were used in the analysis: (1) sown area weighted crop yield (kg/ha); (2) spatially averaged annual precipitation (mm); and (3) spatially averaged mean annual temperature (°C).The partial correlation analysis was carried out separately for each of the three independent variables while the other two independent variables were controlled.A two-tailed Student's t test was carried to test if the partial correlations were significant.We compared the analysis results for the three independent variables to see which combination can explain the most of the observed vegetation change.

Vegetation Change since 2000
We found significant vegetation increases as measured by NDVI in 48.2% (3.02 × 10 5 km 2 ) of the Loess Plateau (Figure 2).The area with positive vegetation change spreads across the Loess Plateau from northeast to southwest with the major part in the north of the Shaanxi province and south of the Gansu and Qinghai provinces.We found significant negative vegetation change in 0.8% (4.9 × 10 3 km 2 ) of the Loess Plateau (Figure 2) and no significant vegetation change in 51% of the study area (3.20 × 10 5 km 2 ).We excluded the region showing no significant vegetation change from the remaining analyses.

Impact of Weather Change on Vegetation Change in the Loess Plateau
We found significant positive correlation between precipitation and lagged NDVI in 82.6% of the Loess Plateau and between NDVI and temperature in 95.2% of the Loess Plateau.Precipitation is the main factor constraining vegetation growth in the far north (mean annual precipitation: 300 mm [18]) but not in the far south, which receives more precipitation (mean annual precipitation: 700 mm [18]).The central area, with moderate precipitation, revealed significant correlation between NDVI and precipitation.Despite the positive correlation between weather and NDVI, the NDVI increases cannot be entirely explained by weather changes.Significantly increasing precipitation (p ≤ 0.1) was only found in the northeastern part, and the far southwestern part of the Loess Plateau, and accounted for just 22.2% of the total area of positive vegetation change (Figure 3).We found significant warming (p ≤ 0.1) in the southwest corner of the Loess Plateau but this warming coincides with positive vegetation change in just 3.8% of the pixels (Figure 3).Based on the error assessments using p-values, the NDVI increases that can be attributed to significant increases in precipitation and temperature were 19.8% to 22.2% and 3.4% to 3.8%, respectively.

Land Cover Change Revealed by MODIS
We compared the land cover classification quality between stable and random classification pixels derived from our three-year grouping strategy during 2003 to 2005 and 2006 to 2008.Results of the two-sample t test show that the mean land cover classification quality for stable pixels was significantly higher than that of random classification pixels during both periods of 2003-2005 (70.9 vs. 61.7,p = 0.01), and 2006-2008 (72.0 vs. 63.1, p = 0.01).The higher land cover classification quality associated with the stable pixels suggests that the three-year average strategy was able to identify high quality classification results and thus can ensure reliable land cover change detection.
A little more than 13% (3.97 × 10 4 km 2 ) of the area with positive vegetation change in the MODIS time series reveals a change in terms of the three broad land cover categories of natural vegetation, cropland, and barren or sparsely vegetated area (Figure 4).Three types of land cover changes can be used to explain NDVI increases: from cropland or sparsely vegetated area to natural vegetation or from sparsely vegetated area to cropland.These three types of land cover changes accounted for 7.6% of the area with significant increases in NDVI.The land cover classification errors for these three types of land cover changes were 53% (from cropland to natural vegetation), 45% (from sparsely vegetated area to natural vegetation), and 50% (from sparsely vegetated area to cropland), respectively.Based on the probability of land cover change analysis errors, the percentage of NDVI increases that were attributable to land cover changes ranged from 3.6% to 7.6%.The remainder of the area with positive vegetation change reveals no land cover change observable with 500 m satellite data.We found that within the area with positive vegetation change based on the MODIS time series analysis, croplands are decreasing while natural vegetation and sparsely vegetated areas are increasing (Table 3).The area identified as random classification accounts for less than 3% of the total area with positive vegetation change.Therefore, it is reasonable to assume that the random classification would not cause significant bias in the land cover change analysis.The conversion between natural vegetation and cropland accounted for the majority of land cover changes detected by the MODIS land cover dataset.This conversion mainly occurred in central Shaanxi and northern Shanxi Province.Although there was a 6.2% increase in the area of barren land, increases in barren land only accounted for about 1% of the area showing land cover change.We found that the distribution of different types of land cover change were population dependent (Table 4).For example, the majority of pixels (58.6% = 1.7% + 49.3% + 7.6%) identified as croplands which were converted to natural vegetation, occurred in anthromes classified as residential anthromes while 37.2% (4.1% + 29.7% + 3.4%) of conversion of croplands to natural vegetation were found in anthromes with higher population density.In contrast, more unchanged croplands can be found in village anthromes (60.5%) compared to residential anthromes (37.1%).The conversion of barren or sparsely vegetated area to natural vegetation that occurred in village, residential and populated anthromes accounted for 11.7%, 47.1%, and 41.2% of the total conversion, respectively.
To summarize, anthromes with higher population densities revealed less cropland conversion than anthromes with lower population densities.The difference in the distribution of cropland conversion and unchanged cropland among anthromes can be explained by the concern for food security by policymakers of GTGP [34].Declining land productivity resulting from severe soil erosion in the Loess Plateau endangers food security in this region [16].The implementation of GTGP was predicted to result in a 10% to 14% loss in grain production with higher losses predicted for some provinces within the Loess Plateau (e.g., 14% to 17% for the Shaanxi Province) [35].In order to achieve the dual goal of vegetation conservation and minimizing grain production loss, GTGP was designed to focus more on marginal croplands with low productivity in remote and mountainous regions hosting lower population density (e.g., croplands in residential anthromes) [7,36,37].By examining the characteristics of 3,397 parcels of land in Shaanxi Province, Kelley and Huo [37] reported that land parcels with steeper slope, lower yield and greater distance from farmers' houses had a significantly higher probability of being enrolled in GTGP.Xu et al. [34] reported that the yield of cropland plots which were transformed under GTGP was about 30% to 50% of the yield by non-GTGP cropland plots based on data acquired by a national survey of GTGP.The targeting strategies of GTGP reported by previous studies agree with what we found by examining the distribution of land cover changes over anthromes: GTGP left croplands with higher productivity in areas with higher population density (e.g., croplands in village anthromes) unchanged to offset the decrease in grain production due to loss of croplands elsewhere.
Conversion of sparsely vegetated areas occurred mainly in anthromes with higher population densities, e.g., residential and populated rangelands instead of remote rangelands.We suspect that since afforestation requires manual labor, a certain population density is necessary, and, thus, more re-vegetated barren lands were found in anthromes with higher population densities [36].

Vegetation Conversion Reported in Statistical Yearbooks
We found that for 68.7% of the area with significant vegetation changes, we could not explain the changes directly by MODIS land cover changes or weather changes.In order to identify conservation activities that may be missing in the MODIS land cover change analysis, we examined the statistical yearbook data for 20 counties located in Northern Shaanxi Province and Western Inner Mongolia.Positive vegetation change within the 20 counties accounts for approximately 21% of the total positive vegetation changes within the Loess Plateau.
We found that accumulated vegetation conservation activities, mean annual temperature and annual precipitation were significantly correlated with maximum annual NDVI (P < 0.05) within five, four, and two of the 20 selected counties, respectively.The distribution of the partial correlation coefficients and p-values for the three independent variables are shown in Figure 5. .Partial correlation analysis results using vegetation conservation data acquired from regional statistical yearbook."Con:r" and "Con:p" refer to the correlation coefficients and p values of the partial correlation between accumulated vegetation conservation and NDVI, with precipitation and temperature as control variables."Precip:r" and "Precip:p" refer to the correlation coefficients and p values of the partial correlation between NDVI and precipitation, with conservation and temperature as control variables."Temp:r" and "Temp:p" represent the correlation coefficients and p values of the partial correlation between NDVI and temperature, with annual precipitation and accumulated conservation as control variables.
There was no significant weather change derived from the SMK trend analysis in 12 of the 20 selected counties.By employing the partial correlation analysis, however, we identified three counties (Yan'an, Yanchuan, and Zhidan) where conservation activities, annual precipitation and mean annual temperature were all found to be significantly correlated with maximum annual NDVI (P < 0.05), when the other two variables were controlled.This demonstrated the importance of using multi-variable analysis in change attribution.
In the remaining eight counties where significant increases in precipitation was detected by SMK trend test, conservation activities were found to be significantly correlated with annual maximum NDVI (P < 0.05) in Shenmu and Zichang.In Zichang county, besides conservation activities, mean annual temperature was also identified to be significantly correlated with annual maximum NDVI.
Conservation program related vegetation increases within the selected counties have been identified by using 30 m Landsat images or field-measured data acquired at very fine scale in previous studies.Zhou et al. [38] investigated the effects of GTGP on land cover changes in Ansai county, Yan'an between 1995 and 2010 using multi-temporal Landsat TM/ETM+ images.They identified a 21.4% increase in newly forested area by 2010 [38].The increases in newly forested area were attributed to conversions of cropland and shrub-grassland to forests by GTGP [38].Cao et al. [12] investigated the effects of afforestation, cropland conversion, and grazing prohibition by GTGP on vegetation dynamics in five counties within the regions of Yan'an and Yulin by sampling vegetation changes in 150 field plots with each plot having a size of 0.5 ha.They reported a 12.5% increase in the average vegetation cover of the five counties between 1998 and 2005 due to GTGP.The unexplained positive NDVI increases in natural vegetation area of the 20 selected counties account for 19.9% of the total area with NDVI increases in the Loess Plateau.We suspect that the vegetation increases within the selected 20 counties were due to vegetation conservation activities that occurred at scales finer than 500 m, thus land cover changes resulting from these conservation activities were missed in the MODIS land cover change analysis.

The Effect of GTGP on Areas with Steep Slopes and Sparsely Vegetated Areas
The goal of GTGP is to convert croplands on steep slopes to forest and grassland and afforest barren and sparsely vegetated areas.Steep slopes are defined as slopes greater or equal to 15° in northwestern China [8] where the major part of the Loess Plateau is located.Our land cover change analyses based on MODIS have confirmed a decrease in the area of cropland and an increase in the area of natural vegetation over the study period (Table 3).We found that 41.6% (8,977.50km 2 ) of the area that changed from cropland to natural vegetation was located in areas with slopes steeper than 15°.The spatially averaged slope on which croplands were converted to natural vegetation is 14.3° (Table 5), which is significantly steeper than the average slope of unchanged croplands (10.9°) (P < 0.01).This may indicate the possible effect of GTGP.The conversion from sparsely vegetated area to natural vegetation was mainly detected in the region of Ordos, which lies between the Hobq Desert and the Mu Us Desert.We speculate that the vegetation increases in this region were due to vegetation conservation activities by GTGP.

The Effect of Conservation Programs on Unchanged Croplands
The pixels identified as unchanged croplands that did not experience significant weather change but did reveal significant NDVI increases accounts for 12.2% of the total area with vegetation increases.We found significant positive relationships between sown area weighted crop yield and annual accumulated NDVI for these areas.Thus, it appears that crop yields are increasing in unchanged croplands in the Loess Plateau.Due to limited availability of crop yield data, we have not obtained crop yield data covering all the unchanged croplands area.
Our results in Figure 6 reveal a strong correlation between increased crop yields and increases in accumulated NDVI in seven of the selected nine counties (P < 0.05), when weather variables were controlled.Croplands with observed NDVI increases in the seven counties accounted for 2.1% of the positive vegetation change area.We did not identify any county where annual precipitation or mean temperature was significantly correlated with accumulated NDVI. Figure 6.Partial correlation analysis using crop yield data acquired from regional statistical yearbook."Yield:r" and "Yield:p" refer to the correlation coefficients and p values of the partial correlation between sown area weighted crop yield and accumulated NDVI, with precipitation and temperature as control variables."Precip:r" and "Precip:p" refer to the correlation coefficients and p values of the partial correlation between NDVI and precipitation, with crop yield and temperature as control variables.'Temp:r' and "Temp:p" represent the correlation coefficients and p values of the partial correlation between NDVI and temperature, with annual precipitation and sown area weighted crop yield as control variables As shown in Tables 4 and 5, the spatially averaged slope in unchanged croplands was significantly lower than the slope of croplands converted to natural vegetation, and the majority of unchanged croplands (60.5%) are found in anthromes with high population densities (e.g., irrigated and rainfed villages).We suggest that these croplands were not participating in GTGP.They were more profitable and located in areas with suitable topography.Positive vegetation change in these croplands may be caused by significant increases in crop yield due to intensified agricultural production to offset economic loss caused by GTGP.Increases in crop yield were achieved through intensified agricultural production on non-GTGP croplands by means such as utilizing improved seeds and increasing crop diversity [34].For example, Xu et al. [34] found increases in crop yield in Delong county of Ningxia Hui Autonomous Region located in the western Loess Plateau from about 2,516 kg/ha to about 4,312 kg/ha after GTGP was initiated in this county.We did not find any articles or other relevant data about the contribution of fertilizer and technology to the increased crop yield in this region.However, it seems reasonable to speculate that increased utilization of fertilizer and advances in agricultural production technology were also potential causes of crop yield increases.

Conclusion
Change attribution is one of the most challenging parts of any change analysis based on satellite data.The observed NDVI increases that we were able to explain accounted for about 34.7% to 37.5% of the total changes observed in the NDVI time series analysis.We were able to attribute about 23.2% to 26.0% of the changes to precipitation and temperature increases.We were able to attribute 3.6% to 7.6% of the NDVI increases to large-scale land cover changes caused by GTGP and NFCP programs.We also attributed 5.8% (1.7% + 1.8% + 2.3%) of the observed NDVI increases to the interacting effects of changes in weather variables and conservation programs.Based on the available county level crop yield data, we were able to identified 2.1% the vegetation increase was a result of crop production intensifications in unchanged cropland areas (Table 6).
Due to the limited availability of crop yield and vegetation conservation data, we used a stepwise approach for the attribution of observed vegetation increases.The lack of data available may have resulted in the omission of some change attribution interaction factors.Approximately 62.5% of the observed NDVI increases were left unexplained, most of which were concentrated in natural vegetation and cropland regions within the Loess Plateau.We speculate that part of these unexplained changes may be the result of fine scale land cover change or changes in agricultural management practices resulting from the large scale conservation programs.The impact of vegetation conservation programs on vegetation change in the Loess Plateau is twofold.On the one hand, vegetation conservation programs target marginal lands (e.g., steep slope croplands or barren lands in remote areas), which are less profitable and more difficult to manage.Thus more vegetation increases due to cropland conversion and afforestation can be found in these regions.On the other hand, intensified agricultural production can be found in non-GTGP/NFCP plots with suitable topography and well-established irrigation systems to offset the agricultural production loss caused by vegetation conservation programs elsewhere.
We expect the results of this study to facilitate the further development of forest transition theory in two ways: (1) the study of the targeting strategy of large scale conservation programs will provide useful knowledge in terms of understanding how policy change would affect the spatial adjustment of agriculture, which was viewed as a key process that would allow forest regrowth [1]; (2) The impacts of displaced domestic crop and forestry product demands by China and Vietnam caused by conservation programs in these countries are a major concern in international literatures [39,40].The study of how conservation programs in China affected crop yield in non-enrolled cropland area may help future studies to assess the possible amounts of displaced crop products demands from China.

Figure 1 .
Figure 1.Upper left inset: Location of the Loess Plateau and the 138 weather stations in China.Main map frame: Location of key case study counties within the Loess Plateau.Case studies counties selected for conservation program analyses are outlined in green (County #1 to #20).Case studies counties selected for crop yield analyses are outlined in black (County #21 to #29).The spatial distribution of anthropogenic biomes is based on the anthropogenic biomes in 2000 as described in Ellis et al. [15].

Figure 2 .
Figure 2. Significant vegetation change between 2000 and 2009 in the Loess Plateau.The significant vegetation change was identified by applying seasonal Mann-Kendall trend test to 500 m MODIS NDVI time series acquired between 2000 and 2009.Green pixels represent significant positive vegetation change, which accounted for 48.2% of the Loess Plateau.Brown pixels indicate significant negative vegetation change, which accounted for 0.8% of the Loess Plateau.Pixels in white represent either areas with no significant vegetation change or areas with high missing data rates due to frequent cloud cover.

Figure 3 .
Figure 3. Significant temperature and precipitation increases from 1999 to 2009 in the Loess Plateau.The significant temperature and precipitation increases were identified by applying seasonal Mann-Kendall trend test to monthly 500 m temperature and precipitation time series interpolated from weather station based monthly observations.Blue pixels represent significant precipitation increases, which accounted for 16.8% of the Loess Plateau.Dark blue pixels indicate the coincidence of significant NDVI and precipitation increases, which accounted for 10.7% and 22.2% of the entire Loess Plateau and positive vegetation change area, respectively.Brown pixels indicate significant temperature increases, which accounted for 4.7% of the Loess Plateau.Dark brown pixels indicate the coincidence of significant NDVI and precipitation increases, which accounted for 1.8% and 3.8% of the entire Loess Plateau and positive vegetation change area, respectively.

Figure 4 .
Figure 4. Land cover change within the positive vegetation change region of the Loess Plateau between 2003-2005 and 2006-2008.The land cover change analysis was conducted among three broad land cover categories: natural vegetation, cropland and barren or sparsely vegetated area.Between 2003-2005 and 2006-2008, there was 4.1% and 6.2% increase in natural vegetation and sparsely vegetated area within the positive vegetation change region, respectively.A 10.6% decrease in cropland area was also identified.

Figure 5
Figure5.Partial correlation analysis results using vegetation conservation data acquired from regional statistical yearbook."Con:r" and "Con:p" refer to the correlation coefficients and p values of the partial correlation between accumulated vegetation conservation and NDVI, with precipitation and temperature as control variables."Precip:r" and "Precip:p" refer to the correlation coefficients and p values of the partial correlation between NDVI and precipitation, with conservation and temperature as control variables."Temp:r" and "Temp:p" represent the correlation coefficients and p values of the partial correlation between NDVI and temperature, with annual precipitation and accumulated conservation as control variables.

Table 1 .
Overview of data sources.

Table 3 .
Land cover change shown by three-year-summarized MODIS land cover data.

Table 4 .
Distribution of land cover changes in major anthromes within the area of positive vegetation change.

Table 5 .
Slope and elevation results.
* Overlapping effects of different causes were removed.