Did Ecological Restoration Hit Its Mark ? Monitoring and Assessing Ecological Changes in the Grain for Green Program Region Using Multi-source Satellite Images

Ecological restoration programs are expected to control environmental deterioration and enhance ecosystem functions under a scenario of increasing human disturbance. The largest ecological restoration program ever implemented in China, the first round of the countrywide Grain for Green Program (GGP), finished in 2010. However, it is not known whether the ecological changes that resulted from the GGP met the restoration goal across the whole implementation region. In this study, we monitored and assessed the ecological changes in the whole GGP region in China over the lifetime of the first round of implementation (2000–2010), by establishing a comprehensive assessment indicator system composed of ecosystem pattern, ecosystem quality (EQ), and key ecosystem services (ESs). Remote sensing interpretation, ecological model simulations based on multi-source images, and trend analysis were used to generate land use and land cover (LULC) datasets and estimate ES and ESs indicators. Results showed that while forest increased by 0.77%, artificial land increased more intensely by 22.38%, and cropland and grassland decreased by 1.81% and 0.68%, respectively. The interconversion of cropland and forest played a primary role in ecosystem pattern change. The increase in ecosystem quality measures, including fractional vegetation cover (0.1459% yr−1), leaf area index (0.0121 yr−1), and net primary productivity (2.6958 gC m−2 yr−1), and the mitigation of ecosystem services deterioration in soil water loss (−0.0841 t ha yr−1) and soil wind loss (−1.0071 t ha yr−1) in the GGP region, indicated the positive ecological change in the GGP region to some extent, while southern GGP subregions improved more than the those in the north on the whole. The GGP implementation other than climate change impacted ecological change, with contributions of 14.23%, 9.94%, 8.23%, 30.45%, and 18.05% in the ecological outputs mentioned above, respectively. However, the water regulation did not improve (−2283 t km−2 yr−1), revealing trade-offs between ecosystem services and inappropriate afforestation in ecological restoration programs. Future GGP implementation should change the practice of large-scale afforestation, and focus more on the restoration of existing forest and cultivation of young plantings, formulating rational and specific plans and designs for afforestation areas through the establishment of near-natural vegetation communities, instead of single-species plantations, guided by regional climate and geographical characteristics. Remote Sens. 2019, 11, 358; doi:10.3390/rs11030358 www.mdpi.com/journal/remotesensing Remote Sens. 2019, 11, 358 2 of 28


Introduction
Humans have displayed a history of deforesting and changing land use and land cover since farming began.China's environmental problems are among the most severe of any major country, because of the need for food security for its large population and the pursuit of high economic growth [1].From the 1950s to the 1980s, China experienced severe losses of forest and grassland because of large-scale deforestation and grassland reclamation for crop production and industrial construction [2][3][4][5].This pattern of changing land use has led to severe soil erosion [6]; land desertification [7,8]; loss of biodiversity [9,10]; drying rivers and shrinking lakes [11]; and the frequent occurrence of dust storms, floods, and drought [1,12,13].
To address these problems and maintain environmental sustainability, several countrywide ecological restoration programs have been implemented [2,14].The Grain for Green Program (GGP), also known as the Sloping Land Conversion Program (SLCP), is regarded as the largest ecological restoration program ever implemented in a developing country [4,15].The GGP, originally planned for 2001-2010 (the first round), aimed to set aside 14.67 million ha of cropland (4.4 million ha with slopes greater than 25 • ) and afforest a roughly equal area of wasteland by 2010 [16].By the end of 2015, the Chinese government had achieved its ambition to increase forest by approximately 29.87 million ha, with a total investment of approximately 405.66 billion RMB (US$58.52 billion) [17].A new round was launched at the end of 2014 [18], and the implementation region was expanded at the end of 2015 [19].An assessment of the impact of the first round of GGP on ecosystems could provide insights into future ecological restoration, and provide input into ecological programs initiated by other countries' governments all around the world.
Previous studies have evaluated the ecological changes brought about by the GGP restoration program from different perspectives, including land use and land cover (LULC) change, vegetation condition, net primary productivity (NPP), and ecosystem services, since the program was designed to convert farmland to forest and grassland, and further reduce soil erosion and desertification [20][21][22][23].Zhu & Wang [24] detected land cover change in Ansai County and found that farmland had been reduced overall from 18.12% to 13.58% in 2004-2010, while forest increased from 14.60% to 37.95%.Li et al. [25] investigated the vegetation in Bijie Prefecture and found that it had been significantly improved, but the trend of degradation had not been reversed in the east of the region.Feng et al. [26] reported that the NPP and evapotranspiration (ET) in Loess Plateau had increased by 9.3 gCm −2 and 4.3 mm per year, respectively, in 2000-2010, but lost soil moisture of 2.4 mm yr −1 and runoff of 0.5 mm yr −1 simultaneously.Wang et al. [19] stated that the transformation from non-irrigated farmland to forestland could potentially improve soil conservation by 24.89%, but the increase extent of GGP implementation decreased NPP and water yield at the sub-watershed scale.With the analysis of different ecological variables, these findings underscore the complexity and trade-offs behind the anthropogenic program implementation.Although previous studies have investigated eco-environmental change in the GGP region, most have focused on the local or regional scale, rather than the national scale.Moreover, the impact of climate change on ecosystems has often been neglected, which may cause unforeseen errors in the assessment.
In this paper, we monitor and assess the ecological changes in the whole region where the GGP was implemented in China from 2000 to 2010 using remote sensing and modelling techniques along with trend analysis, based on a comprehensive assessment indicator system established to reflect changes in ecosystem pattern, ecosystem quality, and key ecosystem services.Our study attempts to answer the following questions: (1) did the first round of the GGP achieve its goal of enhancing local ecological and environmental conditions after investing huge amounts of capital and resources?; (2) how can these ecological changes be measured comprehensively and effectively?; (3) did the whole GGP implementation region show a notable positive trend of ecological change?; and (4) to what extent did the GGP implementation improve the ecological environment after excluding the influence of climate change?

Subregions
The GGP started with a pilot program in three provinces in 1999 and was formally launched in 2002.It expanded to a total of 1897 county units across 25 provinces (autonomous regions and municipalities) and the Xinjiang Production and Construction Corps-a wide area with huge ecological and economic heterogeneity, covering 76.77% of China-and involved tens of millions of rural households [16].Following the principle of hazard prevention, the GGP implementation area was divided into 10 subregions by China's State Forestry Administration (SFA) according to the risk of soil erosion, wind erosion, and desertification; local hydrothermal conditions; and landform characteristics: the Alpines-Gorges Region in Southwestern China (AGRSC), the Mountains-Hills Region in Sichuan, Chongqing, Hubei and Hunan Province (MHRSCHH), the Low Mountains-Hills Region in Middle and Lower Yangtze Plain (LMHRMLYP), the Yunnan-Guizhou Plateau Region (YGPR), the Mountains-Hills Region in Hainan and Guangxi Province (MHRHG), the Alpine Steppe and Meadow Region in the source of Yangtze River and Yellow River (ASMRYY), the Arid Desert Region in Xinjiang Province (ADRX), the Hills-Gullies Region in Loess Plateau (HGRLP), the Arid and Semi-arid Region in Northern China (ASRNC), and the Mountains-Sands Region in Northeastern China (MSRNC) (Figure 1).By taking different conservation measures in each subregion (Table 1), the government hoped to recover and enhance natural ecosystems more efficiently [27].
Remote Sens. 2019, 11, x FOR PEER REVIEW 3 of 27 extent did the GGP implementation improve the ecological environment after excluding the influence of climate change?

Subregions
The GGP started with a pilot program in three provinces in 1999 and was formally launched in 2002.It expanded to a total of 1897 county units across 25 provinces (autonomous regions and municipalities) and the Xinjiang Production and Construction Corps-a wide area with huge ecological and economic heterogeneity, covering 76.77% of China-and involved tens of millions of rural households [16].Following the principle of hazard prevention, the GGP implementation area was divided into 10 subregions by China's State Forestry Administration (SFA) according to the risk of soil erosion, wind erosion, and desertification; local hydrothermal conditions; and landform characteristics: the Alpines-Gorges Region in Southwestern China (AGRSC), the Mountains-Hills Region in Sichuan, Chongqing, Hubei and Hunan Province (MHRSCHH), the Low Mountains-Hills Region in Middle and Lower Yangtze Plain (LMHRMLYP), the Yunnan-Guizhou Plateau Region (YGPR), the Mountains-Hills Region in Hainan and Guangxi Province (MHRHG), the Alpine Steppe and Meadow Region in the source of Yangtze River and Yellow River (ASMRYY), the Arid Desert Region in Xinjiang Province (ADRX), the Hills-Gullies Region in Loess Plateau (HGRLP), the Arid and Semi-arid Region in Northern China (ASRNC), and the Mountains-Sands Region in Northeastern China (MSRNC) (Figure 1).By taking different conservation measures in each subregion (Table 1), the government hoped to recover and enhance natural ecosystems more efficiently [27].

Comprehensive Assessment Indicator System for Assessing Ecological Changes from the GGP
To comprehensively and efficiently monitor and assess the ecological changes in the GGP region, we developed an assessment indicator system from three perspectives: ecosystem pattern, ecosystem quality, and key ecosystem services (Table 2).Ecosystem area and ecosystem conversion area, which have been used in many studies, were selected as indicators to reflect changes in the ecosystem pattern [13,20].As the amount and distribution of vegetation are of prime importance for terrestrial ecosystems [28], the fractional vegetation cover (FVC), leaf area index (LAI), and NPP were chosen to detect ecosystem quality.The FVC is the percentage of vegetation occupying the ground area in a vertical projection [29], while LAI is the ratio of leaf surface area to unit ground surface area [30].The FVC and LAI parameterize the horizontal and vertical density of vegetation, respectively [31].They are two important variables used to describe the condition of land surface vegetation, characterize ecosystem changes, and reflect the related biological and physical processes [31][32][33][34].NPP marks the first visible step of carbon accumulation [35], is a key regulator of ecological processes [36], and is an important component of net ecosystem exchange or net ecosystem productivity.It indicates the metabolism of the ecosystem [37].Because the FVC, LAI, and NPP vary within a year, we chose the annual maximum value of FVC (FVC max ) and LAI (LAI max ) to reflect the instantaneous vegetation growth status, and the annual average value of NPP (NPP avg ) to reflect the accumulated process of carbon accumulation.In addition, owing to the aims of the GGP were to control soil erosion, mitigate flooding, and prevent sandstorms, we selected three key ecosystem services, water regulation (WR), soil conservation (SC), and sandstorm prevention (SP), to evaluate the enhancement of ecosystem services arising from the program implementation.The annual quantity per unit area of water regulation (WR pua ), soil water loss (SWAL pua ), and soil wind loss (SWIL pua ) were selected as assessment indicators to exclude the influence of area.

Generation of Land Use and Land Cover Datasets
Landsat TM/ETM images in 2000 and 2005, and HJ-1 A/B images in 2010 with a 30 m spatial resolution were used to interpret LULC change.Using the FAO Land Cover Classification System (FAO LCCS) tools [38], we classified the images into 38 LULC level II classes [39], and then categorized them into six LULC level I classes, which were forest, grassland, wetland, cropland, artificial land, and other types land, based on the IPCC classification system [40].
We derived the 2010 LULC vector dataset from 2010 HJ-1 A/B images by applying object-oriented automatic classification methods after data preprocessing on a supercomputing platform, with associated data from ground surveys and radar [41].The overall producer's accuracy for the 2010 LULC dataset was validated with a total of 31,675 independent ground survey samples based on random sampling and reached 94% and 86% accuracy on average for LULC classes at level I and level II, respectively.
The 2000 and 2005 LULC vector datasets were generated using object-oriented spectrum eigenvectors similarity change detection.The spectrum eigenvectors of a Landsat TM image patch can be regarded as a vector in a six-dimensional feature space, and the smaller the included angle and closer the norm of two vectors, the more similar two patches are.Based on this principle, we used the cosine of the angle between two vectors and the ratio of the vector norm to measure vector similarity, see Song & Yan [42] for details.In eCognition Developer 8.64, we used the 2010 LULC vector dataset generated above to overlap and segment the 2000 and 2005 Landsat TM images with 2010 HJ-1 images at multiple scales, which ensured that the unchanged areas had the same border and provided a priori information on ground objects, enhancing the classification accuracy.We extracted spectrum eigenvectors from each image patch segmented above and calculated the similarity between the images of target years (2000 and 2005) and the base year (2010) to detect LULC change, and then updated the LULC dataset accordingly [43].A detailed workflow for the LULC datasets in the target years is displayed in Figure 2.
The 2000 and 2005 LULC vector datasets were generated using object-oriented spectrum eigenvectors similarity change detection.The spectrum eigenvectors of a Landsat TM image patch can be regarded as a vector in a six-dimensional feature space, and the smaller the included angle and closer the norm of two vectors, the more similar two patches are.Based on this principle, we used the cosine of the angle between two vectors and the ratio of the vector norm to measure vector similarity, see Song & Yan [42] for details.In eCognition Developer 8.64, we used the 2010 LULC vector dataset generated above to overlap and segment the 2000 and 2005 Landsat TM images with 2010 HJ-1 images at multiple scales, which ensured that the unchanged areas had the same border and provided a priori information on ground objects, enhancing the classification accuracy.We extracted spectrum eigenvectors from each image patch segmented above and calculated the similarity between the images of target years (2000 and 2005) and the base year (2010) to detect LULC change, and then updated the LULC dataset accordingly [43].A detailed workflow for the LULC datasets in the target years is displayed in Figure 2.

Detection of Ecosystem Pattern Changes
Ecosystem types were consistent with the six LULC level I classes in this study, because we had considered the spectrum characteristics, vegetation coverage, and ecosystem plant community component characteristics in our LULC classification system [44].Thus, ecosystem pattern changes can be identified through LULC change detection.Here, we employed an LULC transfer matrix to detect the change between LULC types in 2000-2005, 2005-2010, and 2000-2010.The LULC transfer matrix delivers dynamic information on mutual transfer processes between different LULC types from the beginning to the end of a given period; thus it not only includes the static area data of different LULC types at a point in time, but also contains information about the area transferred out and in during a period [45].The general form of the LULC transfer matrix is: where i and j (i, j = 1, 2, . . ., n) are the LULC type before and after the transformation happened, respectively; n is the total number of LULC types; and s ij is the area of LULC type i converted to type j over the period, or the area of untransformed LUCC type i when i = j.Each row and column of the transfer matrix, respectively, denotes the source and destination of LULC change.
We rasterized the LULC vector datasets for 2000, 2005, and 2010 into 30 m spatial resolution grids.We calculated the LULC transfer matrix pixel by pixel between each two phases of LULC maps using spatial analyst tools in ArcGIS 10.3.

Ecosystem Quality Estimation
The FVC was estimated using the pixel dichotomy model [31], which assumes that image pixels are composed of two parts: vegetation-covered and non-vegetation-covered areas [46], and thus the FVC can be defined as follows: where NDVI veg and NDVI soil are the Normalized Difference Vegetation Index (NDVI) value for pure vegetation pixels and pure bare soil pixels, respectively.We used MOD13A3 data, monthly MODIS NDVI at a 1 km spatial resolution as a gridded level-3 product, as the calculation input after aggregation through the Maximum Value Composite (MVP) approach [47] to obtain the annual maximum value and after time domain processing using a Savitzky-Golay filter [48] to minimize the influence of existing noise.
The LAI was derived from the MODIS MOD15A2 product, which is composited every eight days at a 1 km resolution.We improved our LAI data quality using the temporal spatial filter (TSF) [49] to fill in the spatial and temporal discontinuity of original datasets.We then used MVP to obtain the annual maximum value of LAI.
NPP was evaluated by applying the GLOPEM-CEVSA model [50], which couples a satellite-based model and a process-based model.The model has been validated by carbon flux observations within forest, grassland, and cropland.In the model, NPP is the difference of gross primary productivity (GPP) and autotrophic respiration (R a ), where GPP is estimated by the GLoPEM model, and R a is divided into maintenance respiration (R m ) and growth respiration (R g ) from four plant parts-thick root, thin root, stem, and leaf-with reference to the process-based BEPS model and Biome-BGC model: where PAR (photosynthetically active radiation) was estimated using climatology radiation calculation methods [51] and FPAR (fraction of photosynthetically active radiation) was adapted from satellite remote sensing inversion data based on the canopy radiative transfer model [52].ε g , the vegetation light utilization efficiency in the GPP concept, was the product of the maximum possible light utilization efficiency and environment stresses including vapor pressure deficit (VPD), PAR, air temperature, atmospheric CO 2 , and soil moisture [53].M i and r m,i are the biomass and maintenance respiration coefficients of the four plant organs, respectively.Q 10 represents the change in the maintenance respiration rate when the temperature changes by 10 degrees; T and T b are the actual temperature and base temperature, respectively; and r g is the proportional coefficient of total growth respiration to total plant growth.The input meteorological data, including temperature, precipitation, evaporation, and total radiation, were provided by China's State Meteorological Administration.They were interpolated into grids using thin plate smoothing splines (ANUSPLIN) [54,55].

Estimation of Ecosystem Services
WR refers to the magnitude of runoff, flood, and aquifer recharge influenced by changes in ecosystems [56,57].We estimated WR using the Rainfall Storage Method (RSM), which considers the hydrological regulation effect of ecosystems, including forest, grassland, and cropland, compared with bare land [58,59].RSM can be expressed as follows: where A is the area of the ecosystem; J 0 is the annual average rainfall; K is the ratio of rainfall yielding runoff to total rainfall, which is determined by precipitation and raininess; and R 0 and R g are the runoff yield rates of bare land and the ecosystem.respectively.The thresholds of rainfall to yield runoff at different locations were collected from published literature, and the single rainfall values at the same locations were derived from three-hour TRMM precipitation data (3B42), which were revised by the measured daily precipitation data of neighboring national meteorological stations.By accumulating single rainfall values greater than the corresponding threshold, we obtained the K values at these single points.With R 2 = 0.79, a good linear relationship existed between these K values and annual average runoff coefficient at the same locations, thus a national spatial distribution of K value was derived from the national annual average runoff coefficient [60].The R g of forest was collected from relevant articles using meta-analysis according to factors such as vegetation types, geographic location, climate zone, and measured surface runoff.The R g of grassland was calculated according to grassland coverage [61], except for alpine meadows, which adopted the study result of Li et al. [62].The water conservation amount of floodplain wetlands and inland marshes was calculated at 8100 m 3 per hectare according to Meng et al. [63].SC is the reduction of soil and water loss (SWAL)-the soil erosion resulting from raindrop splash and runoff [64]-by ecosystems through their structures and processes [65,66].It is the difference between potential and actual soil erosion [19,57].Therefore, a decline in SWAL implies the improvement of SC.Using the Revised Universal Soil Loss Equation (RUSLE) [64], SWAL can be expressed as: where R is the rainfall-runoff erosivity factor, K is the soil-erodibility factor, LS is the topographic factor, C is the cover-management factor, and P is the support practice factor.In this study, R was calculated using the Daily Rainfall Erosivity Model [67] based on the daily precipitation from national meteorological stations, K was obtained using the soil-erodibility nomograph [68] with the soil dataset from the Second National Soil Survey of China, C was derived according to its relationship with NDVI [69], P was obtained by reference to the USDA Handbook Number 703 and Cai et al. [69], and the LS was computed using the V4.1 SRTM (Shuttle Radar Topography Mission) DEM product at a 90-m spatial resolution by the integrated method of Sun et al. [70].
Remote Sens. 2019, 11, 358 10 of 28 SP is one of the most important services supplied by ecosystems in arid and semi-arid regions, where vegetation measures help to control desertification and prevent dune movement by fixing dunes and slowing wind velocities [71,72].The Revised Wind Erosion Equation (RWEQ) is a combination of empirical and process-based models and has been extensively tested in simulating soil loss eroded by wind (SWIL) under broad field conditions [57,73].As SP is the difference between potential and actual SWIL, the less SWIL, the more SP.SWIL can be expressed as: where Q max is the actual maximum transport capacity, s is the actual critical field length, WF is the weather factor, EF is the soil erodible fraction, SCF is the soil crust factor, K is the soil surface roughness, and COG is the combined residue factor.For data sources and processing details, see Gong et al. [74].

Trend Analysis
The annual trends for ecological changes from 2000 to 2010 were calculated pixel by pixel using the slope equation fitted with ordinary least-squares regression in ENVI 5.4.The slope indicated strong or weak and positive or negative change for the variables [70], with boundaries defined from a 95% confidence interval, which is expressed as follows: where i is the number of years, j is the number of raster pixels, n is the length of study years, and var i,j stands for the value of the variable (including EQ, ES, and climate indicators) at pixel j in year i.

Changes in Ecosystem Pattern
During the first round of GGP from 2000 to 2010, the ecosystem composition and distribution in the GGP area (Table 3) indicated that the forest ecosystem was the largest ecosystem in the whole GGP area throughout the study period, occupying above 30%, followed by grassland, other types of land, cropland wetland, and artificial land.From 2000 to 2010, forest, wetland, and artificial land showed an increasing trend, while grassland, cropland, and other types of land declined.The artificial land increased dramatically by 22.38%; forest and wetland increased slightly by 0.77% and 1.03%, respectively; and cropland, grassland, and cropland decreased by 1.81%, 0.68%, and 0.41%, respectively.Forest, cropland, and other types of land changed more acutely in the first five years (2000-2005), while grassland, wetland, and artificial land changed more dramatically in the last five years (2005)(2006)(2007)(2008)(2009)(2010).Forest in the GGP area was mainly distributed in the southern and northeastern subregions, while grassland was largely distributed over the northern subregions.Cropland was concentrated in the eastern plains and the Sichuan Basin (southwestern MHRSCHH), while the other types of land class were concentrated in the northwest Gobi desert throughout most of the ADRX and ASMRYY, and the western ASRNC.Wetland and artificial land were scattered across the whole area and were too small to be observed on the map (Figure 3a-c).

Spatial Pattern and Temporal Trend of Ecosystem Quality
The terrestrial ecosystems in the GGP region from 2000 to 2010 had an average of 57.33%, 1.91, and 438.65 gC m -2 in annual FVCmax, LAImax, and NPPavg, respectively (Table 5).The spatial patterns of EQ distribution and change are shown in Figure 4, and varied distinctly across the whole restoration area.Higher FVCmax, LAImax, and NPPavg values were distributed in the MHRHG, YGPR, MHRSCHH, LMHRMLYP, and MSRNC (roughly concentrated between 84%-94%, 3-4, and 700-1000 gC m -2 , respectively, all greater than the corresponding annual average value for the whole GGP region), over southern and northeastern China, where most forests are located.Northwestern China showed poor EQ, especially the HGRLP, ASRNC, ASMRYY, and ADRX, and the EQ and vegetated land descended in turn (Figure 4a-c, Figure 5).The spatial changes of annual FVCmax and LAImax from 2000 to 2010 basically followed the same spatial pattern as distributions-gradually decreasing from southeast to northwest overall-but differed in some areas.For example, the eastern HGRLP, central and southwest of MSRNC, where lower FVCmax and LAImax values were concentrated, had positive changes; the northwestern ADRX, where higher FVCmax and LAImax values were found, however, was getting worse (Figure 4a,b,d,e).In relation to annual NPPavg change, the slope declined from the center The LULC transfer matrix (Table 4) demonstrated that the interconversion of cropland and forest played a primary role in ecosystem pattern change during 2000-2010, with 2.45% of cropland converted to forest and 1.11% of forest transferred to cropland.In addition, the transformations from cropland to artificial land, from grassland to forest, from cropland to grassland, and from grassland to cropland were also significant.The ecosystem transformation mainly took place in 2005-2010, since the transformed area of each ecosystem type was almost twice as much as that from 2000 to 2005.Forest, grassland, and cropland ecosystems were the major transfer direction and sources to each other, while wetland was mainly converted to cropland and grassland and was also mainly replenished from cropland, grassland, and other types of land.Artificial land mostly came from and transferred into cropland.Among all ecosystems, the dynamics of wetland were the most drastic: only 96.54%, 95.Wetland primarily transferred into cropland, which ranked second in ecosystem change intensity (94.10% unchanged in 2000-2010), followed by artificial land (95.43% unchanged), grassland (97.44% unchanged), forest (98.13% unchanged), and other types of land (99.19% unchanged).The ecosystem transformation to forest was scattered throughout the whole GGP, but concentrated in ADRX; while the ecosystem transformation to grassland mainly occurred in the northwestern GGP regions, including western HGRLP, ASMRYY, ADRX, and AGRSC.Artificial land mainly increased in the eastern GGP regions (Figure 3d).Note: Each row denotes the LULC transfer source, and each column denotes the LULC transfer destination.

Spatial Pattern and Temporal Trend of Ecosystem Quality
The terrestrial ecosystems in the GGP region from 2000 to 2010 had an average of 57.33%, 1.91, and 438.65 gC m −2 in annual FVC max , LAI max , and NPP avg , respectively (Table 5).The spatial patterns of EQ distribution and change are shown in Figure 4, and varied distinctly across the whole restoration area.Higher FVC max , LAI max , and NPP avg values were distributed in the MHRHG, YGPR, MHRSCHH, LMHRMLYP, and MSRNC (roughly concentrated between 84%-94%, 3-4, and 700-1000 gC m −2 , respectively, all greater than the corresponding annual average value for the whole GGP region), over southern and northeastern China, where most forests are located.Northwestern China showed poor EQ, especially the HGRLP, ASRNC, ASMRYY, and ADRX, and the EQ and vegetated land descended in turn (Figures 4a-c and 5).The spatial changes of annual FVC max and LAI max from 2000 to 2010 basically followed the same spatial pattern as distributions-gradually decreasing from southeast to northwest overall-but differed in some areas.For example, the eastern HGRLP, central and southwest of MSRNC, where lower FVC max and LAI max values were concentrated, had positive changes; the northwestern ADRX, where higher FVC max and LAI max values were found, however, was getting worse (Figure 4a,b,d,e).In relation to annual NPP avg change, the slope declined from the center (eastern HGRLP and northern MHRSCHH) to the periphery in general, and the southern subregions suffered a more severe loss of NPP avg , whereas the northern subregions had improved (Figure 4c,f).There were significant temporal trends in changes in annual FVC max and LAI max over the whole GGP area from 2000 to 2010, reaching 0.1459% yr −1 (R 2 = 0.3964, P < 0.05) and 0.0121 yr −1 (R 2 = 0.6672, P < 0.01), respectively, although this was not so for the NPP avg , which reached 2.6958 gC m −2 yr −1 (R 2 = 0.1910, P > 0.05) (Figure 5, Table 5).The annual FVC max , LAI max , and NPP avg changes, in general, all showed a slightly increasing trend in the decade, indicating that the EQ of the restoration region had improved.Large variations in the temporal trends of annual FVC max , LAI max , and NPP avg existed across the whole GGP region (Figure 5, Table 5).The annual FVC max and LAI max in most subregions remained basically stable with minor fluctuations and an overall slight increase, except FVC max in ADRX and LAI max in MHRHG.The increase in annual FVC max mainly occurred in HGRLP, with the highest rate of 0.5083% yr −1 (R 2 = 0.7244, P < 0.01), followed by MHRHG, LMHRMLYP, MHRSCHH, and YGPR, with significant slopes of over 0.2% yr −1 (P < 0.05) (Figures 4d and 5a).The annual LAI max mainly increased in MHRHG, with the highest slope of 0.0724 yr −1 (R 2 = 0.6425, P < 0.01), followed by LMHRMLYP and HGRLP, with slopes of over 0.02 yr −1 (P < 0.01).The EQ in these three regions rose very significantly (Figures 4e and 5b).The annual NPP avg change showed relatively large fluctuations in the MHRHG, LMHRMLYP, MHRSCHH, and YGPR, but oscillated with micro amplitudes in other subregions.The annual NPP avg in most subregions rose in general, especially in the HGRLP (with the highest slope of 8.0381 gC m −2 yr −1 , R 2 = 0.7950, P < 0.01), LMHRMLYP, MHRSCHH, and MSRNC, while the annual NPP avg in the ADRX, AGRSC, and YGPR decreased (Figures 4f and 5c).

Spatial Pattern and Temporal Trend of Key Ecosystem Services
The annual WR, SWAL, and SWIL were estimated to be 108.76× 10 4 t km −2 , 15.62 t ha, and 16.58 t ha per unit area per year, respectively, at total quantities of 580.56 billion tons, 11.57 billion tons, and 12.27 billion tons in the whole GGP area from 2000 to 2010 (Table 6).Figure 6 demonstrates that all the three key ESs displayed great spatial heterogeneity in distribution, with the southeast subregions generally performing much better than the northwest.The annual WR pua decreased gradually from southeast to northwest, while the annual SWAL pua increased from northeast to southwest in total, and the SWIL basically occurred in Northwest China (Figure 6a-c).With respect to spatial changes of ESs, the southern subregions, which contained most of the WR owing to abundant precipitation and high vegetation cover, showed weaker WR from 2000 to 2010, while the northern subregions with poor WR showed a positive trend (Figure 6a,d).Soil water erosion was alleviated effectively all over the GGP region, especially in the HGRLP, AGRSC, and YGPR, where the severest SWAL occurred due to the Loess Plateau and alpine valleys, but the situation remained poor in the southern edges of the ADRX and ASMRYY (Figure 6b,e).It was gratifying that positive changes in soil wind erosion had also occurred in the regions with heavy SWIL, i.e., the middle of the ADRX and western ASRNC, but in some regions-the northwest of ASMRYY and HGRLP-the deterioration continued.It should be noted that a slight increase in soil wind erosion appeared in a few areas which seldom had SWIL (eastern ASRNC and parts of south subregions) (Figure 6c,f).
Unlike EQ, the temporal trends of the three key ESs were not very significant.The slopes of annual WR pua , SWAL pua , and SWIL pua changes were −2283 t km −2 yr −1 , −0.0841 t ha yr −1 , and −1.0071 t ha yr −1 , respectively (R 2 = 0.4570, P < 0.05) (Figure 7, Table 6).The declines in SWAL pua and SWIL pua implied that the ESs of SC and SP in the GGP area from 2000 to 2010 had been enhanced overall (especially SP), whereas the improvement of WR was not reflected clearly and tended to exhibit a slight weakening.However, the temporal trend of ESs on the subregion-scale differed from the whole region, especially the WR pua .Half of the subregions, mostly in the north, appeared to show steady rises.In the southern subregions, the annual WR pua decreased overall, with obvious fluctuations (Figure 7a, Table 6).With 56.93 t ha on average, the annual SWAL pua in AGRSC was about two to 24 times higher than any other subregions and had much bigger fluctuations.The annual SWAL pua dropped with minor fluctuations in most areas, and significantly in the YGPR (P < 0.05) and HGRLP (P < 0.01), although there was an increase in three northern subregions: ADRX, MSRNC, and ASMRYY (Figure 7b, Table 6).The annual SWIL pua decreased very significantly in ADRX (P < 0.01), MSRNC (P < 0.01), AGRSC (P < 0.01), and LMHRMLYP (P < 0.05); while that in ASRNC, HGRLP, and ASMRYY showed a downward trend at first, but later an upward trend; and HGRLP and ASMRYY eventually rose slightly on the whole (Figure 7c, Table 6).The annual WR pua , SWAL pua , and SWIL pua in ASMRYY all showed rising trends, which were the opposite to those in the whole GGP area.Table 6.Statistical parameters of ecosystem services indicators in the Grain for Green Program (GGP) area from 2000 to 2010.The units of the mean and slope of WR pua are 10 4 t km −2 and 10 4 t km −2 yr −1 , respectively, while for both SWAL pua and SWIL pua , units are t ha and t ha yr −1 , respectively.Note: the * symbol on the top-right of the P value means the tendency is significant (P < 0.05), and the ** symbol means the tendency is very significant (P < 0.01).

Ecosystem
poor in the southern edges of the ADRX and ASMRYY (Figure 6b,e).It was gratifying that positive changes in soil wind erosion had also occurred in the regions with heavy SWIL, i.e., the middle of the ADRX and western ASRNC, but in some regions-the northwest of ASMRYY and HGRLP-the deterioration continued.It should be noted that a slight increase in soil wind erosion appeared in a few areas which seldom had SWIL (eastern ASRNC and parts of south subregions) (Figure 6c,f).Unlike EQ, the temporal trends of the three key ESs were not very significant.The slopes of annual WRpua, SWALpua, and SWILpua changes were −2283 t km -2 yr -1 , −0.0841 t ha yr -1 , and −1.0071 t ha yr -1 , respectively (R 2 =0.4570,P<0.05) (Figure 7, Table 6).The declines in SWALpua and SWILpua implied that the ESs of SC and SP in the GGP area from 2000 to 2010 had been enhanced overall (especially SP), whereas the improvement of WR was not reflected clearly and tended to exhibit a slight weakening.However, the temporal trend of ESs on the subregion-scale differed from the whole region, especially the WRpua.Half of the subregions, mostly in the north, appeared to show steady rises.In the southern subregions, the annual WRpua decreased overall, with obvious fluctuations (Figure 7a, Table 6).With 56.93 t ha on average, the annual SWALpua in AGRSC was about two to 24 times higher than any other subregions and had much bigger fluctuations.The annual SWALpua dropped with minor fluctuations in most areas, and significantly in the YGPR (P<0.05) and HGRLP (P<0.01) , although there was an increase in three northern subregions: ADRX, MSRNC, and ASMRYY (Figure 7b, Table 6).The annual SWILpua decreased very significantly in ADRX (P<0.01),MSRNC (P<0.01),AGRSC (P<0.01), and LMHRMLYP (P<0.05); while that in ASRNC, HGRLP, and ASMRYY showed a downward trend at first, but later an upward trend; and HGRLP and ASMRYY eventually rose slightly on the whole (Figure 7c, Table 6).The annual WRpua, SWALpua, and SWILpua in ASMRYY all showed rising trends, which were the opposite to those in the whole GGP area.Table 6.Statistical parameters of ecosystem services indicators in the Grain for Green Program (GGP) area from 2000 to 2010.The units of the mean and slope of WRpua are 10 4 t km -2 and 10 4 t km -2 yr -1 , respectively, while for both SWALpua and SWILpua, units are t ha and t ha yr -1 , respectively.

Trade-Offs and Synergies of Ecological Indicators in GGP Area
Large trade-offs have occurred between ecosystem services, largely due to land-use and land management changes [75].Combining multiple ecological outputs under the scenario of GGP demonstrates the extent of the tradeoffs or synergies among these ecological variables.In the whole GGP region, most pair-wise ecological indicators were revealed to be remarkably correlated (Figure 8).FVC, LAI, and NPP showed a strong and significant positive correlation with each other (P < 0.05), presenting synergies among EQ variables.With respect to ESs, WR positively correlated with SWAL (P < 0.05), indicating a remarkable trade-off with SC.The relatively positive correlation between SWIL and WR and between SWIL and SWAL, implied a trade-off between SP and WR and a synergy between SP and SC.While SWAL and SWIL were negatively correlated with relative EQ indicators, presenting synergies of ESs of SC and SP with EQ, WR showed no apparent tendency with EQ.Since the FVC and LAI respectively parameterize the horizontal and vertical density of vegetation [31], and NPP is the rate of net accumulation of organic matter in plants [76], with the growth of vegetation, they increase simultaneously in general.Vegetation generally reduces runoff and soil water erosion and thus enhances soil retention by intercepting rainfall, increasing water infiltration, increasing soil shear strength, reducing raindrop energy and 'splash' effects, and trapping sediment where rainfall is abundant [77][78][79], and halts near-surface sand flow, hence improving sand fixation by stabilizing soil, increasing underlying surface roughness, reducing near-surface wind speed, and weakening sand transportation intensity in arid and semi-arid regions [80][81][82].

Impact of Climate Change on Ecological Change in the GGP Subregions
Climate variability and ecological restoration programs jointly influence vegetation growth [78], thus further impacting EQ and ESs.From 2000 to 2010, the climate in the whole GGP area became warmer and wetter, with average increments of 0.0292°C and 0.4685 mm per year in temperature and Our findings of no apparent tendency between WR and EQ seemed to be slightly different from the discovery of Jia et al. [23], who reported that the water yield in North Shaanxi showed no apparent tendency with NPP, but declined with the increase of biomass after eight years of GGP implementation, and contrary to Cao et al [83] and Feng et al [26], who claimed that there was a trade-off between the water yield and afforestation in the Loess Plateau.These might be attributable to the larger scale of our study area.The relationship between a pair of ESs can differ across different scales and across different socio-ecological systems [84,85].The whole GGP region represented vast spatial heterogeneity of natural characteristics and climate, while the previous work focused on arid or semi-arid areas.Due to the fact that the debate about the balance of benefits and disbenefits of forestry in water regulation is still widespread around the planet [86][87][88][89], the effects of climate and geographic locations on water conservation need further consideration.

Impact of Climate Change on Ecological Change in the GGP Subregions
Climate variability and ecological restoration programs jointly influence vegetation growth [78], thus further impacting EQ and ESs.From 2000 to 2010, the climate in the whole GGP area became warmer and wetter, with average increments of 0.0292 • C and 0.4685 mm per year in temperature and precipitation, respectively.Climate change showed large spatial heterogeneity, with wetter places (mainly southern subregions) becoming drier, and drier places (mainly northern subregions) becoming wetter.Most of the southern and northwestern subregions became warmer, especially in the southwest, while central and northeast subregions became colder (Figure 9).We used the slope of annual FVC max , LAI max , NPP avg , WR pua , SWAL pua , and SWIL pua changes and the slope of annual precipitation and temperature changes to investigate the relationship between the ecological changes and climate change at a subregion scale (Figure 10).We found that the impacts of precipitation and temperature on ecological changes at a subregion level were not significant, except WR pua , which was impacted by precipitation very significantly (R 2 = 0.78, P < 0.01) (Figure 10d).Therefore, the increase in precipitation might explain the rise of WR pua in most northern subregions, while the decrease led to a decline in WR pua in most southern subregions.However, the WR pua and precipitation of subregions did not fluctuate with the same volatility and the response of WR pua to precipitation change in northern subregions was weaker than that in southern subregions, e.g., the slope values for precipitation change in the ASRNC, HGRLP, and MHRSCHH were similar, but the slopes for WR pua change in the ASRNC and HGRLP were far lower than in MHRSCHH (Figure 10d).Considering the abnormal rise of annual FVC max , LAI max , and NPP avg compared with climate change in northern subregions, especially in the HGRLP (Figure 10a-c), the increased forest might lead to a reduction in the water yield in arid and semi-arid regions [78,83,90], whereas the precipitation might become a more constraining factor to water conservation in humid areas.New planting would increase ET and induce a significant decrease in the ratio of river runoff to annual precipitation across hydrological catchments [26], and large-scale afforestation in vulnerable arid and semi-arid regions could increase the severity of water shortages [83].Besides WR pua , the SWIL pua in the ASMRYY and HGRLP still increased slightly under the scenarios of afforestation (Figure 10f), concurring with the discovery of Cao [22], who found that planting trees in arid and semi-arid regions had failed to mitigate desertification.The correlation was relatively high between SWAL pua and precipitation (Figure 10e), indicating that climate change contributes to the enhancement of SC.The fact that the effects of climate change on most ecological variables at a subregion level were not significant implied that the GGP might play a more important role in ecological change, and the difference of ecological outputs among subregions might be attributed to different regional restoration measures (Table 1).The selection and richness of tree species, area, and locations for revegetation or afforestation would significantly influence the ecological effects [91,92].
indicating that climate change contributes to the enhancement of SC.The fact that the effects of climate change on most ecological variables at a subregion level were not significant implied that the GGP might play a more important role in ecological change, and the difference of ecological outputs among subregions might be attributed to different regional restoration measures (Table 1).The selection and richness of tree species, area, and locations for revegetation or afforestation would significantly influence the ecological effects [91,92].

The Contribution of the GGP Implementation to Ecological Change
Climate change can enhance the ecological condition to some extent, but the GGP only affected new vegetation in the GGP area.To find out the contribution of GGP implementation to the ecological benefit, a paired Welch's t-test was conducted to compare the slope of annual FVCmax, LAImax, NPPavg, WRpua, SWALpua, and SWILpua changes between two LULC types: (1) the restored vegetation (including forest and grassland) transferred from another LULC type from 2000 to 2010 (LULC A);

The Contribution of the GGP Implementation to Ecological Change
Climate change can enhance the ecological condition to some extent, but the GGP only affected new vegetation in the GGP area.To find out the contribution of GGP implementation to the ecological benefit, a paired Welch's t-test was conducted to compare the slope of annual FVC max , LAI max , NPP avg , WR pua , SWAL pua , and SWIL pua changes between two LULC types: (1) the restored vegetation (including forest and grassland) transferred from another LULC type from 2000 to 2010 (LULC A); and (2) the retained forest and grassland that existed before 2000 and remained untransformed from 2000 to 2010 (LULC B).The sample points of LULC A (denoted as Sample A) were selected from the pixels that were transferred to forest and grassland in 2000-2010 at a 1 km spatial resolution.The sample points of LULC B (denoted as Sample B) were randomly selected from the intersecting pixels of forest and grassland in the 2000, 2005, and 2010 LULC dataset in the nearest neighborhood of Sample A, to ensure similar climate conditions and further exclude the influence of climate change.Both samples had about 15000 points, involving 56% forest and 44% grassland.The null hypothesis and alternative hypothesis are listed in Table 7. Results showed that the slope of annual FVC max , LAI max , NPP avg , and SWAL pua changes passed the paired Welch's t-test at a very significant level, whereas the slope of annual precipitation and temperature changes did not, indicating that the recovery/improvement rates of EQ and SC in the revegetation area were higher than existing vegetation under similar climate conditions.Therefore, the GGP implementation rather than climatic change was the primary cause for the estimated increase of ecosystem quality and services, which is supported by many studies [15,19,26].The slope of annual FVC max , LAI max , and NPP avg changes in LULC A increased by 14.23%, 9.94%, and 8.23%, respectively, compared with those in LULC B, while the slope of annual SWAL pua change in LULC A decreased by 30.45% compared with LULC B. This change should be considered as the contribution of the GGP implementation.

Alternative Hypothesis (H 1 )
µ Note: the * symbol on the top-right of the P value means the tendency is significant (P < 0.05), for ** symbol means the tendency is very significant (P < 0.01).
As the slope of annual WR pua and SWIL pua changes did not pass the t-test, we further conducted paired Welch's t-tests on them between LULC A and LULC B in forest and grassland separately, and found that the slope of SWIL pua in restored grassland passed the t-test with high significance and it decreased by 18.05% more than that in LULC B grassland (µ A = −0.2260,µ B = −0.1914,P < 0.01).The slope of annual SWIL pua change in restored forest and the slope of annual WR pua change in neither restored forest nor grassland did not pass the t-test.The results further confirmed our findings in the last section.Afforestation with inappropriate species in arid and semi-arid regions can damage local water balances because of the high consumption of plant transpiration and evaporation from the soil surface [93,94], and thus aggravate water shortages and soil degradation [95], leading to erosion and further exacerbating desertification in the long term [89].In this case, introducing stable communities of natural desert steppe and grassland vegetation is a more appropriate way to sustain the natural ecosystem and control the spread of desertification [96].

Implications for the New Round of the GGP
To investigate the gap between new plantings and retained vegetation in improving ecological conditions, a paired Welch's t-test was conducted to compare the annual average FVC max , LAI max , NPP avg , WR pua , SWAL pua , and SWIL pua between LULC A and LULC B from 2000 to 2010.The null hypothesis and alternative hypothesis are listed in Table 8.All the indicators except SWIL pua passed the t-test at a very significant level.The annual average FVC max , LAI max , NPP avg , and WR pua in LULC A were respectively 2.56%, 6.74%, 8.81%, and 22.34% less than those in LULC B, and the annual average SWAL pua and SWIL pua were respectively 7.39% and 5.22% higher than those in LULC B, suggesting that the new planting had not reached its potential maximum in terms of ecosystem condition enhancement.Therefore, as the new round of GGP is initiated, more attention should be paid to the protection and restoration of old-growth original forest, the cultivation of young planted forest, and the management of artificial forests, apart from only afforestation and reforestation.To achieve this goal, besides the institutional arrangements such as legislation, setting up nature reserves and enhancing incentive payments, which the government always does, market tools, the development of new technology, and international programs should be taken into account.Responsible forest management could be promoted by bringing in market-based third-party certification for evaluating forests according to a set of standards [97,98].Improving the biological, physiological, and genetic quality of forest planting stock through plant biotechnology, a high level of environmental control in the growth chambers, and transnational co-operation are also practical ways to carry out forest restoration projects [99][100][101][102].International climate change mitigation programs, which are designed to reduce emissions from deforestation and forest degradation in developing countries (REDD+), should also be taken into consideration as collateral payments for the substantial carbon benefit in maintaining forest [103].Note: the * symbol on the top-right of the P value means the tendency is significant (P < 0.05), and the ** symbol means the tendency is very significant (P < 0.01).
Our studies suggested that not all the ESs increased simultaneously under the scenario of reforestation and afforestation, and significant differences and lags existed between the improvement of each pair of ecological indicators.This implied the existence of trade-offs and synergies between ESs, which complicates the restoration mechanism and must be considered in policy-making [19,23].Considering the more important role GGP played in ecological change and the variation of ecological change among subregions, regional restoration measures need to be concerned carefully.More attention should be paid to re-establishing near-natural forests or the restoration of a region's natural climax vegetation [22,96,[104][105][106].
Owing to the urgent need to address eco-environmental problems, the Chinese government had to act before the achievement of scientific consensus, which needs a step-by-step investigation and does not provide instant solutions [107].Therefore, the massive afforestation program which ignored differences in topography, climate, and hydrology in the first round of the GGP is understandable.In future GGP policy-making and implementation, we suggest that the former practice of aggressive afforestation on a large-scale must be changed.Rational and specific plans and designs for afforestation areas should seek scientific advice before implementation according to regional climate and geographical characteristics, including local temperatures, precipitation, terrain, soil, and hydrology.In addition to the appropriateness of tree species, near-natural vegetation communities which take into account diversity and mixtures of the tree-shrub-herb system should be planted instead of single-species plantations.

Figure 1 .
Figure 1.Distribution of subregions involved in the Grain for Green Program (GGP).For the full names and description of each region, see Table 1.

Figure 1 .
Figure 1.Distribution of subregions involved in the Grain for Green Program (GGP).For the full names and description of each region, see Table 1. 2000

Figure 2 .
Figure 2. The workflow used to generate the 2000 and 2005 land use and land cover datasets.

Figure 2 .
Figure 2. The workflow used to generate the 2000 and 2005 land use and land cover datasets.

Figure 3 .
Figure 3. Spatial pattern of ecosystem distribution and transformation in the Grain for Green Program (GGP) area.(a-c) are the ecosystem distributions in 2010, 2005, and 2010, respectively; (d) is the ecosystem transformation from 2000 to 2010.

Figure 3 .
Figure 3. Spatial pattern of ecosystem distribution and transformation in the Grain for Green Program (GGP) area.(a-c) are the ecosystem distributions in 2010, 2005, and 2010, respectively; (d) is the ecosystem transformation from 2000 to 2010.

Figure 4 .
Figure 4. Spatial pattern of annual average ecosystem quality (EQ) and slope of EQ change in the Grain for Green Program (GGP) area from 2000 to 2010.(a), (b), and (c) are the annual average FVCmax, LAImax, and NPPavg, respectively; (d), (e), and (f) refer to the slope of annual FVCmax, LAImax, and NPPavg changes, respectively.

Figure 5 .
Figure 5. Temporal trends of annual ecosystem quality (EQ) changes in the Grain for Green Program (GGP) area from 2000 to 2010.(a), (b), and (c) are the temporal trends of annual FVCmax, LAImax, and NPPavg changes, respectively.

Figure 4 .
Figure 4. Spatial pattern of annual average ecosystem quality (EQ) and slope of EQ change in the Grain for Green Program (GGP) area from 2000 to 2010.(a-c) are the annual average FVC max , LAI max , and NPP avg , respectively; (d-f) refer to the slope of annual FVC max , LAI max , and NPP avg changes, respectively.

Figure 4 .
Figure 4. Spatial pattern of annual average ecosystem quality (EQ) and slope of EQ change in the Grain for Green Program (GGP) area from 2000 to 2010.(a), (b), and (c) are the annual average FVCmax, LAImax, and NPPavg, respectively; (d), (e), and (f) refer to the slope of annual FVCmax, LAImax, and NPPavg changes, respectively.

Figure 5 .
Figure 5. Temporal trends of annual ecosystem quality (EQ) changes in the Grain for Green Program (GGP) area from 2000 to 2010.(a), (b), and (c) are the temporal trends of annual FVCmax, LAImax, and NPPavg changes, respectively.

Figure 5 .
Figure 5. Temporal trends of annual ecosystem quality (EQ) changes in the Grain for Green Program (GGP) area from 2000 to 2010.(a-c) are the temporal trends of annual FVC max , LAI max , and NPP avg changes, respectively.

Figure 6 .
Figure 6.Spatial pattern of annual average ecosystem services (ESs) and slope of ESs change in the Grain for Green Program (GGP) area from 2000-2010.(a), (b), and (c) are the annual average WRpua, SWALpua, and SWILpua, respectively; (d), (e), and (f) refer to the slope of annual WRpua, SWALpua, and SWILpua changes, respectively.

Figure 6 .
Figure 6.Spatial pattern of annual average ecosystem services (ESs) and slope of ESs change in the Grain for Green Program (GGP) area from 2000-2010.(a-c) are the annual average WR pua , SWAL pua , and SWIL pua , respectively;(d-f) refer to the slope of annual WR pua , SWAL pua , and SWIL pua changes, respectively.

Figure 7 .
Figure 7. Temporal trends of annual ecosystem services (ESs) changes in the Grain for Green Program (GGP) area from 2000 to 2010.(a), (b), and (c) are the temporal trends of annual WRpua, SWALpua, and SWILpua changes, respectively.

Figure 7 .
Figure 7. Temporal trends of annual ecosystem services (ESs) changes in the Grain for Green Program (GGP) area from 2000 to 2010.(a-c) are the temporal trends of annual WR pua , SWAL pua , and SWIL pua changes, respectively.

27 Figure 8 .
Figure 8. Tradeoffs and synergies among variables of EQ and ESs.

Figure 8 .
Figure 8. Tradeoffs and synergies among variables of EQ and ESs.

Figure 9 .
Figure 9. Spatial pattern of annual average climate and slope of climate change in the Grain for Green Program (GGP) area from 2000 to 2010.(a) and (b) are the annual average precipitation and temperature, respectively; (c) and (d) refer to the slope of annual precipitation and temperature changes, respectively.

Figure 9 . 27 Figure 10 .
Figure 9. Spatial pattern of annual average climate and slope of climate change in the Grain for Green Program (GGP) area from 2000 to 2010.(a,b) are the annual average precipitation and temperature, respectively; (c,d) refer to the slope of annual precipitation and temperature changes, respectively.Remote Sens. 2019, 11, x FOR PEER REVIEW 20 of 27

Figure 10 .
Figure 10.Relationship between the difference in ecological changes among Grain for Green Program (GGP) subregions and climate condition from 2000 to 2010.(a-f) are the slope of annual FVC max , LAI max , NPP avg , WR pua , SWAL pua , and SWIL pua changes compared with the slope of annual precipitation and temperature changes from 2000 to 2010, respectively.

Table 7 .
Null hypothesis and alternative hypothesis for the slope of changes of climate, ecosystem quality, and ecosystem services indicator values between Sample A and Sample B from 2000 to 2010.µ A and µ B are, respectively, the mean of the slope value of Sample A and Sample B, and the units from top to bottom are mm yr −1 , • C yr −1 , % yr −1 , yr −1 , gC m −2 yr −1 , 10 4 t km −2 yr −1 , t ha yr −1 , and t ha yr−1.

Table 1 .
Natural characteristics and restoration measures in the Grain for Green Program (GGP) subregions.

Table 2 .
The comprehensive assessment indicator system for ecological changes of the Grain for Green Program (GGP).

Table 3 .
The ecosystem composition and distribution in the Grain for Green Program (GGP) area in 2000, 2005, and 2010.

Table 5 .
Statistical parameter values of ecosystem quality indicators in the Grain for Green Program (GGP) area from 2000 to 2010.The units of the mean and slope of FVC max are % and % yr −1 , respectively; for NPP avg , units are gC m −2 and gC m −2 yr −1 , respectively; and the LAI max is dimensionless.

Table 8 .
Null hypothesis and alternative hypothesis of mean annual values of ecosystem quality and ecosystem services indicators between Sample A and Sample B from 2000-2010.µ A and µ B are, respectively, the mean of annual values of indicators of Sample A and Sample B, and the units from top to bottom are %, dimensionless, gC m −2 , 10 4 t km −2 , t ha, and t ha.